flyopunta/src/src_par/partres5.cpp

/***************************************************************************
                          partres5  -  description
                             -------------------
    begin                : Sat Feb 26 2005
    copyright            : (C) 2005 by Giuseppe Pierazzini
    email                : peppe@peppix
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/
#include "flyoh.h"
#include "fase3max.h"
using namespace std;
Partres5::Partres5()
{
  p1=p2=p3=p4=p5=0;
//  Gout<<"\n Class Partres5  : "<<nome<<" done!";
}


void Partres5::

DecParm()
{
  // rivedere !!!!

// skip.. if  p1 allready defined....
  if ( p1==NULL )
  {
// define the pointer to the four particles;
// and the decay parameters

    Get_PhaseParm();
//  Get_MatrixParm();

  }

//------------------------
// get the invariant values...
  Get_PhaseSpace();

// define the resonace particle parameters..masse...
//
// res 1  == this
  Set_Massa (mres);
  P.Riposo  (mres);
// res 2
  pr2->Set_Massa(mres2);
  pr2->P.Riposo (mres2);
// res 3
  pr3->Set_Massa(mres3);
  pr3->P.Riposo (mres3);
}


void Partres5::Get_PhaseParm()
{
  /*-------------------------------------------
   Five bodies decay  k -> a b c d e  in phase space..
   the decay is  via the Res4 resonance
   generated in the decay chain automatilly...see Make_list
  -----
   the chain is K->a Res4 
                    ->b Res3 
                       ->c Res2 
                          ->d  e
  -----
   Phase space weight
   w = F(padre,Res4,e)*F(Res4,Res,v)*F(Res,Res',pi0)*F(Res,pi0,pi0)
  ------------
  Note: imput pattern  K =  a b c d e 
  gives the following pointer matching
   p1->a, p2->b, p3->c, p4->d, p5->e
   
   pr1->Res4, pr2->Res3, pr3->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!=NULL ) return;
  
    p1=lf;        
    p2=dw;        
    p3=p2->rg->dw;
    p4=p3->rg->dw;
    p5=p4->rg;    
  
 
//resonances... generate come particelle
    pr1=this;
    pr2=dw->rg;
    pr3=pr2->dw->rg;

// 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;
    mchd4 =p4->Get_Massa(); mchd4q=mchd4*mchd4;
    mchd5 =p5->Get_Massa(); mchd5q=mchd5*mchd5;

// massimi e minimi valori delle risonanze
    mmin     =mchd2+mchd3+mchd4+mchd5;  mminq =mmin*mmin;
    mmin2=          mchd3+mchd4+mchd5;  mminq2=mmin2*mmin2;
    mmin3=                mchd4+mchd5;  mminq3=mmin3*mmin3;
    mmax=mpadre-mfrate;                 mmaxq =mmax*mmax;
    mmax2=mmax-mchd2;                   mmaxq2=mmax2*mmax2;
    mmax3=mmax-mchd2-mchd3;             mmaxq3=mmax3*mmax3;
    dmq =mmaxq-mminq;
    dmq2=mmaxq2-mminq2;
    dmq3=mmaxq3-mminq3;
    double pass1=dmq/100.;
    double pass2=dmq2/100.;
    double pass3=dmq3/100.;
    double w1,w2,w3;

//     cout<<"\n mass  "<<mminq<<" "<<mmaxq<<" "<<pr2->mminq<<" "<<pr2->mmaxq<<" "<<pr3->mminq<<" "<<pr3->mmaxq;
   int p_id[6]= {up->Get_Idm(),p1->Get_Idm(),p2->Get_Idm(),p3->Get_Idm(),p4->Get_Idm(),p5->Get_Idm()};
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()<<" "<<p4->Get_Nome()<<" "<<p5->Get_Nome();
      return;
    }
  }


    // define the wfase_max for p.s.  normalizazion
    wfase_max=0.0;
    for ( mresq=mminq; mresq<mmaxq; mresq+=pass1 )
    {
      mres=sqrt ( mresq );
      w1=Fspace ( mpadre,mres,mfrate );
      for ( mresq2=mminq2; mresq2<mmaxq2; mresq2+=pass2 )
      {
        mres2=sqrt ( mresq2 );
        w2=Fspace ( mres,mres2,mchd2 );
        for ( mresq3=mminq3; mresq3<mmaxq3; mresq3+=pass3 )
        {
          mres3=sqrt ( mresq3 );
          w3=Fspace ( mres2,mres3,mchd3 ) *Fspace ( mres3,mchd4,mchd5 );
          w=w1*w2*w3;
          if ( w>wfase_max ) wfase_max=w;
        }
      }
    }

    Gout<<"\n Phase-Space:w="<<wfase_max<<" done for " <<up->Get_Nome()<<" ==> "<<p1->Get_Nome()<<" "<<p2->Get_Nome() 
    <<"  "<<p3->Get_Nome()<<" "<<p4->Get_Nome()<<" "<<p5->Get_Nome();  

    new Fase3Max ( p_id,wfase_max );
}

void  Partres5::Get_PhaseSpace()
{
  iterfs=0; w=0.5;
  while ( iterfs++<300 )
  {
    // resonance mass selection...
    mres     =sqrt ( ( mresq =mminq+Pran() *dmq ) );
    mres2    =sqrt ( ( mresq2=mminq2+Pran() *dmq2 ) );
    if ( mres2+mchd2>mres )  continue;
    mres3    =sqrt ( ( mresq3=mminq3+Pran() *dmq3 ) );
    if ( mres3+mchd3>mres2 ) continue;

    w = Fspace ( mpadre, mres,mfrate ) *Fspace ( mres, mres2,mchd2 );
    w*= Fspace ( mres2,mres3,mchd3 ) *Fspace ( mres3,mchd4,mchd5 );


    if ( w >wfase_max )
    {
      Gout<<"\n Wfas off: Ev "<<evento_.Gen.Event<<" iter="<<iterfs<<" 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(), p4->Get_Idm(),p5->Get_Idm()};
      Fasemx->Set_Wmax ( p_id, wfase_max );
    }

    if ( Pran() <w/wfase_max ) break;

  }
  if ( iterfs>400 )
    Gout<<"\nPartRes5 iter="<<iterfs<<" Risonaze: m2345="<<mres<<" m345="<<mres2<<" m45="<<mres3;
}