estarlight/src/e_narrowResonanceCrossSection.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 //#define _makeGammaPQ2_
0033
0034 #include <iostream>
0035 #include <iomanip>
0036 #include <fstream>
0037 #include <cmath>
0038
0039 #include "starlightconstants.h"
0040 #include "e_narrowResonanceCrossSection.h"
0041
0042 using namespace std;
0043 using namespace starlightConstants;
0044
0045 //______________________________________________________________________________
0046 e_narrowResonanceCrossSection::e_narrowResonanceCrossSection(const inputParameters& inputParametersInstance, const beamBeamSystem&  bbsystem)
0047     :photonNucleusCrossSection(inputParametersInstance, bbsystem)
0048 {
0049     _Ep         = inputParametersInstance.protonEnergy();
0050     //this is in target frame
0051     _electronEnergy = inputParametersInstance.electronEnergy();
0052     //_target_beamLorentz = inputParametersInstance.targetBeamLorentzGamma();
0053     _target_beamLorentz = inputParametersInstance.beamLorentzGamma();
0054     _boost = std::acosh(inputParametersInstance.electronBeamLorentzGamma())
0055       -std::acosh(inputParametersInstance.targetBeamLorentzGamma());
0056     _boost = _boost/2;
0057     // Photon energy limits in the two important frames
0058     _targetMaxPhotonEnergy=inputParametersInstance.targetMaxPhotonEnergy();
0059     _targetMinPhotonEnergy=inputParametersInstance.targetMinPhotonEnergy();
0060     _cmsMaxPhotonEnergy=inputParametersInstance.cmsMaxPhotonEnergy();
0061     _cmsMinPhotonEnergy=inputParametersInstance.cmsMinPhotonEnergy();
0062     //
0063     _VMnumEgamma = inputParametersInstance.nmbEnergyBins();
0064     _useFixedRange = inputParametersInstance.fixedQ2Range();
0065     _gammaMinQ2 = inputParametersInstance.minGammaQ2();
0066     _gammaMaxQ2 = inputParametersInstance.maxGammaQ2();
0067     _targetRadii = inputParametersInstance.targetRadius();
0068     _backwardsProduction = inputParametersInstance.backwardsProduction();
0069 }
0070
0071
0072 //______________________________________________________________________________
0073 e_narrowResonanceCrossSection::~e_narrowResonanceCrossSection()
0074 { }
0075
0076
0077 //______________________________________________________________________________
0078 void
0079 e_narrowResonanceCrossSection::crossSectionCalculation(const double)  // _bwnormsave (unused)
0080 {
0081         // This subroutine calculates the vector meson cross section assuming
0082         // a narrow resonance.  For reference, see STAR Note 386.
0083
0084         double dEgamma, minEgamma;
0085     double ega[3] = {0};
0086     double int_r,dR;
0087     double int_r2, dR2;
0088     int    iEgamma, nEgamma,beam;
0089
0090     //Integration is done with exponential steps, in target frame
0091     //nEgamma = _VMnumEgamma;
0092     nEgamma = 1000;
0093     dEgamma = std::log(_targetMaxPhotonEnergy/_targetMinPhotonEnergy)/nEgamma;
0094     minEgamma = std::log(_targetMinPhotonEnergy);
0095
0096     cout<<" Using Narrow Resonance ..."<<endl;
0097
0098     //
0099         printf(" gamma+nucleon threshold (CMS): %e GeV \n", _cmsMinPhotonEnergy);
0100
0101         int A_1 = getbbs().electronBeam().A();
0102         int A_2 = getbbs().targetBeam().A();
0103
0104     if( A_2 == 0 && A_1 >= 1 ){
0105       // pA, first beam is the nucleus and is in this case the target
0106       beam = 1;
0107     } else if( A_1 ==0 && A_2 >= 1){
0108       // photon energy in Target frame
0109       beam = 2;
0110     } else {
0111       // photon energy in Target frame
0112       beam = 2;
0113     }
0114
0115     int_r=0.;
0116     int_r2= 0;
0117         // Do this first for the case when the first beam is the photon emitter
0118         // The variable beam (=1,2) defines which nucleus is the target
0119     // Target beam ==2 so repidity is negative. Can generalize later
0120     // These might be useful in a future iteration
0121     int nQ2 = 1000;
0122         printf(" gamma+nucleon threshold (Target): %e GeV \n", _targetMinPhotonEnergy);
0123     for(iEgamma = 0 ; iEgamma < nEgamma; ++iEgamma){    // Integral over photon energy
0124       // Target frame energies
0125       ega[0] = exp(minEgamma + iEgamma*dEgamma );
0126       ega[1] = exp(minEgamma + (iEgamma+1)*dEgamma );
0127       ega[2] = 0.5*(ega[0]+ega[1]);
0128
0129       // Integral over Q2
0130       double dndE[3] = {0}; // Effective photon flux
0131       double full_int[3] = {0}; // Full e+X --> e+X+V.M. cross section
0132       //
0133       for( int iEgaInt = 0 ; iEgaInt < 3; ++iEgaInt){    // Loop over the energies for the three points to integrate over Q2
0134
0135         /*
0136         if ( _backwardsProduction)
0137           {
0138         if (_ip->productionMode() != 2 || _ip->prodParticleId() != 223 || _ip->beam1A() == 1 || _ip->beam2A() != 1)
0139           {
0140             std::cout << "Backward production is implemented for ";
0141             std::cout << "PROD_MODE = 2, PROD_PID = 223, ";
0142             std::cout << "BEAM_1_A != 1, and BEAM_2_A == 1. ";
0143             std::cout << "Terminating program." << std::endl;
0144             throw 0;
0145           }
0146           }
0147         */
0148         //These are the physical limits
0149         double Q2_min = std::pow(starlightConstants::mel*ega[iEgaInt],2.0)/(_electronEnergy*(_electronEnergy-ega[iEgaInt]));
0150         double Q2_max = 4.*_electronEnergy*(_electronEnergy-ega[iEgaInt]);
0151         if(_useFixedRange == true){
0152           if( Q2_min < _gammaMinQ2 )
0153         Q2_min = _gammaMinQ2;
0154           if( Q2_max > _gammaMaxQ2 )
0155         Q2_max = _gammaMaxQ2;
0156         }
0157         double lnQ2ratio = std::log(Q2_max/Q2_min)/nQ2;
0158         double lnQ2_min = std::log(Q2_min);
0159         //
0160
0161         for( int iQ2 = 0 ; iQ2 < nQ2; ++iQ2){     // Integral over photon virtuality
0162           //
0163           double q2_1 = exp( lnQ2_min + iQ2*lnQ2ratio);
0164           double q2_2 = exp( lnQ2_min + (iQ2+1)*lnQ2ratio);
0165           double q2_12 = (q2_2+q2_1)/2.;
0166           // Integrating the effective photon flux
0167           dndE[iEgaInt] +=(q2_2-q2_1)*( getcsgA_Q2_dep(q2_1)*photonFlux(ega[iEgaInt],q2_1)
0168                         +getcsgA_Q2_dep(q2_2)*photonFlux(ega[iEgaInt],q2_2)
0169                         +4.*getcsgA_Q2_dep(q2_12)*photonFlux(ega[iEgaInt],q2_12) );
0170           // Integrating cross section
0171           full_int[iEgaInt] += (q2_2-q2_1)*( g(ega[iEgaInt],q2_1)*getcsgA(ega[iEgaInt],q2_1,beam)
0172                          + g(ega[iEgaInt],q2_2)*getcsgA(ega[iEgaInt],q2_2,beam)
0173                          + 4.*g(ega[iEgaInt],q2_12)*getcsgA(ega[iEgaInt],q2_12,beam) );
0174         }
0175
0176         // Finish the Q2 integration for the three end-points (Siumpsons rule)
0177         dndE[iEgaInt] = dndE[iEgaInt]/6.;
0178         full_int[iEgaInt] = full_int[iEgaInt]/6.;
0179       }
0180       // Finishing cross-section integral
0181       dR = full_int[0];
0182       dR += full_int[1];
0183       dR += 4.*full_int[2];
0184       dR = dR*(ega[1]-ega[0])/6.;
0185
0186       // Finishing integral over the effective photon flux
0187       dR2 = dndE[0];
0188       dR2 += dndE[1];
0189       dR2 += 4.*dndE[2];
0190       dR2 = dR2*(ega[1]-ega[0])/6.;
0191       //
0192       int_r = int_r+dR;
0193       int_r2 = int_r2 + dR2;
0194     }
0195     cout<<endl;
0196     if(_useFixedRange == true){
0197       cout<<" Using fixed Q2 range "<<_gammaMinQ2 << " < Q2 < "<<_gammaMaxQ2<<endl;
0198     }
0199
0200     if(_backwardsProduction){cout<<"********************** WARNING ***********************"<<endl;
0201                              cout<<"The total cross section for backward production"<<endl;
0202                              cout<<"should only be trusted for the omega. The exact omega"<<endl;
0203                              cout<<"behavior is implemented for the rho. Other mesons use"<<endl;
0204                              cout<<"the forward cross section dependence"<<endl;}
0205     printCrossSection(" Total cross section: ",int_r);
0206     if(_backwardsProduction)cout<<"******************************************************"<<endl;
0207
0208     //printCrossSection(" gamma+X --> VM+X ", int_r/int_r2);      // commented out for the mean time, not necesary in current implementation
0209     //
0210     //cout<<endl;
0211     setPhotonNucleusSigma(0.01*int_r);
0212     #ifdef _makeGammaPQ2_
0213     makeGammaPQ2dependence();
0214     #endif
0215 }
0216
0217
0218 //______________________________________________________________________________
0219 void
0220 e_narrowResonanceCrossSection::makeGammaPQ2dependence()
0221 {
0222     // This subroutine calculates the Q2 dependence of
0223         // gamma+X -> VM + X cross section for a narrow resonance
0224
0225   /*int const nQ2bins = 19;
0226     double const q2Edge[nQ2bins+1] = { 0.1,1.,2.,3., 4.,5.,
0227                        6.,7.,8.,9.,10.,
0228                        11.,12.,13.,14.,15.,
0229                        20.,30.,40.,50.};
0230   */
0231         int const nQ2bins = 38;
0232     double const q2Edge[nQ2bins+1] = { 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09,
0233                        0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9,
0234                        1, 2, 3, 4, 5, 6, 7, 8, 9,
0235                        10, 20, 30, 40, 50, 60, 70, 80, 90,
0236                        100, 200 };
0237     //
0238     double full_x_section[nQ2bins] = {0};
0239     double effective_flux[nQ2bins] = {0};
0240     //
0241         ofstream gamma_flux,total_xsec;
0242     //
0243     gamma_flux.open("estarlight_gamma_flux.csv");
0244     total_xsec.open("estarlight_total_xsec.csv");
0245     //
0246         double dEgamma, minEgamma;
0247     double ega[3] = {0};
0248     double dR;
0249     double dR2;
0250     int    iEgamma, nEgamma,beam;
0251
0252     //Integration is done with exponential steps, in target frame
0253     //nEgamma = _VMnumEgamma;
0254     nEgamma = 1000;
0255     dEgamma = std::log(_targetMaxPhotonEnergy/_targetMinPhotonEnergy)/nEgamma;
0256     minEgamma = std::log(_targetMinPhotonEnergy);
0257
0258     //cout<<" Using Narrow Resonance ..."<<endl;
0259
0260         //printf(" gamma+nucleon threshold (CMS): %e GeV \n", _cmsMinPhotonEnergy);
0261
0262         int A_1 = getbbs().electronBeam().A();
0263         int A_2 = getbbs().targetBeam().A();
0264
0265     if( A_2 == 0 && A_1 >= 1 ){
0266       // pA, first beam is the nucleus and is in this case the target
0267       beam = 1;
0268     } else if( A_1 ==0 && A_2 >= 1){
0269       // photon energy in Target frame
0270       beam = 2;
0271     } else {
0272       // photon energy in Target frame
0273       beam = 2;
0274     }
0275         // Do this first for the case when the first beam is the photon emitter
0276         // The variable beam (=1,2) defines which nucleus is the target
0277     // Target beam ==2 so repidity is negative. Can generalize later
0278     int nQ2 = 500;
0279         //printf(" gamma+nucleon threshold (Target): %e GeV \n", _targetMinPhotonEnergy);
0280     for( int iQ2Bin  = 0 ; iQ2Bin < nQ2bins; ++iQ2Bin){
0281       for(iEgamma = 0 ; iEgamma < nEgamma; ++iEgamma){    // Integral over photon energy
0282         // Target frame energies
0283         ega[0] = exp(minEgamma + iEgamma*dEgamma );
0284         ega[1] = exp(minEgamma + (iEgamma+1)*dEgamma );
0285         ega[2] = 0.5*(ega[0]+ega[1]);
0286         //
0287         //
0288         // Integral over Q2
0289         double dndE[3] = {0}; // Effective photon flux
0290         double full_int[3] = {0}; // Full e+X --> e+X+V.M. cross section
0291         //
0292         for( int iEgaInt = 0 ; iEgaInt < 3; ++iEgaInt){    // Loop over the energies for the three points to integrate over Q2
0293           //
0294           double Q2_min = q2Edge[iQ2Bin];
0295           double Q2_max = q2Edge[iQ2Bin+1];
0296           double lnQ2ratio = std::log(Q2_max/Q2_min)/nQ2;
0297           double lnQ2_min = std::log(Q2_min);
0298           for( int iQ2 = 0 ; iQ2 < nQ2; ++iQ2){     // Integral over photon virtuality
0299         //
0300         double q2_1 = exp( lnQ2_min + iQ2*lnQ2ratio);
0301         double q2_2 = exp( lnQ2_min + (iQ2+1)*lnQ2ratio);
0302         double q2_12 = (q2_2+q2_1)/2.;
0303         // Integrating the photon flux
0304         dndE[iEgaInt] +=(q2_2-q2_1)*( photonFlux(ega[iEgaInt],q2_1)
0305                           +photonFlux(ega[iEgaInt],q2_2)
0306                           +4.*photonFlux(ega[iEgaInt],q2_12) );
0307
0308         full_int[iEgaInt] += (q2_2-q2_1)*( g(ega[iEgaInt],q2_1)*getcsgA(ega[iEgaInt],q2_1,beam)
0309                            + g(ega[iEgaInt],q2_2)*getcsgA(ega[iEgaInt],q2_2,beam)
0310                            + 4.*g(ega[iEgaInt],q2_12)*getcsgA(ega[iEgaInt],q2_12,beam) );
0311         }
0312           // Finish the Q2 integration for the three end-points (Siumpsons rule)
0313           dndE[iEgaInt] = dndE[iEgaInt]/6.;
0314           full_int[iEgaInt] = full_int[iEgaInt]/6.;
0315         }
0316         // Finishing cross-section integral
0317         dR = full_int[0];
0318         dR += full_int[1];
0319         dR += 4.*full_int[2];
0320         dR = dR*(ega[1]-ega[0])/6.;
0321
0322         // Finishing integral over the effective photon flux
0323         dR2 = dndE[0];
0324         dR2 += dndE[1];
0325         dR2 += 4.*dndE[2];
0326         dR2 = dR2*(ega[1]-ega[0])/6.;
0327         //
0328         full_x_section[iQ2Bin] += dR;
0329         effective_flux[iQ2Bin] += dR2;
0330       }
0331       //cout<<gamma_x_section[iQ2Bin]/effective_flux[iQ2Bin]*1E7<<endl;
0332       gamma_flux<<q2Edge[iQ2Bin]<<","<<q2Edge[iQ2Bin+1]<<","<<effective_flux[iQ2Bin]<<endl;
0333       total_xsec<<q2Edge[iQ2Bin]<<","<<q2Edge[iQ2Bin+1]<<","<<full_x_section[iQ2Bin]<<endl; //No need to remove bin width
0334     }
0335     //
0336     //
0337     gamma_flux.close();
0338     total_xsec.close();
0339 }
0340
0341 void e_narrowResonanceCrossSection::printCrossSection(const string name, const double x_section)
0342 {
0343   if (0.01*x_section > 1.){
0344     cout<< name.c_str() <<0.01*x_section<<" barn."<<endl;
0345   } else if (10.*x_section > 1.){
0346     cout<< name.c_str() <<10.*x_section<<" mb."<<endl;
0347   } else if (10000.*x_section > 1.){
0348     cout<< name.c_str() <<10000.*x_section<<" microb."<<endl;
0349   } else if (10000000.*x_section > 1.){
0350     cout<< name.c_str() <<10000000.*x_section<<" nanob."<<endl;
0351   } else if (1.E10*x_section > 1.){
0352     cout<< name.c_str() <<1.E10*x_section<<" picob."<<endl;
0353   } else {
0354     cout<< name.c_str() <<1.E13*x_section<<" femtob."<<endl;
0355   }
0356   //cout<<endl;
0357 }