FlyoDoc_2011 Pisa 2011 by GmP --- 011
|
00001 /*************************************************************************** 00002 partres5 - description 00003 ------------------- 00004 begin : Sat Feb 26 2005 00005 copyright : (C) 2005 by Giuseppe Pierazzini 00006 email : peppe@peppix 00007 *************************************************************************** 00008 * * 00009 * NA48 simulation program. * 00010 * * 00011 ***************************************************************************/ 00012 #include "flyoh.h" 00013 #include "fase3max.h" 00014 using namespace std; 00021 Partres5::Partres5() 00022 { 00023 p1=p2=p3=p4=p5=0; 00024 // Gout<<"\n Class Partres5 : "<<nome<<" done!"; 00025 } 00026 00027 00028 void Partres5:: 00029 00030 DecParm() 00031 { 00032 // rivedere !!!! 00033 00034 // skip.. if p1 allready defined.... 00035 if ( p1==NULL ) 00036 { 00037 // define the pointer to the four particles; 00038 // and the decay parameters 00039 00040 Get_PhaseParm(); 00041 // Get_MatrixParm(); 00042 00043 } 00044 00045 //------------------------ 00046 // get the invariant values... 00047 Get_PhaseSpace(); 00048 00049 // define the resonace particle parameters..masse... 00050 // 00051 // res 1 == this 00052 Set_Massa (mres); 00053 P.Riposo (mres); 00054 // res 2 00055 pr2->Set_Massa(mres2); 00056 pr2->P.Riposo (mres2); 00057 // res 3 00058 pr3->Set_Massa(mres3); 00059 pr3->P.Riposo (mres3); 00060 } 00061 00062 00063 void Partres5::Get_PhaseParm() 00064 { 00065 /*------------------------------------------- 00066 Five bodies decay k -> a b c d e in phase space.. 00067 the decay is via the Res4 resonance 00068 generated in the decay chain automatilly...see Make_list 00069 ----- 00070 the chain is K->a Res4 00071 ->b Res3 00072 ->c Res2 00073 ->d e 00074 ----- 00075 Phase space weight 00076 w = F(padre,Res4,e)*F(Res4,Res,v)*F(Res,Res',pi0)*F(Res,pi0,pi0) 00077 ------------ 00078 Note: imput pattern K = a b c d e 00079 gives the following pointer matching 00080 p1->a, p2->b, p3->c, p4->d, p5->e 00081 00082 pr1->Res4, pr2->Res3, pr3->Res2 00083 00084 if you change the pattern the chain changes automatically... 00085 00086 -------- here 00087 -define the pointers to the particles; 00088 -the decay parameters 00089 -the wfase_max 00090 --------------------------------------------------*/ 00091 00092 if ( p1!=NULL ) return; 00093 00094 p1=lf; 00095 p2=dw; 00096 p3=p2->rg->dw; 00097 p4=p3->rg->dw; 00098 p5=p4->rg; 00099 00100 00101 //resonances... generate come particelle 00102 pr1=this; 00103 pr2=dw->rg; 00104 pr3=pr2->dw->rg; 00105 00106 // particle parameters 00107 mpadre=up->Get_Massa(); mpadreq=mpadre*mpadre; 00108 mfrate=p1->Get_Massa(); mfrateq=mfrate*mfrate; //mchd1=mfrate 00109 mchd2 =p2->Get_Massa(); mchd2q=mchd2*mchd2; 00110 mchd3 =p3->Get_Massa(); mchd3q=mchd3*mchd3; 00111 mchd4 =p4->Get_Massa(); mchd4q=mchd4*mchd4; 00112 mchd5 =p5->Get_Massa(); mchd5q=mchd5*mchd5; 00113 00114 // massimi e minimi valori delle risonanze 00115 mmin =mchd2+mchd3+mchd4+mchd5; mminq =mmin*mmin; 00116 mmin2= mchd3+mchd4+mchd5; mminq2=mmin2*mmin2; 00117 mmin3= mchd4+mchd5; mminq3=mmin3*mmin3; 00118 mmax=mpadre-mfrate; mmaxq =mmax*mmax; 00119 mmax2=mmax-mchd2; mmaxq2=mmax2*mmax2; 00120 mmax3=mmax-mchd2-mchd3; mmaxq3=mmax3*mmax3; 00121 dmq =mmaxq-mminq; 00122 dmq2=mmaxq2-mminq2; 00123 dmq3=mmaxq3-mminq3; 00124 double pass1=dmq/100.; 00125 double pass2=dmq2/100.; 00126 double pass3=dmq3/100.; 00127 double w1,w2,w3; 00128 00129 // cout<<"\n mass "<<mminq<<" "<<mmaxq<<" "<<pr2->mminq<<" "<<pr2->mmaxq<<" "<<pr3->mminq<<" "<<pr3->mmaxq; 00130 int p_id[6]= {up->Get_Idm(),p1->Get_Idm(),p2->Get_Idm(),p3->Get_Idm(),p4->Get_Idm(),p5->Get_Idm()}; 00131 if ( Fasemx!=0 ) 00132 { 00133 double wmx=Fasemx->Get_Wmax ( p_id); 00134 if ( wmx>0. ) 00135 { 00136 wfase_max=wmx; 00137 // Gout<<"\n Ev "<<Eventi_Fatti<<" F.S. trovato "<<up->Get_Nome() <<" "<<nome<<"=> "<<p1->Get_Nome() <<" "<<p2->Get_Nome() <<" 00138 // "<<p3->Get_Nome()<<" "<<p4->Get_Nome()<<" "<<p5->Get_Nome(); 00139 return; 00140 } 00141 } 00142 00143 00144 // define the wfase_max for p.s. normalizazion 00145 wfase_max=0.0; 00146 for ( mresq=mminq; mresq<mmaxq; mresq+=pass1 ) 00147 { 00148 mres=sqrt ( mresq ); 00149 w1=Fspace ( mpadre,mres,mfrate ); 00150 for ( mresq2=mminq2; mresq2<mmaxq2; mresq2+=pass2 ) 00151 { 00152 mres2=sqrt ( mresq2 ); 00153 w2=Fspace ( mres,mres2,mchd2 ); 00154 for ( mresq3=mminq3; mresq3<mmaxq3; mresq3+=pass3 ) 00155 { 00156 mres3=sqrt ( mresq3 ); 00157 w3=Fspace ( mres2,mres3,mchd3 ) *Fspace ( mres3,mchd4,mchd5 ); 00158 w=w1*w2*w3; 00159 if ( w>wfase_max ) wfase_max=w; 00160 } 00161 } 00162 } 00163 00164 Gout<<"\n Phase-Space:w="<<wfase_max<<" done for " <<up->Get_Nome()<<" ==> "<<p1->Get_Nome()<<" "<<p2->Get_Nome() 00165 <<" "<<p3->Get_Nome()<<" "<<p4->Get_Nome()<<" "<<p5->Get_Nome(); 00166 00167 new Fase3Max ( p_id,wfase_max ); 00168 } 00169 00170 void Partres5::Get_PhaseSpace() 00171 { 00172 iterfs=0; w=0.5; 00173 while ( iterfs++<300 ) 00174 { 00175 // resonance mass selection... 00176 mres =sqrt ( ( mresq =mminq+Pran() *dmq ) ); 00177 mres2 =sqrt ( ( mresq2=mminq2+Pran() *dmq2 ) ); 00178 if ( mres2+mchd2>mres ) continue; 00179 mres3 =sqrt ( ( mresq3=mminq3+Pran() *dmq3 ) ); 00180 if ( mres3+mchd3>mres2 ) continue; 00181 00182 w = Fspace ( mpadre, mres,mfrate ) *Fspace ( mres, mres2,mchd2 ); 00183 w*= Fspace ( mres2,mres3,mchd3 ) *Fspace ( mres3,mchd4,mchd5 ); 00184 00185 00186 if ( w >wfase_max ) 00187 { 00188 Gout<<"\n Wfas off: Ev "<<evento_.Gen.Event<<" iter="<<iterfs<<" prb="<<w 00189 <<" prbmax"<<wfase_max<<" "<<SelRea->Get_Titolo(); 00190 wfase_max=w+0.000001; 00191 // aggiorna wfase_max memorizzato nella lista 00192 int p_id[6]= {up->Get_Idm(),p1->Get_Idm(),p2->Get_Idm(),p3->Get_Idm(), p4->Get_Idm(),p5->Get_Idm()}; 00193 Fasemx->Set_Wmax ( p_id, wfase_max ); 00194 } 00195 00196 if ( Pran() <w/wfase_max ) break; 00197 00198 } 00199 if ( iterfs>400 ) 00200 Gout<<"\nPartRes5 iter="<<iterfs<<" Risonaze: m2345="<<mres<<" m345="<<mres2<<" m45="<<mres3; 00201 }