File indexing completed on 2024-09-27 07:03:40
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
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
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
0118 EQ2lumfile.open(EQ2FileName.c_str());
0119 EQ2lumfile << _nmbGammaQ2Bins <<endl;
0120
0121
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
0138 testint=0.0;
0139
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
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
0164 double dE_target = std::log(_targetMaxPhotonEnergy/_targetMinPhotonEnergy)/_nEBins;
0165
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
0173 Egamma = std::exp( std::log(_targetMinPhotonEnergy)+j*dE_target);
0174 g_E = 0;
0175
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
0182 g_E = Egamma*integrated_Q2_dep(Egamma, Q2min, Q2max);
0183 EQ2lumfile << g_E<<endl;
0184
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