STARLIGHT/starlight/src/gammagammaleptonpair.cpp

   1 ///////////////////////////////////////////////////////////////////////////
   2 //
   3 //    Copyright 2010
   4 //
   5 //    This file is part of starlight.
   6 //
   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.
  11 //
  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.
  16 //
  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/>.
  19 //
  20 ///////////////////////////////////////////////////////////////////////////
  21 //
  22 // File and Version Information:
  23 // $Rev:: 176                         $: revision of last commit
  24 // $Author:: jseger                   $: author of last commit
  25 // $Date:: 2014-06-20 22:15:20 +0200 #$: date of last commit
  26 //
  27 // Description:
  28 //    Nystrand 220710
  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
  33 //    with old starligt.
  34 //
  35 //
  36 ///////////////////////////////////////////////////////////////////////////
  37
  38
  39 #include <iostream>
  40 #include <fstream>
  41 #include <cmath>
  42 #include <vector>
  43
  44 #include "starlightconstants.h"
  45 #include "gammagammaleptonpair.h"
  46
  47
  48 using namespace std;
  49
  50
  51 //_____________________________________________________________________________
  52 Gammagammaleptonpair::Gammagammaleptonpair(beamBeamSystem& bbsystem)
  53 : eventChannel(bbsystem)
  54 {
  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
  63     read();
  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();
  68 }
  69
  70
  71 //______________________________________________________________________________
  72 Gammagammaleptonpair::~Gammagammaleptonpair()
  73 { }
  74
  75
  76 //______________________________________________________________________________
  77 void Gammagammaleptonpair::twoLeptonCrossSection()
  78 {
  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
  84
  85     for(int i=0;i<_GGlepInputnumw;i++)
  86     {
  87         for(int j=0;j<_GGlepInputnumy;j++)
  88         {
  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]);
  91         }
  92     }
  93     //calculate the total two-lepton cross section
  94     double sigmasum =0.;
  95     for(int i=0;i<_GGlepInputnumw-1;i++)
  96     {
  97         for(int j=0;j<_GGlepInputnumy-1;j++)
  98         {
  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]);
 102         }
 103     }
 104     cout << "The total "<<_GGlepInputpidtest<<" cross-section is: "<<sigmasum<<" barns."<<endl;
 105
 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
 108
 109     double sgf=0.;
 110
 111     for(int i=0;i<_ReadInputnumw;i++)
 112     {
 113             _sigofw[i]=0.;
 114             for(int j=0;j<_ReadInputnumy-1;j++)
 115             {
 116                 _sigofw[i] = _sigofw[i]+(_Yarray[j+1]-_Yarray[j])*(_sigmax[i][j+1]+_sigmax[i][j])/2.;
 117             }
 118     }
 119
 120     //calculate the unnormalized sgfint array
 121     _sigfint[0]=0.;
 122     for(int i=0;i<_ReadInputnumw-1;i++)
 123     {
 124         sgf=(_sigofw[i+1]+_sigofw[i])*(_Warray[i+1]-_Warray[i])/2.;
 125         _sigfint[i+1]=_sigfint[i]+sgf;
 126      }
 127
 128      //normalize sgfint array
 129      _signormw=_sigfint[_ReadInputnumw-1];
 130      for(int i=0;i<_ReadInputnumw;i++)
 131      {
 132           _sigfint[i]=_sigfint[i]/_signormw;
 133      }
 134     return;
 135 }
 136
 137
 138 //______________________________________________________________________________
 139 double Gammagammaleptonpair::twoMuonCrossSection(double w)
 140 {
 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.;
 144     s = w*w;
 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);
 151     if(Etest > 1.)
 152         sigmuix = 0.;
 153     return sigmuix;
 154 }
 155
 156
 157 //______________________________________________________________________________
 158 void Gammagammaleptonpair::pickw(double &w)
 159 {
 160 //  This function picks a w for the 2-photon calculation.
 161
 162     double x=0.,remainarea=0.,remainw=0.,a=0.,b=0.,c=0.;
 163     int ivalw=0;
 164
 165     if(_wdelta != 0)
 166     {
 167         w=_wdelta;
 168         ivalw=_ivalwd;
 169         remainw=_remainwd;
 170     }
 171     else{
 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
 175
 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++)
 180         {
 181             if(x > _sigfint[i]) ivalw=i;
 182         }
 183         //remainder above ivalw
 184         remainarea = x - _sigfint[ivalw];
 185
 186         //figure out what point corresponds to the excess area in remainarea
 187         c = -remainarea*_signormw/(_Warray[ivalw+1]-_Warray[ivalw]);
 188         b = _sigofw[ivalw];
 189         a = (_sigofw[ivalw+1]-_sigofw[ivalw])/2.;
 190         if(a==0.)
 191         {
 192             remainw = -c/b;
 193         }
 194         else{
 195             remainw = (-b+sqrt(b*b-4.*a*c))/(2.*a);
 196         }
 197         _ivalwd = ivalw;
 198         _remainwd = remainw;
 199         //calculate the w value
 200         w = _Warray[ivalw]+(_Warray[ivalw+1]-_Warray[ivalw])*remainw;
 201     }
 202 }
 203
 204
 205 //______________________________________________________________________________
 206 void Gammagammaleptonpair::picky(double &y)
 207 {
 208     // This function picks a y given a W
 209
 210     double * sigofy;
 211     double * sgfint;
 212     sigofy = new double[starlightLimits::MAXYBINS];
 213     sgfint = new double[starlightLimits::MAXWBINS];
 214
 215     double remainw =0.,remainarea=0.,remainy=0.,a=0.,b=0.,c=0.,sgf=0.,signorm=0.,x=0.;
 216     int ivalw=0,ivaly=0;
 217
 218     ivalw=_ivalwd;
 219     remainw=_remainwd;
 220     //average over two colums to get y array
 221     for(int j=0;j<_GGlepInputnumy;j++)
 222     {
 223         sigofy[j]=_sigmax[ivalw][j]+(_sigmax[ivalw+1][j]-_sigmax[ivalw][j])*remainw;
 224     }
 225
 226     //calculate the unnormalized sgfint
 227     sgfint[0]=0.;
 228     for(int j=0;j<_GGlepInputnumy-1;j++)
 229     {
 230         sgf = (sigofy[j+1]+sigofy[j])/2.;
 231         sgfint[j+1]=sgfint[j]+sgf*(_Yarray[j+1]-_Yarray[j]);
 232     }
 233
 234     //normalize the sgfint array
 235     signorm = sgfint[_GGlepInputnumy-1];
 236
 237     for(int j=0;j<_GGlepInputnumy;j++)
 238     {
 239         sgfint[j]=sgfint[j]/signorm;
 240     }
 241
 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++)
 246     {
 247         if(x > sgfint[i]) ivaly = i;
 248     }
 249     //remainder above ivaly
 250     remainarea = x - sgfint[ivaly];
 251
 252     //figure what point corresponds to the leftover area in remainarea
 253     c = -remainarea*signorm/(_Yarray[ivaly+1]-_Yarray[ivaly]);
 254     b = sigofy[ivaly];
 255     a = (sigofy[ivaly+1]-sigofy[ivaly])/2.;
 256     if(a==0.)
 257     {
 258         remainy = -c/b;
 259     }
 260     else{
 261         remainy = (-b + sqrt(b*b-4.*a*c))/(2.*a);
 262     }
 263     //calculate the y value
 264     y = _Yarray[ivaly]+(_Yarray[ivaly+1]-_Yarray[ivaly])*remainy;
 265     delete[] sigofy;
 266     delete[] sgfint;
 267 }
 268
 269
 270 //______________________________________________________________________________
 271 void Gammagammaleptonpair::pairMomentum(double w,double y,double &E,double &px,double &py,double &pz)
 272 {
 273     //this function calculates px,py,pz,and E given w and y
 274
 275     double anglepp1=0.,anglepp2=0.,pp1=0.,pp2=0.,E1=0.,E2=0.,signpx=0.,pt=0.;
 276
 277     //E1 and E2 are for the 2 photons in the CM frame
 278     E1 = w*exp(y)/2.;
 279     E2 = w*exp(-y)/2.;
 280
 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);
 285
 286     pp1 = pp_1(E1);
 287     pp2 = pp_2(E2);
 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);
 290
 291     //Compute vector sum Pt=Pt1+Pt2 to find pt for the produced particle
 292     pt = sqrt(px*px+py*py);
 293
 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);
 298
 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;
 302 }
 303
 304
 305 //______________________________________________________________________________
 306 double Gammagammaleptonpair::pp_1(double E)
 307 {
 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.;
 312     int satisfy =0;
 313
 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)));
 320
 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)));
 326
 327     while(satisfy==0){
 328         u = randyInstance.Rndom();
 329         if(u*Coef <= test)
 330         {
 331             satisfy =1;
 332         }
 333         else{
 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));
 338         }
 339     }
 340
 341     return pp;
 342 }
 343
 344 //______________________________________________________________________________
 345 double Gammagammaleptonpair::pp_2(double E)
 346 {
 347
 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.;
 352     int satisfy =0;
 353
 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)));
 360
 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)));
 366
 367     while(satisfy==0){
 368         u = randyInstance.Rndom();
 369         if(u*Coef <= test)
 370         {
 371             satisfy =1;
 372         }
 373         else{
 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));
 378         }
 379     }
 380
 381     return pp;
 382 }
 383
 384
 385
 386 //______________________________________________________________________________
 387 void Gammagammaleptonpair::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
 388                                     double  ,  // E (unused)
 389                                     double  W,
 390                                     double  px0, double  py0, double  pz0,
 391                                     double& px1, double& py1, double& pz1,
 392                                     double& px2, double& py2, double& pz2,
 393                                     int&    iFbadevent)
 394 {
 395     //     This routine decays a particle into two particles of mass mdec,
 396     //     taking spin into account
 397
 398     double mdec=0.,E1=0.,E2=0.;
 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
 404
 405     mdec = getMass();
 406     if(W < 2*mdec)
 407     {
 408         cout<<" ERROR: W="<<W<<endl;
 409         iFbadevent = 1;
 410         return;
 411     }
 412     pmag = sqrt(W*W/4. - mdec*mdec);
 413
 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;
 417
 418     //     find theta, the angle between one of the outgoing particles and
 419     //    the beamline, in the frame of the two photons
 420
 421     if(getSpin()==0.5){
 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.
 424         hiresw = W;
 425
 426         anglelep[0] = 0.;
 427
 428         for(int i =1;i<=20000;i++)
 429         {
 430             hirestheta = starlightConstants::pi * double(i) /20000.;
 431
 432             //  Added sin(theta) phase space factor (not in thetalep) and changed E to W in thetalep call
 433             //  11/9/2000 SRK
 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.
 436
 437             anglelep[i] = anglelep[i-1] + thetalep(hiresw,hirestheta)*sin(hirestheta);
 438         }
 439
 440         hirestheta = 0.;
 441         xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 442         hirestest = xtest;
 443         for(int i =1;i<=20000;i++)
 444         {
 445             if(xtest > (anglelep[i]/anglelep[20000]))
 446                 hirestheta = starlightConstants::pi * double(i) / 20000.;
 447         }
 448         theta=hirestheta;
 449
 450     }
 451
 452     if(getSpin() != 0.5)
 453         cout<<" This model cannot yet handle this spin value for lepton pairs: "<<getSpin()<<endl;
 454
 455
 456     //     compute unboosted momenta
 457     px1 = sin(theta)*cos(phi)*pmag;
 458     py1 = sin(theta)*sin(phi)*pmag;
 459     pz1 = cos(theta)*pmag;
 460     px2 = -px1;
 461     py2 = -py1;
 462     pz2 = -pz1;
 463
 464     //        compute energies
 465     //Changed mass to W 11/9/2000 SRK
 466     Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0);
 467
 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);
 474
 475     //     Lorentz transform into the lab frame
 476     // betax,betay,betaz are the boost of the complete system
 477     betax = -(px0/Ecm);
 478     betay = -(py0/Ecm);
 479     betaz = -(pz0/Ecm);
 480
 481     transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent);
 482     transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent);
 483
 484
 485     if(iFbadevent == 1)
 486         return;
 487
 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;
 496 }
 497
 498
 499 //______________________________________________________________________________
 500 double Gammagammaleptonpair::thetalep(double W,double theta)
 501 {
 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.
 504     //     (or tauons)
 505     //    expression is taken from Brodsky et al. PRD 1971, 1532
 506     //     equation 5.7
 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
 510     //     term is 0
 511
 512     //    real function thetalep (W,theta)
 513
 514     double moverw=0., W1sq=0.;
 515     double thetalep_r=0.,denominator=0.;
 516
 517     W1sq = (W / 2.)*(W/2.);
 518     moverw = getMass()*getMass() / W1sq;
 519     denominator = (1. - (1. - moverw)*(cos(theta)*cos(theta)));
 520
 521     thetalep_r = 2. + 4.*(1.-moverw)*((1.-moverw)*sin(theta)*sin(theta)*cos(theta)*cos(theta) + moverw) / (denominator*denominator);
 522
 523     return thetalep_r;
 524
 525 }
 526
 527
 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.;
 534     double pairE = 0.;
 535     double pairmomx=0.,pairmomy=0.,pairmomz=0.;
 536     int iFbadevent=0;
 537     starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
 538
 539     double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
 540 //this function decays particles and writes events to a file
 541     //zero out the event structure
 542     leptonpair._numberOfTracks=0;
 543     for(int i=0;i<4;i++)
 544     {
 545         leptonpair.px[i]=0.;
 546         leptonpair.py[i]=0.;
 547         leptonpair.pz[i]=0.;
 548         leptonpair._fsParticle[i]=starlightConstants::UNKNOWN;
 549         leptonpair._charge[i]=0;
 550     }
 551
 552     pickw(comenergy);
 553
 554     picky(rapidity);
 555
 556     pairMomentum(comenergy,rapidity,pairE,pairmomx,pairmomy,pairmomz);
 557     twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
 558     if (iFbadevent==0){
 559         int q1=0,q2=0;
 560
 561         double xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 562         if (xtest<0.5)
 563         {
 564             q1=1;
 565             q2=-1;
 566         }
 567         else{
 568             q1=-1;
 569             q2=1;
 570         }
 571         leptonpair._numberOfTracks=2;//leptonpairs are two tracks...
 572         leptonpair.px[0]=px1;
 573         leptonpair.py[0]=py1;
 574         leptonpair.pz[0]=pz1;
 575         leptonpair._fsParticle[0]=ipid;
 576         leptonpair._charge[0]=q1;
 577
 578         leptonpair.px[1]=px2;
 579         leptonpair.py[1]=py2;
 580         leptonpair.pz[1]=pz2;
 581         leptonpair._fsParticle[1]=ipid;
 582         leptonpair._charge[1]=q2;
 583
 584         ievent=ievent+1;
 585     }
 586
 587     return leptonpair;
 588 }
 589
 590
 591 //______________________________________________________________________________
 592 upcEvent Gammagammaleptonpair::produceEvent()
 593 {
 594 //returns the vector with the decay particles inside.
 595
 596    upcEvent event;
 597
 598    double comenergy = 0.;
 599    double rapidity = 0.;
 600    double pairE = 0.;
 601    double pairmomx=0.,pairmomy=0.,pairmomz=0.;
 602    int iFbadevent=0;
 603    starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
 604
 605    double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.,E2=0.,E1=0.;
 606    bool accepted = false;
 607    do{
 608      //this function decays particles and writes events to a file
 609      //zero out the event structure
 610      pickw(comenergy);
 611
 612      picky(rapidity);
 613
 614      pairMomentum(comenergy,rapidity,pairE,pairmomx,pairmomy,pairmomz);
 615
 616
 617      _nmbAttempts++;
 618      twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
 619      double pt1chk = sqrt(px1*px1+py1*py1);
 620      double pt2chk = sqrt(px2*px2+py2*py2);
 621
 622      double eta1 = pseudoRapidity(px1, py1, pz1);
 623      double eta2 = pseudoRapidity(px2, py2, pz2);
 624
 625      if(_ptCutEnabled && !_etaCutEnabled){
 626        if(pt1chk > _ptCutMin && pt1chk < _ptCutMax &&  pt2chk > _ptCutMin && pt2chk < _ptCutMax){
 627          accepted = true;
 628          _nmbAccepted++;
 629        }
 630      }
 631      else if(!_ptCutEnabled && _etaCutEnabled){
 632        if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
 633          accepted = true;
 634          _nmbAccepted++;
 635        }
 636      }
 637      else if(_ptCutEnabled && _etaCutEnabled){
 638        if(pt1chk > _ptCutMin && pt1chk < _ptCutMax &&  pt2chk > _ptCutMin && pt2chk < _ptCutMax){
 639          if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
 640            accepted = true;
 641             _nmbAccepted++;
 642          }
 643        }
 644      }
 645      else if(!_ptCutEnabled && !_etaCutEnabled)
 646         _nmbAccepted++;
 647
 648    }while((_ptCutEnabled || _etaCutEnabled) && !accepted);
 649    //twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
 650
 651    if (iFbadevent==0){
 652      int q1=0,q2=0;
 653
 654      double xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 655      if (xtest<0.5)
 656        {
 657          q1=1;
 658          q2=-1;
 659        }
 660      else{
 661        q1=-1;
 662        q2=1;
 663      }
 664
 665      // The new stuff
 666      double mlepton = getMass();
 667      E1 = sqrt( mlepton*mlepton + px1*px1 + py1*py1 + pz1*pz1 );
 668      E2 = sqrt( mlepton*mlepton + px2*px2 + py2*py2 + pz2*pz2 );
 669
 670      starlightParticle particle1(px1, py1, pz1, E1, starlightConstants::UNKNOWN, -q1*ipid, q1);
 671      event.addParticle(particle1);
 672
 673      starlightParticle particle2(px2, py2, pz2, E2, starlightConstants::UNKNOWN, -q2*ipid, q2);
 674      event.addParticle(particle2);
 675
 676     }
 677    return event;
 678 }
 679
 680
 681 //______________________________________________________________________________
 682 void Gammagammaleptonpair::calculateTable()
 683 {
 684     //     this subroutine calculates the tables that are used
 685     //     elsewhere in the montecarlo
 686     //     the tauon decay is taken from V-A theory, 1 - 1/3 cos(theta)
 687     //     the energy of each of the two leptons in tau decay
 688     //     is calculated using formula 10.35 given
 689     //     in introduction to elementary particles by D. griffiths
 690     //     which assmes that the mass of the electron is 0.
 691     //     the highest energy attainable by the electron in such a system is
 692     //     .5 * mass of the tau
 693
 694     //    subroutine calculateTable
 695
 696
 697     double E,theta;
 698
 699     _tautolangle[0] = 0.;
 700     _dgammade[0] = 0.;
 701
 702
 703     for(int i =1;i<=100;i++)
 704     {
 705         //     calculate energy of tau decay
 706         E = double(i)/100. * .5 * starlightConstants::tauMass;
 707         _dgammade[i] = _dgammade[i-1] + E*E * (1. - 4.*E/(3.*starlightConstants::tauMass));
 708
 709         //     calculate angles for tau
 710         theta = starlightConstants::pi * double(i) / 100.;
 711         _tautolangle[i] = _tautolangle[i-1] + (1 + 0.3333333 * cos(theta));
 712     }
 713
 714
 715 }
 716
 717
 718 //______________________________________________________________________________
 719 void Gammagammaleptonpair::tauDecay(double &px1,double &py1,double &pz1,double &E1,double &px2,double &py2,double &pz2,double &E2)
 720 {
 721     //     this routine assumes that the tauons decay to electrons and
 722     //     calculates the directons of the decays
 723
 724     double Ee1,Ee2,theta1,theta2,phi1,phi2, ran1, ran2 ;
 725     double pmag1,pex1,pey1,pez1,pex2,pey2,pez2,pmag2;
 726     double betax,betay,betaz,dir;
 727
 728     int Tmp_Par=0; // temp variable for the transform function .. kind of starnge - being called with 7 parameter instead of 8
 729
 730     //     the highest energy attainable by the electron in this system is
 731     //     .5 * mass of the tau
 732
 733     //     get two random numbers to compare with
 734
 735
 736     ran1 = randyInstance.Rndom()*_dgammade[100];//(random()/(RAND_MAX+1.0)) * _dgammade[100];
 737     ran2 = randyInstance.Rndom()*_dgammade[100];//(random()/(RAND_MAX+1.0)) * _dgammade[100];
 738
 739     //     compute the energies that correspond to those numbers
 740     Ee1 = 0.;
 741     Ee2 = 0.;
 742
 743     for( int i =1;i<=100;i++)
 744     {
 745         if (ran1 > _dgammade[i])
 746             Ee1 = double(i) /100. * .5 * getMass();
 747         if (ran2 > _dgammade[i])
 748             Ee2 = double(i) /100. * .5 * getMass();
 749     }
 750
 751     //     to find the direction of emmission, first
 752     //     we determine if the tauons have spin of +1 or -1 along the
 753     //     direction of the beam line
 754     dir = 1.;
 755     if ( randyInstance.Rndom() < 0.5 )//(random()/(RAND_MAX+1.0)) < 0.5)
 756         dir = -1.;
 757
 758     //     get two random numbers to compare with
 759     ran1 = randyInstance.Rndom()*_tautolangle[99];//(random()/(RAND_MAX+1.0))  * _tautolangle[100];
 760     ran2 = randyInstance.Rndom()*_tautolangle[99];//(random()/(RAND_MAX+1.0))  * _tautolangle[100];
 761
 762     //     find the angles corrsponding to those numbers
 763     theta1 = 0.;
 764     theta2 = 0.;
 765     for( int i =1;i<=100;i++)
 766     {
 767         if (ran1 > _tautolangle[i]) theta1 = starlightConstants::pi * double(i) /100.;
 768         if (ran2 > _tautolangle[i]) theta2 = starlightConstants::pi * double(i) /100.;
 769     }
 770
 771     //     grab another two random numbers to determine phi's
 772     phi1 = randyInstance.Rndom()*2.*starlightConstants::pi;// (random()/(RAND_MAX+1.0))* 2. * starlightConstants::pi;
 773     phi2 = randyInstance.Rndom()*2.*starlightConstants::pi;// (random()/(RAND_MAX+1.0))* 2. * starlightConstants::pi;
 774     //     figure out the momenta of the electron in the frames of the
 775     //     tauons from which they decayed, that is electron1 is in the
 776     //     rest frame of tauon1 and e2 is in the rest fram of tau2
 777     //     again the electrons are assumed to be massless
 778     pmag1 = Ee1;
 779     pex1 = cos(phi1)*sin(theta1)*pmag1;
 780     pey1 = sin(phi1)*sin(theta1)*pmag1;
 781     pez1 = cos(theta1)*pmag1*dir;
 782     pmag2 = Ee2;
 783     pex2 = cos(phi2)*sin(theta2)*pmag2;
 784     pey2 = sin(phi2)*sin(theta2)*pmag2;
 785     pez2 = cos(theta2)*pmag2*(-1.*dir);
 786     //     now Lorentz transform into the frame of each of the particles
 787     //     do particle one first
 788     betax = -(px1/E1);
 789     betay = -(py1/E1);
 790     betaz = -(pz1/E1);
 791 //cout<<"2decay betax,pex1= "<<betax<<" "<<pex1<<endl;
 792     transform (betax,betay,betaz,Ee1,pex1,pey1,pez1,Tmp_Par);
 793     //     then do particle two
 794     betax = -(px1/E2);
 795     betay = -(py1/E2);
 796     betaz = -(pz1/E2);
 797
 798     transform (betax,betay,betaz,Ee2,pex2,pey2,pez2, Tmp_Par);
 799     //     finally dump the electron values into the approriate
 800     //     variables
 801     E1 = Ee1;
 802     E2 = Ee2;
 803     px1 = pex1;
 804     px2 = pex2;
 805     py1 = pey1;
 806     py2 = pey2;
 807     pz1 = pez1;
 808     pz2 = pez2;
 809 }
 810
 811
 812 //______________________________________________________________________________
 813 double Gammagammaleptonpair::getMass()
 814 {
 815     double leptonmass=0.;
 816     switch(_GGlepInputpidtest){
 817     case starlightConstants::ELECTRON:
 818         leptonmass=starlightConstants::mel;
 819         break;
 820     case starlightConstants::MUON:
 821         leptonmass=starlightConstants::muonMass;
 822         break;
 823     case starlightConstants::TAUON:
 824         leptonmass=starlightConstants::tauMass;
 825         break;
 826     default:
 827         cout<<"Not a recognized lepton, Gammagammaleptonpair::getmass(), mass = 0."<<endl;
 828     }
 829
 830     return leptonmass;
 831 }
 832
 833
 834 //______________________________________________________________________________
 835 double Gammagammaleptonpair::getWidth()
 836 {
 837     double leptonwidth=0.;
 838     return leptonwidth;
 839
 840 }
 841
 842
 843 //______________________________________________________________________________
 844 double Gammagammaleptonpair::getSpin()
 845 {
 846     double leptonspin=0.5;
 847     return leptonspin;
 848 }