EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
photonNucleusCrossSection.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file photonNucleusCrossSection.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:: 276 $: revision of last commit
24 // $Author:: jnystrand $: author of last commit
25 // $Date:: 2016-09-13 19:54:42 +0100 #$: date of last commit
26 //
27 // Description:
28 //
29 //
31 
32 
33 #include <iostream>
34 #include <fstream>
35 #include <cmath>
36 
37 #include "reportingUtils.h"
38 #include "starlightconstants.h"
39 #include "bessel.h"
40 #include "photonNucleusCrossSection.h"
41 
42 
43 using namespace std;
44 using namespace starlightConstants;
45 
46 
47 //______________________________________________________________________________
48 photonNucleusCrossSection::photonNucleusCrossSection(const inputParameters& inputParametersInstance, const beamBeamSystem& bbsystem)
49  : _nWbins (inputParametersInstance.nmbWBins() ),
50  _nYbins (inputParametersInstance.nmbRapidityBins() ),
51  _wMin (inputParametersInstance.minW() ),
52  _wMax (inputParametersInstance.maxW() ),
53  _yMax (inputParametersInstance.maxRapidity() ),
54  _beamLorentzGamma (inputParametersInstance.beamLorentzGamma() ),
55  _bbs (bbsystem ),
56  _protonEnergy (inputParametersInstance.protonEnergy() ),
57  _electronEnergy (inputParametersInstance.electronEnergy() ),
58  _particleType (inputParametersInstance.prodParticleType() ),
59  _beamBreakupMode (inputParametersInstance.beamBreakupMode() ),
60  _backwardsProduction(inputParametersInstance.backwardsProduction()),
61  _productionMode (inputParametersInstance.productionMode() ),
62  _sigmaNucleus (_bbs.targetBeam().A() ),
63  _fixedQ2range (inputParametersInstance.fixedQ2Range() ),
64  _minQ2 (inputParametersInstance.minGammaQ2() ),
65  _maxQ2 (inputParametersInstance.maxGammaQ2() ),
66  _maxPhotonEnergy (inputParametersInstance.cmsMaxPhotonEnergy()),
67  _cmsMinPhotonEnergy(inputParametersInstance.cmsMinPhotonEnergy()),
68  _targetRadii (inputParametersInstance.targetRadius() ),
69  _maxW_GP (inputParametersInstance.maxW_GP() ),
70  _minW_GP (inputParametersInstance.minW_GP() )
71 {
72  // new options - impulse aproximation (per Joakim) and Quantum Glauber (per SK) SKQG
73  _impulseSelected = inputParametersInstance.impulseVM();
74  _quantumGlauber = inputParametersInstance.quantumGlauber();
75  switch(_particleType) {
76  case RHO:
77  _slopeParameter = 11.0; // [(GeV/c)^{-2}]
78  _vmPhotonCoupling = 2.02;
79  _vmQ2Power_c1 = 2.09;
80  _vmQ2Power_c2 = 0.0073; // [ GeV^{-2}]
81  _ANORM = -2.75;
82  _BNORM = 0.0;
83  _defaultC = 1.0;
84  _channelMass = starlightConstants::rho0Mass;
85  _width = starlightConstants::rho0Width;
86  if(_backwardsProduction){
87  _slopeParameter = 21.8; // [(GeV/c)^{-2}]
88  _ANORM = 1.;
89  }
90  break;
91  case RHOZEUS:
92  _slopeParameter =11.0;
93  _vmPhotonCoupling=2.02;
94  _vmQ2Power_c1 = 2.09;
95  _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
96  _ANORM=-2.75;
97  _BNORM=1.84;
98  _defaultC=1.0;
99  _channelMass = starlightConstants::rho0Mass;
100  _width = starlightConstants::rho0Width;
101  break;
102  case FOURPRONG:
103  _slopeParameter = 11.0;
104  _vmPhotonCoupling = 2.02;
105  _vmQ2Power_c1 = 2.09;
106  _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
107  _ANORM = -2.75;
108  _BNORM = 0;
109  _defaultC = 11.0;
110  _channelMass = starlightConstants::rho0PrimeMass;
111  _width = starlightConstants::rho0PrimeWidth;
112  break;
113  case OMEGA:
114  case OMEGA_pi0gamma:
115  _slopeParameter=10.0;
116  _vmPhotonCoupling=23.13;
117  _vmQ2Power_c1 = 2.09;
118  _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
119  _ANORM=-2.75;
120  _BNORM=0.0;
121  _defaultC=1.0;
122  _channelMass = starlightConstants::OmegaMass;
123  _width = starlightConstants::OmegaWidth;
124  if(_backwardsProduction){
125  _slopeParameter = 21.8; // [(GeV/c)^{-2}]
126  _ANORM = 1.;
127  }
128  break;
129  case PHI:
130  _slopeParameter=7.0;
131  _vmPhotonCoupling=13.71;
132  _vmQ2Power_c1 = 2.15;
133  _vmQ2Power_c2 = 0.0074; // [GeV^{-2}]
134  _ANORM=-2.75;
135  _BNORM=0.0;
136  _defaultC=1.0;
137  _channelMass = starlightConstants::PhiMass;
138  _width = starlightConstants::PhiWidth;
139  break;
140  case JPSI:
141  case JPSI_ee:
142  _slopeParameter=4.0;
143  _vmPhotonCoupling=10.45;
144  _vmQ2Power_c1 = 2.45;
145  _vmQ2Power_c2 = 0.00084; // [GeV^{-2}]
146  _ANORM=-2.75;
147  _BNORM=0.0;
148  _defaultC=1.0;
149  _channelMass = starlightConstants::JpsiMass;
150  _width = starlightConstants::JpsiWidth;
151  break;
152  case JPSI_mumu:
153  _slopeParameter=4.0;
154  _vmPhotonCoupling=10.45;
155  _vmQ2Power_c1 = 2.36;
156  _vmQ2Power_c2 = 0.0029; // [GeV^{-2}]
157  _ANORM=-2.75;
158  _BNORM=0.0;
159  _defaultC=1.0;
160  _channelMass = starlightConstants::JpsiMass;
161  _width = starlightConstants::JpsiWidth;
162  break;
163  case JPSI2S:
164  case JPSI2S_ee:
165  case JPSI2S_mumu:
166  _slopeParameter=4.3;
167  _vmPhotonCoupling=26.39;
168  _vmQ2Power_c1 = 2.09;
169  _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
170  _ANORM=-2.75;
171  _BNORM=0.0;
172  _defaultC=1.0;
173  _channelMass = starlightConstants::Psi2SMass;
174  _width = starlightConstants::Psi2SWidth;
175  break;
176  case UPSILON:
177  case UPSILON_ee:
178  case UPSILON_mumu:
179  _slopeParameter=4.0;
180  _vmPhotonCoupling=125.37;
181  _vmQ2Power_c1 = 2.09;
182  _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
183  _ANORM=-2.75;
184  _BNORM=0.0;
185  _defaultC=1.0;
186  _channelMass = starlightConstants::Upsilon1SMass;
187  _width = starlightConstants::Upsilon1SWidth;
188  break;
189  case UPSILON2S:
190  case UPSILON2S_ee:
191  case UPSILON2S_mumu:
192  _slopeParameter=4.0;
193  _vmPhotonCoupling=290.84;
194  _vmQ2Power_c1 = 2.09;
195  _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
196  _ANORM=-2.75;
197  _BNORM=0.0;
198  _defaultC=1.0;
199  _channelMass = starlightConstants::Upsilon2SMass;
200  _width = starlightConstants::Upsilon2SWidth;
201  break;
202  case UPSILON3S:
203  case UPSILON3S_ee:
204  case UPSILON3S_mumu:
205  _slopeParameter=4.0;
206  _vmPhotonCoupling=415.10;
207  _vmQ2Power_c1 = 2.09;
208  _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
209  _ANORM=-2.75;
210  _BNORM=0.0;
211  _defaultC=1.0;
212  _channelMass = starlightConstants::Upsilon3SMass;
213  _width = starlightConstants::Upsilon3SWidth;
214  break;
215  default:
216  cout <<"No sigma constants parameterized for pid: "<<_particleType
217  <<" GammaAcrosssection"<<endl;
218  }
219 
220  //Changed by Lomnitz for e case. Limit is now E_e - 100m_e
221  //_maxPhotonEnergy = 12. * _beamLorentzGamma * hbarc/(_bbs.beam1().nuclearRadius()+_bbs.targetBeam().nuclearRadius());
222  //_maxPhotonEnergy = _electronEnergy - 10.*starlightConstants::mel;
223  /*cout<<" Lomnitz:: max energy in target frame "<< _electronEnergy - 1000.*starlightConstants::mel<<" vs electron energy "<<_electronEnergy<<endl
224  <<" max energy in cms frame "<<_maxPhotonEnergy<<" vs electron energy "<<_beamLorentzGamma*starlightConstants::mel<<endl;
225  cout<<" testing original limit "<< 12. * _beamLorentzGamma * hbarc/(2.*_bbs.targetBeam().nuclearRadius())<<endl;*/
226 
227 
228 }
229 
230 
231 //______________________________________________________________________________
232 photonNucleusCrossSection::~photonNucleusCrossSection()
233 { }
234 
235 
236 //______________________________________________________________________________
237 void
238 photonNucleusCrossSection::crossSectionCalculation(const double)
239 {
240  cout << "Neither narrow/wide resonance cross-section calculation.--Derived" << endl;
241 }
242 
243 //______________________________________________________________________________
244 double
245 photonNucleusCrossSection::getcsgA(const double targetEgamma,
246  const double Q2,
247  const int beam)
248 {
249  //This function returns the cross-section for photon-nucleus interaction
250  //producing vectormesons
251 
252  double Av,Wgp,cs,cvma;
253  double t,tmin,tmax;
254  double csgA,ax,bx;
255  int NGAUSS;
256  //
257  double W = _channelMass; //new change, channel mass center used for the t min integration.
258 
259  // DATA FOR GAUSS INTEGRATION
260  double xg[6] = {0, 0.1488743390, 0.4333953941, 0.6794095683, 0.8650633667, 0.9739065285};
261  double ag[6] = {0, 0.2955242247, 0.2692667193, 0.2190863625, 0.1494513492, 0.0666713443};
262  NGAUSS = 6;
263 
264  // Note: The photon energy passed to this function is now in the target frame. The rest of the calculations are done in the
265  // CMS frame. The next lines boost the photon into the CM frame.
266  double E_prime = _electronEnergy - targetEgamma;
267  double cos_theta_e = 1. - Q2/(2.*_electronEnergy*E_prime);
268  double theta_e = acos(cos_theta_e);
269  double beam_y = acosh(_beamLorentzGamma);
270  double gamma_pt = E_prime*sin(theta_e);
271  double pz_squared = targetEgamma*targetEgamma - Q2 - gamma_pt*gamma_pt;
272  if( pz_squared < 0 || fabs(cos_theta_e) > 1 || 2.*targetEgamma/(Q2+W*W) < _targetRadii)
273  return 0;
274  double temp_pz = sqrt(pz_squared);
275  // Now boost to CM frame
276  double Egamma = targetEgamma*cosh(beam_y) - temp_pz*sinh(beam_y);
277  if( Egamma < _cmsMinPhotonEnergy || Egamma > _maxPhotonEnergy){
278  return 0;
279  }
280  //cout<<" ::: Lomnitz test in photonNucleus ::: pz^2 = "<<pz_squared << " CMS Egamma = "<<Egamma<<endl;
281  // Find gamma-proton CM energy in CMS frame in the limit Q2->0 (this is our assumption, the Q2 dependence is in the factor)
282  Wgp = sqrt( 2.*(protonMass*targetEgamma)+protonMass*protonMass);
283  /*Wgp = sqrt(2. * Egamma * (_protonEnergy
284  + sqrt(_protonEnergy * _protonEnergy - protonMass * protonMass))
285  + protonMass * protonMass);*/
286 
287  //Used for A-A
288  tmin = (W * W / (4. * Egamma * _beamLorentzGamma)) * (W * W / (4. * Egamma * _beamLorentzGamma));
289 
290  if ((_bbs.electronBeam().A() <= 1) && (_bbs.targetBeam().A() <= 1)){
291  // proton-proton, no scaling needed
292  csgA = getcsgA_Q2_dep(Q2)*sigmagp(Wgp);
293  } else {
294  // Check if one or both beams are nuclei
295  int A_1 = _bbs.electronBeam().A();
296  int A_2 = _bbs.targetBeam().A();
297  // coherent AA interactions
298  // Calculate V.M.+proton cross section
299  // cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha);
300  cs = getcsgA_Q2_dep(Q2)*sigma_N(Wgp); //Use member function instead
301 
302  // Calculate V.M.+nucleus cross section
303  cvma = sigma_A(cs,beam);
304 
305  // Do impulse approximation here
306  if( _impulseSelected == 1){
307  if( beam == 1 ){
308  cvma = A_1*cs;
309  } else if ( beam == 2 ){
310  cvma = A_2*cs;
311  }
312  }
313 
314  // Calculate Av = dsigma/dt(t=0) Note Units: fm**s/Gev**2
315  Av = (alpha * cvma * cvma) / (16. * pi * _vmPhotonCoupling * hbarc * hbarc);
316 
317  tmax = tmin + 0.25;
318  ax = 0.5 * (tmax - tmin);
319  bx = 0.5 * (tmax + tmin);
320  csgA = 0.;
321  for (int k = 1; k < NGAUSS; ++k) {
322 
323  t = ax * xg[k] + bx;
324  if( A_1 <= 1 && A_2 != 1){
325  csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
326  }else if(A_2 <=1 && A_1 != 1){
327  csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
328  }else{
329  if( beam==1 ){
330  csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
331  }else if(beam==2){
332  csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
333  }else{
334  cout<<"Something went wrong here, beam= "<<beam<<endl;
335  }
336  }
337 
338  t = ax * (-xg[k]) + bx;
339  if( A_1 <= 1 && A_2 != 1){
340  csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
341  }else if(A_2 <=1 && A_1 != 1){
342  csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
343  }else{
344  if( beam==1 ){
345  csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
346  }else if(beam==2){
347  csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
348  }else{
349  cout<<"Something went wrong here, beam= "<<beam<<endl;
350  }
351  }
352  }
353  csgA = 0.5 * (tmax - tmin) * csgA;
354  csgA = Av * csgA;
355  }
356  return csgA;
357 }
358 
359 
360 //______________________________________________________________________________
361 double
362 photonNucleusCrossSection::e_getcsgA(const double Egamma, double Q2,
363  const double W,
364  const int beam)
365 {
366  //This function returns the cross-section for photon-nucleus interaction
367  //producing vectormesons for e_starlight. Stuff will be returned in CMS frame, but
368  //photon input is take in target frame
369 
370  double Av,Wgp,cs,cvma;
371  double t,tmin,tmax;
372  double csgA,ax,bx;
373  int NGAUSS;
374 
375  // DATA FOR GAUSS INTEGRATION
376  double xg[6] = {0, 0.1488743390, 0.4333953941, 0.6794095683, 0.8650633667, 0.9739065285};
377  double ag[6] = {0, 0.2955242247, 0.2692667193, 0.2190863625, 0.1494513492, 0.0666713443};
378  NGAUSS = 6;
379 
380  // Find gamma-proton CM energy
381  Wgp = sqrt( 2.*(protonMass*Egamma)
382  +protonMass*protonMass + Q2);
383 
384  //Used for A-A
385  tmin = (W * W / (4. * Egamma * _beamLorentzGamma)) * (W * W / (4. * Egamma * _beamLorentzGamma));
386 
387  if ((_bbs.electronBeam().A() <= 1) && (_bbs.targetBeam().A() <= 1)){
388  // proton-proton, no scaling needed
389  csgA = sigmagp(Wgp);
390  } else {
391  // coherent AA interactions
392  // Calculate V.M.+proton cross section
393  // cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha);
394  cs = sigma_N(Wgp); //Use member function instead
395 
396  // Calculate V.M.+nucleus cross section
397  cvma = sigma_A(cs,beam);
398 
399  // Calculate Av = dsigma/dt(t=0) Note Units: fm**s/Gev**2
400  Av = (alpha * cvma * cvma) / (16. * pi * _vmPhotonCoupling * hbarc * hbarc);
401 
402  // Check if one or both beams are nuclei
403  int A_1 = _bbs.electronBeam().A();
404  int A_2 = _bbs.targetBeam().A();
405 
406  tmax = tmin + 0.25;
407  ax = 0.5 * (tmax - tmin);
408  bx = 0.5 * (tmax + tmin);
409  csgA = 0.;
410  for (int k = 1; k < NGAUSS; ++k) {
411 
412  t = ax * xg[k] + bx;
413  if( A_1 <= 1 && A_2 != 1){
414  csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
415  }else if(A_2 <=1 && A_1 != 1){
416  csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
417  }else{
418  if( beam==1 ){
419  csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
420  }else if(beam==2){
421  csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
422  }else{
423  cout<<"Something went wrong here, beam= "<<beam<<endl;
424  }
425  }
426 
427  t = ax * (-xg[k]) + bx;
428  if( A_1 <= 1 && A_2 != 1){
429  csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
430  }else if(A_2 <=1 && A_1 != 1){
431  csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
432  }else{
433  if( beam==1 ){
434  csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
435  }else if(beam==2){
436  csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
437  }else{
438  cout<<"Something went wrong here, beam= "<<beam<<endl;
439  }
440  }
441  }
442  csgA = 0.5 * (tmax - tmin) * csgA;
443  csgA = Av * csgA;
444  }
445  return csgA;
446 }
447 
448 
449 //______________________________________________________________________________
450 double
451 photonNucleusCrossSection::getcsgA_Q2_dep(const double Q2)
452 {
453  double const mv2 = getChannelMass()*getChannelMass();
454  double const n = vmQ2Power(Q2);
455  return std::pow(mv2/(mv2+Q2),n);
456 }
457 
458 
459 //______________________________________________________________________________
460 double
461 photonNucleusCrossSection::photonFlux(const double Egamma, const int beam)
462 {
463  // This routine gives the photon flux as a function of energy Egamma
464  // It works for arbitrary nuclei and gamma; the first time it is
465  // called, it calculates a lookup table which is used on
466  // subsequent calls.
467  // It returns dN_gamma/dE (dimensions 1/E), not dI/dE
468  // energies are in GeV, in the lab frame
469  // rewritten 4/25/2001 by SRK
470 
471  // NOTE: beam (=1,2) defines the photon TARGET
472 
473  double lEgamma,Emin,Emax;
474  static double lnEmax, lnEmin, dlnE;
475  double stepmult,energy,rZ;
476  int nbstep,nrstep,nphistep,nstep;
477  double bmin,bmax,bmult,biter,bold,integratedflux;
478  double fluxelement,deltar,riter;
479  double deltaphi,phiiter,dist;
480  static double dide[401];
481  double lnElt;
482  double flux_r;
483  double Xvar;
484  int Ilt;
485  double RNuc=0.,RSum=0.;
486 
487  RSum=_bbs.electronBeam().nuclearRadius()+_bbs.targetBeam().nuclearRadius();
488  if( beam == 1){
489  rZ=double(_bbs.targetBeam().Z());
490  RNuc = _bbs.electronBeam().nuclearRadius();
491  } else {
492  rZ=double(_bbs.electronBeam().Z());
493  RNuc = _bbs.targetBeam().nuclearRadius();
494  }
495 
496  static int Icheck = 0;
497  static int Ibeam = 0;
498 
499  //Check first to see if pp
500  if( _bbs.electronBeam().A()==1 && _bbs.targetBeam().A()==1 ){
501  int nbsteps = 400;
502  double bmin = 0.5;
503  double bmax = 5.0 + (5.0*_beamLorentzGamma*hbarc/Egamma);
504  double dlnb = (log(bmax)-log(bmin))/(1.*nbsteps);
505 
506  double local_sum=0.0;
507 
508  // Impact parameter loop
509  for (int i = 0; i<=nbsteps;i++){
510 
511  double bnn0 = bmin*exp(i*dlnb);
512  double bnn1 = bmin*exp((i+1)*dlnb);
513  double db = bnn1-bnn0;
514 
515  double ppslope = 19.0;
516  double GammaProfile = exp(-bnn0*bnn0/(2.*hbarc*hbarc*ppslope));
517  double PofB0 = 1. - (1. - GammaProfile)*(1. - GammaProfile);
518  GammaProfile = exp(-bnn1*bnn1/(2.*hbarc*hbarc*ppslope));
519  double PofB1 = 1. - (1. - GammaProfile)*(1. - GammaProfile);
520 
521  double loc_nofe0 = _bbs.electronBeam().photonDensity(bnn0,Egamma);
522  double loc_nofe1 = _bbs.targetBeam().photonDensity(bnn1,Egamma);
523 
524  local_sum += 0.5*loc_nofe0*(1. - PofB0)*2.*starlightConstants::pi*bnn0*db;
525  local_sum += 0.5*loc_nofe1*(1. - PofB1)*2.*starlightConstants::pi*bnn1*db;
526 
527  }
528  // End Impact parameter loop
529  return local_sum;
530  }
531 
532  // first call or new beam? - initialize - calculate photon flux
533  Icheck=Icheck+1;
534 
535  // Do the numerical integration only once for symmetric systems.
536  if( Icheck > 1 && _bbs.electronBeam().A() == _bbs.targetBeam().A() && _bbs.electronBeam().Z() == _bbs.targetBeam().Z() ) goto L1000f;
537  // For asymmetric systems check if we have another beam
538  if( Icheck > 1 && beam == Ibeam ) goto L1000f;
539  Ibeam = beam;
540 
541  // Nuclear breakup is done by PofB
542  // collect number of integration steps here, in one place
543 
544  nbstep=1200;
545  nrstep=60;
546  nphistep=40;
547 
548  // this last one is the number of energy steps
549  nstep=100;
550 
551  // following previous choices, take Emin=10 keV at LHC, Emin = 1 MeV at RHIC
552  Emin=1.E-5;
553  if (_beamLorentzGamma < 500)
554  Emin=1.E-3;
555 
556  Emax=_maxPhotonEnergy;
557  // Emax=12.*hbarc*_beamLorentzGamma/RSum;
558 
559  // >> lnEmin <-> ln(Egamma) for the 0th bin
560  // >> lnEmax <-> ln(Egamma) for the last bin
561  lnEmin=log(Emin);
562  lnEmax=log(Emax);
563  dlnE=(lnEmax-lnEmin)/nstep;
564 
565  printf("Calculating photon flux from Emin = %e GeV to Emax = %e GeV (CM frame) for source with Z = %3.0f \n", Emin, Emax, rZ);
566 
567  stepmult= exp(log(Emax/Emin)/double(nstep));
568  energy=Emin;
569 
570  for (int j = 1; j<=nstep;j++){
571  energy=energy*stepmult;
572 
573  // integrate flux over 2R_A < b < 2R_A+ 6* gamma hbar/energy
574  // use exponential steps
575 
576  bmin=0.8*RSum; //Start slightly below 2*Radius
577  bmax=bmin + 6.*hbarc*_beamLorentzGamma/energy;
578 
579  bmult=exp(log(bmax/bmin)/double(nbstep));
580  biter=bmin;
581  integratedflux=0.;
582 
583  if( (_bbs.electronBeam().A() == 1 && _bbs.targetBeam().A() != 1) || (_bbs.targetBeam().A() == 1 && _bbs.electronBeam().A() != 1) ){
584  // This is pA
585 
586  if( _productionMode == PHOTONPOMERONINCOHERENT ){
587  // This pA incoherent, proton is the target
588 
589  int nbsteps = 400;
590  double bmin = 0.7*RSum;
591  double bmax = 2.0*RSum + (8.0*_beamLorentzGamma*hbarc/energy);
592  double dlnb = (log(bmax)-log(bmin))/(1.*nbsteps);
593 
594  double local_sum=0.0;
595 
596  // Impact parameter loop
597  for (int i = 0; i<=nbsteps; i++){
598 
599  double bnn0 = bmin*exp(i*dlnb);
600  double bnn1 = bmin*exp((i+1)*dlnb);
601  double db = bnn1-bnn0;
602 
603  double PofB0 = _bbs.probabilityOfBreakup(bnn0);
604  double PofB1 = _bbs.probabilityOfBreakup(bnn1);
605 
606  double loc_nofe0 = 0.0;
607  double loc_nofe1 = 0.0;
608  if( _bbs.electronBeam().A() == 1 ){
609  loc_nofe0 = _bbs.targetBeam().photonDensity(bnn0,energy);
610  loc_nofe1 = _bbs.targetBeam().photonDensity(bnn1,energy);
611  }
612  else if( _bbs.targetBeam().A() == 1 ){
613  loc_nofe0 = _bbs.electronBeam().photonDensity(bnn0,energy);
614  loc_nofe1 = _bbs.electronBeam().photonDensity(bnn1,energy);
615  }
616 
617  // cout<<" i: "<<i<<" bnn0: "<<bnn0<<" PofB0: "<<PofB0<<" loc_nofe0: "<<loc_nofe0<<endl;
618 
619  local_sum += 0.5*loc_nofe0*PofB0*2.*starlightConstants::pi*bnn0*db;
620  local_sum += 0.5*loc_nofe1*PofB1*2.*starlightConstants::pi*bnn1*db;
621  } // End Impact parameter loop
622  integratedflux = local_sum;
623  } else if ( _productionMode == PHOTONPOMERONNARROW || _productionMode == PHOTONPOMERONWIDE ){
624  // This is pA coherent, nucleus is the target
625  double localbmin = 0.0;
626  if( _bbs.electronBeam().A() == 1 ){
627  localbmin = _bbs.targetBeam().nuclearRadius() + 0.7;
628  }
629  if( _bbs.targetBeam().A() == 1 ){
630  localbmin = _bbs.electronBeam().nuclearRadius() + 0.7;
631  }
632  integratedflux = nepoint(energy,localbmin);
633  }
634  }else{
635  // This is AA
636  for (int jb = 1; jb<=nbstep;jb++){
637  bold=biter;
638  biter=biter*bmult;
639  // When we get to b>20R_A change methods - just take the photon flux
640  // at the center of the nucleus.
641  if (biter > (10.*RNuc)){
642  // if there is no nuclear breakup or only hadronic breakup, which only
643  // occurs at smaller b, we can analytically integrate the flux from b~20R_A
644  // to infinity, following Jackson (2nd edition), Eq. 15.54
645  Xvar=energy*biter/(hbarc*_beamLorentzGamma);
646  // Here, there is nuclear breakup. So, we can't use the integrated flux
647  // However, we can do a single flux calculation, at the center of the nucleus
648  // Eq. 41 of Vidovic, Greiner and Soff, Phys.Rev.C47,2308(1993), among other places
649  // this is the flux per unit area
650  fluxelement = (rZ*rZ*alpha*energy)*
651  (bessel::dbesk1(Xvar))*(bessel::dbesk1(Xvar))/
652  ((pi*_beamLorentzGamma*hbarc)*
653  (pi*_beamLorentzGamma*hbarc));
654 
655  } else {
656  // integrate over nuclear surface. n.b. this assumes total shadowing -
657  // treat photons hitting the nucleus the same no matter where they strike
658  fluxelement=0.;
659  deltar=RNuc/double(nrstep);
660  riter=-deltar/2.;
661 
662  for (int jr =1; jr<=nrstep;jr++){
663  riter=riter+deltar;
664  // use symmetry; only integrate from 0 to pi (half circle)
665  deltaphi=pi/double(nphistep);
666  phiiter=0.;
667 
668  for( int jphi=1;jphi<= nphistep;jphi++){
669  phiiter=(double(jphi)-0.5)*deltaphi;
670  // dist is the distance from the center of the emitting nucleus
671  // to the point in question
672  dist=sqrt((biter+riter*cos(phiiter))*(biter+riter*
673  cos(phiiter))+(riter*sin(phiiter))*(riter*sin(phiiter)));
674  Xvar=energy*dist/(hbarc*_beamLorentzGamma);
675  flux_r = (rZ*rZ*alpha*energy)*
676  (bessel::dbesk1(Xvar)*bessel::dbesk1(Xvar))/
677  ((pi*_beamLorentzGamma*hbarc)*
678  (pi*_beamLorentzGamma*hbarc));
679 
680  // The surface element is 2.* delta phi* r * delta r
681  // The '2' is because the phi integral only goes from 0 to pi
682  fluxelement=fluxelement+flux_r*2.*deltaphi*riter*deltar;
683  // end phi and r integrations
684  }//for(jphi)
685  }//for(jr)
686  // average fluxelement over the nuclear surface
687  fluxelement=fluxelement/(pi*RNuc*RNuc);
688  }//else
689  // multiply by volume element to get total flux in the volume element
690  fluxelement=fluxelement*2.*pi*biter*(biter-bold);
691  // modulate by the probability of nuclear breakup as f(biter)
692  // cout<<" jb: "<<jb<<" biter: "<<biter<<" fluxelement: "<<fluxelement<<endl;
693  if (_beamBreakupMode > 1){
694  fluxelement=fluxelement*_bbs.probabilityOfBreakup(biter);
695  }
696  // cout<<" jb: "<<jb<<" biter: "<<biter<<" fluxelement: "<<fluxelement<<endl;
697  integratedflux=integratedflux+fluxelement;
698 
699  } //end loop over impact parameter
700  } //end of else (pp, pA, AA)
701 
702  // In lookup table, store k*dN/dk because it changes less
703  // so the interpolation should be better
704  dide[j]=integratedflux*energy;
705  }//end loop over photon energy
706 
707  // for 2nd and subsequent calls, use lookup table immediately
708 
709  L1000f:
710 
711  lEgamma=log(Egamma);
712  if (lEgamma < (lnEmin+dlnE) || lEgamma > lnEmax){
713  flux_r=0.0;
714  cout<<" ERROR: Egamma outside defined range. Egamma= "<<Egamma
715  <<" "<<lnEmax<<" "<<(lnEmin+dlnE)<<endl;
716  }
717  else{
718  // >> Egamma between Ilt and Ilt+1
719  Ilt = int((lEgamma-lnEmin)/dlnE);
720  // >> ln(Egamma) for first point
721  lnElt = lnEmin + Ilt*dlnE;
722  // >> Interpolate
723  flux_r = dide[Ilt] + ((lEgamma-lnElt)/dlnE)*(dide[Ilt+1]- dide[Ilt]);
724  flux_r = flux_r/Egamma;
725  }
726 
727  return flux_r;
728 }
729 
730 
731 //______________________________________________________________________________
732 double
733 photonNucleusCrossSection::integrated_Q2_dep(double const Egamma, double const _min, double const _max)
734 {
735  //Integration over full limits gives more accurate result
736  double Q2_min = std::pow(starlightConstants::mel*Egamma,2.0)/(_electronEnergy*(_electronEnergy-Egamma));
737  double Q2_max = 4.*_electronEnergy*(_electronEnergy-Egamma);
738  //double Q2_max = 2.*Egamma/_targetRadii - _wMax*_wMax;
739  if( _min != 0 || _max !=0){
740  if( _min > Q2_min )
741  Q2_min = _min;
742  if( _max < Q2_max )
743  Q2_max = _max;
744  }
745  // Simpsons rule in using exponential step size and 10000 steps. Tested trapeve rule, linear steps and
746  // nstep = 10,000 100,000 & 10,000,0000
747  int nstep = 1000;
748  double ln_min = std::log(Q2_min);
749  double ratio = std::log(Q2_max/Q2_min)/nstep;
750  double g_int = 0;
751  //double g_int2 = 0 ;
752  //double g_int3 = 0;
753  //cout<<"*** Lomnitz **** Energy "<<Egamma<<" limits "<<Q2_min*1E9<<" x 1E-9 - "<<Q2_max<<endl;
754  for ( int ii = 0 ; ii< nstep; ++ii){
755  double x1 = std::exp(ln_min+(double)ii*ratio);
756  double x3 = std::exp(ln_min+(double)(ii+1)*ratio);
757  double x2 = (x3+x1)/2.;
758  //cout<<"ii : "<<x1<<" "<<x2<<" "<<x3<<endl;
759  g_int += (x3-x1)*( g(Egamma,x3)+g(Egamma,x1) +4.*g(Egamma,x2));
760  //g_int2 += (x3-x1)*( photonFlux(Egamma,x3)+photonFlux(Egamma,x1) +4.*photonFlux(Egamma,x2));
761  //g_int3 += (x3-x1)*( getcsgA_Q2_dep(x3)+getcsgA_Q2_dep(x1) +4.*getcsgA_Q2_dep(x2));
762  }
763  //return g_int2*g_int3/36.;
764  //return g_int2/6.;
765  return g_int/6.;
766 }
767 
768 
769 //______________________________________________________________________________
770 double
771 photonNucleusCrossSection::integrated_x_section(double const Egamma, double const _min, double const _max)
772 {
773  //Integration over full limits gives more accurate result
774  double Q2_min = std::pow(starlightConstants::mel*Egamma,2.0)/(_electronEnergy*(_electronEnergy-Egamma));
775  //double Q2_max = 2.*Egamma/_targetRadii - _wMax*_wMax;
776  double Q2_max = 4.*_electronEnergy*(_electronEnergy-Egamma);
777  // Fixed ranges for plot
778  if( _min != 0 || _max!=0){
779  if( _min > Q2_min)
780  Q2_min = _min;
781  if( _max < Q2_max)
782  Q2_max = _max;
783  }
784  // Simpsons rule in using exponential step size and 10000 steps. Tested trapeve rule, linear steps and
785  // nstep = 10,000 100,000 & 10,000,0000
786  int nstep = 1000;
787  double ln_min = std::log(Q2_min);
788  double ratio = std::log(Q2_max/Q2_min)/nstep;
789  double g_int = 0;
790  for ( int ii = 0 ; ii< nstep; ++ii){
791  double x1 = std::exp(ln_min+(double)ii*ratio);
792  double x3 = std::exp(ln_min+(double)(ii+1)*ratio);
793  double x2 = (x3+x1)/2.;
794  //Tests from HERA https://arxiv.org/pdf/hep-ex/9601009.pdf
795  // g_int += (x3-x1)*( getcsgA_Q2_dep(x3)+getcsgA_Q2_dep(x1) +4.*getcsgA_Q2_dep(x2));
796  g_int += (x3-x1)*( getcsgA_Q2_dep(x3)+getcsgA_Q2_dep(x1) +4.*getcsgA_Q2_dep(x2));
797  }
798  //return g_int2*g_int3/36.;
799  return g_int/6.;
800 }
801 
802 
803 //______________________________________________________________________________
804 pair< double, double >*photonNucleusCrossSection::Q2arraylimits(double const Egamma)
805 {
806  //double Q2max = 2.*Egamma/_targetRadii - _wMax*_wMax;
807  double Q2max = 4.*_electronEnergy*(_electronEnergy-Egamma);
808  double Q2min= std::pow(starlightConstants::mel*Egamma,2.0)/(_electronEnergy*(_electronEnergy-Egamma));
809 
810  if( _fixedQ2range == true){
811  if( Q2min < _minQ2 )
812  Q2min = _minQ2;
813  if( Q2max > _maxQ2 )
814  Q2max = _maxQ2;
815  //cout<<" Imposed limits "<<Q2min<<" - "<<Q2max<<endl;
816  std::pair<double,double>* to_ret = new std::pair<double, double>(Q2min,Q2max);
817  return to_ret;
818  }
819  int Nstep = 1000;
820  //
821  double ratio = std::log(Q2max/Q2min)/(double)Nstep;
822  double ln_min = std::log(Q2min);
823  // -- - -
824  const double limit = 1E-9;
825  std::vector<double>Q2_array;
826  int iNstep = 0;
827  double g_step = 1.;
828  while( g_step>limit ){
829  double Q2 = std::exp(ln_min+iNstep*ratio);
830  if(Q2>Q2max)
831  break;
832  g_step = g(Egamma,Q2);
833  Q2_array.push_back(g_step);
834  iNstep++;
835  }
836  if( std::exp(ln_min+iNstep*ratio) < Q2max)
837  Q2max = std::exp(ln_min+iNstep*ratio);
838  //cout<<Q2max<<" "<<g(Egamma,Q2max)*1E9<<endl;
839  std::pair<double, double>* to_ret = new std::pair<double, double>(Q2min,Q2max);
840 
841  return to_ret;
842 }
843 
844 
845 //______________________________________________________________________________
846 double
847 photonNucleusCrossSection::g(double const Egamma,
848  double const Q2)
849 {
850  //return photonFlux(Egamma,Q2)*getcsgA_Q2_dep(Q2);
851  return photonFlux(Egamma,Q2); //This could be done more elegantly in the future, but this one change should account for everything
852 }
853 
854 //______________________________________________________________________________
855 double
856 photonNucleusCrossSection::photonFlux(double const Egamma,
857  double const Q2)
858 {
859  //Need to check units later
860  //double const hbar = starlightConstants::hbarc / 2.99*pow(10,14); // 6.582 * pow (10,-16) [eVs]
861  //double omega = Egamma/ hbar;
862  //Do we even need a lookup table for this case? This should return N(E,Q2) from dn = N(E,Q2) dE dQ2
863  double const ratio = Egamma/_electronEnergy;
864  double const minQ2 = std::pow( starlightConstants::mel*Egamma,2.0) / (_electronEnergy*(_electronEnergy - Egamma));
865  double to_ret = alpha/(pi) *( 1- ratio + ratio*ratio/2. - (1-ratio)*( fabs(minQ2/Q2)) );
866  //Comparisons:
867  // double temp = pow(2.*_electronEnergy-Egamma,2.)/(Egamma*Egamma + Q2) + 1. - 4.*starlightConstants::mel*starlightConstants::mel/Q2;
868  //temp = alpha*temp*Egamma/(4.*Q2*pi*_electronEnergy*_electronEnergy);
869  // cout<<" *** Lomnitz *** Testing photon flux approximation for electron energy "<<_electronEnergy<<" gamma "<<Egamma<<" Q2 "<<Q2*1E6<<" 1E-6 "<<endl;
870  //cout<<" Full expression "<<temp*1E6<<" vs. apporioximation "<<to_ret/( Egamma*fabs(Q2) )*1E6<<" --- ratio "<<temp/(to_ret/( Egamma*fabs(Q2) ))<<endl;
871  return to_ret/( Egamma*fabs(Q2) );
872  //return temp;
873 }
874 
875 //______________________________________________________________________________
876 double
877 photonNucleusCrossSection::nepoint(const double Egamma,
878  const double bmin)
879 {
880  // Function for the spectrum of virtual photons,
881  // dn/dEgamma, for a point charge q=Ze sweeping
882  // past the origin with velocity gamma
883  // (=1/SQRT(1-(V/c)**2)) integrated over impact
884  // parameter from bmin to infinity
885  // See Jackson eq15.54 Classical Electrodynamics
886  // Declare Local Variables
887  double beta,X,C1,bracket,nepoint_r;
888 
889  beta = sqrt(1.-(1./(_beamLorentzGamma*_beamLorentzGamma)));
890  X = (bmin*Egamma)/(beta*_beamLorentzGamma*hbarc);
891 
892  bracket = -0.5*beta*beta*X*X*(bessel::dbesk1(X)*bessel::dbesk1(X)
893  -bessel::dbesk0(X)*bessel::dbesk0(X));
894 
895  bracket = bracket+X*bessel::dbesk0(X)*bessel::dbesk1(X);
896 
897  // Note: NO Z*Z!!
898  C1=(2.*alpha)/pi;
899 
900  nepoint_r = C1*(1./beta)*(1./beta)*(1./Egamma)*bracket;
901 
902  return nepoint_r;
903 
904 }
905 
906 
907 //______________________________________________________________________________
908 double
909 photonNucleusCrossSection::sigmagp(const double Wgp)
910 {
911  // Function for the gamma-proton --> VectorMeson
912  // cross section. Wgp is the gamma-proton CM energy.
913  // Unit for cross section: fm**2
914  //If W_gp is not in the allowed region, i.e. W_GP_MIN < W_gp < W_GP_MAX, do not sample.
915  if(Wgp < _minW_GP || Wgp > _maxW_GP) return 0;
916  double sigmagp_r=0.;
917 
918  // Near the threshold CM energy (between WgpMin and WgpMax),
919  // we add a linear scaling factor to bring the Xsec down to 0
920  // at the threshold CM value define by WgpMin = m_p + m_vm
921  double WgpMax = 0.;
922  double WgpMin = 0.;
923  double thresholdScaling = 1.0;
924 
925  switch(_particleType)
926  {
927  case RHO:
928  case RHOZEUS:
929  case FOURPRONG:
930  WgpMax = 1.8;
931  WgpMin = 1.60; //this is the cutoff threshold for rho production. But rho has large width so it's lower
932  if(Wgp<WgpMax) thresholdScaling=(Wgp-WgpMin)/(WgpMax-WgpMin);
933  sigmagp_r=thresholdScaling*1.E-4*(5.0*exp(0.22*log(Wgp))+26.0*exp(-1.23*log(Wgp)));
934  //This is based on the omega cross section:
935  //https://arxiv.org/pdf/2107.06748.pdf
936  if(_backwardsProduction)sigmagp_r=thresholdScaling*1.E-4*0.14*pow(Wgp,-2.7);
937  if(_backwardsProduction)sigmagp_r=thresholdScaling*1.E-4*0.206*pow(((Wgp*Wgp-protonMass*protonMass)/(2.0*protonMass)),-2.7);
938  break;
939  case OMEGA:
940  case OMEGA_pi0gamma:
941  WgpMax = 1.8;
942  WgpMin = 1.74; //this is the cutoff threshold for omega production: W > m_p+m_omega = 1.74
943  if(Wgp<WgpMax) thresholdScaling=(Wgp-WgpMin)/(WgpMax-WgpMin);
944  sigmagp_r=thresholdScaling*1.E-4*(0.55*exp(0.22*log(Wgp))+18.0*exp(-1.92*log(Wgp)));
945  //if(_backwardsProduction)sigmagp_r=thresholdScaling*1.E-4*0.14*pow(Wgp,-2.7);//https://arxiv.org/pdf/2107.06748.pdf
946  if(_backwardsProduction)sigmagp_r=thresholdScaling*1.E-4*0.206*pow(((Wgp*Wgp-protonMass*protonMass)/(2.0*protonMass)),-2.7);
947  break;
948  case PHI:
949  sigmagp_r=1.E-4*0.34*exp(0.22*log(Wgp));
950  break;
951  case JPSI:
952  case JPSI_ee:
953  case JPSI_mumu:
954  sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
955  sigmagp_r*=sigmagp_r;
956  sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp));
957  // sigmagp_r=1.E-4*0.0015*exp(0.80*log(Wgp));
958  break;
959  case JPSI2S:
960  case JPSI2S_ee:
961  case JPSI2S_mumu:
962  sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
963  sigmagp_r*=sigmagp_r;
964  sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp));
965  sigmagp_r*=0.166;
966  // sigmagp_r=0.166*(1.E-4*0.0015*exp(0.80*log(Wgp)));
967  break;
968  case UPSILON:
969  case UPSILON_ee:
970  case UPSILON_mumu:
971  // >> This is W**1.7 dependence from QCD calculations
972  // sigmagp_r=1.E-10*(0.060)*exp(1.70*log(Wgp));
973  sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
974  sigmagp_r*=sigmagp_r;
975  sigmagp_r*=1.E-10*6.4*exp(0.74*log(Wgp));
976  break;
977  case UPSILON2S:
978  case UPSILON2S_ee:
979  case UPSILON2S_mumu:
980  // sigmagp_r=1.E-10*(0.0259)*exp(1.70*log(Wgp));
981  sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
982  sigmagp_r*=sigmagp_r;
983  sigmagp_r*=1.E-10*2.9*exp(0.74*log(Wgp));
984  break;
985  case UPSILON3S:
986  case UPSILON3S_ee:
987  case UPSILON3S_mumu:
988  // sigmagp_r=1.E-10*(0.0181)*exp(1.70*log(Wgp));
989  sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
990  sigmagp_r*=sigmagp_r;
991  sigmagp_r*=1.E-10*2.1*exp(0.74*log(Wgp));
992  break;
993  default: cout<< "!!! ERROR: Unidentified Vector Meson: "<< _particleType <<endl;
994  }
995  return sigmagp_r;
996 }
997 
998 
999 //______________________________________________________________________________
1000 double
1001 photonNucleusCrossSection::sigma_A(const double sig_N, const int beam)
1002 {
1003  // Nuclear Cross Section
1004  // sig_N,sigma_A in (fm**2)
1005 
1006  double sum;
1007  double b,bmax,Pint,arg,sigma_A_r;
1008 
1009  int NGAUSS;
1010 
1011  double xg[17]=
1012  {.0,
1013  .0483076656877383162,.144471961582796493,
1014  .239287362252137075, .331868602282127650,
1015  .421351276130635345, .506899908932229390,
1016  .587715757240762329, .663044266930215201,
1017  .732182118740289680, .794483795967942407,
1018  .849367613732569970, .896321155766052124,
1019  .934906075937739689, .964762255587506430,
1020  .985611511545268335, .997263861849481564
1021  };
1022 
1023  double ag[17]=
1024  {.0,
1025  .0965400885147278006, .0956387200792748594,
1026  .0938443990808045654, .0911738786957638847,
1027  .0876520930044038111, .0833119242269467552,
1028  .0781938957870703065, .0723457941088485062,
1029  .0658222227763618468, .0586840934785355471,
1030  .0509980592623761762, .0428358980222266807,
1031  .0342738629130214331, .0253920653092620595,
1032  .0162743947309056706, .00701861000947009660
1033  };
1034 
1035  NGAUSS=16;
1036 
1037  // Check if one or both beams are nuclei
1038  int A_1 = _bbs.electronBeam().A();
1039  int A_2 = _bbs.targetBeam().A();
1040  if( A_1 == 1 && A_2 == 1)cout<<" This is pp, you should not be here..."<<endl;
1041 
1042  // CALCULATE P(int) FOR b=0.0 - bmax (fm)
1043  bmax = 25.0;
1044  sum = 0.;
1045  for(int IB=1;IB<=NGAUSS;IB++){
1046 
1047  b = 0.5*bmax*xg[IB]+0.5*bmax;
1048 
1049  if( A_1 == 1 && A_2 != 1){
1050  arg=-sig_N*_bbs.targetBeam().rho0()*_bbs.targetBeam().thickness(b);
1051  }else if(A_2 == 1 && A_1 != 1){
1052  arg=-sig_N*_bbs.electronBeam().rho0()*_bbs.electronBeam().thickness(b);
1053  }else{
1054  // Check which beam is target
1055  if( beam == 1 ){
1056  arg=-sig_N*_bbs.electronBeam().rho0()*_bbs.electronBeam().thickness(b);
1057  }else if( beam==2 ){
1058  arg=-sig_N*_bbs.targetBeam().rho0()*_bbs.targetBeam().thickness(b);
1059  }else{
1060  cout<<" Something went wrong here, beam= "<<beam<<endl;
1061  }
1062  }
1063 
1064  Pint=1.0-exp(arg);
1065  // If this is a quantum Glauber calculation, use the quantum Glauber formula
1066  if (_quantumGlauber == 1){Pint=2.0*(1.0-exp(arg/2.0));}
1067 
1068  sum=sum+2.*pi*b*Pint*ag[IB];
1069 
1070  b = 0.5*bmax*(-xg[IB])+0.5*bmax;
1071 
1072  if( A_1 == 1 && A_2 != 1){
1073  arg=-sig_N*_bbs.targetBeam().rho0()*_bbs.targetBeam().thickness(b);
1074  }else if(A_2 == 1 && A_1 != 1){
1075  arg=-sig_N*_bbs.electronBeam().rho0()*_bbs.electronBeam().thickness(b);
1076  }else{
1077  // Check which beam is target
1078  if( beam == 1 ){
1079  arg=-sig_N*_bbs.electronBeam().rho0()*_bbs.electronBeam().thickness(b);
1080  }else if(beam==2){
1081  arg=-sig_N*_bbs.targetBeam().rho0()*_bbs.targetBeam().thickness(b);
1082  }else{
1083  cout<<" Something went wrong here, beam= "<<beam<<endl;
1084  }
1085  }
1086 
1087  Pint=1.0-exp(arg);
1088  // If this is a quantum Glauber calculation, use the quantum Glauber formula
1089  if (_quantumGlauber == 1){Pint=2.0*(1.0-exp(arg/2.0));}
1090 
1091  sum=sum+2.*pi*b*Pint*ag[IB];
1092 
1093  }
1094 
1095  sum=0.5*bmax*sum;
1096 
1097  sigma_A_r=sum;
1098 
1099  return sigma_A_r;
1100 }
1101 
1102 //______________________________________________________________________________
1103 double
1104 photonNucleusCrossSection::sigma_N(const double Wgp)
1105 {
1106  // Nucleon Cross Section in (fm**2)
1107  double cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha);
1108  return cs;
1109 }
1110 
1111 
1112 //______________________________________________________________________________
1113 double
1114 photonNucleusCrossSection::breitWigner(const double W,
1115  const double C)
1116 {
1117  // use simple fixed-width s-wave Breit-Wigner without coherent backgorund for rho'
1118  // (PDG '08 eq. 38.56)
1119  if(_particleType==FOURPRONG) {
1120  if (W < 4.01 * pionChargedMass)
1121  return 0;
1122  const double termA = _channelMass * _width;
1123  const double termA2 = termA * termA;
1124  const double termB = W * W - _channelMass * _channelMass;
1125  return C * _ANORM * _ANORM * termA2 / (termB * termB + termA2);
1126  }
1127 
1128  // Relativistic Breit-Wigner according to J.D. Jackson,
1129  // Nuovo Cimento 34, 6692 (1964), with nonresonant term. A is the strength
1130  // of the resonant term and b the strength of the non-resonant
1131  // term. C is an overall normalization.
1132 
1133  double ppi=0.,ppi0=0.,GammaPrim,rat;
1134  double aa,bb,cc;
1135 
1136  double nrbw_r;
1137 
1138  // width depends on energy - Jackson Eq. A.2
1139  // if below threshold, then return 0. Added 5/3/2001 SRK
1140  // 0.5% extra added for safety margin
1141  // omega added here 10/26/2014 SRK
1142  if( _particleType==RHO ||_particleType==RHOZEUS || _particleType==OMEGA){
1143  if (W < 2.01*pionChargedMass){
1144  nrbw_r=0.;
1145  return nrbw_r;
1146  }
1147  ppi=sqrt( ((W/2.)*(W/2.)) - pionChargedMass * pionChargedMass);
1148  ppi0=0.358;
1149  }
1150  if( _particleType==OMEGA_pi0gamma){
1151  if (W < pionNeutralMass){
1152  nrbw_r=0.;
1153  return nrbw_r;
1154  }
1155  ppi=abs((W*W - pionNeutralMass * pionNeutralMass)/(2.0*W));
1156  ppi0=abs((_channelMass*_channelMass - pionNeutralMass * pionNeutralMass)/(2.0*W));
1157  }
1158 
1159  // handle phi-->K+K- properly
1160  if (_particleType == PHI){
1161  if (W < 2.*kaonChargedMass){
1162  nrbw_r=0.;
1163  return nrbw_r;
1164  }
1165  ppi=sqrt( ((W/2.)*(W/2.))- kaonChargedMass*kaonChargedMass);
1166  ppi0=sqrt( ((_channelMass/2.)*(_channelMass/2.))-kaonChargedMass*kaonChargedMass);
1167  }
1168 
1169  //handle J/Psi-->e+e- properly
1170  if (_particleType==JPSI || _particleType==JPSI2S){
1171  if(W<2.*mel){
1172  nrbw_r=0.;
1173  return nrbw_r;
1174  }
1175  ppi=sqrt(((W/2.)*(W/2.))-mel*mel);
1176  ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel);
1177  }
1178  if (_particleType==JPSI_ee){
1179  if(W<2.*mel){
1180  nrbw_r=0.;
1181  return nrbw_r;
1182  }
1183  ppi=sqrt(((W/2.)*(W/2.))-mel*mel);
1184  ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel);
1185  }
1186  if (_particleType==JPSI_mumu){
1187  if(W<2.*muonMass){
1188  nrbw_r=0.;
1189  return nrbw_r;
1190  }
1191  ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass);
1192  ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass);
1193  }
1194  if (_particleType==JPSI2S_ee){
1195  if(W<2.*mel){
1196  nrbw_r=0.;
1197  return nrbw_r;
1198  }
1199  ppi=sqrt(((W/2.)*(W/2.))-mel*mel);
1200  ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel);
1201  }
1202  if (_particleType==JPSI2S_mumu){
1203  if(W<2.*muonMass){
1204  nrbw_r=0.;
1205  return nrbw_r;
1206  }
1207  ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass);
1208  ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass);
1209  }
1210 
1211  if(_particleType==UPSILON || _particleType==UPSILON2S ||_particleType==UPSILON3S ){
1212  if (W<2.*muonMass){
1213  nrbw_r=0.;
1214  return nrbw_r;
1215  }
1216  ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass);
1217  ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass);
1218  }
1219 
1220  if(_particleType==UPSILON_mumu || _particleType==UPSILON2S_mumu ||_particleType==UPSILON3S_mumu ){
1221  if (W<2.*muonMass){
1222  nrbw_r=0.;
1223  return nrbw_r;
1224  }
1225  ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass);
1226  ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass);
1227  }
1228 
1229  if(_particleType==UPSILON_ee || _particleType==UPSILON2S_ee ||_particleType==UPSILON3S_ee ){
1230  if (W<2.*mel){
1231  nrbw_r=0.;
1232  return nrbw_r;
1233  }
1234  ppi=sqrt(((W/2.)*(W/2.))-mel*mel);
1235  ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel);
1236  }
1237 
1238  if(ppi==0.&&ppi0==0.)
1239  cout<<"Improper Gammaacrosssection::breitwigner, ppi&ppi0=0."<<endl;
1240 
1241  rat=ppi/ppi0;
1242  GammaPrim=_width*(_channelMass/W)*rat*rat*rat;
1243 
1244  aa=_ANORM*sqrt(GammaPrim*_channelMass*W);
1245  bb=W*W-_channelMass*_channelMass;
1246  cc=_channelMass*GammaPrim;
1247 
1248  // First real part squared
1249  nrbw_r = (( (aa*bb)/(bb*bb+cc*cc) + _BNORM)*( (aa*bb)/(bb*bb+cc*cc) + _BNORM));
1250 
1251  // Then imaginary part squared
1252  nrbw_r = nrbw_r + (( (aa*cc)/(bb*bb+cc*cc) )*( (aa*cc)/(bb*bb+cc*cc) ));
1253 
1254  // Alternative, a simple, no-background BW, following J. Breitweg et al.
1255  // Eq. 15 of Eur. Phys. J. C2, 247 (1998). SRK 11/10/2000
1256  // nrbw_r = (_ANORM*_mass*GammaPrim/(bb*bb+cc*cc))**2
1257 
1258 
1259  nrbw_r = C*nrbw_r;
1260 
1261  return nrbw_r;
1262 }