1 ///////////////////////////////////////////////////////////////////////////
5 // This file is part of starlight.
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.
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.
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/>.
20 ///////////////////////////////////////////////////////////////////////////
22 // File and Version Information:
23 // $Rev:: $: revision of last commit
24 // $Author:: $: author of last commit
25 // $Date:: $: date of last commit
31 ///////////////////////////////////////////////////////////////////////////
38 #include "inputParameters.h"
39 #include "beambeamsystem.h"
41 #include "starlightconstants.h"
44 #include "gammaaluminosity.h"
50 //______________________________________________________________________________
51 photonNucleusLuminosity::photonNucleusLuminosity(beamBeamSystem& bbsystem)
52 : photonNucleusCrossSection(bbsystem)
53 ,_nPtBinsInterference(inputParametersInstance.nmbPtBinsInterference())
54 ,_ptBinWidthInterference(inputParametersInstance.ptBinWidthInterference())
55 ,_interferenceStrength(inputParametersInstance.interferenceStrength())
57 cout <<"Creating Luminosity Tables."<<endl;
58 photonNucleusDifferentialLuminosity();
59 cout <<"Luminosity Tables created."<<endl;
63 //______________________________________________________________________________
64 photonNucleusLuminosity::~photonNucleusLuminosity()
68 //______________________________________________________________________________
69 void photonNucleusLuminosity::photonNucleusDifferentialLuminosity()
71 // double Av,Wgp,cs,cvma;
74 // double t,tmin,tmax;
75 double testint,dndWdY;
81 wylumfile.precision(15);
85 dW = (_wMax-_wMin)/_nWbins;
86 dY = (_yMax-(-1.0)*_yMax)/_nYbins;
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;
111 // Normalize the Breit-Wigner Distribution and write values of W to slight.txt
113 //Grabbing default value for C in the breit-wigner calculation
115 for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
116 W = _wMin + double(i)*dW + 0.5*dW;
117 testint = testint + breitWigner(W,C)*dW;
118 wylumfile << W << endl;
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;
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)));
133 for(unsigned int i = 0; i <= _nWbins - 1; ++i) {
135 W = _wMin + double(i)*dW + 0.5*dW;
137 for(unsigned int j = 0; j <= _nYbins - 1; ++j) {
139 Y = -1.0*_yMax + double(j)*dY + 0.5*dY;
141 int A_1 = getbbs().beam1().A();
142 int A_2 = getbbs().beam2().A();
143 if( A_2 == 1 && A_1 != 1 ){
144 // pA, first beam is the nucleus
145 Egamma = 0.5*W*exp(-Y);
146 } else if( A_1 ==1 && A_2 != 1){
147 // pA, second beam is the nucleus
148 Egamma = 0.5*W*exp(Y);
150 Egamma = 0.5*W*exp(Y);
155 if( Egamma > Eth && Egamma < maxPhotonEnergy() ){
157 csgA=getcsgA(Egamma,W);
158 dndWdY = Egamma*photonFlux(Egamma)*csgA*breitWigner(W,bwnorm);
162 wylumfile << dndWdY << endl;
169 if(inputParametersInstance.interferenceEnabled()==1)
172 wylumfile.open("slight.txt",ios::app);
173 cout << "bwnorm: "<< bwnorm <<endl;
174 wylumfile << bwnorm << endl;
175 wylumfile << inputParametersInstance.parameterValueKey() << endl;
180 //______________________________________________________________________________
181 void photonNucleusLuminosity::pttablegen()
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
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.
193 // At each y, calculate the photon energies Egamma1 and Egamma2
194 // and the two photon-A cross sections
196 // loop over each p_t entry in the table.
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
203 wylumfile.precision(15);
204 wylumfile.open("slight.txt",ios::app);
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;
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};
234 //Setting input calls to variables/less calls this way.
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();
243 // loop over y from 0 (not -ymax) to yma
246 for(int jy=1;jy<=numy/2;jy++){
247 Yp=(double(jy)-0.5)*dY;
249 // Find the photon energies. Yp >= 0, so Egamma2 is smaller
250 // Use the vector meson mass for W here - neglect the width
252 Egamma1 = 0.5*mass*exp(Yp);
253 Egamma2 = 0.5*mass*exp(-Yp);
255 // Find the sigma(gammaA) for the two directions
256 // Photonuclear Cross Section 1
257 // Gamma-proton CM energy
259 Wgp=sqrt(2.*Egamma1*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass*
260 starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
262 // Calculate V.M.+proton cross section
264 cs=sqrt(16.*starlightConstants::pi*vmPhotonCoupling()*slopeParameter()*
265 starlightConstants::hbarc*starlightConstants::hbarc*sigmagp(Wgp)
266 /starlightConstants::alpha);
268 // Calculate V.M.+nucleus cross section
272 // Calculate Av = dsigma/dt(t=0) Note Units: fm**2/Gev**2
274 Av=(starlightConstants::alpha*cvma*cvma)/(16.*starlightConstants::pi
275 *vmPhotonCoupling()*starlightConstants::hbarc*starlightConstants::hbarc);
277 tmin = ((mass*mass)/(4.*Egamma1*gamma_em)*(mass*mass)/(4.*Egamma1*gamma_em));
279 ax = 0.5*(tmax-tmin);
280 bx = 0.5*(tmax+tmin);
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);
290 csgA = 0.5*(tmax-tmin)*csgA;
294 // Photonuclear Cross Section 2
296 Wgp=sqrt(2.*Egamma2*(Ep+sqrt(Ep*Ep-starlightConstants::protonMass*
297 starlightConstants::protonMass))+starlightConstants::protonMass*starlightConstants::protonMass);
299 cs=sqrt(16.*starlightConstants::pi*vmPhotonCoupling()*slopeParameter()*
300 starlightConstants::hbarc*starlightConstants::hbarc*sigmagp(Wgp)/starlightConstants::alpha);
304 Av=(starlightConstants::alpha*cvma*cvma)/(16.*starlightConstants::pi
305 *vmPhotonCoupling()*starlightConstants::hbarc*starlightConstants::hbarc);
307 tmin = (((mass*mass)/(4.*Egamma2*gamma_em))*((mass*mass)/(4.*Egamma2*gamma_em)));
309 ax = 0.5*(tmax-tmin);
310 bx = 0.5*(tmax+tmin);
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);
320 csgA = 0.5*(tmax-tmin)*csgA;
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
329 ptparam1=vmsigmapt(mass,Egamma1,ptparam1);
330 ptparam2=vmsigmapt(mass,Egamma2,ptparam2);
332 // set bmax according to the smaller photon energy, following flux.f
334 bmax=bmin+6.*starlightConstants::hbarc*gamma_em/Egamma2;
335 bmin = getbbs().beam1().nuclearRadius()+getbbs().beam2().nuclearRadius();
336 // if we allow for nuclear breakup, use a slightly smaller bmin
340 // set number of bins to a reasonable number to start
343 db = (bmax-bmin)/float(NBIN);
345 for(int i=1;i<=NPT;i++){
347 pt = (float(i)-0.5)*_ptBinWidthInterference;
350 for(int j=1;j<=NBIN;j++){
352 b = bmin + (float(j)-0.5)*db;
353 // nofe is the photon flux function
354 A1 = Egamma1*nofe(Egamma1,b)*sig_ga_1*ptparam1[i];
355 A2 = Egamma2*nofe(Egamma2,b)*sig_ga_2*ptparam2[i];
357 // do this as a Gaussian integral, from 0 to pi
358 for(int k=0;k<NGAUSS;k++){
360 theta=xg[k]*starlightConstants::pi;
361 // allow for a linear sum of interfering and non-interfering amplitudes
362 amp_i_2 = A1 + A2 - 2.*_interferenceStrength
363 *sqrt(A1*A2)*cos(pt*b*cos(theta)/starlightConstants::hbarc);
364 sumg = sumg+ag[k]*amp_i_2;
367 // The factor of 2 is because the theta integral is only from 0 to pi
370 sumint=sumint*getbbs().probabilityOfBreakup(b);
371 sum1 = sum1 + sumint;
373 // normalization is done in readDiffLum.f
374 // This is d^2sigma/dpt^2; convert to dsigma/dpt
375 //CHECK THIS OUT----> write(20,*)sum1*pt*dpt
376 wylumfile << sum1*pt*_ptBinWidthInterference <<endl;
385 //______________________________________________________________________________
386 double *photonNucleusLuminosity::vmsigmapt(double W, double Egamma, double *SIGMAPT)
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
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};
419 pxmax = 10.*(starlightConstants::hbarc/getbbs().beam1().nuclearRadius());
420 pymax = 10.*(starlightConstants::hbarc/getbbs().beam1().nuclearRadius());
424 dx = 2.*pxmax/double(Nxbin);
425 Epom = W*W/(4.*Egamma);
427 // >> Loop over total Pt to find distribution
429 for(int k=1;k<=_nPtBinsInterference;k++){
431 pt=_ptBinWidthInterference*(double(k)-0.5);
436 // For each total Pt, integrate over Pt1, , the photon pt
437 // The pt of the Pomeron is the difference
440 for(int i=1;i<=Nxbin;i++){
442 px = -pxmax + (double(i)-0.5)*dx;
444 for(int j=0;j<NGAUSS;j++){
446 py = 0.5*pymax*xg[j]+0.5*pymax;
448 pt1 = sqrt( px*px + py*py );
450 pt2 = sqrt( (px-px0)*(px-px0) + (py-py0)*(py-py0) );
451 q1 = sqrt( ((Egamma/_beamLorentzGamma)*(Egamma/_beamLorentzGamma)) + pt1*pt1 );
452 q2 = sqrt( ((Epom/_beamLorentzGamma)*(Epom/_beamLorentzGamma)) + pt2*pt2 );
454 // photon form factor
455 // add in phase space factor?
456 f1 = (getbbs().beam1().formFactor(q1*q1)*getbbs().beam1().formFactor(q1*q1)*pt1*pt1)/(q1*q1*q1*q1);
458 // Pomeron form factor
459 f2 = getbbs().beam1().formFactor(q2*q2)*getbbs().beam1().formFactor(q2*q2);
460 sumy= sumy + ag[j]*f1*f2;
462 // now consider other half of py phase space - why is this split?
463 py = 0.5*pymax*(-xg[j])+0.5*pymax;
464 pt1 = sqrt( px*px + py*py );
465 pt2 = sqrt( (px-px0)*(px-px0) + (py-py0)*(py-py0) );
466 q1 = sqrt( ((Egamma/_beamLorentzGamma)*Egamma/_beamLorentzGamma) + pt1*pt1 );
467 q2 = sqrt( ((Epom/_beamLorentzGamma)*(Epom/_beamLorentzGamma)) + pt2*pt2 );
468 // add in phase space factor?
469 f1 = (getbbs().beam1().formFactor(q1*q1)*getbbs().beam1().formFactor(q1*q1)*pt1*pt1)/(q1*q1*q1*q1);
470 f2 = getbbs().beam1().formFactor(q2*q2)*getbbs().beam1().formFactor(q2*q2);
471 sumy= sumy + ag[j]*f1*f2;
474 // >> This is to normalize the gaussian integration
475 sumy = 0.5*pymax*sumy;
476 // >> The 2 is to account for py: 0 -- pymax
477 sum = sum + 2.*sumy*dx;
480 if(k==1) norm=1./sum;
487 //______________________________________________________________________________
488 double photonNucleusLuminosity::nofe(double Egamma, double bimp)
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
494 double X=0.,nofex=0.,factor1=0.,factor2=0.,factor3=0.;
496 X = (bimp*Egamma)/(_beamLorentzGamma*starlightConstants::hbarc);
499 cout<<"In nofe, X= "<<X<<endl;
501 factor1 = (double(getbbs().beam1().Z()*getbbs().beam1().Z())
502 *starlightConstants::alpha)/(starlightConstants::pi*starlightConstants::pi);
504 factor2 = 1./(Egamma*bimp*bimp);
505 factor3 = X*X*(bessel::dbesk1(X))*(bessel::dbesk1(X));
506 nofex = factor1*factor2*factor3;