flyopunta/src/src_par/partres4.cpp

/***************************************************************************
                          partres4.cpp  -  description
                             -------------------
    begin                : Mon Jun 18 2001
    copyright            : (C) 2001 by Giuseppe Pierazzini
    email                : peppe@unipi.it
  ***************************************************************************
  *                                                                         *
  *   NA48  simulation program.                                             *
  *                                                                         *
  ***************************************************************************/


#include "flyoh.h"
#include "fase3max.h"
using namespace std;
Partres4::Partres4()
{p1=p2=p3=p4=p5=0;
Gout<<"\n Class Partres4  : "<<nome<<" done!";
}


void Partres4::DecParm()
{     
// skip.. if  p1 allready defined....

  if(p1==0)
{ 
// 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);
  
  
  
 }


void Partres4::Get_PhaseParm()
{
/*-------------------------------------------
 Four bodies decay  k -> a b c d  in phase space..
 the decay is  via the Res3 resonance 
 generated in the decay chain automatilly...see Make_list 
-----
 the chain is K->a Res3 
                    ->b Res
                         ->c d 
-----
 Phase space weight 
 w = F(padre,Res3,a)*F(Res3,Res,b)*F(Res,c,d)
------------
Note: imput pattern like  K = a b c d  ( esempio:e v pi0 pi0) 
gives the following pointer matching
 p1->a, p2->b, p3->c, p4->d
 pr1->Res3, pr2->Res,
 
 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;
     p1=lf;      
     p2=dw;   
     p3=dw->rg->dw;
     p4=p3->rg;      
    
//resonances... generate come particelle 
     pr1=this;
     pr2=dw;
// 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;
     
// massimi e minimi valori delle risonanze 
     mmin     =mchd2+mchd3+mchd4;  mminq =mmin*mmin;
     mmin2=          mchd3+mchd4;  mminq2=mmin2*mmin2;
     mmax=mpadre-mfrate;                 mmaxq =mmax*mmax;
     mmax2=mmax-mchd2;                   mmaxq2=mmax2*mmax2;
     dmq =mmaxq-mminq;
     dmq2=mmaxq2-mminq2;

// nota Fase 

    int p_id[6]= {up->Get_Idm(),p1->Get_Idm(),p2->Get_Idm(),p3->Get_Idm(),p4->Get_Idm(),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()<<" "<<p4->Get_Nome();
      return;
    }
  }




     double pass1=dmq/100.;
     double pass2=dmq2/100.;
     double w1;
        
//     cout<<"\n mass  "<<mminq<<" "<<mmaxq<<" "<<pr2->mminq<<" "<<pr2->mmaxq<<" "<<pr3->mminq<<" "<<pr3->mmaxq;
     
  // 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);
        w=Fspace(mres,mres2,mchd2)*Fspace(mres2,mchd3,mchd4);
         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();  
         
      new Fase3Max ( p_id,wfase_max );
 }

void  Partres4::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;
     w = Fspace(mpadre, mres,mfrate)*Fspace(mres, mres2,mchd2);
     w*= Fspace(mres2,mchd3,mchd4);      
   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(),0};
      Fasemx->Set_Wmax ( p_id, wfase_max );
    }
    
    if(Pran()<w/wfase_max) break;
    
   }
   if(iterfs>300) 
   Gout<<"\nPartRes5 iter="<<iterfs<<" Risonaze: m234="<<mres<<" m34="<<mres2;
   
 }