]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STARLIGHT/starlight/src/photonNucleusCrossSection.cpp
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
updated STARTLIGHT to r191 (http://starlight.hepforge.org/svn/trunk)
[u/mrichter/AliRoot.git] / STARLIGHT / starlight / src / photonNucleusCrossSection.cpp
CommitLineData
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:
be31281a 23// $Rev:: 189 $: revision of last commit
24// $Author:: srklein $: author of last commit
25// $Date:: 2014-10-27 04:29:14 +0100 #$: date of last commit
da32329d
AM		 26//
27// Description:
28//
29//
30///////////////////////////////////////////////////////////////////////////
31
32
33#include <iostream>
34#include <fstream>
35#include <cmath>
36
37#include "reportingUtils.h"
38#include "starlightconstants.h"
39#include "bessel.h"
40#include "photonNucleusCrossSection.h"
41
42
43using namespace std;
44using namespace starlightConstants;
45
46
47//______________________________________________________________________________
48photonNucleusCrossSection::photonNucleusCrossSection(const beamBeamSystem& bbsystem)
49 : _nWbins (inputParametersInstance.nmbWBins() ),
50 _nYbins (inputParametersInstance.nmbRapidityBins() ),
51 _wMin (inputParametersInstance.minW() ),
52 _wMax (inputParametersInstance.maxW() ),
53 _yMax (inputParametersInstance.maxRapidity() ),
54 _beamLorentzGamma (inputParametersInstance.beamLorentzGamma() ),
55 _bbs (bbsystem ),
56 _protonEnergy (inputParametersInstance.protonEnergy() ),
57 _particleType (inputParametersInstance.prodParticleType() ),
58 _beamBreakupMode (inputParametersInstance.beamBreakupMode() ),
59 _coherentProduction(inputParametersInstance.coherentProduction()),
60 _incoherentFactor (inputParametersInstance.incoherentFactor() ),
61 _productionMode (inputParametersInstance.productionMode() ),
62 _sigmaNucleus (_bbs.beam2().A() )
63{
64 // define luminosity for various beam particles in units of 10^{26} cm^{-2} sec^{-1}
65 switch(_bbs.beam1().Z()) {
66 case 1: // proton
67 _luminosity = 1.E8;
68 break;
69 case 8: // O
70 _luminosity = 980.;
71 break;
72 case 14: // Si
73 _luminosity = 440.;
74 break;
75 case 20: // Ca
76 _luminosity = 2000.;
77 break;
78 case 29: // Cu
79 _luminosity = 95.;
80 break;
81 case 49: // Indium, uses same as Iodine
82 _luminosity = 27.;
83 break;
84 case 53: // I
85 _luminosity = 27.;
86 break;
87 case 79: // Au
88 _luminosity = 2.0;
89 break;
90 case 82: // Pb
91 _luminosity = 1.;
92 break;
93 default:
94 printWarn << "luminosity is not defined for beam with Z = " << _bbs.beam1().Z()
95 << ". using " << _luminosity << " 10^{26} cm^{-2} sec^{-1}" << endl;
96 }
97
98 switch(_particleType) {
99 case RHO:
100 _slopeParameter = 11.0; // [(GeV/c)^{-2}]
101 _vmPhotonCoupling = 2.02;
102 _ANORM = -2.75;
103 _BNORM = 0.0;
104 _defaultC = 1.0;
105 _channelMass = 0.7685; // [GeV/c^2]
106 _width = 0.1507; // [GeV/c^2]
107 break;
108 case RHOZEUS:
109 _slopeParameter =11.0;
110 _vmPhotonCoupling=2.02;
111 _ANORM=-2.75;
112 _BNORM=1.84;
113 _defaultC=1.0;
114 _channelMass = 0.7685;
115 _width=0.1507;
116 break;
117 case FOURPRONG:
118 _slopeParameter = 11.0;
119 _vmPhotonCoupling = 2.02;
120 _ANORM = -2.75;
121 _BNORM = 0; // no coherent background component is implemented for four-prong
122 _defaultC = 11.0;
123 _channelMass = 1.350;
124 _width = 0.360;
125 break;
126 case OMEGA:
127 _slopeParameter=10.0;
128 _vmPhotonCoupling=23.13;
129 _ANORM=-2.75;
130 _BNORM=0.0;
131 _defaultC=1.0;
132 _channelMass=0.78194;
133 _width=0.00843;
134 break;
135 case PHI:
136 _slopeParameter=7.0;
137 _vmPhotonCoupling=13.71;
138 _ANORM=-2.75;
139 _BNORM=0.0;
140 _defaultC=1.0;
141 _channelMass=1.019413;
142 _width=0.00443;
143 break;
144 case JPSI:
145 case JPSI_ee:
146 case JPSI_mumu:
147 _slopeParameter=4.0;
148 _vmPhotonCoupling=10.45;
149 _ANORM=-2.75;//Artificial Breit-Wigner parameters--no direct pions
150 _BNORM=0.0;
151 _defaultC=1.0;
152 _channelMass=3.09692;//JN 3.09688
153 _width=0.000091;//JN 0.000087
154 break;
155 case JPSI2S:
156 case JPSI2S_ee:
157 case JPSI2S_mumu:
158 _slopeParameter=4.3;
159 _vmPhotonCoupling=26.39;
160 _ANORM=-2.75;//Artificial
161 _BNORM=0.0;
162 _defaultC=1.0;
163 _channelMass=3.686093;
164 _width=0.000337;
165 break;
166 case UPSILON:
167 case UPSILON_ee:
168 case UPSILON_mumu:
169 _slopeParameter=4.0;
170 _vmPhotonCoupling=125.37;
171 _ANORM=-2.75;//Artificial
172 _BNORM=0.0;
173 _defaultC=1.0;
174 _channelMass=9.46030;
175 _width=0.00005402;
176 break;
177 case UPSILON2S:
178 case UPSILON2S_ee:
179 case UPSILON2S_mumu:
180 _slopeParameter=4.0;
181 _vmPhotonCoupling=290.84;
182 _ANORM=-2.75;
183 _BNORM=0.0;
184 _defaultC=1.0;
185 _channelMass=10.02326;
186 _width=0.00003198;
187 break;
188 case UPSILON3S:
189 case UPSILON3S_ee:
190 case UPSILON3S_mumu:
191 _slopeParameter=4.0;
192 _vmPhotonCoupling=415.10;
193 _ANORM=-2.75;
194 _BNORM=0.0;
195 _defaultC=1.0;
196 _channelMass=10.3552;
197 _width=0.00002032;
198 break;
199 default:
200 cout <<"No sigma constants parameterized for pid: "<<_particleType
201 <<" GammaAcrosssection"<<endl;
202 }
203
204 double maxradius = (_bbs.beam1().nuclearRadius()<_bbs.beam2().nuclearRadius())?_bbs.beam2().nuclearRadius():_bbs.beam1().nuclearRadius();
205 _maxPhotonEnergy = 5. * _beamLorentzGamma * hbarc/maxradius;
206}
207
208
209//______________________________________________________________________________
210photonNucleusCrossSection::~photonNucleusCrossSection()
211{ }
212
213
214//______________________________________________________________________________
215void
216photonNucleusCrossSection::crossSectionCalculation(const double)
217{
218 cout << "Neither narrow/wide resonance cross-section calculation.--Derived" << endl;
219}
220
221
222//______________________________________________________________________________
223double
224photonNucleusCrossSection::getcsgA(const double Egamma,
225 const double W)
226{
227 //This function returns the cross-section for photon-nucleus interaction
228 //producing vectormesons
229
230 double Av,Wgp,cs,cvma;
231 double t,tmin,tmax;
232 double csgA,ax,bx;
233 int NGAUSS;
234
235 // DATA FOR GAUSS INTEGRATION
236 double xg[6] = {0, 0.1488743390, 0.4333953941, 0.6794095683, 0.8650633667, 0.9739065285};
237 double ag[6] = {0, 0.2955242247, 0.2692667193, 0.2190863625, 0.1494513492, 0.0666713443};
238 NGAUSS = 6;
239
240 // Find gamma-proton CM energy
241 Wgp = sqrt(2. * Egamma * (_protonEnergy
242 + sqrt(_protonEnergy * _protonEnergy - protonMass * protonMass))
243 + protonMass * protonMass);
244
245 //Used for d-A and A-A
246 tmin = (W * W / (4. * Egamma * _beamLorentzGamma)) * (W * W / (4. * Egamma * _beamLorentzGamma));
247
248 if ((_bbs.beam1().A() == 1) && (_bbs.beam2().A() == 1)) // proton-proton, no scaling needed
249 csgA = sigmagp(Wgp);
250 else if ((_bbs.beam2().Z() == 1) && (_bbs.beam2().A() == 2)) { // deuteron-A interaction
251 Av = _slopeParameter * sigmagp(Wgp);
252
253 tmax = tmin + 0.64; //0.64
254 ax = 0.5 * (tmax - tmin);
255 bx = 0.5 * (tmax + tmin);
256 csgA = 0.;
257
258 for (int k = 1; k < NGAUSS; ++k) {
259 t = ax * xg[k] + bx;
260 // We use beam2 here since the input stores the deuteron as nucleus 2
261 // and nucleus 2 is the pomeron field source
262 // Also this is the way sergey formatted the formfactor.
263 csgA = csgA + ag[k] * _bbs.beam2().formFactor(t);
264 t = ax * (-xg[k]) + bx;
265 csgA = csgA + ag[k] * _bbs.beam2().formFactor(t);
266 }
267 csgA = 0.5 * (tmax - tmin) * csgA;
268 csgA = Av * csgA;
269 } else if (!_coherentProduction &&
270 (!((_bbs.beam2().Z() == 1) && (_bbs.beam2().A() == 2)))) { // incoherent AA interactions
271 // For incoherent AA interactions, since incoherent treating it as gamma-p
272 // Calculate the differential V.M.+proton cross section
273 csgA = 1.E-4 * _incoherentFactor * _sigmaNucleus * _slopeParameter * sigmagp(Wgp);
274
275 } else { // coherent AA interactions
276 // For typical AA interactions.
277 // Calculate V.M.+proton cross section
278 cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha);
279
280 // Calculate V.M.+nucleus cross section
281 // cvma = _bbs.beam1().A()*cs;
282 cvma = sigma_A(cs);
283
284 // Calculate Av = dsigma/dt(t=0) Note Units: fm**s/Gev**2
285 Av = (alpha * cvma * cvma) / (16. * pi * _vmPhotonCoupling * hbarc * hbarc);
286
287 // Check if one or both beams are nuclei
288 int A_1 = _bbs.beam1().A();
289 int A_2 = _bbs.beam2().A();
290
291 tmax = tmin + 0.25;
292 ax = 0.5 * (tmax - tmin);
293 bx = 0.5 * (tmax + tmin);
294 csgA = 0.;
295 for (int k = 1; k < NGAUSS; ++k) {
296
297 t = ax * xg[k] + bx;
298 if( A_1 == 1 && A_2 != 1){
299 csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
300 }else if(A_2 ==1 && A_1 != 1){
301 csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
302 }else{
303 csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
304 }
305
306 t = ax * (-xg[k]) + bx;
307 if( A_1 == 1 && A_2 != 1){
308 csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
309 }else if(A_2 ==1 && A_1 != 1){
310 csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
311 }else{
312 csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
313 }
314 }
315 csgA = 0.5 * (tmax - tmin) * csgA;
316 csgA = Av * csgA;
317 }
318
319 return csgA;
320}
321
322
323//______________________________________________________________________________
324double
325photonNucleusCrossSection::photonFlux(const double Egamma)
326{
327 // This routine gives the photon flux as a function of energy Egamma
328 // It works for arbitrary nuclei and gamma; the first time it is
329 // called, it calculates a lookup table which is used on
330 // subsequent calls.
331 // It returns dN_gamma/dE (dimensions 1/E), not dI/dE
332 // energies are in GeV, in the lab frame
333 // rewritten 4/25/2001 by SRK
334
335 double lEgamma,Emin,Emax;
336 static double lnEmax, lnEmin, dlnE;
337 double stepmult,energy,rZ,rA;
338 int nbstep,nrstep,nphistep,nstep;
339 double bmin,bmax,bmult,biter,bold,integratedflux;
340 double fluxelement,deltar,riter;
341 double deltaphi,phiiter,dist;
342 static double dide[401];
343 double lnElt;
344 double rA2, rZ2;
345 double flux_r;
346 double Xvar;
347 int Ilt;
348 double RNuc=0.,RNuc2=0.,maxradius=0.;
349
350 RNuc=_bbs.beam1().nuclearRadius();
351 RNuc2=_bbs.beam2().nuclearRadius();
352 maxradius = (RNuc<RNuc2)?RNuc2:RNuc;
353 // static ->>> dide,lnEMax,lnEmin,dlnE
354 static int Icheck = 0;
355
356 //Check first to see if pp
357 if( _bbs.beam1().A()==1 && _bbs.beam2().A()==1 ){
358 int nbsteps = 400;
359 double bmin = 0.5;
360 double bmax = 5.0 + (5.0*_beamLorentzGamma*hbarc/Egamma);
361 double dlnb = (log(bmax)-log(bmin))/(1.*nbsteps);
362
363 double local_sum=0.0;
364
365 // Impact parameter loop
366 for (int i = 0; i<=nbsteps;i++){
367
368 double bnn0 = bmin*exp(i*dlnb);
369 double bnn1 = bmin*exp((i+1)*dlnb);
370 double db = bnn1-bnn0;
371
372 // double PofB0 = 1.0;
373 // if( bnn0 > 1.4 )PofB0=0.0;
374 // double PofB1 = 1.0;
375 // if( bnn1 > 1.4 )PofB1=0.0;
376
377 double ppslope = 19.0;
378 double GammaProfile = exp(-bnn0*bnn0/(2.*hbarc*hbarc*ppslope));
379 double PofB0 = 1. - (1. - GammaProfile)*(1. - GammaProfile);
380 GammaProfile = exp(-bnn1*bnn1/(2.*hbarc*hbarc*ppslope));
381 double PofB1 = 1. - (1. - GammaProfile)*(1. - GammaProfile);
382
383 double Xarg = Egamma*bnn0/(hbarc*_beamLorentzGamma);
384 double loc_nofe0 = (_bbs.beam1().Z()*_bbs.beam1().Z()*alpha)/
385 (pi*pi);
386 loc_nofe0 *= (1./(Egamma*bnn0*bnn0));
387 loc_nofe0 *= Xarg*Xarg*(bessel::dbesk1(Xarg))*(bessel::dbesk1(Xarg));
388
389 Xarg = Egamma*bnn1/(hbarc*_beamLorentzGamma);
390 double loc_nofe1 = (_bbs.beam1().Z()*_bbs.beam1().Z()*alpha)/
391 (pi*pi);
392 loc_nofe1 *= (1./(Egamma*bnn1*bnn1));
393 loc_nofe1 *= Xarg*Xarg*(bessel::dbesk1(Xarg))*(bessel::dbesk1(Xarg));
394
395 local_sum += loc_nofe0*(1. - PofB0)*bnn0*db;
396 local_sum += loc_nofe1*(1. - PofB1)*bnn1*db;
397
398 }
399 // End Impact parameter loop
400
401 // Note: 2*pi --> pi because of no factor 2 above
402 double flux_r=local_sum*pi;
403 return flux_r;
404
405 // bmin = nuclearRadius+nuclearRadius;
406 // flux_r = nepoint(Egamma,bmin);
407 // return flux_r;
408 }
409
410 // first call? - initialize - calculate photon flux
411 Icheck=Icheck+1;
412 if(Icheck > 1) goto L1000f;
413
414 rZ=double(_bbs.beam1().Z());
415 rA=double(_bbs.beam1().A());
416 rZ2=double(_bbs.beam2().Z()); //Sergey--dAu
417 rA2=double(_bbs.beam2().A()); //Sergey
418
419 // Nuclear breakup is done by PofB
420 // collect number of integration steps here, in one place
421
422 nbstep=1200;
423 nrstep=60;
424 nphistep=40;
425
426 // this last one is the number of energy steps
427 nstep=100;
428
429 // following previous choices, take Emin=10 keV at LHC, Emin = 1 MeV at RHIC
430 Emin=1.E-5;
431 if (_beamLorentzGamma < 500)
432 Emin=1.E-3;
433
434 // maximum energy is 6 times the cutoff
435 // Emax=12.*hbarc*_beamLorentzGamma/RNuc;
436 Emax=6.*hbarc*_beamLorentzGamma/maxradius;
437
438 // >> lnEmin <-> ln(Egamma) for the 0th bin
439 // >> lnEmax <-> ln(Egamma) for the last bin
440
441 lnEmin=log(Emin);
442 lnEmax=log(Emax);
443 dlnE=(lnEmax-lnEmin)/nstep;
444
445 cout<<" Calculating flux for photon energies from E= "<<Emin
446 <<" to "<<Emax<<" GeV (CM frame) "<<endl;
447
448
449 stepmult= exp(log(Emax/Emin)/double(nstep));
450 energy=Emin;
451
452 for (int j = 1; j<=nstep;j++){
453 energy=energy*stepmult;
454
455 // integrate flux over 2R_A < b < 2R_A+ 6* gamma hbar/energy
456 // use exponential steps
457
45d54d9a 458 bmin=0.8*(RNuc+RNuc2); //Start slightly below 2*Radius
da32329d
AM		 459 bmax=bmin + 6.*hbarc*_beamLorentzGamma/energy;
460
461 bmult=exp(log(bmax/bmin)/double(nbstep));
462 biter=bmin;
463 integratedflux=0.;
464
465 if (_bbs.beam2().Z()==1&&_bbs.beam1().A()==2){
466 //This is for deuteron-gold
467 Xvar = (RNuc+RNuc2)*energy/(hbarc*(_beamLorentzGamma));
468
469 fluxelement = (2.0/pi)*rZ*rZ*alpha/
470 energy*(Xvar*bessel::dbesk0(Xvar)*bessel::dbesk1(Xvar)-(1/2)*Xvar*Xvar*
471 (bessel::dbesk1(Xvar)*bessel::dbesk1(Xvar)-bessel::dbesk0(Xvar)*bessel::dbesk0(Xvar)));
472
473 integratedflux=integratedflux+fluxelement;
474 }else if( (_bbs.beam1().A() == 1 && _bbs.beam2().A() != 1) || (_bbs.beam2().A() == 1 && _bbs.beam1().A() != 1) ){
475 // This is pA
476 if( _productionMode == PHOTONPOMERONINCOHERENT ){
477 // This is pA incoherent
da32329d 478 double zproj = 0.0;
da32329d
AM		 479 if( _bbs.beam1().A() == 1 ){
480 zproj = (_bbs.beam2().Z());
45d54d9a 481 }
482 else if( _bbs.beam2().A() == 1 ){
da32329d 483 zproj = (_bbs.beam1().Z());
da32329d 484 }
45d54d9a 485
486 int nbsteps = 400;
487 double bmin = 0.7*(RNuc+RNuc2);
488 double bmax = 2.0*(RNuc+RNuc2) + (8.0*_beamLorentzGamma*hbarc/energy);
489 double dlnb = (log(bmax)-log(bmin))/(1.*nbsteps);
490
491 double local_sum=0.0;
492
493 // Impact parameter loop
494 for (int i = 0; i<=nbsteps; i++){
495
496 double bnn0 = bmin*exp(i*dlnb);
497 double bnn1 = bmin*exp((i+1)*dlnb);
498 double db = bnn1-bnn0;
499
500 // double PofB0 = 1.0;
501 // if( bnn0 > _bbs.beam2().nuclearRadius() + 0.7 )PofB0=0.0;
502 // double PofB1 = 1.0;
503 // if( bnn1 > _bbs.beam2().nuclearRadius() + 0.7 )PofB1=0.0;
504
505 double PofB0 = _bbs.probabilityOfBreakup(bnn0);
506 double PofB1 = _bbs.probabilityOfBreakup(bnn1);
507
508 double Xarg = energy*bnn0/(hbarc*_beamLorentzGamma);
509 double loc_nofe0 = (zproj*zproj*alpha)/
510 (pi*pi);
511 loc_nofe0 *= (1./(energy*bnn0*bnn0));
512 loc_nofe0 *= Xarg*Xarg*(bessel::dbesk1(Xarg))*(bessel::dbesk1(Xarg));
513
514 Xarg = energy*bnn1/(hbarc*_beamLorentzGamma);
515 double loc_nofe1 = (zproj*zproj*alpha)/
516 (pi*pi);
517 loc_nofe1 *= (1./(energy*bnn1*bnn1));
518 loc_nofe1 *= Xarg*Xarg*(bessel::dbesk1(Xarg))*(bessel::dbesk1(Xarg));
519 // cout<<" i: "<<i<<" bnn0: "<<bnn0<<" PofB0: "<<PofB0<<" loc_nofe0: "<<loc_nofe0<<endl;
520
521 local_sum += loc_nofe0*PofB0*bnn0*db;
522 local_sum += loc_nofe1*PofB1*bnn1*db;
523 } // End Impact parameter loop
524
525 // Note: 2*pi --> pi because of no factor 2 above
526 integratedflux = local_sum*pi;
da32329d
AM		 527 } else if ( _productionMode == PHOTONPOMERONNARROW || _productionMode == PHOTONPOMERONWIDE ){
528 // cout<<" This is pA coherent "<<" j= "<<j<<endl;
529 double localbmin = 0.0;
530 if( _bbs.beam1().A() == 1 ){
531 localbmin = _bbs.beam2().nuclearRadius() + 0.7;
532 }
533 if( _bbs.beam2().A() == 1 ){
534 localbmin = _bbs.beam1().nuclearRadius() + 0.7;
535 }
536 integratedflux = nepoint(energy,localbmin);
537 }
538 }else{
539 for (int jb = 1; jb<=nbstep;jb++){
540 bold=biter;
541 biter=biter*bmult;
542 // When we get to b>20R_A change methods - just take the photon flux
543 // at the center of the nucleus.
544 if (biter > (10.*RNuc))
545 {
546 // if there is no nuclear breakup or only hadronic breakup, which only
547 // occurs at smaller b, we can analytically integrate the flux from b~20R_A
548 // to infinity, following Jackson (2nd edition), Eq. 15.54
549 Xvar=energy*biter/(hbarc*_beamLorentzGamma);
550 // Here, there is nuclear breakup. So, we can't use the integrated flux
551 // However, we can do a single flux calculation, at the center of the nucleus
552 // Eq. 41 of Vidovic, Greiner and Soff, Phys.Rev.C47,2308(1993), among other places
553 // this is the flux per unit area
554 fluxelement = (rZ*rZ*alpha*energy)*
555 (bessel::dbesk1(Xvar))*(bessel::dbesk1(Xvar))/
556 ((pi*_beamLorentzGamma*hbarc)*
557 (pi*_beamLorentzGamma*hbarc));
558
559 }//if biter>10
560 else{
561 // integrate over nuclear surface. n.b. this assumes total shadowing -
562 // treat photons hitting the nucleus the same no matter where they strike
563 fluxelement=0.;
564 deltar=RNuc/double(nrstep);
565 riter=-deltar/2.;
566
567 for (int jr =1; jr<=nrstep;jr++){
568 riter=riter+deltar;
569 // use symmetry; only integrate from 0 to pi (half circle)
570 deltaphi=pi/double(nphistep);
571 phiiter=0.;
572
573 for( int jphi=1;jphi<= nphistep;jphi++){
574 phiiter=(double(jphi)-0.5)*deltaphi;
575 // dist is the distance from the center of the emitting nucleus to the point in question
576 dist=sqrt((biter+riter*cos(phiiter))*(biter+riter*
577 cos(phiiter))+(riter*sin(phiiter))*(riter*sin(phiiter)));
578 Xvar=energy*dist/(hbarc*_beamLorentzGamma);
579 flux_r = (rZ*rZ*alpha*energy)*
580 (bessel::dbesk1(Xvar)*bessel::dbesk1(Xvar))/
581 ((pi*_beamLorentzGamma*hbarc)*
582 (pi*_beamLorentzGamma*hbarc));
583
584 // The surface element is 2.* delta phi* r * delta r
585 // The '2' is because the phi integral only goes from 0 to pi
586 fluxelement=fluxelement+flux_r*2.*deltaphi*riter*deltar;
587 // end phi and r integrations
588 }//for(jphi)
589 }//for(jr)
590 // average fluxelement over the nuclear surface
591 fluxelement=fluxelement/(pi*RNuc*RNuc);
592 }//else
593 // multiply by volume element to get total flux in the volume element
594 fluxelement=fluxelement*2.*pi*biter*(biter-bold);
595 // modulate by the probability of nuclear breakup as f(biter)
596 if (_beamBreakupMode > 1){
597 fluxelement=fluxelement*_bbs.probabilityOfBreakup(biter);
598 }
599 integratedflux=integratedflux+fluxelement;
600
601 } //end loop over impact parameter
602 } //end of else (pp, pA, AA)
603
604 // In lookup table, store k*dN/dk because it changes less
605 // so the interpolation should be better
606 dide[j]=integratedflux*energy;
607
608 }//end loop over photon energy
609
610 // for 2nd and subsequent calls, use lookup table immediately
611
612 L1000f:
613
614 lEgamma=log(Egamma);
615 if (lEgamma < (lnEmin+dlnE) || lEgamma > lnEmax){
616 flux_r=0.0;
617 cout<<" ERROR: Egamma outside defined range. Egamma= "<<Egamma
618 <<" "<<lnEmax<<" "<<(lnEmin+dlnE)<<endl;
619 }
620 else{
621 // >> Egamma between Ilt and Ilt+1
622 Ilt = int((lEgamma-lnEmin)/dlnE);
623 // >> ln(Egamma) for first point
624 lnElt = lnEmin + Ilt*dlnE;
625 // >> Interpolate
626 flux_r = dide[Ilt] + ((lEgamma-lnElt)/dlnE)*(dide[Ilt+1]- dide[Ilt]);
627 flux_r = flux_r/Egamma;
628 }
629
630 return flux_r;
631}
632
633
634//______________________________________________________________________________
635double
636photonNucleusCrossSection::nepoint(const double Egamma,
637 const double bmin)
638{
639 // Function for the spectrum of virtual photons,
640 // dn/dEgamma, for a point charge q=Ze sweeping
641 // past the origin with velocity gamma
642 // (=1/SQRT(1-(V/c)**2)) integrated over impact
643 // parameter from bmin to infinity
644 // See Jackson eq15.54 Classical Electrodynamics
645 // Declare Local Variables
646 double beta,X,C1,bracket,nepoint_r;
647
648 beta = sqrt(1.-(1./(_beamLorentzGamma*_beamLorentzGamma)));
649 X = (bmin*Egamma)/(beta*_beamLorentzGamma*hbarc);
650
651 bracket = -0.5*beta*beta*X*X*(bessel::dbesk1(X)*bessel::dbesk1(X)
652 -bessel::dbesk0(X)*bessel::dbesk0(X));
653
654 bracket = bracket+X*bessel::dbesk0(X)*bessel::dbesk1(X);
655
656 // C1=(2.*double((_bbs.beam1().Z())*(_bbs.beam1().Z()))*
657 // alpha)/pi;
658
659 // Note: NO Z*Z!!
660 C1=(2.*alpha)/pi;
661
662 nepoint_r = C1*(1./beta)*(1./beta)*(1./Egamma)*bracket;
663
664 return nepoint_r;
665
666}
667
668
669//______________________________________________________________________________
670double
671photonNucleusCrossSection::sigmagp(const double Wgp)
672{
673 // >> Function for the gamma-proton --> VectorMeson
674 // >> cross section. Wgp is the gamma-proton CM energy.
675 // >> Unit for cross section: fm**2
676
677 double sigmagp_r=0.;
678
679 switch(_particleType)
680 {
681 case RHO:
682 case RHOZEUS:
683 case FOURPRONG:
684 sigmagp_r=1.E-4*(5.0*exp(0.22*log(Wgp))+26.0*exp(-1.23*log(Wgp)));
685 break;
686 case OMEGA:
687 sigmagp_r=1.E-4*(0.55*exp(0.22*log(Wgp))+18.0*exp(-1.92*log(Wgp)));
688 break;
689 case PHI:
690 sigmagp_r=1.E-4*0.34*exp(0.22*log(Wgp));
691 break;
692 case JPSI:
693 case JPSI_ee:
694 case JPSI_mumu:
695 sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
696 sigmagp_r*=sigmagp_r;
697 sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp));
698 // sigmagp_r=1.E-4*0.0015*exp(0.80*log(Wgp));
699 break;
700 case JPSI2S:
701 case JPSI2S_ee:
702 case JPSI2S_mumu:
703 sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp));
704 sigmagp_r*=sigmagp_r;
705 sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp));
706 sigmagp_r*=0.166;
707 // sigmagp_r=0.166*(1.E-4*0.0015*exp(0.80*log(Wgp)));
708 break;
709 case UPSILON:
710 case UPSILON_ee:
711 case UPSILON_mumu:
712 // >> This is W**1.7 dependence from QCD calculations
713 sigmagp_r=1.E-10*(0.060)*exp(1.70*log(Wgp));
714 break;
715 case UPSILON2S:
716 case UPSILON2S_ee:
717 case UPSILON2S_mumu:
718 sigmagp_r=1.E-10*(0.0259)*exp(1.70*log(Wgp));
719 break;
720 case UPSILON3S:
721 case UPSILON3S_ee:
722 case UPSILON3S_mumu:
723 sigmagp_r=1.E-10*(0.0181)*exp(1.70*log(Wgp));
724 break;
725 default: cout<< "!!! ERROR: Unidentified Vector Meson: "<< _particleType <<endl;
726 }
727 return sigmagp_r;
728}
729
730
731//______________________________________________________________________________
732double
733photonNucleusCrossSection::sigma_A(const double sig_N)
734{
735 // Nuclear Cross Section
736 // sig_N,sigma_A in (fm**2)
737
738 double sum;
739 double b,bmax,Pint,arg,sigma_A_r;
740
741 int NGAUSS;
742
743 double xg[17]=
744 {.0,
745 .0483076656877383162,.144471961582796493,
746 .239287362252137075, .331868602282127650,
747 .421351276130635345, .506899908932229390,
748 .587715757240762329, .663044266930215201,
749 .732182118740289680, .794483795967942407,
750 .849367613732569970, .896321155766052124,
751 .934906075937739689, .964762255587506430,
752 .985611511545268335, .997263861849481564
753 };
754
755 double ag[17]=
756 {.0,
757 .0965400885147278006, .0956387200792748594,
758 .0938443990808045654, .0911738786957638847,
759 .0876520930044038111, .0833119242269467552,
760 .0781938957870703065, .0723457941088485062,
761 .0658222227763618468, .0586840934785355471,
762 .0509980592623761762, .0428358980222266807,
763 .0342738629130214331, .0253920653092620595,
764 .0162743947309056706, .00701861000947009660
765 };
766
767 NGAUSS=16;
768
769 // Check if one or both beams are nuclei
770 int A_1 = _bbs.beam1().A();
771 int A_2 = _bbs.beam2().A();
772 if( A_1 == 1 && A_2 == 1)cout<<" This is pp, you should not be here..."<<endl;
773
774 // CALCULATE P(int) FOR b=0.0 - bmax (fm)
775 bmax = 25.0;
776 sum = 0.;
777 for(int IB=1;IB<=NGAUSS;IB++){
778
779 b = 0.5*bmax*xg[IB]+0.5*bmax;
780
781 if( A_1 == 1 && A_2 != 1){
782 arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b);
783 }else if(A_2 == 1 && A_1 != 1){
784 arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
785 }else{
786 arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
787 }
788
789 Pint=1.0-exp(arg);
790 sum=sum+2.*pi*b*Pint*ag[IB];
791
792
793 b = 0.5*bmax*(-xg[IB])+0.5*bmax;
794
795 if( A_1 == 1 && A_2 != 1){
796 arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b);
797 }else if(A_2 == 1 && A_1 != 1){
798 arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
799 }else{
800 arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
801 }
802
803 Pint=1.0-exp(arg);
804 sum=sum+2.*pi*b*Pint*ag[IB];
805
806 }
807
808 sum=0.5*bmax*sum;
809
810 sigma_A_r=sum;
811
812 return sigma_A_r;
813}
814
815//______________________________________________________________________________
816double
817photonNucleusCrossSection::sigma_N(const double Wgp)
818{
819 // Nucleon Cross Section in (fm**2)
820 double cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha);
821 return cs;
822}
823
824
825//______________________________________________________________________________
826double
827photonNucleusCrossSection::breitWigner(const double W,
828 const double C)
829{
830 // use simple fixed-width s-wave Breit-Wigner without coherent backgorund for rho'
831 // (PDG '08 eq. 38.56)
832 if(_particleType==FOURPRONG) {
833 if (W < 4.01 * pionChargedMass)
834 return 0;
835 const double termA = _channelMass * _width;
836 const double termA2 = termA * termA;
837 const double termB = W * W - _channelMass * _channelMass;
838 return C * _ANORM * _ANORM * termA2 / (termB * termB + termA2);
839 }
840
841 // Relativistic Breit-Wigner according to J.D. Jackson,
842 // Nuovo Cimento 34, 6692 (1964), with nonresonant term. A is the strength
843 // of the resonant term and b the strength of the non-resonant
844 // term. C is an overall normalization.
845
846 double ppi=0.,ppi0=0.,GammaPrim,rat;
847 double aa,bb,cc;
848
849 double nrbw_r;
850
851 // width depends on energy - Jackson Eq. A.2
852 // if below threshold, then return 0. Added 5/3/2001 SRK
853 // 0.5% extra added for safety margin
be31281a 854 // omega added here 10/26/2014 SRK
855 if( _particleType==RHO ||_particleType==RHOZEUS || _particleType==OMEGA){
da32329d
AM		 856 if (W < 2.01*pionChargedMass){
857 nrbw_r=0.;
858 return nrbw_r;
859 }
860 ppi=sqrt( ((W/2.)*(W/2.)) - pionChargedMass * pionChargedMass);
861 ppi0=0.358;
862 }
863
864 // handle phi-->K+K- properly
865 if (_particleType == PHI){
866 if (W < 2.*kaonChargedMass){
867 nrbw_r=0.;
868 return nrbw_r;
869 }
870 ppi=sqrt( ((W/2.)*(W/2.))- kaonChargedMass*kaonChargedMass);
871 ppi0=sqrt( ((_channelMass/2.)*(_channelMass/2.))-kaonChargedMass*kaonChargedMass);
872 }
873
874 //handle J/Psi-->e+e- properly
875 if (_particleType==JPSI || _particleType==JPSI2S){
876 if(W<2.*mel){
877 nrbw_r=0.;
878 return nrbw_r;
879 }
880 ppi=sqrt(((W/2.)*(W/2.))-mel*mel);
881 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel);
882 }
883 if (_particleType==JPSI_ee){
884 if(W<2.*mel){
885 nrbw_r=0.;
886 return nrbw_r;
887 }
888 ppi=sqrt(((W/2.)*(W/2.))-mel*mel);
889 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel);
890 }
891 if (_particleType==JPSI_mumu){
892 if(W<2.*muonMass){
893 nrbw_r=0.;
894 return nrbw_r;
895 }
896 ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass);
897 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass);
898 }
899 if (_particleType==JPSI2S_ee){
900 if(W<2.*mel){
901 nrbw_r=0.;
902 return nrbw_r;
903 }
904 ppi=sqrt(((W/2.)*(W/2.))-mel*mel);
905 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel);
906 }
907 if (_particleType==JPSI2S_mumu){
908 if(W<2.*muonMass){
909 nrbw_r=0.;
910 return nrbw_r;
911 }
912 ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass);
913 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass);
914 }
915
916 if(_particleType==UPSILON || _particleType==UPSILON2S ||_particleType==UPSILON3S ){
917 if (W<2.*muonMass){
918 nrbw_r=0.;
919 return nrbw_r;
920 }
921 ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass);
922 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass);
923 }
924
925 if(_particleType==UPSILON_mumu || _particleType==UPSILON2S_mumu ||_particleType==UPSILON3S_mumu ){
926 if (W<2.*muonMass){
927 nrbw_r=0.;
928 return nrbw_r;
929 }
930 ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass);
931 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass);
932 }
933
934 if(_particleType==UPSILON_ee || _particleType==UPSILON2S_ee ||_particleType==UPSILON3S_ee ){
935 if (W<2.*mel){
936 nrbw_r=0.;
937 return nrbw_r;
938 }
939 ppi=sqrt(((W/2.)*(W/2.))-mel*mel);
940 ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel);
941 }
942
943 if(ppi==0.&&ppi0==0.)
944 cout<<"Improper Gammaacrosssection::breitwigner, ppi&ppi0=0."<<endl;
945
946 rat=ppi/ppi0;
947 GammaPrim=_width*(_channelMass/W)*rat*rat*rat;
948
949 aa=_ANORM*sqrt(GammaPrim*_channelMass*W);
950 bb=W*W-_channelMass*_channelMass;
951 cc=_channelMass*GammaPrim;
952
953 // First real part squared
954 nrbw_r = (( (aa*bb)/(bb*bb+cc*cc) + _BNORM)*( (aa*bb)/(bb*bb+cc*cc) + _BNORM));
955
956 // Then imaginary part squared
957 nrbw_r = nrbw_r + (( (aa*cc)/(bb*bb+cc*cc) )*( (aa*cc)/(bb*bb+cc*cc) ));
958
959 // Alternative, a simple, no-background BW, following J. Breitweg et al.
960 // Eq. 15 of Eur. Phys. J. C2, 247 (1998). SRK 11/10/2000
961 // nrbw_r = (_ANORM*_mass*GammaPrim/(bb*bb+cc*cc))**2
962
963 nrbw_r = C*nrbw_r;
964
965 return nrbw_r;
966}