estarlight/src/gammaaluminosity.cpp

0001 ///////////////////////////////////////////////////////////////////////////
0002 //
0003 //    Copyright 2010
0004 //
0005 //    This file is part of starlight.
0006 //
0007 //    starlight is free software: you can redistribute it and/or modify
0008 //    it under the terms of the GNU General Public License as published by
0009 //    the Free Software Foundation, either version 3 of the License, or
0010 //    (at your option) any later version.
0011 //
0012 //    starlight is distributed in the hope that it will be useful,
0013 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0015 //    GNU General Public License for more details.
0016 //
0017 //    You should have received a copy of the GNU General Public License
0018 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
0019 //
0020 ///////////////////////////////////////////////////////////////////////////
0021 //
0022 // File and Version Information:
0023 // $Rev:: 264                         $: revision of last commit
0024 // $Author:: jseger                   $: author of last commit
0025 // $Date:: 2016-06-06 21:05:12 +0100 #$: date of last commit
0026 //
0027 // Description:
0028 //
0029 //
0030 //
0031 ///////////////////////////////////////////////////////////////////////////
0032
0033
0034 #include <iostream>
0035 #include <fstream>
0036 #include <cmath>
0037
0038 #include "inputParameters.h"
0039 #include "beambeamsystem.h"
0040 #include "beam.h"
0041 #include "starlightconstants.h"
0042 #include "nucleus.h"
0043 #include "bessel.h"
0044 #include "gammaaluminosity.h"
0045
0046
0047 using namespace std;
0048 using namespace starlightConstants;
0049
0050
0051 //______________________________________________________________________________
0052 photonNucleusLuminosity::photonNucleusLuminosity(const inputParameters& inputParametersInstance, beamBeamSystem& bbsystem)
0053   : photonNucleusCrossSection(inputParametersInstance, bbsystem)
0054   ,_ptBinWidthInterference(inputParametersInstance.ptBinWidthInterference())
0055   ,_interferenceStrength(inputParametersInstance.interferenceStrength())
0056   ,_protonEnergy(inputParametersInstance.protonEnergy())
0057   ,_beamLorentzGamma(inputParametersInstance.beamLorentzGamma())
0058   ,_baseFileName(inputParametersInstance.baseFileName())
0059   ,_maxW(inputParametersInstance.maxW())
0060   ,_minW(inputParametersInstance.minW())
0061   ,_nmbWBins(inputParametersInstance.nmbWBins())
0062   ,_maxRapidity(inputParametersInstance.maxRapidity())
0063   ,_nmbRapidityBins(inputParametersInstance.nmbRapidityBins())
0064   ,_productionMode(inputParametersInstance.productionMode())
0065   ,_beamBreakupMode(inputParametersInstance.beamBreakupMode())
0066   ,_interferenceEnabled(inputParametersInstance.interferenceEnabled())
0067   ,_maxPtInterference(inputParametersInstance.maxPtInterference())
0068   ,_nmbPtBinsInterference(inputParametersInstance.nmbPtBinsInterference())
0069 {
0070   cout <<"Creating Luminosity Tables."<<endl;
0071   photonNucleusDifferentialLuminosity();
0072   cout <<"Luminosity Tables created."<<endl;
0073 }
0074
0075
0076 //______________________________________________________________________________
0077 photonNucleusLuminosity::~photonNucleusLuminosity()
0078 { }
0079
0080
0081 //______________________________________________________________________________
0082 void photonNucleusLuminosity::photonNucleusDifferentialLuminosity()
0083 {
0084   double W,dW,dY;
0085   double Egamma,Y;
0086   double testint,dndWdY;
0087   double csgA;
0088   double C;
0089   int beam;
0090
0091   std::string wyFileName;
0092   wyFileName = _baseFileName +".txt";
0093
0094   ofstream wylumfile;
0095   wylumfile.precision(15);
0096
0097   double  bwnorm,Eth;
0098
0099   dW = (_wMax-_wMin)/_nWbins;
0100   dY  = (_yMax-(-1.0)*_yMax)/_nYbins;
0101
0102   // Write the values of W used in the calculation to slight.txt.
0103   wylumfile.open(wyFileName.c_str());
0104   wylumfile << getbbs().targetBeam().Z() <<endl;
0105   wylumfile << getbbs().targetBeam().A() <<endl;
0106   wylumfile << _beamLorentzGamma <<endl;
0107   wylumfile << _maxW <<endl;
0108   wylumfile << _minW <<endl;
0109   wylumfile << _nmbWBins <<endl;
0110   wylumfile << _maxRapidity <<endl;
0111   wylumfile << _nmbRapidityBins <<endl;
0112   wylumfile << _productionMode <<endl;
0113   wylumfile << _beamBreakupMode <<endl;
0114   wylumfile << _interferenceEnabled <<endl;
0115   wylumfile << _interferenceStrength <<endl;
0116   wylumfile << starlightConstants::deuteronSlopePar <<endl;
0117   wylumfile << _maxPtInterference <<endl;
0118   wylumfile << _nmbPtBinsInterference <<endl;
0119
0120   //     Normalize the Breit-Wigner Distribution and write values of W to slight.txt
0121   testint=0.0;
0122   //Grabbing default value for C in the breit-wigner calculation
0123   C=getDefaultC();
0124   for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
0125     W = _wMin + double(i)*dW + 0.5*dW;
0126     testint = testint + breitWigner(W,C)*dW;
0127     wylumfile << W << endl;
0128   }
0129   bwnorm = 1./testint;
0130
0131   //     Write the values of Y used in the calculation to slight.txt.
0132   for(unsigned int i = 0; i <= _nYbins - 1; ++i) {
0133     Y = -1.0*_yMax + double(i)*dY + 0.5*dY;
0134     wylumfile << Y << endl;
0135   }
0136
0137   Eth=0.5*(((_wMin+starlightConstants::protonMass)*(_wMin +starlightConstants::protonMass)-starlightConstants::protonMass*starlightConstants::protonMass)/(_protonEnergy+sqrt(_protonEnergy*_protonEnergy-starlightConstants::protonMass*starlightConstants::protonMass)));
0138
0139   int A_1 = 0;
0140   int A_2 = getbbs().targetBeam().A();
0141
0142   // Do this first for the case when the first beam is the photon emitter
0143   // Treat pA separately with defined beams
0144   // The variable beam (=1,2) defines which nucleus is the target
0145   for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
0146
0147     W = _wMin + double(i)*dW + 0.5*dW;
0148
0149     for(unsigned int j = 0; j <= _nYbins - 1; ++j) {
0150
0151       Y = -1.0*_yMax + double(j)*dY + 0.5*dY;
0152
0153       if( A_2 == 1 && A_1 != 1 ){
0154         // pA, first beam is the nucleus and is in this case the target
0155         Egamma = 0.5*W*exp(-Y);
0156         beam = 1;
0157       } else if( A_1 ==1 && A_2 != 1){
0158         // pA, second beam is the nucleus and is in this case the target
0159         Egamma = 0.5*W*exp(Y);
0160         beam = 2;
0161       } else {
0162         Egamma = 0.5*W*exp(Y);
0163         beam = 2;
0164       }
0165
0166       dndWdY = 0.;
0167
0168       if( Egamma > Eth && Egamma < maxPhotonEnergy() ){
0169
0170     csgA=getcsgA(Egamma,W,beam);
0171         dndWdY = Egamma*photonFlux(Egamma,beam)*csgA*breitWigner(W,bwnorm);
0172
0173       }
0174
0175       wylumfile << dndWdY << endl;
0176
0177     }
0178   }
0179
0180   // Repeat the loop for the case when the second beam is the photon emitter.
0181   // Don't repeat for pA
0182   if( !( (A_2 == 1 && A_1 != 1) || (A_1 == 1 && A_2 != 1) ) ){
0183
0184     for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
0185
0186       W = _wMin + double(i)*dW + 0.5*dW;
0187
0188       for(unsigned int j = 0; j <= _nYbins - 1; ++j) {
0189
0190         Y = -1.0*_yMax + double(j)*dY + 0.5*dY;
0191
0192         beam = 1;
0193         Egamma = 0.5*W*exp(-Y);
0194
0195         dndWdY = 0.;
0196
0197         if( Egamma > Eth && Egamma < maxPhotonEnergy() ){
0198
0199       csgA=getcsgA(Egamma,W,beam);
0200           dndWdY = Egamma*photonFlux(Egamma,beam)*csgA*breitWigner(W,bwnorm);
0201
0202         }
0203
0204         wylumfile << dndWdY << endl;
0205
0206       }
0207     }
0208   }
0209
0210   wylumfile << bwnorm << endl;
0211   wylumfile.close();
0212
0213   if(_interferenceEnabled==1)
0214     pttablegen();
0215
0216 }
0217
0218
0219 //______________________________________________________________________________
0220 void photonNucleusLuminosity::pttablegen()
0221 {
0222   //  Calculates the pt spectra for VM production with interference
0223   //  Follows S. Klein and J. Nystrand, Phys. Rev Lett. 84, 2330 (2000).
0224   //  Written by S. Klein, 8/2002
0225
0226   //  fill in table pttable in one call
0227   //  Integrate over all y (using the same y values as in table yarray
0228   //  note that the cross section goes from ymin (<0) to ymax (>0), in numy points
0229   //  here,  we go from 0 to ymax in (numy/2)+1 points
0230   //  numy must be even.
0231
0232   //  At each y, calculate the photon energies Egamma1 and Egamma2
0233   //  and the two photon-A cross sections
0234
0235   //  loop over each p_t entry in the table.
0236
0237   //  Then, loop over b and phi (the angle between the VM \vec{p_t} and \vec{b}
0238   //  and calculate the cross section at each step.
0239   //  Put the results in pttable
0240
0241   std::string wyFileName;
0242   wyFileName = _baseFileName +".txt";
0243
0244   ofstream wylumfile;
0245   wylumfile.precision(15);
0246
0247   wylumfile.open(wyFileName.c_str(),ios::app);
0248
0249   double param1pt[500],param2pt[500];
0250   double  *ptparam1=param1pt;
0251   double  *ptparam2=param2pt;
0252   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.;
0253   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.;
0254   double sumg=0.,theta=0.,amp_i_2=0.,sumint=0.;
0255   int NGAUSS=0,NBIN=0;
0256
0257   double xg[16]={.0483076656877383162E0,.144471961582796493E0,
0258          .239287362252137075E0, .331868602282127650E0,
0259          .421351276130635345E0, .506899908932229390E0,
0260          .587715757240762329E0, .663044266930215201E0,
0261          .732182118740289680E0, .794483795967942407E0,
0262          .849367613732569970E0, .896321155766052124E0,
0263          .934906075937739689E0, .964762255587506430E0,
0264          .985611511545268335E0, .997263861849481564E0};
0265   double ag[16]={.0965400885147278006E0, .0956387200792748594E0,
0266          .0938443990808045654E0, .0911738786957638847E0,
0267          .0876520930044038111E0, .0833119242269467552E0,
0268          .0781938957870703065E0, .0723457941088485062E0,
0269          .0658222227763618468E0, .0586840934785355471E0,
0270          .0509980592623761762E0, .0428358980222266807E0,
0271          .0342738629130214331E0, .0253920653092620595E0,
0272          .0162743947309056706E0, .00701861000947009660E0};
0273
0274   NGAUSS=16;
0275
0276   //Setting input calls to variables/less calls this way.
0277   double Ymax=_yMax;
0278   int numy = _nYbins;
0279   double Ep = _protonEnergy;
0280   int ibreakup = _beamBreakupMode;
0281   double NPT = _nmbPtBinsInterference;
0282   double gamma_em = _beamLorentzGamma;
0283   double mass = getChannelMass();
0284   int beam;
0285
0286   //  loop over y from 0 (not -ymax) to yma
0287   // changed this to go from -ymax to ymax to aid asymmetric collisions
0288
0289   dY=(2.*Ymax)/numy;
0290   for(int jy=1;jy<=numy;jy++){
0291     Yp=-Ymax+((double(jy)-0.5)*dY);
0292
0293     // Find the photon energies.  Yp >= 0, so Egamma2 is smaller (no longer true if we integrate over all Y)
0294     // Use the vector meson mass for W here - neglect the width
0295
0296     Egamma1 = 0.5*mass*exp(Yp);
0297     Egamma2 = 0.5*mass*exp(-Yp);
0298
0299     //  Find the sigma(gammaA) for the two directions
0300     //  Photonuclear Cross Section 1
0301     //  Gamma-proton CM energy
0302     beam=2;
0303
0304     Wgp=sqrt(2.*Egamma1*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass*
0305                  starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
0306
0307     // Calculate V.M.+proton cross section
0308
0309     cs=sqrt(16.*starlightConstants::pi*vmPhotonCoupling()*slopeParameter()*
0310         starlightConstants::hbarc*starlightConstants::hbarc*sigmagp(Wgp)
0311         /starlightConstants::alpha);
0312     // Calculate V.M.+nucleus cross section
0313     cvma=sigma_A(cs,beam);
0314
0315     // Calculate Av = dsigma/dt(t=0) Note Units: fm**2/Gev**2
0316
0317     Av=(starlightConstants::alpha*cvma*cvma)/(16.*starlightConstants::pi
0318                           *vmPhotonCoupling()*starlightConstants::hbarc*starlightConstants::hbarc);
0319
0320     tmin  = ((mass*mass)/(4.*Egamma1*gamma_em)*(mass*mass)/(4.*Egamma1*gamma_em));
0321     tmax  = tmin + 0.25;
0322     ax    = 0.5*(tmax-tmin);
0323     bx    = 0.5*(tmax+tmin);
0324     csgA  = 0.;
0325
0326     for(int k=0;k<NGAUSS;k++){
0327       t     = sqrt(ax*xg[k]+bx);
0328       csgA  = csgA + ag[k]*getbbs().targetBeam().formFactor(t)*getbbs().targetBeam().formFactor(t);
0329       t     = sqrt(ax*(-xg[k])+bx);
0330       csgA  = csgA + ag[k]*getbbs().targetBeam().formFactor(t)*getbbs().targetBeam().formFactor(t);
0331     }
0332
0333     csgA = 0.5*(tmax-tmin)*csgA;
0334     csgA = Av*csgA;
0335     sig_ga_1 = csgA;
0336
0337     // Photonuclear Cross Section 2
0338     beam=1;
0339
0340     Wgp=sqrt(2.*Egamma2*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass*
0341                  starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
0342
0343     cs=sqrt(16.*starlightConstants::pi*vmPhotonCoupling()*slopeParameter()*
0344         starlightConstants::hbarc*starlightConstants::hbarc*sigmagp(Wgp)/starlightConstants::alpha);
0345
0346     cvma=sigma_A(cs,beam);
0347
0348     Av=(starlightConstants::alpha*cvma*cvma)/(16.*starlightConstants::pi
0349                           *vmPhotonCoupling()*starlightConstants::hbarc*starlightConstants::hbarc);
0350
0351     tmin  = (((mass*mass)/(4.*Egamma2*gamma_em))*((mass*mass)/(4.*Egamma2*gamma_em)));
0352     tmax  = tmin + 0.25;
0353     ax    = 0.5*(tmax-tmin);
0354     bx    = 0.5*(tmax+tmin);
0355     csgA  = 0.;
0356
0357     for(int k=0;k<NGAUSS;k++){
0358       t     = sqrt(ax*xg[k]+bx);
0359       csgA  = csgA + ag[k]*getbbs().electronBeam().formFactor(t)*getbbs().electronBeam().formFactor(t);
0360       t     = sqrt(ax*(-xg[k])+bx);
0361       csgA  = csgA + ag[k]*getbbs().electronBeam().formFactor(t)*getbbs().electronBeam().formFactor(t);
0362     }
0363
0364     csgA = 0.5*(tmax-tmin)*csgA;
0365     csgA = Av*csgA;
0366     sig_ga_2 = csgA;
0367
0368     //  Set up pttables - they find the reduction in sigma(pt)
0369     //  due to the nuclear form factors.
0370     //  Use the vector meson mass for W here - neglect width in
0371     //  interference calculation
0372
0373     ptparam1=vmsigmapt(mass,Egamma1,ptparam1, 2);
0374     ptparam2=vmsigmapt(mass,Egamma2,ptparam2, 1);
0375
0376     bmin = getbbs().electronBeam().nuclearRadius()+getbbs().targetBeam().nuclearRadius();
0377     //  if we allow for nuclear breakup, use a slightly smaller bmin
0378     if (ibreakup != 1)
0379       bmin=0.95*bmin;
0380
0381     //  set  bmax according to the smaller photon energy, following flux.f
0382     if (Egamma1 >=Egamma2) {
0383     bmax=bmin+6.*starlightConstants::hbarc*gamma_em/Egamma2;
0384     }
0385     else {
0386     bmax=bmin+6.*starlightConstants::hbarc*gamma_em/Egamma1;
0387     }
0388     //  set number of bins to a reasonable number to start
0389     NBIN = 2000;
0390     db   = (bmax-bmin)/float(NBIN);
0391     // loop over pt
0392     for(int i=1;i<=NPT;i++){
0393
0394       pt = (float(i)-0.5)*_ptBinWidthInterference;
0395       sum1=0.0;
0396       // loop over b
0397       for(int j=1;j<=NBIN;j++){
0398
0399     b = bmin + (float(j)-0.5)*db;
0400     A1 = Egamma1*getbbs().electronBeam().photonDensity(Egamma1,b)*sig_ga_1*ptparam1[i];
0401     A2 = Egamma2*getbbs().targetBeam().photonDensity(Egamma2,b)*sig_ga_2*ptparam2[i];
0402     sumg=0.0;
0403
0404     //  do this as a Gaussian integral, from 0 to pi
0405     for(int k=0;k<NGAUSS;k++){
0406
0407       theta=xg[k]*starlightConstants::pi;
0408       //  allow for a linear sum of interfering and non-interfering amplitudes
0409       amp_i_2 = A1 + A2 - 2.*_interferenceStrength
0410         *sqrt(A1*A2)*cos(pt*b*cos(theta)/starlightConstants::hbarc);
0411       sumg  = sumg+ag[k]*amp_i_2;
0412     }
0413     //  this is dn/dpt^2
0414     //  The factor of 2 is because the theta integral is only from 0 to pi
0415     sumint=2.*sumg*b*db;
0416     if (ibreakup > 1)
0417       sumint=sumint*getbbs().probabilityOfBreakup(b);
0418     sum1 = sum1 + sumint;
0419
0420       }
0421       //  normalization is done in readDiffLum.f
0422       //  This is d^2sigma/dpt^2; convert to dsigma/dpt
0423
0424
0425       wylumfile << sum1*pt*_ptBinWidthInterference <<endl;
0426       //  end of pt loop
0427     }
0428     //  end of y loop
0429   }
0430   wylumfile.close();
0431 }
0432
0433
0434 //______________________________________________________________________________
0435 double *photonNucleusLuminosity::vmsigmapt(double W, double Egamma, double *SIGMAPT, int beam)
0436 {
0437   //
0438   //  This subroutine calculates the effect of the nuclear form factor
0439   // on the pt spectrum, for use in interference calculations
0440   // For an interaction with mass W and photon energy Egamma,
0441   // it calculates the cross section suppression SIGMAPT(PT)
0442   // as a function of pt.
0443   // The input pt values come from pttable.inc
0444
0445
0446   double pxmax=0.,pymax=0.,dx=0.,Epom=0.,pt=0.,px0=0.,py0=0.,sum=0.,sumy=0.;
0447   double py=0.,px=0.,pt1=0.,pt2=0.,f1=0.,f2=0.,q1=0.,q2=0.,norm=0.;
0448   int NGAUSS =0,Nxbin=0;
0449   double xg[16]={.0483076656877383162e0,.144471961582796493e0,
0450          .239287362252137075e0, .331868602282127650e0,
0451          .421351276130635345e0, .506899908932229390e0,
0452          .587715757240762329e0, .663044266930215201e0,
0453          .732182118740289680e0, .794483795967942407e0,
0454          .849367613732569970e0, .896321155766052124e0,
0455          .934906075937739689e0, .964762255587506430e0,
0456          .985611511545268335e0, .997263861849481564e0};
0457   double ag[16]={.0965400885147278006e0, .0956387200792748594e0,
0458          .0938443990808045654e0, .0911738786957638847e0,
0459          .0876520930044038111e0, .0833119242269467552e0,
0460          .0781938957870703065e0, .0723457941088485062e0,
0461          .0658222227763618468e0, .0586840934785355471e0,
0462          .0509980592623761762e0, .0428358980222266807e0,
0463          .0342738629130214331e0, .0253920653092620595e0,
0464          .0162743947309056706e0, .00701861000947009660e0};
0465   NGAUSS=16;
0466
0467   //     >> Initialize
0468   if (beam == 1) {
0469   pxmax = 10.*(starlightConstants::hbarc/getbbs().targetBeam().nuclearRadius());
0470   pymax = 10.*(starlightConstants::hbarc/getbbs().targetBeam().nuclearRadius());
0471   }
0472   else {
0473   pxmax = 10.*(starlightConstants::hbarc/getbbs().electronBeam().nuclearRadius());
0474   pymax = 10.*(starlightConstants::hbarc/getbbs().electronBeam().nuclearRadius());
0475   }
0476
0477   Nxbin = 500;
0478
0479   dx = 2.*pxmax/double(Nxbin);
0480   Epom   = W*W/(4.*Egamma);
0481
0482   //     >> Loop over total Pt to find distribution
0483
0484       for(int k=1;k<=_nmbPtBinsInterference;k++){
0485
0486               pt=_ptBinWidthInterference*(double(k)-0.5);
0487
0488               px0 = pt;
0489               py0 = 0.0;
0490
0491               //  For each total Pt, integrate over Pt1, , the photon pt
0492               //  The pt of the Pomeron  is the difference
0493               //  pt1 is
0494               sum=0.;
0495               for(int i=1;i<=Nxbin;i++){
0496
0497         px = -pxmax + (double(i)-0.5)*dx;
0498         sumy=0.0;
0499         for(int j=0;j<NGAUSS;j++){
0500
0501           py = 0.5*pymax*xg[j]+0.5*pymax;
0502           //  photon pt
0503           pt1 = sqrt( px*px + py*py );
0504           //  pomeron pt
0505           pt2 = sqrt( (px-px0)*(px-px0) + (py-py0)*(py-py0) );
0506           q1  = sqrt( ((Egamma/_beamLorentzGamma)*(Egamma/_beamLorentzGamma)) + pt1*pt1 );
0507           q2  = sqrt( ((Epom/_beamLorentzGamma)*(Epom/_beamLorentzGamma))   + pt2*pt2 );
0508
0509           //  photon form factor
0510           // add in phase space factor?
0511           if (beam ==2) {
0512           f1  = (getbbs().electronBeam().formFactor(q1*q1)*getbbs().electronBeam().formFactor(q1*q1)*pt1*pt1)/(q1*q1*q1*q1);
0513
0514           //  Pomeron form factor
0515           f2  = getbbs().targetBeam().formFactor(q2*q2)*getbbs().targetBeam().formFactor(q2*q2);
0516           }
0517           else {
0518           f1  = (getbbs().targetBeam().formFactor(q1*q1)*getbbs().targetBeam().formFactor(q1*q1)*pt1*pt1)/(q1*q1*q1*q1);
0519
0520           //  Pomeron form factor
0521           f2  = getbbs().electronBeam().formFactor(q2*q2)*getbbs().electronBeam().formFactor(q2*q2);
0522           }
0523           sumy= sumy + ag[j]*f1*f2;
0524
0525           //  now consider other half of py phase space - why is this split?
0526           py = 0.5*pymax*(-xg[j])+0.5*pymax;
0527           pt1 = sqrt( px*px + py*py );
0528           pt2 = sqrt( (px-px0)*(px-px0) + (py-py0)*(py-py0) );
0529           q1  = sqrt( ((Egamma/_beamLorentzGamma)*Egamma/_beamLorentzGamma) + pt1*pt1 );
0530           q2  = sqrt( ((Epom/_beamLorentzGamma)*(Epom/_beamLorentzGamma))   + pt2*pt2 );
0531           //  add in phase space factor?
0532           if (beam ==2) {
0533             f1  = (getbbs().electronBeam().formFactor(q1*q1)*getbbs().electronBeam().formFactor(q1*q1)*pt1*pt1)/(q1*q1*q1*q1);
0534
0535           //  Pomeron form factor
0536           f2  = getbbs().targetBeam().formFactor(q2*q2)*getbbs().targetBeam().formFactor(q2*q2);
0537           }
0538           else {
0539           f1  = (getbbs().targetBeam().formFactor(q1*q1)*getbbs().targetBeam().formFactor(q1*q1)*pt1*pt1)/(q1*q1*q1*q1);
0540
0541           //  Pomeron form factor
0542           f2  = getbbs().electronBeam().formFactor(q2*q2)*getbbs().electronBeam().formFactor(q2*q2);
0543           }
0544
0545           sumy= sumy + ag[j]*f1*f2;
0546
0547         }
0548         //         >> This is to normalize the gaussian integration
0549         sumy = 0.5*pymax*sumy;
0550         //         >> The 2 is to account for py: 0 -- pymax
0551         sum  = sum + 2.*sumy*dx;
0552           }
0553
0554           if(k==1) norm=1./sum;
0555           SIGMAPT[k]=sum*norm;
0556       }
0557       return (SIGMAPT);
0558 }
0559