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