estarlight/src/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
0033
0034 #include <iostream>
0035 #include <iomanip>
0036 #include <fstream>
0037 #include <cmath>
0038
0039 #include "starlightconstants.h"
0040 #include "narrowResonanceCrossSection.h"
0041
0042
0043 using namespace std;
0044 using namespace starlightConstants;
0045
0046
0047 //______________________________________________________________________________
0048 narrowResonanceCrossSection::narrowResonanceCrossSection(const inputParameters& inputParametersInstance, const beamBeamSystem&  bbsystem)
0049     :photonNucleusCrossSection(inputParametersInstance, bbsystem)
0050 {
0051     _narrowYmax = inputParametersInstance.maxRapidity();
0052     _narrowYmin = -1.0*_narrowYmax;
0053     _narrowNumY = inputParametersInstance.nmbRapidityBins();
0054     _Ep         = inputParametersInstance.protonEnergy();
0055 }
0056
0057
0058 //______________________________________________________________________________
0059 narrowResonanceCrossSection::~narrowResonanceCrossSection()
0060 { }
0061
0062
0063 //______________________________________________________________________________
0064 void
0065 narrowResonanceCrossSection::crossSectionCalculation(const double)  // _bwnormsave (unused)
0066 {
0067     // This subroutine calculates the vector meson cross section assuming
0068     // a narrow resonance.  For reference, see STAR Note 386.
0069
0070     double W,dY;
0071     double y1,y2,y12,ega1,ega2,ega12;
0072     double csgA1,csgA2,csgA12,int_r,dR;
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+nuclexon  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         y1  = _narrowYmin + double(J)*dY;
0099         y2  = _narrowYmin + double(J+1)*dY;
0100         y12 = 0.5*(y1+y2);
0101
0102                 if( A_2 == 1 && A_1 != 1 ){
0103                   // pA, first beam is the nucleus and is in this case the target
0104                   ega1  = 0.5*W*exp(-y1);
0105                   ega2  = 0.5*W*exp(-y2);
0106                   ega12 = 0.5*W*exp(-y12);
0107                   beam = 1;
0108                 } else if( A_1 ==1 && A_2 != 1){
0109                   // pA, second beam is the nucleus and is in this case the target
0110                   ega1  = 0.5*W*exp(y1);
0111                   ega2  = 0.5*W*exp(y2);
0112                   ega12 = 0.5*W*exp(y12);
0113                   beam = 2;
0114                 } else {
0115                   ega1  = 0.5*W*exp(y1);
0116                   ega2  = 0.5*W*exp(y2);
0117                   ega12 = 0.5*W*exp(y12);
0118                   beam = 2;
0119                 }
0120
0121         if(ega1 < Eth || ega2 < Eth)
0122             continue;
0123         if(ega2 > maxPhotonEnergy() || ega1 > maxPhotonEnergy() )
0124             continue;
0125
0126         csgA1=getcsgA(ega1,W,beam);
0127
0128         // Middle Point                      =====>>>
0129         csgA12=getcsgA(ega12,W,beam);
0130
0131         // Second Point                      =====>>>
0132         csgA2=getcsgA(ega2,W,beam);
0133
0134         // Sum the contribution for this W,Y.
0135         dR  = ega1*photonFlux(ega1,beam)*csgA1;
0136         dR  = dR + 4.*ega12*photonFlux(ega12,beam)*csgA12;
0137         dR  = dR + ega2*photonFlux(ega2,beam)*csgA2;
0138         dR  = dR*(dY/6.);
0139
0140         // cout<<" y: "<<y12<<" egamma: "<<ega12<<" flux: "<<photonFlux(ega12,beam)<<" sigma_gA: "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*dR/dY<<endl;
0141
0142         int_r = int_r+dR;
0143     }
0144
0145         // Repeat the loop for the case when the second beam is the photon emitter.
0146         // Don't repeat for pA
0147         if( !( (A_2 == 1 && A_1 != 1) || (A_1 == 1 && A_2 != 1) ) ){
0148       for(J=0;J<=(NY-1);J++){
0149
0150         y1  = _narrowYmin + double(J)*dY;
0151         y2  = _narrowYmin + double(J+1)*dY;
0152         y12 = 0.5*(y1+y2);
0153
0154                 beam = 1;
0155         ega1  = 0.5*W*exp(-y1);
0156         ega2  = 0.5*W*exp(-y2);
0157         ega12 = 0.5*W*exp(-y12);
0158
0159         if(ega2 < Eth || ega1 < Eth)
0160             continue;
0161         if(ega1 > maxPhotonEnergy() || ega2 > maxPhotonEnergy() )
0162             continue;
0163
0164         csgA1=getcsgA(ega1,W,beam);
0165
0166         // Middle Point                      =====>>>
0167         csgA12=getcsgA(ega12,W,beam);
0168
0169         // Second Point                      =====>>>
0170         csgA2=getcsgA(ega2,W,beam);
0171
0172         // Sum the contribution for this W,Y.
0173         dR  = ega1*photonFlux(ega1,beam)*csgA1;
0174         dR  = dR + 4.*ega12*photonFlux(ega12,beam)*csgA12;
0175         dR  = dR + ega2*photonFlux(ega2,beam)*csgA2;
0176         dR  = dR*(dY/6.);
0177
0178         // cout<<" y: "<<y12<<" egamma: "<<ega12<<" flux: "<<photonFlux(ega12,beam)<<" sigma_gA: "<<10000000.*csgA12<<" dsig/dy (microb): "<<10000.*dR/dY<<endl;
0179
0180         int_r = int_r+dR;
0181       }
0182         }
0183
0184     cout<<endl;
0185     if (0.01*int_r > 1.){
0186       cout<< " Total cross section: "<<0.01*int_r<<" barn."<<endl;
0187     } else if (10.*int_r > 1.){
0188       cout<< " Total cross section: " <<10.*int_r<<" mb."<<endl;
0189         } else if (10000.*int_r > 1.){
0190       cout<< " Total cross section: " <<10000.*int_r<<" microb."<<endl;
0191         } else if (10000000.*int_r > 1.){
0192       cout<< " Total cross section: " <<10000000.*int_r<<" nanob."<<endl;
0193         } else if (1.E10*int_r > 1.){
0194       cout<< " Total cross section: "<<1.E10*int_r<<" picob."<<endl;
0195         } else {
0196       cout<< " Total cross section: " <<1.E13*int_r<<" femtob."<<endl;
0197         }
0198     cout<<endl;
0199     setPhotonNucleusSigma(0.01*int_r);
0200 }