[u/mrichter/AliRoot.git] / STARLIGHT / starlight / src / gammaaluminosity.cpp

///////////////////////////////////////////////////////////////////////////
//
//    Copyright 2010
//
//    This file is part of starlight.
//
//    starlight is free software: you can redistribute it and/or modify
//    it under the terms of the GNU General Public License as published by
//    the Free Software Foundation, either version 3 of the License, or
//    (at your option) any later version.
//
//    starlight is distributed in the hope that it will be useful,
//    but WITHOUT ANY WARRANTY; without even the implied warranty of
//    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
//    GNU General Public License for more details.
//
//    You should have received a copy of the GNU General Public License
//    along with starlight. If not, see <http://www.gnu.org/licenses/>.
//
///////////////////////////////////////////////////////////////////////////
//
// File and Version Information:
// $Rev:: 164                         $: revision of last commit
// $Author:: odjuvsla                 $: author of last commit
// $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit
//
// Description:
//
//
//
///////////////////////////////////////////////////////////////////////////


#include <iostream>
#include <fstream>
#include <cmath>

#include "inputParameters.h"
#include "beambeamsystem.h"
#include "beam.h"
#include "starlightconstants.h"
#include "nucleus.h"
#include "bessel.h"
#include "gammaaluminosity.h"


using namespace std;


//______________________________________________________________________________
photonNucleusLuminosity::photonNucleusLuminosity(beamBeamSystem& bbsystem)
  : photonNucleusCrossSection(bbsystem)
  ,_nPtBinsInterference(inputParametersInstance.nmbPtBinsInterference())
  ,_ptBinWidthInterference(inputParametersInstance.ptBinWidthInterference())
  ,_interferenceStrength(inputParametersInstance.interferenceStrength())
{
  cout <<"Creating Luminosity Tables."<<endl;
  photonNucleusDifferentialLuminosity();
  cout <<"Luminosity Tables created."<<endl;
}


//______________________________________________________________________________
photonNucleusLuminosity::~photonNucleusLuminosity()
{ }


//______________________________________________________________________________
void photonNucleusLuminosity::photonNucleusDifferentialLuminosity()
{
	// double Av,Wgp,cs,cvma;
  double W,dW,dY;
  double Egamma,Y;
  // double t,tmin,tmax;
  double testint,dndWdY;
  double csgA;
  // double ax,bx;
  double C;  

  ofstream wylumfile;
  wylumfile.precision(15);
  
  double  bwnorm,Eth;

  dW = (_wMax-_wMin)/_nWbins;
  dY  = (_yMax-(-1.0)*_yMax)/_nYbins;
    
  // Write the values of W used in the calculation to slight.txt.  
  wylumfile.open("slight.txt");
  wylumfile << inputParametersInstance.parameterValueKey() << endl;
  wylumfile << getbbs().beam1().Z() <<endl;
  wylumfile << getbbs().beam1().A() <<endl;
  wylumfile << getbbs().beam2().Z() <<endl;
  wylumfile << getbbs().beam2().A() <<endl;
  wylumfile << inputParametersInstance.beamLorentzGamma() <<endl;
  wylumfile << inputParametersInstance.maxW() <<endl;
  wylumfile << inputParametersInstance.minW() <<endl;
  wylumfile << inputParametersInstance.nmbWBins() <<endl;
  wylumfile << inputParametersInstance.maxRapidity() <<endl;
  wylumfile << inputParametersInstance.nmbRapidityBins() <<endl;
  wylumfile << inputParametersInstance.productionMode() <<endl;
  wylumfile << inputParametersInstance.beamBreakupMode() <<endl;
  wylumfile << inputParametersInstance.interferenceEnabled() <<endl;
  wylumfile << inputParametersInstance.interferenceStrength() <<endl;
  wylumfile << inputParametersInstance.coherentProduction() <<endl;
  wylumfile << inputParametersInstance.incoherentFactor() <<endl;
  wylumfile << inputParametersInstance.deuteronSlopePar() <<endl;
  wylumfile << inputParametersInstance.maxPtInterference() <<endl;
  wylumfile << inputParametersInstance.nmbPtBinsInterference() <<endl;
  
  //     Normalize the Breit-Wigner Distribution and write values of W to slight.txt
  testint=0.0;
  //Grabbing default value for C in the breit-wigner calculation
  C=getDefaultC();
  for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
    W = _wMin + double(i)*dW + 0.5*dW;
    testint = testint + breitWigner(W,C)*dW;
    wylumfile << W << endl;
  }
  bwnorm = 1./testint;
  
  //     Write the values of Y used in the calculation to slight.txt.
  for(unsigned int i = 0; i <= _nYbins - 1; ++i) {
    Y = -1.0*_yMax + double(i)*dY + 0.5*dY;
    wylumfile << Y << endl;
  }
    
  Eth=0.5*(((_wMin+starlightConstants::protonMass)*(_wMin
							    +starlightConstants::protonMass)-starlightConstants::protonMass*starlightConstants::protonMass)/
	   (inputParametersInstance.protonEnergy()+sqrt(inputParametersInstance.protonEnergy()*
					       inputParametersInstance.protonEnergy()-starlightConstants::protonMass*starlightConstants::protonMass)));
  
  for(unsigned int i = 0; i <= _nWbins - 1; ++i) {

    W = _wMin + double(i)*dW + 0.5*dW;
    
    for(unsigned int j = 0; j <= _nYbins - 1; ++j) {

      Y = -1.0*_yMax + double(j)*dY + 0.5*dY;

      int A_1 = getbbs().beam1().A(); 
      int A_2 = getbbs().beam2().A();
      if( A_2 == 1 && A_1 != 1 ){
        // pA, first beam is the nucleus 
        Egamma = 0.5*W*exp(-Y);
      } else if( A_1 ==1 && A_2 != 1){
        // pA, second beam is the nucleus 
        Egamma = 0.5*W*exp(Y); 
      } else {
        Egamma = 0.5*W*exp(Y);        
      }

      dndWdY = 0.; 

      if( Egamma > Eth && Egamma < maxPhotonEnergy() ){

       csgA=getcsgA(Egamma,W);
       dndWdY = Egamma*photonFlux(Egamma)*csgA*breitWigner(W,bwnorm);

      }

      wylumfile << dndWdY << endl;

    }
  }

  wylumfile.close();
  
  if(inputParametersInstance.interferenceEnabled()==1) 
    pttablegen();
 
  wylumfile.open("slight.txt",ios::app);
  cout << "bwnorm: "<< bwnorm <<endl;
  wylumfile << bwnorm << endl;
  wylumfile << inputParametersInstance.parameterValueKey() << endl;
  wylumfile.close();
}


//______________________________________________________________________________
void photonNucleusLuminosity::pttablegen()
{
  //  Calculates the pt spectra for VM production with interference
  //  Follows S. Klein and J. Nystrand, Phys. Rev Lett. 84, 2330 (2000).
  //  Written by S. Klein, 8/2002
  
  //  fill in table pttable in one call
  //  Integrate over all y (using the same y values as in table yarray
  //  note that the cross section goes from ymin (<0) to ymax (>0), in numy points
  //  here,  we go from 0 to ymax in (numy/2)+1 points
  //  numy must be even.
  
  //  At each y, calculate the photon energies Egamma1 and Egamma2
  //  and the two photon-A cross sections
  
  //  loop over each p_t entry in the table.
  
  //  Then, loop over b and phi (the angle between the VM \vec{p_t} and \vec{b}
  //  and calculate the cross section at each step.
  //  Put the results in pttable
  
  ofstream wylumfile;
  wylumfile.precision(15);
  wylumfile.open("slight.txt",ios::app);
  
  
  double param1pt[500],param2pt[500];
  double  *ptparam1=param1pt;
  double  *ptparam2=param2pt;
  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.;
  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.;
  double sumg=0.,theta=0.,amp_i_2=0.,sumint=0.;
  int NGAUSS=0,NBIN=0,NThetaBIN=0;
  
  double xg[16]={.0483076656877383162E0,.144471961582796493E0,
		 .239287362252137075E0, .331868602282127650E0,
		 .421351276130635345E0, .506899908932229390E0,
		 .587715757240762329E0, .663044266930215201E0,
		 .732182118740289680E0, .794483795967942407E0,
		 .849367613732569970E0, .896321155766052124E0,
		 .934906075937739689E0, .964762255587506430E0,
		 .985611511545268335E0, .997263861849481564E0};
  double ag[16]={.0965400885147278006E0, .0956387200792748594E0,
		 .0938443990808045654E0, .0911738786957638847E0,
		 .0876520930044038111E0, .0833119242269467552E0,
		 .0781938957870703065E0, .0723457941088485062E0,
		 .0658222227763618468E0, .0586840934785355471E0,
		 .0509980592623761762E0, .0428358980222266807E0,
		 .0342738629130214331E0, .0253920653092620595E0,
		 .0162743947309056706E0, .00701861000947009660E0};

  NGAUSS=16;

  //Setting input calls to variables/less calls this way.
  double Ymax=_yMax;
  int numy = _nYbins;
  double Ep = inputParametersInstance.protonEnergy();
  int ibreakup = inputParametersInstance.beamBreakupMode();
  double NPT = inputParametersInstance.nmbPtBinsInterference();
  double gamma_em=inputParametersInstance.beamLorentzGamma();
  double mass= getChannelMass();
  
  //  loop over y from 0 (not -ymax) to yma
  
  dY=(2.*Ymax)/numy;
  for(int jy=1;jy<=numy/2;jy++){
    Yp=(double(jy)-0.5)*dY;
    
    // Find the photon energies.  Yp >= 0, so Egamma2 is smaller
    // Use the vector meson mass for W here - neglect the width
    
    Egamma1 = 0.5*mass*exp(Yp);
    Egamma2 = 0.5*mass*exp(-Yp);
    
    //  Find the sigma(gammaA) for the two directions
    //  Photonuclear Cross Section 1
    //  Gamma-proton CM energy
    
    Wgp=sqrt(2.*Egamma1*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass*
				 starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
    
    // Calculate V.M.+proton cross section
    
    cs=sqrt(16.*starlightConstants::pi*vmPhotonCoupling()*slopeParameter()*
	    starlightConstants::hbarc*starlightConstants::hbarc*sigmagp(Wgp)
	    /starlightConstants::alpha);
    
    // Calculate V.M.+nucleus cross section
    
    cvma=sigma_A(cs);
    
    // Calculate Av = dsigma/dt(t=0) Note Units: fm**2/Gev**2
    
    Av=(starlightConstants::alpha*cvma*cvma)/(16.*starlightConstants::pi
					      *vmPhotonCoupling()*starlightConstants::hbarc*starlightConstants::hbarc);
    
    tmin  = ((mass*mass)/(4.*Egamma1*gamma_em)*(mass*mass)/(4.*Egamma1*gamma_em));
    tmax  = tmin + 0.25;
    ax    = 0.5*(tmax-tmin);
    bx    = 0.5*(tmax+tmin);
    csgA  = 0.;
    
    for(int k=0;k<NGAUSS;k++){
      t     = sqrt(ax*xg[k]+bx);
      csgA  = csgA + ag[k]*getbbs().beam2().formFactor(t)*getbbs().beam2().formFactor(t);
      t     = sqrt(ax*(-xg[k])+bx);
      csgA  = csgA + ag[k]*getbbs().beam2().formFactor(t)*getbbs().beam2().formFactor(t);
    }
    
    csgA = 0.5*(tmax-tmin)*csgA;
    csgA = Av*csgA;
    sig_ga_1 = csgA;
	   
    // Photonuclear Cross Section 2
    
    Wgp=sqrt(2.*Egamma2*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass*
				 starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
    
    cs=sqrt(16.*starlightConstants::pi*vmPhotonCoupling()*slopeParameter()*
	    starlightConstants::hbarc*starlightConstants::hbarc*sigmagp(Wgp)/starlightConstants::alpha);
    
    cvma=sigma_A(cs);
    
    Av=(starlightConstants::alpha*cvma*cvma)/(16.*starlightConstants::pi
					      *vmPhotonCoupling()*starlightConstants::hbarc*starlightConstants::hbarc);
    
    tmin  = (((mass*mass)/(4.*Egamma2*gamma_em))*((mass*mass)/(4.*Egamma2*gamma_em)));
    tmax  = tmin + 0.25;
    ax    = 0.5*(tmax-tmin);
    bx    = 0.5*(tmax+tmin);
    csgA  = 0.;
    
    for(int k=0;k<NGAUSS;k++){
      t     = sqrt(ax*xg[k]+bx);
      csgA  = csgA + ag[k]*getbbs().beam2().formFactor(t)*getbbs().beam2().formFactor(t);
      t     = sqrt(ax*(-xg[k])+bx);
      csgA  = csgA + ag[k]*getbbs().beam2().formFactor(t)*getbbs().beam2().formFactor(t);
    }
	   
    csgA = 0.5*(tmax-tmin)*csgA;
    csgA = Av*csgA;
    sig_ga_2 = csgA;
    
    //  Set up pttables - they find the reduction in sigma(pt)
    //  due to the nuclear form factors.
    //  Use the vector meson mass for W here - neglect width in
    //  interference calculation
    
    ptparam1=vmsigmapt(mass,Egamma1,ptparam1);
    ptparam2=vmsigmapt(mass,Egamma2,ptparam2);
    
    //  set  bmax according to the smaller photon energy, following flux.f
    
    bmax=bmin+6.*starlightConstants::hbarc*gamma_em/Egamma2;
    bmin = getbbs().beam1().nuclearRadius()+getbbs().beam2().nuclearRadius();
    //  if we allow for nuclear breakup, use a slightly smaller bmin
    
    if (ibreakup != 1) 
      bmin=0.95*bmin;
    //  set number of bins to a reasonable number to start
    NBIN = 2000;
    NThetaBIN = 1000;
    db   = (bmax-bmin)/float(NBIN);
    // loop over pt
    for(int i=1;i<=NPT;i++){
      
      pt = (float(i)-0.5)*_ptBinWidthInterference;
      sum1=0.0;
      // loop over b
      for(int j=1;j<=NBIN;j++){
	
	b = bmin + (float(j)-0.5)*db;
	//  nofe is the photon flux function
	A1 = Egamma1*nofe(Egamma1,b)*sig_ga_1*ptparam1[i];
	A2 = Egamma2*nofe(Egamma2,b)*sig_ga_2*ptparam2[i];
	sumg=0.0;
	//  do this as a Gaussian integral, from 0 to pi
	for(int k=0;k<NGAUSS;k++){
	  
	  theta=xg[k]*starlightConstants::pi;
	  //  allow for a linear sum of interfering and non-interfering amplitudes
	  amp_i_2 = A1 + A2 - 2.*_interferenceStrength
	    *sqrt(A1*A2)*cos(pt*b*cos(theta)/starlightConstants::hbarc);
	  sumg  = sumg+ag[k]*amp_i_2;
	}
	//  this is dn/dpt^2
	//  The factor of 2 is because the theta integral is only from 0 to pi
	sumint=2.*sumg*b*db;
	if (ibreakup > 1)
	  sumint=sumint*getbbs().probabilityOfBreakup(b);
	sum1 = sum1 + sumint;
      }
      //  normalization is done in readDiffLum.f
      //  This is d^2sigma/dpt^2; convert to dsigma/dpt
      //CHECK THIS OUT---->      write(20,*)sum1*pt*dpt
      wylumfile << sum1*pt*_ptBinWidthInterference <<endl;
      //  end of pt loop
    }
    //  end of y loop
  }
  wylumfile.close();
}


//______________________________________________________________________________
double *photonNucleusLuminosity::vmsigmapt(double W, double Egamma, double *SIGMAPT)
{
  //
  //  This subroutine calculates the effect of the nuclear form factor
  // on the pt spectrum, for use in interference calculations
  // For an interaction with mass W and photon energy Egamma,
  // it calculates the cross section suppression SIGMAPT(PT)
  // as a function of pt.
  // The input pt values come from pttable.inc

  
  double pxmax=0.,pymax=0.,dx=0.,Epom=0.,pt=0.,px0=0.,py0=0.,sum=0.,sumy=0.;
  double py=0.,px=0.,pt1=0.,pt2=0.,f1=0.,f2=0.,q1=0.,q2=0.,norm=0.;
  int NGAUSS =0,Nxbin=0;
  double xg[16]={.0483076656877383162e0,.144471961582796493e0,
		 .239287362252137075e0, .331868602282127650e0,
		 .421351276130635345e0, .506899908932229390e0,
		 .587715757240762329e0, .663044266930215201e0,
		 .732182118740289680e0, .794483795967942407e0,
		 .849367613732569970e0, .896321155766052124e0,
		 .934906075937739689e0, .964762255587506430e0,
		 .985611511545268335e0, .997263861849481564e0};
  double ag[16]={.0965400885147278006e0, .0956387200792748594e0,
		 .0938443990808045654e0, .0911738786957638847e0,
		 .0876520930044038111e0, .0833119242269467552e0,
		 .0781938957870703065e0, .0723457941088485062e0,
		 .0658222227763618468e0, .0586840934785355471e0,
		 .0509980592623761762e0, .0428358980222266807e0,
		 .0342738629130214331e0, .0253920653092620595e0,
		 .0162743947309056706e0, .00701861000947009660e0};
  NGAUSS=16;

  //     >> Initialize
  pxmax = 10.*(starlightConstants::hbarc/getbbs().beam1().nuclearRadius());
  pymax = 10.*(starlightConstants::hbarc/getbbs().beam1().nuclearRadius());
  
  Nxbin = 500;
  
  dx = 2.*pxmax/double(Nxbin);
  Epom   = W*W/(4.*Egamma);
  
  //     >> Loop over total Pt to find distribution
      
      for(int k=1;k<=_nPtBinsInterference;k++){
	
              pt=_ptBinWidthInterference*(double(k)-0.5);
              
              px0 = pt;
              py0 = 0.0;
              
              //  For each total Pt, integrate over Pt1, , the photon pt
              //  The pt of the Pomeron  is the difference
              //  pt1 is
              sum=0.;
              for(int i=1;i<=Nxbin;i++){
		
		px = -pxmax + (double(i)-0.5)*dx;
		sumy=0.0;
		for(int j=0;j<NGAUSS;j++){

		  py = 0.5*pymax*xg[j]+0.5*pymax;
		  //  photon pt
		  pt1 = sqrt( px*px + py*py );
		  //  pomeron pt
		  pt2 = sqrt( (px-px0)*(px-px0) + (py-py0)*(py-py0) );
		  q1  = sqrt( ((Egamma/_beamLorentzGamma)*(Egamma/_beamLorentzGamma)) + pt1*pt1 );        
		  q2  = sqrt( ((Epom/_beamLorentzGamma)*(Epom/_beamLorentzGamma))   + pt2*pt2 );
		  
		  //  photon form factor
		  // add in phase space factor?
		  f1  = (getbbs().beam1().formFactor(q1*q1)*getbbs().beam1().formFactor(q1*q1)*pt1*pt1)/(q1*q1*q1*q1);
		  
		  //  Pomeron form factor
		  f2  = getbbs().beam1().formFactor(q2*q2)*getbbs().beam1().formFactor(q2*q2);
		  sumy= sumy + ag[j]*f1*f2;
		  
		  //  now consider other half of py phase space - why is this split?
		  py = 0.5*pymax*(-xg[j])+0.5*pymax;
		  pt1 = sqrt( px*px + py*py );
		  pt2 = sqrt( (px-px0)*(px-px0) + (py-py0)*(py-py0) );
		  q1  = sqrt( ((Egamma/_beamLorentzGamma)*Egamma/_beamLorentzGamma) + pt1*pt1 );
		  q2  = sqrt( ((Epom/_beamLorentzGamma)*(Epom/_beamLorentzGamma))   + pt2*pt2 );
		  //  add in phase space factor?
		  f1  = (getbbs().beam1().formFactor(q1*q1)*getbbs().beam1().formFactor(q1*q1)*pt1*pt1)/(q1*q1*q1*q1);
		  f2  = getbbs().beam1().formFactor(q2*q2)*getbbs().beam1().formFactor(q2*q2);
		  sumy= sumy + ag[j]*f1*f2;
      
		}
		//         >> This is to normalize the gaussian integration
		sumy = 0.5*pymax*sumy;
		//         >> The 2 is to account for py: 0 -- pymax
		sum  = sum + 2.*sumy*dx;
	      }
	      
	      if(k==1) norm=1./sum;
	      SIGMAPT[k]=sum*norm;
      }
      return (SIGMAPT);
}


//______________________________________________________________________________
double photonNucleusLuminosity::nofe(double Egamma, double bimp)
{
  //Function for the calculation of the "photon density".
  //nofe=numberofgammas/(energy*area)
  //Assume beta=1.0 and gamma>>1, i.e. neglect the (1/gamma^2)*K0(x) term
  
  double X=0.,nofex=0.,factor1=0.,factor2=0.,factor3=0.;
  
  X = (bimp*Egamma)/(_beamLorentzGamma*starlightConstants::hbarc);
  
  if( X <= 0.0) 
    cout<<"In nofe, X= "<<X<<endl;
  
  factor1 = (double(getbbs().beam1().Z()*getbbs().beam1().Z())
	     *starlightConstants::alpha)/(starlightConstants::pi*starlightConstants::pi);

  factor2 = 1./(Egamma*bimp*bimp);
  factor3 = X*X*(bessel::dbesk1(X))*(bessel::dbesk1(X));
  nofex    = factor1*factor2*factor3;
  return nofex;
}
Commit	Line	Data
da32329d AM	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.5dW;
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.5dW;
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.5Wexp(-Y);
146	} else if( A_1 ==1 && A_2 != 1){
147	// pA, second beam is the nucleus
148	Egamma = 0.5Wexp(Y);
149	} else {
150	Egamma = 0.5Wexp(Y);
151	}
152
153	dndWdY = 0.;
154
155	if( Egamma > Eth && Egamma < maxPhotonEnergy() ){
156
157	csgA=getcsgA(Egamma,W);
158	dndWdY = EgammaphotonFlux(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.5massexp(Yp);
253	Egamma2 = 0.5massexp(-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(EpEp-starlightConstants::protonMass
260	starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
261
262	// Calculate V.M.+proton cross section
263
264	cs=sqrt(16.starlightConstants::pivmPhotonCoupling()slopeParameter()
265	starlightConstants::hbarcstarlightConstants::hbarcsigmagp(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: fm2/Gev2
273
274	Av=(starlightConstants::alphacvmacvma)/(16.*starlightConstants::pi
275	vmPhotonCoupling()starlightConstants::hbarc*starlightConstants::hbarc);
276
277	tmin = ((massmass)/(4.Egamma1gamma_em)(massmass)/(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(EpEp-starlightConstants::protonMass
297	starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
298
299	cs=sqrt(16.starlightConstants::pivmPhotonCoupling()slopeParameter()
300	starlightConstants::hbarcstarlightConstants::hbarcsigmagp(Wgp)/starlightConstants::alpha);
301
302	cvma=sigma_A(cs);
303
304	Av=(starlightConstants::alphacvmacvma)/(16.*starlightConstants::pi
305	vmPhotonCoupling()starlightConstants::hbarc*starlightConstants::hbarc);
306
307	tmin = (((massmass)/(4.Egamma2gamma_em))((massmass)/(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::hbarcgamma_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 = Egamma1nofe(Egamma1,b)sig_ga_1*ptparam1[i];
355	A2 = Egamma2nofe(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(A1A2)cos(ptb*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.sumgb*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,)sum1pt*dpt
376	wylumfile << sum1pt_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 = WW/(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.5pymaxxg[j]+0.5*pymax;
447	// photon pt
448	pt1 = sqrt( pxpx + pypy );
449	// pomeron pt
450	pt2 = sqrt( (px-px0)(px-px0) + (py-py0)(py-py0) );
451	q1 = sqrt( ((Egamma/_beamLorentzGamma)(Egamma/_beamLorentzGamma)) + pt1pt1 );
452	q2 = sqrt( ((Epom/_beamLorentzGamma)(Epom/_beamLorentzGamma)) + pt2pt2 );
453
454	// photon form factor
455	// add in phase space factor?
456	f1 = (getbbs().beam1().formFactor(q1q1)getbbs().beam1().formFactor(q1q1)pt1pt1)/(q1q1q1q1);
457
458	// Pomeron form factor
459	f2 = getbbs().beam1().formFactor(q2q2)getbbs().beam1().formFactor(q2*q2);
460	sumy= sumy + ag[j]f1f2;
461
462	// now consider other half of py phase space - why is this split?
463	py = 0.5pymax(-xg[j])+0.5*pymax;
464	pt1 = sqrt( pxpx + pypy );
465	pt2 = sqrt( (px-px0)(px-px0) + (py-py0)(py-py0) );
466	q1 = sqrt( ((Egamma/_beamLorentzGamma)Egamma/_beamLorentzGamma) + pt1pt1 );
467	q2 = sqrt( ((Epom/_beamLorentzGamma)(Epom/_beamLorentzGamma)) + pt2pt2 );
468	// add in phase space factor?
469	f1 = (getbbs().beam1().formFactor(q1q1)getbbs().beam1().formFactor(q1q1)pt1pt1)/(q1q1q1q1);
470	f2 = getbbs().beam1().formFactor(q2q2)getbbs().beam1().formFactor(q2*q2);
471	sumy= sumy + ag[j]f1f2;
472
473	}
474	// >> This is to normalize the gaussian integration
475	sumy = 0.5pymaxsumy;
476	// >> The 2 is to account for py: 0 -- pymax
477	sum = sum + 2.sumydx;
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 = (bimpEgamma)/(_beamLorentzGammastarlightConstants::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::pistarlightConstants::pi);
503
504	factor2 = 1./(Egammabimpbimp);
505	factor3 = XX(bessel::dbesk1(X))*(bessel::dbesk1(X));
506	nofex = factor1factor2factor3;
507	return nofex;
508	}