]>
Commit | Line | Data |
---|---|---|
da32329d AM |
1 | /////////////////////////////////////////////////////////////////////////// |
2 | // | |
3 | // Copyright 2010 | |
4 | // | |
5 | // This file is part of starlight. | |
6 | // | |
7 | // starlight is free software: you can redistribute it and/or modify | |
8 | // it under the terms of the GNU General Public License as published by | |
9 | // the Free Software Foundation, either version 3 of the License, or | |
10 | // (at your option) any later version. | |
11 | // | |
12 | // starlight is distributed in the hope that it will be useful, | |
13 | // but WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
15 | // GNU General Public License for more details. | |
16 | // | |
17 | // You should have received a copy of the GNU General Public License | |
18 | // along with starlight. If not, see <http://www.gnu.org/licenses/>. | |
19 | // | |
20 | /////////////////////////////////////////////////////////////////////////// | |
21 | // | |
22 | // File and Version Information: | |
23 | // $Rev:: 164 $: revision of last commit | |
24 | // $Author:: odjuvsla $: author of last commit | |
25 | // $Date:: 2013-10-06 16:18:08 +0200 #$: date of last commit | |
26 | // | |
27 | // Description: | |
28 | // | |
29 | // | |
30 | /////////////////////////////////////////////////////////////////////////// | |
31 | ||
32 | ||
33 | #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 | ||
43 | using namespace std; | |
44 | using namespace starlightConstants; | |
45 | ||
46 | ||
47 | //______________________________________________________________________________ | |
48 | photonNucleusCrossSection::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 | //______________________________________________________________________________ | |
210 | photonNucleusCrossSection::~photonNucleusCrossSection() | |
211 | { } | |
212 | ||
213 | ||
214 | //______________________________________________________________________________ | |
215 | void | |
216 | photonNucleusCrossSection::crossSectionCalculation(const double) | |
217 | { | |
218 | cout << "Neither narrow/wide resonance cross-section calculation.--Derived" << endl; | |
219 | } | |
220 | ||
221 | ||
222 | //______________________________________________________________________________ | |
223 | double | |
224 | photonNucleusCrossSection::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 | //______________________________________________________________________________ | |
324 | double | |
325 | photonNucleusCrossSection::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 | ||
458 | bmin=RNuc+RNuc2; //2.*nuclearRadius; Sergey | |
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 | |
478 | // cout<<" This is incoherent! "<<" j = "<<j<<endl; | |
479 | double zproj = 0.0; | |
480 | double localbmin = 0.0; | |
481 | if( _bbs.beam1().A() == 1 ){ | |
482 | zproj = (_bbs.beam2().Z()); | |
483 | localbmin = _bbs.beam2().nuclearRadius() + 0.7; | |
484 | } | |
485 | if( _bbs.beam2().A() == 1 ){ | |
486 | zproj = (_bbs.beam1().Z()); | |
487 | localbmin = _bbs.beam1().nuclearRadius() + 0.7; | |
488 | } | |
489 | integratedflux = zproj*zproj*nepoint(energy,localbmin); | |
490 | } else if ( _productionMode == PHOTONPOMERONNARROW || _productionMode == PHOTONPOMERONWIDE ){ | |
491 | // cout<<" This is pA coherent "<<" j= "<<j<<endl; | |
492 | double localbmin = 0.0; | |
493 | if( _bbs.beam1().A() == 1 ){ | |
494 | localbmin = _bbs.beam2().nuclearRadius() + 0.7; | |
495 | } | |
496 | if( _bbs.beam2().A() == 1 ){ | |
497 | localbmin = _bbs.beam1().nuclearRadius() + 0.7; | |
498 | } | |
499 | integratedflux = nepoint(energy,localbmin); | |
500 | } | |
501 | }else{ | |
502 | for (int jb = 1; jb<=nbstep;jb++){ | |
503 | bold=biter; | |
504 | biter=biter*bmult; | |
505 | // When we get to b>20R_A change methods - just take the photon flux | |
506 | // at the center of the nucleus. | |
507 | if (biter > (10.*RNuc)) | |
508 | { | |
509 | // if there is no nuclear breakup or only hadronic breakup, which only | |
510 | // occurs at smaller b, we can analytically integrate the flux from b~20R_A | |
511 | // to infinity, following Jackson (2nd edition), Eq. 15.54 | |
512 | Xvar=energy*biter/(hbarc*_beamLorentzGamma); | |
513 | // Here, there is nuclear breakup. So, we can't use the integrated flux | |
514 | // However, we can do a single flux calculation, at the center of the nucleus | |
515 | // Eq. 41 of Vidovic, Greiner and Soff, Phys.Rev.C47,2308(1993), among other places | |
516 | // this is the flux per unit area | |
517 | fluxelement = (rZ*rZ*alpha*energy)* | |
518 | (bessel::dbesk1(Xvar))*(bessel::dbesk1(Xvar))/ | |
519 | ((pi*_beamLorentzGamma*hbarc)* | |
520 | (pi*_beamLorentzGamma*hbarc)); | |
521 | ||
522 | }//if biter>10 | |
523 | else{ | |
524 | // integrate over nuclear surface. n.b. this assumes total shadowing - | |
525 | // treat photons hitting the nucleus the same no matter where they strike | |
526 | fluxelement=0.; | |
527 | deltar=RNuc/double(nrstep); | |
528 | riter=-deltar/2.; | |
529 | ||
530 | for (int jr =1; jr<=nrstep;jr++){ | |
531 | riter=riter+deltar; | |
532 | // use symmetry; only integrate from 0 to pi (half circle) | |
533 | deltaphi=pi/double(nphistep); | |
534 | phiiter=0.; | |
535 | ||
536 | for( int jphi=1;jphi<= nphistep;jphi++){ | |
537 | phiiter=(double(jphi)-0.5)*deltaphi; | |
538 | // dist is the distance from the center of the emitting nucleus to the point in question | |
539 | dist=sqrt((biter+riter*cos(phiiter))*(biter+riter* | |
540 | cos(phiiter))+(riter*sin(phiiter))*(riter*sin(phiiter))); | |
541 | Xvar=energy*dist/(hbarc*_beamLorentzGamma); | |
542 | flux_r = (rZ*rZ*alpha*energy)* | |
543 | (bessel::dbesk1(Xvar)*bessel::dbesk1(Xvar))/ | |
544 | ((pi*_beamLorentzGamma*hbarc)* | |
545 | (pi*_beamLorentzGamma*hbarc)); | |
546 | ||
547 | // The surface element is 2.* delta phi* r * delta r | |
548 | // The '2' is because the phi integral only goes from 0 to pi | |
549 | fluxelement=fluxelement+flux_r*2.*deltaphi*riter*deltar; | |
550 | // end phi and r integrations | |
551 | }//for(jphi) | |
552 | }//for(jr) | |
553 | // average fluxelement over the nuclear surface | |
554 | fluxelement=fluxelement/(pi*RNuc*RNuc); | |
555 | }//else | |
556 | // multiply by volume element to get total flux in the volume element | |
557 | fluxelement=fluxelement*2.*pi*biter*(biter-bold); | |
558 | // modulate by the probability of nuclear breakup as f(biter) | |
559 | if (_beamBreakupMode > 1){ | |
560 | fluxelement=fluxelement*_bbs.probabilityOfBreakup(biter); | |
561 | } | |
562 | integratedflux=integratedflux+fluxelement; | |
563 | ||
564 | } //end loop over impact parameter | |
565 | } //end of else (pp, pA, AA) | |
566 | ||
567 | // In lookup table, store k*dN/dk because it changes less | |
568 | // so the interpolation should be better | |
569 | dide[j]=integratedflux*energy; | |
570 | ||
571 | }//end loop over photon energy | |
572 | ||
573 | // for 2nd and subsequent calls, use lookup table immediately | |
574 | ||
575 | L1000f: | |
576 | ||
577 | lEgamma=log(Egamma); | |
578 | if (lEgamma < (lnEmin+dlnE) || lEgamma > lnEmax){ | |
579 | flux_r=0.0; | |
580 | cout<<" ERROR: Egamma outside defined range. Egamma= "<<Egamma | |
581 | <<" "<<lnEmax<<" "<<(lnEmin+dlnE)<<endl; | |
582 | } | |
583 | else{ | |
584 | // >> Egamma between Ilt and Ilt+1 | |
585 | Ilt = int((lEgamma-lnEmin)/dlnE); | |
586 | // >> ln(Egamma) for first point | |
587 | lnElt = lnEmin + Ilt*dlnE; | |
588 | // >> Interpolate | |
589 | flux_r = dide[Ilt] + ((lEgamma-lnElt)/dlnE)*(dide[Ilt+1]- dide[Ilt]); | |
590 | flux_r = flux_r/Egamma; | |
591 | } | |
592 | ||
593 | return flux_r; | |
594 | } | |
595 | ||
596 | ||
597 | //______________________________________________________________________________ | |
598 | double | |
599 | photonNucleusCrossSection::nepoint(const double Egamma, | |
600 | const double bmin) | |
601 | { | |
602 | // Function for the spectrum of virtual photons, | |
603 | // dn/dEgamma, for a point charge q=Ze sweeping | |
604 | // past the origin with velocity gamma | |
605 | // (=1/SQRT(1-(V/c)**2)) integrated over impact | |
606 | // parameter from bmin to infinity | |
607 | // See Jackson eq15.54 Classical Electrodynamics | |
608 | // Declare Local Variables | |
609 | double beta,X,C1,bracket,nepoint_r; | |
610 | ||
611 | beta = sqrt(1.-(1./(_beamLorentzGamma*_beamLorentzGamma))); | |
612 | X = (bmin*Egamma)/(beta*_beamLorentzGamma*hbarc); | |
613 | ||
614 | bracket = -0.5*beta*beta*X*X*(bessel::dbesk1(X)*bessel::dbesk1(X) | |
615 | -bessel::dbesk0(X)*bessel::dbesk0(X)); | |
616 | ||
617 | bracket = bracket+X*bessel::dbesk0(X)*bessel::dbesk1(X); | |
618 | ||
619 | // C1=(2.*double((_bbs.beam1().Z())*(_bbs.beam1().Z()))* | |
620 | // alpha)/pi; | |
621 | ||
622 | // Note: NO Z*Z!! | |
623 | C1=(2.*alpha)/pi; | |
624 | ||
625 | nepoint_r = C1*(1./beta)*(1./beta)*(1./Egamma)*bracket; | |
626 | ||
627 | return nepoint_r; | |
628 | ||
629 | } | |
630 | ||
631 | ||
632 | //______________________________________________________________________________ | |
633 | double | |
634 | photonNucleusCrossSection::sigmagp(const double Wgp) | |
635 | { | |
636 | // >> Function for the gamma-proton --> VectorMeson | |
637 | // >> cross section. Wgp is the gamma-proton CM energy. | |
638 | // >> Unit for cross section: fm**2 | |
639 | ||
640 | double sigmagp_r=0.; | |
641 | ||
642 | switch(_particleType) | |
643 | { | |
644 | case RHO: | |
645 | case RHOZEUS: | |
646 | case FOURPRONG: | |
647 | sigmagp_r=1.E-4*(5.0*exp(0.22*log(Wgp))+26.0*exp(-1.23*log(Wgp))); | |
648 | break; | |
649 | case OMEGA: | |
650 | sigmagp_r=1.E-4*(0.55*exp(0.22*log(Wgp))+18.0*exp(-1.92*log(Wgp))); | |
651 | break; | |
652 | case PHI: | |
653 | sigmagp_r=1.E-4*0.34*exp(0.22*log(Wgp)); | |
654 | break; | |
655 | case JPSI: | |
656 | case JPSI_ee: | |
657 | case JPSI_mumu: | |
658 | sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp)); | |
659 | sigmagp_r*=sigmagp_r; | |
660 | sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp)); | |
661 | // sigmagp_r=1.E-4*0.0015*exp(0.80*log(Wgp)); | |
662 | break; | |
663 | case JPSI2S: | |
664 | case JPSI2S_ee: | |
665 | case JPSI2S_mumu: | |
666 | sigmagp_r=(1.0-((_channelMass+protonMass)*(_channelMass+protonMass))/(Wgp*Wgp)); | |
667 | sigmagp_r*=sigmagp_r; | |
668 | sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp)); | |
669 | sigmagp_r*=0.166; | |
670 | // sigmagp_r=0.166*(1.E-4*0.0015*exp(0.80*log(Wgp))); | |
671 | break; | |
672 | case UPSILON: | |
673 | case UPSILON_ee: | |
674 | case UPSILON_mumu: | |
675 | // >> This is W**1.7 dependence from QCD calculations | |
676 | sigmagp_r=1.E-10*(0.060)*exp(1.70*log(Wgp)); | |
677 | break; | |
678 | case UPSILON2S: | |
679 | case UPSILON2S_ee: | |
680 | case UPSILON2S_mumu: | |
681 | sigmagp_r=1.E-10*(0.0259)*exp(1.70*log(Wgp)); | |
682 | break; | |
683 | case UPSILON3S: | |
684 | case UPSILON3S_ee: | |
685 | case UPSILON3S_mumu: | |
686 | sigmagp_r=1.E-10*(0.0181)*exp(1.70*log(Wgp)); | |
687 | break; | |
688 | default: cout<< "!!! ERROR: Unidentified Vector Meson: "<< _particleType <<endl; | |
689 | } | |
690 | return sigmagp_r; | |
691 | } | |
692 | ||
693 | ||
694 | //______________________________________________________________________________ | |
695 | double | |
696 | photonNucleusCrossSection::sigma_A(const double sig_N) | |
697 | { | |
698 | // Nuclear Cross Section | |
699 | // sig_N,sigma_A in (fm**2) | |
700 | ||
701 | double sum; | |
702 | double b,bmax,Pint,arg,sigma_A_r; | |
703 | ||
704 | int NGAUSS; | |
705 | ||
706 | double xg[17]= | |
707 | {.0, | |
708 | .0483076656877383162,.144471961582796493, | |
709 | .239287362252137075, .331868602282127650, | |
710 | .421351276130635345, .506899908932229390, | |
711 | .587715757240762329, .663044266930215201, | |
712 | .732182118740289680, .794483795967942407, | |
713 | .849367613732569970, .896321155766052124, | |
714 | .934906075937739689, .964762255587506430, | |
715 | .985611511545268335, .997263861849481564 | |
716 | }; | |
717 | ||
718 | double ag[17]= | |
719 | {.0, | |
720 | .0965400885147278006, .0956387200792748594, | |
721 | .0938443990808045654, .0911738786957638847, | |
722 | .0876520930044038111, .0833119242269467552, | |
723 | .0781938957870703065, .0723457941088485062, | |
724 | .0658222227763618468, .0586840934785355471, | |
725 | .0509980592623761762, .0428358980222266807, | |
726 | .0342738629130214331, .0253920653092620595, | |
727 | .0162743947309056706, .00701861000947009660 | |
728 | }; | |
729 | ||
730 | NGAUSS=16; | |
731 | ||
732 | // Check if one or both beams are nuclei | |
733 | int A_1 = _bbs.beam1().A(); | |
734 | int A_2 = _bbs.beam2().A(); | |
735 | if( A_1 == 1 && A_2 == 1)cout<<" This is pp, you should not be here..."<<endl; | |
736 | ||
737 | // CALCULATE P(int) FOR b=0.0 - bmax (fm) | |
738 | bmax = 25.0; | |
739 | sum = 0.; | |
740 | for(int IB=1;IB<=NGAUSS;IB++){ | |
741 | ||
742 | b = 0.5*bmax*xg[IB]+0.5*bmax; | |
743 | ||
744 | if( A_1 == 1 && A_2 != 1){ | |
745 | arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b); | |
746 | }else if(A_2 == 1 && A_1 != 1){ | |
747 | arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b); | |
748 | }else{ | |
749 | arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b); | |
750 | } | |
751 | ||
752 | Pint=1.0-exp(arg); | |
753 | sum=sum+2.*pi*b*Pint*ag[IB]; | |
754 | ||
755 | ||
756 | b = 0.5*bmax*(-xg[IB])+0.5*bmax; | |
757 | ||
758 | if( A_1 == 1 && A_2 != 1){ | |
759 | arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b); | |
760 | }else if(A_2 == 1 && A_1 != 1){ | |
761 | arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b); | |
762 | }else{ | |
763 | arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b); | |
764 | } | |
765 | ||
766 | Pint=1.0-exp(arg); | |
767 | sum=sum+2.*pi*b*Pint*ag[IB]; | |
768 | ||
769 | } | |
770 | ||
771 | sum=0.5*bmax*sum; | |
772 | ||
773 | sigma_A_r=sum; | |
774 | ||
775 | return sigma_A_r; | |
776 | } | |
777 | ||
778 | //______________________________________________________________________________ | |
779 | double | |
780 | photonNucleusCrossSection::sigma_N(const double Wgp) | |
781 | { | |
782 | // Nucleon Cross Section in (fm**2) | |
783 | double cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha); | |
784 | return cs; | |
785 | } | |
786 | ||
787 | ||
788 | //______________________________________________________________________________ | |
789 | double | |
790 | photonNucleusCrossSection::breitWigner(const double W, | |
791 | const double C) | |
792 | { | |
793 | // use simple fixed-width s-wave Breit-Wigner without coherent backgorund for rho' | |
794 | // (PDG '08 eq. 38.56) | |
795 | if(_particleType==FOURPRONG) { | |
796 | if (W < 4.01 * pionChargedMass) | |
797 | return 0; | |
798 | const double termA = _channelMass * _width; | |
799 | const double termA2 = termA * termA; | |
800 | const double termB = W * W - _channelMass * _channelMass; | |
801 | return C * _ANORM * _ANORM * termA2 / (termB * termB + termA2); | |
802 | } | |
803 | ||
804 | // Relativistic Breit-Wigner according to J.D. Jackson, | |
805 | // Nuovo Cimento 34, 6692 (1964), with nonresonant term. A is the strength | |
806 | // of the resonant term and b the strength of the non-resonant | |
807 | // term. C is an overall normalization. | |
808 | ||
809 | double ppi=0.,ppi0=0.,GammaPrim,rat; | |
810 | double aa,bb,cc; | |
811 | ||
812 | double nrbw_r; | |
813 | ||
814 | // width depends on energy - Jackson Eq. A.2 | |
815 | // if below threshold, then return 0. Added 5/3/2001 SRK | |
816 | // 0.5% extra added for safety margin | |
817 | if( _particleType==RHO ||_particleType==RHOZEUS){ | |
818 | if (W < 2.01*pionChargedMass){ | |
819 | nrbw_r=0.; | |
820 | return nrbw_r; | |
821 | } | |
822 | ppi=sqrt( ((W/2.)*(W/2.)) - pionChargedMass * pionChargedMass); | |
823 | ppi0=0.358; | |
824 | } | |
825 | ||
826 | // handle phi-->K+K- properly | |
827 | if (_particleType == PHI){ | |
828 | if (W < 2.*kaonChargedMass){ | |
829 | nrbw_r=0.; | |
830 | return nrbw_r; | |
831 | } | |
832 | ppi=sqrt( ((W/2.)*(W/2.))- kaonChargedMass*kaonChargedMass); | |
833 | ppi0=sqrt( ((_channelMass/2.)*(_channelMass/2.))-kaonChargedMass*kaonChargedMass); | |
834 | } | |
835 | ||
836 | //handle J/Psi-->e+e- properly | |
837 | if (_particleType==JPSI || _particleType==JPSI2S){ | |
838 | if(W<2.*mel){ | |
839 | nrbw_r=0.; | |
840 | return nrbw_r; | |
841 | } | |
842 | ppi=sqrt(((W/2.)*(W/2.))-mel*mel); | |
843 | ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel); | |
844 | } | |
845 | if (_particleType==JPSI_ee){ | |
846 | if(W<2.*mel){ | |
847 | nrbw_r=0.; | |
848 | return nrbw_r; | |
849 | } | |
850 | ppi=sqrt(((W/2.)*(W/2.))-mel*mel); | |
851 | ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel); | |
852 | } | |
853 | if (_particleType==JPSI_mumu){ | |
854 | if(W<2.*muonMass){ | |
855 | nrbw_r=0.; | |
856 | return nrbw_r; | |
857 | } | |
858 | ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass); | |
859 | ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass); | |
860 | } | |
861 | if (_particleType==JPSI2S_ee){ | |
862 | if(W<2.*mel){ | |
863 | nrbw_r=0.; | |
864 | return nrbw_r; | |
865 | } | |
866 | ppi=sqrt(((W/2.)*(W/2.))-mel*mel); | |
867 | ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel); | |
868 | } | |
869 | if (_particleType==JPSI2S_mumu){ | |
870 | if(W<2.*muonMass){ | |
871 | nrbw_r=0.; | |
872 | return nrbw_r; | |
873 | } | |
874 | ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass); | |
875 | ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass); | |
876 | } | |
877 | ||
878 | if(_particleType==UPSILON || _particleType==UPSILON2S ||_particleType==UPSILON3S ){ | |
879 | if (W<2.*muonMass){ | |
880 | nrbw_r=0.; | |
881 | return nrbw_r; | |
882 | } | |
883 | ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass); | |
884 | ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass); | |
885 | } | |
886 | ||
887 | if(_particleType==UPSILON_mumu || _particleType==UPSILON2S_mumu ||_particleType==UPSILON3S_mumu ){ | |
888 | if (W<2.*muonMass){ | |
889 | nrbw_r=0.; | |
890 | return nrbw_r; | |
891 | } | |
892 | ppi=sqrt(((W/2.)*(W/2.))-muonMass*muonMass); | |
893 | ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-muonMass*muonMass); | |
894 | } | |
895 | ||
896 | if(_particleType==UPSILON_ee || _particleType==UPSILON2S_ee ||_particleType==UPSILON3S_ee ){ | |
897 | if (W<2.*mel){ | |
898 | nrbw_r=0.; | |
899 | return nrbw_r; | |
900 | } | |
901 | ppi=sqrt(((W/2.)*(W/2.))-mel*mel); | |
902 | ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-mel*mel); | |
903 | } | |
904 | ||
905 | if(ppi==0.&&ppi0==0.) | |
906 | cout<<"Improper Gammaacrosssection::breitwigner, ppi&ppi0=0."<<endl; | |
907 | ||
908 | rat=ppi/ppi0; | |
909 | GammaPrim=_width*(_channelMass/W)*rat*rat*rat; | |
910 | ||
911 | aa=_ANORM*sqrt(GammaPrim*_channelMass*W); | |
912 | bb=W*W-_channelMass*_channelMass; | |
913 | cc=_channelMass*GammaPrim; | |
914 | ||
915 | // First real part squared | |
916 | nrbw_r = (( (aa*bb)/(bb*bb+cc*cc) + _BNORM)*( (aa*bb)/(bb*bb+cc*cc) + _BNORM)); | |
917 | ||
918 | // Then imaginary part squared | |
919 | nrbw_r = nrbw_r + (( (aa*cc)/(bb*bb+cc*cc) )*( (aa*cc)/(bb*bb+cc*cc) )); | |
920 | ||
921 | // Alternative, a simple, no-background BW, following J. Breitweg et al. | |
922 | // Eq. 15 of Eur. Phys. J. C2, 247 (1998). SRK 11/10/2000 | |
923 | // nrbw_r = (_ANORM*_mass*GammaPrim/(bb*bb+cc*cc))**2 | |
924 | ||
925 | nrbw_r = C*nrbw_r; | |
926 | ||
927 | return nrbw_r; | |
928 | } |