BOSS_7.0.2: EvtGenKine Class Reference

double EvtGenKine::PhaseSpace	(	int	ndaug,
		double	mass[30],
		EvtVector4R	p4[30],
		double	mp
	)			`[static]`

Definition at line 50 of file EvtGenKine.cc.

References alpha, EvtVector4R::applyRotateEuler(), cos(), energy, calibUtil::ERROR, EvtPawt(), EvtRandom::Flat(), genRecEmupikp::i, report(), eformat::write::set(), EvtVector4R::set(), sin(), subSeperate::temp, and EvtConst::twoPi.

Referenced by EvtBtoXsll::decay(), EvtBtoXsgamma::decay(), EvtBtoXsEtap::decay(), and EvtParticle::initializePhaseSpace().

00057 {
00058 
00059   double energy, p3, alpha, beta;
00060 
00061   if ( ndaug == 1 ) {
00062      p4[0].set(mass[0],0.0,0.0,0.0);
00063      return 1.0;
00064   }
00065 
00066 
00067   if ( ndaug == 2 ) {
00068 
00069 //  Two body phase space
00070 
00071      energy = ( mp*mp + mass[0]*mass[0] -
00072               mass[1]*mass[1] ) / ( 2.0 * mp );
00073 
00074      p3 = sqrt( energy*energy - mass[0]*mass[0] );
00075 
00076      p4[0].set( energy, 0.0, 0.0, p3 );
00077 
00078      energy = mp - energy;
00079      p3 = -1.0*p3;
00080      p4[1].set( energy, 0.0, 0.0, p3 );
00081 
00082 //   Now rotate four vectors.
00083 
00084      alpha = EvtRandom::Flat( EvtConst::twoPi );
00085      beta = acos(EvtRandom::Flat( -1.0, 1.0 ));
00086    
00087      p4[0].applyRotateEuler( alpha, beta, -alpha );
00088      p4[1].applyRotateEuler( alpha, beta, -alpha );
00089 
00090      return 1.0;
00091   }
00092 
00093   if ( ndaug != 2 ) {
00094 
00095     double wtmax=0.0;
00096     double pm[5][30],to[4],pmin,pmax,psum,rnd[30];
00097     double ran,wt,pa,costh,sinth,phi,p[4][30],be[4],bep,temp;
00098     int i,il,ilr,i1,il1u,il1,il2r,ilu;
00099     int il2=0;
00100     
00101      for(i=0;i<ndaug;i++){
00102        pm[4][i]=0.0;
00103        rnd[i]=0.0;
00104      }     
00105     
00106      pm[0][0]=mp;
00107      pm[1][0]=0.0;
00108      pm[2][0]=0.0;
00109      pm[3][0]=0.0;
00110      pm[4][0]=mp;
00111 
00112      to[0]=mp;
00113      to[1]=0.0;
00114      to[2]=0.0;
00115      to[3]=0.0;
00116 
00117      psum=0.0;
00118      for(i=1;i<ndaug+1;i++){
00119        psum=psum+mass[i-1];
00120      }
00121 
00122      pm[4][ndaug-1]=mass[ndaug-1];
00123 
00124      switch (ndaug) {
00125 
00126      case 1:
00127        wtmax=1.0/16.0;
00128        break;
00129      case 2:
00130        wtmax=1.0/150.0;
00131        break;
00132      case 3:
00133        wtmax=1.0/2.0;
00134        break;
00135      case 4:
00136        wtmax=1.0/5.0;
00137        break;
00138      case 5:
00139        wtmax=1.0/15.0;
00140        break;
00141      case 6:
00142        wtmax=1.0/15.0;
00143        break;
00144      case 7:
00145        wtmax=1.0/15.0;
00146        break;
00147      case 8:
00148        wtmax=1.0/15.0;
00149        break;
00150      case 9:
00151        wtmax=1.0/15.0;
00152        break;
00153      case 10:
00154        wtmax=1.0/15.0;
00155        break;
00156      case 11:
00157        wtmax=1.0/15.0;
00158        break;
00159      case 12:
00160        wtmax=1.0/15.0;
00161        break;
00162      case 13:
00163        wtmax=1.0/15.0;
00164        break;
00165      case 14:
00166        wtmax=1.0/15.0;
00167        break;
00168      case 15:
00169        wtmax=1.0/15.0;
00170        break;
00171      default:
00172        report(ERROR,"EvtGen") << "too many daughters for phase space..." << ndaug << " "<< mp <<endl;;
00173        break;
00174      }
00175 
00176      pmax=mp-psum+mass[ndaug-1];
00177 
00178      pmin=0.0;
00179 
00180      for(ilr=2;ilr<ndaug+1;ilr++){
00181        il=ndaug+1-ilr;
00182        pmax=pmax+mass[il-1];
00183        pmin=pmin+mass[il+1-1];
00184        wtmax=wtmax*EvtPawt(pmax,pmin,mass[il-1]);
00185      }
00186 
00187      //report(INFO,"EvtGen") << "wtmax:"<<wtmax<<endl;
00188 
00189      do{
00190 
00191        rnd[0]=1.0;
00192        il1u=ndaug-1;
00193 
00194        for (il1=2;il1<il1u+1;il1++){
00195          ran=EvtRandom::Flat();
00196          for (il2r=2;il2r<il1+1;il2r++){
00197            il2=il1+1-il2r;
00198            if (ran<=rnd[il2-1]) goto two39;
00199            rnd[il2+1-1]=rnd[il2-1];
00200          }
00201        two39:
00202          rnd[il2+1-1]=ran;
00203        }
00204 
00205        rnd[ndaug-1]=0.0;
00206        wt=1.0;
00207        for(ilr=2;ilr<ndaug+1;ilr++){
00208          il=ndaug+1-ilr;
00209          pm[4][il-1]=pm[4][il+1-1]+mass[il-1]+
00210            (rnd[il-1]-rnd[il+1-1])*(mp-psum);
00211          wt=wt*EvtPawt(pm[4][il-1],pm[4][il+1-1],mass[il-1]);
00212        }
00213        if (wt>wtmax) {
00214          report(ERROR,"EvtGen") << "wtmax to small in EvtPhaseSpace with "
00215                                 << ndaug <<" daughters"<<endl;;
00216        }
00217      } while (wt<EvtRandom::Flat(wtmax));
00218      
00219      ilu=ndaug-1;
00220 
00221      for (il=1;il<ilu+1;il++){
00222        pa=EvtPawt(pm[4][il-1],pm[4][il+1-1],mass[il-1]);
00223        costh=EvtRandom::Flat(-1.0,1.0);
00224        sinth=sqrt(1.0-costh*costh);
00225        phi=EvtRandom::Flat(EvtConst::twoPi);
00226        p[1][il-1]=pa*sinth*cos(phi);
00227        p[2][il-1]=pa*sinth*sin(phi);
00228        p[3][il-1]=pa*costh;
00229        pm[1][il+1-1]=-p[1][il-1];
00230        pm[2][il+1-1]=-p[2][il-1];
00231        pm[3][il+1-1]=-p[3][il-1];
00232        p[0][il-1]=sqrt(pa*pa+mass[il-1]*mass[il-1]);
00233        pm[0][il+1-1]=sqrt(pa*pa+pm[4][il+1-1]*pm[4][il+1-1]);
00234      }
00235 
00236       p[0][ndaug-1]=pm[0][ndaug-1];
00237       p[1][ndaug-1]=pm[1][ndaug-1];
00238       p[2][ndaug-1]=pm[2][ndaug-1];
00239       p[3][ndaug-1]=pm[3][ndaug-1];
00240 
00241       for (ilr=2;ilr<ndaug+1;ilr++){
00242         il=ndaug+1-ilr;
00243         be[0]=pm[0][il-1]/pm[4][il-1];
00244         be[1]=pm[1][il-1]/pm[4][il-1];
00245         be[2]=pm[2][il-1]/pm[4][il-1];
00246         be[3]=pm[3][il-1]/pm[4][il-1];
00247 
00248         for (i1=il;i1<ndaug+1;i1++){
00249           bep=be[1]*p[1][i1-1]+be[2]*p[2][i1-1]+
00250               be[3]*p[3][i1-1]+be[0]*p[0][i1-1];
00251           temp=(p[0][i1-1]+bep)/(be[0]+1.0);
00252           p[1][i1-1]=p[1][i1-1]+temp*be[1];
00253           p[2][i1-1]=p[2][i1-1]+temp*be[2];
00254           p[3][i1-1]=p[3][i1-1]+temp*be[3];
00255           p[0][i1-1]=bep;
00256 
00257         }
00258       }
00259 
00260      for (ilr=0;ilr<ndaug;ilr++){
00261        p4[ilr].set(p[0][ilr],p[1][ilr],p[2][ilr],p[3][ilr]);
00262      }
00263 
00264      return 1.0;
00265   }
00266 
00267   return 1.0;
00268 
00269 }

double EvtGenKine::PhaseSpacePole	(	double	M,
		double	m1,
		double	m2,
		double	m3,
		double	a,
		EvtVector4R	p4[10]
	)			`[static]`

Definition at line 272 of file EvtGenKine.cc.

References alpha, EvtVector4R::applyRotateEuler(), EvtRandom::Flat(), EvtVector4R::set(), and EvtConst::twoPi.

Referenced by EvtParticle::initializePhaseSpace().

00277 {
00278 
00279   //f1   = 1  (phasespace)
00280   //f2   = a*(1/m12sq)^2
00281 
00282   double m12sqmax=(M-m3)*(M-m3);
00283   double m12sqmin=(m1+m2)*(m1+m2);
00284 
00285   double m13sqmax=(M-m2)*(M-m2);
00286   double m13sqmin=(m1+m3)*(m1+m3);
00287 
00288   double v1=(m12sqmax-m12sqmin)*(m13sqmax-m13sqmin);
00289   double v2= a*(1.0/m12sqmin-1.0/m12sqmax)*(m13sqmax-m13sqmin);
00290 
00291   double m12sq,m13sq;
00292 
00293   double r=v1/(v1+v2);
00294 
00295   //report(INFO,"EvtGen") << "v1,v2:"<<v1<<","<<v2<<endl;
00296 
00297   double m13min,m13max;
00298 
00299   do{
00300 
00301     m13sq=EvtRandom::Flat(m13sqmin,m13sqmax);
00302   
00303     if (r>EvtRandom::Flat()){
00304       m12sq=EvtRandom::Flat(m12sqmin,m12sqmax);
00305     }
00306     else{
00307       m12sq=1.0/(1.0/m12sqmin-EvtRandom::Flat()*(1.0/m12sqmin-1.0/m12sqmax));
00308     }
00309 
00310     //kinematically allowed?
00311     double E3star=(M*M-m12sq-m3*m3)/sqrt(4*m12sq);
00312     double E1star=(m12sq+m1*m1-m2*m2)/sqrt(4*m12sq);
00313     double p3star=sqrt(E3star*E3star-m3*m3);
00314     double p1star=sqrt(E1star*E1star-m1*m1);
00315     m13max=(E3star+E1star)*(E3star+E1star)-
00316       (p3star-p1star)*(p3star-p1star);
00317     m13min=(E3star+E1star)*(E3star+E1star)-
00318       (p3star+p1star)*(p3star+p1star);
00319 
00320   }while(m13sq<m13min||m13sq>m13max);
00321 
00322   
00323   double E2=(M*M+m2*m2-m13sq)/(2.0*M);
00324   double E3=(M*M+m3*m3-m12sq)/(2.0*M);
00325   double E1=M-E2-E3;
00326   double p1mom=sqrt(E1*E1-m1*m1);      
00327   double p3mom=sqrt(E3*E3-m3*m3);
00328   double cost13=(2.0*E1*E3+m1*m1+m3*m3-m13sq)/(2.0*p1mom*p3mom);
00329 
00330   //report(INFO,"EvtGen") << m13sq << endl;
00331   //report(INFO,"EvtGen") << m12sq << endl;
00332   //report(INFO,"EvtGen") << E1 << endl;
00333   //report(INFO,"EvtGen") << E2 << endl;
00334   //report(INFO,"EvtGen") << E3 << endl;
00335   //report(INFO,"EvtGen") << p1mom << endl;
00336   //report(INFO,"EvtGen") << p3mom << endl;
00337   //report(INFO,"EvtGen") << cost13 << endl;
00338   
00339 
00340   p4[2].set(E3,0.0,0.0,p3mom);
00341   p4[0].set(E1,p1mom*sqrt(1.0-cost13*cost13),0.0,p1mom*cost13);
00342   p4[1].set(E2,-p1mom*sqrt(1.0-cost13*cost13),0.0,-p1mom*cost13-p3mom);
00343 
00344   //report(INFO,"EvtGen") << "p4:"<<p4[0]<<p4[1]<<p4[2]<<endl;
00345 
00346   double alpha = EvtRandom::Flat( EvtConst::twoPi );
00347   double beta = acos(EvtRandom::Flat( -1.0, 1.0 ));
00348   double gamma = EvtRandom::Flat( EvtConst::twoPi );
00349   
00350   p4[0].applyRotateEuler( alpha, beta, gamma );
00351   p4[1].applyRotateEuler( alpha, beta, gamma );
00352   p4[2].applyRotateEuler( alpha, beta, gamma );
00353   
00354   return 1.0+a/(m12sq*m12sq);
00355 
00356 }


Static Public Member Functions
static double	PhaseSpace (int ndaug, double mass[30], EvtVector4R p4[30], double mp)
static double	PhaseSpacePole (double M, double m1, double m2, double m3, double a, EvtVector4R p4[10])

EvtGenKine Class Reference

Static Public Member Functions

Detailed Description

Member Function Documentation