Line data Source code
1 : #include "crpropa/module/PhotoPionProduction.h"
2 : #include "crpropa/Units.h"
3 : #include "crpropa/ParticleID.h"
4 : #include "crpropa/Random.h"
5 :
6 : #include "kiss/convert.h"
7 : #include "kiss/logger.h"
8 : #include "sophia.h"
9 :
10 : #include <limits>
11 : #include <cmath>
12 : #include <sstream>
13 : #include <fstream>
14 : #include <stdexcept>
15 :
16 : namespace crpropa {
17 :
18 11 : PhotoPionProduction::PhotoPionProduction(ref_ptr<PhotonField> field, bool photons, bool neutrinos, bool electrons, bool antiNucleons, double l, bool redshift) {
19 11 : havePhotons = photons;
20 11 : haveNeutrinos = neutrinos;
21 11 : haveElectrons = electrons;
22 11 : haveAntiNucleons = antiNucleons;
23 11 : haveRedshiftDependence = redshift;
24 11 : limit = l;
25 11 : setPhotonField(field);
26 11 : }
27 :
28 22 : void PhotoPionProduction::setPhotonField(ref_ptr<PhotonField> field) {
29 22 : photonField = field;
30 22 : std::string fname = photonField->getFieldName();
31 22 : if (haveRedshiftDependence) {
32 0 : if (photonField->hasRedshiftDependence() == false){
33 0 : std::cout << "PhotoPionProduction: tabulated redshift dependence not needed for " + fname + ", switching off" << std::endl;
34 0 : haveRedshiftDependence = false;
35 : }
36 : else {
37 0 : KISS_LOG_WARNING << "PhotoPionProduction: You are using the 2-dimensional tabulated redshift evolution, which is not available for other interactions. To be consistent across all interactions you may deactivate this <setHaveRedshiftDependence(False)>.";
38 : }
39 : }
40 :
41 22 : setDescription("PhotoPionProduction: " + fname);
42 22 : if (haveRedshiftDependence){
43 0 : initRate(getDataPath("PhotoPionProduction/rate_" + fname.replace(0, 3, "IRBz") + ".txt"));
44 : }
45 : else
46 66 : initRate(getDataPath("PhotoPionProduction/rate_" + fname + ".txt"));
47 22 : }
48 :
49 1 : void PhotoPionProduction::setHavePhotons(bool b) {
50 1 : havePhotons = b;
51 1 : }
52 :
53 0 : void PhotoPionProduction::setHaveElectrons(bool b) {
54 0 : haveElectrons = b;
55 0 : }
56 :
57 0 : void PhotoPionProduction::setHaveNeutrinos(bool b) {
58 0 : haveNeutrinos = b;
59 0 : }
60 :
61 0 : void PhotoPionProduction::setHaveAntiNucleons(bool b) {
62 0 : haveAntiNucleons = b;
63 0 : }
64 :
65 0 : void PhotoPionProduction::setHaveRedshiftDependence(bool b) {
66 0 : haveRedshiftDependence = b;
67 0 : setPhotonField(photonField);
68 0 : }
69 :
70 0 : void PhotoPionProduction::setLimit(double l) {
71 0 : limit = l;
72 0 : }
73 :
74 22 : void PhotoPionProduction::initRate(std::string filename) {
75 : // clear previously loaded tables
76 : tabLorentz.clear();
77 : tabRedshifts.clear();
78 : tabProtonRate.clear();
79 : tabNeutronRate.clear();
80 :
81 22 : std::ifstream infile(filename.c_str());
82 22 : if (!infile.good())
83 0 : throw std::runtime_error("PhotoPionProduction: could not open file " + filename);
84 :
85 22 : if (haveRedshiftDependence) {
86 : double zOld = -1, aOld = -1;
87 0 : while (infile.good()) {
88 0 : if (infile.peek() == '#') {
89 0 : infile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
90 0 : continue;
91 : }
92 : double z, a, b, c;
93 : infile >> z >> a >> b >> c;
94 0 : if (!infile)
95 : break;
96 0 : if (z > zOld) {
97 0 : tabRedshifts.push_back(z);
98 0 : zOld = z;
99 : }
100 0 : if (a > aOld) {
101 0 : tabLorentz.push_back(pow(10, a));
102 0 : aOld = a;
103 : }
104 0 : tabProtonRate.push_back(b / Mpc);
105 0 : tabNeutronRate.push_back(c / Mpc);
106 : }
107 : } else {
108 5610 : while (infile.good()) {
109 5610 : if (infile.peek() == '#') {
110 66 : infile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
111 66 : continue;
112 : }
113 : double a, b, c;
114 : infile >> a >> b >> c;
115 5544 : if (!infile)
116 : break;
117 5522 : tabLorentz.push_back(pow(10, a));
118 5522 : tabProtonRate.push_back(b / Mpc);
119 5522 : tabNeutronRate.push_back(c / Mpc);
120 : }
121 : }
122 :
123 22 : infile.close();
124 22 : }
125 :
126 995 : double PhotoPionProduction::nucleonMFP(double gamma, double z, bool onProton) const {
127 995 : const std::vector<double> &tabRate = (onProton)? tabProtonRate : tabNeutronRate;
128 :
129 : // scale nucleus energy instead of background photon energy
130 995 : gamma *= (1 + z);
131 995 : if (gamma < tabLorentz.front() or (gamma > tabLorentz.back()))
132 : return std::numeric_limits<double>::max();
133 :
134 : double rate;
135 995 : if (haveRedshiftDependence)
136 0 : rate = interpolate2d(z, gamma, tabRedshifts, tabLorentz, tabRate);
137 : else
138 995 : rate = interpolate(gamma, tabLorentz, tabRate) * photonField->getRedshiftScaling(z);
139 :
140 : // cosmological scaling
141 995 : rate *= pow_integer<2>(1 + z);
142 :
143 995 : return 1. / rate;
144 : }
145 :
146 995 : double PhotoPionProduction::nucleiModification(int A, int X) const {
147 995 : if (A == 1)
148 : return 1.;
149 209 : if (A <= 8)
150 209 : return 0.85 * pow(X, 2. / 3.);
151 0 : return 0.85 * X;
152 : }
153 :
154 873 : void PhotoPionProduction::process(Candidate *candidate) const {
155 873 : double step = candidate->getCurrentStep();
156 873 : double z = candidate->getRedshift();
157 : // the loop is processed at least once for limiting the next step
158 : do {
159 : // check if nucleus
160 892 : int id = candidate->current.getId();
161 892 : if (!isNucleus(id))
162 : return;
163 :
164 : // find interaction with minimum random distance
165 891 : Random &random = Random::instance();
166 : double randDistance = std::numeric_limits<double>::max();
167 : double meanFreePath;
168 : double totalRate = 0;
169 : bool onProton = true; // interacting particle: proton or neutron
170 :
171 891 : int A = massNumber(id);
172 891 : int Z = chargeNumber(id);
173 891 : int N = A - Z;
174 891 : double gamma = candidate->current.getLorentzFactor();
175 :
176 : // check for interaction on protons
177 891 : if (Z > 0) {
178 867 : meanFreePath = nucleonMFP(gamma, z, true) / nucleiModification(A, Z);
179 867 : randDistance = -log(random.rand()) * meanFreePath;
180 867 : totalRate += 1. / meanFreePath;
181 : }
182 : // check for interaction on neutrons
183 891 : if (N > 0) {
184 128 : meanFreePath = nucleonMFP(gamma, z, false) / nucleiModification(A, N);
185 128 : totalRate += 1. / meanFreePath;
186 128 : double d = -log(random.rand()) * meanFreePath;
187 128 : if (d < randDistance) {
188 : randDistance = d;
189 : onProton = false;
190 : }
191 : }
192 :
193 : // check if interaction does not happen
194 891 : if (step < randDistance) {
195 872 : if (totalRate > 0.)
196 872 : candidate->limitNextStep(limit / totalRate);
197 872 : return;
198 : }
199 :
200 : // interact and repeat with remaining step
201 19 : performInteraction(candidate, onProton);
202 19 : step -= randDistance;
203 19 : } while (step > 0);
204 : }
205 :
206 29 : void PhotoPionProduction::performInteraction(Candidate *candidate, bool onProton) const {
207 29 : int id = candidate->current.getId();
208 29 : int A = massNumber(id);
209 29 : int Z = chargeNumber(id);
210 29 : double E = candidate->current.getEnergy();
211 29 : double EpA = E / A;
212 29 : double z = candidate->getRedshift();
213 :
214 : // SOPHIA simulates interactions only for protons / neutrons.
215 : // For anti-protons / neutrons assume charge symmetry and change all
216 : // interaction products from particle <--> anti-particle (sign)
217 29 : int sign = (id > 0) ? 1 : -1;
218 :
219 : // check if below SOPHIA's energy threshold
220 29 : double E_threshold = (photonField->getFieldName() == "CMB") ? 3.72e18 * eV : 5.83e15 * eV;
221 29 : if (EpA * (1 + z) < E_threshold)
222 0 : return;
223 :
224 : // SOPHIA - input:
225 29 : int nature = 1 - static_cast<int>(onProton); // 0=proton, 1=neutron
226 29 : double Ein = EpA / GeV; // GeV is the SOPHIA standard unit
227 29 : double eps = sampleEps(onProton, EpA, z) / GeV; // GeV for SOPHIA
228 :
229 : // SOPHIA - output:
230 : double outputEnergy[5][2000]; // [GeV/c, GeV/c, GeV/c, GeV, GeV/c^2]
231 : int outPartID[2000];
232 : int nParticles;
233 :
234 58 : #pragma omp critical(SophiaEvent)
235 : {
236 29 : sophiaevent_(nature, Ein, eps, outputEnergy, outPartID, nParticles);
237 : }
238 :
239 29 : Random &random = Random::instance();
240 29 : Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
241 : std::vector<int> pnType; // filled with either 13 (proton) or 14 (neutron)
242 : std::vector<double> pnEnergy; // corresponding energies of proton or neutron
243 29 : if (nParticles == 0)
244 : return;
245 170 : for (int i = 0; i < nParticles; i++) { // loop over out-going particles
246 141 : double Eout = outputEnergy[3][i] * GeV; // only the energy is used; could be changed for more detail
247 141 : int pType = outPartID[i];
248 141 : switch (pType) {
249 29 : case 13: // proton
250 : case 14: // neutron
251 : // proton and neutron data is taken to determine primary particle in a later step
252 29 : pnType.push_back(pType);
253 29 : pnEnergy.push_back(Eout);
254 : break;
255 0 : case -13: // anti-proton
256 : case -14: // anti-neutron
257 0 : if (haveAntiNucleons)
258 : try
259 : {
260 0 : candidate->addSecondary(-sign * nucleusId(1, 14 + pType), Eout, pos, 1., interactionTag);
261 : }
262 0 : catch (std::runtime_error &e)
263 : {
264 0 : KISS_LOG_ERROR<< "Something went wrong in the PhotoPionProduction (anti-nucleon production)\n" << "Something went wrong in the PhotoPionProduction\n"<< "Please report this error on https://github.com/CRPropa/CRPropa3/issues including your simulation setup and the following random seed:\n" << Random::instance().getSeed_base64();
265 0 : throw;
266 0 : }
267 : break;
268 40 : case 1: // photon
269 40 : if (havePhotons)
270 18 : candidate->addSecondary(22, Eout, pos, 1., interactionTag);
271 : break;
272 12 : case 2: // positron
273 12 : if (haveElectrons)
274 1 : candidate->addSecondary(sign * -11, Eout, pos, 1., interactionTag);
275 : break;
276 6 : case 3: // electron
277 6 : if (haveElectrons)
278 0 : candidate->addSecondary(sign * 11, Eout, pos, 1., interactionTag);
279 : break;
280 12 : case 15: // nu_e
281 12 : if (haveNeutrinos)
282 1 : candidate->addSecondary(sign * 12, Eout, pos, 1., interactionTag);
283 : break;
284 6 : case 16: // anti-nu_e
285 6 : if (haveNeutrinos)
286 0 : candidate->addSecondary(sign * -12, Eout, pos, 1., interactionTag);
287 : break;
288 18 : case 17: // nu_mu
289 18 : if (haveNeutrinos)
290 1 : candidate->addSecondary(sign * 14, Eout, pos, 1., interactionTag);
291 : break;
292 18 : case 18: // anti-nu_mu
293 18 : if (haveNeutrinos)
294 1 : candidate->addSecondary(sign * -14, Eout, pos, 1., interactionTag);
295 : break;
296 0 : default:
297 0 : throw std::runtime_error("PhotoPionProduction: unexpected particle " + kiss::str(pType));
298 : }
299 : }
300 29 : double maxEnergy = *std::max_element(pnEnergy.begin(), pnEnergy.end()); // criterion for being declared primary
301 58 : for (int i = 0; i < pnEnergy.size(); ++i) {
302 29 : if (pnEnergy[i] == maxEnergy) { // nucleon is primary particle
303 29 : if (A == 1) {
304 : // single interacting nucleon
305 26 : candidate->current.setEnergy(pnEnergy[i]);
306 : try
307 : {
308 26 : candidate->current.setId(sign * nucleusId(1, 14 - pnType[i]));
309 : }
310 0 : catch (std::runtime_error &e)
311 : {
312 0 : KISS_LOG_ERROR<< "Something went wrong in the PhotoPionProduction (primary particle, A==1)\n" << "Please report this error on https://github.com/CRPropa/CRPropa3/issues including your simulation setup and the following random seed:\n" << Random::instance().getSeed_base64();
313 0 : throw;
314 0 : }
315 : } else {
316 : // interacting nucleon is part of nucleus: it is emitted from the nucleus
317 3 : candidate->current.setEnergy(E - EpA);
318 : try
319 : {
320 3 : candidate->current.setId(sign * nucleusId(A - 1, Z - int(onProton)));
321 3 : candidate->addSecondary(sign * nucleusId(1, 14 - pnType[i]), pnEnergy[i], pos, 1., interactionTag);
322 : }
323 0 : catch (std::runtime_error &e)
324 : {
325 0 : KISS_LOG_ERROR<< "Something went wrong in the PhotoPionProduction (primary particle, A!=1)\n" << "Please report this error on https://github.com/CRPropa/CRPropa3/issues including your simulation setup and the following random seed:\n" << Random::instance().getSeed_base64();
326 0 : throw;
327 0 : }
328 : }
329 : } else { // nucleon is secondary proton or neutron
330 0 : candidate->addSecondary(sign * nucleusId(1, 14 - pnType[i]), pnEnergy[i], pos, 1., interactionTag);
331 : }
332 : }
333 29 : }
334 :
335 0 : double PhotoPionProduction::lossLength(int id, double gamma, double z) {
336 0 : int A = massNumber(id);
337 0 : int Z = chargeNumber(id);
338 0 : int N = A - Z;
339 :
340 : double lossRate = 0;
341 0 : if (Z > 0)
342 0 : lossRate += 1 / nucleonMFP(gamma, z, true) * nucleiModification(A, Z);
343 0 : if (N > 0)
344 0 : lossRate += 1 / nucleonMFP(gamma, z, false) * nucleiModification(A, N);
345 :
346 : // approximate the relative energy loss
347 : // - nucleons keep the fraction of mass to delta-resonance mass
348 : // - nuclei lose the energy 1/A the interacting nucleon is carrying
349 0 : double relativeEnergyLoss = (A == 1) ? 1 - 938. / 1232. : 1. / A;
350 0 : lossRate *= relativeEnergyLoss;
351 :
352 : // scaling factor: interaction rate --> energy loss rate
353 0 : lossRate *= (1 + z);
354 :
355 0 : return 1. / lossRate;
356 : }
357 :
358 0 : SophiaEventOutput PhotoPionProduction::sophiaEvent(bool onProton, double Ein, double eps) const {
359 : // SOPHIA - input:
360 0 : int nature = 1 - static_cast<int>(onProton); // 0=proton, 1=neutron
361 0 : Ein /= GeV; // GeV is the SOPHIA standard unit
362 0 : eps /= GeV; // GeV for SOPHIA
363 :
364 : // SOPHIA - output:
365 : double outputEnergy[5][2000]; // [Px GeV/c, Py GeV/c, Pz GeV/c, E GeV, m0 GeV/c^2]
366 : int outPartID[2000];
367 : int nParticles;
368 :
369 0 : sophiaevent_(nature, Ein, eps, outputEnergy, outPartID, nParticles);
370 :
371 : // convert SOPHIA IDs to PDG naming convention & create particles
372 : SophiaEventOutput output;
373 0 : output.nParticles = nParticles;
374 0 : for (int i = 0; i < nParticles; ++i) {
375 0 : int id = 0;
376 0 : int partType = outPartID[i];
377 0 : switch (partType) {
378 0 : case 13: // proton
379 : case 14: // neutron
380 0 : id = nucleusId(1, 14 - partType);
381 0 : break;
382 0 : case -13: // anti-proton
383 : case -14: // anti-neutron
384 0 : id = -nucleusId(1, 14 + partType);
385 0 : break;
386 0 : case 1: // photon
387 0 : id = 22;
388 0 : break;
389 0 : case 2: // positron
390 0 : id = -11;
391 0 : break;
392 0 : case 3: // electron
393 0 : id = 11;
394 0 : break;
395 0 : case 15: // nu_e
396 0 : id = 12;
397 0 : break;
398 0 : case 16: // anti-nu_e
399 0 : id = -12;
400 0 : break;
401 0 : case 17: // nu_mu
402 0 : id = 14;
403 0 : break;
404 0 : case 18: // anti-nu_mu
405 0 : id = -14;
406 0 : break;
407 0 : default:
408 0 : throw std::runtime_error("PhotoPionProduction: unexpected particle " + kiss::str(partType));
409 : }
410 0 : output.energy.push_back(outputEnergy[3][i] * GeV); // only the energy is used; could be changed for more detail
411 0 : output.id.push_back(id);
412 : }
413 0 : return output;
414 : }
415 :
416 29 : double PhotoPionProduction::sampleEps(bool onProton, double E, double z) const {
417 : // sample eps between epsMin ... epsMax
418 29 : double Ein = E / GeV;
419 58 : double epsMin = std::max(photonField -> getMinimumPhotonEnergy(z) / eV, epsMinInteraction(onProton, Ein));
420 29 : double epsMax = photonField -> getMaximumPhotonEnergy(z) / eV;
421 29 : double pEpsMax = probEpsMax(onProton, Ein, z, epsMin, epsMax);
422 :
423 29 : Random &random = Random::instance();
424 6111 : for (int i = 0; i < 1000000; i++) {
425 6111 : double eps = epsMin + random.rand() * (epsMax - epsMin);
426 6111 : double pEps = probEps(eps, onProton, Ein, z);
427 6111 : if (random.rand() * pEpsMax < pEps)
428 29 : return eps * eV;
429 : }
430 0 : throw std::runtime_error("error: no photon found in sampleEps, please make sure that photon field provides photons for the interaction by adapting the energy range of the tabulated photon field.");
431 : }
432 :
433 29 : double PhotoPionProduction::epsMinInteraction(bool onProton, double Ein) const {
434 : // labframe energy of least energetic photon where PPP can occur
435 : // this kind-of ties samplingEps to the PPP and SOPHIA
436 29 : const double m = mass(onProton);
437 29 : const double p = momentum(onProton, Ein);
438 29 : double epsMin = 1.e9 * (1.1646 - m * m) / 2. / (Ein + p); // eV
439 29 : return epsMin;
440 : }
441 :
442 30 : double PhotoPionProduction::probEpsMax(bool onProton, double Ein, double z, double epsMin, double epsMax) const {
443 : // find pEpsMax by testing photon energies (eps) for their interaction
444 : // probabilities (p) in order to find the maximum (max) probability
445 : const int nrSteps = 100;
446 : double pEpsMaxTested = 0.;
447 : double step = 0.;
448 30 : if (sampleLog){
449 : // sample in logspace with stepsize that is at max Δlog(E/eV) = 0.01 or otherwise dep. on size of energy range with nrSteps+1 steps log. equidis. spaced
450 30 : step = std::min(0.01, std::log10(epsMax / epsMin) / nrSteps);
451 : } else
452 0 : step = (epsMax - epsMin) / nrSteps;
453 :
454 : double epsDummy = 0.;
455 : int i = 0;
456 5917 : while (epsDummy < epsMax) {
457 5887 : if (sampleLog)
458 5887 : epsDummy = epsMin * pow(10, step * i);
459 : else
460 0 : epsDummy = epsMin + step * i;
461 5887 : double p = probEps(epsDummy, onProton, Ein, z);
462 5887 : if(p > pEpsMaxTested)
463 : pEpsMaxTested = p;
464 5887 : i++;
465 : }
466 : // the following factor corrects for only trying to find the maximum on nrIteration photon energies
467 : // the factor should be determined in convergence tests
468 30 : double pEpsMax = pEpsMaxTested * correctionFactor;
469 :
470 30 : if(pEpsMax == 0) {
471 0 : KISS_LOG_WARNING << "pEpsMax is 0 in the following configuration: \n"
472 0 : << "\t" << "onProton: " << onProton << "\n"
473 0 : << "\t" << "Ein: " << Ein << " [GeV] \n"
474 0 : << "\t" << "epsRange [eV] " << epsMin << "\t" << epsMax << "\n"
475 0 : << "\t" << "redshift: " << z << "\n"
476 0 : << "\t" << "sample Log " << sampleLog << " with step " << step << " [eV] \n";
477 : }
478 :
479 30 : return pEpsMax;
480 : }
481 :
482 11998 : double PhotoPionProduction::probEps(double eps, bool onProton, double Ein, double z) const {
483 : // probEps returns "probability to encounter a photon of energy eps", given a primary nucleon
484 : // note, probEps does not return a normalized probability [0,...,1]
485 11998 : double photonDensity = photonField->getPhotonDensity(eps * eV, z) * ccm / eps;
486 11998 : if (photonDensity != 0.) {
487 11997 : const double p = momentum(onProton, Ein);
488 11997 : const double sMax = mass(onProton) * mass(onProton) + 2. * eps * (Ein + p) / 1.e9;
489 11997 : if (sMax <= sMin())
490 : return 0;
491 107703 : double sIntegr = gaussInt([this, onProton](double s) { return this->functs(s, onProton); }, sMin(), sMax);
492 11967 : return photonDensity * sIntegr / eps / eps / p / 8. * 1.e18 * 1.e6;
493 : }
494 : return 0;
495 : }
496 :
497 12026 : double PhotoPionProduction::momentum(bool onProton, double Ein) const {
498 12026 : const double m = mass(onProton);
499 12026 : const double momentumHadron = sqrt(Ein * Ein - m * m); // GeV/c
500 12026 : return momentumHadron;
501 : }
502 :
503 191472 : double PhotoPionProduction::crossection(double eps, bool onProton) const {
504 191472 : const double m = mass(onProton);
505 191472 : const double s = m * m + 2. * m * eps;
506 191472 : if (s < sMin())
507 : return 0.;
508 : double cross_res = 0.;
509 : double cross_dir = 0.;
510 : double cross_dir1 = 0.;
511 : double cross_dir2 = 0.;
512 : double sig_res[9];
513 :
514 : // first half of array: 9x proton resonance data | second half of array 9x neutron resonance data
515 : static const double AMRES[18] = {1.231, 1.440, 1.515, 1.525, 1.675, 1.680, 1.690, 1.895, 1.950, 1.231, 1.440, 1.515, 1.525, 1.675, 1.675, 1.690, 1.895, 1.950};
516 : static const double BGAMMA[18] = {5.6, 0.5, 4.6, 2.5, 1.0, 2.1, 2.0, 0.2, 1.0, 6.1, 0.3, 4.0, 2.5, 0.0, 0.2, 2.0, 0.2, 1.0};
517 : static const double WIDTH[18] = {0.11, 0.35, 0.11, 0.1, 0.16, 0.125, 0.29, 0.35, 0.3, 0.11, 0.35, 0.11, 0.1, 0.16, 0.150, 0.29, 0.35, 0.3};
518 : static const double RATIOJ[18] = {1., 0.5, 1., 0.5, 0.5, 1.5, 1., 1.5, 2., 1., 0.5, 1., 0.5, 0.5, 1.5, 1., 1.5, 2.};
519 : static const double AM2[2] = {0.882792, 0.880351};
520 :
521 191472 : const int idx = onProton? 0 : 9;
522 : double SIG0[9];
523 1914720 : for (int i = 0; i < 9; ++i) {
524 1723248 : SIG0[i] = 4.893089117 / AM2[int(onProton)] * RATIOJ[i + idx] * BGAMMA[i + idx];
525 : }
526 191472 : if (eps <= 10.) {
527 153453 : cross_res = breitwigner(SIG0[0], WIDTH[0 + idx], AMRES[0 + idx], eps, onProton) * Ef(eps, 0.152, 0.17);
528 : sig_res[0] = cross_res;
529 1381077 : for (int i = 1; i < 9; ++i) {
530 1227624 : sig_res[i] = breitwigner(SIG0[i], WIDTH[i + idx], AMRES[i + idx], eps, onProton) * Ef(eps, 0.15, 0.38);
531 1227624 : cross_res += sig_res[i];
532 : }
533 : // direct channel
534 153453 : if ((eps > 0.1) && (eps < 0.6)) {
535 64299 : cross_dir1 = 92.7 * Pl(eps, 0.152, 0.25, 2.0) // single pion production
536 64299 : + 40. * std::exp(-(eps - 0.29) * (eps - 0.29) / 0.002)
537 64299 : - 15. * std::exp(-(eps - 0.37) * (eps - 0.37) / 0.002);
538 : } else {
539 89154 : cross_dir1 = 92.7 * Pl(eps, 0.152, 0.25, 2.0); // single pion production
540 : }
541 153453 : cross_dir2 = 37.7 * Pl(eps, 0.4, 0.6, 2.0); // double pion production
542 153453 : cross_dir = cross_dir1 + cross_dir2;
543 : }
544 : // fragmentation 2:
545 191472 : double cross_frag2 = onProton? 80.3 : 60.2;
546 191472 : cross_frag2 *= Ef(eps, 0.5, 0.1) * std::pow(s, -0.34);
547 : // multipion production/fragmentation 1 cross section
548 : double cs_multidiff = 0.;
549 : double cs_multi = 0.;
550 : double cross_diffr1 = 0.;
551 : double cross_diffr2 = 0.;
552 : double cross_diffr = 0.;
553 191472 : if (eps > 0.85) {
554 116083 : double ss1 = (eps - 0.85) / 0.69;
555 116083 : double ss2 = onProton? 29.3 : 26.4;
556 116083 : ss2 *= std::pow(s, -0.34) + 59.3 * std::pow(s, 0.095);
557 116083 : cs_multidiff = (1. - std::exp(-ss1)) * ss2;
558 116083 : cs_multi = 0.89 * cs_multidiff;
559 : // diffractive scattering:
560 : cross_diffr1 = 0.099 * cs_multidiff;
561 : cross_diffr2 = 0.011 * cs_multidiff;
562 116083 : cross_diffr = 0.11 * cs_multidiff;
563 : // **************************************
564 116083 : ss1 = std::pow(eps - 0.85, 0.75) / 0.64;
565 116083 : ss2 = 74.1 * std::pow(eps, -0.44) + 62. * std::pow(s, 0.08);
566 116083 : double cs_tmp = 0.96 * (1. - std::exp(-ss1)) * ss2;
567 116083 : cross_diffr1 = 0.14 * cs_tmp;
568 116083 : cross_diffr2 = 0.013 * cs_tmp;
569 116083 : double cs_delta = cross_frag2 - (cross_diffr1 + cross_diffr2 - cross_diffr);
570 116083 : if (cs_delta < 0.) {
571 : cross_frag2 = 0.;
572 0 : cs_multi += cs_delta;
573 : } else {
574 : cross_frag2 = cs_delta;
575 : }
576 : cross_diffr = cross_diffr1 + cross_diffr2;
577 116083 : cs_multidiff = cs_multi + cross_diffr;
578 : // in the original SOPHIA code, here is a switch for the return argument.
579 : // Here, only one case (compare in SOPHIA: NDIR=3) is needed.
580 : }
581 191472 : return cross_res + cross_dir + cs_multidiff + cross_frag2;
582 : }
583 :
584 306906 : double PhotoPionProduction::Pl(double eps, double epsTh, double epsMax, double alpha) const {
585 306906 : if (epsTh > eps)
586 : return 0.;
587 254660 : const double a = alpha * epsMax / epsTh;
588 254660 : const double prod1 = std::pow((eps - epsTh) / (epsMax - epsTh), a - alpha);
589 254660 : const double prod2 = std::pow(eps / epsMax, -a);
590 254660 : return prod1 * prod2;
591 : }
592 :
593 1572549 : double PhotoPionProduction::Ef(double eps, double epsTh, double w) const {
594 1572549 : const double wTh = w + epsTh;
595 1572549 : if (eps <= epsTh) {
596 : return 0.;
597 1512860 : } else if ((eps > epsTh) && (eps < wTh)) {
598 530142 : return (eps - epsTh) / w;
599 982718 : } else if (eps >= wTh) {
600 : return 1.;
601 : } else {
602 0 : throw std::runtime_error("error in function Ef");
603 : }
604 : }
605 :
606 1381077 : double PhotoPionProduction::breitwigner(double sigma0, double gamma, double DMM, double epsPrime, bool onProton) const {
607 1381077 : const double m = mass(onProton);
608 1381077 : const double s = m * m + 2. * m * epsPrime;
609 1381077 : const double gam2s = gamma * gamma * s;
610 1381077 : return sigma0 * (s / epsPrime / epsPrime) * gam2s / ((s - DMM * DMM) * (s - DMM * DMM) + gam2s);
611 : }
612 :
613 191472 : double PhotoPionProduction::functs(double s, bool onProton) const {
614 191472 : const double m = mass(onProton);
615 191472 : const double factor = s - m * m;
616 191472 : const double epsPrime = factor / 2. / m;
617 191472 : const double sigmaPg = crossection(epsPrime, onProton);
618 191472 : return factor * sigmaPg;
619 : }
620 :
621 1800070 : double PhotoPionProduction::mass(bool onProton) const {
622 1800070 : const double m = onProton ? mass_proton : mass_neutron;
623 1800070 : return m / GeV * c_squared;
624 : }
625 :
626 215436 : double PhotoPionProduction::sMin() const {
627 215436 : return 1.1646; // [GeV^2] head-on collision
628 : }
629 :
630 0 : void PhotoPionProduction::setSampleLog(bool b) {
631 0 : sampleLog = b;
632 0 : }
633 :
634 0 : void PhotoPionProduction::setCorrectionFactor(double factor) {
635 0 : correctionFactor = factor;
636 0 : }
637 :
638 1 : ref_ptr<PhotonField> PhotoPionProduction::getPhotonField() const {
639 1 : return photonField;
640 : }
641 :
642 0 : bool PhotoPionProduction::getHavePhotons() const {
643 0 : return havePhotons;
644 : }
645 :
646 0 : bool PhotoPionProduction::getHaveNeutrinos() const {
647 0 : return haveNeutrinos;
648 : }
649 :
650 0 : bool PhotoPionProduction::getHaveElectrons() const {
651 0 : return haveElectrons;
652 : }
653 :
654 0 : bool PhotoPionProduction::getHaveAntiNucleons() const {
655 0 : return haveAntiNucleons;
656 : }
657 :
658 0 : bool PhotoPionProduction::getHaveRedshiftDependence() const {
659 0 : return haveRedshiftDependence;
660 : }
661 :
662 0 : double PhotoPionProduction::getLimit() const {
663 0 : return limit;
664 : }
665 :
666 0 : bool PhotoPionProduction::getSampleLog() const {
667 0 : return sampleLog;
668 : }
669 :
670 1 : double PhotoPionProduction::getCorrectionFactor() const {
671 1 : return correctionFactor;
672 : }
673 :
674 1 : void PhotoPionProduction::setInteractionTag(std::string tag) {
675 1 : interactionTag = tag;
676 1 : }
677 :
678 2 : std::string PhotoPionProduction::getInteractionTag() const {
679 2 : return interactionTag;
680 : }
681 :
682 : } // namespace crpropa
|
