Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2024-09-27 07:03:41 

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 
0033 #include <iostream>
0034 #include <fstream>
0035 #include <cmath>
0036 
0037 #include "reportingUtils.h"
0038 #include "starlightconstants.h"
0039 #include "bessel.h"
0040 #include "photonNucleusCrossSection.h"
0041 
0042 
0043 using namespace std;
0044 using namespace starlightConstants;
0045 
0046 
0047 //______________________________________________________________________________
0048 photonNucleusCrossSection::photonNucleusCrossSection(const inputParameters& inputParametersInstance, const beamBeamSystem& bbsystem)
0049     : _nWbins            (inputParametersInstance.nmbWBins()          ),
0050       _nYbins            (inputParametersInstance.nmbRapidityBins()   ),
0051       _wMin              (inputParametersInstance.minW()              ),
0052       _wMax              (inputParametersInstance.maxW()              ),
0053       _yMax              (inputParametersInstance.maxRapidity()       ),
0054       _beamLorentzGamma  (inputParametersInstance.beamLorentzGamma()  ),
0055       _bbs               (bbsystem                                    ),
0056       _protonEnergy      (inputParametersInstance.protonEnergy()      ),
0057       _electronEnergy    (inputParametersInstance.electronEnergy()    ),
0058       _particleType      (inputParametersInstance.prodParticleType()  ),
0059       _beamBreakupMode   (inputParametersInstance.beamBreakupMode()   ),
0060       _backwardsProduction(inputParametersInstance.backwardsProduction()),
0061       _productionMode    (inputParametersInstance.productionMode()    ),
0062       _sigmaNucleus      (_bbs.targetBeam().A()          ),
0063       _fixedQ2range      (inputParametersInstance.fixedQ2Range()      ),
0064       _minQ2             (inputParametersInstance.minGammaQ2()        ),
0065       _maxQ2             (inputParametersInstance.maxGammaQ2()        ),
0066       _maxPhotonEnergy   (inputParametersInstance.cmsMaxPhotonEnergy()),
0067       _cmsMinPhotonEnergy(inputParametersInstance.cmsMinPhotonEnergy()),
0068       _targetRadii       (inputParametersInstance.targetRadius()      ),
0069       _maxW_GP           (inputParametersInstance.maxW_GP()              ),
0070       _minW_GP           (inputParametersInstance.minW_GP()              )
0071 {
0072         // new options - impulse aproximation (per Joakim) and Quantum Glauber (per SK) SKQG
0073         _impulseSelected = inputParametersInstance.impulseVM();
0074     _quantumGlauber = inputParametersInstance.quantumGlauber();
0075     switch(_particleType) {
0076     case RHO:
0077       _slopeParameter = 11.0;  // [(GeV/c)^{-2}]
0078       _vmPhotonCoupling = 2.02;
0079       _vmQ2Power_c1 = 2.09;
0080       _vmQ2Power_c2 = 0.0073; // [ GeV^{-2}]
0081       _ANORM       = -2.75;
0082       _BNORM       = 0.0;
0083       _defaultC    = 1.0;
0084       _channelMass = starlightConstants::rho0Mass; 
0085       _width       = starlightConstants::rho0Width; 
0086       if(_backwardsProduction){
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      = 11.0;
0104       _vmPhotonCoupling      = 2.02;
0105       _vmQ2Power_c1 = 2.09;
0106       _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
0107       _ANORM       = -2.75;
0108       _BNORM       = 0;  
0109       _defaultC    = 11.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         _slopeParameter = 21.8;  // [(GeV/c)^{-2}]
0126         _ANORM          = 1.;
0127       }
0128       break;
0129     case PHI:
0130       _slopeParameter=7.0;
0131       _vmPhotonCoupling=13.71;
0132       _vmQ2Power_c1 = 2.15;
0133       _vmQ2Power_c2 = 0.0074; // [GeV^{-2}]
0134       _ANORM=-2.75;
0135       _BNORM=0.0;
0136       _defaultC=1.0;
0137       _channelMass  = starlightConstants::PhiMass;
0138       _width        = starlightConstants::PhiWidth;
0139       break;
0140     case JPSI:
0141     case JPSI_ee:
0142       _slopeParameter=4.0;
0143       _vmPhotonCoupling=10.45;
0144       _vmQ2Power_c1 = 2.45;
0145       _vmQ2Power_c2 = 0.00084; // [GeV^{-2}]
0146       _ANORM=-2.75; 
0147       _BNORM=0.0;
0148       _defaultC=1.0;
0149       _channelMass  = starlightConstants::JpsiMass; 
0150       _width        = starlightConstants::JpsiWidth; 
0151       break;
0152     case JPSI_mumu:
0153       _slopeParameter=4.0;
0154       _vmPhotonCoupling=10.45;
0155       _vmQ2Power_c1 = 2.36;
0156       _vmQ2Power_c2 = 0.0029; // [GeV^{-2}]
0157       _ANORM=-2.75; 
0158       _BNORM=0.0;
0159       _defaultC=1.0;
0160       _channelMass  = starlightConstants::JpsiMass; 
0161       _width        = starlightConstants::JpsiWidth; 
0162       break;
0163     case JPSI2S:
0164     case JPSI2S_ee:
0165     case JPSI2S_mumu:
0166       _slopeParameter=4.3;
0167       _vmPhotonCoupling=26.39;
0168       _vmQ2Power_c1 = 2.09;
0169       _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
0170       _ANORM=-2.75; 
0171       _BNORM=0.0;
0172       _defaultC=1.0;
0173       _channelMass  = starlightConstants::Psi2SMass;
0174       _width        = starlightConstants::Psi2SWidth;
0175       break;
0176     case UPSILON:
0177     case UPSILON_ee:
0178     case UPSILON_mumu:
0179       _slopeParameter=4.0;
0180       _vmPhotonCoupling=125.37;
0181       _vmQ2Power_c1 = 2.09;
0182       _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
0183       _ANORM=-2.75; 
0184       _BNORM=0.0;
0185       _defaultC=1.0;
0186       _channelMass  = starlightConstants::Upsilon1SMass;
0187       _width        = starlightConstants::Upsilon1SWidth;
0188       break;
0189     case UPSILON2S:
0190     case UPSILON2S_ee:
0191     case UPSILON2S_mumu:
0192       _slopeParameter=4.0;
0193       _vmPhotonCoupling=290.84;
0194       _vmQ2Power_c1 = 2.09;
0195       _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
0196       _ANORM=-2.75;
0197       _BNORM=0.0;
0198       _defaultC=1.0;
0199       _channelMass  = starlightConstants::Upsilon2SMass;
0200       _width        = starlightConstants::Upsilon2SWidth;       
0201       break;
0202     case UPSILON3S:
0203     case UPSILON3S_ee:
0204     case UPSILON3S_mumu:
0205       _slopeParameter=4.0;
0206       _vmPhotonCoupling=415.10;
0207       _vmQ2Power_c1 = 2.09;
0208       _vmQ2Power_c2 = 0.0073; // [GeV^{-2}]
0209       _ANORM=-2.75;
0210       _BNORM=0.0;
0211       _defaultC=1.0;
0212       _channelMass  = starlightConstants::Upsilon3SMass;
0213       _width        = starlightConstants::Upsilon3SWidth;
0214       break;
0215     default:
0216         cout <<"No sigma constants parameterized for pid: "<<_particleType
0217              <<" GammaAcrosssection"<<endl;
0218     }
0219 
0220     //Changed by Lomnitz for e case. Limit is now E_e - 100m_e
0221     //_maxPhotonEnergy = 12. * _beamLorentzGamma * hbarc/(_bbs.beam1().nuclearRadius()+_bbs.targetBeam().nuclearRadius());
0222     //_maxPhotonEnergy = _electronEnergy - 10.*starlightConstants::mel;
0223     /*cout<<" Lomnitz:: max energy in target frame "<< _electronEnergy - 1000.*starlightConstants::mel<<" vs electron energy "<<_electronEnergy<<endl
0224         <<"           max energy in cms frame    "<<_maxPhotonEnergy<<"  vs electron energy "<<_beamLorentzGamma*starlightConstants::mel<<endl;
0225         cout<<" testing original limit "<< 12. * _beamLorentzGamma * hbarc/(2.*_bbs.targetBeam().nuclearRadius())<<endl;*/
0226     
0227       
0228 }
0229 
0230 
0231 //______________________________________________________________________________
0232 photonNucleusCrossSection::~photonNucleusCrossSection()
0233 { }
0234 
0235 
0236 //______________________________________________________________________________
0237 void
0238 photonNucleusCrossSection::crossSectionCalculation(const double)
0239 {
0240     cout << "Neither narrow/wide resonance cross-section calculation.--Derived" << endl;
0241 }
0242 
0243 //______________________________________________________________________________
0244 double
0245 photonNucleusCrossSection::getcsgA(const double targetEgamma,
0246                                    const double Q2, 
0247                                    const int beam)
0248 {
0249     //This function returns the cross-section for photon-nucleus interaction 
0250     //producing vectormesons
0251   
0252     double Av,Wgp,cs,cvma;
0253     double t,tmin,tmax;
0254     double csgA,ax,bx;
0255     int NGAUSS; 
0256     //
0257     double W = _channelMass; //new change, channel mass center used for the t min integration.
0258       
0259     //     DATA FOR GAUSS INTEGRATION
0260     double xg[6] = {0, 0.1488743390, 0.4333953941, 0.6794095683, 0.8650633667, 0.9739065285};
0261     double ag[6] = {0, 0.2955242247, 0.2692667193, 0.2190863625, 0.1494513492, 0.0666713443};
0262     NGAUSS = 6;
0263     
0264     //       Note: The photon energy passed to this function is now in the target frame. The rest of the calculations are done in the
0265     //       CMS frame. The next lines boost the photon into the CM frame.
0266     double E_prime = _electronEnergy - targetEgamma;
0267     double cos_theta_e = 1. - Q2/(2.*_electronEnergy*E_prime);
0268     double theta_e = acos(cos_theta_e);
0269     double beam_y = acosh(_beamLorentzGamma);   
0270     double gamma_pt = E_prime*sin(theta_e);
0271     double pz_squared = targetEgamma*targetEgamma - Q2 - gamma_pt*gamma_pt;
0272     if( pz_squared < 0 || fabs(cos_theta_e) > 1 || 2.*targetEgamma/(Q2+W*W) < _targetRadii)
0273       return 0;
0274     double temp_pz = sqrt(pz_squared);
0275     // Now boost to CM frame
0276     double Egamma = targetEgamma*cosh(beam_y) - temp_pz*sinh(beam_y);
0277     if( Egamma < _cmsMinPhotonEnergy || Egamma > _maxPhotonEnergy){
0278       return 0;
0279     }
0280     //cout<<" ::: Lomnitz test in photonNucleus ::: pz^2 = "<<pz_squared << " CMS Egamma = "<<Egamma<<endl;
0281     //       Find gamma-proton CM energy in CMS frame in the limit Q2->0 (this is our assumption, the Q2 dependence is in the factor)
0282     Wgp = sqrt( 2.*(protonMass*targetEgamma)+protonMass*protonMass);
0283     /*Wgp = sqrt(2. * Egamma * (_protonEnergy
0284                               + sqrt(_protonEnergy * _protonEnergy - protonMass * protonMass))
0285                   + protonMass * protonMass);*/
0286     
0287     //Used for A-A
0288     tmin = (W * W / (4. * Egamma * _beamLorentzGamma)) * (W * W / (4. * Egamma * _beamLorentzGamma));
0289   
0290     if ((_bbs.electronBeam().A() <= 1) && (_bbs.targetBeam().A() <= 1)){
0291        // proton-proton, no scaling needed
0292       csgA = getcsgA_Q2_dep(Q2)*sigmagp(Wgp);
0293     } else {
0294        // Check if one or both beams are nuclei 
0295        int A_1 = _bbs.electronBeam().A(); 
0296        int A_2 = _bbs.targetBeam().A(); 
0297        // coherent AA interactions
0298        // Calculate V.M.+proton cross section
0299            // cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha); 
0300            cs = getcsgA_Q2_dep(Q2)*sigma_N(Wgp); //Use member function instead 
0301     
0302        // Calculate V.M.+nucleus cross section
0303        cvma = sigma_A(cs,beam); 
0304 
0305            // Do impulse approximation here
0306            if( _impulseSelected == 1){
0307              if( beam == 1 ){
0308            cvma = A_1*cs;
0309          } else if ( beam == 2 ){
0310                cvma = A_2*cs;
0311          }   
0312            }       
0313 
0314        // Calculate Av = dsigma/dt(t=0) Note Units: fm**s/Gev**2
0315        Av = (alpha * cvma * cvma) / (16. * pi * _vmPhotonCoupling * hbarc * hbarc);
0316    
0317        tmax   = tmin + 0.25;
0318        ax     = 0.5 * (tmax - tmin);
0319        bx     = 0.5 * (tmax + tmin);
0320        csgA   = 0.;
0321        for (int k = 1; k < NGAUSS; ++k) { 
0322 
0323            t    = ax * xg[k] + bx;
0324                if( A_1 <= 1 && A_2 != 1){ 
0325           csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
0326                }else if(A_2 <=1 && A_1 != 1){
0327           csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0328                }else{     
0329                   if( beam==1 ){
0330              csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0331                   }else if(beam==2){
0332              csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);   
0333           }else{
0334              cout<<"Something went wrong here, beam= "<<beam<<endl; 
0335                   }
0336                }
0337 
0338            t    = ax * (-xg[k]) + bx;
0339                if( A_1 <= 1 && A_2 != 1){ 
0340               csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
0341                }else if(A_2 <=1 && A_1 != 1){
0342               csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0343                }else{     
0344                   if( beam==1 ){
0345                 csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0346                   }else if(beam==2){
0347                 csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);    
0348           }else{
0349                 cout<<"Something went wrong here, beam= "<<beam<<endl; 
0350                   }
0351            }
0352        }
0353        csgA = 0.5 * (tmax - tmin) * csgA;
0354        csgA = Av * csgA;
0355     }
0356     return csgA;    
0357 }
0358 
0359 
0360 //______________________________________________________________________________
0361 double
0362 photonNucleusCrossSection::e_getcsgA(const double Egamma, double Q2,
0363                                    const double W, 
0364                                    const int beam)
0365 {
0366     //This function returns the cross-section for photon-nucleus interaction 
0367     //producing vectormesons for e_starlight. Stuff will be returned in CMS frame, but
0368         //photon input is take in target frame
0369   
0370     double Av,Wgp,cs,cvma;
0371     double t,tmin,tmax;
0372     double csgA,ax,bx;
0373     int NGAUSS; 
0374   
0375     //     DATA FOR GAUSS INTEGRATION
0376     double xg[6] = {0, 0.1488743390, 0.4333953941, 0.6794095683, 0.8650633667, 0.9739065285};
0377     double ag[6] = {0, 0.2955242247, 0.2692667193, 0.2190863625, 0.1494513492, 0.0666713443};
0378     NGAUSS = 6;
0379   
0380     //       Find gamma-proton CM energy
0381     Wgp = sqrt( 2.*(protonMass*Egamma)
0382             +protonMass*protonMass + Q2);
0383     
0384     //Used for A-A
0385     tmin = (W * W / (4. * Egamma * _beamLorentzGamma)) * (W * W / (4. * Egamma * _beamLorentzGamma));
0386   
0387     if ((_bbs.electronBeam().A() <= 1) && (_bbs.targetBeam().A() <= 1)){
0388        // proton-proton, no scaling needed
0389        csgA = sigmagp(Wgp);
0390     } else {
0391        // coherent AA interactions
0392        // Calculate V.M.+proton cross section
0393            // cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha); 
0394            cs = sigma_N(Wgp); //Use member function instead 
0395     
0396        // Calculate V.M.+nucleus cross section
0397        cvma = sigma_A(cs,beam); 
0398 
0399        // Calculate Av = dsigma/dt(t=0) Note Units: fm**s/Gev**2
0400        Av = (alpha * cvma * cvma) / (16. * pi * _vmPhotonCoupling * hbarc * hbarc);
0401 
0402            // Check if one or both beams are nuclei 
0403            int A_1 = _bbs.electronBeam().A(); 
0404            int A_2 = _bbs.targetBeam().A(); 
0405    
0406        tmax   = tmin + 0.25;
0407        ax     = 0.5 * (tmax - tmin);
0408        bx     = 0.5 * (tmax + tmin);
0409        csgA   = 0.;
0410        for (int k = 1; k < NGAUSS; ++k) { 
0411 
0412            t    = ax * xg[k] + bx;
0413                if( A_1 <= 1 && A_2 != 1){ 
0414           csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
0415                }else if(A_2 <=1 && A_1 != 1){
0416           csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0417                }else{     
0418                   if( beam==1 ){
0419              csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0420                   }else if(beam==2){
0421              csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);   
0422           }else{
0423              cout<<"Something went wrong here, beam= "<<beam<<endl; 
0424                   }
0425                }
0426 
0427            t    = ax * (-xg[k]) + bx;
0428                if( A_1 <= 1 && A_2 != 1){ 
0429               csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);
0430                }else if(A_2 <=1 && A_1 != 1){
0431               csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0432                }else{     
0433                   if( beam==1 ){
0434                 csgA = csgA + ag[k] * _bbs.electronBeam().formFactor(t) * _bbs.electronBeam().formFactor(t);
0435                   }else if(beam==2){
0436                 csgA = csgA + ag[k] * _bbs.targetBeam().formFactor(t) * _bbs.targetBeam().formFactor(t);    
0437           }else{
0438                 cout<<"Something went wrong here, beam= "<<beam<<endl; 
0439                   }
0440            }
0441        }
0442        csgA = 0.5 * (tmax - tmin) * csgA;
0443        csgA = Av * csgA;
0444     }
0445     return csgA;    
0446 }
0447 
0448 
0449 //______________________________________________________________________________
0450 double
0451 photonNucleusCrossSection::getcsgA_Q2_dep(const double Q2)
0452 {
0453   double const mv2 = getChannelMass()*getChannelMass();
0454   double const n = vmQ2Power(Q2);
0455   return std::pow(mv2/(mv2+Q2),n);
0456 }
0457 
0458 
0459 //______________________________________________________________________________
0460 double
0461 photonNucleusCrossSection::photonFlux(const double Egamma, const int beam)
0462 {
0463     // This routine gives the photon flux as a function of energy Egamma
0464     // It works for arbitrary nuclei and gamma; the first time it is
0465     // called, it calculates a lookup table which is used on
0466     // subsequent calls.
0467     // It returns dN_gamma/dE (dimensions 1/E), not dI/dE
0468     // energies are in GeV, in the lab frame
0469     // rewritten 4/25/2001 by SRK
0470 
0471         // NOTE: beam (=1,2) defines the photon TARGET
0472 
0473     double lEgamma,Emin,Emax;
0474     static double lnEmax, lnEmin, dlnE;
0475     double stepmult,energy,rZ;
0476     int nbstep,nrstep,nphistep,nstep;
0477     double bmin,bmax,bmult,biter,bold,integratedflux;
0478     double fluxelement,deltar,riter;
0479     double deltaphi,phiiter,dist;
0480     static double dide[401];
0481     double lnElt;
0482     double flux_r; 
0483     double Xvar;
0484     int Ilt;
0485     double RNuc=0.,RSum=0.;
0486 
0487     RSum=_bbs.electronBeam().nuclearRadius()+_bbs.targetBeam().nuclearRadius();
0488         if( beam == 1){
0489           rZ=double(_bbs.targetBeam().Z());
0490           RNuc = _bbs.electronBeam().nuclearRadius();
0491         } else { 
0492       rZ=double(_bbs.electronBeam().Z());
0493           RNuc = _bbs.targetBeam().nuclearRadius();
0494         }
0495 
0496     static int  Icheck = 0;
0497     static int  Ibeam  = 0; 
0498   
0499     //Check first to see if pp 
0500     if( _bbs.electronBeam().A()==1 && _bbs.targetBeam().A()==1 ){
0501         int nbsteps = 400;
0502         double bmin = 0.5;
0503         double bmax = 5.0 + (5.0*_beamLorentzGamma*hbarc/Egamma);
0504         double dlnb = (log(bmax)-log(bmin))/(1.*nbsteps);
0505  
0506         double local_sum=0.0;
0507 
0508         // Impact parameter loop 
0509         for (int i = 0; i<=nbsteps;i++){
0510 
0511             double bnn0 = bmin*exp(i*dlnb);
0512             double bnn1 = bmin*exp((i+1)*dlnb);
0513             double db   = bnn1-bnn0;
0514         
0515             double ppslope = 19.0; 
0516             double GammaProfile = exp(-bnn0*bnn0/(2.*hbarc*hbarc*ppslope));  
0517             double PofB0 = 1. - (1. - GammaProfile)*(1. - GammaProfile);   
0518             GammaProfile = exp(-bnn1*bnn1/(2.*hbarc*hbarc*ppslope));  
0519             double PofB1 = 1. - (1. - GammaProfile)*(1. - GammaProfile);   
0520 
0521                         double loc_nofe0 = _bbs.electronBeam().photonDensity(bnn0,Egamma);
0522             double loc_nofe1 = _bbs.targetBeam().photonDensity(bnn1,Egamma);
0523 
0524             local_sum += 0.5*loc_nofe0*(1. - PofB0)*2.*starlightConstants::pi*bnn0*db; 
0525             local_sum += 0.5*loc_nofe1*(1. - PofB1)*2.*starlightConstants::pi*bnn1*db; 
0526 
0527         }
0528         // End Impact parameter loop 
0529         return local_sum;
0530     }
0531 
0532     //   first call or new beam?  - initialize - calculate photon flux
0533     Icheck=Icheck+1;
0534 
0535     // Do the numerical integration only once for symmetric systems. 
0536         if( Icheck > 1 && _bbs.electronBeam().A() == _bbs.targetBeam().A() && _bbs.electronBeam().Z() == _bbs.targetBeam().Z() ) goto L1000f;
0537         // For asymmetric systems check if we have another beam 
0538     if( Icheck > 1 && beam == Ibeam ) goto L1000f; 
0539         Ibeam = beam; 
0540   
0541     //  Nuclear breakup is done by PofB
0542     //  collect number of integration steps here, in one place
0543   
0544     nbstep=1200;
0545     nrstep=60;
0546     nphistep=40;
0547   
0548     //  this last one is the number of energy steps
0549     nstep=100;
0550  
0551     // following previous choices, take Emin=10 keV at LHC, Emin = 1 MeV at RHIC
0552     Emin=1.E-5;
0553     if (_beamLorentzGamma < 500) 
0554         Emin=1.E-3;
0555   
0556         Emax=_maxPhotonEnergy; 
0557     // Emax=12.*hbarc*_beamLorentzGamma/RSum;
0558  
0559     //     >> lnEmin <-> ln(Egamma) for the 0th bin
0560     //     >> lnEmax <-> ln(Egamma) for the last bin
0561     lnEmin=log(Emin);
0562     lnEmax=log(Emax);
0563     dlnE=(lnEmax-lnEmin)/nstep; 
0564 
0565         printf("Calculating photon flux from Emin = %e GeV to Emax = %e GeV (CM frame) for source with Z = %3.0f \n", Emin, Emax, rZ);
0566     
0567     stepmult= exp(log(Emax/Emin)/double(nstep));
0568     energy=Emin;
0569   
0570     for (int j = 1; j<=nstep;j++){
0571         energy=energy*stepmult;
0572     
0573         //  integrate flux over 2R_A < b < 2R_A+ 6* gamma hbar/energy
0574         //  use exponential steps
0575     
0576         bmin=0.8*RSum; //Start slightly below 2*Radius 
0577         bmax=bmin + 6.*hbarc*_beamLorentzGamma/energy;
0578     
0579         bmult=exp(log(bmax/bmin)/double(nbstep));
0580         biter=bmin;
0581         integratedflux=0.;
0582 
0583         if( (_bbs.electronBeam().A() == 1 && _bbs.targetBeam().A() != 1) || (_bbs.targetBeam().A() == 1 && _bbs.electronBeam().A() != 1) ){
0584             // This is pA 
0585 
0586           if( _productionMode == PHOTONPOMERONINCOHERENT ){
0587               // This pA incoherent, proton is the target
0588 
0589               int nbsteps = 400;
0590               double bmin = 0.7*RSum;
0591               double bmax = 2.0*RSum + (8.0*_beamLorentzGamma*hbarc/energy);
0592               double dlnb = (log(bmax)-log(bmin))/(1.*nbsteps);
0593 
0594               double local_sum=0.0;
0595 
0596               // Impact parameter loop 
0597               for (int i = 0; i<=nbsteps; i++){
0598 
0599                   double bnn0 = bmin*exp(i*dlnb);
0600               double bnn1 = bmin*exp((i+1)*dlnb);
0601               double db   = bnn1-bnn0;
0602 
0603                           double PofB0 = _bbs.probabilityOfBreakup(bnn0); 
0604                           double PofB1 = _bbs.probabilityOfBreakup(bnn1); 
0605       
0606               double loc_nofe0 = 0.0;
0607               double loc_nofe1 = 0.0; 
0608                           if( _bbs.electronBeam().A() == 1 ){
0609                 loc_nofe0 = _bbs.targetBeam().photonDensity(bnn0,energy);
0610                 loc_nofe1 = _bbs.targetBeam().photonDensity(bnn1,energy);
0611                           }
0612                   else if( _bbs.targetBeam().A() == 1 ){
0613                 loc_nofe0 = _bbs.electronBeam().photonDensity(bnn0,energy);
0614                 loc_nofe1 = _bbs.electronBeam().photonDensity(bnn1,energy);             
0615                           }
0616 
0617                           // cout<<" i: "<<i<<" bnn0: "<<bnn0<<" PofB0: "<<PofB0<<" loc_nofe0: "<<loc_nofe0<<endl; 
0618 
0619               local_sum += 0.5*loc_nofe0*PofB0*2.*starlightConstants::pi*bnn0*db; 
0620               local_sum += 0.5*loc_nofe1*PofB1*2.*starlightConstants::pi*bnn1*db;  
0621               }  // End Impact parameter loop 
0622                   integratedflux = local_sum; 
0623                     } else if ( _productionMode == PHOTONPOMERONNARROW ||  _productionMode == PHOTONPOMERONWIDE ){
0624                       // This is pA coherent, nucleus is the target 
0625                       double localbmin = 0.0;   
0626                       if( _bbs.electronBeam().A() == 1 ){
0627             localbmin = _bbs.targetBeam().nuclearRadius() + 0.7; 
0628               }
0629                       if( _bbs.targetBeam().A() == 1 ){ 
0630             localbmin = _bbs.electronBeam().nuclearRadius() + 0.7; 
0631                       }
0632                       integratedflux = nepoint(energy,localbmin); 
0633             }
0634         }else{ 
0635         // This is AA
0636         for (int jb = 1; jb<=nbstep;jb++){
0637             bold=biter;
0638             biter=biter*bmult;
0639             // When we get to b>20R_A change methods - just take the photon flux
0640             //  at the center of the nucleus.
0641             if (biter > (10.*RNuc)){
0642                // if there is no nuclear breakup or only hadronic breakup, which only
0643                // occurs at smaller b, we can analytically integrate the flux from b~20R_A
0644                // to infinity, following Jackson (2nd edition), Eq. 15.54
0645                Xvar=energy*biter/(hbarc*_beamLorentzGamma);
0646                // Here, there is nuclear breakup.  So, we can't use the integrated flux
0647                // However, we can do a single flux calculation, at the center of the nucleus
0648                // Eq. 41 of Vidovic, Greiner and Soff, Phys.Rev.C47,2308(1993), among other places
0649                // this is the flux per unit area
0650                fluxelement  = (rZ*rZ*alpha*energy)*
0651                        (bessel::dbesk1(Xvar))*(bessel::dbesk1(Xvar))/
0652                        ((pi*_beamLorentzGamma*hbarc)*
0653                        (pi*_beamLorentzGamma*hbarc));
0654         
0655                    } else {
0656                // integrate over nuclear surface. n.b. this assumes total shadowing -
0657                // treat photons hitting the nucleus the same no matter where they strike
0658                fluxelement=0.;
0659                deltar=RNuc/double(nrstep);
0660                riter=-deltar/2.;
0661           
0662                for (int jr =1; jr<=nrstep;jr++){
0663                riter=riter+deltar;
0664                // use symmetry;  only integrate from 0 to pi (half circle)
0665                deltaphi=pi/double(nphistep);
0666                phiiter=0.;
0667             
0668                for( int jphi=1;jphi<= nphistep;jphi++){
0669                    phiiter=(double(jphi)-0.5)*deltaphi;
0670                    // dist is the distance from the center of the emitting nucleus 
0671                    // to the point in question
0672                    dist=sqrt((biter+riter*cos(phiiter))*(biter+riter*
0673                      cos(phiiter))+(riter*sin(phiiter))*(riter*sin(phiiter)));
0674                    Xvar=energy*dist/(hbarc*_beamLorentzGamma);  
0675                    flux_r = (rZ*rZ*alpha*energy)*
0676                      (bessel::dbesk1(Xvar)*bessel::dbesk1(Xvar))/
0677                      ((pi*_beamLorentzGamma*hbarc)*
0678                      (pi*_beamLorentzGamma*hbarc));
0679           
0680                      //  The surface  element is 2.* delta phi* r * delta r
0681                      //  The '2' is because the phi integral only goes from 0 to pi
0682                      fluxelement=fluxelement+flux_r*2.*deltaphi*riter*deltar;
0683                      //  end phi and r integrations
0684                }//for(jphi)
0685                }//for(jr)
0686               //  average fluxelement over the nuclear surface
0687               fluxelement=fluxelement/(pi*RNuc*RNuc);
0688            }//else
0689               //  multiply by volume element to get total flux in the volume element
0690               fluxelement=fluxelement*2.*pi*biter*(biter-bold);
0691               //  modulate by the probability of nuclear breakup as f(biter)
0692                           // cout<<" jb: "<<jb<<" biter: "<<biter<<" fluxelement: "<<fluxelement<<endl; 
0693               if (_beamBreakupMode > 1){
0694                   fluxelement=fluxelement*_bbs.probabilityOfBreakup(biter);
0695               }
0696                           // cout<<" jb: "<<jb<<" biter: "<<biter<<" fluxelement: "<<fluxelement<<endl; 
0697               integratedflux=integratedflux+fluxelement;
0698       
0699         } //end loop over impact parameter 
0700         }  //end of else (pp, pA, AA) 
0701     
0702         //  In lookup table, store k*dN/dk because it changes less
0703         //  so the interpolation should be better    
0704         dide[j]=integratedflux*energy;                                     
0705     }//end loop over photon energy 
0706        
0707     //  for 2nd and subsequent calls, use lookup table immediately
0708   
0709  L1000f:
0710   
0711     lEgamma=log(Egamma);
0712     if (lEgamma < (lnEmin+dlnE) ||  lEgamma  > lnEmax){
0713         flux_r=0.0;
0714         cout<<"  ERROR: Egamma outside defined range. Egamma= "<<Egamma
0715             <<"   "<<lnEmax<<" "<<(lnEmin+dlnE)<<endl;
0716     }
0717     else{
0718         //       >> Egamma between Ilt and Ilt+1
0719         Ilt = int((lEgamma-lnEmin)/dlnE);
0720         //       >> ln(Egamma) for first point 
0721         lnElt = lnEmin + Ilt*dlnE; 
0722         //       >> Interpolate
0723         flux_r = dide[Ilt] + ((lEgamma-lnElt)/dlnE)*(dide[Ilt+1]- dide[Ilt]);
0724         flux_r = flux_r/Egamma;
0725     }
0726   
0727     return flux_r;
0728 }
0729 
0730 
0731 //______________________________________________________________________________
0732 double 
0733 photonNucleusCrossSection::integrated_Q2_dep(double const Egamma, double const _min, double const _max)
0734 {
0735   //Integration over full limits gives more accurate result
0736   double Q2_min =  std::pow(starlightConstants::mel*Egamma,2.0)/(_electronEnergy*(_electronEnergy-Egamma));
0737   double Q2_max = 4.*_electronEnergy*(_electronEnergy-Egamma);
0738   //double Q2_max = 2.*Egamma/_targetRadii - _wMax*_wMax;
0739   if( _min != 0 || _max !=0){
0740     if( _min > Q2_min )
0741       Q2_min = _min;
0742     if( _max < Q2_max )
0743       Q2_max = _max;
0744   }
0745   // Simpsons rule in using exponential step size and 10000 steps. Tested trapeve rule, linear steps and
0746   // nstep = 10,000 100,000 & 10,000,0000
0747   int nstep = 1000;
0748   double ln_min = std::log(Q2_min);
0749   double ratio = std::log(Q2_max/Q2_min)/nstep;
0750   double g_int = 0;
0751   //double g_int2 = 0 ;
0752   //double g_int3 = 0;
0753   //cout<<"*** Lomnitz **** Energy "<<Egamma<<" limits "<<Q2_min*1E9<<" x 1E-9 -  "<<Q2_max<<endl;
0754   for ( int ii = 0 ; ii< nstep; ++ii){
0755     double x1 =  std::exp(ln_min+(double)ii*ratio);
0756     double x3 =  std::exp(ln_min+(double)(ii+1)*ratio);
0757     double x2 =  (x3+x1)/2.;
0758     //cout<<"ii : "<<x1<<" "<<x2<<" "<<x3<<endl;
0759     g_int += (x3-x1)*( g(Egamma,x3)+g(Egamma,x1) +4.*g(Egamma,x2));
0760     //g_int2 += (x3-x1)*( photonFlux(Egamma,x3)+photonFlux(Egamma,x1) +4.*photonFlux(Egamma,x2));
0761     //g_int3 += (x3-x1)*( getcsgA_Q2_dep(x3)+getcsgA_Q2_dep(x1) +4.*getcsgA_Q2_dep(x2));
0762   }
0763   //return g_int2*g_int3/36.; 
0764   //return g_int2/6.;
0765   return g_int/6.;
0766 }
0767 
0768 
0769 //______________________________________________________________________________
0770 double 
0771 photonNucleusCrossSection::integrated_x_section(double const Egamma, double const _min, double const _max)
0772 {
0773   //Integration over full limits gives more accurate result
0774   double Q2_min =  std::pow(starlightConstants::mel*Egamma,2.0)/(_electronEnergy*(_electronEnergy-Egamma));
0775   //double Q2_max = 2.*Egamma/_targetRadii - _wMax*_wMax;
0776   double Q2_max  = 4.*_electronEnergy*(_electronEnergy-Egamma);
0777   // Fixed ranges for plot
0778   if( _min != 0 || _max!=0){
0779     if( _min > Q2_min)
0780       Q2_min = _min;
0781     if( _max < Q2_max)
0782       Q2_max = _max;
0783   }
0784   // Simpsons rule in using exponential step size and 10000 steps. Tested trapeve rule, linear steps and
0785   // nstep = 10,000 100,000 & 10,000,0000
0786   int nstep = 1000;
0787   double ln_min = std::log(Q2_min);
0788   double ratio = std::log(Q2_max/Q2_min)/nstep;
0789   double g_int = 0;
0790   for ( int ii = 0 ; ii< nstep; ++ii){
0791     double x1 =  std::exp(ln_min+(double)ii*ratio);
0792     double x3 =  std::exp(ln_min+(double)(ii+1)*ratio);
0793     double x2 =  (x3+x1)/2.;
0794     //Tests from HERA https://arxiv.org/pdf/hep-ex/9601009.pdf
0795     //    g_int += (x3-x1)*( getcsgA_Q2_dep(x3)+getcsgA_Q2_dep(x1) +4.*getcsgA_Q2_dep(x2));
0796     g_int += (x3-x1)*( getcsgA_Q2_dep(x3)+getcsgA_Q2_dep(x1) +4.*getcsgA_Q2_dep(x2));
0797   }
0798   //return g_int2*g_int3/36.; 
0799   return g_int/6.;
0800 }
0801 
0802 
0803 //______________________________________________________________________________
0804 pair< double, double >*photonNucleusCrossSection::Q2arraylimits(double const Egamma)
0805 {
0806   //double Q2max = 2.*Egamma/_targetRadii - _wMax*_wMax;
0807   double Q2max = 4.*_electronEnergy*(_electronEnergy-Egamma);
0808   double Q2min= std::pow(starlightConstants::mel*Egamma,2.0)/(_electronEnergy*(_electronEnergy-Egamma));
0809 
0810   if( _fixedQ2range == true){
0811     if( Q2min < _minQ2 )
0812       Q2min = _minQ2;
0813     if( Q2max > _maxQ2 )
0814       Q2max = _maxQ2;
0815     //cout<<" Imposed limits "<<Q2min<<" - "<<Q2max<<endl;
0816     std::pair<double,double>* to_ret = new std::pair<double, double>(Q2min,Q2max);
0817     return to_ret;
0818   }
0819   int Nstep = 1000;
0820   //
0821   double ratio = std::log(Q2max/Q2min)/(double)Nstep;
0822   double ln_min = std::log(Q2min);
0823   // -- - -
0824   const double limit = 1E-9;
0825   std::vector<double>Q2_array;
0826   int iNstep = 0;
0827   double g_step = 1.;
0828   while( g_step>limit ){
0829     double Q2 = std::exp(ln_min+iNstep*ratio);
0830     if(Q2>Q2max) 
0831       break;
0832     g_step = g(Egamma,Q2);
0833     Q2_array.push_back(g_step);
0834     iNstep++;
0835   }
0836   if( std::exp(ln_min+iNstep*ratio) < Q2max)
0837     Q2max = std::exp(ln_min+iNstep*ratio);
0838   //cout<<Q2max<<" "<<g(Egamma,Q2max)*1E9<<endl;
0839   std::pair<double, double>* to_ret = new std::pair<double, double>(Q2min,Q2max);
0840 
0841   return to_ret;
0842 }
0843 
0844 
0845 //______________________________________________________________________________
0846 double 
0847 photonNucleusCrossSection::g(double const Egamma,
0848                  double const Q2)
0849 {
0850   //return photonFlux(Egamma,Q2)*getcsgA_Q2_dep(Q2);
0851   return photonFlux(Egamma,Q2); //This could be done more elegantly in the future, but this one change should account for everything
0852 }
0853 
0854 //______________________________________________________________________________
0855 double 
0856 photonNucleusCrossSection::photonFlux(double const Egamma, 
0857                       double const Q2)
0858 {
0859   //Need to check units later
0860   //double const hbar = starlightConstants::hbarc / 2.99*pow(10,14); // 6.582 * pow (10,-16) [eVs]
0861   //double omega = Egamma/ hbar;
0862   //Do we even need a lookup table for this case? This should return N(E,Q2) from dn = N(E,Q2) dE dQ2
0863   double const ratio = Egamma/_electronEnergy;
0864   double const minQ2 = std::pow( starlightConstants::mel*Egamma,2.0) / (_electronEnergy*(_electronEnergy - Egamma));
0865   double to_ret = alpha/(pi) *( 1- ratio + ratio*ratio/2. - (1-ratio)*( fabs(minQ2/Q2)) );
0866   //Comparisons:
0867   //  double temp = pow(2.*_electronEnergy-Egamma,2.)/(Egamma*Egamma + Q2) + 1. - 4.*starlightConstants::mel*starlightConstants::mel/Q2;
0868   //temp = alpha*temp*Egamma/(4.*Q2*pi*_electronEnergy*_electronEnergy);
0869   //  cout<<" *** Lomnitz *** Testing photon flux approximation for electron energy "<<_electronEnergy<<" gamma "<<Egamma<<" Q2 "<<Q2*1E6<<" 1E-6 "<<endl;
0870   //cout<<" Full expression "<<temp*1E6<<" vs. apporioximation "<<to_ret/( Egamma*fabs(Q2) )*1E6<<" --- ratio "<<temp/(to_ret/( Egamma*fabs(Q2) ))<<endl;
0871   return to_ret/( Egamma*fabs(Q2) );
0872   //return temp;
0873 }
0874 
0875 //______________________________________________________________________________
0876 double
0877 photonNucleusCrossSection::nepoint(const double Egamma,
0878                                    const double bmin)
0879 {
0880     // Function for the spectrum of virtual photons,
0881     // dn/dEgamma, for a point charge q=Ze sweeping
0882     // past the origin with velocity gamma
0883     // (=1/SQRT(1-(V/c)**2)) integrated over impact
0884     // parameter from bmin to infinity
0885     // See Jackson eq15.54 Classical Electrodynamics
0886     // Declare Local Variables
0887     double beta,X,C1,bracket,nepoint_r;
0888   
0889     beta = sqrt(1.-(1./(_beamLorentzGamma*_beamLorentzGamma)));
0890     X = (bmin*Egamma)/(beta*_beamLorentzGamma*hbarc);
0891   
0892     bracket = -0.5*beta*beta*X*X*(bessel::dbesk1(X)*bessel::dbesk1(X)
0893                                   -bessel::dbesk0(X)*bessel::dbesk0(X));
0894 
0895     bracket = bracket+X*bessel::dbesk0(X)*bessel::dbesk1(X);
0896   
0897     // Note: NO  Z*Z!!
0898     C1=(2.*alpha)/pi;
0899   
0900     nepoint_r = C1*(1./beta)*(1./beta)*(1./Egamma)*bracket;
0901   
0902     return nepoint_r;
0903   
0904 }
0905 
0906 
0907 //______________________________________________________________________________
0908 double
0909 photonNucleusCrossSection::sigmagp(const double Wgp)
0910 {
0911     // Function for the gamma-proton --> VectorMeson
0912     // cross section. Wgp is the gamma-proton CM energy.
0913     // Unit for cross section: fm**2
0914     //If W_gp is not in the allowed region, i.e. W_GP_MIN < W_gp < W_GP_MAX, do not sample.
0915     if(Wgp < _minW_GP || Wgp > _maxW_GP) return 0;
0916     double sigmagp_r=0.;
0917 
0918     // Near the threshold CM energy (between WgpMin and WgpMax), 
0919     // we add a linear scaling factor to bring the Xsec down to 0 
0920     // at the threshold CM value define by WgpMin = m_p + m_vm
0921     double WgpMax = 0.;
0922     double WgpMin = 0.;
0923     double thresholdScaling = 1.0;
0924   
0925     switch(_particleType)
0926         { 
0927         case RHO:
0928         case RHOZEUS:
0929         case FOURPRONG:
0930             WgpMax = 1.8;
0931             WgpMin = 1.60; //this is the cutoff threshold for rho production. But rho has large width so it's lower
0932             if(Wgp<WgpMax) thresholdScaling=(Wgp-WgpMin)/(WgpMax-WgpMin);
0933             sigmagp_r=thresholdScaling*1.E-4*(5.0*exp(0.22*log(Wgp))+26.0*exp(-1.23*log(Wgp)));
0934             //This is based on the omega cross section:
0935             //https://arxiv.org/pdf/2107.06748.pdf
0936             if(_backwardsProduction)sigmagp_r=thresholdScaling*1.E-4*0.14*pow(Wgp,-2.7);
0937             if(_backwardsProduction)sigmagp_r=thresholdScaling*1.E-4*0.206*pow(((Wgp*Wgp-protonMass*protonMass)/(2.0*protonMass)),-2.7);
0938             break;
0939         case OMEGA:
0940         case OMEGA_pi0gamma:
0941             WgpMax = 1.8;
0942             WgpMin = 1.74; //this is the cutoff threshold for omega production: W > m_p+m_omega = 1.74
0943             if(Wgp<WgpMax) thresholdScaling=(Wgp-WgpMin)/(WgpMax-WgpMin);
0944             sigmagp_r=thresholdScaling*1.E-4*(0.55*exp(0.22*log(Wgp))+18.0*exp(-1.92*log(Wgp)));
0945             //if(_backwardsProduction)sigmagp_r=thresholdScaling*1.E-4*0.14*pow(Wgp,-2.7);//https://arxiv.org/pdf/2107.06748.pdf
0946             if(_backwardsProduction)sigmagp_r=thresholdScaling*1.E-4*0.206*pow(((Wgp*Wgp-protonMass*protonMass)/(2.0*protonMass)),-2.7);
0947             break;                                                      
0948         case PHI:
0949             sigmagp_r=1.E-4*0.34*exp(0.22*log(Wgp));
0950             break;
0951         case JPSI:
0952         case JPSI_ee:
0953         case JPSI_mumu:
0954             sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
0955             sigmagp_r*=sigmagp_r;
0956             sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp));
0957             // sigmagp_r=1.E-4*0.0015*exp(0.80*log(Wgp));
0958             break;
0959         case JPSI2S:
0960         case JPSI2S_ee:
0961         case JPSI2S_mumu:
0962             sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
0963             sigmagp_r*=sigmagp_r;
0964             sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp));
0965             sigmagp_r*=0.166;  
0966             //      sigmagp_r=0.166*(1.E-4*0.0015*exp(0.80*log(Wgp)));
0967             break;
0968         case UPSILON:
0969         case UPSILON_ee:
0970         case UPSILON_mumu:
0971             //       >> This is W**1.7 dependence from QCD calculations
0972             //  sigmagp_r=1.E-10*(0.060)*exp(1.70*log(Wgp));
0973             sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
0974             sigmagp_r*=sigmagp_r;
0975             sigmagp_r*=1.E-10*6.4*exp(0.74*log(Wgp));
0976             break;
0977         case UPSILON2S:
0978         case UPSILON2S_ee:
0979         case UPSILON2S_mumu:
0980                 // sigmagp_r=1.E-10*(0.0259)*exp(1.70*log(Wgp));
0981                 sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
0982             sigmagp_r*=sigmagp_r;
0983             sigmagp_r*=1.E-10*2.9*exp(0.74*log(Wgp)); 
0984             break;
0985         case UPSILON3S:
0986         case UPSILON3S_ee:
0987         case UPSILON3S_mumu:
0988                 // sigmagp_r=1.E-10*(0.0181)*exp(1.70*log(Wgp));
0989                 sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
0990             sigmagp_r*=sigmagp_r;
0991             sigmagp_r*=1.E-10*2.1*exp(0.74*log(Wgp)); 
0992             break;
0993         default: cout<< "!!!  ERROR: Unidentified Vector Meson: "<< _particleType <<endl;
0994         }                                                                  
0995     return sigmagp_r;
0996 }
0997 
0998 
0999 //______________________________________________________________________________
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!