Line data    Source code

       1             : #include "crpropa/module/EMPairProduction.h"
       2             : #include "crpropa/Units.h"
       3             : #include "crpropa/Random.h"
       4             : 
       5             : #include <fstream>
       6             : #include <limits>
       7             : #include <stdexcept>
       8             : 
       9             : 
      10             : namespace crpropa {
      11             : 
      12             : static const double mec2 = mass_electron * c_squared;
      13             : 
      14           9 : EMPairProduction::EMPairProduction(ref_ptr<PhotonField> photonField, bool haveElectrons, double thinning, double limit) {
      15           9 :         setPhotonField(photonField);
      16           9 :         setThinning(thinning);
      17           9 :         setLimit(limit);
      18           9 :         setHaveElectrons(haveElectrons);
      19           9 : }
      20             : 
      21          33 : void EMPairProduction::setPhotonField(ref_ptr<PhotonField> photonField) {
      22          33 :         this->photonField = photonField;
      23          33 :         std::string fname = photonField->getFieldName();
      24          33 :         setDescription("EMPairProduction: " + fname);
      25          99 :         initRate(getDataPath("EMPairProduction/rate_" + fname + ".txt"));
      26          99 :         initCumulativeRate(getDataPath("EMPairProduction/cdf_" + fname + ".txt"));
      27          33 : }
      28             : 
      29          14 : void EMPairProduction::setHaveElectrons(bool haveElectrons) {
      30          14 :         this->haveElectrons = haveElectrons;
      31          14 : }
      32             : 
      33           9 : void EMPairProduction::setLimit(double limit) {
      34           9 :         this->limit = limit;
      35           9 : }
      36             : 
      37          13 : void EMPairProduction::setThinning(double thinning) {
      38          13 :         this->thinning = thinning;
      39          13 : }
      40             : 
      41          33 : void EMPairProduction::initRate(std::string filename) {
      42          33 :         std::ifstream infile(filename.c_str());
      43             : 
      44          33 :         if (!infile.good())
      45           0 :                 throw std::runtime_error("EMPairProduction: could not open file " + filename);
      46             : 
      47             :         // clear previously loaded interaction rates
      48          33 :         tabEnergy.clear();
      49          33 :         tabRate.clear();
      50             : 
      51        7510 :         while (infile.good()) {
      52        7477 :                 if (infile.peek() != '#') {
      53             :                         double a, b;
      54             :                         infile >> a >> b;
      55        7345 :                         if (infile) {
      56        7312 :                                 tabEnergy.push_back(pow(10, a) * eV);
      57        7312 :                                 tabRate.push_back(b / Mpc);
      58             :                         }
      59             :                 }
      60        7477 :                 infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
      61             :         }
      62          33 :         infile.close();
      63          33 : }
      64             : 
      65          33 : void EMPairProduction::initCumulativeRate(std::string filename) {
      66          33 :         std::ifstream infile(filename.c_str());
      67             : 
      68          33 :         if (!infile.good())
      69           0 :                 throw std::runtime_error("EMPairProduction: could not open file " + filename);
      70             : 
      71             :         // clear previously loaded tables
      72          33 :         tabE.clear();
      73          33 :         tabs.clear();
      74          33 :         tabCDF.clear();
      75             :         
      76             :         // skip header
      77         165 :         while (infile.peek() == '#')
      78         132 :                 infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
      79             : 
      80             :         // read s values in first line
      81             :         double a;
      82             :         infile >> a; // skip first value
      83        3663 :         while (infile.good() and (infile.peek() != '\n')) {
      84             :                 infile >> a;
      85        3630 :                 tabs.push_back(pow(10, a) * eV * eV);
      86             :         }
      87             : 
      88             :         // read all following lines: E, cdf values
      89        7345 :         while (infile.good()) {
      90             :                 infile >> a;
      91        7345 :                 if (!infile)
      92             :                         break;  // end of file
      93        7312 :                 tabE.push_back(pow(10, a) * eV);
      94             :                 std::vector<double> cdf;
      95      811632 :                 for (int i = 0; i < tabs.size(); i++) {
      96             :                         infile >> a;
      97      804320 :                         cdf.push_back(a / Mpc);
      98             :                 }
      99        7312 :                 tabCDF.push_back(cdf);
     100             :         }
     101          33 :         infile.close();
     102          33 : }
     103             : 
     104             : // Hold an data array to interpolate the energy distribution on
     105             : class PPSecondariesEnergyDistribution {
     106             :         private:
     107             :                 std::vector<double> tab_s;
     108             :                 std::vector< std::vector<double> > data;
     109             :                 size_t N;
     110             : 
     111             :         public:
     112             :                 // differential cross section for pair production for x = Epositron/Egamma, compare Lee 96 arXiv:9604098
     113             :                 double dSigmadE_PPx(double x, double beta) {
     114     1000000 :                         double A = (x / (1. - x) + (1. - x) / x );
     115     1000000 :                         double B =  (1. / x + 1. / (1. - x) );
     116        1000 :                         double y = (1 - beta * beta);
     117     1000000 :                         return A + y * B - y * y / 4 * B * B;
     118             :                 }
     119             : 
     120           1 :                 PPSecondariesEnergyDistribution() {
     121           1 :                         N = 1000;
     122             :                         size_t Ns = 1000;
     123             :                         double s_min = 4 * mec2 * mec2;
     124             :                         double s_max = 1e23 * eV * eV;
     125             :                         double dls = log(s_max / s_min) / Ns;
     126           1 :                         data = std::vector< std::vector<double> >(Ns, std::vector<double>(N));
     127           1 :                         tab_s = std::vector<double>(Ns + 1);
     128             : 
     129        1002 :                         for (size_t i = 0; i < Ns + 1; ++i)
     130        1001 :                                 tab_s[i] = s_min * exp(i*dls); // tabulate s bin borders
     131             : 
     132        1001 :                         for (size_t i = 0; i < Ns; i++) {
     133        1000 :                                 double s = s_min * exp(i*dls + 0.5*dls);
     134        1000 :                                 double beta = sqrt(1 - s_min/s);
     135        1000 :                                 double x0 = (1 - beta) / 2;
     136        1000 :                                 double dx = log((1 + beta) / (1 - beta)) / N;
     137             : 
     138             :                                 // cumulative midpoint integration
     139        1000 :                                 std::vector<double> data_i(1000);
     140        1000 :                                 data_i[0] = dSigmadE_PPx(x0, beta) * expm1(dx);
     141     1000000 :                                 for (size_t j = 1; j < N; j++) {
     142      999000 :                                         double x = x0 * exp(j*dx + 0.5*dx);
     143      999000 :                                         double binWidth = exp((j+1)*dx)-exp(j*dx);
     144      999000 :                                         data_i[j] = dSigmadE_PPx(x, beta) * binWidth + data_i[j-1];
     145             :                                 }
     146        1000 :                                 data[i] = data_i;
     147             :                         }
     148           1 :                 }
     149             : 
     150             :                 // sample positron energy from cdf(E, s_kin)
     151         559 :                 double sample(double E0, double s) {
     152             :                         // get distribution for given s
     153         559 :                         size_t idx = std::lower_bound(tab_s.begin(), tab_s.end(), s) - tab_s.begin();
     154         559 :                         if (idx > data.size())
     155             :                                 return NAN;
     156             :                                 
     157         559 :                         std::vector<double> s0 = data[idx];
     158             : 
     159             :                         // draw random bin
     160         559 :                         Random &random = Random::instance();
     161         559 :                         size_t j = random.randBin(s0) + 1;
     162             : 
     163             :                         double s_min = 4. * mec2 * mec2;
     164         559 :                         double beta = sqrtl(1. - s_min / s);
     165         559 :                         double x0 = (1. - beta) / 2.;
     166         559 :                         double dx = log((1 + beta) / (1 - beta)) / N;
     167         559 :                         double binWidth = x0 * (exp(j*dx) - exp((j-1)*dx));
     168         559 :                         if (random.rand() < 0.5)
     169         282 :                                 return E0 * (x0 * exp((j-1) * dx) + binWidth);
     170             :                         else
     171         277 :                                 return E0 * (1 - (x0 * exp((j-1) * dx) + binWidth));
     172             :                 }
     173             : };
     174             : 
     175         559 : void EMPairProduction::performInteraction(Candidate *candidate) const {
     176             :         // scale particle energy instead of background photon energy
     177         559 :         double z = candidate->getRedshift();
     178         559 :         double E = candidate->current.getEnergy() * (1 + z);
     179             : 
     180             :         // check if secondary electron pair needs to be produced
     181         559 :         if (not haveElectrons)
     182           0 :                 return;
     183             : 
     184             :         // check if in tabulated energy range
     185         559 :         if (E < tabE.front() or (E > tabE.back()))
     186             :                 return;
     187             : 
     188             :         // sample the value of s
     189         559 :         Random &random = Random::instance();
     190         559 :         size_t i = closestIndex(E, tabE);  // find closest tabulation point
     191         559 :         size_t j = random.randBin(tabCDF[i]);
     192        1076 :         double lo = std::max(4 * mec2 * mec2, tabs[j-1]);  // first s-tabulation point below min(s_kin) = (2 me c^2)^2; ensure physical value
     193         559 :         double hi = tabs[j];
     194         559 :         double s = lo + random.rand() * (hi - lo);
     195             : 
     196             :         // sample electron / positron energy
     197         559 :         static PPSecondariesEnergyDistribution interpolation;
     198         559 :         double Ee = interpolation.sample(E, s);
     199         559 :         double Ep = E - Ee;
     200         559 :         double f = Ep / E;
     201             : 
     202             :         // for some backgrounds Ee=nan due to precision limitations.
     203         559 :         if (not std::isfinite(Ee) || not std::isfinite(Ep))
     204             :                 return;
     205             : 
     206             :         // photon is lost after interacting
     207         559 :         candidate->setActive(false);
     208             : 
     209             :         // sample random position along current step
     210         559 :         Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
     211             :         // apply sampling
     212         559 :         if (random.rand() < pow(f, thinning)) {
     213         559 :                 double w = 1. / pow(f, thinning);
     214         559 :                 candidate->addSecondary(11, Ep / (1 + z), pos, w, interactionTag);
     215             :         }
     216         559 :         if (random.rand() < pow(1 - f, thinning)){
     217         559 :                 double w = 1. / pow(1 - f, thinning);
     218         559 :                 candidate->addSecondary(-11, Ee / (1 + z), pos, w, interactionTag);  
     219             :         }
     220             : }
     221             : 
     222       66881 : void EMPairProduction::process(Candidate *candidate) const {
     223             :         // check if photon
     224       66881 :         if (candidate->current.getId() != 22)
     225             :                 return;
     226             : 
     227             :         // scale particle energy instead of background photon energy
     228         841 :         double z = candidate->getRedshift();
     229         841 :         double E = candidate->current.getEnergy() * (1 + z);
     230             : 
     231             :         // check if in tabulated energy range
     232         841 :         if ((E < tabEnergy.front()) or (E > tabEnergy.back()))
     233             :                 return;
     234             : 
     235             :         // interaction rate
     236         651 :         double rate = interpolate(E, tabEnergy, tabRate);
     237         651 :         rate *= pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z);
     238             : 
     239             :         // run this loop at least once to limit the step size 
     240         651 :         double step = candidate->getCurrentStep();
     241         651 :         Random &random = Random::instance();
     242             :         do {
     243         651 :                 double randDistance = -log(random.rand()) / rate;
     244             :                 // check for interaction; if it doesn't ocurr, limit next step
     245         651 :                 if (step < randDistance) { 
     246          93 :                         candidate->limitNextStep(limit / rate);
     247             :                 } else {
     248         558 :                         performInteraction(candidate);
     249         558 :                         return;
     250             :                 }
     251          93 :                 step -= randDistance; 
     252          93 :         } while (step > 0.);
     253             : 
     254             : }
     255             : 
     256           1 : void EMPairProduction::setInteractionTag(std::string tag) {
     257           1 :         interactionTag = tag;
     258           1 : }
     259             : 
     260           2 : std::string EMPairProduction::getInteractionTag() const {
     261           2 :         return interactionTag;
     262             : }
     263             : 
     264             : } // namespace crpropa