[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:: 176                         $: revision of last commit
// $Author:: jseger                   $: author of last commit
// $Date:: 2014-06-20 22:15:20 +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;
using namespace starlightConstants;


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