EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
e_wideResonanceCrossSection.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file e_wideResonanceCrossSection.cpp
1 
2 //
3 // Copyright 2010
4 //
5 // This file is part of starlight.
6 //
7 // starlight is free software: you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation, either version 3 of the License, or
10 // (at your option) any later version.
11 //
12 // starlight is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // You should have received a copy of the GNU General Public License
18 // along with starlight. If not, see <http://www.gnu.org/licenses/>.
19 //
21 //
22 // File and Version Information:
23 // $Rev:: 264 $: revision of last commit
24 // $Author:: mlomnitz $: author of last commit
25 // $Date:: 2017-03-14 21:05:12 +0100 #$: date of last commit
26 //
27 // Description:
28 //
29 //
30 //
32 //#define _makeGammaPQ2_
33 
34 #include <iostream>
35 #include <fstream>
36 #include <cmath>
37 
38 #include "starlightconstants.h"
39 #include "e_wideResonanceCrossSection.h"
40 
41 
42 using namespace std;
43 using namespace starlightConstants;
44 
45 
46 //______________________________________________________________________________
47 e_wideResonanceCrossSection::e_wideResonanceCrossSection(const inputParameters& inputParametersInstance, const beamBeamSystem& bbsystem)
48  : photonNucleusCrossSection(inputParametersInstance, bbsystem)//hrm
49 {
50  _wideWmax = _wMax;
51  _wideWmin = _wMin;
52  _targetMaxPhotonEnergy=inputParametersInstance.targetMaxPhotonEnergy();
53  _targetMinPhotonEnergy=inputParametersInstance.targetMinPhotonEnergy();
54  _Ep = inputParametersInstance.protonEnergy();
55  _electronEnergy = inputParametersInstance.electronEnergy();
56  _target_beamLorentz = inputParametersInstance.beamLorentzGamma();
57  _VMnumEgamma = inputParametersInstance.nmbEnergyBins();
58  _useFixedRange = inputParametersInstance.fixedQ2Range();
59  _gammaMinQ2 = inputParametersInstance.minGammaQ2();
60  _gammaMaxQ2 = inputParametersInstance.maxGammaQ2();
61  _targetRadii = inputParametersInstance.targetRadius();
62 }
63 
64 
65 //______________________________________________________________________________
66 e_wideResonanceCrossSection::~e_wideResonanceCrossSection()
67 {
68 
69 }
70 
71 
72 //______________________________________________________________________________
73 void
74 e_wideResonanceCrossSection::crossSectionCalculation(const double bwnormsave)
75 {
76  // This subroutine calculates the cross-section assuming a wide
77  // (Breit-Wigner) resonance.
78 
79  double W,dW, dEgamma, minEgamma;
80  double ega[3] = {0};
81  double int_r,dR;
82  double int_r2,dR2;
83  int iW,nW,iEgamma,nEgamma,beam;
84 
85  double bwnorm = bwnormsave; //used to transfer the bwnorm from the luminosity tables
86 
87  // For W integration
88  nW = 100;
89  dW = (_wideWmax-_wideWmin)/double(nW);
90  // For Egamma integration
91  nEgamma = 1000;
92  dEgamma = std::log(_targetMaxPhotonEnergy/_targetMinPhotonEnergy)/nEgamma;
93  minEgamma = std::log(_targetMinPhotonEnergy);
94 
95 
96  if (getBNORM() == 0.){
97  cout<<" Using Breit-Wigner Resonance Profile."<<endl;
98  }
99  else{
100  cout<<" Using Breit-Wigner plus direct pi+pi- profile."<<endl;
101  }
102 
103  cout<<" Integrating over W from "<<_wideWmin<<" to "<<_wideWmax<<endl;
104 
105  int A_1 = getbbs().electronBeam().A();
106  int A_2 = getbbs().targetBeam().A();
107 
108  if( A_2 == 0 && A_1 >= 1 ){
109  // eA, first beam is the nucleus and is in this case the target
110  beam = 1;
111  } else if( A_1 ==0 && A_2 >= 1){
112  // eA, second beam is the nucleus and is in this case the target
113  beam = 2;
114  } else {
115  beam = 2;
116  }
117 
118  int_r=0.;
119  int_r2 = 0.;
120  // Do this first for the case when the first beam is the photon emitter
121  // The variable beam (=1,2) defines which nucleus is the target
122  // Integration done using Simpson's rule
123  for(iW=0;iW<=nW-1;iW++){
124 
125  W = _wideWmin + double(iW)*dW + 0.5*dW;
126  int nQ2 = 1000;
127  for(iEgamma = 0 ; iEgamma < nEgamma; ++iEgamma){ // Integral over photon energy
128  // Target frame photon energies
129  ega[0] = exp(minEgamma + iEgamma*dEgamma );
130  ega[1] = exp(minEgamma + (iEgamma+1)*dEgamma );
131  ega[2] = 0.5*(ega[0]+ega[1]);
132  // Integral over Q2
133  double full_int[3] = {0}; // Full e+X --> e+X+V.M. cross section
134  double dndE[3] = {0}; // Full e+X --> e+X+V.M. cross section
135  //
136  for( int iEgaInt = 0 ; iEgaInt < 3; ++iEgaInt){ // Loop over the energies for the three points to integrate over Q2
137  //These are the physical limits
138  double Q2_min = std::pow(starlightConstants::mel*ega[iEgaInt],2.0)/(_electronEnergy*(_electronEnergy-ega[iEgaInt]));
139  double Q2_max = 4.*_electronEnergy*(_electronEnergy-ega[iEgaInt]);
140  if(_useFixedRange == true){
141  if( Q2_min < _gammaMinQ2 )
142  Q2_min = _gammaMinQ2;
143  if( Q2_max > _gammaMaxQ2 )
144  Q2_max = _gammaMaxQ2;
145  }
146  double lnQ2ratio = std::log(Q2_max/Q2_min)/nQ2;
147  double lnQ2_min = std::log(Q2_min);
148  //
149  for( int iQ2 = 0 ; iQ2 < nQ2; ++iQ2){ // Integral over photon virtuality
150  //
151  double q2_1 = exp( lnQ2_min + iQ2*lnQ2ratio);
152  double q2_2 = exp( lnQ2_min + (iQ2+1)*lnQ2ratio);
153  double q2_12 = (q2_2+q2_1)/2.;
154  //
155  // Integrating cross section
156  full_int[iEgaInt] += (q2_2-q2_1)*( g(ega[iEgaInt],q2_1)*getcsgA(ega[iEgaInt],q2_1,beam)
157  + g(ega[iEgaInt],q2_2)*getcsgA(ega[iEgaInt],q2_2,beam)
158  + 4.*g(ega[iEgaInt],q2_12)*getcsgA(ega[iEgaInt],q2_12,beam) );
159  // Effective flux
160  dndE[iEgaInt] +=(q2_2-q2_1)*( getcsgA_Q2_dep(q2_1)*photonFlux(ega[iEgaInt],q2_1)
161  +getcsgA_Q2_dep(q2_2)*photonFlux(ega[iEgaInt],q2_2)
162  +4.*getcsgA_Q2_dep(q2_12)*photonFlux(ega[iEgaInt],q2_12) );
163  }
164  full_int[iEgaInt] = full_int[iEgaInt]/6.;
165  dndE[iEgaInt] = dndE[iEgaInt]/6.;
166  }
167  // Finishing cross-section integral
168  dR = full_int[0];
169  dR += full_int[1];
170  dR += 4.*full_int[2];
171  dR = dR*(ega[1]-ega[0])/6.;
172  int_r = int_r + dR*breitWigner(W,bwnorm)*dW;
173  // Finishing effective flux integral
174  // Finishing integral over the effective photon flux
175  dR2 = dndE[0];
176  dR2 += dndE[1];
177  dR2 += 4.*dndE[2];
178  dR2 = dR2*(ega[1]-ega[0])/6.;
179  int_r2 = int_r2 + dR2*breitWigner(W,bwnorm)*dW;
180  }
181  }
182  cout<<endl;
183  if(_useFixedRange == true){
184  cout<<" Using fixed Q2 range "<<_gammaMinQ2 << " < Q2 < "<<_gammaMaxQ2<<endl;
185  }
186  printCrossSection(" Total cross section: ",int_r);
187  //printCrossSection(" gamma+X --> VM+X ", int_r/int_r2);
188  setPhotonNucleusSigma(0.01*int_r);
189  //
190 #ifdef _makeGammaPQ2_
191  makeGammaPQ2dependence(bwnormsave);
192 #endif
193 }
194 
195 #ifdef _makeGammaPQ2_
196 //______________________________________________________________________________
197 void
198 e_wideResonanceCrossSection::makeGammaPQ2dependence( double bwnormsave)
199 {
200  // This subroutine calculates the Q2 dependence of
201  // gamma+X -> VM + X cross section for a narrow resonance
202 
203  int const nQ2bins = 19;
204  double const q2Edge[nQ2bins+1] = { 0.,1.,2.,3., 4.,5.,
205  6.,7.,8.,9.,10.,
206  11.,12.,13.,14.,15.,
207  20.,30.,40.,50.};
208  double W = 0,dW,dY;
209  double y1,y2,y12,ega1,ega2,ega12;
210  double int_r,dR,dR2;
211  double csgA1, csgA2, csgA12;
212  double Eth;
213  int I,J,NW,NY,beam;
214 
215  double bwnorm = bwnormsave; //used to transfer the bwnorm from the luminosity tables
216 
217  NW = 100;
218  dW = (_wideWmax-_wideWmin)/double(NW);
219 
220  if (getBNORM() == 0.){
221  cout<<" Using Breit-Wigner Resonance Profile."<<endl;
222  }
223  else{
224  cout<<" Using Breit-Wigner plus direct pi+pi- profile."<<endl;
225  }
226 
227  cout<<" Integrating over W from "<<_wideWmin<<" to "<<_wideWmax<<endl;
228 
229  //Lomnitz old used for XX
230  Eth=0.5*(((W+protonMass)*(W+protonMass)-
231  protonMass*protonMass)/(_Ep+sqrt(_Ep*_Ep-protonMass*protonMass)));
232  // Adapted for eX
233  //Eth=0.5*(((W+starlightConstants::mel)*(W +starlightConstants::mel)-
234  // starlightConstants::mel*starlightConstants::mel)/(_electronEnergy+sqrt(_electronEnergy*_electronEnergy-starlightConstants::mel*starlightConstants::mel)));
235 
236  printf(" gamma+nucleon threshold: %e GeV \n", Eth);
237 
238  int A_1 = getbbs().electronBeam().A();
239  int A_2 = getbbs().targetBeam().A();
240 
241 
242  // Do this first for the case when the first beam is the photon emitter
243  // The variable beam (=1,2) defines which nucleus is the target
244  // Target beam ==2 so repidity is negative. Can generalize later
246  cout<<" Lomnitz debug :: sigma_gamma_p --> VM_p "<<endl;
247  cout<<" Q2+MV2 \t \t"<<" sigma_gamma_p --> VM_p (nanob)"<<endl;
248  double target_cm = acosh(_target_beamLorentz);
249  // another - sign from subraction in addition rule
250  double exp_target_cm = exp(-target_cm);
251  double int_r2;
252  for( int iQ2 = 0 ; iQ2 < nQ2bins; ++iQ2){
253  int_r=0.;
254  int_r2=0.;
255  //
256  double q2_cor = getcsgA_Q2_dep( (q2Edge[iQ2+1] + q2Edge[iQ2])/2. );
257  for(I=0;I<=NW-1;I++){
258 
259  W = _wideWmin + double(I)*dW + 0.5*dW;
260  for(J=0;J<=(NY-1);J++){
261  y1 = _wideYmin + double(J)*dY;
262  y2 = _wideYmin + double(J+1)*dY;
263  y12 = 0.5*(y1+y2);
264  double target_ega1, target_ega2, target_ega12;
265  if( A_2 == 0 && A_1 >= 1 ){
266  // pA, first beam is the nucleus and is in this case the target
267  ega1 = 0.5*W*exp(-y1);
268  ega2 = 0.5*W*exp(-y2);
269  ega12 = 0.5*W*exp(-y12);
270  beam = 1;
271  } else if( A_1 ==0 && A_2 >= 1){
272  // pA, second beam is the nucleus and is in this case the target
273  ega1 = 0.5*W*exp(y1);
274  ega2 = 0.5*W*exp(y2);
275  ega12 = 0.5*W*exp(y12);
276  // photon energy in Target frame
277  beam = 2;
278  } else {
279  ega1 = 0.5*W*exp(y1);
280  ega2 = 0.5*W*exp(y2);
281  ega12 = 0.5*W*exp(y12);
282  // photon energy in Target frame
283  beam = 2;
284  }
285  //
286  if(ega1 < Eth || ega2 < Eth)
287  continue;
288  if(ega2 > maxPhotonEnergy() || ega1 > maxPhotonEnergy() )
289  continue;
290  target_ega1 = ega1*exp_target_cm;
291  target_ega12 = ega12*exp_target_cm;
292  target_ega2 = ega2*exp_target_cm;
293  //cout<<"Nortmalizations "<<integrated_x_section(ega1,0,50)<<endl;
294  //
295  csgA1=getcsgA(ega1,W,beam);
296  double full_range_1 = integrated_x_section(target_ega1);
297  // >> Middle Point =====>>>
298  csgA12=getcsgA(ega12,W,beam);
299  double full_range_12 = integrated_x_section(target_ega12);
300  // >> Second Point =====>>>
301  csgA2=getcsgA(ega2,W,beam);
302  double full_range_2 = integrated_x_section(target_ega2);
303  //
304 
305  //>> Sum the contribution for this W,Y. The 2 accounts for the 2 beams
306  dR = q2_cor*csgA1;
307  dR = dR + 4.*q2_cor*csgA12;
308  dR = dR + q2_cor*csgA2;
309  dR = dR*(dY/6.)*breitWigner(W,bwnorm)*dW;
310  //
311  dR2 = full_range_1*csgA1;
312  dR2 = dR2 + 4.*full_range_12*csgA12;
313  dR2 = dR2 + full_range_2*csgA2;
314  dR2 = dR2*(dY/6.)*breitWigner(W,bwnorm)*dW;
315  //
316  int_r = int_r+dR;
317  int_r2 = int_r2 +dR2;
318  }
319  //
320  }
321  if( iQ2 ==0 )
322  cout<<"Full range "<<int_r2*10000000<<endl;
323  cout<<(q2Edge[iQ2+1]+q2Edge[iQ2])/2.+W*W<<" , "<<10000000.*int_r<<endl;
324  }
325 
326 }
327 #endif
328 void e_wideResonanceCrossSection::printCrossSection(const string name, const double x_section)
329 {
330  if (0.01*x_section > 1.){
331  cout<< name.c_str() <<0.01*x_section<<" barn."<<endl;
332  } else if (10.*x_section > 1.){
333  cout<< name.c_str() <<10.*x_section<<" mb."<<endl;
334  } else if (10000.*x_section > 1.){
335  cout<< name.c_str() <<10000.*x_section<<" microb."<<endl;
336  } else if (10000000.*x_section > 1.){
337  cout<< name.c_str() <<10000000.*x_section<<" nanob."<<endl;
338  } else if (1.E10*x_section > 1.){
339  cout<< name.c_str() <<1.E10*x_section<<" picob."<<endl;
340  } else {
341  cout<< name.c_str() <<1.E13*x_section<<" femtob."<<endl;
342  }
343  //cout<<endl;
344 }