estarlight/src/photonNucleusCrossSection.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:: 276                         $: revision of last commit
0024 // $Author:: jnystrand                $: author of last commit
0025 // $Date:: 2016-09-13 19:54:42 +0100 #$: date of last commit
0026 //
0027 // Description:
0028 //
0029 //
0030 ///////////////////////////////////////////////////////////////////////////
0031
0032 #include <iostream>
0033 #include <fstream>
0034 #include <cmath>
0035
0036 #include "reportingUtils.h"
0037 #include "starlightconstants.h"
0038 #include "bessel.h"
0039 #include "photonNucleusCrossSection.h"
0040
0041 using namespace std;
0042 using namespace starlightConstants;
0043
0044 //______________________________________________________________________________
0045 photonNucleusCrossSection::photonNucleusCrossSection(const inputParameters &inputParametersInstance, const beamBeamSystem &bbsystem)
0046     : _nWbins(inputParametersInstance.nmbWBins()),
0047       _nYbins(inputParametersInstance.nmbRapidityBins()),
0048       _wMin(inputParametersInstance.minW()),
0049       _wMax(inputParametersInstance.maxW()),
0050       _yMax(inputParametersInstance.maxRapidity()),
0051       _beamLorentzGamma(inputParametersInstance.beamLorentzGamma()),
0052       _bbs(bbsystem),
0053       _protonEnergy(inputParametersInstance.protonEnergy()),
0054       _electronEnergy(inputParametersInstance.electronEnergy()),
0055       _particleType(inputParametersInstance.prodParticleType()),
0056       _beamBreakupMode(inputParametersInstance.beamBreakupMode()),
0057       _backwardsProduction(inputParametersInstance.backwardsProduction()),
0058       _rho0PrimeBrFourProng(inputParametersInstance.rho0PrimeBrFourProng()),
0059       _productionMode(inputParametersInstance.productionMode()),
0060       _sigmaNucleus(_bbs.targetBeam().A()),
0061       _fixedQ2range(inputParametersInstance.fixedQ2Range()),
0062       _minQ2(inputParametersInstance.minGammaQ2()),
0063       _maxQ2(inputParametersInstance.maxGammaQ2()),
0064       _maxPhotonEnergy(inputParametersInstance.cmsMaxPhotonEnergy()),
0065       _cmsMinPhotonEnergy(inputParametersInstance.cmsMinPhotonEnergy()),
0066       _targetRadii(inputParametersInstance.targetRadius()),
0067       _maxW_GP(inputParametersInstance.maxW_GP()),
0068       _minW_GP(inputParametersInstance.minW_GP())
0069 {
0070     // new options - impulse aproximation (per Joakim) and Quantum Glauber (per SK) SKQG
0071     _impulseSelected = inputParametersInstance.impulseVM();
0072     _quantumGlauber = inputParametersInstance.quantumGlauber();
0073     switch (_particleType)
0074     {
0075     case RHO:
0076         _slopeParameter = 11.0; // [(GeV/c)^{-2}]
0077         _vmPhotonCoupling = 2.02;
0078         _vmQ2Power_c1 = 2.09;
0079         _vmQ2Power_c2 = 0.0073; // [ GeV^{-2}]
0080         _ANORM = -2.75;
0081         _BNORM = 0.0;
0082         _defaultC = 1.0;
0083         _channelMass = starlightConstants::rho0Mass;
0084         _width = starlightConstants::rho0Width;
0085         if (_backwardsProduction)
0086         {
0087             _slopeParameter = 21.8; // [(GeV/c)^{-2}]
0088             _ANORM = 1.;
0089         }
0090         break;
0091     case RHOZEUS:
0092         _slopeParameter = 11.0;
0093         _vmPhotonCoupling = 2.02;
0094         _vmQ2Power_c1 = 2.09;
0095         _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
0096         _ANORM = -2.75;
0097         _BNORM = 1.84;
0098         _defaultC = 1.0;
0099         _channelMass = starlightConstants::rho0Mass;
0100         _width = starlightConstants::rho0Width;
0101         break;
0102     case FOURPRONG:
0103         _slopeParameter = 9.4;
0104         _vmPhotonCoupling = inputParametersInstance.rho0PrimeCoupling();
0105         _vmQ2Power_c1 = 2.09;
0106         _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
0107         _ANORM = -2.75;
0108         _BNORM = 0;
0109         _defaultC = 1.0;
0110         _channelMass = starlightConstants::rho0PrimeMass;
0111         _width = starlightConstants::rho0PrimeWidth;
0112         break;
0113     case OMEGA:
0114     case OMEGA_pi0gamma:
0115         _slopeParameter = 10.0;
0116         _vmPhotonCoupling = 23.13;
0117         _vmQ2Power_c1 = 2.09;
0118         _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
0119         _ANORM = -2.75;
0120         _BNORM = 0.0;
0121         _defaultC = 1.0;
0122         _channelMass = starlightConstants::OmegaMass;
0123         _width = starlightConstants::OmegaWidth;
0124         if (_backwardsProduction)
0125         {
0126             _slopeParameter = 21.8; // [(GeV/c)^{-2}]
0127             _ANORM = 1.;
0128         }
0129         break;
0130     case PHI:
0131         _slopeParameter = 7.0;
0132         _vmPhotonCoupling = 13.71;
0133         _vmQ2Power_c1 = 2.15;
0134         _vmQ2Power_c2 = 0.0074; // [GeV^{-2}]
0135         _ANORM = -2.75;
0136         _BNORM = 0.0;
0137         _defaultC = 1.0;
0138         _channelMass = starlightConstants::PhiMass;
0139         _width = starlightConstants::PhiWidth;
0140         break;
0141     case JPSI:
0142     case JPSI_ee:
0143         _slopeParameter = 4.0;
0144         _vmPhotonCoupling = 10.45;
0145         _vmQ2Power_c1 = 2.45;
0146         _vmQ2Power_c2 = 0.00084; // [GeV^{-2}]
0147         _ANORM = -2.75;
0148         _BNORM = 0.0;
0149         _defaultC = 1.0;
0150         _channelMass = starlightConstants::JpsiMass;
0151         _width = starlightConstants::JpsiWidth;
0152         break;
0153     case JPSI_mumu:
0154         _slopeParameter = 4.0;
0155         _vmPhotonCoupling = 10.45;
0156         _vmQ2Power_c1 = 2.36;
0157         _vmQ2Power_c2 = 0.0029; // [GeV^{-2}]
0158         _ANORM = -2.75;
0159         _BNORM = 0.0;
0160         _defaultC = 1.0;
0161         _channelMass = starlightConstants::JpsiMass;
0162         _width = starlightConstants::JpsiWidth;
0163         break;
0164     case JPSI2S:
0165     case JPSI2S_ee:
0166     case JPSI2S_mumu:
0167         _slopeParameter = 4.3;
0168         _vmPhotonCoupling = 26.39;
0169         _vmQ2Power_c1 = 2.09;
0170         _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
0171         _ANORM = -2.75;
0172         _BNORM = 0.0;
0173         _defaultC = 1.0;
0174         _channelMass = starlightConstants::Psi2SMass;
0175         _width = starlightConstants::Psi2SWidth;
0176         break;
0177     case UPSILON:
0178     case UPSILON_ee:
0179     case UPSILON_mumu:
0180         _slopeParameter = 4.0;
0181         _vmPhotonCoupling = 125.37;
0182         _vmQ2Power_c1 = 2.09;
0183         _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
0184         _ANORM = -2.75;
0185         _BNORM = 0.0;
0186         _defaultC = 1.0;
0187         _channelMass = starlightConstants::Upsilon1SMass;
0188         _width = starlightConstants::Upsilon1SWidth;
0189         break;
0190     case UPSILON2S:
0191     case UPSILON2S_ee:
0192     case UPSILON2S_mumu:
0193         _slopeParameter = 4.0;
0194         _vmPhotonCoupling = 290.84;
0195         _vmQ2Power_c1 = 2.09;
0196         _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
0197         _ANORM = -2.75;
0198         _BNORM = 0.0;
0199         _defaultC = 1.0;
0200         _channelMass = starlightConstants::Upsilon2SMass;
0201         _width = starlightConstants::Upsilon2SWidth;
0202         break;
0203     case UPSILON3S:
0204     case UPSILON3S_ee:
0205     case UPSILON3S_mumu:
0206         _slopeParameter = 4.0;
0207         _vmPhotonCoupling = 415.10;
0208         _vmQ2Power_c1 = 2.09;
0209         _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
0210         _ANORM = -2.75;
0211         _BNORM = 0.0;
0212         _defaultC = 1.0;
0213         _channelMass = starlightConstants::Upsilon3SMass;
0214         _width = starlightConstants::Upsilon3SWidth;
0215         break;
0216     default:
0217         cout << "No sigma constants parameterized for pid: " << _particleType
0218              << " GammaAcrosssection" << endl;
0219     }
0220
0221     // Changed by Lomnitz for e case. Limit is now E_e - 100m_e
0222     //_maxPhotonEnergy = 12. * _beamLorentzGamma * hbarc/(_bbs.beam1().nuclearRadius()+_bbs.targetBeam().nuclearRadius());
0223     //_maxPhotonEnergy = _electronEnergy - 10.*starlightConstants::mel;
0224     /*cout<<" Lomnitz:: max energy in target frame "<< _electronEnergy - 1000.*starlightConstants::mel<<" vs electron energy "<<_electronEnergy<<endl
0225         <<"           max energy in cms frame    "<<_maxPhotonEnergy<<"  vs electron energy "<<_beamLorentzGamma*starlightConstants::mel<<endl;
0226         cout<<" testing original limit "<< 12. * _beamLorentzGamma * hbarc/(2.*_bbs.targetBeam().nuclearRadius())<<endl;*/
0227 }
0228
0229 //______________________________________________________________________________
0230 photonNucleusCrossSection::~photonNucleusCrossSection()
0231 {
0232 }
0233
0234 //______________________________________________________________________________
0235 void photonNucleusCrossSection::crossSectionCalculation(const double)
0236 {
0237     cout << "Neither narrow/wide resonance cross-section calculation.--Derived" << endl;
0238 }
0239
0240 //______________________________________________________________________________
0241 double
0242 photonNucleusCrossSection::getcsgA(const double targetEgamma,
0243                                    const double Q2,
0244                                    const int beam)
0245 {
0246     // This function returns the cross-section for photon-nucleus interaction
0247     // producing vectormesons
0248
0249     double Av, Wgp, cs, cvma;
0250     double t, tmin, tmax;
0251     double csgA, ax, bx;
0252     int NGAUSS;
0253     //
0254     double W = _channelMass; // new change, channel mass center used for the t min integration.
0255
0256     //     DATA FOR GAUSS INTEGRATION
0257     double xg[6] = {0, 0.1488743390, 0.4333953941, 0.6794095683, 0.8650633667, 0.9739065285};
0258     double ag[6] = {0, 0.2955242247, 0.2692667193, 0.2190863625, 0.1494513492, 0.0666713443};
0259     NGAUSS = 6;
0260
0261     //       Note: The photon energy passed to this function is now in the target frame. The rest of the calculations are done in the
0262     //       CMS frame. The next lines boost the photon into the CM frame.
0263     double E_prime = _electronEnergy - targetEgamma;
0264     double cos_theta_e = 1. - Q2 / (2. * _electronEnergy * E_prime);
0265     double theta_e = acos(cos_theta_e);
0266     double beam_y = acosh(_beamLorentzGamma);
0267     double gamma_pt = E_prime * sin(theta_e);
0268     double pz_squared = targetEgamma * targetEgamma - Q2 - gamma_pt * gamma_pt;
0269     if (pz_squared < 0 || fabs(cos_theta_e) > 1 || 2. * targetEgamma / (Q2 + W * W) < _targetRadii)
0270         return 0;
0271     double temp_pz = sqrt(pz_squared);
0272     // Now boost to CM frame
0273     double Egamma = targetEgamma * cosh(beam_y) - temp_pz * sinh(beam_y);
0274     if (Egamma < _cmsMinPhotonEnergy || Egamma > _maxPhotonEnergy)
0275     {
0276         return 0;
0277     }
0278     // cout<<" ::: Lomnitz test in photonNucleus ::: pz^2 = "<<pz_squared << " CMS Egamma = "<<Egamma<<endl;
0279     //        Find gamma-proton CM energy in CMS frame in the limit Q2->0 (this is our assumption, the Q2 dependence is in the factor)
0280     Wgp = sqrt(2. * (protonMass * targetEgamma) + protonMass * protonMass);
0281     /*Wgp = sqrt(2. * Egamma * (_protonEnergy
0282                               + sqrt(_protonEnergy * _protonEnergy - protonMass * protonMass))
0283                   + protonMass * protonMass);*/
0284
0285     // Used for A-A
0286     tmin = (W * W / (4. * Egamma * _beamLorentzGamma)) * (W * W / (4. * Egamma * _beamLorentzGamma));
0287
0288     if ((_bbs.electronBeam().A() <= 1) && (_bbs.targetBeam().A() <= 1))
0289     {
0290         // proton-proton, no scaling needed
0291         csgA = getcsgA_Q2_dep(Q2) * sigmagp(Wgp);
0292     }
0293     else
0294     {
0295         // Check if one or both beams are nuclei
0296         int A_1 = _bbs.electronBeam().A();
0297         int A_2 = _bbs.targetBeam().A();
0298         // coherent AA interactions
0299         // Calculate V.M.+proton cross section
0300         // cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha);
0301         cs = getcsgA_Q2_dep(Q2) * sigma_N(Wgp); // Use member function instead
0302
0303         // Calculate V.M.+nucleus cross section
0304         cvma = sigma_A(cs, beam);
0305
0306         // Do impulse approximation here
0307         if (_impulseSelected == 1)
0308         {
0309             if (beam == 1)
0310             {
0311                 cvma = A_1 * cs;
0312             }
0313             else if (beam == 2)
0314             {
0315                 cvma = A_2 * cs;
0316             }
0317         }
0318
0319         // Calculate Av = dsigma/dt(t=0) Note Units: fm**s/Gev**2
0320         Av = (alpha * cvma * cvma) / (16. * pi * _vmPhotonCoupling * hbarc * hbarc);
0321
0322         tmax = tmin + 0.25;
0323         ax = 0.5 * (tmax - tmin);
0324         bx = 0.5 * (tmax + tmin);
0325         csgA = 0.;
0326         for (int k = 1; k < NGAUSS; ++k)
0327         {
0328
0329             t = ax * xg[k] + bx;
0330             if (A_1 <= 1 && A_2 != 1)
0331             {
0332                 csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
0333             }
0334             else if (A_2 <= 1 && A_1 != 1)
0335             {
0336                 csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0337             }
0338             else
0339             {
0340                 if (beam == 1)
0341                 {
0342                     csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0343                 }
0344                 else if (beam == 2)
0345                 {
0346                     csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
0347                 }
0348                 else
0349                 {
0350                     cout << "Something went wrong here, beam= " << beam << endl;
0351                 }
0352             }
0353
0354             t = ax * (-xg[k]) + bx;
0355             if (A_1 <= 1 && A_2 != 1)
0356             {
0357                 csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
0358             }
0359             else if (A_2 <= 1 && A_1 != 1)
0360             {
0361                 csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0362             }
0363             else
0364             {
0365                 if (beam == 1)
0366                 {
0367                     csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0368                 }
0369                 else if (beam == 2)
0370                 {
0371                     csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
0372                 }
0373                 else
0374                 {
0375                     cout << "Something went wrong here, beam= " << beam << endl;
0376                 }
0377             }
0378         }
0379         csgA = 0.5 * (tmax - tmin) * csgA;
0380         csgA = Av * csgA;
0381     }
0382     return csgA;
0383 }
0384
0385 //______________________________________________________________________________
0386 double
0387 photonNucleusCrossSection::e_getcsgA(const double Egamma, double Q2,
0388                                      const double W,
0389                                      const int beam)
0390 {
0391     // This function returns the cross-section for photon-nucleus interaction
0392     // producing vectormesons for e_starlight. Stuff will be returned in CMS frame, but
0393     // photon input is take in target frame
0394
0395     double Av, Wgp, cs, cvma;
0396     double t, tmin, tmax;
0397     double csgA, ax, bx;
0398     int NGAUSS;
0399
0400     //     DATA FOR GAUSS INTEGRATION
0401     double xg[6] = {0, 0.1488743390, 0.4333953941, 0.6794095683, 0.8650633667, 0.9739065285};
0402     double ag[6] = {0, 0.2955242247, 0.2692667193, 0.2190863625, 0.1494513492, 0.0666713443};
0403     NGAUSS = 6;
0404
0405     //       Find gamma-proton CM energy
0406     Wgp = sqrt(2. * (protonMass * Egamma) + protonMass * protonMass + Q2);
0407
0408     // Used for A-A
0409     tmin = (W * W / (4. * Egamma * _beamLorentzGamma)) * (W * W / (4. * Egamma * _beamLorentzGamma));
0410
0411     if ((_bbs.electronBeam().A() <= 1) && (_bbs.targetBeam().A() <= 1))
0412     {
0413         // proton-proton, no scaling needed
0414         csgA = sigmagp(Wgp);
0415     }
0416     else
0417     {
0418         // coherent AA interactions
0419         // Calculate V.M.+proton cross section
0420         // cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha);
0421         cs = sigma_N(Wgp); // Use member function instead
0422
0423         // Calculate V.M.+nucleus cross section
0424         cvma = sigma_A(cs, beam);
0425
0426         // Calculate Av = dsigma/dt(t=0) Note Units: fm**s/Gev**2
0427         Av = (alpha * cvma * cvma) / (16. * pi * _vmPhotonCoupling * hbarc * hbarc);
0428
0429         // Check if one or both beams are nuclei
0430         int A_1 = _bbs.electronBeam().A();
0431         int A_2 = _bbs.targetBeam().A();
0432
0433         tmax = tmin + 0.25;
0434         ax = 0.5 * (tmax - tmin);
0435         bx = 0.5 * (tmax + tmin);
0436         csgA = 0.;
0437         for (int k = 1; k < NGAUSS; ++k)
0438         {
0439
0440             t = ax * xg[k] + bx;
0441             if (A_1 <= 1 && A_2 != 1)
0442             {
0443                 csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
0444             }
0445             else if (A_2 <= 1 && A_1 != 1)
0446             {
0447                 csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0448             }
0449             else
0450             {
0451                 if (beam == 1)
0452                 {
0453                     csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0454                 }
0455                 else if (beam == 2)
0456                 {
0457                     csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
0458                 }
0459                 else
0460                 {
0461                     cout << "Something went wrong here, beam= " << beam << endl;
0462                 }
0463             }
0464
0465             t = ax * (-xg[k]) + bx;
0466             if (A_1 <= 1 && A_2 != 1)
0467             {
0468                 csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
0469             }
0470             else if (A_2 <= 1 && A_1 != 1)
0471             {
0472                 csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0473             }
0474             else
0475             {
0476                 if (beam == 1)
0477                 {
0478                     csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0479                 }
0480                 else if (beam == 2)
0481                 {
0482                     csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
0483                 }
0484                 else
0485                 {
0486                     cout << "Something went wrong here, beam= " << beam << endl;
0487                 }
0488             }
0489         }
0490         csgA = 0.5 * (tmax - tmin) * csgA;
0491         csgA = Av * csgA;
0492     }
0493     return csgA;
0494 }
0495
0496 //______________________________________________________________________________
0497 double
0498 photonNucleusCrossSection::getcsgA_Q2_dep(const double Q2)
0499 {
0500     double const mv2 = getChannelMass() * getChannelMass();
0501     double const n = vmQ2Power(Q2);
0502     return std::pow(mv2 / (mv2 + Q2), n);
0503 }
0504
0505 //______________________________________________________________________________
0506 double
0507 photonNucleusCrossSection::photonFlux(const double Egamma, const int beam)
0508 {
0509     // This routine gives the photon flux as a function of energy Egamma
0510     // It works for arbitrary nuclei and gamma; the first time it is
0511     // called, it calculates a lookup table which is used on
0512     // subsequent calls.
0513     // It returns dN_gamma/dE (dimensions 1/E), not dI/dE
0514     // energies are in GeV, in the lab frame
0515     // rewritten 4/25/2001 by SRK
0516
0517     // NOTE: beam (=1,2) defines the photon TARGET
0518
0519     double lEgamma, Emin, Emax;
0520     static double lnEmax, lnEmin, dlnE;
0521     double stepmult, energy, rZ;
0522     int nbstep, nrstep, nphistep, nstep;
0523     double bmin, bmax, bmult, biter, bold, integratedflux;
0524     double fluxelement, deltar, riter;
0525     double deltaphi, phiiter, dist;
0526     static double dide[401];
0527     double lnElt;
0528     double flux_r;
0529     double Xvar;
0530     int Ilt;
0531     double RNuc = 0., RSum = 0.;
0532
0533     RSum = _bbs.electronBeam().nuclearRadius() + _bbs.targetBeam().nuclearRadius();
0534     if (beam == 1)
0535     {
0536         rZ = double(_bbs.targetBeam().Z());
0537         RNuc = _bbs.electronBeam().nuclearRadius();
0538     }
0539     else
0540     {
0541         rZ = double(_bbs.electronBeam().Z());
0542         RNuc = _bbs.targetBeam().nuclearRadius();
0543     }
0544
0545     static int Icheck = 0;
0546     static int Ibeam = 0;
0547
0548     // Check first to see if pp
0549     if (_bbs.electronBeam().A() == 1 && _bbs.targetBeam().A() == 1)
0550     {
0551         int nbsteps = 400;
0552         double bmin = 0.5;
0553         double bmax = 5.0 + (5.0 * _beamLorentzGamma * hbarc / Egamma);
0554         double dlnb = (log(bmax) - log(bmin)) / (1. * nbsteps);
0555
0556         double local_sum = 0.0;
0557
0558         // Impact parameter loop
0559         for (int i = 0; i <= nbsteps; i++)
0560         {
0561
0562             double bnn0 = bmin * exp(i * dlnb);
0563             double bnn1 = bmin * exp((i + 1) * dlnb);
0564             double db = bnn1 - bnn0;
0565
0566             double ppslope = 19.0;
0567             double GammaProfile = exp(-bnn0 * bnn0 / (2. * hbarc * hbarc * ppslope));
0568             double PofB0 = 1. - (1. - GammaProfile) * (1. - GammaProfile);
0569             GammaProfile = exp(-bnn1 * bnn1 / (2. * hbarc * hbarc * ppslope));
0570             double PofB1 = 1. - (1. - GammaProfile) * (1. - GammaProfile);
0571
0572             double loc_nofe0 = _bbs.electronBeam().photonDensity(bnn0, Egamma);
0573             double loc_nofe1 = _bbs.targetBeam().photonDensity(bnn1, Egamma);
0574
0575             local_sum += 0.5 * loc_nofe0 * (1. - PofB0) * 2. * starlightConstants::pi * bnn0 * db;
0576             local_sum += 0.5 * loc_nofe1 * (1. - PofB1) * 2. * starlightConstants::pi * bnn1 * db;
0577         }
0578         // End Impact parameter loop
0579         return local_sum;
0580     }
0581
0582     //   first call or new beam?  - initialize - calculate photon flux
0583     Icheck = Icheck + 1;
0584
0585     // Do the numerical integration only once for symmetric systems.
0586     if (Icheck > 1 && _bbs.electronBeam().A() == _bbs.targetBeam().A() && _bbs.electronBeam().Z() == _bbs.targetBeam().Z())
0587         goto L1000f;
0588     // For asymmetric systems check if we have another beam
0589     if (Icheck > 1 && beam == Ibeam)
0590         goto L1000f;
0591     Ibeam = beam;
0592
0593     //  Nuclear breakup is done by PofB
0594     //  collect number of integration steps here, in one place
0595
0596     nbstep = 1200;
0597     nrstep = 60;
0598     nphistep = 40;
0599
0600     //  this last one is the number of energy steps
0601     nstep = 100;
0602
0603     // following previous choices, take Emin=10 keV at LHC, Emin = 1 MeV at RHIC
0604     Emin = 1.E-5;
0605     if (_beamLorentzGamma < 500)
0606         Emin = 1.E-3;
0607
0608     Emax = _maxPhotonEnergy;
0609     // Emax=12.*hbarc*_beamLorentzGamma/RSum;
0610
0611     //     >> lnEmin <-> ln(Egamma) for the 0th bin
0612     //     >> lnEmax <-> ln(Egamma) for the last bin
0613     lnEmin = log(Emin);
0614     lnEmax = log(Emax);
0615     dlnE = (lnEmax - lnEmin) / nstep;
0616
0617     printf("Calculating photon flux from Emin = %e GeV to Emax = %e GeV (CM frame) for source with Z = %3.0f \n", Emin, Emax, rZ);
0618
0619     stepmult = exp(log(Emax / Emin) / double(nstep));
0620     energy = Emin;
0621
0622     for (int j = 1; j <= nstep; j++)
0623     {
0624         energy = energy * stepmult;
0625
0626         //  integrate flux over 2R_A < b < 2R_A+ 6* gamma hbar/energy
0627         //  use exponential steps
0628
0629         bmin = 0.8 * RSum; // Start slightly below 2*Radius
0630         bmax = bmin + 6. * hbarc * _beamLorentzGamma / energy;
0631
0632         bmult = exp(log(bmax / bmin) / double(nbstep));
0633         biter = bmin;
0634         integratedflux = 0.;
0635
0636         if ((_bbs.electronBeam().A() == 1 && _bbs.targetBeam().A() != 1) || (_bbs.targetBeam().A() == 1 && _bbs.electronBeam().A() != 1))
0637         {
0638             // This is pA
0639
0640             if (_productionMode == PHOTONPOMERONINCOHERENT)
0641             {
0642                 // This pA incoherent, proton is the target
0643
0644                 int nbsteps = 400;
0645                 double bmin = 0.7 * RSum;
0646                 double bmax = 2.0 * RSum + (8.0 * _beamLorentzGamma * hbarc / energy);
0647                 double dlnb = (log(bmax) - log(bmin)) / (1. * nbsteps);
0648
0649                 double local_sum = 0.0;
0650
0651                 // Impact parameter loop
0652                 for (int i = 0; i <= nbsteps; i++)
0653                 {
0654
0655                     double bnn0 = bmin * exp(i * dlnb);
0656                     double bnn1 = bmin * exp((i + 1) * dlnb);
0657                     double db = bnn1 - bnn0;
0658
0659                     double PofB0 = _bbs.probabilityOfBreakup(bnn0);
0660                     double PofB1 = _bbs.probabilityOfBreakup(bnn1);
0661
0662                     double loc_nofe0 = 0.0;
0663                     double loc_nofe1 = 0.0;
0664                     if (_bbs.electronBeam().A() == 1)
0665                     {
0666                         loc_nofe0 = _bbs.targetBeam().photonDensity(bnn0, energy);
0667                         loc_nofe1 = _bbs.targetBeam().photonDensity(bnn1, energy);
0668                     }
0669                     else if (_bbs.targetBeam().A() == 1)
0670                     {
0671                         loc_nofe0 = _bbs.electronBeam().photonDensity(bnn0, energy);
0672                         loc_nofe1 = _bbs.electronBeam().photonDensity(bnn1, energy);
0673                     }
0674
0675                     // cout<<" i: "<<i<<" bnn0: "<<bnn0<<" PofB0: "<<PofB0<<" loc_nofe0: "<<loc_nofe0<<endl;
0676
0677                     local_sum += 0.5 * loc_nofe0 * PofB0 * 2. * starlightConstants::pi * bnn0 * db;
0678                     local_sum += 0.5 * loc_nofe1 * PofB1 * 2. * starlightConstants::pi * bnn1 * db;
0679                 } // End Impact parameter loop
0680                 integratedflux = local_sum;
0681             }
0682             else if (_productionMode == PHOTONPOMERONNARROW || _productionMode == PHOTONPOMERONWIDE)
0683             {
0684                 // This is pA coherent, nucleus is the target
0685                 double localbmin = 0.0;
0686                 if (_bbs.electronBeam().A() == 1)
0687                 {
0688                     localbmin = _bbs.targetBeam().nuclearRadius() + 0.7;
0689                 }
0690                 if (_bbs.targetBeam().A() == 1)
0691                 {
0692                     localbmin = _bbs.electronBeam().nuclearRadius() + 0.7;
0693                 }
0694                 integratedflux = nepoint(energy, localbmin);
0695             }
0696         }
0697         else
0698         {
0699             // This is AA
0700             for (int jb = 1; jb <= nbstep; jb++)
0701             {
0702                 bold = biter;
0703                 biter = biter * bmult;
0704                 // When we get to b>20R_A change methods - just take the photon flux
0705                 //  at the center of the nucleus.
0706                 if (biter > (10. * RNuc))
0707                 {
0708                     // if there is no nuclear breakup or only hadronic breakup, which only
0709                     // occurs at smaller b, we can analytically integrate the flux from b~20R_A
0710                     // to infinity, following Jackson (2nd edition), Eq. 15.54
0711                     Xvar = energy * biter / (hbarc * _beamLorentzGamma);
0712                     // Here, there is nuclear breakup.  So, we can't use the integrated flux
0713                     // However, we can do a single flux calculation, at the center of the nucleus
0714                     // Eq. 41 of Vidovic, Greiner and Soff, Phys.Rev.C47,2308(1993), among other places
0715                     // this is the flux per unit area
0716                     fluxelement = (rZ * rZ * alpha * energy) *
0717                                   (bessel::dbesk1(Xvar)) * (bessel::dbesk1(Xvar)) /
0718                                   ((pi * _beamLorentzGamma * hbarc) *
0719                                    (pi * _beamLorentzGamma * hbarc));
0720                 }
0721                 else
0722                 {
0723                     // integrate over nuclear surface. n.b. this assumes total shadowing -
0724                     // treat photons hitting the nucleus the same no matter where they strike
0725                     fluxelement = 0.;
0726                     deltar = RNuc / double(nrstep);
0727                     riter = -deltar / 2.;
0728
0729                     for (int jr = 1; jr <= nrstep; jr++)
0730                     {
0731                         riter = riter + deltar;
0732                         // use symmetry;  only integrate from 0 to pi (half circle)
0733                         deltaphi = pi / double(nphistep);
0734                         phiiter = 0.;
0735
0736                         for (int jphi = 1; jphi <= nphistep; jphi++)
0737                         {
0738                             phiiter = (double(jphi) - 0.5) * deltaphi;
0739                             // dist is the distance from the center of the emitting nucleus
0740                             // to the point in question
0741                             dist = sqrt((biter + riter * cos(phiiter)) * (biter + riter *
0742                                                                                       cos(phiiter)) +
0743                                         (riter * sin(phiiter)) * (riter * sin(phiiter)));
0744                             Xvar = energy * dist / (hbarc * _beamLorentzGamma);
0745                             flux_r = (rZ * rZ * alpha * energy) *
0746                                      (bessel::dbesk1(Xvar) * bessel::dbesk1(Xvar)) /
0747                                      ((pi * _beamLorentzGamma * hbarc) *
0748                                       (pi * _beamLorentzGamma * hbarc));
0749
0750                             //  The surface  element is 2.* delta phi* r * delta r
0751                             //  The '2' is because the phi integral only goes from 0 to pi
0752                             fluxelement = fluxelement + flux_r * 2. * deltaphi * riter * deltar;
0753                             //  end phi and r integrations
0754                         } // for(jphi)
0755                     } // for(jr)
0756                     //  average fluxelement over the nuclear surface
0757                     fluxelement = fluxelement / (pi * RNuc * RNuc);
0758                 } // else
0759                 //  multiply by volume element to get total flux in the volume element
0760                 fluxelement = fluxelement * 2. * pi * biter * (biter - bold);
0761                 //  modulate by the probability of nuclear breakup as f(biter)
0762                 // cout<<" jb: "<<jb<<" biter: "<<biter<<" fluxelement: "<<fluxelement<<endl;
0763                 if (_beamBreakupMode > 1)
0764                 {
0765                     fluxelement = fluxelement * _bbs.probabilityOfBreakup(biter);
0766                 }
0767                 // cout<<" jb: "<<jb<<" biter: "<<biter<<" fluxelement: "<<fluxelement<<endl;
0768                 integratedflux = integratedflux + fluxelement;
0769
0770             } // end loop over impact parameter
0771         } // end of else (pp, pA, AA)
0772
0773         //  In lookup table, store k*dN/dk because it changes less
0774         //  so the interpolation should be better
0775         dide[j] = integratedflux * energy;
0776     } // end loop over photon energy
0777
0778     //  for 2nd and subsequent calls, use lookup table immediately
0779
0780 L1000f:
0781
0782     lEgamma = log(Egamma);
0783     if (lEgamma < (lnEmin + dlnE) || lEgamma > lnEmax)
0784     {
0785         flux_r = 0.0;
0786         cout << "  ERROR: Egamma outside defined range. Egamma= " << Egamma
0787              << "   " << lnEmax << " " << (lnEmin + dlnE) << endl;
0788     }
0789     else
0790     {
0791         //       >> Egamma between Ilt and Ilt+1
0792         Ilt = int((lEgamma - lnEmin) / dlnE);
0793         //       >> ln(Egamma) for first point
0794         lnElt = lnEmin + Ilt * dlnE;
0795         //       >> Interpolate
0796         flux_r = dide[Ilt] + ((lEgamma - lnElt) / dlnE) * (dide[Ilt + 1] - dide[Ilt]);
0797         flux_r = flux_r / Egamma;
0798     }
0799
0800     return flux_r;
0801 }
0802
0803 //______________________________________________________________________________
0804 double
0805 photonNucleusCrossSection::integrated_Q2_dep(double const Egamma, double const _min, double const _max)
0806 {
0807     // Integration over full limits gives more accurate result
0808     double Q2_min = std::pow(starlightConstants::mel * Egamma, 2.0) / (_electronEnergy * (_electronEnergy - Egamma));
0809     double Q2_max = 4. * _electronEnergy * (_electronEnergy - Egamma);
0810     // double Q2_max = 2.*Egamma/_targetRadii - _wMax*_wMax;
0811     if (_min != 0 || _max != 0)
0812     {
0813         if (_min > Q2_min)
0814             Q2_min = _min;
0815         if (_max < Q2_max)
0816             Q2_max = _max;
0817     }
0818     // Simpsons rule in using exponential step size and 10000 steps. Tested trapeve rule, linear steps and
0819     // nstep = 10,000 100,000 & 10,000,0000
0820     int nstep = 1000;
0821     double ln_min = std::log(Q2_min);
0822     double ratio = std::log(Q2_max / Q2_min) / nstep;
0823     double g_int = 0;
0824     // double g_int2 = 0 ;
0825     // double g_int3 = 0;
0826     // cout<<"*** Lomnitz **** Energy "<<Egamma<<" limits "<<Q2_min*1E9<<" x 1E-9 -  "<<Q2_max<<endl;
0827     for (int ii = 0; ii < nstep; ++ii)
0828     {
0829         double x1 = std::exp(ln_min + (double)ii * ratio);
0830         double x3 = std::exp(ln_min + (double)(ii + 1) * ratio);
0831         double x2 = (x3 + x1) / 2.;
0832         // cout<<"ii : "<<x1<<" "<<x2<<" "<<x3<<endl;
0833         g_int += (x3 - x1) * (g(Egamma, x3) + g(Egamma, x1) + 4. * g(Egamma, x2));
0834         // g_int2 += (x3-x1)*( photonFlux(Egamma,x3)+photonFlux(Egamma,x1) +4.*photonFlux(Egamma,x2));
0835         // g_int3 += (x3-x1)*( getcsgA_Q2_dep(x3)+getcsgA_Q2_dep(x1) +4.*getcsgA_Q2_dep(x2));
0836     }
0837     // return g_int2*g_int3/36.;
0838     // return g_int2/6.;
0839     return g_int / 6.;
0840 }
0841
0842 //______________________________________________________________________________
0843 double
0844 photonNucleusCrossSection::integrated_x_section(double const Egamma, double const _min, double const _max)
0845 {
0846     // Integration over full limits gives more accurate result
0847     double Q2_min = std::pow(starlightConstants::mel * Egamma, 2.0) / (_electronEnergy * (_electronEnergy - Egamma));
0848     // double Q2_max = 2.*Egamma/_targetRadii - _wMax*_wMax;
0849     double Q2_max = 4. * _electronEnergy * (_electronEnergy - Egamma);
0850     // Fixed ranges for plot
0851     if (_min != 0 || _max != 0)
0852     {
0853         if (_min > Q2_min)
0854             Q2_min = _min;
0855         if (_max < Q2_max)
0856             Q2_max = _max;
0857     }
0858     // Simpsons rule in using exponential step size and 10000 steps. Tested trapeve rule, linear steps and
0859     // nstep = 10,000 100,000 & 10,000,0000
0860     int nstep = 1000;
0861     double ln_min = std::log(Q2_min);
0862     double ratio = std::log(Q2_max / Q2_min) / nstep;
0863     double g_int = 0;
0864     for (int ii = 0; ii < nstep; ++ii)
0865     {
0866         double x1 = std::exp(ln_min + (double)ii * ratio);
0867         double x3 = std::exp(ln_min + (double)(ii + 1) * ratio);
0868         double x2 = (x3 + x1) / 2.;
0869         // Tests from HERA https://arxiv.org/pdf/hep-ex/9601009.pdf
0870         //     g_int += (x3-x1)*( getcsgA_Q2_dep(x3)+getcsgA_Q2_dep(x1) +4.*getcsgA_Q2_dep(x2));
0871         g_int += (x3 - x1) * (getcsgA_Q2_dep(x3) + getcsgA_Q2_dep(x1) + 4. * getcsgA_Q2_dep(x2));
0872     }
0873     // return g_int2*g_int3/36.;
0874     return g_int / 6.;
0875 }
0876
0877 //______________________________________________________________________________
0878 pair<double, double> *photonNucleusCrossSection::Q2arraylimits(double const Egamma)
0879 {
0880     // double Q2max = 2.*Egamma/_targetRadii - _wMax*_wMax;
0881     double Q2max = 4. * _electronEnergy * (_electronEnergy - Egamma);
0882     double Q2min = std::pow(starlightConstants::mel * Egamma, 2.0) / (_electronEnergy * (_electronEnergy - Egamma));
0883
0884     if (_fixedQ2range == true)
0885     {
0886         if (Q2min < _minQ2)
0887             Q2min = _minQ2;
0888         if (Q2max > _maxQ2)
0889             Q2max = _maxQ2;
0890         // cout<<" Imposed limits at E_gamma = "<<Egamma<<" Q2min: "<<Q2min<<" - Q2max: "<<Q2max<<endl;
0891         std::pair<double, double> *to_ret = new std::pair<double, double>(Q2min, Q2max);
0892         return to_ret;
0893     }
0894     int Nstep = 1000;
0895     //
0896     double ratio = std::log(Q2max / Q2min) / (double)Nstep;
0897     double ln_min = std::log(Q2min);
0898     // -- - -
0899     const double limit = 1E-9;
0900     std::vector<double> Q2_array;
0901     int iNstep = 0;
0902     double g_step = 1.;
0903     while (g_step > limit)
0904     {
0905         double Q2 = std::exp(ln_min + iNstep * ratio);
0906         if (Q2 > Q2max)
0907             break;
0908         g_step = g(Egamma, Q2);
0909         Q2_array.push_back(g_step);
0910         iNstep++;
0911     }
0912     if (std::exp(ln_min + iNstep * ratio) < Q2max)
0913         Q2max = std::exp(ln_min + iNstep * ratio);
0914     // cout<<Q2max<<" "<<g(Egamma,Q2max)*1E9<<endl;
0915     std::pair<double, double> *to_ret = new std::pair<double, double>(Q2min, Q2max);
0916
0917     return to_ret;
0918 }
0919
0920 //______________________________________________________________________________
0921 double
0922 photonNucleusCrossSection::g(double const Egamma,
0923                              double const Q2)
0924 {
0925     // return photonFlux(Egamma,Q2)*getcsgA_Q2_dep(Q2);
0926     return photonFlux(Egamma, Q2); // This could be done more elegantly in the future, but this one change should account for everything
0927 }
0928
0929 //______________________________________________________________________________
0930 double
0931 photonNucleusCrossSection::photonFlux(double const Egamma,
0932                                       double const Q2)
0933 {
0934     // Need to check units later
0935     // double const hbar = starlightConstants::hbarc / 2.99*pow(10,14); // 6.582 * pow (10,-16) [eVs]
0936     // double omega = Egamma/ hbar;
0937     // Do we even need a lookup table for this case? This should return N(E,Q2) from dn = N(E,Q2) dE dQ2
0938     double const ratio = Egamma / _electronEnergy;
0939     double const minQ2 = std::pow(starlightConstants::mel * Egamma, 2.0) / (_electronEnergy * (_electronEnergy - Egamma));
0940     double to_ret = alpha / (pi) * (1 - ratio + ratio * ratio / 2. - (1 - ratio) * (fabs(minQ2 / Q2)));
0941     // Comparisons:
0942     //   double temp = pow(2.*_electronEnergy-Egamma,2.)/(Egamma*Egamma + Q2) + 1. - 4.*starlightConstants::mel*starlightConstants::mel/Q2;
0943     // temp = alpha*temp*Egamma/(4.*Q2*pi*_electronEnergy*_electronEnergy);
0944     //   cout<<" *** Lomnitz *** Testing photon flux approximation for electron energy "<<_electronEnergy<<" gamma "<<Egamma<<" Q2 "<<Q2*1E6<<" 1E-6 "<<endl;
0945     // cout<<" Full expression "<<temp*1E6<<" vs. apporioximation "<<to_ret/( Egamma*fabs(Q2) )*1E6<<" --- ratio "<<temp/(to_ret/( Egamma*fabs(Q2) ))<<endl;
0946     return to_ret / (Egamma * fabs(Q2));
0947     // return temp;
0948 }
0949
0950 //______________________________________________________________________________
0951 double
0952 photonNucleusCrossSection::nepoint(const double Egamma,
0953                                    const double bmin)
0954 {
0955     // Function for the spectrum of virtual photons,
0956     // dn/dEgamma, for a point charge q=Ze sweeping
0957     // past the origin with velocity gamma
0958     // (=1/SQRT(1-(V/c)**2)) integrated over impact
0959     // parameter from bmin to infinity
0960     // See Jackson eq15.54 Classical Electrodynamics
0961     // Declare Local Variables
0962     double beta, X, C1, bracket, nepoint_r;
0963
0964     beta = sqrt(1. - (1. / (_beamLorentzGamma * _beamLorentzGamma)));
0965     X = (bmin * Egamma) / (beta * _beamLorentzGamma * hbarc);
0966
0967     bracket = -0.5 * beta * beta * X * X * (bessel::dbesk1(X) * bessel::dbesk1(X) - bessel::dbesk0(X) * bessel::dbesk0(X));
0968
0969     bracket = bracket + X * bessel::dbesk0(X) * bessel::dbesk1(X);
0970
0971     // Note: NO  Z*Z!!
0972     C1 = (2. * alpha) / pi;
0973
0974     nepoint_r = C1 * (1. / beta) * (1. / beta) * (1. / Egamma) * bracket;
0975
0976     return nepoint_r;
0977 }
0978
0979 //______________________________________________________________________________
0980 double
0981 photonNucleusCrossSection::sigmagp(const double Wgp)
0982 {
0983     // Function for the gamma-proton --> VectorMeson
0984     // cross section. Wgp is the gamma-proton CM energy.
0985     // Unit for cross section: fm**2
0986     // If W_gp is not in the allowed region, i.e. W_GP_MIN < W_gp < W_GP_MAX, do not sample.
0987     if (Wgp < _minW_GP || Wgp > _maxW_GP)
0988         return 0;
0989     double sigmagp_r = 0.;
0990
0991     // Near the threshold CM energy (between WgpMin and WgpMax),
0992     // we add a linear scaling factor to bring the Xsec down to 0
0993     // at the threshold CM value define by WgpMin = m_p + m_vm
0994     double WgpMax = 0.;
0995     double WgpMin = 0.;
0996     double thresholdScaling = 1.0;
0997
0998     switch (_particleType)
0999     {
1000     case RHO:
1001     case RHOZEUS:
1002         WgpMax = 1.8;
1003         WgpMin = 1.60; // this is the cutoff threshold for rho production. But rho has large width so it's lower
1004         if (Wgp < WgpMax)
1005             thresholdScaling = (Wgp - WgpMin) / (WgpMax - WgpMin);
1006         sigmagp_r = thresholdScaling * 1.E-4 * (5.0 * exp(0.22 * log(Wgp)) + 26.0 * exp(-1.23 * log(Wgp)));
1007         // This is based on the omega cross section:
1008         // https://arxiv.org/pdf/2107.06748.pdf
1009         if (_backwardsProduction)
1010             sigmagp_r = thresholdScaling * 1.E-4 * 0.14 * pow(Wgp, -2.7);
1011         if (_backwardsProduction)
1012             sigmagp_r = thresholdScaling * 1.E-4 * 0.206 * pow(((Wgp * Wgp - protonMass * protonMass) / (2.0 * protonMass)), -2.7);
1013         break;
1014     case FOURPRONG:
1015         sigmagp_r = (1. / _rho0PrimeBrFourProng) * 1.E-4 * (0.26 * exp(0.28 * log(Wgp)) + 20.78 * exp(-1.21 * log(Wgp)));
1016         break;
1017     case OMEGA:
1018     case OMEGA_pi0gamma:
1019         WgpMax = 1.8;
1020         WgpMin = 1.74; // this is the cutoff threshold for omega production: W > m_p+m_omega = 1.74
1021         if (Wgp < WgpMax)
1022             thresholdScaling = (Wgp - WgpMin) / (WgpMax - WgpMin);
1023         sigmagp_r = thresholdScaling * 1.E-4 * (0.55 * exp(0.22 * log(Wgp)) + 18.0 * exp(-1.92 * log(Wgp)));
1024         // if(_backwardsProduction)sigmagp_r=thresholdScaling*1.E-4*0.14*pow(Wgp,-2.7);//https://arxiv.org/pdf/2107.06748.pdf
1025         if (_backwardsProduction)
1026             sigmagp_r = thresholdScaling * 1.E-4 * 0.206 * pow(((Wgp * Wgp - protonMass * protonMass) / (2.0 * protonMass)), -2.7);
1027         break;
1028     case PHI:
1029         sigmagp_r = 1.E-4 * 0.34 * exp(0.22 * log(Wgp));
1030         break;
1031     case JPSI:
1032     case JPSI_ee:
1033     case JPSI_mumu:
1034         sigmagp_r = (1.0 - ((_channelMass + protonMass) * (_channelMass + protonMass)) / (Wgp * Wgp));
1035         sigmagp_r *= sigmagp_r;
1036         sigmagp_r *= 1.E-4 * 0.00406 * exp(0.65 * log(Wgp));
1037         // sigmagp_r=1.E-4*0.0015*exp(0.80*log(Wgp));
1038         break;
1039     case JPSI2S:
1040     case JPSI2S_ee:
1041     case JPSI2S_mumu:
1042         sigmagp_r = (1.0 - ((_channelMass + protonMass) * (_channelMass + protonMass)) / (Wgp * Wgp));
1043         sigmagp_r *= sigmagp_r;
1044         sigmagp_r *= 1.E-4 * 0.00406 * exp(0.65 * log(Wgp));
1045         sigmagp_r *= 0.166;
1046         //      sigmagp_r=0.166*(1.E-4*0.0015*exp(0.80*log(Wgp)));
1047         break;
1048     case UPSILON:
1049     case UPSILON_ee:
1050     case UPSILON_mumu:
1051         //       >> This is W**1.7 dependence from QCD calculations
1052         //  sigmagp_r=1.E-10*(0.060)*exp(1.70*log(Wgp));
1053         sigmagp_r = (1.0 - ((_channelMass + protonMass) * (_channelMass + protonMass)) / (Wgp * Wgp));
1054         sigmagp_r *= sigmagp_r;
1055         sigmagp_r *= 1.E-10 * 6.4 * exp(0.74 * log(Wgp));
1056         break;
1057     case UPSILON2S:
1058     case UPSILON2S_ee:
1059     case UPSILON2S_mumu:
1060         // sigmagp_r=1.E-10*(0.0259)*exp(1.70*log(Wgp));
1061         sigmagp_r = (1.0 - ((_channelMass + protonMass) * (_channelMass + protonMass)) / (Wgp * Wgp));
1062         sigmagp_r *= sigmagp_r;
1063         sigmagp_r *= 1.E-10 * 2.9 * exp(0.74 * log(Wgp));
1064         break;
1065     case UPSILON3S:
1066     case UPSILON3S_ee:
1067     case UPSILON3S_mumu:
1068         // sigmagp_r=1.E-10*(0.0181)*exp(1.70*log(Wgp));
1069         sigmagp_r = (1.0 - ((_channelMass + protonMass) * (_channelMass + protonMass)) / (Wgp * Wgp));
1070         sigmagp_r *= sigmagp_r;
1071         sigmagp_r *= 1.E-10 * 2.1 * exp(0.74 * log(Wgp));
1072         break;
1073     default:
1074         cout << "!!!  ERROR: Unidentified Vector Meson: " << _particleType << endl;
1075     }
1076     return sigmagp_r;
1077 }
1078
1079 //______________________________________________________________________________
1080 double
1081 photonNucleusCrossSection::sigma_A(const double sig_N, const int beam)
1082 {
1083     // Nuclear Cross Section
1084     // sig_N,sigma_A in (fm**2)
1085
1086     double sum;
1087     double b, bmax, Pint, arg, sigma_A_r;
1088
1089     int NGAUSS;
1090
1091     double xg[17] =
1092         {.0,
1093          .0483076656877383162, .144471961582796493,
1094          .239287362252137075, .331868602282127650,
1095          .421351276130635345, .506899908932229390,
1096          .587715757240762329, .663044266930215201,
1097          .732182118740289680, .794483795967942407,
1098          .849367613732569970, .896321155766052124,
1099          .934906075937739689, .964762255587506430,
1100          .985611511545268335, .997263861849481564};
1101
1102     double ag[17] =
1103         {.0,
1104          .0965400885147278006, .0956387200792748594,
1105          .0938443990808045654, .0911738786957638847,
1106          .0876520930044038111, .0833119242269467552,
1107          .0781938957870703065, .0723457941088485062,
1108          .0658222227763618468, .0586840934785355471,
1109          .0509980592623761762, .0428358980222266807,
1110          .0342738629130214331, .0253920653092620595,
1111          .0162743947309056706, .00701861000947009660};
1112
1113     NGAUSS = 16;
1114
1115     // Check if one or both beams are nuclei
1116     int A_1 = _bbs.electronBeam().A();
1117     int A_2 = _bbs.targetBeam().A();
1118     if (A_1 == 1 && A_2 == 1)
1119         cout << " This is pp, you should not be here..." << endl;
1120
1121     // CALCULATE P(int) FOR b=0.0 - bmax (fm)
1122     bmax = 25.0;
1123     sum = 0.;
1124     for (int IB = 1; IB <= NGAUSS; IB++)
1125     {
1126
1127         b = 0.5 * bmax * xg[IB] + 0.5 * bmax;
1128
1129         if (A_1 == 1 && A_2 != 1)
1130         {
1131             arg = -sig_N * _bbs.targetBeam().rho0() * _bbs.targetBeam().thickness(b);
1132         }
1133         else if (A_2 == 1 && A_1 != 1)
1134         {
1135             arg = -sig_N * _bbs.electronBeam().rho0() * _bbs.electronBeam().thickness(b);
1136         }
1137         else
1138         {
1139             // Check which beam is target
1140             if (beam == 1)
1141             {
1142                 arg = -sig_N * _bbs.electronBeam().rho0() * _bbs.electronBeam().thickness(b);
1143             }
1144             else if (beam == 2)
1145             {
1146                 arg = -sig_N * _bbs.targetBeam().rho0() * _bbs.targetBeam().thickness(b);
1147             }
1148             else
1149             {
1150                 cout << " Something went wrong here, beam= " << beam << endl;
1151             }
1152         }
1153
1154         Pint = 1.0 - exp(arg);
1155         // If this is a quantum Glauber calculation, use the quantum Glauber formula
1156         if (_quantumGlauber == 1)
1157         {
1158             Pint = 2.0 * (1.0 - exp(arg / 2.0));
1159         }
1160
1161         sum = sum + 2. * pi * b * Pint * ag[IB];
1162
1163         b = 0.5 * bmax * (-xg[IB]) + 0.5 * bmax;
1164
1165         if (A_1 == 1 && A_2 != 1)
1166         {
1167             arg = -sig_N * _bbs.targetBeam().rho0() * _bbs.targetBeam().thickness(b);
1168         }
1169         else if (A_2 == 1 && A_1 != 1)
1170         {
1171             arg = -sig_N * _bbs.electronBeam().rho0() * _bbs.electronBeam().thickness(b);
1172         }
1173         else
1174         {
1175             // Check which beam is target
1176             if (beam == 1)
1177             {
1178                 arg = -sig_N * _bbs.electronBeam().rho0() * _bbs.electronBeam().thickness(b);
1179             }
1180             else if (beam == 2)
1181             {
1182                 arg = -sig_N * _bbs.targetBeam().rho0() * _bbs.targetBeam().thickness(b);
1183             }
1184             else
1185             {
1186                 cout << " Something went wrong here, beam= " << beam << endl;
1187             }
1188         }
1189
1190         Pint = 1.0 - exp(arg);
1191         // If this is a quantum Glauber calculation, use the quantum Glauber formula
1192         if (_quantumGlauber == 1)
1193         {
1194             Pint = 2.0 * (1.0 - exp(arg / 2.0));
1195         }
1196
1197         sum = sum + 2. * pi * b * Pint * ag[IB];
1198     }
1199
1200     sum = 0.5 * bmax * sum;
1201
1202     sigma_A_r = sum;
1203
1204     return sigma_A_r;
1205 }
1206
1207 //______________________________________________________________________________
1208 double
1209 photonNucleusCrossSection::sigma_N(const double Wgp)
1210 {
1211     // Nucleon Cross Section in (fm**2)
1212     double cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha);
1213     return cs;
1214 }
1215
1216 //______________________________________________________________________________
1217 double
1218 photonNucleusCrossSection::breitWigner(const double W,
1219                                        const double C)
1220 {
1221     // use simple fixed-width s-wave Breit-Wigner without coherent backgorund for rho'
1222     // (PDG '08 eq. 38.56)
1223     if (_particleType == FOURPRONG)
1224     {
1225         if (W < 4.01 * pionChargedMass)
1226             return 0;
1227         const double termA = _channelMass * _width;
1228         const double termA2 = termA * termA;
1229         const double termB = W * W - _channelMass * _channelMass;
1230         return C * _ANORM * _ANORM * termA2 / (termB * termB + termA2);
1231     }
1232
1233     // Relativistic Breit-Wigner according to J.D. Jackson,
1234     // Nuovo Cimento 34, 6692 (1964), with nonresonant term. A is the strength
1235     // of the resonant term and b the strength of the non-resonant
1236     // term. C is an overall normalization.
1237
1238     double ppi = 0., ppi0 = 0., GammaPrim, rat;
1239     double aa, bb, cc;
1240
1241     double nrbw_r;
1242
1243     // width depends on energy - Jackson Eq. A.2
1244     // if below threshold, then return 0.  Added 5/3/2001 SRK
1245     // 0.5% extra added for safety margin
1246     // omega added here 10/26/2014 SRK
1247     if (_particleType == RHO || _particleType == RHOZEUS || _particleType == OMEGA)
1248     {
1249         if (W < 2.01 * pionChargedMass)
1250         {
1251             nrbw_r = 0.;
1252             return nrbw_r;
1253         }
1254         ppi = sqrt(((W / 2.) * (W / 2.)) - pionChargedMass * pionChargedMass);
1255         ppi0 = 0.358;
1256     }
1257     if (_particleType == OMEGA_pi0gamma)
1258     {
1259         if (W < pionNeutralMass)
1260         {
1261             nrbw_r = 0.;
1262             return nrbw_r;
1263         }
1264         ppi = abs((W * W - pionNeutralMass * pionNeutralMass) / (2.0 * W));
1265         ppi0 = abs((_channelMass * _channelMass - pionNeutralMass * pionNeutralMass) / (2.0 * W));
1266     }
1267
1268     // handle phi-->K+K- properly
1269     if (_particleType == PHI)
1270     {
1271         if (W < 2. * kaonChargedMass)
1272         {
1273             nrbw_r = 0.;
1274             return nrbw_r;
1275         }
1276         ppi = sqrt(((W / 2.) * (W / 2.)) - kaonChargedMass * kaonChargedMass);
1277         ppi0 = sqrt(((_channelMass / 2.) * (_channelMass / 2.)) - kaonChargedMass * kaonChargedMass);
1278     }
1279
1280     // handle J/Psi-->e+e- properly
1281     if (_particleType == JPSI || _particleType == JPSI2S)
1282     {
1283         if (W < 2. * mel)
1284         {
1285             nrbw_r = 0.;
1286             return nrbw_r;
1287         }
1288         ppi = sqrt(((W / 2.) * (W / 2.)) - mel * mel);
1289         ppi0 = sqrt(((_channelMass / 2.) * (_channelMass / 2.)) - mel * mel);
1290     }
1291     if (_particleType == JPSI_ee)
1292     {
1293         if (W < 2. * mel)
1294         {
1295             nrbw_r = 0.;
1296             return nrbw_r;
1297         }
1298         ppi = sqrt(((W / 2.) * (W / 2.)) - mel * mel);
1299         ppi0 = sqrt(((_channelMass / 2.) * (_channelMass / 2.)) - mel * mel);
1300     }
1301     if (_particleType == JPSI_mumu)
1302     {
1303         if (W < 2. * muonMass)
1304         {
1305             nrbw_r = 0.;
1306             return nrbw_r;
1307         }
1308         ppi = sqrt(((W / 2.) * (W / 2.)) - muonMass * muonMass);
1309         ppi0 = sqrt(((_channelMass / 2.) * (_channelMass / 2.)) - muonMass * muonMass);
1310     }
1311     if (_particleType == JPSI2S_ee)
1312     {
1313         if (W < 2. * mel)
1314         {
1315             nrbw_r = 0.;
1316             return nrbw_r;
1317         }
1318         ppi = sqrt(((W / 2.) * (W / 2.)) - mel * mel);
1319         ppi0 = sqrt(((_channelMass / 2.) * (_channelMass / 2.)) - mel * mel);
1320     }
1321     if (_particleType == JPSI2S_mumu)
1322     {
1323         if (W < 2. * muonMass)
1324         {
1325             nrbw_r = 0.;
1326             return nrbw_r;
1327         }
1328         ppi = sqrt(((W / 2.) * (W / 2.)) - muonMass * muonMass);
1329         ppi0 = sqrt(((_channelMass / 2.) * (_channelMass / 2.)) - muonMass * muonMass);
1330     }
1331
1332     if (_particleType == UPSILON || _particleType == UPSILON2S || _particleType == UPSILON3S)
1333     {
1334         if (W < 2. * muonMass)
1335         {
1336             nrbw_r = 0.;
1337             return nrbw_r;
1338         }
1339         ppi = sqrt(((W / 2.) * (W / 2.)) - muonMass * muonMass);
1340         ppi0 = sqrt(((_channelMass / 2.) * (_channelMass / 2.)) - muonMass * muonMass);
1341     }
1342
1343     if (_particleType == UPSILON_mumu || _particleType == UPSILON2S_mumu || _particleType == UPSILON3S_mumu)
1344     {
1345         if (W < 2. * muonMass)
1346         {
1347             nrbw_r = 0.;
1348             return nrbw_r;
1349         }
1350         ppi = sqrt(((W / 2.) * (W / 2.)) - muonMass * muonMass);
1351         ppi0 = sqrt(((_channelMass / 2.) * (_channelMass / 2.)) - muonMass * muonMass);
1352     }
1353
1354     if (_particleType == UPSILON_ee || _particleType == UPSILON2S_ee || _particleType == UPSILON3S_ee)
1355     {
1356         if (W < 2. * mel)
1357         {
1358             nrbw_r = 0.;
1359             return nrbw_r;
1360         }
1361         ppi = sqrt(((W / 2.) * (W / 2.)) - mel * mel);
1362         ppi0 = sqrt(((_channelMass / 2.) * (_channelMass / 2.)) - mel * mel);
1363     }
1364
1365     if (ppi == 0. && ppi0 == 0.)
1366         cout << "Improper Gammaacrosssection::breitwigner, ppi&ppi0=0." << endl;
1367
1368     rat = ppi / ppi0;
1369     GammaPrim = _width * (_channelMass / W) * rat * rat * rat;
1370
1371     aa = _ANORM * sqrt(GammaPrim * _channelMass * W);
1372     bb = W * W - _channelMass * _channelMass;
1373     cc = _channelMass * GammaPrim;
1374
1375     // First real part squared
1376     nrbw_r = (((aa * bb) / (bb * bb + cc * cc) + _BNORM) * ((aa * bb) / (bb * bb + cc * cc) + _BNORM));
1377
1378     // Then imaginary part squared
1379     nrbw_r = nrbw_r + (((aa * cc) / (bb * bb + cc * cc)) * ((aa * cc) / (bb * bb + cc * cc)));
1380
1381     //  Alternative, a simple, no-background BW, following J. Breitweg et al.
1382     //  Eq. 15 of Eur. Phys. J. C2, 247 (1998).  SRK 11/10/2000
1383     //      nrbw_r = (_ANORM*_mass*GammaPrim/(bb*bb+cc*cc))**2
1384
1385     nrbw_r = C * nrbw_r;
1386
1387     return nrbw_r;
1388 }