STARLIGHT/starlight/src/incoherentPhotonNucleusLuminosity.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:: 45                          $: revision of last commit
  24 // $Author:: bgrube                   $: author of last commit
  25 // $Date:: 2011-02-27 20:52:35 +0100 #$: date of last commit
  26 //
  27 // Description:
  28 //
  29 //
  30 //
  31 ///////////////////////////////////////////////////////////////////////////
  32
  33
  34 #include <iostream>
  35 #include <fstream>
  36 #include <cmath>
  37
  38 #include "inputParameters.h"
  39 #include "beambeamsystem.h"
  40 #include "beam.h"
  41 #include "starlightconstants.h"
  42 #include "nucleus.h"
  43 #include "bessel.h"
  44 #include "incoherentPhotonNucleusLuminosity.h"
  45
  46
  47 using namespace std;
  48
  49
  50 //______________________________________________________________________________
  51 incoherentPhotonNucleusLuminosity::incoherentPhotonNucleusLuminosity(beamBeamSystem& bbsystem)
  52   : photonNucleusCrossSection(bbsystem)
  53 {
  54   cout <<"Creating Luminosity Tables for incoherent vector meson production."<<endl;
  55   incoherentPhotonNucleusDifferentialLuminosity();
  56   cout <<"Luminosity Tables created."<<endl;
  57 }
  58
  59
  60 //______________________________________________________________________________
  61 incoherentPhotonNucleusLuminosity::~incoherentPhotonNucleusLuminosity()
  62 { }
  63
  64
  65 //______________________________________________________________________________
  66 void incoherentPhotonNucleusLuminosity::incoherentPhotonNucleusDifferentialLuminosity()
  67 {
  68         // double Av,Wgp,cs,cvma;
  69   double W,dW,dY;
  70   double Egamma,Y;
  71   // double t,tmin,tmax;
  72   double testint,dndWdY;
  73   // double ax,bx;
  74   double C;
  75
  76   ofstream wylumfile;
  77   wylumfile.precision(15);
  78
  79   double  bwnorm,Eth;
  80
  81   dW = (_wMax - _wMin)/_nWbins;
  82   dY  = (_yMax-(-1.0)*_yMax)/_nYbins;
  83
  84   // Write the values of W used in the calculation to slight.txt.
  85   wylumfile.open("slight.txt");
  86   wylumfile << inputParametersInstance.parameterValueKey() << endl;
  87   wylumfile << getbbs().beam1().Z() <<endl;
  88   wylumfile << getbbs().beam1().A() <<endl;
  89   wylumfile << getbbs().beam2().Z() <<endl;
  90   wylumfile << getbbs().beam2().A() <<endl;
  91   wylumfile << inputParametersInstance.beamLorentzGamma() <<endl;
  92   wylumfile << inputParametersInstance.maxW() <<endl;
  93   wylumfile << inputParametersInstance.minW() <<endl;
  94   wylumfile << inputParametersInstance.nmbWBins() <<endl;
  95   wylumfile << inputParametersInstance.maxRapidity() <<endl;
  96   wylumfile << inputParametersInstance.nmbRapidityBins() <<endl;
  97   wylumfile << inputParametersInstance.productionMode() <<endl;
  98   wylumfile << inputParametersInstance.beamBreakupMode() <<endl;
  99   wylumfile << inputParametersInstance.interferenceEnabled() <<endl;
 100   wylumfile << inputParametersInstance.interferenceStrength() <<endl;
 101   wylumfile << inputParametersInstance.coherentProduction() <<endl;
 102   wylumfile << inputParametersInstance.incoherentFactor() <<endl;
 103   wylumfile << inputParametersInstance.deuteronSlopePar() <<endl;
 104   wylumfile << inputParametersInstance.maxPtInterference() <<endl;
 105   wylumfile << inputParametersInstance.nmbPtBinsInterference() <<endl;
 106
 107   //     Normalize the Breit-Wigner Distribution and write values of W to slight.txt
 108   testint=0.0;
 109   //Grabbing default value for C in the breit-wigner calculation
 110   C=getDefaultC();
 111   for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
 112     W = _wMin + double(i)*dW + 0.5*dW;
 113     testint = testint + breitWigner(W,C)*dW;
 114     wylumfile << W << endl;
 115   }
 116   bwnorm = 1./testint;
 117
 118   //     Write the values of Y used in the calculation to slight.txt.
 119   for(unsigned int i = 0; i <= _nYbins - 1; ++i) {
 120     Y = -1.0*_yMax + double(i)*dY + 0.5*dY;
 121     wylumfile << Y << endl;
 122   }
 123
 124   //  Eth=0.5*(((_wMin+starlightConstants::protonMass)*(_wMin
 125   //                                                        +starlightConstants::protonMass)-starlightConstants::protonMass*starlightConstants::protonMass)/
 126   //       (Ep + sqrt(Ep*Ep-starlightConstants::protonMass*starlightConstants::protonMass)));
 127
 128   for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
 129
 130     W = _wMin + double(i)*dW + 0.5*dW;
 131
 132     double Ep = inputParametersInstance.protonEnergy();
 133
 134     Eth=0.5*(((W+starlightConstants::protonMass)*(W+starlightConstants::protonMass)-starlightConstants::protonMass*starlightConstants::protonMass)/
 135            (Ep + sqrt(Ep*Ep-starlightConstants::protonMass*starlightConstants::protonMass)));
 136
 137     for(unsigned int j = 0; j <= _nYbins - 1; ++j) {
 138
 139       Y = -1.0*_yMax + double(j)*dY + 0.5*dY;
 140
 141       int A_1 = getbbs().beam1().A();
 142       int A_2 = getbbs().beam2().A();
 143       if( A_2 == 1 && A_1 != 1 ){
 144         // pA, first beam is the nucleus
 145         Egamma = 0.5*W*exp(Y);
 146       } else if( A_1 ==1 && A_2 != 1){
 147         // pA, second beam is the nucleus
 148         Egamma = 0.5*W*exp(-Y);
 149       } else {
 150         Egamma = 0.5*W*exp(Y);
 151       }
 152
 153       dndWdY = 0.;
 154
 155       if(Egamma > Eth){
 156         if(Egamma > maxPhotonEnergy())Egamma = maxPhotonEnergy();
 157         double Wgp = sqrt(2.*Egamma*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass*
 158                                  starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
 159
 160         double localsig = sigmagp(Wgp);
 161         // int localz = 0;
 162         // double localbmin = 0;
 163         if( A_1 == 1 && A_2 != 1 ){
 164           // localbmin = getbbs().beam2().nuclearRadius() + 0.7;
 165           // localz = getbbs().beam2().Z();
 166           //   dndWdY = Egamma*localz*localz*nepoint(Egamma,localbmin)*localsig*breitWigner(W,bwnorm);
 167           dndWdY = Egamma*photonFlux(Egamma)*localsig*breitWigner(W,bwnorm);
 168         }else if (A_2 ==1 && A_1 !=1){
 169           // localbmin = getbbs().beam1().nuclearRadius() + 0.7;
 170           // localz = getbbs().beam1().Z();
 171           //   dndWdY = Egamma*localz*localz*nepoint(Egamma,localbmin)*localsig*breitWigner(W,bwnorm);
 172           dndWdY = Egamma*photonFlux(Egamma)*localsig*breitWigner(W,bwnorm);
 173         }else{
 174           double csVN = sigma_N(Wgp);
 175           double csVA = sigma_A(csVN);
 176           double csgA= (csVA/csVN)*sigmagp(Wgp);
 177           dndWdY = Egamma*photonFlux(Egamma)*csgA*breitWigner(W,bwnorm);
 178         }
 179       }
 180
 181       wylumfile << dndWdY << endl;
 182     }
 183   }
 184
 185   wylumfile << bwnorm << endl;
 186   wylumfile << inputParametersInstance.parameterValueKey() << endl;
 187   wylumfile.close();
 188
 189 //   wylumfile.open("slight.txt",ios::app);
 190   cout << "bwnorm: "<< bwnorm <<endl;
 191 //   wylumfile << bwnorm << endl;
 192 //   wylumfile << inputParametersInstance.parameterValueKey() << endl;
 193 //   wylumfile.close();
 194 }
 195
 196
 197 //______________________________________________________________________________
 198 double incoherentPhotonNucleusLuminosity::nofe(double Egamma, double bimp)
 199 {
 200   //Function for the calculation of the "photon density".
 201   //nofe=numberofgammas/(energy*area)
 202   //Assume beta=1.0 and gamma>>1, i.e. neglect the (1/gamma^2)*K0(x) term
 203
 204   double X=0.,nofex=0.,factor1=0.,factor2=0.,factor3=0.;
 205
 206   X = (bimp*Egamma)/(_beamLorentzGamma*starlightConstants::hbarc);
 207
 208   if( X <= 0.0)
 209     cout<<"In nofe, X= "<<X<<endl;
 210
 211   factor1 = (double(getbbs().beam1().Z()*getbbs().beam1().Z())
 212              *starlightConstants::alpha)/(starlightConstants::pi*starlightConstants::pi);
 213
 214   factor2 = 1./(Egamma*bimp*bimp);
 215   factor3 = X*X*(bessel::dbesk1(X))*(bessel::dbesk1(X));
 216   nofex    = factor1*factor2*factor3;
 217   return nofex;
 218 }