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          15 :         tabEnergy.clear();
      45          15 :         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 (haveElectrons) {
      82           1 :                 if (random.rand() < pow(1 - f, thinning)) {
      83           1 :                         double w = 1. / pow(1 - f, thinning);
      84           1 :                         candidate->addSecondary( 11, Ee / (1 + z), pos, w, interactionTag);
      85             :                 } 
      86           1 :                 if (random.rand() < pow(f, thinning)) {
      87           1 :                         double w = 1. / pow(f, thinning);
      88           1 :                         candidate->addSecondary(-11, Ee / (1 + z), pos, w, interactionTag);
      89             :                 }
      90             :         }
      91             : }
      92             : 
      93           1 : void EMDoublePairProduction::process(Candidate *candidate) const {
      94             :         // check if photon
      95           1 :         if (candidate->current.getId() != 22)
      96             :                 return;
      97             : 
      98             :         // scale the electron energy instead of background photons
      99           1 :         double z = candidate->getRedshift();
     100           1 :         double E = (1 + z) * candidate->current.getEnergy();
     101             : 
     102             :         // check if in tabulated energy range
     103           1 :         if (E < tabEnergy.front() or (E > tabEnergy.back()))
     104             :                 return;
     105             : 
     106             :         // interaction rate
     107           1 :         double rate = interpolate(E, tabEnergy, tabRate);
     108           1 :         rate *= pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z);
     109             : 
     110             :         // check for interaction
     111           1 :         Random &random = Random::instance();
     112           1 :         double randDistance = -log(random.rand()) / rate;
     113           1 :         double step = candidate->getCurrentStep();
     114           1 :         if (step < randDistance) {
     115           1 :                 candidate->limitNextStep(limit / rate);
     116           1 :                 return;
     117             :         } else { // after performing interaction photon ceases to exist (hence return)
     118           0 :                 performInteraction(candidate);
     119           0 :                 return;
     120             :         }
     121             : 
     122             : }
     123             : 
     124           1 : void EMDoublePairProduction::setInteractionTag(std::string tag) {
     125           1 :         interactionTag = tag;
     126           1 : }
     127             : 
     128           2 : std::string EMDoublePairProduction::getInteractionTag() const {
     129           2 :         return interactionTag;
     130             : }
     131             : 
     132             : 
     133             : } // namespace crpropa