estarlight/src/gammagammasingle.cpp

0001 ///////////////////////////////////////////////////////////////////////////
0002 //
0003 //    Copyright 2010
0004 //
0005 //    This file is part of starlight.
0006 //
0007 //    starlight is free software: you can redistribute it and/or modify
0008 //    it under the terms of the GNU General Public License as published by
0009 //    the Free Software Foundation, either version 3 of the License, or
0010 //    (at your option) any later version.
0011 //
0012 //    starlight is distributed in the hope that it will be useful,
0013 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0015 //    GNU General Public License for more details.
0016 //
0017 //    You should have received a copy of the GNU General Public License
0018 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
0019 //
0020 ///////////////////////////////////////////////////////////////////////////
0021 //
0022 // File and Version Information:
0023 // $Rev:: 274                         $: revision of last commit
0024 // $Author:: butter                   $: author of last commit
0025 // $Date:: 2016-09-11 23:40:25 +0100 #$: date of last commit
0026 //
0027 // Description:
0028 //
0029 //
0030 //
0031 ///////////////////////////////////////////////////////////////////////////
0032
0033 #include <iostream>
0034 #include <fstream>
0035 #include <cmath>
0036 #include <vector>
0037
0038
0039 #include "starlightconstants.h"
0040 #include "gammagammasingle.h"
0041 #include "starlightconfig.h"
0042
0043 using namespace std;
0044
0045
0046 //______________________________________________________________________________
0047 Gammagammasingle::Gammagammasingle(const inputParameters& inputParametersInstance, beamBeamSystem& bbsystem)
0048 : eventChannel(inputParametersInstance, bbsystem)
0049 #ifdef ENABLE_PYTHIA
0050 ,_pyDecayer()
0051 #endif
0052 {
0053
0054 #ifdef ENABLE_PYTHIA
0055     _pyDecayer.init();
0056 #endif
0057
0058   //Storing inputparameters into protected members for use
0059   _GGsingInputnumw=inputParametersInstance.nmbWBins();
0060   _GGsingInputnumy=inputParametersInstance.nmbRapidityBins();
0061   _GGsingInputpidtest=inputParametersInstance.prodParticleType();
0062   _GGsingInputGamma_em=inputParametersInstance.beamLorentzGamma();
0063   _axionMass=inputParametersInstance.axionMass(); // AXION HACK
0064   cout<<"SINGLE MESON pid test: "<<_GGsingInputpidtest<<endl;
0065   //reading in luminosity tables
0066   read();
0067   //Now calculating crosssection
0068   singleCrossSection();
0069 }
0070
0071
0072 //______________________________________________________________________________
0073 Gammagammasingle::~Gammagammasingle()
0074 { }
0075
0076
0077 //______________________________________________________________________________
0078 void Gammagammasingle::singleCrossSection()
0079 {
0080   //This function carries out a delta function cross-section calculation. For reference, see STAR Note 243, Eq. 8
0081   //Multiply all _Farray[] by _f_max
0082   double _sigmaSum=0.,remainw=0.;//_remainwd=0.;
0083   int ivalw =0;//_ivalwd;
0084   //calculate the differential cross section and place in the sigma table
0085   cout << "MASS  " << getMass() << "\n"; // AXION HACK, optional
0086   cout << "WIDTH  " << getWidth() << "\n";// AXION HACK, optional
0087   _wdelta=getMass();
0088   for(int i=0;i<_GGsingInputnumw;i++){
0089     for(int j=0;j<_GGsingInputnumy;j++){
0090       // Eq. 1 of starnote 347
0091       _sigmax[i][j]=(getSpin()*2.+1.)*4*starlightConstants::pi*starlightConstants::pi*getWidth()/
0092     (getMass()*getMass()*getMass())*_f_max*_Farray[i][j]*starlightConstants::hbarc*starlightConstants::hbarc/100.;
0093     }
0094   }
0095   //find the index, i,for the value of w just less than the mass because we want to use the value from the sigma table that has w=mass
0096
0097   for(int i=0;i<_GGsingInputnumw;i++){
0098     if(getMass()>_Warray[i]) ivalw=i;
0099   }
0100
0101   remainw = (getMass()-_Warray[ivalw])/(_Warray[ivalw+1]-_Warray[ivalw]);
0102   _ivalwd = ivalw;
0103   _remainwd = remainw;
0104   //if we are interested rho pairs at threshold, the just set sigma to 100nb
0105   switch(_GGsingInputpidtest){
0106   case starlightConstants::ZOVERZ03:
0107     _sigmaSum =0.;
0108     for(int j=0;j<_GGsingInputnumy-1;j++){
0109                         _sigmaSum = _sigmaSum +(_Yarray[j+1]-_Yarray[j])*
0110               100.0E-9*(.1/getMass())*((1.-remainw)*_f_max*
0111                            (_Farray[ivalw][j]+_Farray[ivalw][j])/2.+remainw*_f_max*
0112                            (_Farray[ivalw+1][j]+_Farray[ivalw+1][j+1])/2.);
0113     }
0114     break;
0115   default:
0116     //Sum to find the total cross-section
0117     _sigmaSum =0.;
0118     for(int j =0;j<_GGsingInputnumy-1;j++){
0119                         _sigmaSum = _sigmaSum+
0120               (_Yarray[j+1]-_Yarray[j])*((1.-remainw)*
0121                            (_sigmax[ivalw][j]+_sigmax[ivalw][j+1])/2.+remainw*
0122                            (_sigmax[ivalw+1][j]+_sigmax[ivalw+1][j+1])/2.);
0123     }
0124   }
0125   // if(_sigmaSum > 0.1) cout <<"The total cross-section is: "<<_sigmaSum<<" barns."<<endl;
0126   // else if(_sigmaSum > 0.0001)cout <<"The total cross-section is: "<<_sigmaSum*1000<<" mb."<<endl;
0127   // else cout <<"The total cross-section is: "<<_sigmaSum*1000000<<" ub."<<endl;
0128   cout<<endl;
0129   if (_sigmaSum > 1.){
0130      cout << "Total cross section: "<<_sigmaSum<<" barn."<<endl;
0131   } else if (1000.*_sigmaSum > 1.){
0132      cout << "Total cross section: "<<1000.*_sigmaSum<<" mb."<<endl;
0133   } else if (1000000.*_sigmaSum > 1.){
0134     cout << "Total cross section: "<<1000000.*_sigmaSum<<" microbarn."<<endl;
0135   } else if (1.E9*_sigmaSum > 1.){
0136     cout << "Total cross section: "<<1.E9*_sigmaSum<<" nanobarn."<<endl;
0137   } else if (1.E12*_sigmaSum > 1.){
0138     cout << "Total cross section: "<<1.E12*_sigmaSum<<" picobarn."<<endl;
0139   } else {
0140     cout << "Total cross section: "<<1.E15*_sigmaSum<<" femtobarn."<<endl;
0141   }
0142   cout<<endl;
0143   setTotalChannelCrossSection(_sigmaSum);
0144
0145   return;
0146 }
0147
0148
0149 //______________________________________________________________________________
0150 void Gammagammasingle::pickw(double &w)
0151 {
0152   //This function picks a w for the 2-photon calculation.
0153   double sgf=0.,signorm=0.,x=0.,remainarea=0.,remainw=0.,a=0.,b=0.,c=0.;
0154   int ivalw=0;
0155
0156   double * _sigofw;
0157   double * sgfint;
0158   _sigofw = new double[starlightLimits::MAXWBINS];
0159   sgfint = new double[starlightLimits::MAXYBINS];
0160
0161   if(_wdelta != 0){
0162     w=_wdelta;
0163     ivalw=_ivalwd;
0164     remainw=_remainwd;
0165   }
0166   else{
0167     //deal with the case where sigma is an array
0168     //_sigofw is simga integrated over y using a linear interpolation
0169     //sigint is the integral of sgfint, normalized
0170
0171     //integrate sigma down to a function of just w
0172     for(int i=0;i<_GGsingInputnumw;i++){
0173       _sigofw[i]=0.;
0174       for(int j=0;j<_GGsingInputnumy-1;j++){
0175     _sigofw[i] = _sigofw[i]+(_Yarray[j+1]-_Yarray[j])*(_sigmax[i][j+1]+_sigmax[i][j])/2.;
0176       }
0177     }
0178     //calculate the unnormalized sgfint array
0179     sgfint[0]=0.;
0180     for(int i=0;i<_GGsingInputnumw-1;i++){
0181       sgf=(_sigofw[i+1]+_sigofw[i])*(_Warray[i+1]-_Warray[i])/2.;
0182       sgfint[i+1]=sgfint[i]+sgf;
0183     }
0184     //normalize sgfint array
0185     signorm=sgfint[_GGsingInputnumw-1];
0186
0187     for(int i=0;i<_GGsingInputnumw;i++){
0188       sgfint[i]=sgfint[i]/signorm;
0189     }
0190     //pick a random number
0191     x = _randy.Rndom();
0192     //compare x and sgfint to find the ivalue which is just less than the random number x
0193     for(int i=0;i<_GGsingInputnumw;i++){
0194       if(x > sgfint[i]) ivalw=i;
0195     }
0196     //remainder above ivalw
0197     remainarea = x - sgfint[ivalw];
0198
0199     //figure out what point corresponds to the excess area in remainarea
0200     c = -remainarea*signorm/(_Warray[ivalw+1]-_Warray[ivalw]);
0201     b = _sigofw[ivalw];
0202     a = (_sigofw[ivalw+1]-_sigofw[ivalw])/2.;
0203     if(a==0.){
0204       remainw = -c/b;
0205     }
0206     else{
0207       remainw = (-b+sqrt(b*b-4.*a*c))/(2.*a);
0208     }
0209     _ivalwd = ivalw;
0210     _remainwd = remainw;
0211     //calculate the w value
0212     w = _Warray[ivalw]+(_Warray[ivalw+1]-_Warray[ivalw])*remainw;
0213     }
0214
0215   delete[] _sigofw;
0216   delete[] sgfint;
0217 }
0218
0219
0220 //______________________________________________________________________________
0221 void Gammagammasingle::picky(double &y)
0222 {
0223   double * sigofy;
0224   double * sgfint;
0225   sigofy = new double[starlightLimits::MAXYBINS];
0226   sgfint = new double[starlightLimits::MAXYBINS];
0227
0228   double remainw =0.,remainarea=0.,remainy=0.,a=0.,b=0.,c=0.,sgf=0.,signorm=0.,x=0.;
0229   int ivalw=0,ivaly=0;
0230
0231   ivalw=_ivalwd;
0232   remainw=_remainwd;
0233   //average over two colums to get y array
0234   for(int j=0;j<_GGsingInputnumy;j++){
0235     sigofy[j]=_sigmax[ivalw][j]+(_sigmax[ivalw+1][j]-_sigmax[ivalw][j])*remainw;
0236   }
0237   //calculate the unnormalized sgfint
0238
0239   sgfint[0]=0.;
0240   for(int j=0;j<_GGsingInputnumy-1;j++){
0241     sgf = (sigofy[j+1]+sigofy[j])/2.;
0242     sgfint[j+1]=sgfint[j]+sgf*(_Yarray[j+1]-_Yarray[j]);
0243   }
0244
0245   //normalize the sgfint array
0246   signorm = sgfint[_GGsingInputnumy-1];
0247
0248   for(int j=0;j<_GGsingInputnumy;j++){
0249     sgfint[j]=sgfint[j]/signorm;
0250   }
0251   //pick a random number
0252   x = _randy.Rndom();
0253   //compare x and sgfint to find the ivalue which is just less then the random number x
0254   for(int i=0;i<_GGsingInputnumy;i++){
0255     if(x > sgfint[i])
0256       ivaly = i;
0257   }
0258   //remainder above ivaly
0259   remainarea = x - sgfint[ivaly];
0260   //figure what point corresponds to the leftover area in remainarea
0261   c = -remainarea*signorm/(_Yarray[ivaly+1]-_Yarray[ivaly]);
0262   b = sigofy[ivaly];
0263   a = (sigofy[ivaly+1]-sigofy[ivaly])/2.;
0264   if(a==0.){
0265     remainy = -c/b;
0266   }
0267   else{
0268     remainy = (-b + sqrt(b*b-4.*a*c))/(2.*a);
0269   }
0270   //calculate the y value
0271   y = _Yarray[ivaly]+(_Yarray[ivaly+1]-_Yarray[ivaly])*remainy;
0272   delete[] sigofy;
0273   delete[] sgfint;
0274 }
0275
0276
0277 //______________________________________________________________________________
0278 void Gammagammasingle::parentMomentum(double w,double y,double &E,double &px,double &py,double &pz)
0279 {
0280   //this function calculates px,py,pz,and E given w and y
0281   double anglepp1=0.,anglepp2=0.,ppp1=0.,ppp2=0.,E1=0.,E2=0.,signpx=0.,pt=0.;
0282
0283   //E1 and E2 are for the 2 photons in the CM frame
0284   E1 = w*exp(y)/2.;
0285   E2 = w*exp(-y)/2.;
0286   //calculate px and py
0287   //to get x and y components-- phi is random between 0 and 2*pi
0288   anglepp1 = _randy.Rndom();
0289   anglepp2 = _randy.Rndom();
0290
0291   ppp1 = pp1(E1);
0292   ppp2 = pp2(E2);
0293   px = ppp1*cos(2.*starlightConstants::pi*anglepp1)+ppp2*cos(2.*starlightConstants::pi*anglepp2);
0294   py = ppp1*sin(2.*starlightConstants::pi*anglepp1)+ppp2*sin(2.*starlightConstants::pi*anglepp2);
0295   //Compute vector sum Pt=Pt1+Pt2 to find pt for the produced particle
0296   pt = sqrt(px*px+py*py);
0297   //W is the mass of the produced particle (not necessarily on-mass-shell).Now compute its energy and pz
0298   E = sqrt(w*w+pt*pt)*cosh(y);
0299   pz= sqrt(w*w+pt*pt)*sinh(y);
0300   signpx = _randy.Rndom();
0301   //pick the z direction
0302   if(signpx > 0.5)
0303     pz = -pz;
0304 }
0305
0306
0307 //______________________________________________________________________________
0308 double Gammagammasingle::pp1(double E)
0309 {
0310   // First 'copy' of pp, for nucleus 1 form factor.  The split was needed to handle asymmetric beams.  SRK 4/2015
0311   //  will probably have to pass in beambeamsys? that way we can do beam1.formFactor(t) or beam2..., careful with the way sergey did it for asymmetry
0312   //  returns on random draw from pp(E) distribution
0313
0314   double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
0315   double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
0316   int satisfy =0;
0317
0318   ereds = (E/_GGsingInputGamma_em)*(E/_GGsingInputGamma_em);
0319   Cm = sqrt(3.)*E/_GGsingInputGamma_em;
0320   //the amplitude of the p_t spectrum at the maximum
0321   singleformfactorCm=_bbs.electronBeam().formFactor(Cm*Cm+ereds);
0322   //Doing this once and then storing it as a double, for the beam 1 form factor.
0323   Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
0324
0325   //pick a test value pp, and find the amplitude there
0326   x = _randy.Rndom();
0327   pp = x*5.*starlightConstants::hbarc/_bbs.electronBeam().nuclearRadius(); //Will use nucleus #1, there should be two for symmetry//nextline
0328   singleformfactorpp1=_bbs.electronBeam().formFactor(pp*pp+ereds);
0329   test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
0330
0331   while(satisfy==0){
0332     u = _randy.Rndom();
0333     if(u*Coef <= test){
0334       satisfy =1;
0335     }
0336     else{
0337       x =_randy.Rndom();
0338       pp = 5*starlightConstants::hbarc/_bbs.electronBeam().nuclearRadius()*x;
0339       singleformfactorpp2=_bbs.electronBeam().formFactor(pp*pp+ereds);//Symmetry
0340       test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
0341     }
0342   }
0343   return pp;
0344 }
0345
0346 //______________________________________________________________________________
0347 double Gammagammasingle::pp2(double E)
0348 {
0349   // Second 'copy' of pp, for nucleus 1 form factor.  The split was needed to handle asymmetric beams.  SRK 4/2015
0350   //  will probably have to pass in beambeamsys? that way we can do electronBeam.formFactor(t) or beam2..., careful with the way sergey did it for asymmetry
0351   //  returns on random draw from pp(E) distribution
0352
0353   double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
0354   double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
0355   int satisfy =0;
0356
0357   ereds = (E/_GGsingInputGamma_em)*(E/_GGsingInputGamma_em);
0358   Cm = sqrt(3.)*E/_GGsingInputGamma_em;
0359   //the amplitude of the p_t spectrum at the maximum
0360   singleformfactorCm=_bbs.targetBeam().formFactor(Cm*Cm+ereds);
0361   //Doing this once and then storing it as a double, which we square later...SYMMETRY?using electronBeam for now.
0362   Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
0363
0364   //pick a test value pp, and find the amplitude there
0365   x = _randy.Rndom();
0366   pp = x*5.*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius(); //Will use nucleus #1, there should be two for symmetry//nextline
0367   singleformfactorpp1=_bbs.targetBeam().formFactor(pp*pp+ereds);
0368   test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
0369
0370   while(satisfy==0){
0371     u = _randy.Rndom();
0372     if(u*Coef <= test){
0373       satisfy =1;
0374     }
0375     else{
0376       x =_randy.Rndom();
0377       pp = 5*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius()*x;
0378       singleformfactorpp2=_bbs.targetBeam().formFactor(pp*pp+ereds);//Symmetry
0379       test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
0380     }
0381   }
0382   return pp;
0383 }
0384
0385
0386 //______________________________________________________________________________
0387 void Gammagammasingle::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,double W,double px0,double py0,double pz0,double &px1,double &py1,double &pz1,double &px2,double &py2,double &pz2,int &iFbadevent)
0388 {
0389   //     This routine decays a particle into two particles of mass mdec,
0390   //     taking spin into account
0391
0392   double mdec=0.,E1=0.,E2=0.;
0393   double pmag,ytest=0.;
0394   double phi,theta,xtest,dndtheta,Ecm;
0395   double  betax,betay,betaz;
0396
0397   //    set the mass of the daughter particles
0398   switch(_GGsingInputpidtest){
0399   case starlightConstants::ZOVERZ03:
0400   case starlightConstants::F2:
0401     mdec = starlightConstants::pionChargedMass;
0402     break;
0403   case starlightConstants::AXION:       // AXION HACK
0404     mdec = 0;//axion decays to two photons, set mass of decay products to zero
0405     break;
0406   case starlightConstants::F2PRIME:
0407     //  decays 50% to K+/K-, 50% to K_0's
0408     ytest = _randy.Rndom();
0409     if(ytest >= 0.5){
0410       mdec = starlightConstants::kaonChargedMass;
0411     }
0412     else{
0413       mdec = 0.493677;
0414     }
0415     break;
0416   default :
0417     cout<<"No default mass selected for single photon-photon particle, expect errant results"<<endl;
0418   }
0419
0420   //Calculating the momentum's magnitude
0421     if(W < 2*mdec){
0422       cout<<" ERROR: W="<<W<<endl;
0423       iFbadevent = 1;
0424       return;
0425     }
0426     pmag = sqrt(W*W/4. - mdec*mdec);
0427 //  }
0428   //     pick an orientation, based on the spin
0429   //      phi has a flat distribution in 2*pi
0430   phi = _randy.Rndom()*2.*starlightConstants::pi;
0431
0432   //     find theta, the angle between one of the outgoing particles and
0433   //    the beamline, in the frame of the two photons
0434   //this will depend on spin, F2,F2' and z/z03 all have spin 2, all other photonphoton-single mesons are handled by jetset/pythia
0435   //Applies to spin2 mesons.
0436  L300td:
0437   theta = starlightConstants::pi*_randy.Rndom();
0438   xtest = _randy.Rndom();
0439   dndtheta = sin(theta)*sin(theta)*sin(theta)*sin(theta)*sin(theta);
0440   if(xtest > dndtheta)
0441     goto L300td;
0442
0443   //     compute unboosted momenta
0444   px1 = sin(theta)*cos(phi)*pmag;
0445   py1 = sin(theta)*sin(phi)*pmag;
0446   pz1 = cos(theta)*pmag;
0447   px2 = -px1;
0448   py2 = -py1;
0449   pz2 = -pz1;
0450   //        compute energies
0451   //Changed mass to W 11/9/2000 SRK
0452   Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0);
0453   E1 = sqrt(mdec*mdec+px1*px1+py1*py1+pz1*pz1);
0454   E2 = sqrt(mdec*mdec+px2*px2+py2*py2+pz2*pz2);
0455
0456   //     Lorentz transform into the lab frame
0457   // betax,betay,betaz are the boost of the complete system
0458   betax = -(px0/Ecm);
0459   betay = -(py0/Ecm);
0460   betaz = -(pz0/Ecm);
0461
0462   transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent);
0463   transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent);
0464
0465
0466   if(iFbadevent == 1)
0467     return;
0468
0469   //       change particle id from that of parent to that of daughters
0470
0471   switch(_GGsingInputpidtest){
0472     //These decay into a pi+ pi- pair
0473   case starlightConstants::ZOVERZ03:
0474   case starlightConstants::F2:
0475     ipid=starlightConstants::PION;
0476     break;
0477   case starlightConstants::AXION:// AXION HACK
0478     ipid=starlightConstants::PHOTON;  // AXION HACK
0479     break;      // AXION HACK
0480   case starlightConstants::F2PRIME:
0481     if( ytest >= 0.5 )
0482       {
0483     //Decays 50/50 into k+ k- or k_s k_l
0484     ipid=starlightConstants::KAONCHARGE;
0485       }
0486     else
0487       {
0488     ipid=starlightConstants::KAONNEUTRAL;
0489       }
0490     break;
0491   default:
0492     cout<<"Rethink the daughter particles"<<endl;
0493   }
0494 }
0495
0496
0497 //______________________________________________________________________________
0498 starlightConstants::event Gammagammasingle::produceEvent(int &/*ievent*/)
0499 {
0500   // Not in use anymore, default event struct returned
0501   return starlightConstants::event();
0502 }
0503
0504
0505 //______________________________________________________________________________
0506 double Gammagammasingle::getMass()
0507 {
0508   using namespace starlightConstants;
0509   double singlemass=0.;
0510   switch(_GGsingInputpidtest){
0511   case starlightConstants::ETA:
0512     singlemass = starlightConstants::etaMass;
0513     break;
0514   case starlightConstants::ETAPRIME:
0515     singlemass = starlightConstants::etaPrimeMass;
0516     break;
0517   case starlightConstants::ETAC:
0518     singlemass = starlightConstants::etaCMass;
0519     break;
0520   case starlightConstants::F0:
0521     singlemass = starlightConstants::f0Mass;
0522     break;
0523   case starlightConstants::F2:
0524     singlemass = starlightConstants::f2Mass;
0525     break;
0526   case starlightConstants::A2:
0527     singlemass = starlightConstants::a2Mass;
0528     break;
0529   case starlightConstants::F2PRIME:
0530     singlemass = starlightConstants::f2PrimeMass;
0531     break;
0532   case starlightConstants::ZOVERZ03:
0533     singlemass = starlightConstants::zoverz03Mass;
0534     break;
0535   case starlightConstants::AXION: // AXION HACK
0536     singlemass = _axionMass;      // AXION HACK
0537     break; // AXION HACK
0538   default:
0539     cout<<"Not a recognized single particle, Gammagammasingle::getmass(), mass = 0."<<endl;
0540   }
0541   return singlemass;
0542 }
0543
0544
0545 //______________________________________________________________________________
0546 double Gammagammasingle::getWidth()
0547 {
0548
0549   /* Partial widths(GAMMA(gammgamma)) taken from PDG 2014- Chinese Physics C 38, no 9, Sept. 2014.*/
0550   double singlewidth=0.;
0551   switch(_GGsingInputpidtest){
0552   case starlightConstants::ETA:
0553     singlewidth = starlightConstants::etaPartialggWidth;
0554     break;
0555   case starlightConstants::ETAPRIME:
0556     singlewidth = starlightConstants::etaPrimePartialggWidth;
0557     break;
0558   case starlightConstants::ETAC:
0559     singlewidth = starlightConstants::etaCPartialggWidth;
0560     break;
0561   case starlightConstants::F0:
0562     singlewidth = starlightConstants::f0PartialggWidth;
0563     break;
0564   case starlightConstants::F2:
0565     singlewidth = starlightConstants::f2PartialggWidth;
0566     break;
0567   case starlightConstants::A2:
0568     singlewidth = starlightConstants::a2PartialggWidth;
0569     break;
0570   case starlightConstants::F2PRIME:
0571     singlewidth = starlightConstants::f2PrimePartialggWidth;
0572     break;
0573   case starlightConstants::ZOVERZ03:
0574     singlewidth = starlightConstants::zoverz03PartialggWidth;
0575     break;
0576   case starlightConstants::AXION: // AXION HACK
0577     singlewidth = 1/(64*starlightConstants::pi)*_axionMass*_axionMass*_axionMass/(1000*1000);//Fix Lambda=1000 GeV,rescaling is trivial.    // AXION HACK
0578     break;
0579   default:
0580     cout<<"Not a recognized single particle, Gammagammasingle::getwidth(), width = 0."<<endl;
0581   }
0582   return singlewidth;
0583 }
0584
0585
0586 //______________________________________________________________________________
0587 double Gammagammasingle::getSpin()
0588 {
0589   double singlespin=0.5;
0590   switch(_GGsingInputpidtest){
0591   case starlightConstants::ETA:
0592     singlespin = starlightConstants::etaSpin;
0593     break;
0594   case starlightConstants::ETAPRIME:
0595     singlespin = starlightConstants::etaPrimeSpin;
0596     break;
0597   case starlightConstants::ETAC:
0598     singlespin = starlightConstants::etaCSpin;
0599     break;
0600   case starlightConstants::F0:
0601     singlespin = starlightConstants::f0Spin;
0602     break;
0603   case starlightConstants::F2:
0604     singlespin = starlightConstants::f2Spin;
0605     break;
0606   case starlightConstants::A2:
0607     singlespin = starlightConstants::a2Spin;
0608     break;
0609   case starlightConstants::F2PRIME:
0610     singlespin = starlightConstants::f2PrimeSpin;
0611     break;
0612   case starlightConstants::ZOVERZ03:
0613     singlespin = starlightConstants::zoverz03Spin;
0614     break;
0615   case starlightConstants::AXION:// AXION HACK
0616     singlespin = starlightConstants::axionSpin;// AXION HACK
0617     break;// AXION HACK
0618   default:
0619     cout<<"Not a recognized single particle, Gammagammasingle::getspin(), spin = 0."<<endl;
0620   }
0621   return singlespin;
0622 }
0623 eXEvent Gammagammasingle::e_produceEvent()
0624 {
0625   eXEvent event;
0626   cout<< " Gamma+Gamma -> single particle is not implemente in current eSTARlight build. REturning empty event"<<endl;
0627   return event;
0628 }
0629
0630