+++ /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:: $: revision of last commit
-// $Author:: $: author of last commit
-// $Date:: $: 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