File indexing completed on 2024-09-27 07:03:42
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 "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)
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
0070
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;
0079
0080
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
0105
0106
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
0119
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
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
0144 csgA12=getcsgA(ega12,W,beam);
0145
0146
0147 csgA2=getcsgA(ega2,W,beam);
0148
0149
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
0156
0157
0158
0159
0160
0161 int_r = int_r+dR;
0162 }
0163 }
0164
0165
0166
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
0189 csgA12=getcsgA(ega12,W,beam);
0190
0191
0192 csgA2=getcsgA(ega2,W,beam);
0193
0194
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
0201
0202
0203
0204
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 }