Line data    Source code

       1             : #include "crpropa/module/EMTripletPairProduction.h"
       2             : #include "crpropa/Units.h"
       3             : #include "crpropa/Random.h"
       4             : 
       5             : #include <fstream>
       6             : #include <limits>
       7             : #include <stdexcept>
       8             : 
       9             : namespace crpropa {
      10             : 
      11             : static const double mec2 = mass_electron * c_squared;
      12             : 
      13           3 : EMTripletPairProduction::EMTripletPairProduction(ref_ptr<PhotonField> photonField, bool haveElectrons, double thinning, double limit) {
      14           3 :         setPhotonField(photonField);
      15           3 :         setHaveElectrons(haveElectrons);
      16           3 :         setLimit(limit);
      17           3 :         setThinning(thinning);
      18           3 : }
      19             : 
      20          15 : void EMTripletPairProduction::setPhotonField(ref_ptr<PhotonField> photonField) {
      21          15 :         this->photonField = photonField;
      22          15 :         std::string fname = photonField->getFieldName();
      23          15 :         setDescription("EMTripletPairProduction: " + fname);
      24          45 :         initRate(getDataPath("EMTripletPairProduction/rate_" + fname + ".txt"));
      25          45 :         initCumulativeRate(getDataPath("EMTripletPairProduction/cdf_" + fname + ".txt"));
      26          15 : }
      27             : 
      28           4 : void EMTripletPairProduction::setHaveElectrons(bool haveElectrons) {
      29           4 :         this->haveElectrons = haveElectrons;
      30           4 : }
      31             : 
      32           3 : void EMTripletPairProduction::setLimit(double limit) {
      33           3 :         this->limit = limit;
      34           3 : }
      35             : 
      36           3 : void EMTripletPairProduction::setThinning(double thinning) {
      37           3 :         this->thinning = thinning;
      38           3 : }
      39             : 
      40          15 : void EMTripletPairProduction::initRate(std::string filename) {
      41          15 :         std::ifstream infile(filename.c_str());
      42             : 
      43          15 :         if (!infile.good())
      44           0 :                 throw std::runtime_error("EMTripletPairProduction: could not open file " + filename);
      45             : 
      46             :         // clear previously loaded interaction rates
      47          15 :         tabEnergy.clear();
      48          15 :         tabRate.clear();
      49             : 
      50        3340 :         while (infile.good()) {
      51        3325 :                 if (infile.peek() != '#') {
      52             :                         double a, b;
      53             :                         infile >> a >> b;
      54        3265 :                         if (infile) {
      55        3250 :                                 tabEnergy.push_back(pow(10, a) * eV);
      56        3250 :                                 tabRate.push_back(b / Mpc);
      57             :                         }
      58             :                 }
      59        3325 :                 infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
      60             :         }
      61          15 :         infile.close();
      62          15 : }
      63             : 
      64          15 : void EMTripletPairProduction::initCumulativeRate(std::string filename) {
      65          15 :         std::ifstream infile(filename.c_str());
      66             : 
      67          15 :         if (!infile.good())
      68           0 :                 throw std::runtime_error(
      69           0 :                                 "EMTripletPairProduction: could not open file " + filename);
      70             : 
      71             :         // clear previously loaded tables
      72          15 :         tabE.clear();
      73          15 :         tabs.clear();
      74          15 :         tabCDF.clear();
      75             :         
      76             :         // skip header
      77          75 :         while (infile.peek() == '#')
      78          60 :                 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        1590 :         while (infile.good() and (infile.peek() != '\n')) {
      84             :                 infile >> a;
      85        1575 :                 tabs.push_back(pow(10, a) * eV * eV);
      86             :         }
      87             : 
      88             :         // read all following lines: E, cdf values
      89        3265 :         while (infile.good()) {
      90             :                 infile >> a;
      91        3265 :                 if (!infile)
      92             :                         break;  // end of file
      93        3250 :                 tabE.push_back(pow(10, a) * eV);
      94             :                 std::vector<double> cdf;
      95      344500 :                 for (int i = 0; i < tabs.size(); i++) {
      96             :                         infile >> a;
      97      341250 :                         cdf.push_back(a / Mpc);
      98             :                 }
      99        3250 :                 tabCDF.push_back(cdf);
     100             :         }
     101          15 :         infile.close();
     102          15 : }
     103             : 
     104           1 : void EMTripletPairProduction::performInteraction(Candidate *candidate) const {
     105           1 :         int id = candidate->current.getId();
     106           1 :         if  (abs(id) != 11)
     107             :                 return;
     108             : 
     109             :         // scale the particle energy instead of background photons
     110           1 :         double z = candidate->getRedshift();
     111           1 :         double E = candidate->current.getEnergy() * (1 + z);
     112             : 
     113           1 :         if (E < tabE.front() or E > tabE.back())
     114             :                 return;
     115             : 
     116             :         // sample the value of eps
     117           1 :         Random &random = Random::instance();
     118           1 :         size_t i = closestIndex(E, tabE);
     119           1 :         size_t j = random.randBin(tabCDF[i]);
     120           1 :         double s_kin = pow(10, log10(tabs[j]) + (random.rand() - 0.5) * 0.1);
     121           1 :         double eps = s_kin / 4. / E; // random background photon energy
     122             : 
     123             :         // Use approximation from A. Mastichiadis et al., Astroph. Journ. 300:178-189 (1986), eq. 30.
     124             :         // This approx is valid only for alpha >=100 where alpha = p0*eps*costheta - E0*eps
     125             :         // For our purposes, me << E0 --> p0~E0 --> alpha = E0*eps*(costheta - 1) >= 100
     126           1 :         double Epp = 5.7e-1 * pow(eps / mec2, -0.56) * pow(E / mec2, 0.44) * mec2;
     127             : 
     128           1 :         double f = Epp / E;
     129             : 
     130           1 :         if (haveElectrons) {
     131           1 :                 Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
     132           1 :                 if (random.rand() < pow(1 - f, thinning)) {
     133           1 :                         double w = 1. / pow(1 - f, thinning);
     134           1 :                         candidate->addSecondary(11, Epp / (1 + z), pos, w, interactionTag);
     135             :                 }
     136           1 :                 if (random.rand() < pow(f, thinning)) {
     137           1 :                         double w = 1. / pow(f, thinning);
     138           1 :                         candidate->addSecondary(-11, Epp / (1 + z), pos, w, interactionTag);
     139             :                 }
     140             :         }
     141             :         // Update the primary particle energy.
     142             :         // This is done after adding the secondaries to correctly set the secondaries parent
     143           1 :         candidate->current.setEnergy((E - 2 * Epp) / (1. + z));
     144             : }
     145             : 
     146           1 : void EMTripletPairProduction::process(Candidate *candidate) const {
     147             :         // check if electron / positron
     148           1 :         int id = candidate->current.getId();
     149           1 :         if (abs(id) != 11)
     150             :                 return;
     151             : 
     152             :         // scale the particle energy instead of background photons
     153           1 :         double z = candidate->getRedshift();
     154           1 :         double E = (1 + z) * candidate->current.getEnergy();
     155             : 
     156             :         // check if in tabulated energy range
     157           1 :         if ((E < tabEnergy.front()) or (E > tabEnergy.back()))
     158             :                 return;
     159             : 
     160             :         // cosmological scaling of interaction distance (comoving)
     161           1 :         double scaling = pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z);
     162           1 :         double rate = scaling * interpolate(E, tabEnergy, tabRate);
     163             : 
     164             :         // run this loop at least once to limit the step size
     165           1 :         double step = candidate->getCurrentStep();
     166           1 :         Random &random = Random::instance();
     167             :         do {
     168           1 :                 double randDistance = -log(random.rand()) / rate;
     169             :                 // check for interaction; if it doesn't occur, limit next step
     170           1 :                 if (step < randDistance) { 
     171           1 :                         candidate->limitNextStep(limit / rate);
     172           1 :                         return;
     173             :                 }
     174           0 :                 performInteraction(candidate);
     175           0 :                 step -= randDistance; 
     176           0 :         } while (step > 0.);
     177             : }
     178             : 
     179           1 : void EMTripletPairProduction::setInteractionTag(std::string tag) {
     180           1 :         interactionTag = tag;
     181           1 : }
     182             : 
     183           2 : std::string EMTripletPairProduction::getInteractionTag() const {
     184           2 :         return interactionTag;
     185             : }
     186             : 
     187             : } // namespace crpropa