Line data    Source code

       1              : #include "crpropa/module/EMDoublePairProduction.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            3 : EMDoublePairProduction::EMDoublePairProduction(ref_ptr<PhotonField> photonField, bool haveElectrons, double thinning, double limit) {
      12            3 :         setPhotonField(photonField);
      13            3 :         setHaveElectrons(haveElectrons);
      14            3 :         setLimit(limit);
      15            3 :         setThinning(thinning);
      16            3 : }
      17              : 
      18           15 : void EMDoublePairProduction::setPhotonField(ref_ptr<PhotonField> photonField) {
      19           15 :         this->photonField = photonField;
      20           15 :         std::string fname = photonField->getFieldName();
      21           15 :         setDescription("EMDoublePairProduction: " + fname);
      22           45 :         initRate(getDataPath("EMDoublePairProduction/rate_" + fname + ".txt"));
      23           15 : }
      24              : 
      25            4 : void EMDoublePairProduction::setHaveElectrons(bool haveElectrons) {
      26            4 :         this->haveElectrons = haveElectrons;
      27            4 : }
      28              : 
      29            3 : void EMDoublePairProduction::setLimit(double limit) {
      30            3 :         this->limit = limit;
      31            3 : }
      32              : 
      33            3 : void EMDoublePairProduction::setThinning(double thinning) {
      34            3 :         this->thinning = thinning;
      35            3 : }
      36              : 
      37           15 : void EMDoublePairProduction::initRate(std::string filename) {
      38           15 :         std::ifstream infile(filename.c_str());
      39              : 
      40           15 :         if (!infile.good())
      41            0 :                 throw std::runtime_error("EMDoublePairProduction: could not open file " + filename);
      42              : 
      43              :         // clear previously loaded interaction rates
      44              :         tabEnergy.clear();
      45              :         tabRate.clear();
      46              : 
      47         3306 :         while (infile.good()) {
      48         3291 :                 if (infile.peek() != '#') {
      49              :                         double a, b;
      50              :                         infile >> a >> b;
      51         3231 :                         if (infile) {
      52         3216 :                                 tabEnergy.push_back(pow(10, a) * eV);
      53         3216 :                                 tabRate.push_back(b / Mpc);
      54              :                         }
      55              :                 }
      56         3291 :                 infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
      57              :         }
      58           15 :         infile.close();
      59           15 : }
      60              : 
      61              : 
      62            1 : void EMDoublePairProduction::performInteraction(Candidate *candidate) const {
      63              :         // the photon is lost after the interaction
      64            1 :         candidate->setActive(false);
      65              : 
      66            1 :         if (not haveElectrons)
      67            0 :                 return;
      68              : 
      69              :         // Use assumption of Lee 96 arXiv:9604098
      70              :         // Energy is equally shared between one e+e- pair, but take mass of second e+e- pair into account.
      71              :         // This approximation has been shown to be valid within -1.5%.
      72            1 :         double z = candidate->getRedshift();
      73            1 :         double E = candidate->current.getEnergy() * (1 + z);
      74            1 :         double Ee = (E - 2 * mass_electron * c_squared) / 2;
      75              : 
      76            1 :         Random &random = Random::instance();
      77            1 :         Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
      78              : 
      79            1 :         double f = Ee / E;
      80              : 
      81            1 :                 if (random.rand() < pow(1 - f, thinning)) {
      82            1 :                         double w = 1. / pow(1 - f, thinning);
      83            1 :                         candidate->addSecondary( 11, Ee / (1 + z), pos, w, interactionTag);
      84              :                 } 
      85            1 :                 if (random.rand() < pow(f, thinning)) {
      86            1 :                         double w = 1. / pow(f, thinning);
      87            1 :                         candidate->addSecondary(-11, Ee / (1 + z), pos, w, interactionTag);
      88              :                 }
      89              : }
      90              : 
      91            1 : void EMDoublePairProduction::process(Candidate *candidate) const {
      92              :         // check if photon
      93            1 :         if (candidate->current.getId() != 22)
      94              :                 return;
      95              : 
      96              :         // scale the electron energy instead of background photons
      97            1 :         double z = candidate->getRedshift();
      98            1 :         double E = (1 + z) * candidate->current.getEnergy();
      99              : 
     100              :         // check if in tabulated energy range
     101            1 :         if (E < tabEnergy.front() or (E > tabEnergy.back()))
     102              :                 return;
     103              : 
     104              :         // interaction rate
     105            1 :         double rate = interpolate(E, tabEnergy, tabRate);
     106            1 :         rate *= pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z);
     107              : 
     108              :         // check for interaction
     109            1 :         Random &random = Random::instance();
     110            1 :         double randDistance = -log(random.rand()) / rate;
     111            1 :         double step = candidate->getCurrentStep();
     112            1 :         if (step < randDistance) {
     113            1 :                 candidate->limitNextStep(limit / rate);
     114            1 :                 return;
     115              :         } else { // after performing interaction photon ceases to exist (hence return)
     116            0 :                 performInteraction(candidate);
     117            0 :                 return;
     118              :         }
     119              : 
     120              : }
     121              : 
     122            1 : void EMDoublePairProduction::setInteractionTag(std::string tag) {
     123            1 :         interactionTag = tag;
     124            1 : }
     125              : 
     126            2 : std::string EMDoublePairProduction::getInteractionTag() const {
     127            2 :         return interactionTag;
     128              : }
     129              : 
     130              : 
     131              : } // namespace crpropa