STARLIGHT/starlight/src/.svn/text-base/photonNucleusCrossSection.cpp.svn-base

   1 ///////////////////////////////////////////////////////////////////////////
   2 //
   3 //    Copyright 2010
   4 //
   5 //    This file is part of starlight.
   6 //
   7 //    starlight is free software: you can redistribute it and/or modify
   8 //    it under the terms of the GNU General Public License as published by
   9 //    the Free Software Foundation, either version 3 of the License, or
  10 //    (at your option) any later version.
  11 //
  12 //    starlight is distributed in the hope that it will be useful,
  13 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 //    GNU General Public License for more details.
  16 //
  17 //    You should have received a copy of the GNU General Public License
  18 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
  19 //
  20 ///////////////////////////////////////////////////////////////////////////
  21 //
  22 // File and Version Information:
  23 // $Rev::                             $: revision of last commit
  24 // $Author::                          $: author of last commit
  25 // $Date::                            $: date of last commit
  26 //
  27 // Description:
  28 //
  29 //
  30 ///////////////////////////////////////////////////////////////////////////
  31
  32
  33 #include <iostream>
  34 #include <fstream>
  35 #include <cmath>
  36
  37 #include "reportingUtils.h"
  38 #include "starlightconstants.h"
  39 #include "bessel.h"
  40 #include "photonNucleusCrossSection.h"
  41
  42
  43 using namespace std;
  44 using namespace starlightConstants;
  45
  46
  47 //______________________________________________________________________________
  48 photonNucleusCrossSection::photonNucleusCrossSection(const beamBeamSystem& bbsystem)
  49         : _nWbins            (inputParametersInstance.nmbWBins()          ),
  50           _nYbins            (inputParametersInstance.nmbRapidityBins()   ),
  51           _wMin              (inputParametersInstance.minW()              ),
  52           _wMax              (inputParametersInstance.maxW()              ),
  53           _yMax              (inputParametersInstance.maxRapidity()       ),
  54           _beamLorentzGamma  (inputParametersInstance.beamLorentzGamma()  ),
  55           _bbs               (bbsystem                                    ),
  56           _protonEnergy      (inputParametersInstance.protonEnergy()      ),
  57           _particleType      (inputParametersInstance.prodParticleType()  ),
  58           _beamBreakupMode   (inputParametersInstance.beamBreakupMode()   ),
  59           _coherentProduction(inputParametersInstance.coherentProduction()),
  60           _incoherentFactor  (inputParametersInstance.incoherentFactor()  ),
  61           _productionMode    (inputParametersInstance.productionMode()    ),
  62           _sigmaNucleus      (_bbs.beam2().A()          )
  63 {
  64         // define luminosity for various beam particles in units of 10^{26} cm^{-2} sec^{-1}
  65         switch(_bbs.beam1().Z()) {
  66         case 1:   // proton
  67                 _luminosity = 1.E8;
  68                 break;
  69         case 8:   // O
  70                 _luminosity = 980.;
  71                 break;
  72         case 14:  // Si
  73                 _luminosity = 440.;
  74                 break;
  75         case 20:  // Ca
  76                 _luminosity = 2000.;
  77                 break;
  78         case 29:  // Cu
  79                 _luminosity = 95.;
  80                 break;
  81         case 49:  // Indium, uses same as Iodine
  82                 _luminosity = 27.;
  83                 break;
  84         case 53:  // I
  85                 _luminosity = 27.;
  86                 break;
  87         case 79:  // Au
  88                 _luminosity = 2.0;
  89                 break;
  90         case 82:  // Pb
  91                 _luminosity = 1.;
  92                 break;
  93         default:
  94                 printWarn << "luminosity is not defined for beam with Z = " << _bbs.beam1().Z()
  95                           << ". using " << _luminosity << " 10^{26} cm^{-2} sec^{-1}" << endl;
  96         }
  97
  98         switch(_particleType) {
  99         case RHO:
 100                 _slopeParameter = 11.0;  // [(GeV/c)^{-2}]
 101                 _vmPhotonCoupling = 2.02;
 102                 _ANORM       = -2.75;
 103                 _BNORM       = 0.0;
 104                 _defaultC    = 1.0;
 105                 _channelMass = 0.7685;  // [GeV/c^2]
 106                 _width       = 0.1507;  // [GeV/c^2]
 107                 break;
 108         case RHOZEUS:
 109                 _slopeParameter =11.0;
 110                 _vmPhotonCoupling=2.02;
 111                 _ANORM=-2.75;
 112                 _BNORM=1.84;
 113                 _defaultC=1.0;
 114                 _channelMass = 0.7685;
 115                 _width=0.1507;
 116                 break;
 117         case FOURPRONG:
 118                 _slopeParameter      = 11.0;
 119                 _vmPhotonCoupling      = 2.02;
 120                 _ANORM       = -2.75;
 121                 _BNORM       = 0;  // no coherent background component is implemented for four-prong
 122                 _defaultC    = 11.0;
 123                 _channelMass = 1.350;
 124                 _width       = 0.360;
 125                 break;
 126         case OMEGA:
 127                 _slopeParameter=10.0;
 128                 _vmPhotonCoupling=23.13;
 129                 _ANORM=-2.75;
 130                 _BNORM=0.0;
 131                 _defaultC=1.0;
 132                 _channelMass=0.78194;
 133                 _width=0.00843;
 134                 break;
 135         case PHI:
 136                 _slopeParameter=7.0;
 137                 _vmPhotonCoupling=13.71;
 138                 _ANORM=-2.75;
 139                 _BNORM=0.0;
 140                 _defaultC=1.0;
 141                 _channelMass=1.019413;
 142                 _width=0.00443;
 143                 break;
 144         case JPSI:
 145         case JPSI_ee:
 146         case JPSI_mumu:
 147                 _slopeParameter=4.0;
 148                 _vmPhotonCoupling=10.45;
 149                 _ANORM=-2.75;//Artificial Breit-Wigner parameters--no direct pions
 150                 _BNORM=0.0;
 151                 _defaultC=1.0;
 152                 _channelMass=3.09692;//JN 3.09688
 153                 _width=0.000091;//JN 0.000087
 154                 break;
 155         case JPSI2S:
 156         case JPSI2S_ee:
 157         case JPSI2S_mumu:
 158                 _slopeParameter=4.3;
 159                 _vmPhotonCoupling=26.39;
 160                 _ANORM=-2.75;//Artificial
 161                 _BNORM=0.0;
 162                 _defaultC=1.0;
 163                 _channelMass=3.686093;
 164                 _width=0.000337;
 165                 break;
 166         case UPSILON:
 167         case UPSILON_ee:
 168         case UPSILON_mumu:
 169                 _slopeParameter=4.0;
 170                 _vmPhotonCoupling=125.37;
 171                 _ANORM=-2.75;//Artificial
 172                 _BNORM=0.0;
 173                 _defaultC=1.0;
 174                 _channelMass=9.46030;
 175                 _width=0.00005402;
 176                 break;
 177         case UPSILON2S:
 178         case UPSILON2S_ee:
 179         case UPSILON2S_mumu:
 180                 _slopeParameter=4.0;
 181                 _vmPhotonCoupling=290.84;
 182                 _ANORM=-2.75;
 183                 _BNORM=0.0;
 184                 _defaultC=1.0;
 185                 _channelMass=10.02326;
 186                 _width=0.00003198;
 187                 break;
 188         case UPSILON3S:
 189         case UPSILON3S_ee:
 190         case UPSILON3S_mumu:
 191                 _slopeParameter=4.0;
 192                 _vmPhotonCoupling=415.10;
 193                 _ANORM=-2.75;
 194                 _BNORM=0.0;
 195                 _defaultC=1.0;
 196                 _channelMass=10.3552;
 197                 _width=0.00002032;
 198                 break;
 199         default:
 200                 cout <<"No sigma constants parameterized for pid: "<<_particleType
 201                      <<" GammaAcrosssection"<<endl;
 202         }
 203
 204         double maxradius = (_bbs.beam1().nuclearRadius()<_bbs.beam2().nuclearRadius())?_bbs.beam2().nuclearRadius():_bbs.beam1().nuclearRadius();
 205         _maxPhotonEnergy = 5. * _beamLorentzGamma * hbarc/maxradius;
 206 }
 207
 208
 209 //______________________________________________________________________________
 210 photonNucleusCrossSection::~photonNucleusCrossSection()
 211 { }
 212
 213
 214 //______________________________________________________________________________
 215 void
 216 photonNucleusCrossSection::crossSectionCalculation(const double)
 217 {
 218         cout << "Neither narrow/wide resonance cross-section calculation.--Derived" << endl;
 219 }
 220
 221
 222 //______________________________________________________________________________
 223 double
 224 photonNucleusCrossSection::getcsgA(const double Egamma,
 225                                    const double W)
 226 {
 227         //This function returns the cross-section for photon-nucleus interaction
 228         //producing vectormesons
 229
 230         double Av,Wgp,cs,cvma;
 231         double t,tmin,tmax;
 232         double csgA,ax,bx;
 233         int NGAUSS;
 234
 235         //     DATA FOR GAUSS INTEGRATION
 236         double xg[6] = {0, 0.1488743390, 0.4333953941, 0.6794095683, 0.8650633667, 0.9739065285};
 237         double ag[6] = {0, 0.2955242247, 0.2692667193, 0.2190863625, 0.1494513492, 0.0666713443};
 238         NGAUSS = 6;
 239
 240         //       Find gamma-proton CM energy
 241         Wgp = sqrt(2. * Egamma * (_protonEnergy
 242                                   + sqrt(_protonEnergy * _protonEnergy - protonMass * protonMass))
 243                    + protonMass * protonMass);
 244
 245         //Used for d-A and A-A
 246         tmin = (W * W / (4. * Egamma * _beamLorentzGamma)) * (W * W / (4. * Egamma * _beamLorentzGamma));
 247
 248         if ((_bbs.beam1().A() == 1) && (_bbs.beam2().A() == 1))  // proton-proton, no scaling needed
 249                 csgA = sigmagp(Wgp);
 250         else if ((_bbs.beam2().Z() == 1) && (_bbs.beam2().A() == 2)) {  // deuteron-A interaction
 251                 Av = _slopeParameter * sigmagp(Wgp);
 252
 253                 tmax   = tmin + 0.64;   //0.64
 254                 ax     = 0.5 * (tmax - tmin);
 255                 bx     = 0.5 * (tmax + tmin);
 256                 csgA   = 0.;
 257
 258                 for (int k = 1; k < NGAUSS; ++k) {
 259                         t    = ax * xg[k] + bx;
 260                         // We use beam2 here since the input stores the deuteron as nucleus 2
 261                         // and nucleus 2 is the pomeron field source
 262                         // Also this is the way sergey formatted the formfactor.
 263                         csgA = csgA + ag[k] * _bbs.beam2().formFactor(t);
 264                         t    = ax * (-xg[k]) + bx;
 265                         csgA = csgA + ag[k] * _bbs.beam2().formFactor(t);
 266                 }
 267                 csgA = 0.5 * (tmax - tmin) * csgA;
 268                 csgA = Av * csgA;
 269         }       else if (!_coherentProduction &&
 270                          (!((_bbs.beam2().Z() == 1) && (_bbs.beam2().A() == 2)))) {  // incoherent AA interactions
 271                 // For incoherent AA interactions, since incoherent treating it as gamma-p
 272                 // Calculate the differential V.M.+proton cross section
 273                 csgA = 1.E-4 * _incoherentFactor * _sigmaNucleus * _slopeParameter * sigmagp(Wgp);
 274
 275         }       else {  // coherent AA interactions
 276                 // For typical AA interactions.
 277                 // Calculate V.M.+proton cross section
 278                 cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha);
 279
 280                 // Calculate V.M.+nucleus cross section
 281                 // cvma = _bbs.beam1().A()*cs;
 282                 cvma = sigma_A(cs);
 283
 284                 // Calculate Av = dsigma/dt(t=0) Note Units: fm**s/Gev**2
 285                 Av = (alpha * cvma * cvma) / (16. * pi * _vmPhotonCoupling * hbarc * hbarc);
 286
 287                 // Check if one or both beams are nuclei
 288                 int A_1 = _bbs.beam1().A();
 289                 int A_2 = _bbs.beam2().A();
 290
 291                 tmax   = tmin + 0.25;
 292                 ax     = 0.5 * (tmax - tmin);
 293                 bx     = 0.5 * (tmax + tmin);
 294                 csgA   = 0.;
 295                 for (int k = 1; k < NGAUSS; ++k) {
 296
 297                         t    = ax * xg[k] + bx;
 298                         if( A_1 == 1 && A_2 != 1){
 299                           csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
 300                         }else if(A_2 ==1 && A_1 != 1){
 301                           csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
 302                         }else{
 303                           csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
 304                         }
 305
 306                         t    = ax * (-xg[k]) + bx;
 307                         if( A_1 == 1 && A_2 != 1){
 308                           csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
 309                         }else if(A_2 ==1 && A_1 != 1){
 310                           csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
 311                         }else{
 312                           csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
 313                         }
 314                 }
 315                 csgA = 0.5 * (tmax - tmin) * csgA;
 316                 csgA = Av * csgA;
 317         }
 318
 319         return csgA;
 320 }
 321
 322
 323 //______________________________________________________________________________
 324 double
 325 photonNucleusCrossSection::photonFlux(const double Egamma)
 326 {
 327         // This routine gives the photon flux as a function of energy Egamma
 328         // It works for arbitrary nuclei and gamma; the first time it is
 329         // called, it calculates a lookup table which is used on
 330         // subsequent calls.
 331         // It returns dN_gamma/dE (dimensions 1/E), not dI/dE
 332         // energies are in GeV, in the lab frame
 333         // rewritten 4/25/2001 by SRK
 334
 335         double lEgamma,Emin,Emax;
 336         static double lnEmax, lnEmin, dlnE;
 337         double stepmult,energy,rZ,rA;
 338         int nbstep,nrstep,nphistep,nstep;
 339         double bmin,bmax,bmult,biter,bold,integratedflux;
 340         double fluxelement,deltar,riter;
 341         double deltaphi,phiiter,dist;
 342         static double dide[401];
 343         double lnElt;
 344         double rA2, rZ2;
 345         double flux_r;
 346         double Xvar;
 347         int Ilt;
 348         double RNuc=0.,RNuc2=0.,maxradius=0.;
 349
 350         RNuc=_bbs.beam1().nuclearRadius();
 351         RNuc2=_bbs.beam2().nuclearRadius();
 352         maxradius =  (RNuc<RNuc2)?RNuc2:RNuc;
 353         // static ->>> dide,lnEMax,lnEmin,dlnE
 354         static int  Icheck = 0;
 355
 356         //Check first to see if pp
 357         if( _bbs.beam1().A()==1 && _bbs.beam2().A()==1 ){
 358                 int nbsteps = 400;
 359                 double bmin = 0.5;
 360                 double bmax = 5.0 + (5.0*_beamLorentzGamma*hbarc/Egamma);
 361                 double dlnb = (log(bmax)-log(bmin))/(1.*nbsteps);
 362
 363                 double local_sum=0.0;
 364
 365                 // Impact parameter loop
 366                 for (int i = 0; i<=nbsteps;i++){
 367
 368                         double bnn0 = bmin*exp(i*dlnb);
 369                         double bnn1 = bmin*exp((i+1)*dlnb);
 370                         double db   = bnn1-bnn0;
 371
 372                         //      double PofB0 = 1.0;
 373                         //      if( bnn0 > 1.4 )PofB0=0.0;
 374                         //      double PofB1 = 1.0;
 375                         //      if( bnn1 > 1.4 )PofB1=0.0;
 376
 377                         double ppslope = 19.0;
 378                         double GammaProfile = exp(-bnn0*bnn0/(2.*hbarc*hbarc*ppslope));
 379                         double PofB0 = 1. - (1. - GammaProfile)*(1. - GammaProfile);
 380                         GammaProfile = exp(-bnn1*bnn1/(2.*hbarc*hbarc*ppslope));
 381                         double PofB1 = 1. - (1. - GammaProfile)*(1. - GammaProfile);
 382
 383                         double Xarg = Egamma*bnn0/(hbarc*_beamLorentzGamma);
 384                         double loc_nofe0 = (_bbs.beam1().Z()*_bbs.beam1().Z()*alpha)/
 385                                 (pi*pi);
 386                         loc_nofe0 *= (1./(Egamma*bnn0*bnn0));
 387                         loc_nofe0 *= Xarg*Xarg*(bessel::dbesk1(Xarg))*(bessel::dbesk1(Xarg));
 388
 389                         Xarg = Egamma*bnn1/(hbarc*_beamLorentzGamma);
 390                         double loc_nofe1 = (_bbs.beam1().Z()*_bbs.beam1().Z()*alpha)/
 391                                 (pi*pi);
 392                         loc_nofe1 *= (1./(Egamma*bnn1*bnn1));
 393                         loc_nofe1 *= Xarg*Xarg*(bessel::dbesk1(Xarg))*(bessel::dbesk1(Xarg));
 394
 395                         local_sum += loc_nofe0*(1. - PofB0)*bnn0*db;
 396                         local_sum += loc_nofe1*(1. - PofB1)*bnn1*db;
 397
 398                 }
 399                 // End Impact parameter loop
 400
 401                 // Note: 2*pi --> pi because of no factor 2 above
 402                 double flux_r=local_sum*pi;
 403                 return flux_r;
 404
 405                 //    bmin = nuclearRadius+nuclearRadius;
 406                 //    flux_r = nepoint(Egamma,bmin);
 407                 //    return flux_r;
 408         }
 409
 410         //   first call?  - initialize - calculate photon flux
 411         Icheck=Icheck+1;
 412         if(Icheck > 1) goto L1000f;
 413
 414         rZ=double(_bbs.beam1().Z());
 415         rA=double(_bbs.beam1().A());
 416         rZ2=double(_bbs.beam2().Z());  //Sergey--dAu
 417         rA2=double(_bbs.beam2().A());  //Sergey
 418
 419         //  Nuclear breakup is done by PofB
 420         //  collect number of integration steps here, in one place
 421
 422         nbstep=1200;
 423         nrstep=60;
 424         nphistep=40;
 425
 426         //  this last one is the number of energy steps
 427         nstep=100;
 428
 429         // following previous choices, take Emin=10 keV at LHC, Emin = 1 MeV at RHIC
 430         Emin=1.E-5;
 431         if (_beamLorentzGamma < 500)
 432                 Emin=1.E-3;
 433
 434         //  maximum energy is 6 times the cutoff
 435         //      Emax=12.*hbarc*_beamLorentzGamma/RNuc;
 436         Emax=6.*hbarc*_beamLorentzGamma/maxradius;
 437
 438         //     >> lnEmin <-> ln(Egamma) for the 0th bin
 439         //     >> lnEmax <-> ln(Egamma) for the last bin
 440
 441         lnEmin=log(Emin);
 442         lnEmax=log(Emax);
 443         dlnE=(lnEmax-lnEmin)/nstep;
 444
 445         cout<<" Calculating flux for photon energies from E= "<<Emin
 446             <<" to  "<<Emax<<"  GeV (CM frame) "<<endl;
 447
 448
 449         stepmult= exp(log(Emax/Emin)/double(nstep));
 450         energy=Emin;
 451
 452         for (int j = 1; j<=nstep;j++){
 453                 energy=energy*stepmult;
 454
 455                 //  integrate flux over 2R_A < b < 2R_A+ 6* gamma hbar/energy
 456                 //  use exponential steps
 457
 458                 bmin=RNuc+RNuc2; //2.*nuclearRadius; Sergey
 459                 bmax=bmin + 6.*hbarc*_beamLorentzGamma/energy;
 460
 461                 bmult=exp(log(bmax/bmin)/double(nbstep));
 462                 biter=bmin;
 463                 integratedflux=0.;
 464
 465                 if (_bbs.beam2().Z()==1&&_bbs.beam1().A()==2){
 466                      //This is for deuteron-gold
 467                      Xvar = (RNuc+RNuc2)*energy/(hbarc*(_beamLorentzGamma));
 468
 469                      fluxelement = (2.0/pi)*rZ*rZ*alpha/
 470                          energy*(Xvar*bessel::dbesk0(Xvar)*bessel::dbesk1(Xvar)-(1/2)*Xvar*Xvar*
 471                                 (bessel::dbesk1(Xvar)*bessel::dbesk1(Xvar)-bessel::dbesk0(Xvar)*bessel::dbesk0(Xvar)));
 472
 473                      integratedflux=integratedflux+fluxelement;
 474                 }else if( (_bbs.beam1().A() == 1 && _bbs.beam2().A() != 1) || (_bbs.beam2().A() == 1 && _bbs.beam1().A() != 1) ){
 475                     // This is pA
 476                     if( _productionMode == PHOTONPOMERONINCOHERENT ){
 477                       // This is pA incoherent
 478                       // cout<<" This is incoherent! "<<" j = "<<j<<endl;
 479                       double zproj = 0.0;
 480                       double localbmin = 0.0;
 481                       if( _bbs.beam1().A() == 1 ){
 482                         zproj = (_bbs.beam2().Z());
 483                         localbmin = _bbs.beam2().nuclearRadius() + 0.7;
 484                       }
 485                       if( _bbs.beam2().A() == 1 ){
 486                         zproj = (_bbs.beam1().Z());
 487                         localbmin = _bbs.beam1().nuclearRadius() + 0.7;
 488                       }
 489                       integratedflux = zproj*zproj*nepoint(energy,localbmin);
 490                     } else if ( _productionMode == PHOTONPOMERONNARROW ||  _productionMode == PHOTONPOMERONWIDE ){
 491                       // cout<<" This is pA coherent "<<" j= "<<j<<endl;
 492                       double localbmin = 0.0;
 493                       if( _bbs.beam1().A() == 1 ){
 494                         localbmin = _bbs.beam2().nuclearRadius() + 0.7;
 495                       }
 496                       if( _bbs.beam2().A() == 1 ){
 497                         localbmin = _bbs.beam1().nuclearRadius() + 0.7;
 498                       }
 499                       integratedflux = nepoint(energy,localbmin);
 500                     }
 501                 }else{
 502                         for (int jb = 1; jb<=nbstep;jb++){
 503                                 bold=biter;
 504                                 biter=biter*bmult;
 505                                 // When we get to b>20R_A change methods - just take the photon flux
 506                                 //  at the center of the nucleus.
 507                                 if (biter > (10.*RNuc))
 508                                         {
 509                                                 // if there is no nuclear breakup or only hadronic breakup, which only
 510                                                 // occurs at smaller b, we can analytically integrate the flux from b~20R_A
 511                                                 // to infinity, following Jackson (2nd edition), Eq. 15.54
 512                                                 Xvar=energy*biter/(hbarc*_beamLorentzGamma);
 513                                                 // Here, there is nuclear breakup.  So, we can't use the integrated flux
 514                                                 // However, we can do a single flux calculation, at the center of the nucleus
 515                                                 // Eq. 41 of Vidovic, Greiner and Soff, Phys.Rev.C47,2308(1993), among other places
 516                                                 // this is the flux per unit area
 517                                                 fluxelement  = (rZ*rZ*alpha*energy)*
 518                                                         (bessel::dbesk1(Xvar))*(bessel::dbesk1(Xvar))/
 519                                                         ((pi*_beamLorentzGamma*hbarc)*
 520                                                          (pi*_beamLorentzGamma*hbarc));
 521
 522                                         }//if biter>10
 523                                 else{
 524                                         // integrate over nuclear surface. n.b. this assumes total shadowing -
 525                                         // treat photons hitting the nucleus the same no matter where they strike
 526                                         fluxelement=0.;
 527                                         deltar=RNuc/double(nrstep);
 528                                         riter=-deltar/2.;
 529
 530                                         for (int jr =1; jr<=nrstep;jr++){
 531                                                 riter=riter+deltar;
 532                                                 // use symmetry;  only integrate from 0 to pi (half circle)
 533                                                 deltaphi=pi/double(nphistep);
 534                                                 phiiter=0.;
 535
 536                                                 for( int jphi=1;jphi<= nphistep;jphi++){
 537                                                         phiiter=(double(jphi)-0.5)*deltaphi;
 538                                                         //  dist is the distance from the center of the emitting nucleus to the point in question
 539                                                         dist=sqrt((biter+riter*cos(phiiter))*(biter+riter*
 540                                                                    cos(phiiter))+(riter*sin(phiiter))*(riter*sin(phiiter)));
 541                                                         Xvar=energy*dist/(hbarc*_beamLorentzGamma);
 542                                                         flux_r = (rZ*rZ*alpha*energy)*
 543                                                                 (bessel::dbesk1(Xvar)*bessel::dbesk1(Xvar))/
 544                                                                 ((pi*_beamLorentzGamma*hbarc)*
 545                                                                  (pi*_beamLorentzGamma*hbarc));
 546
 547                                                         //  The surface  element is 2.* delta phi* r * delta r
 548                                                         //  The '2' is because the phi integral only goes from 0 to pi
 549                                                         fluxelement=fluxelement+flux_r*2.*deltaphi*riter*deltar;
 550                                                         //  end phi and r integrations
 551                                                 }//for(jphi)
 552                                         }//for(jr)
 553                                         //  average fluxelement over the nuclear surface
 554                                         fluxelement=fluxelement/(pi*RNuc*RNuc);
 555                                 }//else
 556                                 //  multiply by volume element to get total flux in the volume element
 557                                 fluxelement=fluxelement*2.*pi*biter*(biter-bold);
 558                                 //  modulate by the probability of nuclear breakup as f(biter)
 559                                 if (_beamBreakupMode > 1){
 560                                         fluxelement=fluxelement*_bbs.probabilityOfBreakup(biter);
 561                                 }
 562                                 integratedflux=integratedflux+fluxelement;
 563
 564                         } //end loop over impact parameter
 565                 }  //end of else (pp, pA, AA)
 566
 567                 //  In lookup table, store k*dN/dk because it changes less
 568                 //  so the interpolation should be better
 569                 dide[j]=integratedflux*energy;
 570
 571         }//end loop over photon energy
 572
 573         //  for 2nd and subsequent calls, use lookup table immediately
 574
 575  L1000f:
 576
 577         lEgamma=log(Egamma);
 578         if (lEgamma < (lnEmin+dlnE) ||  lEgamma  > lnEmax){
 579                 flux_r=0.0;
 580                 cout<<"  ERROR: Egamma outside defined range. Egamma= "<<Egamma
 581                     <<"   "<<lnEmax<<" "<<(lnEmin+dlnE)<<endl;
 582         }
 583         else{
 584                 //       >> Egamma between Ilt and Ilt+1
 585                 Ilt = int((lEgamma-lnEmin)/dlnE);
 586                 //       >> ln(Egamma) for first point
 587                 lnElt = lnEmin + Ilt*dlnE;
 588                 //       >> Interpolate
 589                 flux_r = dide[Ilt] + ((lEgamma-lnElt)/dlnE)*(dide[Ilt+1]- dide[Ilt]);
 590                 flux_r = flux_r/Egamma;
 591         }
 592
 593         return flux_r;
 594 }
 595
 596
 597 //______________________________________________________________________________
 598 double
 599 photonNucleusCrossSection::nepoint(const double Egamma,
 600                                    const double bmin)
 601 {
 602         // Function for the spectrum of virtual photons,
 603         // dn/dEgamma, for a point charge q=Ze sweeping
 604         // past the origin with velocity gamma
 605         // (=1/SQRT(1-(V/c)**2)) integrated over impact
 606         // parameter from bmin to infinity
 607         // See Jackson eq15.54 Classical Electrodynamics
 608         // Declare Local Variables
 609         double beta,X,C1,bracket,nepoint_r;
 610
 611         beta = sqrt(1.-(1./(_beamLorentzGamma*_beamLorentzGamma)));
 612         X = (bmin*Egamma)/(beta*_beamLorentzGamma*hbarc);
 613
 614         bracket = -0.5*beta*beta*X*X*(bessel::dbesk1(X)*bessel::dbesk1(X)
 615                                       -bessel::dbesk0(X)*bessel::dbesk0(X));
 616
 617         bracket = bracket+X*bessel::dbesk0(X)*bessel::dbesk1(X);
 618
 619         //      C1=(2.*double((_bbs.beam1().Z())*(_bbs.beam1().Z()))*
 620         //    alpha)/pi;
 621
 622         // Note: NO  Z*Z!!
 623         C1=(2.*alpha)/pi;
 624
 625         nepoint_r = C1*(1./beta)*(1./beta)*(1./Egamma)*bracket;
 626
 627         return nepoint_r;
 628
 629 }
 630
 631
 632 //______________________________________________________________________________
 633 double
 634 photonNucleusCrossSection::sigmagp(const double Wgp)
 635 {
 636         //     >> Function for the gamma-proton --> VectorMeson
 637         //     >> cross section. Wgp is the gamma-proton CM energy.
 638         //     >> Unit for cross section: fm**2
 639
 640         double sigmagp_r=0.;
 641
 642         switch(_particleType)
 643                 {
 644                 case RHO:
 645                 case RHOZEUS:
 646                 case FOURPRONG:
 647                         sigmagp_r=1.E-4*(5.0*exp(0.22*log(Wgp))+26.0*exp(-1.23*log(Wgp)));
 648                         break;
 649                 case OMEGA:
 650                         sigmagp_r=1.E-4*(0.55*exp(0.22*log(Wgp))+18.0*exp(-1.92*log(Wgp)));
 651                         break;
 652                 case PHI:
 653                         sigmagp_r=1.E-4*0.34*exp(0.22*log(Wgp));
 654                         break;
 655                 case JPSI:
 656                 case JPSI_ee:
 657                 case JPSI_mumu:
 658                         sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
 659                         sigmagp_r*=sigmagp_r;
 660                         sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp));
 661                         // sigmagp_r=1.E-4*0.0015*exp(0.80*log(Wgp));
 662                         break;
 663                 case JPSI2S:
 664                 case JPSI2S_ee:
 665                 case JPSI2S_mumu:
 666                         sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
 667                         sigmagp_r*=sigmagp_r;
 668                         sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp));
 669                         sigmagp_r*=0.166;
 670                         //      sigmagp_r=0.166*(1.E-4*0.0015*exp(0.80*log(Wgp)));
 671                         break;
 672                 case UPSILON:
 673                 case UPSILON_ee:
 674                 case UPSILON_mumu:
 675                         //       >> This is W**1.7 dependence from QCD calculations
 676                         sigmagp_r=1.E-10*(0.060)*exp(1.70*log(Wgp));
 677                         break;
 678                 case UPSILON2S:
 679                 case UPSILON2S_ee:
 680                 case UPSILON2S_mumu:
 681                         sigmagp_r=1.E-10*(0.0259)*exp(1.70*log(Wgp));
 682                         break;
 683                 case UPSILON3S:
 684                 case UPSILON3S_ee:
 685                 case UPSILON3S_mumu:
 686                         sigmagp_r=1.E-10*(0.0181)*exp(1.70*log(Wgp));
 687                         break;
 688                 default: cout<< "!!!  ERROR: Unidentified Vector Meson: "<< _particleType <<endl;
 689                 }
 690         return sigmagp_r;
 691 }
 692
 693
 694 //______________________________________________________________________________
 695 double
 696 photonNucleusCrossSection::sigma_A(const double sig_N)
 697 {
 698         // Nuclear Cross Section
 699         // sig_N,sigma_A in (fm**2)
 700
 701         double sum;
 702         double b,bmax,Pint,arg,sigma_A_r;
 703
 704         int NGAUSS;
 705
 706         double xg[17]=
 707                 {.0,
 708                  .0483076656877383162,.144471961582796493,
 709                  .239287362252137075, .331868602282127650,
 710                  .421351276130635345, .506899908932229390,
 711                  .587715757240762329, .663044266930215201,
 712                  .732182118740289680, .794483795967942407,
 713                  .849367613732569970, .896321155766052124,
 714                  .934906075937739689, .964762255587506430,
 715                  .985611511545268335, .997263861849481564
 716                 };
 717
 718         double ag[17]=
 719                 {.0,
 720                  .0965400885147278006, .0956387200792748594,
 721                  .0938443990808045654, .0911738786957638847,
 722                  .0876520930044038111, .0833119242269467552,
 723                  .0781938957870703065, .0723457941088485062,
 724                  .0658222227763618468, .0586840934785355471,
 725                  .0509980592623761762, .0428358980222266807,
 726                  .0342738629130214331, .0253920653092620595,
 727                  .0162743947309056706, .00701861000947009660
 728                 };
 729
 730         NGAUSS=16;
 731
 732         // Check if one or both beams are nuclei
 733         int A_1 = _bbs.beam1().A();
 734         int A_2 = _bbs.beam2().A();
 735         if( A_1 == 1 && A_2 == 1)cout<<" This is pp, you should not be here..."<<endl;
 736
 737         // CALCULATE P(int) FOR b=0.0 - bmax (fm)
 738         bmax = 25.0;
 739         sum  = 0.;
 740         for(int IB=1;IB<=NGAUSS;IB++){
 741
 742                 b = 0.5*bmax*xg[IB]+0.5*bmax;
 743
 744                 if( A_1 == 1 && A_2 != 1){
 745                   arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b);
 746                 }else if(A_2 == 1 && A_1 != 1){
 747                   arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
 748                 }else{
 749                   arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
 750                 }
 751
 752                 Pint=1.0-exp(arg);
 753                 sum=sum+2.*pi*b*Pint*ag[IB];
 754
 755
 756                 b = 0.5*bmax*(-xg[IB])+0.5*bmax;
 757
 758                 if( A_1 == 1 && A_2 != 1){
 759                   arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b);
 760                 }else if(A_2 == 1 && A_1 != 1){
 761                   arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
 762                 }else{
 763                   arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
 764                 }
 765
 766                 Pint=1.0-exp(arg);
 767                 sum=sum+2.*pi*b*Pint*ag[IB];
 768
 769         }
 770
 771         sum=0.5*bmax*sum;
 772
 773         sigma_A_r=sum;
 774
 775         return sigma_A_r;
 776 }
 777
 778 //______________________________________________________________________________
 779 double
 780 photonNucleusCrossSection::sigma_N(const double Wgp)
 781 {
 782         // Nucleon Cross Section in (fm**2)
 783         double cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha);
 784         return cs;
 785 }
 786
 787
 788 //______________________________________________________________________________
 789 double
 790 photonNucleusCrossSection::breitWigner(const double W,
 791                                        const double C)
 792 {
 793         // use simple fixed-width s-wave Breit-Wigner without coherent backgorund for rho'
 794         // (PDG '08 eq. 38.56)
 795         if(_particleType==FOURPRONG) {
 796                 if (W < 4.01 * pionChargedMass)
 797                         return 0;
 798                 const double termA  = _channelMass * _width;
 799                 const double termA2 = termA * termA;
 800                 const double termB  = W * W - _channelMass * _channelMass;
 801                 return C * _ANORM * _ANORM * termA2 / (termB * termB + termA2);
 802         }
 803
 804         // Relativistic Breit-Wigner according to J.D. Jackson,
 805         // Nuovo Cimento 34, 6692 (1964), with nonresonant term. A is the strength
 806         // of the resonant term and b the strength of the non-resonant
 807         // term. C is an overall normalization.
 808
 809         double ppi=0.,ppi0=0.,GammaPrim,rat;
 810         double aa,bb,cc;
 811
 812         double nrbw_r;
 813
 814         // width depends on energy - Jackson Eq. A.2
 815         // if below threshold, then return 0.  Added 5/3/2001 SRK
 816         // 0.5% extra added for safety margin
 817         if( _particleType==RHO ||_particleType==RHOZEUS){
 818                 if (W < 2.01*pionChargedMass){
 819                         nrbw_r=0.;
 820                         return nrbw_r;
 821                 }
 822                 ppi=sqrt( ((W/2.)*(W/2.)) - pionChargedMass * pionChargedMass);
 823                 ppi0=0.358;
 824         }
 825
 826         // handle phi-->K+K- properly
 827         if (_particleType  ==  PHI){
 828                 if (W < 2.*kaonChargedMass){
 829                         nrbw_r=0.;
 830                         return nrbw_r;
 831                 }
 832                 ppi=sqrt( ((W/2.)*(W/2.))- kaonChargedMass*kaonChargedMass);
 833                 ppi0=sqrt( ((_channelMass/2.)*(_channelMass/2.))-kaonChargedMass*kaonChargedMass);
 834         }
 835
 836         //handle J/Psi-->e+e- properly
 837         if (_particleType==JPSI || _particleType==JPSI2S){
 838                 if(W<2.*mel){
 839                         nrbw_r=0.;
 840                         return nrbw_r;
 841                 }
 842                 ppi=sqrt(((W/2.)*(W/2.))-mel*mel);
 843                 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel);
 844         }
 845         if (_particleType==JPSI_ee){
 846                 if(W<2.*mel){
 847                         nrbw_r=0.;
 848                         return nrbw_r;
 849                 }
 850                 ppi=sqrt(((W/2.)*(W/2.))-mel*mel);
 851                 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel);
 852         }
 853         if (_particleType==JPSI_mumu){
 854                 if(W<2.*muonMass){
 855                         nrbw_r=0.;
 856                         return nrbw_r;
 857                 }
 858                 ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass);
 859                 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass);
 860         }
 861         if (_particleType==JPSI2S_ee){
 862                 if(W<2.*mel){
 863                         nrbw_r=0.;
 864                         return nrbw_r;
 865                 }
 866                 ppi=sqrt(((W/2.)*(W/2.))-mel*mel);
 867                 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel);
 868         }
 869         if (_particleType==JPSI2S_mumu){
 870                 if(W<2.*muonMass){
 871                         nrbw_r=0.;
 872                         return nrbw_r;
 873                 }
 874                 ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass);
 875                 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass);
 876         }
 877
 878         if(_particleType==UPSILON || _particleType==UPSILON2S ||_particleType==UPSILON3S ){
 879                 if (W<2.*muonMass){
 880                         nrbw_r=0.;
 881                         return nrbw_r;
 882                 }
 883                 ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass);
 884                 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass);
 885         }
 886
 887         if(_particleType==UPSILON_mumu || _particleType==UPSILON2S_mumu ||_particleType==UPSILON3S_mumu ){
 888                 if (W<2.*muonMass){
 889                         nrbw_r=0.;
 890                         return nrbw_r;
 891                 }
 892                 ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass);
 893                 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass);
 894         }
 895
 896         if(_particleType==UPSILON_ee || _particleType==UPSILON2S_ee ||_particleType==UPSILON3S_ee ){
 897                 if (W<2.*mel){
 898                         nrbw_r=0.;
 899                         return nrbw_r;
 900                 }
 901                 ppi=sqrt(((W/2.)*(W/2.))-mel*mel);
 902                 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel);
 903         }
 904
 905         if(ppi==0.&&ppi0==0.)
 906                 cout<<"Improper Gammaacrosssection::breitwigner, ppi&ppi0=0."<<endl;
 907
 908         rat=ppi/ppi0;
 909         GammaPrim=_width*(_channelMass/W)*rat*rat*rat;
 910
 911         aa=_ANORM*sqrt(GammaPrim*_channelMass*W);
 912         bb=W*W-_channelMass*_channelMass;
 913         cc=_channelMass*GammaPrim;
 914
 915         // First real part squared
 916         nrbw_r = (( (aa*bb)/(bb*bb+cc*cc) + _BNORM)*( (aa*bb)/(bb*bb+cc*cc) + _BNORM));
 917
 918         // Then imaginary part squared
 919         nrbw_r = nrbw_r + (( (aa*cc)/(bb*bb+cc*cc) )*( (aa*cc)/(bb*bb+cc*cc) ));
 920
 921         //  Alternative, a simple, no-background BW, following J. Breitweg et al.
 922         //  Eq. 15 of Eur. Phys. J. C2, 247 (1998).  SRK 11/10/2000
 923         //      nrbw_r = (_ANORM*_mass*GammaPrim/(bb*bb+cc*cc))**2
 924
 925         nrbw_r = C*nrbw_r;
 926
 927         return nrbw_r;
 928 }