estarlight/src/wideResonanceCrossSection.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 "starlightconstants.h"
0039 #include "wideResonanceCrossSection.h"
0040
0041
0042 using namespace std;
0043 using namespace starlightConstants;
0044
0045
0046 //______________________________________________________________________________
0047 wideResonanceCrossSection::wideResonanceCrossSection(const inputParameters& inputParametersInstance, const beamBeamSystem&  bbsystem)
0048     : photonNucleusCrossSection(inputParametersInstance, bbsystem)//hrm
0049 {
0050     _wideWmax = _wMax;
0051     _wideWmin = _wMin;
0052     _wideYmax = _yMax;
0053     _wideYmin = -1.0 * _wideYmax;
0054     _Ep       = inputParametersInstance.protonEnergy();
0055 }
0056
0057
0058 //______________________________________________________________________________
0059 wideResonanceCrossSection::~wideResonanceCrossSection()
0060 {
0061
0062 }
0063
0064
0065 //______________________________________________________________________________
0066 void
0067 wideResonanceCrossSection::crossSectionCalculation(const double bwnormsave)
0068 {
0069     //     This subroutine calculates the cross-section assuming a wide
0070     //     (Breit-Wigner) resonance.
0071
0072     double W,dW,dY;
0073     double y1,y2,y12,ega1,ega2,ega12;
0074     double csgA1,csgA2,csgA12,int_r,dR;
0075     double Eth;
0076     int    I,J,NW,NY,beam;
0077
0078     double bwnorm = bwnormsave; //used to transfer the bwnorm from the luminosity tables
0079
0080     //gamma+nucleon threshold.
0081     Eth=0.5*(((_wideWmin+protonMass)*(_wideWmin+protonMass)
0082               -protonMass*protonMass)/(_Ep+sqrt(_Ep*_Ep-protonMass*protonMass)));
0083
0084     NW   = 100;
0085     dW   = (_wideWmax-_wideWmin)/double(NW);
0086
0087     NY   =  1200;
0088     dY   = (_wideYmax-_wideYmin)/double(NY);
0089
0090     if (getBNORM()  ==  0.){
0091         cout<<" Using Breit-Wigner Resonance Profile."<<endl;
0092     }
0093     else{
0094         cout<<" Using Breit-Wigner plus direct pi+pi- profile."<<endl;
0095     }
0096
0097     cout<<" Integrating over W from "<<_wideWmin<<" to "<<_wideWmax<<endl;
0098
0099         int A_1 = getbbs().electronBeam().A();
0100         int A_2 = getbbs().targetBeam().A();
0101
0102     int_r=0.;
0103
0104         // Do this first for the case when the first beam is the photon emitter
0105         // Treat pA separately with defined beams
0106         // The variable beam (=1,2) defines which nucleus is the target
0107     for(I=0;I<=NW-1;I++){
0108
0109         W = _wideWmin + double(I)*dW + 0.5*dW;
0110
0111         for(J=0;J<=NY-1;J++){
0112
0113             y1  = _wideYmin + double(J)*dY;
0114             y2  = _wideYmin + double(J+1)*dY;
0115             y12 = 0.5*(y1+y2);
0116
0117                         if( A_2 == 1 && A_1 != 1 ){
0118                           // pA, first beam is the nucleus and is in this case the target
0119                           // Egamma = 0.5*W*exp(-Y);
0120                           ega1  = 0.5*W*exp(-y1);
0121                           ega2  = 0.5*W*exp(-y2);
0122                           ega12 = 0.5*W*exp(-y12);
0123                           beam = 1;
0124                         } else if( A_1 ==1 && A_2 != 1){
0125                           // pA, second beam is the nucleus and is in this case the target
0126                           ega1  = 0.5*W*exp(y1);
0127                           ega2  = 0.5*W*exp(y2);
0128                           ega12 = 0.5*W*exp(y12);
0129                           beam = 2;
0130                         } else {
0131                           ega1  = 0.5*W*exp(y1);
0132                           ega2  = 0.5*W*exp(y2);
0133                           ega12 = 0.5*W*exp(y12);
0134                           beam = 2;
0135                         }
0136
0137
0138             if(ega1 < Eth || ega2 < Eth) continue;
0139             if(ega2 > maxPhotonEnergy()) continue;
0140
0141             csgA1=getcsgA(ega1,W,beam);
0142
0143             //         >> Middle Point                      =====>>>
0144             csgA12=getcsgA(ega12,W,beam);
0145
0146             //         >> Second Point                      =====>>>
0147             csgA2=getcsgA(ega2,W,beam);
0148
0149             //>> Sum the contribution for this W,Y. The 2 accounts for the 2 beams
0150             dR  = ega1*photonFlux(ega1,beam)*csgA1;
0151             dR  = dR + 4.*ega12*photonFlux(ega12,beam)*csgA12;
0152             dR  = dR + ega2*photonFlux(ega2,beam)*csgA2;
0153             dR  = dR*(dY/6.)*breitWigner(W,bwnorm)*dW;
0154
0155             //For identical beams, we double.  Either may emit photon/scatter
0156             //For large differences in Z, we approx, that only electronBeam emits photon
0157             //and targetBeam scatters, eg d-Au where electronBeam=au and targetBeam=d
0158             //if(getbbs().electronBeam().A()==getbbs().targetBeam().A()){
0159             //  dR  = 2.*dR;
0160             //}
0161             int_r = int_r+dR;
0162         }
0163     }
0164
0165         // Repeat the loop for the case when the second beam is the photon emitter.
0166         // Don't repeat for pA
0167         if( !( (A_2 == 1 && A_1 != 1) || (A_1 == 1 && A_2 != 1) ) ){
0168       for(I=0;I<=NW-1;I++){
0169
0170         W = _wideWmin + double(I)*dW + 0.5*dW;
0171
0172         for(J=0;J<=NY-1;J++){
0173
0174             y1  = _wideYmin + double(J)*dY;
0175             y2  = _wideYmin + double(J+1)*dY;
0176             y12 = 0.5*(y1+y2);
0177
0178                         beam = 1;
0179             ega1  = 0.5*W*exp(-y1);
0180             ega2  = 0.5*W*exp(-y2);
0181             ega12 = 0.5*W*exp(-y12);
0182
0183             if(ega1< Eth || ega2 < Eth) continue;
0184             if(ega1 > maxPhotonEnergy()) continue;
0185
0186             csgA1=getcsgA(ega1,W,beam);
0187
0188             //         >> Middle Point                      =====>>>
0189             csgA12=getcsgA(ega12,W,beam);
0190
0191             //         >> Second Point                      =====>>>
0192             csgA2=getcsgA(ega2,W,beam);
0193
0194             //>> Sum the contribution for this W,Y. The 2 accounts for the 2 beams
0195             dR  = ega1*photonFlux(ega1,beam)*csgA1;
0196             dR  = dR + 4.*ega12*photonFlux(ega12,beam)*csgA12;
0197             dR  = dR + ega2*photonFlux(ega2,beam)*csgA2;
0198             dR  = dR*(dY/6.)*breitWigner(W,bwnorm)*dW;
0199
0200             //For identical beams, we double.  Either may emit photon/scatter
0201             //For large differences in Z, we approx, that only electronBeam emits photon
0202             //and targetBeam scatters, eg d-Au where electronBeam=au and targetBeam=d
0203             // if(getbbs().electronBeam().A()==getbbs().targetBeam().A()){
0204             //  dR  = 2.*dR;
0205             // }
0206             int_r = int_r+dR;
0207         }
0208       }
0209         }
0210
0211     cout<<endl;
0212     if (0.01*int_r > 1.){
0213       cout<< " Total cross section: "<<0.01*int_r<<" barn."<<endl;
0214     } else if (10.*int_r > 1.){
0215       cout<< " Total cross section: " <<10.*int_r<<" mb."<<endl;
0216         } else if (10000.*int_r > 1.){
0217       cout<< " Total cross section: " <<10000.*int_r<<" microb."<<endl;
0218         } else if (10000000.*int_r > 1.){
0219       cout<< " Total cross section: " <<10000000.*int_r<<" nanob."<<endl;
0220         } else if (1.E10*int_r > 1.){
0221       cout<< " Total cross section: "<<1.E10*int_r<<" picob."<<endl;
0222         } else {
0223       cout<< " Total cross section: " <<1.E13*int_r<<" femtob."<<endl;
0224         }
0225     cout<<endl;
0226     setPhotonNucleusSigma(0.01*int_r);
0227
0228 }