estarlight/src/incoherentVMCrossSection.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 "starlightconstants.h"
0039 #include "incoherentVMCrossSection.h"
0040
0041
0042 using namespace std;
0043 using namespace starlightConstants;
0044
0045
0046 //______________________________________________________________________________
0047 incoherentVMCrossSection::incoherentVMCrossSection(const inputParameters& inputParametersInstance, const beamBeamSystem&  bbsystem)
0048     :photonNucleusCrossSection(inputParametersInstance, bbsystem)
0049 {
0050     _narrowYmax = inputParametersInstance.maxRapidity();
0051     _narrowYmin = -1.0*_narrowYmax;
0052     _narrowNumY = inputParametersInstance.nmbRapidityBins();
0053     _Ep         = inputParametersInstance.protonEnergy();
0054 }
0055
0056
0057 //______________________________________________________________________________
0058 incoherentVMCrossSection::~incoherentVMCrossSection()
0059 { }
0060
0061
0062 //______________________________________________________________________________
0063 void
0064 incoherentVMCrossSection::crossSectionCalculation(const double)  // _bwnormsave (unused)
0065 {
0066     // This subroutine calculates the vector meson cross section assuming
0067     // a narrow resonance.  For reference, see STAR Note 386.
0068
0069     double W,dY;
0070     double y1,y2,y12,ega1,ega2,ega12;
0071     double csgA1,csgA2,csgA12,int_r,dR;
0072         double Wgp,csVN,csVA;
0073     double Eth;
0074     int    J,NY,beam;
0075
0076     NY   =  _narrowNumY;
0077     dY   = (_narrowYmax-_narrowYmin)/double(NY);
0078
0079     cout<<" Using Narrow Resonance ..."<<endl;
0080
0081     W = getChannelMass();
0082     Eth=0.5*(((W+protonMass)*(W+protonMass)-
0083               protonMass*protonMass)/(_Ep+sqrt(_Ep*_Ep-protonMass*protonMass)));
0084
0085     // cout<<" gamma+nucleon  Threshold: "<<Eth<<endl;
0086         printf(" gamma+nucleon threshold: %e GeV \n", Eth);
0087
0088         int A_1 = getbbs().electronBeam().A();
0089         int A_2 = getbbs().targetBeam().A();
0090
0091     int_r=0.;
0092
0093         // Do this first for the case when the first beam is the photon emitter
0094         // Treat pA separately with defined beams
0095         // The variable beam (=1,2) defines which nucleus is the target
0096     for(J=0;J<=(NY-1);J++){
0097
0098             // This is the fdefault
0099         y1  = _narrowYmin + double(J)*dY;
0100         y2  = _narrowYmin + double(J+1)*dY;
0101         y12 = 0.5*(y1+y2);
0102
0103                 if( A_2 == 1 && A_1 != 1 ){
0104                   // pA, first beam is the nucleus and photon emitter
0105           ega1  = 0.5*W*exp(y1);
0106           ega2  = 0.5*W*exp(y2);
0107           ega12 = 0.5*W*exp(y12);
0108                   beam = 2;
0109                 } else if( A_1 ==1 && A_2 != 1){
0110                   // pA, second beam is the nucleus and photon emitter
0111           ega1  = 0.5*W*exp(-y1);
0112           ega2  = 0.5*W*exp(-y2);
0113           ega12 = 0.5*W*exp(-y12);
0114                   beam = 1;
0115                 } else {
0116           ega1  = 0.5*W*exp(y1);
0117           ega2  = 0.5*W*exp(y2);
0118           ega12 = 0.5*W*exp(y12);
0119                   beam = 2;
0120                 }
0121
0122                 // This is for checking things in the lab frame
0123                 // y1lab  = _narrowYmin + double(J)*dY;
0124                 // y2lab  = _narrowYmin + double(J+1)*dY;
0125                 // y12lab = 0.5*(y1lab+y2lab);
0126
0127                 // p+Pb
0128                 // y1 = y1lab + 0.465;
0129                 // y2 = y2lab + 0.465;
0130                 // y12 = y12lab + 0.465;
0131                 // ega1  = 0.5*W*exp(y1);
0132                 // ega2  = 0.5*W*exp(y2);
0133                 // ega12 = 0.5*W*exp(y12);
0134
0135                 // Pb+p
0136                 // y1 = y1lab - 0.465;
0137                 // y2 = y2lab - 0.465;
0138                 // y12 = y12lab - 0.465;
0139                 // ega1  = 0.5*W*exp(-y1);
0140                 // ega2  = 0.5*W*exp(-y2);
0141                 // ega12 = 0.5*W*exp(-y12);
0142
0143
0144         if(ega1 < Eth)
0145             continue;
0146         if(ega2 > maxPhotonEnergy())
0147             continue;
0148
0149         // First point
0150                 Wgp = sqrt(2.*ega1*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
0151                       +starlightConstants::protonMass*starlightConstants::protonMass);
0152                 csVN = sigma_N(Wgp);
0153                 csVA = sigma_A(csVN,beam);
0154                 csgA1 = (csVA/csVN)*sigmagp(Wgp);
0155                 if( getbbs().electronBeam().A() == 1 || getbbs().targetBeam().A()==1 ){
0156                   csgA1 = sigmagp(Wgp);
0157                 }
0158
0159         // Middle point
0160                 Wgp = sqrt(2.*ega12*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
0161                       +starlightConstants::protonMass*starlightConstants::protonMass);
0162                 csVN = sigma_N(Wgp);
0163                 csVA = sigma_A(csVN,beam);
0164                 csgA12 = (csVA/csVN)*sigmagp(Wgp);
0165                 if( getbbs().electronBeam().A() == 1 || getbbs().targetBeam().A()==1 ){
0166                   csgA12 = sigmagp(Wgp);
0167                 }
0168
0169         // Last point
0170                 Wgp = sqrt(2.*ega2*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
0171                       +starlightConstants::protonMass*starlightConstants::protonMass);
0172                 csVN = sigma_N(Wgp);
0173                 csVA = sigma_A(csVN,beam);
0174                 csgA2 = (csVA/csVN)*sigmagp(Wgp);
0175                 if( getbbs().electronBeam().A() == 1 || getbbs().targetBeam().A()==1 ){
0176                   csgA2 = sigmagp(Wgp);
0177                 }
0178
0179                 dR = ega1*photonFlux(ega1,beam)*csgA1;
0180                 dR = dR + 4*ega12*photonFlux(ega12,beam)*csgA12;
0181                 dR = dR + ega2*photonFlux(ega2,beam)*csgA2;
0182                 dR = dR*(dY/6.);
0183
0184         // cout<<" y: "<<y12<<" egamma: "<<ega12<<" flux: "<<photonFlux(ega12)<<" sigma_gA: "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*dR/dY<<endl;
0185                 // cout<<" y: "<<y12lab<<" egamma: "<<ega12<<" flux: "<<ega12*photonFlux(ega12)<<" W: "<<Wgpm<<" Wflux: "<<Wgpm*photonFlux(ega12)<<" sigma_gA (nb): "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*ega12*photonFlux(ega12)*csgA12<<endl;
0186
0187         int_r = int_r+dR;
0188
0189     }
0190
0191         // Repeat the loop for the case when the second beam is the photon emitter.
0192         // Don't repeat for pA
0193         if( !( (A_2 == 1 && A_1 != 1) || (A_1 == 1 && A_2 != 1) ) ){
0194       for(J=0;J<=(NY-1);J++){
0195
0196             // This is the fdefault
0197         y1  = _narrowYmin + double(J)*dY;
0198         y2  = _narrowYmin + double(J+1)*dY;
0199         y12 = 0.5*(y1+y2);
0200
0201                 beam = 1;
0202         ega1  = 0.5*W*exp(-y1);
0203         ega2  = 0.5*W*exp(-y2);
0204         ega12 = 0.5*W*exp(-y12);
0205
0206                 // This is for checking things in the lab frame
0207                 // y1lab  = _narrowYmin + double(J)*dY;
0208                 // y2lab  = _narrowYmin + double(J+1)*dY;
0209                 // y12lab = 0.5*(y1lab+y2lab);
0210
0211                 // p+Pb
0212                 // y1 = y1lab + 0.465;
0213                 // y2 = y2lab + 0.465;
0214                 // y12 = y12lab + 0.465;
0215                 // ega1  = 0.5*W*exp(y1);
0216                 // ega2  = 0.5*W*exp(y2);
0217                 // ega12 = 0.5*W*exp(y12);
0218
0219                 // Pb+p
0220                 // y1 = y1lab - 0.465;
0221                 // y2 = y2lab - 0.465;
0222                 // y12 = y12lab - 0.465;
0223                 // ega1  = 0.5*W*exp(-y1);
0224                 // ega2  = 0.5*W*exp(-y2);
0225                 // ega12 = 0.5*W*exp(-y12);
0226
0227
0228         if(ega2 < Eth)
0229             continue;
0230         if(ega1 > maxPhotonEnergy())
0231             continue;
0232
0233         // First point
0234                 Wgp = sqrt(2.*ega1*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
0235                       +starlightConstants::protonMass*starlightConstants::protonMass);
0236                 csVN = sigma_N(Wgp);
0237                 csVA = sigma_A(csVN,beam);
0238                 csgA1 = (csVA/csVN)*sigmagp(Wgp);
0239                 if( getbbs().electronBeam().A() == 1 || getbbs().targetBeam().A()==1 ){
0240                   csgA1 = sigmagp(Wgp);
0241                 }
0242
0243         // Middle point
0244                 Wgp = sqrt(2.*ega12*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
0245                       +starlightConstants::protonMass*starlightConstants::protonMass);
0246                 csVN = sigma_N(Wgp);
0247                 csVA = sigma_A(csVN,beam);
0248                 csgA12 = (csVA/csVN)*sigmagp(Wgp);
0249                 if( getbbs().electronBeam().A() == 1 || getbbs().targetBeam().A()==1 ){
0250                   csgA12 = sigmagp(Wgp);
0251                 }
0252
0253         // Last point
0254                 Wgp = sqrt(2.*ega2*(_Ep+sqrt(_Ep*_Ep-starlightConstants::protonMass*starlightConstants::protonMass))
0255                       +starlightConstants::protonMass*starlightConstants::protonMass);
0256                 csVN = sigma_N(Wgp);
0257                 csVA = sigma_A(csVN,beam);
0258                 csgA2 = (csVA/csVN)*sigmagp(Wgp);
0259                 if( getbbs().electronBeam().A() == 1 || getbbs().targetBeam().A()==1 ){
0260                   csgA2 = sigmagp(Wgp);
0261                 }
0262
0263                 dR = ega1*photonFlux(ega1,beam)*csgA1;
0264                 dR = dR + 4*ega12*photonFlux(ega12,beam)*csgA12;
0265                 dR = dR + ega2*photonFlux(ega2,beam)*csgA2;
0266                 dR = dR*(dY/6.);
0267
0268         // cout<<" y: "<<y12<<" egamma: "<<ega12<<" flux: "<<photonFlux(ega12)<<" sigma_gA: "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*dR/dY<<endl;
0269                 // cout<<" y: "<<y12lab<<" egamma: "<<ega12<<" flux: "<<ega12*photonFlux(ega12)<<" W: "<<Wgpm<<" Wflux: "<<Wgpm*photonFlux(ega12)<<" sigma_gA (nb): "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*ega12*photonFlux(ega12)*csgA12<<endl;
0270
0271         int_r = int_r+dR;
0272
0273       }
0274         }
0275
0276     cout<<endl;
0277     if (0.01*int_r > 1.){
0278       cout<< " Total cross section: "<<0.01*int_r<<" barn."<<endl;
0279     } else if (10.*int_r > 1.){
0280       cout<< " Total cross section: " <<10.*int_r<<" mb."<<endl;
0281         } else if (10000.*int_r > 1.){
0282       cout<< " Total cross section: " <<10000.*int_r<<" microb."<<endl;
0283         } else if (10000000.*int_r > 1.){
0284       cout<< " Total cross section: " <<10000000.*int_r<<" nanob."<<endl;
0285         } else if (1.E10*int_r > 1.){
0286       cout<< " Total cross section: "<<1.E10*int_r<<" picob."<<endl;
0287         } else {
0288       cout<< " Total cross section: " <<1.E13*int_r<<" femtob."<<endl;
0289         }
0290     cout<<endl;
0291     setPhotonNucleusSigma(0.01*int_r);
0292 }