1 ///////////////////////////////////////////////////////////////////////////
5 // This file is part of starlight.
7 // starlight is free software: you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation, either version 3 of the License, or
10 // (at your option) any later version.
12 // starlight is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
17 // You should have received a copy of the GNU General Public License
18 // along with starlight. If not, see <http://www.gnu.org/licenses/>.
20 ///////////////////////////////////////////////////////////////////////////
22 // File and Version Information:
23 // $Rev:: 184 $: revision of last commit
24 // $Author:: jnystrand $: author of last commit
25 // $Date:: 2014-09-12 00:59:43 +0200 #$: date of last commit
29 // Fixed bug which gave incorrect minv distribution in gammagammaleptonpair.
30 // Moved loop over W and Y from pickw to twoLeptonCrossSection in
31 // gammagammaleptonpair to speed up event generation.
32 // Changed to Woods-Saxon radius in twophotonluminosity to be consistent
36 ///////////////////////////////////////////////////////////////////////////
44 #include "starlightconstants.h"
45 #include "gammagammaleptonpair.h"
51 //_____________________________________________________________________________
52 Gammagammaleptonpair::Gammagammaleptonpair(beamBeamSystem& bbsystem)
53 : eventChannel(bbsystem)
55 //Initialize randomgenerator with our seed.
56 cout<<"Randy in leptonpair construction: "<<randyInstance.Rndom()<<endl;
57 //Storing inputparameters into protected members for use
58 _GGlepInputnumw=inputParametersInstance.nmbWBins();
59 _GGlepInputnumy=inputParametersInstance.nmbRapidityBins();
60 _GGlepInputpidtest=inputParametersInstance.prodParticleType();
61 _GGlepInputGamma_em=inputParametersInstance.beamLorentzGamma();
62 //Let us read in the luminosity tables
64 //Now we will calculate the crosssection
65 twoLeptonCrossSection();
66 //If it is a tauon, calculate its tables
67 if(inputParametersInstance.prodParticleId()==starlightConstants::TAUONDECAY) calculateTable();
71 //______________________________________________________________________________
72 Gammagammaleptonpair::~Gammagammaleptonpair()
76 //______________________________________________________________________________
77 void Gammagammaleptonpair::twoLeptonCrossSection()
79 //This function calculates section for 2-particle decay. For reference, see STAR Note 243, Eq. 9.
80 //calculate the 2-lepton differential cross section
81 //the 100 is to convert to barns
82 //the 2 is to account for the fact that we integrate only over one half of the rapidity range
83 //Multiply all _Farray[] by _f_max
85 for(int i=0;i<_GGlepInputnumw;i++)
87 for(int j=0;j<_GGlepInputnumy;j++)
89 // _sigmax[i][j]=2.*Gammagammaleptonpair::twoMuonCrossSection(_Warray[i])*_f_max*_Farray[i][j]/100.;
90 _sigmax[i][j]=twoMuonCrossSection(_Warray[i])*_f_max*_Farray[i][j]/(100.*_Warray[i]);
93 //calculate the total two-lepton cross section
95 for(int i=0;i<_GGlepInputnumw-1;i++)
97 for(int j=0;j<_GGlepInputnumy-1;j++)
99 // _sigmaSum = _sigmaSum +2.*((_sigmax[i][j]+_sigmax[i+1][j]+_sigmax[i][j+1]+_sigmax[i+1][j+1])/4.*(_Yarray[j+1]-_Yarray[j])*(_Warray[i+1]-_Warray[i])/((_Warray[i+1]+_Warray[i])/2.));
100 // _sigmaSum = _sigmaSum +((_sigmax[i][j]+_sigmax[i+1][j]+_sigmax[i][j+1]+_sigmax[i+1][j+1])/4.*(_Yarray[j+1]-_Yarray[j])*(_Warray[i+1]-_Warray[i])/((_Warray[i+1]+_Warray[i])/2.));
101 sigmasum = sigmasum +(_sigmax[i][j]+_sigmax[i+1][j]+_sigmax[i][j+1]+_sigmax[i+1][j+1])/4.*(_Yarray[j+1]-_Yarray[j])*(_Warray[i+1]-_Warray[i]);
104 cout << "The total "<<_GGlepInputpidtest<<" cross-section is: "<<sigmasum<<" barns."<<endl;
106 // Do this integration here, once per run rather than once per event (JN 220710)
107 //integrate sigma down to a function of just w
111 for(int i=0;i<_ReadInputnumw;i++)
114 for(int j=0;j<_ReadInputnumy-1;j++)
116 _sigofw[i] = _sigofw[i]+(_Yarray[j+1]-_Yarray[j])*(_sigmax[i][j+1]+_sigmax[i][j])/2.;
120 //calculate the unnormalized sgfint array
122 for(int i=0;i<_ReadInputnumw-1;i++)
124 sgf=(_sigofw[i+1]+_sigofw[i])*(_Warray[i+1]-_Warray[i])/2.;
125 _sigfint[i+1]=_sigfint[i]+sgf;
128 //normalize sgfint array
129 _signormw=_sigfint[_ReadInputnumw-1];
130 for(int i=0;i<_ReadInputnumw;i++)
132 _sigfint[i]=_sigfint[i]/_signormw;
138 //______________________________________________________________________________
139 double Gammagammaleptonpair::twoMuonCrossSection(double w)
141 //This function gives the two muon cross section as a function of Y&W.
142 //Using the formula given in G.Soff et. al Nuclear Equation of State, part B, 579
143 double s=0.,Etest=0.,deltat=0.,xL=0.,sigmuix=0.,alphasquared=0.,hbarcsquared=0.;
145 Etest = 4.*getMass()*getMass()/s;
146 deltat = s * sqrt(1.-Etest);
147 xL = 2.*log(sqrt(s)/(2.*getMass())+sqrt(1./Etest-1.));
148 alphasquared = starlightConstants::alpha*starlightConstants::alpha;
149 hbarcsquared = starlightConstants::hbarc*starlightConstants::hbarc;
150 sigmuix = 4.*starlightConstants::pi*alphasquared/s*hbarcsquared*((1+Etest-0.5*Etest*Etest)*xL-(1./s+Etest/s)*deltat);
157 //______________________________________________________________________________
158 void Gammagammaleptonpair::pickw(double &w)
160 // This function picks a w for the 2-photon calculation.
162 double x=0.,remainarea=0.,remainw=0.,a=0.,b=0.,c=0.;
172 //deal with the case where sigma is an array
173 //_sigofw is simga integrated over y using a linear interpolation
174 //sigint is the integral of sgfint, normalized
176 //pick a random number
177 x = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
178 //compare x and sgfint to find the ivalue which is just less than the random number x
179 for(int i=0;i<_GGlepInputnumw;i++)
181 if(x > _sigfint[i]) ivalw=i;
183 //remainder above ivalw
184 remainarea = x - _sigfint[ivalw];
186 //figure out what point corresponds to the excess area in remainarea
187 c = -remainarea*_signormw/(_Warray[ivalw+1]-_Warray[ivalw]);
189 a = (_sigofw[ivalw+1]-_sigofw[ivalw])/2.;
195 remainw = (-b+sqrt(b*b-4.*a*c))/(2.*a);
199 //calculate the w value
200 w = _Warray[ivalw]+(_Warray[ivalw+1]-_Warray[ivalw])*remainw;
205 //______________________________________________________________________________
206 void Gammagammaleptonpair::picky(double &y)
208 // This function picks a y given a W
212 sigofy = new double[starlightLimits::MAXYBINS];
213 sgfint = new double[starlightLimits::MAXWBINS];
215 double remainw =0.,remainarea=0.,remainy=0.,a=0.,b=0.,c=0.,sgf=0.,signorm=0.,x=0.;
220 //average over two colums to get y array
221 for(int j=0;j<_GGlepInputnumy;j++)
223 sigofy[j]=_sigmax[ivalw][j]+(_sigmax[ivalw+1][j]-_sigmax[ivalw][j])*remainw;
226 //calculate the unnormalized sgfint
228 for(int j=0;j<_GGlepInputnumy-1;j++)
230 sgf = (sigofy[j+1]+sigofy[j])/2.;
231 sgfint[j+1]=sgfint[j]+sgf*(_Yarray[j+1]-_Yarray[j]);
234 //normalize the sgfint array
235 signorm = sgfint[_GGlepInputnumy-1];
237 for(int j=0;j<_GGlepInputnumy;j++)
239 sgfint[j]=sgfint[j]/signorm;
242 //pick a random number
243 x = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
244 //compare x and sgfint to find the ivalue which is just less then the random number x
245 for(int i=0;i<_GGlepInputnumy;i++)
247 if(x > sgfint[i]) ivaly = i;
249 //remainder above ivaly
250 remainarea = x - sgfint[ivaly];
252 //figure what point corresponds to the leftover area in remainarea
253 c = -remainarea*signorm/(_Yarray[ivaly+1]-_Yarray[ivaly]);
255 a = (sigofy[ivaly+1]-sigofy[ivaly])/2.;
261 remainy = (-b + sqrt(b*b-4.*a*c))/(2.*a);
263 //calculate the y value
264 y = _Yarray[ivaly]+(_Yarray[ivaly+1]-_Yarray[ivaly])*remainy;
270 //______________________________________________________________________________
271 void Gammagammaleptonpair::pairMomentum(double w,double y,double &E,double &px,double &py,double &pz)
273 //this function calculates px,py,pz,and E given w and y
275 double anglepp1=0.,anglepp2=0.,pp1=0.,pp2=0.,E1=0.,E2=0.,signpx=0.,pt=0.;
277 //E1 and E2 are for the 2 photons in the CM frame
281 //calculate px and py
282 //to get x and y components-- phi is random between 0 and 2*pi
283 anglepp1 = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
284 anglepp2 = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
288 px = pp1*cos(2.*starlightConstants::pi*anglepp1)+pp2*cos(2.*starlightConstants::pi*anglepp2);
289 py = pp1*sin(2.*starlightConstants::pi*anglepp1)+pp2*sin(2.*starlightConstants::pi*anglepp2);
291 //Compute vector sum Pt=Pt1+Pt2 to find pt for the produced particle
292 pt = sqrt(px*px+py*py);
294 //W is the mass of the produced particle (not necessarily on-mass-shell).Now compute its energy and pz
295 E = sqrt(w*w+pt*pt)*cosh(y);
296 pz= sqrt(w*w+pt*pt)*sinh(y);
297 signpx = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
299 //pick the z direction
300 //Don't do this anymore since y goes from -ymax to +ymax (JN 15-02-2013)
301 //if(signpx > 0.5) pz = -pz;
305 //______________________________________________________________________________
306 double Gammagammaleptonpair::pp_1(double E)
308 // This is for beam 1
309 // returns on random draw from pp(E) distribution
310 double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
311 double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
314 ereds = (E/_GGlepInputGamma_em)*(E/_GGlepInputGamma_em);
315 //sqrt(3)*E/gamma_em is p_t where the distribution is a maximum
316 Cm = sqrt(3.)*E/_GGlepInputGamma_em;
317 //the amplitude of the p_t spectrum at the maximum
318 singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
319 Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
321 //pick a test value pp, and find the amplitude there
322 x = randyInstance.Rndom();
323 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
324 singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
325 test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
328 u = randyInstance.Rndom();
334 x =randyInstance.Rndom();
335 pp = 5*starlightConstants::hbarc/_bbs.beam1().nuclearRadius()*x;
336 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
337 test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
344 //______________________________________________________________________________
345 double Gammagammaleptonpair::pp_2(double E)
348 // This is for beam 2
349 //returns on random draw from pp(E) distribution
350 double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
351 double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
354 ereds = (E/_GGlepInputGamma_em)*(E/_GGlepInputGamma_em);
355 //sqrt(3)*E/gamma_em is p_t where the distribution is a maximum
356 Cm = sqrt(3.)*E/_GGlepInputGamma_em;
357 //the amplitude of the p_t spectrum at the maximum
358 singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
359 Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
361 //pick a test value pp, and find the amplitude there
362 x = randyInstance.Rndom();
363 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius(); //Will use nucleus #1
364 singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
365 test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
368 u = randyInstance.Rndom();
374 x =randyInstance.Rndom();
375 pp = 5*starlightConstants::hbarc/_bbs.beam2().nuclearRadius()*x;
376 singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds);
377 test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
386 //______________________________________________________________________________
387 void Gammagammaleptonpair::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
388 double , // E (unused)
390 double px0, double py0, double pz0,
391 double& px1, double& py1, double& pz1, double& E1,
392 double& px2, double& py2, double& pz2, double& E2,
396 // This routine decays a particle into two particles of mass mdec,
397 // taking spin into account
399 double pmag, anglelep[20001];
400 // double ytest=0.,dndtheta;
401 double phi,theta,xtest,Ecm;
402 double betax,betay,betaz;
403 double hirestheta,hirestest,hiresw; //added from JN->needed precision
408 cout<<" ERROR: W="<<W<<endl;
412 pmag = sqrt(W*W/4. - mdec*mdec);
414 // pick an orientation, based on the spin
415 // phi has a flat distribution in 2*pi
416 phi = randyInstance.Rndom()*2.*starlightConstants::pi; //(random()/(RAND_MAX+1.0))* 2.*starlightConstants::pi;
418 // find theta, the angle between one of the outgoing particles and
419 // the beamline, in the frame of the two photons
422 // calculate a table of integrated angle values for leptons
423 // JN05: Go from 1000->20000bins, use double precision for anglelep and thetalep. needed when W >6Gev.
428 for(int i =1;i<=20000;i++)
430 hirestheta = starlightConstants::pi * double(i) /20000.;
432 // Added sin(theta) phase space factor (not in thetalep) and changed E to W in thetalep call
434 // Note that thetalep is form invariant, so it could be called for E, theta_lab just
435 // as well as W,theta_cm. Since there is a boost from cm to lab below, the former is fine.
437 anglelep[i] = anglelep[i-1] + thetalep(hiresw,hirestheta)*sin(hirestheta);
441 xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
443 for(int i =1;i<=20000;i++)
445 if(xtest > (anglelep[i]/anglelep[20000]))
446 hirestheta = starlightConstants::pi * double(i) / 20000.;
453 cout<<" This model cannot yet handle this spin value for lepton pairs: "<<getSpin()<<endl;
456 // compute unboosted momenta
457 px1 = sin(theta)*cos(phi)*pmag;
458 py1 = sin(theta)*sin(phi)*pmag;
459 pz1 = cos(theta)*pmag;
465 //Changed mass to W 11/9/2000 SRK
466 Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0);
468 E1 = sqrt(mdec*mdec+px1*px1+py1*py1+pz1*pz1);
469 E2 = sqrt(mdec*mdec+px2*px2+py2*py2+pz2*pz2);
470 // decay tau to electrons
471 // note that after this routine px1, etc., refer to the electrons
472 if(_GGlepInputpidtest == starlightConstants::TAUONDECAY)
473 tauDecay(px1,py1,pz1,E1,px2,py2,pz2,E2);
475 // Lorentz transform into the lab frame
476 // betax,betay,betaz are the boost of the complete system
481 transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent);
482 transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent);
488 // change particle id from that of parent to that of daughters
489 // change taoun id into electron id, already switched particles in tau decay
490 if(_GGlepInputpidtest == starlightConstants::TAUONDECAY)
491 ipid = starlightConstants::ELECTRON;
492 // electrons remain electrons; muons remain muons
493 if ( (_GGlepInputpidtest == starlightConstants::ELECTRON) || (_GGlepInputpidtest == starlightConstants::MUON) ||
494 (_GGlepInputpidtest == starlightConstants::TAUON) )
495 ipid = _GGlepInputpidtest;
499 //______________________________________________________________________________
500 double Gammagammaleptonpair::thetalep(double W,double theta)
502 // This function calculates the cross-section as a function of
503 // angle for a given W and Y, for the production of two muons.
505 // expression is taken from Brodsky et al. PRD 1971, 1532
507 // factor that are not dependant on theta are scrapped, so the
508 // the absolute crosssections given by this function are inaccurate
509 // here we are working in the CM frame of the photons and the last
512 // real function thetalep (W,theta)
514 double moverw=0., W1sq=0.;
515 double thetalep_r=0.,denominator=0.;
517 W1sq = (W / 2.)*(W/2.);
518 moverw = getMass()*getMass() / W1sq;
519 denominator = (1. - (1. - moverw)*(cos(theta)*cos(theta)));
521 thetalep_r = 2. + 4.*(1.-moverw)*((1.-moverw)*sin(theta)*sin(theta)*cos(theta)*cos(theta) + moverw) / (denominator*denominator);
528 //______________________________________________________________________________
529 starlightConstants::event Gammagammaleptonpair::produceEvent(int &ievent)
530 {//returns the vector with the decay particles inside.
531 starlightConstants::event leptonpair; //This object will store all the tracks for a single event
532 double comenergy = 0.;
533 double rapidity = 0.;
535 double pairmomx=0.,pairmomy=0.,pairmomz=0.;
537 starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
539 double px2=0.,px1=0.,py2=0.,pt2=0.,py1=0.,pz2=0.,pz1=0.,pt1=0.;
540 double E1=0.,E2=0.,mdec=0.;
541 //this function decays particles and writes events to a file
542 //zero out the event structure
543 leptonpair._numberOfTracks=0;
549 leptonpair._fsParticle[i]=starlightConstants::UNKNOWN;
550 leptonpair._charge[i]=0;
557 pairMomentum(comenergy,rapidity,pairE,pairmomx,pairmomy,pairmomz);
558 twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,E1,px1,px2,py2,E2,mdec,iFbadevent);
560 const bool xtest(randyInstance.Rndom() < 0.5);
561 const int charges[2] = {
565 leptonpair._numberOfTracks=2;//leptonpairs are two tracks...
566 leptonpair.px[0]=px1;
567 leptonpair.py[0]=py1;
568 leptonpair.pz[0]=pz1;
569 leptonpair._fsParticle[0]=ipid;
570 leptonpair._charge[0]=charges[0];
572 leptonpair.px[1]=px2;
573 leptonpair.py[1]=py2;
574 leptonpair.pz[1]=pz2;
575 leptonpair._fsParticle[1]=ipid;
576 leptonpair._charge[1]=charges[1];
585 //______________________________________________________________________________
586 upcEvent Gammagammaleptonpair::produceEvent()
588 //returns the vector with the decay particles inside.
592 double comenergy = 0.;
593 double rapidity = 0.;
595 double pairmomx=0.,pairmomy=0.,pairmomz=0.;
597 starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
599 double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.,E2=0.,E1=0.,mdec=0.;
600 bool accepted = false;
602 //this function decays particles and writes events to a file
603 //zero out the event structure
608 pairMomentum(comenergy,rapidity,pairE,pairmomx,pairmomy,pairmomz);
612 twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,E1,px2,py2,pz2,E2,mdec,iFbadevent);
613 double pt1chk = sqrt(px1*px1+py1*py1);
614 double pt2chk = sqrt(px2*px2+py2*py2);
616 double eta1 = pseudoRapidity(px1, py1, pz1);
617 double eta2 = pseudoRapidity(px2, py2, pz2);
619 if(_ptCutEnabled && !_etaCutEnabled){
620 if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){
625 else if(!_ptCutEnabled && _etaCutEnabled){
626 if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
631 else if(_ptCutEnabled && _etaCutEnabled){
632 if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){
633 if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
639 else if(!_ptCutEnabled && !_etaCutEnabled)
642 }while((_ptCutEnabled || _etaCutEnabled) && !accepted);
643 //twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
646 const bool xtest(randyInstance.Rndom() < 0.5);
647 const int charges[2] = {
653 starlightParticle particle1(px1, py1, pz1, E1, mdec, -charges[0]*ipid, charges[0]);
654 event.addParticle(particle1);
656 starlightParticle particle2(px2, py2, pz2, E2, mdec, -charges[1]*ipid, charges[1]);
657 event.addParticle(particle2);
663 //______________________________________________________________________________
664 void Gammagammaleptonpair::calculateTable()
666 // this subroutine calculates the tables that are used
667 // elsewhere in the montecarlo
668 // the tauon decay is taken from V-A theory, 1 - 1/3 cos(theta)
669 // the energy of each of the two leptons in tau decay
670 // is calculated using formula 10.35 given
671 // in introduction to elementary particles by D. griffiths
672 // which assmes that the mass of the electron is 0.
673 // the highest energy attainable by the electron in such a system is
674 // .5 * mass of the tau
676 // subroutine calculateTable
681 _tautolangle[0] = 0.;
685 for(int i =1;i<=100;i++)
687 // calculate energy of tau decay
688 E = double(i)/100. * .5 * starlightConstants::tauMass;
689 _dgammade[i] = _dgammade[i-1] + E*E * (1. - 4.*E/(3.*starlightConstants::tauMass));
691 // calculate angles for tau
692 theta = starlightConstants::pi * double(i) / 100.;
693 _tautolangle[i] = _tautolangle[i-1] + (1 + 0.3333333 * cos(theta));
700 //______________________________________________________________________________
701 void Gammagammaleptonpair::tauDecay(double &px1,double &py1,double &pz1,double &E1,double &px2,double &py2,double &pz2,double &E2)
703 // this routine assumes that the tauons decay to electrons and
704 // calculates the directons of the decays
706 double Ee1,Ee2,theta1,theta2,phi1,phi2, ran1, ran2 ;
707 double pmag1,pex1,pey1,pez1,pex2,pey2,pez2,pmag2;
708 double betax,betay,betaz,dir;
710 int Tmp_Par=0; // temp variable for the transform function .. kind of starnge - being called with 7 parameter instead of 8
712 // the highest energy attainable by the electron in this system is
713 // .5 * mass of the tau
715 // get two random numbers to compare with
718 ran1 = randyInstance.Rndom()*_dgammade[100];//(random()/(RAND_MAX+1.0)) * _dgammade[100];
719 ran2 = randyInstance.Rndom()*_dgammade[100];//(random()/(RAND_MAX+1.0)) * _dgammade[100];
721 // compute the energies that correspond to those numbers
725 for( int i =1;i<=100;i++)
727 if (ran1 > _dgammade[i])
728 Ee1 = double(i) /100. * .5 * getMass();
729 if (ran2 > _dgammade[i])
730 Ee2 = double(i) /100. * .5 * getMass();
733 // to find the direction of emmission, first
734 // we determine if the tauons have spin of +1 or -1 along the
735 // direction of the beam line
737 if ( randyInstance.Rndom() < 0.5 )//(random()/(RAND_MAX+1.0)) < 0.5)
740 // get two random numbers to compare with
741 ran1 = randyInstance.Rndom()*_tautolangle[99];//(random()/(RAND_MAX+1.0)) * _tautolangle[100];
742 ran2 = randyInstance.Rndom()*_tautolangle[99];//(random()/(RAND_MAX+1.0)) * _tautolangle[100];
744 // find the angles corrsponding to those numbers
747 for( int i =1;i<=100;i++)
749 if (ran1 > _tautolangle[i]) theta1 = starlightConstants::pi * double(i) /100.;
750 if (ran2 > _tautolangle[i]) theta2 = starlightConstants::pi * double(i) /100.;
753 // grab another two random numbers to determine phi's
754 phi1 = randyInstance.Rndom()*2.*starlightConstants::pi;// (random()/(RAND_MAX+1.0))* 2. * starlightConstants::pi;
755 phi2 = randyInstance.Rndom()*2.*starlightConstants::pi;// (random()/(RAND_MAX+1.0))* 2. * starlightConstants::pi;
756 // figure out the momenta of the electron in the frames of the
757 // tauons from which they decayed, that is electron1 is in the
758 // rest frame of tauon1 and e2 is in the rest fram of tau2
759 // again the electrons are assumed to be massless
761 pex1 = cos(phi1)*sin(theta1)*pmag1;
762 pey1 = sin(phi1)*sin(theta1)*pmag1;
763 pez1 = cos(theta1)*pmag1*dir;
765 pex2 = cos(phi2)*sin(theta2)*pmag2;
766 pey2 = sin(phi2)*sin(theta2)*pmag2;
767 pez2 = cos(theta2)*pmag2*(-1.*dir);
768 // now Lorentz transform into the frame of each of the particles
769 // do particle one first
773 //cout<<"2decay betax,pex1= "<<betax<<" "<<pex1<<endl;
774 transform (betax,betay,betaz,Ee1,pex1,pey1,pez1,Tmp_Par);
775 // then do particle two
780 transform (betax,betay,betaz,Ee2,pex2,pey2,pez2, Tmp_Par);
781 // finally dump the electron values into the approriate
794 //______________________________________________________________________________
795 double Gammagammaleptonpair::getMass()
797 double leptonmass=0.;
798 switch(_GGlepInputpidtest){
799 case starlightConstants::ELECTRON:
800 leptonmass=starlightConstants::mel;
802 case starlightConstants::MUON:
803 leptonmass=starlightConstants::muonMass;
805 case starlightConstants::TAUON:
806 leptonmass=starlightConstants::tauMass;
808 case starlightConstants::TAUONDECAY:
809 leptonmass=starlightConstants::tauMass;
812 cout<<"Not a recognized lepton, Gammagammaleptonpair::getmass(), mass = 0."<<endl;
819 //______________________________________________________________________________
820 double Gammagammaleptonpair::getWidth()
822 double leptonwidth=0.;
828 //______________________________________________________________________________
829 double Gammagammaleptonpair::getSpin()
831 double leptonspin=0.5;