flyopunta/src/src_par/partfase.cpp

/***************************************************************************
                          partfase.cpp  -  description
                             -------------------
    begin                : Mon Dec 25 2000
    copyright            : (C) 2000 by Giuseppe M Pierazzini
    email                : pierazzini@unipi.it
 ***************************************************************************/

#define ESTERNO  1
#define INTERNO  0
#define INBUCO   -1

#define EPSR     1.0E-200
#include "flyoh.h"
#include "flyoana.h"
#include "fase3max.h"
using namespace GmpConst;

//using namespace std;


//-------------------------------------------------

//--------------------------------------------------------------
//             F s p a c e
//--------------------------------------------------------------
//  Fspace    routine

double
Fspace ( double res, double rm1, double rm2 )
{
  double dm,sm;
  sm = ( rm1 + rm2 ) ;
  if ( sm==0.0 )  return 1.0;
  if ( sm>res )   return 0.0;
  sm/=res;
  dm = ( rm1 - rm2 ) /res ;
  return  sqrt ( ( 1.0 - dm * dm ) * ( 1.0 - sm*sm ) );
}
//-----------------------------------------------------
//-----------------------------------------------------
PartFase::PartFase()
{
  type=Res;
  this->res_next =0;
  mpadre=0.;
  p1=p2=p3=p4=p5=0;

  if ( evento_.Gen.Event< 1 ) Gout<<"\n Classe PartFase  : "<<nome<<" done!"<<flush;

}

//--------------------------------------------------------------*/
//                  D e c P a r m                               */
//--------------------------------------------------------------*/
// Get the resonance mass   and decay angles                     Get_mass       */
//

void PartFase::
DecParm()
{
// K -> mfrate mres
//                 ->   mchd1 mchd2
// solo spazio fasi.....

  if ( p1==0 )  // inizialization... only first time
    {
      Get_PhaseParm();   // initialize links and phase parms
      Get_MatrixParm();  // Matrix parms if needed

    }
// -------------------------------

// qui sempre...
// get the fase_space point: e0,e1,e2,mres,mresq,cst
  int itermtr=0;
  while ( itermtr++<200 )
    {
// nota: normalmente in Get_PhaseSpace ci sono 1...6 iterazioni!
      Get_PhaseSpace(); //  Attenzione!  qui ora hai dipsonibili e0,e1,e2, mres,mresq,cst
//   if the matrix prob is ok break to iter !!
      if ( Get_MatrixValue() ) break;  // get matrix value if needed ... here does nothing.
    }
  if ( itermtr>100 )
    Gout<<"\n Ev "<< evento_.Gen.Event
    <<"  Matrix warning: too many iters... chek!  "<<itermtr<<"  "<< SelRea->Get_Titolo();


//-----------------------------------------------------

// define the mass and the quadrimomentum (qvet) of the decay particle.(mres)...
  Set_Massa ( mres );
  P.Riposo ( mres );  //Definizione di P a riposo nel cm di mres

// calcolo la direzione di decadimento dei figli.
// uso i puntatori definiti in Get_PhaseParm() per non sbagliare le particelle...
  fi = DuePiGreco * Pran();
  p2->cstar=cst;
  p2->fistar=fi;
  p3->cstar=-cst;
  p3->fistar=fi+PiGreco;

}

// ---------------------------------

void  PartFase::Get_PhaseSpace()
{
  iterfs=0;

// attenzione nel caso di un urto con produzione di tre corpi, per esempio
// cms -> Kl res2 ; con res2 ->pip pip
// la massa del centro di massa può variare (anche se di poco.. dipende dallo spettro del beam)
// quindi il limte superiore della risonanza   res2 va ricalcolata!
// attenzione  wfase_max, calcolato inizialmente per un valore  della massa invariante di cms,viene via via aggiornato.....

  if ( up->Get_Nome() ==Cms ) //C'è stato un urto!
    {
      mpadre=up->Get_Massa();
      mmax=mpadre-mfrate;
      mmaxq =mmax*mmax;
      dmq =mmaxq-mminq;
    }
//---

  while ( iterfs++<155 )
    {
      // resonance mass selection...
      mres=sqrt ( ( mresq =mminq+Pran() *dmq ) );
      if ( mres+mfrate>mpadre )  continue;
      w= Fspace ( mres, mchd2, mchd3 ) * Fspace ( mpadre, mfrate, mres );
      if ( w >wfase_max )
        {
          Gout<<setprecision ( 5 ) <<"\nEv "<<evento_.Gen.Event<<" iter="<<iterfs<<" W-Ph.Sp. off: prb="<<w
          <<" prbmax="<<wfase_max<<" "<<SelRea->Get_Titolo();
          wfase_max=w+0.000001;
          // aggiorna  wfase_max memorizzato nella lista
          int p_id[6]= {up->Get_Idm(),p1->Get_Idm(),p2->Get_Idm(),p3->Get_Idm(), 0,0};
          Fasemx->Set_Wmax ( p_id, wfase_max );
        }
      if ( Pran() <w/wfase_max ) break;
    }
  if ( iterfs>150 ) Gout<<"\nPh.Sp. warning: Too many iters...="<<iterfs<<" Res2  "<<mres;
// byproduct  mresq,mres,eres,e0,e1,e2,cstù
// Attenzione: remapping 1,2,3 == 0,1,2 ATTENZIONE...

  e0=0.5* ( mpadreq+mfrateq-mresq ) /mpadre;
  cst= 1.0 -  2.0 * Pran();
  eres=mpadre-e0;
  gam = eres/mres;
  eta=sqrt ( e0*e0-mfrateq ) /mres;
  estar1 =0.5* ( mresq+mchd2q-mchd3q ) /mres;
  pstar1=sqrt ( estar1*estar1-mchd2q );
  e1 = gam*estar1 + eta*pstar1*cst;
  e2=eres-e1;

// memorizzo i valori cinematici per decadimenti a tre corpi
// se la particella padre è il beam

  if ( up==SelRea->avo )   Svalues();


}
//---------------------------------------------------

void PartFase::Get_PhaseParm()
{
  /*-------------------------------------------
   Three bodies decay  k -> a b c   in phase space..
   the decay is  via the Res2 resonance
   generated in the decay chain automatilly...see Make_list
  -----
   the chain is K->a Res2
                      ->b c

  -----
   Phase space weight
   w = F(padre,Res2,a)*F(Res2,b,c)
  ------------
  Note: imput pattern like  K = a b c
  gives the following pointer matching
   p1->a, p2->b, p3->c
   note this corrisonde a Res2

   if you change the pattern the chain changes automatically...

   -------- here
   -define the pointers to the  particles;
   -the decay parameters
   -the wfase_max
  --------------------------------------------------*/

  if ( p1!=0 ) return;
  if ( rg !=0 ) p1=rg;
  else p1=lf;   //Gout<<"\n PartFase: "<<p1->Get_Nome()<<endl;
  p2=dw;        //Gout<<"\n PartFase: "<<p2->Get_Nome()<<endl;
  p3=dw->rg;    //Gout<<"\n PartFase: "<<p3->Get_Nome()<<endl;

//resonances... generate come particelle
//     pr1=this;
// particle parameters
  mpadre=up->Get_Massa();
  mpadreq=mpadre*mpadre;
  mfrate=p1->Get_Massa();
  mfrateq=mfrate*mfrate;
  mchd1=mfrate;
  mchd2 =p2->Get_Massa();
  mchd2q=mchd2*mchd2;
  mchd3 =p3->Get_Massa();
  mchd3q=mchd3*mchd3;

// massimi e minimi valori della risonanza
  mmin     =mchd2+mchd3;
  mminq =mmin*mmin;
  mmax=mpadre-mfrate;
  mmaxq =mmax*mmax;
  dmq =mmaxq-mminq;

// controlla se wfase_max è già stato calcolato
  int p_id[6]= {up->Get_Idm(),p1->Get_Idm(),p2->Get_Idm(),p3->Get_Idm(), 0,0};

  if ( Fasemx!=0 )
    {
      double  wmx=Fasemx->Get_Wmax ( p_id );
      if ( wmx>0. )
        {
          wfase_max=wmx;
//      Gout<<"\n Ev "<<Eventi_Fatti<<" F.S. trovato "<<up->Get_Nome() <<" "<<nome<<"=> "<<p1->Get_Nome()
//      <<" "<<p2->Get_Nome() <<" "<<p3->Get_Nome();
          return;
        }
    }
  double pass1=dmq/500.;
//     double w1;
  // define the wfase_max for p.s.  normalizazion
  wfase_max=0.0;
  for ( mresq=mminq; mresq<mmaxq; mresq+=pass1 )
    {
      mres=sqrt ( mresq );
      w=Fspace ( mpadre,mres,mfrate ) *Fspace ( mres,mchd2,mchd3 );
      if ( w>wfase_max ) wfase_max=w;
    }

//  Gout<<"\n Ev "<<Eventi_Fatti<<" ***** Ph.S3: "<<setprecision ( 3 ) <<up->Get_Nome() <<" ["<< nome <<"]->  "
//  <<p1->Get_Nome() <<" "<<p2->Get_Nome() <<" "<<p3->Get_Nome() <<" "<<setw ( 8 ) <<wfase_max;
  Gout<<"\n Phase-Space:w="<<wfase_max<<" done for " <<up->Get_Nome() <<" ==> "<<p1->Get_Nome() <<" "<<p2->Get_Nome()
  <<"  "<<p3->Get_Nome() <<"  Masse range: "<<sqrt ( mminq ) <<" " <<sqrt ( mmaxq );

  new Fase3Max ( p_id,wfase_max );

}


//----------------------------------------------------
//----------- Svalues -----
int PartFase:: Svalues()
{

  // solo nel caso del decadimento a tre corpi...
// get the sister pointers
  Pt=0.0;
  if ( p1!=0 && p1->Get_Nata() !=0 )
    {
      // calcolo il Ptq della particella p1
      //(spesso quella di carica opposta.. attenzione nella posizione in epc)
      gvet Vk=SelRea->avo->Get_Gp().Verso();
      P1=p1->Get_Gp();
      P1.Norma();
      double pro=P1%Vk;
      Pt=sqrt ( P1.normaq-pro*pro );
    }
// Attenzione: remapping 1,2,3 == 0,1,2 ATTENZIONE...
// e0,e1,e2=e0 , s0=msq12,s1=msq20,s2=msq01
// and ... save data in reaction class
  s0=SelRea->ss0=mresq;  //==mpadreq+mchd1q-2.*mpadre*e0;
  s1=SelRea->ss1=mpadreq+mchd2q-2.*mpadre*e1;
  s2=SelRea->ss2=mpadreq+mchd3q-2.*mpadre*e2;
  SelRea->ss_done=1;
// definizione di u,v generati nella struttura evento
//per futura memoria in case of Pipeline on

  evento_.Gen.gss0=s0;
  evento_.Gen.gss1=s1;
  evento_.Gen.gss2=s2;

  evento_.Gen.gvss0=fabs ( s1-s2 ) *0.577350269;//1./sqrt(3.)
  evento_.Gen.gvss1=fabs ( s2-s0 ) *0.577350269;
  evento_.Gen.gvss2=fabs ( s0-s1 ) *0.577350269;

  // randomizzazione degli u,v   ========
  int Irand=int ( 3*Pran() ) %3;
  evento_.Gen.Irand_ss=Irand;

  if ( Irand==0 )
    {
      evento_.Gen.gss_sort=s0;
      evento_.Gen.gvss_sort= evento_.Gen.gvss0;
    }
  else if ( Irand==1 )
    {
      evento_.Gen.gss_sort=s1;
      evento_.Gen.gvss_sort=evento_.Gen.gvss1;
    }
  else if ( Irand==2 )
    {
      evento_.Gen.gss_sort=s2;
      evento_.Gen.gvss_sort=evento_.Gen.gvss2;
    }

  Scrivi_Dati ( 3 );
  return 1;
}