TCrossPart Class Reference

#include <TCrossPart.h>

Inheritance diagram for TCrossPart:


Public Member Functions
	TCrossPart (double e, double de, double nth0)
virtual	~TCrossPart ()
double	GetValue (const unsigned int npart)
void	MakeParts (double err)
size_t	GenUnWeightedEvent ()
size_t	GenWeightedEvent (double &w)
TLorentzVector **	GetParticles ()
void	SetNRad (size_t x=25000)
size_t	GetNRad ()
size_t	GetNfinal ()
int	GetPid (size_t i)
void	BeamSpread ()
void	SetBeamSpread (double x=1)
Protected Member Functions
double	radB (const TLorentzVector &k)
double	rad (double zp, double zm)
void	SetFinalParticles (size_t, const double , const int )
double	Generate ()
virtual void	SetJ ()=0
virtual bool	Accepted ()=0
void	J3PseudoScalars ()
Protected Attributes
TGenPhaseSpace	fevent
rb::TPhoton *	fphot
rb::TPhotonD *	fsphot
rb::TDFun *	fD
double	fe
double	fs
double	flogs
double	fbeta2
double	fdelta
double	fK
double	fep
double	fem
double	fsigmap
double	fsigmam
double	fppz
double	fpmz
bool	fBeamSpread
double	fsum [NPARTMAX]
double	fsum2 [NPARTMAX]
double	fmin [NPARTMAX]
double	fmax [NPARTMAX]
size_t	fN [NPARTMAX]
double	fcs [NPARTMAX]
double	fecs [NPARTMAX]
double	facs [NPARTMAX]
double	fxm
double	fxp
bool	inc [NPARTMAX]
size_t	fNRad
TLorentzVector	fpp
TLorentzVector	fpm
TLorentzVector	fW
TLorentzVector	fq
TLorentzVectorC	fJc
size_t	fNfinal
TLorentzVector **	fres
double	fm [7]
int	fpid [7+2]
double	fq2

Detailed Description

Definition at line 11 of file TCrossPart.h.

Constructor & Destructor Documentation

TCrossPart::TCrossPart	(	double	e,
		double	de,
		double	nth0
	)

Definition at line 53 of file TCrossPart.C.

References cos(), fBeamSpread, fbeta2, fcs, fD, fdelta, fe, fecs, fem, fep, fK, flogs, fmax, fmin, fN, fNRad, fphot, fpm, fpmz, fpp, fppz, fs, fsphot, fsum, fsum2, fW, rb::TDFun::GetBeta2(), genRecEmupikp::i, inc, rb::TPhotonD::Init(), M_PI, me, NPARTMAX, rb::TPhoton::SetThetaRange(), and theta0.

00053                                                       {
00054   fe = e;
00055   fbeta2 = 1 - pow(me/fe,2);
00056   double beta = sqrt(fbeta2);
00057   fppz =  fe*beta;
00058   fpmz = -fe*beta;
00059   double theta0 = nth0*sqrt(me/fe);
00060   //    flogs = log(theta0*theta0/4);
00061   long double c = beta*cos(theta0);
00062   flogs = -log((1l+c)/(1l-c));
00063   
00064   fdelta = de/fe;
00065   fK = M_PI*M_PI/3 - 1./2.;
00066   //  fphot = new TPhoton(e,de,420);
00067   fphot = new TPhoton(e);
00068   fphot->SetThetaRange(theta0, M_PI - theta0);
00069   fD = new TDFun(fe);
00070   fsphot = new TPhotonD();
00071   fsphot->Init(fdelta,2*fD->GetBeta2());
00072   for(size_t i=0;i<NPARTMAX;i++){
00073     fN[i] = 0;
00074     fsum[i] = 0;
00075     fsum2[i] = 0;
00076     fmin[i] = 1e100;
00077     fmax[i] = -1e100;
00078     fcs[i] = 0;
00079     fecs[i] = 0;
00080   }
00081   inc[0] = true;
00082   inc[1] = true;
00083   inc[2] = true;
00084   inc[3] = true;
00085   inc[4] = true;
00086 
00087   fNRad = 25000;
00088 
00089   fep = fe;
00090   fem = fe;
00091   fpp.SetPxPyPzE(0,0,fppz,fep);
00092   fpm.SetPxPyPzE(0,0,fpmz,fem);
00093   fW.SetPxPyPzE(0,0,0,2*fe);
00094   fs = fW.M2();
00095 
00096   fBeamSpread = false;
00097 }

TCrossPart::~TCrossPart ( ) [virtual]

Definition at line 99 of file TCrossPart.C.

References fD, fNfinal, fphot, fres, and fsphot.

00099                        {
00100   delete fsphot;
00101   delete fD;
00102   delete fphot;
00103   delete fres[fNfinal+1];
00104   delete [] fres;
00105 }

Member Function Documentation

virtual bool TCrossPart::Accepted ( ) [protected, pure virtual]

Implemented in T2piCrossPart, T3piCrossPart, T3piEtaCrossPart, T4piCrossPart, T5piCrossPart, TKKEtaCrossPart, TKKPCrossPart, and TRhoEtaCrossPart.

Referenced by rad(), and radB().

void TCrossPart::BeamSpread ( )

Definition at line 107 of file TCrossPart.C.

References fe, fem, fep, fpm, fpmz, fpp, fppz, fs, fsigmam, fsigmap, fW, and me.

Referenced by GetValue().

00107                            {
00108   fep = fe + gRandom->Gaus(0,fsigmap);
00109   fem = fe + gRandom->Gaus(0,fsigmam);
00110   double betap2 = 1 - pow(me/fep,2);
00111   double betap = sqrt(betap2);
00112   double betam2 = 1 - pow(me/fem,2);
00113   double betam = sqrt(betam2);
00114   fppz =  fep*betap;
00115   fpmz = -fem*betam;
00116   fpp.SetPxPyPzE(0,0,fppz,fep);
00117   fpm.SetPxPyPzE(0,0,fpmz,fem);
00118   fW.SetPxPyPzE(0,0,0,fep+fem);
00119   fs = fW.M2();
00120 }

double TCrossPart::Generate ( ) [protected]

Definition at line 32 of file TCrossPart.C.

References fevent, fm, fNfinal, and fq.

Referenced by rad(), and radB().

00032                            {
00033   // Attention!!! Before invoking this method you must set fq !
00034   if(!fevent.SetDecay(fq, fNfinal, fm, "Fermi")) return 0;
00035   /* 
00036      for N body decay dGamma has dimension of [MeV^(2*(N-2))] from Fermi
00037      dependence. Meaning dGamma see paper by Kuraev et al
00038   */
00039   const double i2piN[] = 
00040     {1,
00041      2.533029591058444286096986e-2,  // = 1/pow(2*M_PI,2) 2 final particles
00042      1.021176138454182960196328e-4,  // = 1/pow(2*M_PI,5) 3 final particles
00043      4.116812173964596341083153e-7,  // = 1/pow(2*M_PI,8) 4 final particles
00044      1.659668869795425306504901e-9,  // = 1/pow(2*M_PI,11) 5 final particles
00045      6.690858462739554865983570e-12, // = 1/pow(2*M_PI,14) 6 final particles
00046      2.697380651233861939614714e-14  // = 1/pow(2*M_PI,17) 7 final particles
00047      };
00048 
00049   double dGamma = fevent.Generate()*i2piN[fNfinal-1];
00050   return dGamma;
00051 }

size_t TCrossPart::GenUnWeightedEvent ( )

Definition at line 253 of file TCrossPart.C.

References facs, fem, fep, fmax, fNfinal, fres, fxm, fxp, and GetValue().

Referenced by Mcgpj::execute().

00253                                      {
00254   double r = gRandom->Rndm();
00255   size_t ipart = 0;
00256   while(r>facs[ipart]) ipart++;
00257   while (true){
00258     double m = GetValue(ipart);
00259     if ( m != 0  &&  m>fmax[ipart]*gRandom->Rndm() ) break;
00260   }
00261 
00262   if(ipart==4){
00263     fres[fNfinal+1]->SetPxPyPzE(0,0,0,0);
00264   } else {
00265     fres[fNfinal]->SetPxPyPzE(0,0,-fxm*fem,fxm*fem);
00266     fres[fNfinal+1]->SetPxPyPzE(0,0, fxp*fep,fxp*fep);
00267   }
00268   return ipart;
00269 }

size_t TCrossPart::GenWeightedEvent ( double & w )

Definition at line 271 of file TCrossPart.C.

References facs, fcs, fem, fep, fNfinal, fres, fxm, fxp, and GetValue().

00271                                             {
00272   double m = 0;
00273   size_t ipart = 0;
00274   do {
00275     double r = gRandom->Rndm();
00276     ipart = 0;
00277     while(r>facs[ipart]) ipart++;
00278     m = GetValue(ipart);
00279     if ( m > 0 ) break;
00280   } while (m<=0);
00281 
00282   if(ipart==4){
00283     fres[fNfinal+1]->SetPxPyPzE(0,0,0,0);
00284   } else {
00285     fres[fNfinal]->SetPxPyPzE(0,0,-fxm*fem,fxm*fem);
00286     fres[fNfinal+1]->SetPxPyPzE(0,0, fxp*fep,fxp*fep);
00287   }
00288 
00289   w = m/fcs[ipart];
00290   return ipart;
00291 }

size_t TCrossPart::GetNfinal ( ) [inline]

Definition at line 70 of file TCrossPart.h.

References fNfinal.

Referenced by Mcgpj::execute().

00070 {return fNfinal;}

size_t TCrossPart::GetNRad ( ) [inline]

Definition at line 69 of file TCrossPart.h.

References fNRad.

00069 {return fNRad;}

TLorentzVector** TCrossPart::GetParticles ( ) [inline]

Definition at line 67 of file TCrossPart.h.

References fres.

Referenced by Mcgpj::execute().

00067 {return fres;}

int TCrossPart::GetPid ( size_t i ) [inline]

Definition at line 71 of file TCrossPart.h.

References fpid.

Referenced by Mcgpj::execute().

00071 {return fpid[i];}

double TCrossPart::GetValue ( const unsigned int npart )

Definition at line 122 of file TCrossPart.C.

References alpha, BeamSpread(), rb::TDFun::Eval(), rb::TDFun::EvalSoft(), fBeamSpread, fD, fe, fK, flogs, fphot, fsphot, fxm, fxp, rb::TPhotonD::GetEnergy(), rb::TPhoton::GetEnergy(), rb::TPhotonD::GetENorm(), rb::TPhoton::GetENorm(), rb::TPhoton::GetNewPhoton(), rb::TPhoton::GetPhotNorm(), inc, M_PI, rad(), and radB().

Referenced by GenUnWeightedEvent(), GenWeightedEvent(), and MakeParts().

00122                                                    {
00123   double m2 = 0;
00124   const double AP = alpha/M_PI;
00125   if(!inc[npart]) return 0;
00126   if(fBeamSpread) BeamSpread();
00127   switch ( npart ){
00128   case 4:
00129     {
00130       const TLorentzVector &k = fphot->GetNewPhoton();
00131       double d3k_om = k.E()*fphot->GetPhotNorm()*fe;
00132       m2 = radB(k)*d3k_om;
00133     }
00134     break;
00135   case 1:
00136     {
00137       fxp = fphot->GetEnergy();
00138       double normp = fphot->GetENorm();
00139       fxm = fsphot->GetEnergy();
00140       double normm = fsphot->GetENorm();
00141       
00142       double comp = ((1-fxp+0.5*fxp*fxp)*flogs+0.5*fxp*fxp)/fxp;
00143       //        double comp = 0;
00144       m2 = rad(1-fxp, 1-fxm)*
00145         ((1+AP*fK)*
00146          fD->Eval(1-fxp)*
00147          fD->EvalSoft(fxm)*normm
00148          + AP*comp)*normp;
00149     }
00150     break;
00151   case 2:
00152     {
00153       fxm = fphot->GetEnergy();
00154       double normm = fphot->GetENorm();
00155       fxp = fsphot->GetEnergy();
00156       double normp = fsphot->GetENorm();
00157       
00158       double comp = ((1-fxm+0.5*fxm*fxm)*flogs+0.5*fxm*fxm)/fxm;
00159       //        double comp = 0;
00160       m2 = rad(1-fxp, 1-fxm)*
00161         ((1+AP*fK)*
00162          fD->Eval(1-fxm)*
00163          fD->EvalSoft(fxp)*normp
00164          + AP*comp)*normm;
00165     }
00166     break;
00167   case 0:
00168     {
00169       //      return 0;
00170       fxm = fsphot->GetEnergy();
00171       double normm = fsphot->GetENorm();
00172       fxp = fsphot->GetEnergy();
00173       double normp = fsphot->GetENorm();
00174       
00175       double cs = rad(1-fxp, 1-fxm);
00176       m2 = cs*((1+AP*fK)*
00177                fD->EvalSoft(fxm)*normm*
00178                fD->EvalSoft(fxp)*normp);
00179     }
00180     break;
00181   case 3:
00182     {
00183       fxm = fphot->GetEnergy();
00184       double normm = fphot->GetENorm();
00185       fxp = fphot->GetEnergy();
00186       double normp = fphot->GetENorm();
00187       
00188       double cs = rad(1-fxp, 1-fxm);
00189       m2 = cs*(1+AP*fK)*
00190         fD->Eval(1-fxm)*normm*
00191         fD->Eval(1-fxp)*normp;
00192     }
00193     break;
00194   default:
00195     break;
00196   }
00197   return m2;
00198 }

void TCrossPart::J3PseudoScalars ( ) [protected]

Definition at line 340 of file TCrossPart.C.

References fJc, fres, and TLorentzVectorC::SetPxPyPzE().

Referenced by TKKPCrossPart::SetJ(), TKKEtaCrossPart::SetJ(), and T3piCrossPart::SetJ().

00340                                 {
00341   // The hadronic current 
00342   // here for 3pi J_\mu=\epsilon_{a b g \mu}q+_{a}q-_{b}q0_{g}
00343 
00344   TLorentzVector &qp = *fres[0];  //p1 momentum
00345   TLorentzVector &qm = *fres[1];  //p2 momentum
00346   TLorentzVector &q0 = *fres[2];  //p3 momentum
00347   
00348   double Jx = qp.Y()*qm.Z()*q0.T()
00349     -    qp.Y()*qm.T()*q0.Z()
00350     +    qp.T()*qm.Y()*q0.Z()
00351     -    qp.T()*qm.Z()*q0.Y()
00352     +    qp.Z()*qm.T()*q0.Y()
00353     -    qp.Z()*qm.Y()*q0.T();
00354   
00355   double Jy = qp.X()*qm.Z()*q0.T()
00356     -    qp.X()*qm.T()*q0.Z()
00357     +    qp.T()*qm.X()*q0.Z()
00358     -    qp.T()*qm.Z()*q0.X()
00359     +    qp.Z()*qm.T()*q0.X()
00360     -    qp.Z()*qm.X()*q0.T();
00361 
00362   double Jz = qp.Y()*qm.X()*q0.T()
00363     -    qp.Y()*qm.T()*q0.X()
00364     +    qp.T()*qm.Y()*q0.X()
00365     -    qp.T()*qm.X()*q0.Y()
00366     +    qp.X()*qm.T()*q0.Y()
00367     -    qp.X()*qm.Y()*q0.T();
00368 
00369   double Jt = qp.Y()*qm.Z()*q0.X()
00370     -    qp.Y()*qm.X()*q0.Z()
00371     +    qp.X()*qm.Y()*q0.Z()
00372     -    qp.X()*qm.Z()*q0.Y()
00373     +    qp.Z()*qm.X()*q0.Y()
00374     -    qp.Z()*qm.Y()*q0.X();
00375 
00376   fJc.SetPxPyPzE(Jx,Jy,Jz,Jt);
00377 }

void TCrossPart::MakeParts ( double err )

Definition at line 200 of file TCrossPart.C.

References facs, fcs, fecs, fmax, fmin, fN, fNRad, fsum, fsum2, GetValue(), genRecEmupikp::i, and NPARTMAX.

Referenced by Mcgpj::initialize().

00200                                     {
00201   double dcs = pow(err,2)/NPARTMAX;
00202   for(size_t npart = 0; npart<NPARTMAX; npart++){
00203     //      if(inc[npart]) continue;
00204     double err2;
00205     do {
00206       for(size_t i=0;i<fNRad;i++){
00207         double m2 = GetValue(npart);
00208         fsum[npart] += m2; 
00209         fsum2[npart] += m2*m2; 
00210         fN[npart]++;
00211         if(m2>fmax[npart]) fmax[npart] = m2;
00212         if(m2<fmin[npart]) fmin[npart] = m2;
00213       }
00214       double iN = 1./fN[npart];
00215       err2 = (fsum2[npart]*iN - pow(fsum[npart]*iN,2))*iN;
00216       std::cout<<npart<<" "<<fsum[npart]/fN[npart]<<" "<<sqrt(err2)<<std::endl;
00217     } while(err2 > dcs);
00218   }
00219   
00220   for(size_t i=0;i<NPARTMAX;i++){
00221     if(fN[i]>0) fcs[i] =  fsum[i]/fN[i];
00222     if(fN[i]>0) fecs[i] = (fsum2[i]/fN[i] - pow(fcs[i],2))/fN[i];
00223   }
00224   
00225   double cs_sum = 0, ecs_sum = 0;
00226   for(size_t i=0;i<NPARTMAX;i++){
00227     cs_sum += fcs[i];
00228     ecs_sum += fecs[i];
00229     std::cout<<"Part = "<<i<<" "
00230         <<"CS = "<<fcs[i]<<" +- "<<sqrt(fecs[i])<<" nb "
00231         <<"Min = "<<fmin[i]<<" "
00232         <<"Max = "<<fmax[i]<<" "
00233         <<"N = "<<fN[i]<<" "
00234         <<std::endl;
00235   }
00236   
00237   facs[0] = fcs[0];
00238   for(size_t i=1;i<NPARTMAX;i++){
00239     facs[i] = facs[i-1] + fcs[i];
00240   }
00241   for(size_t i=0;i<NPARTMAX;i++){
00242     facs[i] /= cs_sum;
00243     std::cout<<facs[i]<<" ";
00244   }
00245   std::cout<<std::endl;
00246   
00247   ecs_sum = sqrt(ecs_sum);
00248   std::cout<<"sigma = "<<cs_sum<<" +- "<<ecs_sum<<" nb, "
00249            <<"relative error is "<<ecs_sum/cs_sum*100<<" %"
00250            <<std::endl;
00251 }

double TCrossPart::rad	(	double	zp,
		double	zm
	)			`[protected]`

Definition at line 319 of file TCrossPart.C.

References Accepted(), alpha, conj(), fem, fep, fJc, fpm, fpmz, fpp, fppz, fq, fq2, Generate(), M_PI, rb::MeV2nb, SetJ(), and T.

Referenced by GetValue().

00319                                           {
00320   fq.SetPxPyPzE(0,0,zp*fppz+zm*fpmz,zp*fep+zm*fem); // photon virtuality 
00321   fq2 = fq*fq;//fq.M2(); // s' after photon emission
00322   //  TLorentzVector q = zp*fpp + zm*fpm;
00323   double dGamma = Generate();
00324   if(!Accepted()) return 0;             // cut which can be removed
00325   
00326   SetJ(); // Set hadronic current J
00327   
00328   std::complex<double> ppJ = fJc*fpp; ppJ *= zp;
00329   std::complex<double> pmJ = fJc*fpm; pmJ *= zm;
00330   double J2 = std::real(fJc*fJc);
00331 
00332   double T = 2*std::real(ppJ*std::conj(pmJ)) - 0.5*fq2*J2; // folded leptonic and hadronic tensors
00333 
00334   const double AP28 = 8*alpha*M_PI*alpha*M_PI;
00335   double cs = AP28/(fq2*fq2*fq2)*T*dGamma*MeV2nb;
00336 
00337   return cs;
00338 }

double TCrossPart::radB ( const TLorentzVector & k ) [protected]

Definition at line 294 of file TCrossPart.C.

References Accepted(), alpha, conj(), fem, fep, fJc, fpm, fpp, fq, fq2, fs, Generate(), rb::MeV2nb, SetJ(), and T.

Referenced by GetValue().

00294                                               { // photon momentum
00295   fq.SetPxPyPzE(-k.Px(),-k.Py(),-k.Pz(),fep+fem-k.E()); // photon virtuality
00296   fq2 = fq*fq;//fq.M2();
00297   double dGamma = Generate();
00298   if(!Accepted()) return 0;              // cut which can be removed
00299 
00300   SetJ(); // Set hadronic current J
00301 
00302   std::complex<double> ppJ = fJc*fpp;
00303   std::complex<double> pmJ = fJc*fpm;
00304   double J2 = std::real(fJc*fJc);
00305 
00306   double chip = k*fpp; //chi1 in Kuraev's paper
00307   double chim = k*fpm; //chi2 in Kuraev's paper
00308 
00309   double T = -0.5*(fs*fq2 + 2*chip*chip + 2*chim*chim)*J2 -
00310     fq2*(std::real(pmJ*std::conj(pmJ))+std::real(ppJ*std::conj(ppJ))); // folded bremsstrahlung and hadronic tensors
00311 
00312   const double alpha3 = alpha*alpha*alpha;
00313   double cs = 2*alpha3/(fs*fq2*fq2*chim*chip)*T*dGamma*MeV2nb;
00314   return cs;
00315 }

void TCrossPart::SetBeamSpread ( double x = 1 ) [inline]

Definition at line 73 of file TCrossPart.h.

References fBeamSpread, fsigmam, fsigmap, and x.

Referenced by Mcgpj::initialize().

00073                                 { // beam spread in MeV
00074     fBeamSpread = true;
00075     fsigmap = x;
00076     fsigmam = x;
00077   }

void TCrossPart::SetFinalParticles	(	size_t	,
		const double *	,
		const int *
	)			`[protected]`

Definition at line 6 of file TCrossPart.C.

References fevent, fm, fNfinal, fphot, fpid, fres, fW, rb::TPhoton::GetPhoton(), and genRecEmupikp::i.

Referenced by T2piCrossPart::T2piCrossPart(), T3piCrossPart::T3piCrossPart(), T3piEtaCrossPart::T3piEtaCrossPart(), T4piCrossPart::T4piCrossPart(), T5piCrossPart::T5piCrossPart(), TKKEtaCrossPart::TKKEtaCrossPart(), TKKPCrossPart::TKKPCrossPart(), and TRhoEtaCrossPart::TRhoEtaCrossPart().

00006                                                                            {
00007   // set number and masses of final particles
00008   if(n>7){
00009     std::cout<<"Too many particles"<<std::endl;
00010     exit(1);
00011   }
00012   fNfinal = n;
00013   fres = new TLorentzVector*[fNfinal+2];
00014   if(fevent.SetDecay(fW, fNfinal, (double*)m, "Fermi")){
00015     for(size_t i=0;i<fNfinal;i++){
00016       fm[i] = m[i];
00017       fpid[i] = pid[i];
00018       fres[i] = fevent.GetDecay(i);
00019     }
00020   } else {
00021     std::cout<<"Energy is not enough to generate process!"<<std::endl;
00022     exit(1);
00023   }
00024   // include to the list of final particles two hard photons
00025   // (collinear and not collinear)
00026   fres[fNfinal] = (TLorentzVector*)&fphot->GetPhoton();
00027   fpid[fNfinal] = 22;
00028   fres[fNfinal+1] = new TLorentzVector(0,0,0,0);
00029   fpid[fNfinal+1] = 22;
00030 }