STARLIGHT/starlight/src/gammaaluminosity.cpp

   1 ///////////////////////////////////////////////////////////////////////////
   2 //
   3 //    Copyright 2010
   4 //
   5 //    This file is part of starlight.
   6 //
   7 //    starlight is free software: you can redistribute it and/or modify
   8 //    it under the terms of the GNU General Public License as published by
   9 //    the Free Software Foundation, either version 3 of the License, or
  10 //    (at your option) any later version.
  11 //
  12 //    starlight is distributed in the hope that it will be useful,
  13 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 //    GNU General Public License for more details.
  16 //
  17 //    You should have received a copy of the GNU General Public License
  18 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
  19 //
  20 ///////////////////////////////////////////////////////////////////////////
  21 //
  22 // File and Version Information:
  23 // $Rev:: 164                         $: revision of last commit
  24 // $Author:: odjuvsla                 $: author of last commit
  25 // $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
  26 //
  27 // Description:
  28 //
  29 //
  30 //
  31 ///////////////////////////////////////////////////////////////////////////
  32
  33
  34 #include <iostream>
  35 #include <fstream>
  36 #include <cmath>
  37
  38 #include "inputParameters.h"
  39 #include "beambeamsystem.h"
  40 #include "beam.h"
  41 #include "starlightconstants.h"
  42 #include "nucleus.h"
  43 #include "bessel.h"
  44 #include "gammaaluminosity.h"
  45
  46
  47 using namespace std;
  48
  49
  50 //______________________________________________________________________________
  51 photonNucleusLuminosity::photonNucleusLuminosity(beamBeamSystem& bbsystem)
  52   : photonNucleusCrossSection(bbsystem)
  53   ,_nPtBinsInterference(inputParametersInstance.nmbPtBinsInterference())
  54   ,_ptBinWidthInterference(inputParametersInstance.ptBinWidthInterference())
  55   ,_interferenceStrength(inputParametersInstance.interferenceStrength())
  56 {
  57   cout <<"Creating Luminosity Tables."<<endl;
  58   photonNucleusDifferentialLuminosity();
  59   cout <<"Luminosity Tables created."<<endl;
  60 }
  61
  62
  63 //______________________________________________________________________________
  64 photonNucleusLuminosity::~photonNucleusLuminosity()
  65 { }
  66
  67
  68 //______________________________________________________________________________
  69 void photonNucleusLuminosity::photonNucleusDifferentialLuminosity()
  70 {
  71         // double Av,Wgp,cs,cvma;
  72   double W,dW,dY;
  73   double Egamma,Y;
  74   // double t,tmin,tmax;
  75   double testint,dndWdY;
  76   double csgA;
  77   // double ax,bx;
  78   double C;
  79
  80   ofstream wylumfile;
  81   wylumfile.precision(15);
  82
  83   double  bwnorm,Eth;
  84
  85   dW = (_wMax-_wMin)/_nWbins;
  86   dY  = (_yMax-(-1.0)*_yMax)/_nYbins;
  87
  88   // Write the values of W used in the calculation to slight.txt.
  89   wylumfile.open("slight.txt");
  90   wylumfile << inputParametersInstance.parameterValueKey() << endl;
  91   wylumfile << getbbs().beam1().Z() <<endl;
  92   wylumfile << getbbs().beam1().A() <<endl;
  93   wylumfile << getbbs().beam2().Z() <<endl;
  94   wylumfile << getbbs().beam2().A() <<endl;
  95   wylumfile << inputParametersInstance.beamLorentzGamma() <<endl;
  96   wylumfile << inputParametersInstance.maxW() <<endl;
  97   wylumfile << inputParametersInstance.minW() <<endl;
  98   wylumfile << inputParametersInstance.nmbWBins() <<endl;
  99   wylumfile << inputParametersInstance.maxRapidity() <<endl;
 100   wylumfile << inputParametersInstance.nmbRapidityBins() <<endl;
 101   wylumfile << inputParametersInstance.productionMode() <<endl;
 102   wylumfile << inputParametersInstance.beamBreakupMode() <<endl;
 103   wylumfile << inputParametersInstance.interferenceEnabled() <<endl;
 104   wylumfile << inputParametersInstance.interferenceStrength() <<endl;
 105   wylumfile << inputParametersInstance.coherentProduction() <<endl;
 106   wylumfile << inputParametersInstance.incoherentFactor() <<endl;
 107   wylumfile << inputParametersInstance.deuteronSlopePar() <<endl;
 108   wylumfile << inputParametersInstance.maxPtInterference() <<endl;
 109   wylumfile << inputParametersInstance.nmbPtBinsInterference() <<endl;
 110
 111   //     Normalize the Breit-Wigner Distribution and write values of W to slight.txt
 112   testint=0.0;
 113   //Grabbing default value for C in the breit-wigner calculation
 114   C=getDefaultC();
 115   for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
 116     W = _wMin + double(i)*dW + 0.5*dW;
 117     testint = testint + breitWigner(W,C)*dW;
 118     wylumfile << W << endl;
 119   }
 120   bwnorm = 1./testint;
 121
 122   //     Write the values of Y used in the calculation to slight.txt.
 123   for(unsigned int i = 0; i <= _nYbins - 1; ++i) {
 124     Y = -1.0*_yMax + double(i)*dY + 0.5*dY;
 125     wylumfile << Y << endl;
 126   }
 127
 128   Eth=0.5*(((_wMin+starlightConstants::protonMass)*(_wMin
 129                                                             +starlightConstants::protonMass)-starlightConstants::protonMass*starlightConstants::protonMass)/
 130            (inputParametersInstance.protonEnergy()+sqrt(inputParametersInstance.protonEnergy()*
 131                                                inputParametersInstance.protonEnergy()-starlightConstants::protonMass*starlightConstants::protonMass)));
 132
 133   for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
 134
 135     W = _wMin + double(i)*dW + 0.5*dW;
 136
 137     for(unsigned int j = 0; j <= _nYbins - 1; ++j) {
 138
 139       Y = -1.0*_yMax + double(j)*dY + 0.5*dY;
 140
 141       int A_1 = getbbs().beam1().A();
 142       int A_2 = getbbs().beam2().A();
 143       if( A_2 == 1 && A_1 != 1 ){
 144         // pA, first beam is the nucleus
 145         Egamma = 0.5*W*exp(-Y);
 146       } else if( A_1 ==1 && A_2 != 1){
 147         // pA, second beam is the nucleus
 148         Egamma = 0.5*W*exp(Y);
 149       } else {
 150         Egamma = 0.5*W*exp(Y);
 151       }
 152
 153       dndWdY = 0.;
 154
 155       if( Egamma > Eth && Egamma < maxPhotonEnergy() ){
 156
 157        csgA=getcsgA(Egamma,W);
 158        dndWdY = Egamma*photonFlux(Egamma)*csgA*breitWigner(W,bwnorm);
 159
 160       }
 161
 162       wylumfile << dndWdY << endl;
 163
 164     }
 165   }
 166
 167   wylumfile.close();
 168
 169   if(inputParametersInstance.interferenceEnabled()==1)
 170     pttablegen();
 171
 172   wylumfile.open("slight.txt",ios::app);
 173   cout << "bwnorm: "<< bwnorm <<endl;
 174   wylumfile << bwnorm << endl;
 175   wylumfile << inputParametersInstance.parameterValueKey() << endl;
 176   wylumfile.close();
 177 }
 178
 179
 180 //______________________________________________________________________________
 181 void photonNucleusLuminosity::pttablegen()
 182 {
 183   //  Calculates the pt spectra for VM production with interference
 184   //  Follows S. Klein and J. Nystrand, Phys. Rev Lett. 84, 2330 (2000).
 185   //  Written by S. Klein, 8/2002
 186
 187   //  fill in table pttable in one call
 188   //  Integrate over all y (using the same y values as in table yarray
 189   //  note that the cross section goes from ymin (<0) to ymax (>0), in numy points
 190   //  here,  we go from 0 to ymax in (numy/2)+1 points
 191   //  numy must be even.
 192
 193   //  At each y, calculate the photon energies Egamma1 and Egamma2
 194   //  and the two photon-A cross sections
 195
 196   //  loop over each p_t entry in the table.
 197
 198   //  Then, loop over b and phi (the angle between the VM \vec{p_t} and \vec{b}
 199   //  and calculate the cross section at each step.
 200   //  Put the results in pttable
 201
 202   ofstream wylumfile;
 203   wylumfile.precision(15);
 204   wylumfile.open("slight.txt",ios::app);
 205
 206
 207   double param1pt[500],param2pt[500];
 208   double  *ptparam1=param1pt;
 209   double  *ptparam2=param2pt;
 210   double dY=0.,Yp=0.,Egamma1=0.,Egamma2=0.,Wgp=0.,cs=0.,cvma=0.,Av=0.,tmin=0.,tmax=0.,ax=0.,bx=0.;
 211   double csgA=0.,t=0.,sig_ga_1=0.,sig_ga_2=0.,bmax=0.,bmin=0.,db=0.,pt=0.,sum1=0.,b=0.,A1=0.,A2=0.;
 212   double sumg=0.,theta=0.,amp_i_2=0.,sumint=0.;
 213   int NGAUSS=0,NBIN=0,NThetaBIN=0;
 214
 215   double xg[16]={.0483076656877383162E0,.144471961582796493E0,
 216                  .239287362252137075E0, .331868602282127650E0,
 217                  .421351276130635345E0, .506899908932229390E0,
 218                  .587715757240762329E0, .663044266930215201E0,
 219                  .732182118740289680E0, .794483795967942407E0,
 220                  .849367613732569970E0, .896321155766052124E0,
 221                  .934906075937739689E0, .964762255587506430E0,
 222                  .985611511545268335E0, .997263861849481564E0};
 223   double ag[16]={.0965400885147278006E0, .0956387200792748594E0,
 224                  .0938443990808045654E0, .0911738786957638847E0,
 225                  .0876520930044038111E0, .0833119242269467552E0,
 226                  .0781938957870703065E0, .0723457941088485062E0,
 227                  .0658222227763618468E0, .0586840934785355471E0,
 228                  .0509980592623761762E0, .0428358980222266807E0,
 229                  .0342738629130214331E0, .0253920653092620595E0,
 230                  .0162743947309056706E0, .00701861000947009660E0};
 231
 232   NGAUSS=16;
 233
 234   //Setting input calls to variables/less calls this way.
 235   double Ymax=_yMax;
 236   int numy = _nYbins;
 237   double Ep = inputParametersInstance.protonEnergy();
 238   int ibreakup = inputParametersInstance.beamBreakupMode();
 239   double NPT = inputParametersInstance.nmbPtBinsInterference();
 240   double gamma_em=inputParametersInstance.beamLorentzGamma();
 241   double mass= getChannelMass();
 242
 243   //  loop over y from 0 (not -ymax) to yma
 244
 245   dY=(2.*Ymax)/numy;
 246   for(int jy=1;jy<=numy/2;jy++){
 247     Yp=(double(jy)-0.5)*dY;
 248
 249     // Find the photon energies.  Yp >= 0, so Egamma2 is smaller
 250     // Use the vector meson mass for W here - neglect the width
 251
 252     Egamma1 = 0.5*mass*exp(Yp);
 253     Egamma2 = 0.5*mass*exp(-Yp);
 254
 255     //  Find the sigma(gammaA) for the two directions
 256     //  Photonuclear Cross Section 1
 257     //  Gamma-proton CM energy
 258
 259     Wgp=sqrt(2.*Egamma1*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass*
 260                                  starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
 261
 262     // Calculate V.M.+proton cross section
 263
 264     cs=sqrt(16.*starlightConstants::pi*vmPhotonCoupling()*slopeParameter()*
 265             starlightConstants::hbarc*starlightConstants::hbarc*sigmagp(Wgp)
 266             /starlightConstants::alpha);
 267
 268     // Calculate V.M.+nucleus cross section
 269
 270     cvma=sigma_A(cs);
 271
 272     // Calculate Av = dsigma/dt(t=0) Note Units: fm**2/Gev**2
 273
 274     Av=(starlightConstants::alpha*cvma*cvma)/(16.*starlightConstants::pi
 275                                               *vmPhotonCoupling()*starlightConstants::hbarc*starlightConstants::hbarc);
 276
 277     tmin  = ((mass*mass)/(4.*Egamma1*gamma_em)*(mass*mass)/(4.*Egamma1*gamma_em));
 278     tmax  = tmin + 0.25;
 279     ax    = 0.5*(tmax-tmin);
 280     bx    = 0.5*(tmax+tmin);
 281     csgA  = 0.;
 282
 283     for(int k=0;k<NGAUSS;k++){
 284       t     = sqrt(ax*xg[k]+bx);
 285       csgA  = csgA + ag[k]*getbbs().beam2().formFactor(t)*getbbs().beam2().formFactor(t);
 286       t     = sqrt(ax*(-xg[k])+bx);
 287       csgA  = csgA + ag[k]*getbbs().beam2().formFactor(t)*getbbs().beam2().formFactor(t);
 288     }
 289
 290     csgA = 0.5*(tmax-tmin)*csgA;
 291     csgA = Av*csgA;
 292     sig_ga_1 = csgA;
 293
 294     // Photonuclear Cross Section 2
 295
 296     Wgp=sqrt(2.*Egamma2*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass*
 297                                  starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
 298
 299     cs=sqrt(16.*starlightConstants::pi*vmPhotonCoupling()*slopeParameter()*
 300             starlightConstants::hbarc*starlightConstants::hbarc*sigmagp(Wgp)/starlightConstants::alpha);
 301
 302     cvma=sigma_A(cs);
 303
 304     Av=(starlightConstants::alpha*cvma*cvma)/(16.*starlightConstants::pi
 305                                               *vmPhotonCoupling()*starlightConstants::hbarc*starlightConstants::hbarc);
 306
 307     tmin  = (((mass*mass)/(4.*Egamma2*gamma_em))*((mass*mass)/(4.*Egamma2*gamma_em)));
 308     tmax  = tmin + 0.25;
 309     ax    = 0.5*(tmax-tmin);
 310     bx    = 0.5*(tmax+tmin);
 311     csgA  = 0.;
 312
 313     for(int k=0;k<NGAUSS;k++){
 314       t     = sqrt(ax*xg[k]+bx);
 315       csgA  = csgA + ag[k]*getbbs().beam2().formFactor(t)*getbbs().beam2().formFactor(t);
 316       t     = sqrt(ax*(-xg[k])+bx);
 317       csgA  = csgA + ag[k]*getbbs().beam2().formFactor(t)*getbbs().beam2().formFactor(t);
 318     }
 319
 320     csgA = 0.5*(tmax-tmin)*csgA;
 321     csgA = Av*csgA;
 322     sig_ga_2 = csgA;
 323
 324     //  Set up pttables - they find the reduction in sigma(pt)
 325     //  due to the nuclear form factors.
 326     //  Use the vector meson mass for W here - neglect width in
 327     //  interference calculation
 328
 329     ptparam1=vmsigmapt(mass,Egamma1,ptparam1);
 330     ptparam2=vmsigmapt(mass,Egamma2,ptparam2);
 331
 332     //  set  bmax according to the smaller photon energy, following flux.f
 333
 334     bmax=bmin+6.*starlightConstants::hbarc*gamma_em/Egamma2;
 335     bmin = getbbs().beam1().nuclearRadius()+getbbs().beam2().nuclearRadius();
 336     //  if we allow for nuclear breakup, use a slightly smaller bmin
 337
 338     if (ibreakup != 1)
 339       bmin=0.95*bmin;
 340     //  set number of bins to a reasonable number to start
 341     NBIN = 2000;
 342     NThetaBIN = 1000;
 343     db   = (bmax-bmin)/float(NBIN);
 344     // loop over pt
 345     for(int i=1;i<=NPT;i++){
 346
 347       pt = (float(i)-0.5)*_ptBinWidthInterference;
 348       sum1=0.0;
 349       // loop over b
 350       for(int j=1;j<=NBIN;j++){
 351
 352         b = bmin + (float(j)-0.5)*db;
 353         //  nofe is the photon flux function
 354         A1 = Egamma1*nofe(Egamma1,b)*sig_ga_1*ptparam1[i];
 355         A2 = Egamma2*nofe(Egamma2,b)*sig_ga_2*ptparam2[i];
 356         sumg=0.0;
 357         //  do this as a Gaussian integral, from 0 to pi
 358         for(int k=0;k<NGAUSS;k++){
 359
 360           theta=xg[k]*starlightConstants::pi;
 361           //  allow for a linear sum of interfering and non-interfering amplitudes
 362           amp_i_2 = A1 + A2 - 2.*_interferenceStrength
 363             *sqrt(A1*A2)*cos(pt*b*cos(theta)/starlightConstants::hbarc);
 364           sumg  = sumg+ag[k]*amp_i_2;
 365         }
 366         //  this is dn/dpt^2
 367         //  The factor of 2 is because the theta integral is only from 0 to pi
 368         sumint=2.*sumg*b*db;
 369         if (ibreakup > 1)
 370           sumint=sumint*getbbs().probabilityOfBreakup(b);
 371         sum1 = sum1 + sumint;
 372       }
 373       //  normalization is done in readDiffLum.f
 374       //  This is d^2sigma/dpt^2; convert to dsigma/dpt
 375       //CHECK THIS OUT---->      write(20,*)sum1*pt*dpt
 376       wylumfile << sum1*pt*_ptBinWidthInterference <<endl;
 377       //  end of pt loop
 378     }
 379     //  end of y loop
 380   }
 381   wylumfile.close();
 382 }
 383
 384
 385 //______________________________________________________________________________
 386 double *photonNucleusLuminosity::vmsigmapt(double W, double Egamma, double *SIGMAPT)
 387 {
 388   //
 389   //  This subroutine calculates the effect of the nuclear form factor
 390   // on the pt spectrum, for use in interference calculations
 391   // For an interaction with mass W and photon energy Egamma,
 392   // it calculates the cross section suppression SIGMAPT(PT)
 393   // as a function of pt.
 394   // The input pt values come from pttable.inc
 395
 396
 397   double pxmax=0.,pymax=0.,dx=0.,Epom=0.,pt=0.,px0=0.,py0=0.,sum=0.,sumy=0.;
 398   double py=0.,px=0.,pt1=0.,pt2=0.,f1=0.,f2=0.,q1=0.,q2=0.,norm=0.;
 399   int NGAUSS =0,Nxbin=0;
 400   double xg[16]={.0483076656877383162e0,.144471961582796493e0,
 401                  .239287362252137075e0, .331868602282127650e0,
 402                  .421351276130635345e0, .506899908932229390e0,
 403                  .587715757240762329e0, .663044266930215201e0,
 404                  .732182118740289680e0, .794483795967942407e0,
 405                  .849367613732569970e0, .896321155766052124e0,
 406                  .934906075937739689e0, .964762255587506430e0,
 407                  .985611511545268335e0, .997263861849481564e0};
 408   double ag[16]={.0965400885147278006e0, .0956387200792748594e0,
 409                  .0938443990808045654e0, .0911738786957638847e0,
 410                  .0876520930044038111e0, .0833119242269467552e0,
 411                  .0781938957870703065e0, .0723457941088485062e0,
 412                  .0658222227763618468e0, .0586840934785355471e0,
 413                  .0509980592623761762e0, .0428358980222266807e0,
 414                  .0342738629130214331e0, .0253920653092620595e0,
 415                  .0162743947309056706e0, .00701861000947009660e0};
 416   NGAUSS=16;
 417
 418   //     >> Initialize
 419   pxmax = 10.*(starlightConstants::hbarc/getbbs().beam1().nuclearRadius());
 420   pymax = 10.*(starlightConstants::hbarc/getbbs().beam1().nuclearRadius());
 421
 422   Nxbin = 500;
 423
 424   dx = 2.*pxmax/double(Nxbin);
 425   Epom   = W*W/(4.*Egamma);
 426
 427   //     >> Loop over total Pt to find distribution
 428
 429       for(int k=1;k<=_nPtBinsInterference;k++){
 430
 431               pt=_ptBinWidthInterference*(double(k)-0.5);
 432
 433               px0 = pt;
 434               py0 = 0.0;
 435
 436               //  For each total Pt, integrate over Pt1, , the photon pt
 437               //  The pt of the Pomeron  is the difference
 438               //  pt1 is
 439               sum=0.;
 440               for(int i=1;i<=Nxbin;i++){
 441
 442                 px = -pxmax + (double(i)-0.5)*dx;
 443                 sumy=0.0;
 444                 for(int j=0;j<NGAUSS;j++){
 445
 446                   py = 0.5*pymax*xg[j]+0.5*pymax;
 447                   //  photon pt
 448                   pt1 = sqrt( px*px + py*py );
 449                   //  pomeron pt
 450                   pt2 = sqrt( (px-px0)*(px-px0) + (py-py0)*(py-py0) );
 451                   q1  = sqrt( ((Egamma/_beamLorentzGamma)*(Egamma/_beamLorentzGamma)) + pt1*pt1 );
 452                   q2  = sqrt( ((Epom/_beamLorentzGamma)*(Epom/_beamLorentzGamma))   + pt2*pt2 );
 453
 454                   //  photon form factor
 455                   // add in phase space factor?
 456                   f1  = (getbbs().beam1().formFactor(q1*q1)*getbbs().beam1().formFactor(q1*q1)*pt1*pt1)/(q1*q1*q1*q1);
 457
 458                   //  Pomeron form factor
 459                   f2  = getbbs().beam1().formFactor(q2*q2)*getbbs().beam1().formFactor(q2*q2);
 460                   sumy= sumy + ag[j]*f1*f2;
 461
 462                   //  now consider other half of py phase space - why is this split?
 463                   py = 0.5*pymax*(-xg[j])+0.5*pymax;
 464                   pt1 = sqrt( px*px + py*py );
 465                   pt2 = sqrt( (px-px0)*(px-px0) + (py-py0)*(py-py0) );
 466                   q1  = sqrt( ((Egamma/_beamLorentzGamma)*Egamma/_beamLorentzGamma) + pt1*pt1 );
 467                   q2  = sqrt( ((Epom/_beamLorentzGamma)*(Epom/_beamLorentzGamma))   + pt2*pt2 );
 468                   //  add in phase space factor?
 469                   f1  = (getbbs().beam1().formFactor(q1*q1)*getbbs().beam1().formFactor(q1*q1)*pt1*pt1)/(q1*q1*q1*q1);
 470                   f2  = getbbs().beam1().formFactor(q2*q2)*getbbs().beam1().formFactor(q2*q2);
 471                   sumy= sumy + ag[j]*f1*f2;
 472
 473                 }
 474                 //         >> This is to normalize the gaussian integration
 475                 sumy = 0.5*pymax*sumy;
 476                 //         >> The 2 is to account for py: 0 -- pymax
 477                 sum  = sum + 2.*sumy*dx;
 478               }
 479
 480               if(k==1) norm=1./sum;
 481               SIGMAPT[k]=sum*norm;
 482       }
 483       return (SIGMAPT);
 484 }
 485
 486
 487 //______________________________________________________________________________
 488 double photonNucleusLuminosity::nofe(double Egamma, double bimp)
 489 {
 490   //Function for the calculation of the "photon density".
 491   //nofe=numberofgammas/(energy*area)
 492   //Assume beta=1.0 and gamma>>1, i.e. neglect the (1/gamma^2)*K0(x) term
 493
 494   double X=0.,nofex=0.,factor1=0.,factor2=0.,factor3=0.;
 495
 496   X = (bimp*Egamma)/(_beamLorentzGamma*starlightConstants::hbarc);
 497
 498   if( X <= 0.0)
 499     cout<<"In nofe, X= "<<X<<endl;
 500
 501   factor1 = (double(getbbs().beam1().Z()*getbbs().beam1().Z())
 502              *starlightConstants::alpha)/(starlightConstants::pi*starlightConstants::pi);
 503
 504   factor2 = 1./(Egamma*bimp*bimp);
 505   factor3 = X*X*(bessel::dbesk1(X))*(bessel::dbesk1(X));
 506   nofex    = factor1*factor2*factor3;
 507   return nofex;
 508 }