]>
Commit | Line | Data |
---|---|---|
da32329d AM |
1 | /////////////////////////////////////////////////////////////////////////// |
2 | // | |
3 | // Copyright 2010 | |
4 | // | |
5 | // This file is part of starlight. | |
6 | // | |
7 | // starlight is free software: you can redistribute it and/or modify | |
8 | // it under the terms of the GNU General Public License as published by | |
9 | // the Free Software Foundation, either version 3 of the License, or | |
10 | // (at your option) any later version. | |
11 | // | |
12 | // starlight is distributed in the hope that it will be useful, | |
13 | // but WITHOUT ANY WARRANTY; without even the implied warranty of | |
14 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
15 | // GNU General Public License for more details. | |
16 | // | |
17 | // You should have received a copy of the GNU General Public License | |
18 | // along with starlight. If not, see <http://www.gnu.org/licenses/>. | |
19 | // | |
20 | /////////////////////////////////////////////////////////////////////////// | |
21 | // | |
22 | // File and Version Information: | |
45d54d9a | 23 | // $Rev:: 169 $: revision of last commit |
24 | // $Author:: jnystrand $: author of last commit | |
25 | // $Date:: 2014-03-28 17:17:37 +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 | ||
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 | ||
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 | //______________________________________________________________________________ | |
635 | double | |
636 | photonNucleusCrossSection::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 | //______________________________________________________________________________ | |
670 | double | |
671 | photonNucleusCrossSection::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 | //______________________________________________________________________________ | |
732 | double | |
733 | photonNucleusCrossSection::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 | //______________________________________________________________________________ | |
816 | double | |
817 | photonNucleusCrossSection::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 | //______________________________________________________________________________ | |
826 | double | |
827 | photonNucleusCrossSection::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 | |
854 | if( _particleType==RHO ||_particleType==RHOZEUS){ | |
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 | } |