estarlight/src/incoherentPhotonNucleusLuminosity.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:: 45                          $: revision of last commit
0024 // $Author:: bgrube                   $: author of last commit
0025 // $Date:: 2011-02-27 20:52:35 +0100 #$: date of last commit
0026 //
0027 // Description:
0028 //
0029 //
0030 //
0031 ///////////////////////////////////////////////////////////////////////////
0032
0033
0034 #include <iostream>
0035 #include <fstream>
0036 #include <cmath>
0037
0038 #include "inputParameters.h"
0039 #include "beambeamsystem.h"
0040 #include "beam.h"
0041 #include "starlightconstants.h"
0042 #include "nucleus.h"
0043 #include "bessel.h"
0044 #include "incoherentPhotonNucleusLuminosity.h"
0045
0046
0047 using namespace std;
0048 using namespace starlightConstants;
0049
0050
0051 //______________________________________________________________________________
0052 incoherentPhotonNucleusLuminosity::incoherentPhotonNucleusLuminosity(const inputParameters& inputParametersInstance, beamBeamSystem& bbsystem)
0053   : photonNucleusCrossSection(inputParametersInstance, bbsystem)
0054   ,_baseFileName(inputParametersInstance.baseFileName())
0055   ,_beamLorentzGamma(inputParametersInstance.beamLorentzGamma())
0056   ,_maxW(inputParametersInstance.maxW())
0057   ,_minW(inputParametersInstance.minW())
0058   ,_nmbWBins(inputParametersInstance.nmbWBins())
0059   ,_maxRapidity(inputParametersInstance.maxRapidity())
0060   ,_nmbRapidityBins(inputParametersInstance.nmbRapidityBins())
0061   ,_productionMode(inputParametersInstance.productionMode())
0062   ,_beamBreakupMode(inputParametersInstance.beamBreakupMode())
0063   ,_interferenceEnabled(inputParametersInstance.interferenceEnabled())
0064   ,_interferenceStrength(inputParametersInstance.interferenceStrength())
0065   ,_maxPtInterference(inputParametersInstance.maxPtInterference())
0066   ,_nmbPtBinsInterference(inputParametersInstance.nmbPtBinsInterference())
0067   ,_protonEnergy(inputParametersInstance.protonEnergy())
0068   ,_parameterValueKey(inputParametersInstance.parameterValueKey())
0069 {
0070   cout <<"Creating Luminosity Tables for incoherent vector meson production."<<endl;
0071   incoherentPhotonNucleusDifferentialLuminosity();
0072   cout <<"Luminosity Tables created."<<endl;
0073 }
0074
0075
0076 //______________________________________________________________________________
0077 incoherentPhotonNucleusLuminosity::~incoherentPhotonNucleusLuminosity()
0078 { }
0079
0080
0081 //______________________________________________________________________________
0082 void incoherentPhotonNucleusLuminosity::incoherentPhotonNucleusDifferentialLuminosity()
0083 {
0084   double W,dW,dY;
0085   double Egamma,Y;
0086   double testint,dndWdY;
0087   double C;
0088   int beam;
0089
0090   std::string wyFileName;
0091   wyFileName = _baseFileName +".txt";
0092
0093   ofstream wylumfile;
0094   wylumfile.precision(15);
0095
0096   double  bwnorm,Eth;
0097
0098   dW = (_wMax - _wMin)/_nWbins;
0099   dY  = (_yMax-(-1.0)*_yMax)/_nYbins;
0100
0101   // Write the values of W used in the calculation to slight.txt.
0102   wylumfile.open(wyFileName.c_str());
0103   wylumfile << getbbs().electronBeam().Z() <<endl;
0104   wylumfile << getbbs().electronBeam().A() <<endl;
0105   wylumfile << getbbs().targetBeam().Z() <<endl;
0106   wylumfile << getbbs().targetBeam().A() <<endl;
0107   wylumfile << _beamLorentzGamma <<endl;
0108   wylumfile << _maxW <<endl;
0109   wylumfile << _minW <<endl;
0110   wylumfile << _nmbWBins <<endl;
0111   wylumfile << _maxRapidity <<endl;
0112   wylumfile << _nmbRapidityBins <<endl;
0113   wylumfile << _productionMode <<endl;
0114   wylumfile << _beamBreakupMode <<endl;
0115   wylumfile << _interferenceEnabled <<endl;
0116   wylumfile << _interferenceStrength <<endl;
0117   wylumfile << starlightConstants::deuteronSlopePar <<endl;
0118   wylumfile << _maxPtInterference <<endl;
0119   wylumfile << _nmbPtBinsInterference <<endl;
0120
0121   //     Normalize the Breit-Wigner Distribution and write values of W to slight.txt
0122   testint=0.0;
0123   //Grabbing default value for C in the breit-wigner calculation
0124   C=getDefaultC();
0125   for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
0126     W = _wMin + double(i)*dW + 0.5*dW;
0127     testint = testint + breitWigner(W,C)*dW;
0128     wylumfile << W << endl;
0129   }
0130   bwnorm = 1./testint;
0131
0132   //     Write the values of Y used in the calculation to slight.txt.
0133   for(unsigned int i = 0; i <= _nYbins - 1; ++i) {
0134     Y = -1.0*_yMax + double(i)*dY + 0.5*dY;
0135     wylumfile << Y << endl;
0136   }
0137
0138   int A_1 = getbbs().electronBeam().A();
0139   int A_2 = getbbs().targetBeam().A();
0140
0141   // Do this first for the case when the first beam is the photon emitter
0142   // Treat pA separately with defined beams
0143   // The variable beam (=1,2) defines which nucleus is the target
0144   for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
0145
0146     W = _wMin + double(i)*dW + 0.5*dW;
0147
0148     double Ep = _protonEnergy;
0149
0150     Eth=0.5*(((W+starlightConstants::protonMass)*(W+starlightConstants::protonMass)-starlightConstants::protonMass*starlightConstants::protonMass)/
0151        (Ep + sqrt(Ep*Ep-starlightConstants::protonMass*starlightConstants::protonMass)));
0152
0153     for(unsigned int j = 0; j <= _nYbins - 1; ++j) {
0154
0155       Y = -1.0*_yMax + double(j)*dY + 0.5*dY;
0156
0157       if( A_2 == 1 && A_1 != 1 ){
0158         // pA, first beam is the nucleus and photon emitter
0159         Egamma = 0.5*W*exp(Y);
0160         beam = 2;
0161       } else if( A_1 ==1 && A_2 != 1){
0162         // pA, second beam is the nucleus and photon emitter
0163         Egamma = 0.5*W*exp(-Y);
0164         beam = 1;
0165       } else {
0166         Egamma = 0.5*W*exp(Y);
0167         beam = 2;
0168       }
0169
0170       dndWdY = 0.;
0171
0172       if(Egamma > Eth){
0173     if(Egamma > maxPhotonEnergy())Egamma = maxPhotonEnergy();
0174         double Wgp = sqrt(2.*Egamma*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass*
0175                                  starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
0176
0177         double localsig = sigmagp(Wgp);
0178         if( A_1 == 1 && A_2 != 1 ){
0179           dndWdY = Egamma*photonFlux(Egamma,beam)*localsig*breitWigner(W,bwnorm);
0180         }else if (A_2 ==1 && A_1 !=1){
0181           dndWdY = Egamma*photonFlux(Egamma,beam)*localsig*breitWigner(W,bwnorm);
0182         }else{
0183           double csVN = sigma_N(Wgp);
0184           double csVA = sigma_A(csVN,beam);
0185           double csgA= (csVA/csVN)*sigmagp(Wgp);
0186           dndWdY = Egamma*photonFlux(Egamma,beam)*csgA*breitWigner(W,bwnorm);
0187         }
0188       }
0189
0190       wylumfile << dndWdY << endl;
0191
0192     }
0193   }
0194
0195   // Repeat the loop for the case when the second beam is the photon emitter.
0196   // Don't repeat for pA
0197   if( !( (A_2 == 1 && A_1 != 1) || (A_1 == 1 && A_2 != 1) ) ){
0198     for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
0199
0200       W = _wMin + double(i)*dW + 0.5*dW;
0201
0202       double Ep = _protonEnergy;
0203
0204       Eth=0.5*(((W+starlightConstants::protonMass)*(W+starlightConstants::protonMass)-starlightConstants::protonMass*starlightConstants::protonMass)/
0205        (Ep + sqrt(Ep*Ep-starlightConstants::protonMass*starlightConstants::protonMass)));
0206
0207       for(unsigned int j = 0; j <= _nYbins - 1; ++j) {
0208
0209         Y = -1.0*_yMax + double(j)*dY + 0.5*dY;
0210
0211         beam = 1;
0212         Egamma = 0.5*W*exp(-Y);
0213
0214         dndWdY = 0.;
0215
0216         if(Egamma > Eth){
0217       if(Egamma > maxPhotonEnergy())Egamma = maxPhotonEnergy();
0218           double Wgp = sqrt(2.*Egamma*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass*
0219                                  starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
0220
0221           double csVN = sigma_N(Wgp);
0222           double csVA = sigma_A(csVN,beam);
0223           double csgA= (csVA/csVN)*sigmagp(Wgp);
0224           dndWdY = Egamma*photonFlux(Egamma,beam)*csgA*breitWigner(W,bwnorm);
0225
0226         }
0227
0228         wylumfile << dndWdY << endl;
0229
0230       }
0231     }
0232   }
0233
0234   wylumfile << bwnorm << endl;
0235   wylumfile << _parameterValueKey << endl;
0236   wylumfile.close();
0237
0238 }
0239
0240