--- /dev/null
+///////////////////////////////////////////////////////////////////////////
+//
+// Copyright 2010
+//
+// This file is part of starlight.
+//
+// starlight is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or
+// (at your option) any later version.
+//
+// starlight is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with starlight. If not, see <http://www.gnu.org/licenses/>.
+//
+///////////////////////////////////////////////////////////////////////////
+//
+// File and Version Information:
+// $Rev:: 164 $: revision of last commit
+// $Author:: odjuvsla $: author of last commit
+// $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
+//
+// Description:
+// Nystrand 220710
+// Fixed bug which gave incorrect minv distribution in gammagammaleptonpair.
+// Moved loop over W and Y from pickw to twoLeptonCrossSection in
+// gammagammaleptonpair to speed up event generation.
+// Changed to Woods-Saxon radius in twophotonluminosity to be consistent
+// with old starligt.
+//
+//
+///////////////////////////////////////////////////////////////////////////
+
+
+#include <iostream>
+#include <fstream>
+#include <cmath>
+#include <vector>
+
+#include "starlightconstants.h"
+#include "gammagammaleptonpair.h"
+
+
+using namespace std;
+
+
+//_____________________________________________________________________________
+Gammagammaleptonpair::Gammagammaleptonpair(beamBeamSystem& bbsystem)
+: eventChannel(bbsystem)
+{
+ //Initialize randomgenerator with our seed.
+ cout<<"Randy in leptonpair construction: "<<randyInstance.Rndom()<<endl;
+ //Storing inputparameters into protected members for use
+ _GGlepInputnumw=inputParametersInstance.nmbWBins();
+ _GGlepInputnumy=inputParametersInstance.nmbRapidityBins();
+ _GGlepInputpidtest=inputParametersInstance.prodParticleType();
+ _GGlepInputGamma_em=inputParametersInstance.beamLorentzGamma();
+ //Let us read in the luminosity tables
+ read();
+ //Now we will calculate the crosssection
+ twoLeptonCrossSection();
+ //If it is a tauon, calculate its tables
+ if(inputParametersInstance.prodParticleId()==starlightConstants::TAUON) calculateTable();
+}
+
+
+//______________________________________________________________________________
+Gammagammaleptonpair::~Gammagammaleptonpair()
+{ }
+
+
+//______________________________________________________________________________
+void Gammagammaleptonpair::twoLeptonCrossSection()
+{
+ //This function calculates section for 2-particle decay. For reference, see STAR Note 243, Eq. 9.
+ //calculate the 2-lepton differential cross section
+ //the 100 is to convert to barns
+ //the 2 is to account for the fact that we integrate only over one half of the rapidity range
+ //Multiply all _Farray[] by _f_max
+
+ for(int i=0;i<_GGlepInputnumw;i++)
+ {
+ for(int j=0;j<_GGlepInputnumy;j++)
+ {
+ // _sigmax[i][j]=2.*Gammagammaleptonpair::twoMuonCrossSection(_Warray[i])*_f_max*_Farray[i][j]/100.;
+ _sigmax[i][j]=twoMuonCrossSection(_Warray[i])*_f_max*_Farray[i][j]/(100.*_Warray[i]);
+ }
+ }
+ //calculate the total two-lepton cross section
+ double sigmasum =0.;
+ for(int i=0;i<_GGlepInputnumw-1;i++)
+ {
+ for(int j=0;j<_GGlepInputnumy-1;j++)
+ {
+ // _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.));
+ // _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.));
+ 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]);
+ }
+ }
+ cout << "The total "<<_GGlepInputpidtest<<" cross-section is: "<<sigmasum<<" barns."<<endl;
+
+ // Do this integration here, once per run rather than once per event (JN 220710)
+ //integrate sigma down to a function of just w
+
+ double sgf=0.;
+
+ for(int i=0;i<_ReadInputnumw;i++)
+ {
+ _sigofw[i]=0.;
+ for(int j=0;j<_ReadInputnumy-1;j++)
+ {
+ _sigofw[i] = _sigofw[i]+(_Yarray[j+1]-_Yarray[j])*(_sigmax[i][j+1]+_sigmax[i][j])/2.;
+ }
+ }
+
+ //calculate the unnormalized sgfint array
+ _sigfint[0]=0.;
+ for(int i=0;i<_ReadInputnumw-1;i++)
+ {
+ sgf=(_sigofw[i+1]+_sigofw[i])*(_Warray[i+1]-_Warray[i])/2.;
+ _sigfint[i+1]=_sigfint[i]+sgf;
+ }
+
+ //normalize sgfint array
+ _signormw=_sigfint[_ReadInputnumw-1];
+ for(int i=0;i<_ReadInputnumw;i++)
+ {
+ _sigfint[i]=_sigfint[i]/_signormw;
+ }
+ return;
+}
+
+
+//______________________________________________________________________________
+double Gammagammaleptonpair::twoMuonCrossSection(double w)
+{
+ //This function gives the two muon cross section as a function of Y&W.
+ //Using the formula given in G.Soff et. al Nuclear Equation of State, part B, 579
+ double s=0.,Etest=0.,deltat=0.,xL=0.,sigmuix=0.,alphasquared=0.,hbarcsquared=0.;
+ s = w*w;
+ Etest = 4.*getMass()*getMass()/s;
+ deltat = s * sqrt(1.-Etest);
+ xL = 2.*log(sqrt(s)/(2.*getMass())+sqrt(1./Etest-1.));
+ alphasquared = starlightConstants::alpha*starlightConstants::alpha;
+ hbarcsquared = starlightConstants::hbarc*starlightConstants::hbarc;
+ sigmuix = 4.*starlightConstants::pi*alphasquared/s*hbarcsquared*((1+Etest-0.5*Etest*Etest)*xL-(1./s+Etest/s)*deltat);
+ if(Etest > 1.)
+ sigmuix = 0.;
+ return sigmuix;
+}
+
+
+//______________________________________________________________________________
+void Gammagammaleptonpair::pickw(double &w)
+{
+// This function picks a w for the 2-photon calculation.
+
+ double x=0.,remainarea=0.,remainw=0.,a=0.,b=0.,c=0.;
+ int ivalw=0;
+
+ if(_wdelta != 0)
+ {
+ w=_wdelta;
+ ivalw=_ivalwd;
+ remainw=_remainwd;
+ }
+ else{
+ //deal with the case where sigma is an array
+ //_sigofw is simga integrated over y using a linear interpolation
+ //sigint is the integral of sgfint, normalized
+
+ //pick a random number
+ x = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
+ //compare x and sgfint to find the ivalue which is just less than the random number x
+ for(int i=0;i<_GGlepInputnumw;i++)
+ {
+ if(x > _sigfint[i]) ivalw=i;
+ }
+ //remainder above ivalw
+ remainarea = x - _sigfint[ivalw];
+
+ //figure out what point corresponds to the excess area in remainarea
+ c = -remainarea*_signormw/(_Warray[ivalw+1]-_Warray[ivalw]);
+ b = _sigofw[ivalw];
+ a = (_sigofw[ivalw+1]-_sigofw[ivalw])/2.;
+ if(a==0.)
+ {
+ remainw = -c/b;
+ }
+ else{
+ remainw = (-b+sqrt(b*b-4.*a*c))/(2.*a);
+ }
+ _ivalwd = ivalw;
+ _remainwd = remainw;
+ //calculate the w value
+ w = _Warray[ivalw]+(_Warray[ivalw+1]-_Warray[ivalw])*remainw;
+
+ }
+}
+
+
+//______________________________________________________________________________
+void Gammagammaleptonpair::picky(double &y)
+{
+ // This function picks a y given a W
+
+ double * sigofy;
+ double * sgfint;
+ sigofy = new double[starlightLimits::MAXYBINS];
+ sgfint = new double[starlightLimits::MAXWBINS];
+
+ double remainw =0.,remainarea=0.,remainy=0.,a=0.,b=0.,c=0.,sgf=0.,signorm=0.,x=0.;
+ int ivalw=0,ivaly=0;
+
+ ivalw=_ivalwd;
+ remainw=_remainwd;
+ //average over two colums to get y array
+ for(int j=0;j<_GGlepInputnumy;j++)
+ {
+ sigofy[j]=_sigmax[ivalw][j]+(_sigmax[ivalw+1][j]-_sigmax[ivalw][j])*remainw;
+ }
+
+ //calculate the unnormalized sgfint
+ sgfint[0]=0.;
+ for(int j=0;j<_GGlepInputnumy-1;j++)
+ {
+ sgf = (sigofy[j+1]+sigofy[j])/2.;
+ sgfint[j+1]=sgfint[j]+sgf*(_Yarray[j+1]-_Yarray[j]);
+ }
+
+ //normalize the sgfint array
+ signorm = sgfint[_GGlepInputnumy-1];
+
+ for(int j=0;j<_GGlepInputnumy;j++)
+ {
+ sgfint[j]=sgfint[j]/signorm;
+ }
+
+ //pick a random number
+ x = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
+ //compare x and sgfint to find the ivalue which is just less then the random number x
+ for(int i=0;i<_GGlepInputnumy;i++)
+ {
+ if(x > sgfint[i]) ivaly = i;
+ }
+ //remainder above ivaly
+ remainarea = x - sgfint[ivaly];
+
+ //figure what point corresponds to the leftover area in remainarea
+ c = -remainarea*signorm/(_Yarray[ivaly+1]-_Yarray[ivaly]);
+ b = sigofy[ivaly];
+ a = (sigofy[ivaly+1]-sigofy[ivaly])/2.;
+ if(a==0.)
+ {
+ remainy = -c/b;
+ }
+ else{
+ remainy = (-b + sqrt(b*b-4.*a*c))/(2.*a);
+ }
+ //calculate the y value
+ y = _Yarray[ivaly]+(_Yarray[ivaly+1]-_Yarray[ivaly])*remainy;
+ delete[] sigofy;
+ delete[] sgfint;
+}
+
+
+//______________________________________________________________________________
+void Gammagammaleptonpair::pairMomentum(double w,double y,double &E,double &px,double &py,double &pz)
+{
+ //this function calculates px,py,pz,and E given w and y
+
+ double anglepp1=0.,anglepp2=0.,pp1=0.,pp2=0.,E1=0.,E2=0.,signpx=0.,pt=0.;
+
+ //E1 and E2 are for the 2 photons in the CM frame
+ E1 = w*exp(y)/2.;
+ E2 = w*exp(-y)/2.;
+
+ //calculate px and py
+ //to get x and y components-- phi is random between 0 and 2*pi
+ anglepp1 = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
+ anglepp2 = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
+
+ pp1 = pp_1(E1);
+ pp2 = pp_2(E2);
+ px = pp1*cos(2.*starlightConstants::pi*anglepp1)+pp2*cos(2.*starlightConstants::pi*anglepp2);
+ py = pp1*sin(2.*starlightConstants::pi*anglepp1)+pp2*sin(2.*starlightConstants::pi*anglepp2);
+
+ //Compute vector sum Pt=Pt1+Pt2 to find pt for the produced particle
+ pt = sqrt(px*px+py*py);
+
+ //W is the mass of the produced particle (not necessarily on-mass-shell).Now compute its energy and pz
+ E = sqrt(w*w+pt*pt)*cosh(y);
+ pz= sqrt(w*w+pt*pt)*sinh(y);
+ signpx = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
+
+ //pick the z direction
+ //Don't do this anymore since y goes from -ymax to +ymax (JN 15-02-2013)
+ //if(signpx > 0.5) pz = -pz;
+}
+
+
+//______________________________________________________________________________
+double Gammagammaleptonpair::pp_1(double E)
+{
+ // This is for beam 1
+ // returns on random draw from pp(E) distribution
+ double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
+ double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
+ int satisfy =0;
+
+ ereds = (E/_GGlepInputGamma_em)*(E/_GGlepInputGamma_em);
+ //sqrt(3)*E/gamma_em is p_t where the distribution is a maximum
+ Cm = sqrt(3.)*E/_GGlepInputGamma_em;
+ //the amplitude of the p_t spectrum at the maximum
+ singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
+ Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
+
+ //pick a test value pp, and find the amplitude there
+ x = randyInstance.Rndom();
+ pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
+ singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
+ test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
+
+ while(satisfy==0){
+ u = randyInstance.Rndom();
+ if(u*Coef <= test)
+ {
+ satisfy =1;
+ }
+ else{
+ x =randyInstance.Rndom();
+ pp = 5*starlightConstants::hbarc/_bbs.beam1().nuclearRadius()*x;
+ singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
+ test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
+ }
+ }
+
+ return pp;
+}
+
+//______________________________________________________________________________
+double Gammagammaleptonpair::pp_2(double E)
+{
+
+ // This is for beam 2
+ //returns on random draw from pp(E) distribution
+ double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
+ double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
+ int satisfy =0;
+
+ ereds = (E/_GGlepInputGamma_em)*(E/_GGlepInputGamma_em);
+ //sqrt(3)*E/gamma_em is p_t where the distribution is a maximum
+ Cm = sqrt(3.)*E/_GGlepInputGamma_em;
+ //the amplitude of the p_t spectrum at the maximum
+ singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
+ Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
+
+ //pick a test value pp, and find the amplitude there
+ x = randyInstance.Rndom();
+ pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius(); //Will use nucleus #1
+ singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
+ test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
+
+ while(satisfy==0){
+ u = randyInstance.Rndom();
+ if(u*Coef <= test)
+ {
+ satisfy =1;
+ }
+ else{
+ x =randyInstance.Rndom();
+ pp = 5*starlightConstants::hbarc/_bbs.beam2().nuclearRadius()*x;
+ singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds);
+ test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
+ }
+ }
+
+ return pp;
+}
+
+
+
+//______________________________________________________________________________
+void Gammagammaleptonpair::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
+ double , // E (unused)
+ double W,
+ double px0, double py0, double pz0,
+ double& px1, double& py1, double& pz1,
+ double& px2, double& py2, double& pz2,
+ int& iFbadevent)
+{
+ // This routine decays a particle into two particles of mass mdec,
+ // taking spin into account
+
+ double mdec=0.,E1=0.,E2=0.;
+ double pmag, anglelep[20001];
+ // double ytest=0.,dndtheta;
+ double phi,theta,xtest,Ecm;
+ double betax,betay,betaz;
+ double hirestheta,hirestest,hiresw; //added from JN->needed precision
+
+ // set the mass of the daughter particles
+
+ mdec = getMass();
+ if(W < 2*mdec)
+ {
+ cout<<" ERROR: W="<<W<<endl;
+ iFbadevent = 1;
+ return;
+ }
+ pmag = sqrt(W*W/4. - mdec*mdec);
+
+ // pick an orientation, based on the spin
+ // phi has a flat distribution in 2*pi
+ phi = randyInstance.Rndom()*2.*starlightConstants::pi; //(random()/(RAND_MAX+1.0))* 2.*starlightConstants::pi;
+
+ // find theta, the angle between one of the outgoing particles and
+ // the beamline, in the frame of the two photons
+
+ if(getSpin()==0.5){
+ // calculate a table of integrated angle values for leptons
+ // JN05: Go from 1000->20000bins, use double precision for anglelep and thetalep. needed when W >6Gev.
+ hiresw = W;
+
+ anglelep[0] = 0.;
+
+ for(int i =1;i<=20000;i++)
+ {
+ hirestheta = starlightConstants::pi * double(i) /20000.;
+
+ // Added sin(theta) phase space factor (not in thetalep) and changed E to W in thetalep call
+ // 11/9/2000 SRK
+ // Note that thetalep is form invariant, so it could be called for E, theta_lab just
+ // as well as W,theta_cm. Since there is a boost from cm to lab below, the former is fine.
+
+ anglelep[i] = anglelep[i-1] + thetalep(hiresw,hirestheta)*sin(hirestheta);
+ }
+
+ hirestheta = 0.;
+ xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
+ hirestest = xtest;
+ for(int i =1;i<=20000;i++)
+ {
+ if(xtest > (anglelep[i]/anglelep[20000]))
+ hirestheta = starlightConstants::pi * double(i) / 20000.;
+ }
+ theta=hirestheta;
+
+ }
+
+ if(getSpin() != 0.5)
+ cout<<" This model cannot yet handle this spin value for lepton pairs: "<<getSpin()<<endl;
+
+
+ // compute unboosted momenta
+ px1 = sin(theta)*cos(phi)*pmag;
+ py1 = sin(theta)*sin(phi)*pmag;
+ pz1 = cos(theta)*pmag;
+ px2 = -px1;
+ py2 = -py1;
+ pz2 = -pz1;
+
+ // compute energies
+ //Changed mass to W 11/9/2000 SRK
+ Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0);
+
+ E1 = sqrt(mdec*mdec+px1*px1+py1*py1+pz1*pz1);
+ E2 = sqrt(mdec*mdec+px2*px2+py2*py2+pz2*pz2);
+ // decay tau to electrons
+ // note that after this routine px1, etc., refer to the electrons
+ if(_GGlepInputpidtest == starlightConstants::TAUON)
+ tauDecay(px1,py1,pz1,E1,px2,py2,pz2,E2);
+
+ // Lorentz transform into the lab frame
+ // betax,betay,betaz are the boost of the complete system
+ betax = -(px0/Ecm);
+ betay = -(py0/Ecm);
+ betaz = -(pz0/Ecm);
+
+ transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent);
+ transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent);
+
+
+ if(iFbadevent == 1)
+ return;
+
+ // change particle id from that of parent to that of daughters
+ // change taoun id into electron id, already switched particles in tau decay
+ if(_GGlepInputpidtest == starlightConstants::TAUON)
+ ipid = starlightConstants::ELECTRON;
+ // electrons remain electrons; muons remain muons
+ if ((_GGlepInputpidtest == starlightConstants::ELECTRON) || (_GGlepInputpidtest == starlightConstants::MUON))
+ ipid = _GGlepInputpidtest;
+}
+
+
+//______________________________________________________________________________
+double Gammagammaleptonpair::thetalep(double W,double theta)
+{
+ // This function calculates the cross-section as a function of
+ // angle for a given W and Y, for the production of two muons.
+ // (or tauons)
+ // expression is taken from Brodsky et al. PRD 1971, 1532
+ // equation 5.7
+ // factor that are not dependant on theta are scrapped, so the
+ // the absolute crosssections given by this function are inaccurate
+ // here we are working in the CM frame of the photons and the last
+ // term is 0
+
+ // real function thetalep (W,theta)
+
+ double moverw=0., W1sq=0.;
+ double thetalep_r=0.,denominator=0.;
+
+ W1sq = (W / 2.)*(W/2.);
+ moverw = getMass()*getMass() / W1sq;
+ denominator = (1. - (1. - moverw)*(cos(theta)*cos(theta)));
+
+ thetalep_r = 2. + 4.*(1.-moverw)*((1.-moverw)*sin(theta)*sin(theta)*cos(theta)*cos(theta) + moverw) / (denominator*denominator);
+
+ return thetalep_r;
+
+}
+
+
+//______________________________________________________________________________
+starlightConstants::event Gammagammaleptonpair::produceEvent(int &ievent)
+{//returns the vector with the decay particles inside.
+ starlightConstants::event leptonpair; //This object will store all the tracks for a single event
+ double comenergy = 0.;
+ double rapidity = 0.;
+ double pairE = 0.;
+ double pairmomx=0.,pairmomy=0.,pairmomz=0.;
+ int iFbadevent=0;
+ starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
+
+ double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
+//this function decays particles and writes events to a file
+ //zero out the event structure
+ leptonpair._numberOfTracks=0;
+ for(int i=0;i<4;i++)
+ {
+ leptonpair.px[i]=0.;
+ leptonpair.py[i]=0.;
+ leptonpair.pz[i]=0.;
+ leptonpair._fsParticle[i]=starlightConstants::UNKNOWN;
+ leptonpair._charge[i]=0;
+ }
+
+ pickw(comenergy);
+
+ picky(rapidity);
+
+ pairMomentum(comenergy,rapidity,pairE,pairmomx,pairmomy,pairmomz);
+ twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
+ if (iFbadevent==0){
+ int q1=0,q2=0;
+
+ double xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
+ if (xtest<0.5)
+ {
+ q1=1;
+ q2=-1;
+ }
+ else{
+ q1=-1;
+ q2=1;
+ }
+ leptonpair._numberOfTracks=2;//leptonpairs are two tracks...
+ leptonpair.px[0]=px1;
+ leptonpair.py[0]=py1;
+ leptonpair.pz[0]=pz1;
+ leptonpair._fsParticle[0]=ipid;
+ leptonpair._charge[0]=q1;
+
+ leptonpair.px[1]=px2;
+ leptonpair.py[1]=py2;
+ leptonpair.pz[1]=pz2;
+ leptonpair._fsParticle[1]=ipid;
+ leptonpair._charge[1]=q2;
+
+ ievent=ievent+1;
+ }
+
+ return leptonpair;
+}
+
+
+//______________________________________________________________________________
+upcEvent Gammagammaleptonpair::produceEvent()
+{
+//returns the vector with the decay particles inside.
+
+ upcEvent event;
+
+ double comenergy = 0.;
+ double rapidity = 0.;
+ double pairE = 0.;
+ double pairmomx=0.,pairmomy=0.,pairmomz=0.;
+ int iFbadevent=0;
+ starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
+
+ double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.,E2=0.,E1=0.;
+ bool accepted = false;
+ do{
+ //this function decays particles and writes events to a file
+ //zero out the event structure
+ pickw(comenergy);
+
+ picky(rapidity);
+
+ pairMomentum(comenergy,rapidity,pairE,pairmomx,pairmomy,pairmomz);
+
+
+ _nmbAttempts++;
+ twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
+ double pt1chk = sqrt(px1*px1+py1*py1);
+ double pt2chk = sqrt(px2*px2+py2*py2);
+
+ double eta1 = pseudoRapidity(px1, py1, pz1);
+ double eta2 = pseudoRapidity(px2, py2, pz2);
+
+ if(_ptCutEnabled && !_etaCutEnabled){
+ if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){
+ accepted = true;
+ _nmbAccepted++;
+ }
+ }
+ else if(!_ptCutEnabled && _etaCutEnabled){
+ if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
+ accepted = true;
+ _nmbAccepted++;
+ }
+ }
+ else if(_ptCutEnabled && _etaCutEnabled){
+ if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){
+ if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
+ accepted = true;
+ _nmbAccepted++;
+ }
+ }
+ }
+
+ }while((_ptCutEnabled || _etaCutEnabled) && !accepted);
+ //twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
+
+ if (iFbadevent==0){
+ int q1=0,q2=0;
+
+ double xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
+ if (xtest<0.5)
+ {
+ q1=1;
+ q2=-1;
+ }
+ else{
+ q1=-1;
+ q2=1;
+ }
+
+ // The new stuff
+ double mlepton = getMass();
+ E1 = sqrt( mlepton*mlepton + px1*px1 + py1*py1 + pz1*pz1 );
+ E2 = sqrt( mlepton*mlepton + px2*px2 + py2*py2 + pz2*pz2 );
+
+ starlightParticle particle1(px1, py1, pz1, E1, starlightConstants::UNKNOWN, -q1*ipid, q1);
+ event.addParticle(particle1);
+
+ starlightParticle particle2(px2, py2, pz2, E2, starlightConstants::UNKNOWN, -q2*ipid, q2);
+ event.addParticle(particle2);
+
+ }
+ return event;
+}
+
+
+//______________________________________________________________________________
+void Gammagammaleptonpair::calculateTable()
+{
+ // this subroutine calculates the tables that are used
+ // elsewhere in the montecarlo
+ // the tauon decay is taken from V-A theory, 1 - 1/3 cos(theta)
+ // the energy of each of the two leptons in tau decay
+ // is calculated using formula 10.35 given
+ // in introduction to elementary particles by D. griffiths
+ // which assmes that the mass of the electron is 0.
+ // the highest energy attainable by the electron in such a system is
+ // .5 * mass of the tau
+
+ // subroutine calculateTable
+
+
+ double E,theta;
+
+ _tautolangle[0] = 0.;
+ _dgammade[0] = 0.;
+
+
+ for(int i =1;i<=100;i++)
+ {
+ // calculate energy of tau decay
+ E = double(i)/100. * .5 * starlightConstants::tauMass;
+ _dgammade[i] = _dgammade[i-1] + E*E * (1. - 4.*E/(3.*starlightConstants::tauMass));
+
+ // calculate angles for tau
+ theta = starlightConstants::pi * double(i) / 100.;
+ _tautolangle[i] = _tautolangle[i-1] + (1 + 0.3333333 * cos(theta));
+ }
+
+
+}
+
+
+//______________________________________________________________________________
+void Gammagammaleptonpair::tauDecay(double &px1,double &py1,double &pz1,double &E1,double &px2,double &py2,double &pz2,double &E2)
+{
+ // this routine assumes that the tauons decay to electrons and
+ // calculates the directons of the decays
+
+ double Ee1,Ee2,theta1,theta2,phi1,phi2, ran1, ran2 ;
+ double pmag1,pex1,pey1,pez1,pex2,pey2,pez2,pmag2;
+ double betax,betay,betaz,dir;
+
+ int Tmp_Par=0; // temp variable for the transform function .. kind of starnge - being called with 7 parameter instead of 8
+
+ // the highest energy attainable by the electron in this system is
+ // .5 * mass of the tau
+
+ // get two random numbers to compare with
+
+
+ ran1 = randyInstance.Rndom()*_dgammade[100];//(random()/(RAND_MAX+1.0)) * _dgammade[100];
+ ran2 = randyInstance.Rndom()*_dgammade[100];//(random()/(RAND_MAX+1.0)) * _dgammade[100];
+
+ // compute the energies that correspond to those numbers
+ Ee1 = 0.;
+ Ee2 = 0.;
+
+ for( int i =1;i<=100;i++)
+ {
+ if (ran1 > _dgammade[i])
+ Ee1 = double(i) /100. * .5 * getMass();
+ if (ran2 > _dgammade[i])
+ Ee2 = double(i) /100. * .5 * getMass();
+ }
+
+ // to find the direction of emmission, first
+ // we determine if the tauons have spin of +1 or -1 along the
+ // direction of the beam line
+ dir = 1.;
+ if ( randyInstance.Rndom() < 0.5 )//(random()/(RAND_MAX+1.0)) < 0.5)
+ dir = -1.;
+
+ // get two random numbers to compare with
+ ran1 = randyInstance.Rndom()*_tautolangle[99];//(random()/(RAND_MAX+1.0)) * _tautolangle[100];
+ ran2 = randyInstance.Rndom()*_tautolangle[99];//(random()/(RAND_MAX+1.0)) * _tautolangle[100];
+
+ // find the angles corrsponding to those numbers
+ theta1 = 0.;
+ theta2 = 0.;
+ for( int i =1;i<=100;i++)
+ {
+ if (ran1 > _tautolangle[i]) theta1 = starlightConstants::pi * double(i) /100.;
+ if (ran2 > _tautolangle[i]) theta2 = starlightConstants::pi * double(i) /100.;
+ }
+
+ // grab another two random numbers to determine phi's
+ phi1 = randyInstance.Rndom()*2.*starlightConstants::pi;// (random()/(RAND_MAX+1.0))* 2. * starlightConstants::pi;
+ phi2 = randyInstance.Rndom()*2.*starlightConstants::pi;// (random()/(RAND_MAX+1.0))* 2. * starlightConstants::pi;
+ // figure out the momenta of the electron in the frames of the
+ // tauons from which they decayed, that is electron1 is in the
+ // rest frame of tauon1 and e2 is in the rest fram of tau2
+ // again the electrons are assumed to be massless
+ pmag1 = Ee1;
+ pex1 = cos(phi1)*sin(theta1)*pmag1;
+ pey1 = sin(phi1)*sin(theta1)*pmag1;
+ pez1 = cos(theta1)*pmag1*dir;
+ pmag2 = Ee2;
+ pex2 = cos(phi2)*sin(theta2)*pmag2;
+ pey2 = sin(phi2)*sin(theta2)*pmag2;
+ pez2 = cos(theta2)*pmag2*(-1.*dir);
+ // now Lorentz transform into the frame of each of the particles
+ // do particle one first
+ betax = -(px1/E1);
+ betay = -(py1/E1);
+ betaz = -(pz1/E1);
+//cout<<"2decay betax,pex1= "<<betax<<" "<<pex1<<endl;
+ transform (betax,betay,betaz,Ee1,pex1,pey1,pez1,Tmp_Par);
+ // then do particle two
+ betax = -(px1/E2);
+ betay = -(py1/E2);
+ betaz = -(pz1/E2);
+
+ transform (betax,betay,betaz,Ee2,pex2,pey2,pez2, Tmp_Par);
+ // finally dump the electron values into the approriate
+ // variables
+ E1 = Ee1;
+ E2 = Ee2;
+ px1 = pex1;
+ px2 = pex2;
+ py1 = pey1;
+ py2 = pey2;
+ pz1 = pez1;
+ pz2 = pez2;
+}
+
+
+//______________________________________________________________________________
+double Gammagammaleptonpair::getMass()
+{
+ double leptonmass=0.;
+ switch(_GGlepInputpidtest){
+ case starlightConstants::ELECTRON:
+ leptonmass=starlightConstants::mel;
+ break;
+ case starlightConstants::MUON:
+ leptonmass=starlightConstants::muonMass;
+ break;
+ case starlightConstants::TAUON:
+ leptonmass=starlightConstants::tauMass;
+ break;
+ default:
+ cout<<"Not a recognized lepton, Gammagammaleptonpair::getmass(), mass = 0."<<endl;
+ }
+
+ return leptonmass;
+}
+
+
+//______________________________________________________________________________
+double Gammagammaleptonpair::getWidth()
+{
+ double leptonwidth=0.;
+ return leptonwidth;
+
+}
+
+
+//______________________________________________________________________________
+double Gammagammaleptonpair::getSpin()
+{
+ double leptonspin=0.5;
+ return leptonspin;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff