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          22 :         tabLorentz.clear();
      77          22 :         tabRedshifts.clear();
      78          22 :         tabProtonRate.clear();
      79          22 :         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       65838 : double PhotoPionProduction::nucleonMFP(double gamma, double z, bool onProton) const {
     127       65838 :         const std::vector<double> &tabRate = (onProton)? tabProtonRate : tabNeutronRate;
     128             : 
     129             :         // scale nucleus energy instead of background photon energy
     130       65838 :         gamma *= (1 + z);
     131       65838 :         if (gamma < tabLorentz.front() or (gamma > tabLorentz.back()))
     132             :                 return std::numeric_limits<double>::max();
     133             : 
     134             :         double rate;
     135       65838 :         if (haveRedshiftDependence)
     136           0 :                 rate = interpolate2d(z, gamma, tabRedshifts, tabLorentz, tabRate);
     137             :         else
     138       65838 :                 rate = interpolate(gamma, tabLorentz, tabRate) * photonField->getRedshiftScaling(z);
     139             : 
     140             :         // cosmological scaling
     141       65838 :         rate *= pow_integer<2>(1 + z);
     142             : 
     143       65838 :         return 1. / rate;
     144             : }
     145             : 
     146       65838 : double PhotoPionProduction::nucleiModification(int A, int X) const {
     147       65838 :         if (A == 1)
     148             :                 return 1.;
     149        1161 :         if (A <= 8)
     150         345 :                 return 0.85 * pow(X, 2. / 3.);
     151         816 :         return 0.85 * X;
     152             : }
     153             : 
     154       65205 : void PhotoPionProduction::process(Candidate *candidate) const {
     155       65205 :         double step = candidate->getCurrentStep();
     156       65205 :         double z = candidate->getRedshift();
     157             :         // the loop is processed at least once for limiting the next step
     158             :         do {
     159             :                 // check if nucleus
     160       65260 :                 int id = candidate->current.getId();
     161       65260 :                 if (!isNucleus(id))
     162             :                         return;
     163             : 
     164             :                 // find interaction with minimum random distance
     165       65259 :                 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       65259 :                 int A = massNumber(id);
     172       65259 :                 int Z = chargeNumber(id);
     173       65259 :                 int N = A - Z;
     174       65259 :                 double gamma = candidate->current.getLorentzFactor();
     175             : 
     176             :                 // check for interaction on protons
     177       65259 :                 if (Z > 0) {
     178       64707 :                         meanFreePath = nucleonMFP(gamma, z, true) / nucleiModification(A, Z);
     179       64707 :                         randDistance = -log(random.rand()) * meanFreePath;
     180       64707 :                         totalRate += 1. / meanFreePath;
     181             :                 }
     182             :                 // check for interaction on neutrons
     183       65259 :                 if (N > 0) {
     184        1131 :                         meanFreePath = nucleonMFP(gamma, z, false) / nucleiModification(A, N);
     185        1131 :                         totalRate += 1. / meanFreePath;
     186        1131 :                         double d = -log(random.rand()) * meanFreePath;
     187        1131 :                         if (d < randDistance) {
     188             :                                 randDistance = d;
     189             :                                 onProton = false;
     190             :                         }
     191             :                 }
     192             : 
     193             :                 // check if interaction does not happen
     194       65259 :                 if (step < randDistance) {
     195       65204 :                         if (totalRate > 0.)
     196       65204 :                                 candidate->limitNextStep(limit / totalRate);
     197       65204 :                         return;
     198             :                 }
     199             : 
     200             :                 // interact and repeat with remaining step
     201          55 :                 performInteraction(candidate, onProton);
     202          55 :                 step -= randDistance;
     203          55 :         } while (step > 0);
     204             : }
     205             : 
     206          65 : void PhotoPionProduction::performInteraction(Candidate *candidate, bool onProton) const {
     207          65 :         int id = candidate->current.getId();
     208          65 :         int A = massNumber(id);
     209          65 :         int Z = chargeNumber(id);
     210          65 :         double E = candidate->current.getEnergy();
     211          65 :         double EpA = E / A;
     212          65 :         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          65 :         int sign = (id > 0) ? 1 : -1;
     218             : 
     219             :         // check if below SOPHIA's energy threshold
     220         130 :         double E_threshold = (photonField->getFieldName() == "CMB") ? 3.72e18 * eV : 5.83e15 * eV;
     221          65 :         if (EpA * (1 + z) < E_threshold)
     222           0 :                 return;
     223             : 
     224             :         // SOPHIA - input:
     225          65 :         int nature = 1 - static_cast<int>(onProton);  // 0=proton, 1=neutron
     226          65 :         double Ein = EpA / GeV;  // GeV is the SOPHIA standard unit
     227          65 :         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         130 : #pragma omp critical
     235             :         {
     236          65 :                 sophiaevent_(nature, Ein, eps, outputEnergy, outPartID, nParticles);
     237             :         }
     238             : 
     239          65 :         Random &random = Random::instance();
     240          65 :         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          65 :         if (nParticles == 0)
     244             :                 return;
     245         482 :         for (int i = 0; i < nParticles; i++) { // loop over out-going particles
     246         417 :                 double Eout = outputEnergy[3][i] * GeV; // only the energy is used; could be changed for more detail
     247         417 :                 int pType = outPartID[i];
     248         417 :                 switch (pType) {
     249          66 :                 case 13: // proton
     250             :                 case 14: // neutron
     251             :                         // proton and neutron data is taken to determine primary particle in a later step
     252          66 :                         pnType.push_back(pType);
     253          66 :                         pnEnergy.push_back(Eout);
     254             :                         break;
     255           1 :                 case -13: // anti-proton
     256             :                 case -14: // anti-neutron
     257           1 :                         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          90 :                 case 1: // photon
     269          90 :                         if (havePhotons)
     270          14 :                                 candidate->addSecondary(22, Eout, pos, 1., interactionTag);
     271             :                         break;
     272          48 :                 case 2: // positron
     273          48 :                         if (haveElectrons)
     274           2 :                                 candidate->addSecondary(sign * -11, Eout, pos, 1., interactionTag);
     275             :                         break;
     276          17 :                 case 3: // electron
     277          17 :                         if (haveElectrons)
     278           2 :                                 candidate->addSecondary(sign * 11, Eout, pos, 1., interactionTag);
     279             :                         break;
     280          48 :                 case 15: // nu_e
     281          48 :                         if (haveNeutrinos)
     282           2 :                                 candidate->addSecondary(sign * 12, Eout, pos, 1., interactionTag);
     283             :                         break;
     284          17 :                 case 16: // anti-nu_e
     285          17 :                         if (haveNeutrinos)
     286           2 :                                 candidate->addSecondary(sign * -12, Eout, pos, 1., interactionTag);
     287             :                         break;
     288          65 :                 case 17: // nu_mu
     289          65 :                         if (haveNeutrinos)
     290           4 :                                 candidate->addSecondary(sign * 14, Eout, pos, 1., interactionTag);
     291             :                         break;
     292          65 :                 case 18: // anti-nu_mu
     293          65 :                         if (haveNeutrinos)
     294           4 :                                 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          65 :         double maxEnergy = *std::max_element(pnEnergy.begin(), pnEnergy.end());  // criterion for being declared primary
     301         131 :         for (int i = 0; i < pnEnergy.size(); ++i) {
     302          66 :                 if (pnEnergy[i] == maxEnergy) {  // nucleon is primary particle
     303          65 :                         if (A == 1) {
     304             :                                 // single interacting nucleon
     305          61 :                                 candidate->current.setEnergy(pnEnergy[i]);
     306             :                                 try
     307             :                                 {
     308          61 :                                         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           4 :                                 candidate->current.setEnergy(E - EpA);
     318             :                                 try
     319             :                                 {
     320           4 :                                         candidate->current.setId(sign * nucleusId(A - 1, Z - int(onProton)));
     321           4 :                                         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           1 :                         candidate->addSecondary(sign * nucleusId(1, 14 - pnType[i]), pnEnergy[i], pos, 1., interactionTag);
     331             :                 }
     332             :         }
     333             : }
     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           0 : }
     415             : 
     416          65 : double PhotoPionProduction::sampleEps(bool onProton, double E, double z) const {
     417             :         // sample eps between epsMin ... epsMax
     418          65 :         double Ein = E / GeV;
     419         130 :         double epsMin = std::max(photonField -> getMinimumPhotonEnergy(z) / eV, epsMinInteraction(onProton, Ein));
     420          65 :         double epsMax = photonField -> getMaximumPhotonEnergy(z) / eV;
     421          65 :         double pEpsMax = probEpsMax(onProton, Ein, z, epsMin, epsMax);
     422             : 
     423          65 :         Random &random = Random::instance();
     424       47452 :         for (int i = 0; i < 1000000; i++) {
     425       47452 :                 double eps = epsMin + random.rand() * (epsMax - epsMin);
     426       47452 :                 double pEps = probEps(eps, onProton, Ein, z);
     427       47452 :                 if (random.rand() * pEpsMax < pEps)
     428          65 :                         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          65 : 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          65 :         const double m = mass(onProton);
     437          65 :         const double p = momentum(onProton, Ein);
     438          65 :         double epsMin = 1.e9 * (1.1646 - m * m) / 2. / (Ein + p); // eV
     439          65 :         return epsMin;
     440             : }
     441             : 
     442          66 : 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          66 :         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          66 :                 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       16509 :         while (epsDummy < epsMax) {
     457       16443 :                 if (sampleLog)
     458       16443 :                         epsDummy = epsMin * pow(10, step * i);
     459             :                 else
     460           0 :                         epsDummy = epsMin + step * i;
     461       16443 :                 double p = probEps(epsDummy, onProton, Ein, z);
     462       16443 :                 if(p > pEpsMaxTested)
     463             :                         pEpsMaxTested = p;
     464       16443 :                 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          66 :         double pEpsMax = pEpsMaxTested * correctionFactor;
     469             : 
     470          66 :         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          66 :         return pEpsMax;
     480             : }
     481             : 
     482       63895 : 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       63895 :         double photonDensity = photonField->getPhotonDensity(eps * eV, z) * ccm / eps;
     486       63895 :         if (photonDensity != 0.) {
     487       63873 :                 const double p = momentum(onProton, Ein);
     488       63873 :                 const double sMax = mass(onProton) * mass(onProton) + 2. * eps * (Ein + p) / 1.e9;
     489       63873 :                 if (sMax <= sMin())
     490             :                         return 0;
     491      574263 :                 double sIntegr = gaussInt([this, onProton](double s) { return this->functs(s, onProton); }, sMin(), sMax);
     492       63807 :                 return photonDensity * sIntegr / eps / eps / p / 8. * 1.e18 * 1.e6;
     493             :         }
     494             :         return 0;
     495             : }
     496             : 
     497       63938 : double PhotoPionProduction::momentum(bool onProton, double Ein) const {
     498       63938 :         const double m = mass(onProton);
     499       63938 :         const double momentumHadron = sqrt(Ein * Ein - m * m);  // GeV/c
     500       63938 :         return momentumHadron;
     501             : }
     502             : 
     503     1020912 : double PhotoPionProduction::crossection(double eps, bool onProton) const {
     504     1020912 :         const double m = mass(onProton);
     505     1020912 :         const double s = m * m + 2. * m * eps;
     506     1020912 :         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     1020912 :         const int idx = onProton? 0 : 9;
     522             :         double SIG0[9];
     523    10209120 :         for (int i = 0; i < 9; ++i) {
     524     9188208 :                 SIG0[i] = 4.893089117 / AM2[int(onProton)] * RATIOJ[i + idx] * BGAMMA[i + idx];
     525             :         }
     526     1020912 :         if (eps <= 10.) {
     527      389793 :                 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     3508137 :                 for (int i = 1; i < 9; ++i) {
     530     3118344 :                         sig_res[i] = breitwigner(SIG0[i], WIDTH[i + idx], AMRES[i + idx], eps, onProton) * Ef(eps, 0.15, 0.38);
     531     3118344 :                         cross_res += sig_res[i];
     532             :                 }
     533             :                 // direct channel
     534      389793 :                 if ((eps > 0.1) && (eps < 0.6)) {
     535      141887 :                         cross_dir1 = 92.7 * Pl(eps, 0.152, 0.25, 2.0)  // single pion production
     536      141887 :                                            + 40. * std::exp(-(eps - 0.29) * (eps - 0.29) / 0.002)
     537      141887 :                                            - 15. * std::exp(-(eps - 0.37) * (eps - 0.37) / 0.002);
     538             :                 } else {
     539      247906 :                         cross_dir1 = 92.7 * Pl(eps, 0.152, 0.25, 2.0);  // single pion production
     540             :                 }
     541      389793 :                 cross_dir2 = 37.7 * Pl(eps, 0.4, 0.6, 2.0);  // double pion production
     542      389793 :                 cross_dir = cross_dir1 + cross_dir2;
     543             :         }
     544             :         // fragmentation 2:
     545     1020912 :         double cross_frag2 = onProton? 80.3 : 60.2;
     546     1020912 :         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     1020912 :         if (eps > 0.85) {
     554      852363 :                 double ss1 = (eps - 0.85) / 0.69;
     555      852363 :                 double ss2 = onProton? 29.3 : 26.4;
     556      852363 :                 ss2 *= std::pow(s, -0.34) + 59.3 * std::pow(s, 0.095);
     557      852363 :                 cs_multidiff = (1. - std::exp(-ss1)) * ss2;
     558      852363 :                 cs_multi = 0.89 * cs_multidiff;
     559             :                 // diffractive scattering:
     560             :                 cross_diffr1 = 0.099 * cs_multidiff;
     561             :                 cross_diffr2 = 0.011 * cs_multidiff;
     562      852363 :                 cross_diffr = 0.11 * cs_multidiff;
     563             :                 // **************************************
     564      852363 :                 ss1 = std::pow(eps - 0.85, 0.75) / 0.64;
     565      852363 :                 ss2 = 74.1 * std::pow(eps, -0.44) + 62. * std::pow(s, 0.08);
     566      852363 :                 double cs_tmp = 0.96 * (1. - std::exp(-ss1)) * ss2;
     567      852363 :                 cross_diffr1 = 0.14 * cs_tmp;
     568      852363 :                 cross_diffr2 = 0.013 * cs_tmp;
     569      852363 :                 double cs_delta = cross_frag2 - (cross_diffr1 + cross_diffr2 - cross_diffr);
     570      852363 :                 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      852363 :                 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     1020912 :         return cross_res + cross_dir + cs_multidiff + cross_frag2;
     582             : }
     583             : 
     584      779586 : double PhotoPionProduction::Pl(double eps, double epsTh, double epsMax, double alpha) const {
     585      779586 :         if (epsTh > eps)
     586             :                 return 0.;
     587      665380 :         const double a = alpha * epsMax / epsTh;
     588      665380 :         const double prod1 = std::pow((eps - epsTh) / (epsMax - epsTh), a - alpha);
     589      665380 :         const double prod2 = std::pow(eps / epsMax, -a);
     590      665380 :         return prod1 * prod2;
     591             : }
     592             : 
     593     4529049 : double PhotoPionProduction::Ef(double eps, double epsTh, double w) const {
     594     4529049 :         const double wTh = w + epsTh;
     595     4529049 :         if (eps <= epsTh) {
     596             :                 return 0.;
     597     4398330 :         } else if ((eps > epsTh) && (eps < wTh)) {
     598     1166849 :                 return (eps - epsTh) / w;
     599     3231481 :         } else if (eps >= wTh) {
     600             :                 return 1.;
     601             :         } else {
     602           0 :                 throw std::runtime_error("error in function Ef");
     603             :         }
     604             : }
     605             : 
     606     3508137 : double PhotoPionProduction::breitwigner(double sigma0, double gamma, double DMM, double epsPrime, bool onProton) const {
     607     3508137 :         const double m = mass(onProton);
     608     3508137 :         const double s = m * m + 2. * m * epsPrime;
     609     3508137 :         const double gam2s = gamma * gamma * s;
     610     3508137 :         return sigma0 * (s / epsPrime / epsPrime) * gam2s / ((s - DMM * DMM) * (s - DMM * DMM) + gam2s);
     611             : }
     612             : 
     613     1020912 : double PhotoPionProduction::functs(double s, bool onProton) const {
     614     1020912 :         const double m = mass(onProton);
     615     1020912 :         const double factor = s - m * m;
     616     1020912 :         const double epsPrime = factor / 2. / m;
     617     1020912 :         const double sigmaPg = crossection(epsPrime, onProton);
     618     1020912 :         return factor * sigmaPg;
     619             : }
     620             : 
     621     5741710 : double PhotoPionProduction::mass(bool onProton) const {
     622     5741710 :         const double m =  onProton ? mass_proton : mass_neutron;
     623     5741710 :         return m / GeV * c_squared;
     624             : }
     625             : 
     626     1148592 : double PhotoPionProduction::sMin() const {
     627     1148592 :         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