Line data    Source code

       1              : #include "crpropa/module/EMDoublePairProduction.h"
       2              : #include "crpropa/Units.h"
       3              : #include "crpropa/Random.h"
       4              : #include "crpropa/Common.h"
       5              : 
       6              : #include <vector>
       7              : #include <cmath>
       8              : 
       9              : namespace crpropa {
      10              : 
      11            3 : EMDoublePairProduction::EMDoublePairProduction(ref_ptr<PhotonField> photonField, bool haveElectrons, double thinning, double limit, ref_ptr<Surface> surface) {
      12            6 :         setSurface(surface);
      13            3 :         setPhotonField(photonField);
      14            3 :         setHaveElectrons(haveElectrons);
      15            3 :         setLimit(limit);
      16            3 :         setThinning(thinning);
      17            3 : }
      18              : 
      19           15 : void EMDoublePairProduction::setPhotonField(ref_ptr<PhotonField> photonField) {
      20              :         
      21           15 :         this->photonField = photonField;
      22           15 :         std::string fname = photonField->getFieldName();
      23           30 :         setDescription("EMDoublePairProduction: " + fname);
      24              :         
      25              :         // choose the right interaction rates for the used photon field
      26           15 :         if (!this->photonField->hasPositionDependence()) {
      27              :                 this->interactionRates = new InteractionRatesHomogeneous(
      28           45 :                         getDataPath("EMDoublePairProduction/rate_" + fname + ".txt")
      29           45 :                 );
      30              :         } else {
      31              :                 this->interactionRates = new InteractionRatesPositionDependent(
      32            0 :                         getDataPath("EMDoublePairProduction/"+fname+"/Rate/"),
      33              :                         "",
      34              :                         this->surface
      35            0 :                 );
      36              :         }
      37           15 : }
      38              : 
      39            4 : void EMDoublePairProduction::setHaveElectrons(bool haveElectrons) {
      40            4 :         this->haveElectrons = haveElectrons;
      41            4 : }
      42              : 
      43            3 : void EMDoublePairProduction::setLimit(double limit) {
      44            3 :         this->limit = limit;
      45            3 : }
      46              : 
      47            3 : void EMDoublePairProduction::setThinning(double thinning) {
      48            3 :         this->thinning = thinning;
      49            3 : }
      50              : 
      51            3 : void EMDoublePairProduction::setSurface(ref_ptr<Surface> surface) {
      52            3 :         this->surface = surface;
      53            3 : }
      54              : 
      55            0 : ref_ptr<Surface> EMDoublePairProduction::getSurface() const {
      56            0 :         return this->surface;
      57              : }
      58              : 
      59            0 : void EMDoublePairProduction::setInteractionRates(ref_ptr<InteractionRates> intRates) {
      60            0 :         this->interactionRates = intRates;
      61            0 : }
      62              : 
      63            0 : ref_ptr<InteractionRates> EMDoublePairProduction::getInteractionRates() const {
      64            0 :         return this->interactionRates;
      65              : }
      66              : 
      67            0 : void EMDoublePairProduction::initRate(std::string path) {
      68            0 :         this->interactionRates->initRate(path);
      69            0 : }
      70              : 
      71            1 : double EMDoublePairProduction::getRate(double E, const Vector3d &position, double z) const {
      72            1 :         return this->interactionRates->getProcessRate(E, position) * pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z);
      73              : }
      74              : 
      75            1 : void EMDoublePairProduction::performInteraction(Candidate *candidate) const {
      76              : 
      77              :         // Use assumption of Lee 96 arXiv:9604098
      78              :         // Energy is equally shared between one e+e- pair, but take mass of second e+e- pair into account.
      79              :         // This approximation has been shown to be valid within -1.5%.
      80            1 :         double z = candidate->getRedshift();
      81            1 :         double E = candidate->current.getEnergy() * (1 + z);
      82            1 :         double Ee = (E - 2 * mass_electron * c_squared) / 2;
      83              : 
      84              :         // the photon is lost after the interaction
      85            1 :         candidate->setActive(false);
      86              : 
      87            1 :         if (not haveElectrons)
      88            0 :                 return;
      89              :         
      90            1 :         Random &random = Random::instance();
      91            1 :         Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
      92              : 
      93            1 :         double f = Ee / E;
      94              : 
      95            1 :                 if (random.rand() < pow(1 - f, thinning)) {
      96            1 :                         double w = 1. / pow(1 - f, thinning);
      97            1 :                         candidate->addSecondary( 11, Ee / (1 + z), pos, w, interactionTag);
      98              :                 } 
      99            1 :                 if (random.rand() < pow(f, thinning)) {
     100            1 :                         double w = 1. / pow(f, thinning);
     101            1 :                         candidate->addSecondary(-11, Ee / (1 + z), pos, w, interactionTag);
     102              :                 }
     103              : }
     104              : 
     105            1 : void EMDoublePairProduction::process(Candidate *candidate) const {
     106              :         
     107              :         // check if photon
     108            1 :         if (candidate->current.getId() != 22)
     109              :                 return;
     110              : 
     111              :         // scale the electron energy instead of background photons
     112            1 :         double z = candidate->getRedshift();
     113            1 :         double E = (1 + z) * candidate->current.getEnergy();
     114            1 :         Vector3d position = candidate->current.getPosition();
     115              : 
     116              :         // interaction rate
     117            1 :         double rate = getRate(E, position, z);
     118            1 :         if (rate < 0)
     119              :                 return;
     120              :                 
     121              :         // check for interaction
     122            1 :         Random &random = Random::instance();
     123            1 :         double randDistance = -log(random.rand()) / rate;
     124            1 :         double step = candidate->getCurrentStep();
     125            1 :         if (step < randDistance) {
     126            1 :                 candidate->limitNextStep(limit / rate);
     127            1 :                 return;
     128              :         } else { // after performing interaction photon ceases to exist (hence return)
     129            0 :                 performInteraction(candidate);
     130            0 :                 return;
     131              :         }
     132              : 
     133              : }
     134              : 
     135            1 : void EMDoublePairProduction::setInteractionTag(std::string tag) {
     136            1 :         interactionTag = tag;
     137            1 : }
     138              : 
     139            2 : std::string EMDoublePairProduction::getInteractionTag() const {
     140            2 :         return interactionTag;
     141              : }
     142              : 
     143              : 
     144              : } // namespace crpropa