estarlight/src/gammaeluminosity.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:: 264                         $: revision of last commit
0024 // $Author:: jseger                   $: author of last commit
0025 // $Date:: 2016-06-06 21:05:12 +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 "gammaeluminosity.h"
0045
0046
0047 using namespace std;
0048 using namespace starlightConstants;
0049
0050
0051 //______________________________________________________________________________
0052 photonElectronLuminosity::photonElectronLuminosity(const inputParameters& inputParametersInstance, beamBeamSystem& bbsystem)
0053   : photonNucleusCrossSection(inputParametersInstance, bbsystem)
0054   ,_protonEnergy(inputParametersInstance.protonEnergy())
0055   ,_electronEnergy(inputParametersInstance.electronEnergy())
0056   ,_beamLorentzGamma(inputParametersInstance.beamLorentzGamma())
0057   ,_baseFileName(inputParametersInstance.baseFileName())
0058   ,_maxW(inputParametersInstance.maxW())
0059   ,_minW(inputParametersInstance.minW())
0060   ,_nmbWBins(inputParametersInstance.nmbWBins())
0061   ,_maxRapidity(inputParametersInstance.maxRapidity())
0062   ,_nmbRapidityBins(inputParametersInstance.nmbRapidityBins())
0063   ,_nEBins(inputParametersInstance.nmbEnergyBins())
0064   ,_minGammaQ2(inputParametersInstance.minGammaQ2())
0065   ,_maxGammaQ2(inputParametersInstance.maxGammaQ2())
0066   ,_nmbGammaQ2Bins(inputParametersInstance.nmbGammaQ2Bins())
0067   ,_cmsMaxPhotonEnergy(inputParametersInstance.cmsMaxPhotonEnergy())
0068   ,_cmsMinPhotonEnergy(inputParametersInstance.cmsMinPhotonEnergy())
0069   ,_targetMaxPhotonEnergy(inputParametersInstance.targetMaxPhotonEnergy())
0070   ,_targetMinPhotonEnergy(inputParametersInstance.targetMinPhotonEnergy())
0071   ,_productionMode(inputParametersInstance.productionMode())
0072   ,_beamBreakupMode(inputParametersInstance.beamBreakupMode())
0073 {
0074   cout <<"Creating Luminosity Tables."<<endl;
0075   photonNucleusDifferentialLuminosity();
0076   cout <<"Luminosity Tables created."<<endl;
0077 }
0078
0079
0080 //______________________________________________________________________________
0081 photonElectronLuminosity::~photonElectronLuminosity()
0082 { }
0083
0084
0085 //______________________________________________________________________________
0086 void photonElectronLuminosity::photonNucleusDifferentialLuminosity()
0087 {
0088   double W,dW, dY;
0089   double Egamma,Y;
0090   //
0091   double testint;
0092   //
0093   double g_E;
0094   double csgA;
0095   double C;
0096   int beam;
0097   //
0098   //int nQ2steps =100;
0099
0100   std::string wyFileName;
0101   wyFileName = _baseFileName +".txt";
0102
0103   ofstream wylumfile;
0104   wylumfile.precision(15);
0105
0106   std::string EQ2FileName;
0107   EQ2FileName = "e_"+_baseFileName +".txt";
0108
0109   ofstream EQ2lumfile;
0110   EQ2lumfile.precision(15);
0111
0112   double  bwnorm;
0113
0114   dW = (_wMax-_wMin)/_nWbins;
0115   dY  = (_yMax-(-1.0)*_yMax)/_nYbins;
0116
0117   // Look-up table storing double and single photon flux differentials
0118   EQ2lumfile.open(EQ2FileName.c_str());
0119   EQ2lumfile << _nmbGammaQ2Bins <<endl;
0120   // Look-up table storing photo-nuclear cross section
0121   // Write the values of W used in the calculation to slight.txt.
0122   wylumfile.open(wyFileName.c_str());
0123   wylumfile << getbbs().electronBeam().Z() <<endl;
0124   wylumfile << getbbs().electronBeam().A() <<endl;
0125   wylumfile << getbbs().targetBeam().Z() <<endl;
0126   wylumfile << getbbs().targetBeam().A() <<endl;
0127   wylumfile << _beamLorentzGamma <<endl;
0128   wylumfile << _maxW <<endl;
0129   wylumfile << _minW <<endl;
0130   wylumfile << _nmbWBins <<endl;
0131   wylumfile << _maxRapidity <<endl;
0132   wylumfile << _nmbRapidityBins <<endl;
0133   wylumfile << _productionMode <<endl;
0134   wylumfile << _beamBreakupMode <<endl;
0135   wylumfile << starlightConstants::deuteronSlopePar <<endl;
0136
0137   // Normalize the Breit-Wigner Distribution and write values of W to slight.txt
0138   testint=0.0;
0139   // Grabbing default value for C in the breit-wigner calculation
0140   C=getDefaultC();
0141   for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
0142     W = _wMin + double(i)*dW + 0.5*dW;
0143     testint = testint + breitWigner(W,C)*dW;
0144     wylumfile << W << endl;
0145   }
0146   bwnorm = 1./testint;
0147
0148   // Write the values of Y used in the calculation to slight.txt.
0149   for(unsigned int i = 0; i <= _nYbins - 1; ++i) {
0150     Y = -1.0*_yMax + double(i)*dY + 0.5*dY;
0151     wylumfile << Y << endl;
0152   }
0153
0154   int A_1 = getbbs().electronBeam().A();
0155   int A_2 = getbbs().targetBeam().A();
0156   if( A_2 == 0 && A_1 != 0 ){
0157     beam = 1;
0158   } else if( A_1 ==0 && A_2 != 0){
0159     beam = 2;
0160   } else{
0161     beam = 2;
0162   }
0163   // Exponential steps for both tables: Target frame for photon generation and CMS frame for final state generation
0164   double dE_target = std::log(_targetMaxPhotonEnergy/_targetMinPhotonEnergy)/_nEBins;
0165   // - - - - Calculations in Target frame
0166   for(unsigned int i = 0; i < _nWbins; ++i) {
0167
0168     W = _wMin + double(i)*dW + 0.5*dW;
0169     wylumfile << breitWigner(W,bwnorm) <<endl;
0170   }
0171   for(int j = 0; j < _nEBins ; ++j) {
0172     // Used for calculation of g(Egamma) which must be done in the ion (target) rest frame
0173     Egamma = std::exp( std::log(_targetMinPhotonEnergy)+j*dE_target);
0174     g_E = 0;
0175     // Photon energy limits are determnined in target frame - multiply by Egamma to speed up event generation
0176     std::pair< double, double >* this_energy = Q2arraylimits(Egamma);
0177     double Q2min = this_energy->first;
0178     double Q2max = this_energy->second;
0179     if( Q2min > Q2max)
0180       continue;
0181     //Accounts for the phase space factor
0182     g_E = Egamma*integrated_Q2_dep(Egamma, Q2min, Q2max);
0183     EQ2lumfile << g_E<<endl;
0184     // Write out Q2 range used for generation
0185     EQ2lumfile << Q2min << endl;
0186     EQ2lumfile << Q2max << endl;
0187     for( unsigned int iQ2 =0 ;iQ2<_nmbGammaQ2Bins; ++iQ2){
0188       double Q2 = std::exp( std::log(Q2min)+iQ2*std::log(Q2max/Q2min)/_nmbGammaQ2Bins );
0189       EQ2lumfile<< g(Egamma,Q2) <<endl;
0190       csgA=getcsgA(Egamma,Q2,beam);
0191       wylumfile << csgA << endl;
0192     }
0193   }
0194
0195   EQ2lumfile.close();
0196
0197   wylumfile << bwnorm << endl;
0198   wylumfile.close();
0199
0200 }
0201 string photonElectronLuminosity::gammaTableParse(int ii, int jj)
0202 {
0203   ostringstream tag1, tag2;
0204   tag1<<ii;
0205   tag2<<jj;
0206   string to_ret = tag1.str()+","+tag2.str();
0207   return to_ret;
0208 }
0209