Line data    Source code

       1             : #include "crpropa/module/ElectronPairProduction.h"
       2             : #include "crpropa/Units.h"
       3             : #include "crpropa/ParticleID.h"
       4             : #include "crpropa/ParticleMass.h"
       5             : #include "crpropa/Random.h"
       6             : 
       7             : #include <fstream>
       8             : #include <limits>
       9             : #include <stdexcept>
      10             : 
      11             : namespace crpropa {
      12             : 
      13          10 : ElectronPairProduction::ElectronPairProduction(ref_ptr<PhotonField> photonField,
      14          10 :                 bool haveElectrons, double limit) {
      15          10 :         this->haveElectrons = haveElectrons;
      16          10 :         this->limit = limit;
      17          10 :         setPhotonField(photonField);
      18          10 : }
      19             : 
      20          19 : void ElectronPairProduction::setPhotonField(ref_ptr<PhotonField> photonField) {
      21          19 :         this->photonField = photonField;
      22          19 :         std::string fname = photonField->getFieldName();
      23          19 :         setDescription("ElectronPairProduction: " + fname);
      24          57 :         initRate(getDataPath("ElectronPairProduction/lossrate_" + fname + ".txt"));
      25          19 :         if (haveElectrons) { // Load secondary spectra only if electrons should be produced
      26           0 :                 initSpectrum(getDataPath("ElectronPairProduction/spectrum_" + fname.substr(0,3) + ".txt"));
      27             :         }
      28          19 : }
      29             : 
      30           1 : void ElectronPairProduction::setHaveElectrons(bool haveElectrons) {
      31           1 :         this->haveElectrons = haveElectrons;
      32           1 :         if (haveElectrons) { // Load secondary spectra in case haveElectrons was changed to true
      33           1 :                 std::string fname = photonField->getFieldName();
      34           3 :                 initSpectrum(getDataPath("ElectronPairProduction/spectrum_" + fname.substr(0,3) + ".txt"));
      35             :         }
      36           1 : }
      37             : 
      38           0 : void ElectronPairProduction::setLimit(double limit) {
      39           0 :         this->limit = limit;
      40           0 : }
      41             : 
      42          19 : void ElectronPairProduction::initRate(std::string filename) {
      43          19 :         std::ifstream infile(filename.c_str());
      44             : 
      45          19 :         if (!infile.good())
      46           0 :                 throw std::runtime_error("ElectronPairProduction: could not open file " + filename);
      47             : 
      48             :         // clear previously loaded interaction rates
      49          19 :         tabLorentzFactor.clear();
      50          19 :         tabLossRate.clear();
      51             : 
      52        2853 :         while (infile.good()) {
      53        2834 :                 if (infile.peek() != '#') {
      54             :                         double a, b;
      55             :                         infile >> a >> b;
      56        2777 :                         if (infile) {
      57        2758 :                                 tabLorentzFactor.push_back(pow(10, a));
      58        2758 :                                 tabLossRate.push_back(b / Mpc);
      59             :                         }
      60             :                 }
      61        2834 :                 infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
      62             :         }
      63          19 :         infile.close();
      64          19 : }
      65             : 
      66           1 : void ElectronPairProduction::initSpectrum(std::string filename) {
      67           1 :         std::ifstream infile(filename.c_str());
      68           1 :         if (!infile.good())
      69           0 :                 throw std::runtime_error("ElectronPairProduction: could not open file " + filename);
      70             : 
      71             :         double dNdE;
      72           1 :         tabSpectrum.resize(70);
      73          71 :         for (size_t i = 0; i < 70; i++) {
      74          70 :                 tabSpectrum[i].resize(170);
      75       11970 :                 for (size_t j = 0; j < 170; j++) {
      76             :                         infile >> dNdE;
      77       11900 :                         tabSpectrum[i][j] = dNdE * pow(10, (7 + 0.1 * j)); // read electron distribution pdf(Ee) ~ dN/dEe * Ee
      78             :                 }
      79       11900 :                 for (size_t j = 1; j < 170; j++) {
      80       11830 :                         tabSpectrum[i][j] += tabSpectrum[i][j - 1]; // cdf(Ee), unnormalized
      81             :                 }
      82             :         }
      83           1 :         infile.close();
      84           1 : }
      85             : 
      86       65362 : double ElectronPairProduction::lossLength(int id, double lf, double z) const {
      87       65362 :         double Z = chargeNumber(id);
      88       65362 :         if (Z == 0)
      89             :                 return std::numeric_limits<double>::max(); // no pair production on uncharged particles
      90             : 
      91       64916 :         lf *= (1 + z);
      92       64916 :         if (lf < tabLorentzFactor.front())
      93             :                 return std::numeric_limits<double>::max(); // below energy threshold
      94             : 
      95             :         double rate;
      96       64893 :         if (lf < tabLorentzFactor.back())
      97       64893 :                 rate = interpolate(lf, tabLorentzFactor, tabLossRate); // interpolation
      98             :         else
      99           0 :                 rate = tabLossRate.back() * pow(lf / tabLorentzFactor.back(), -0.6); // extrapolation
     100             : 
     101       64893 :         double A = nuclearMass(id) / mass_proton; // more accurate than massNumber(Id)
     102       64893 :         rate *= Z * Z / A * pow_integer<3>(1 + z) * photonField->getRedshiftScaling(z);
     103       64893 :         return 1. / rate;
     104             : }
     105             : 
     106       65363 : void ElectronPairProduction::process(Candidate *c) const {
     107       65363 :         int id = c->current.getId();
     108       65363 :         if (not (isNucleus(id)))
     109             :                 return; // only nuclei
     110             : 
     111       65362 :         double lf = c->current.getLorentzFactor();
     112       65362 :         double z = c->getRedshift();
     113       65362 :         double losslen = lossLength(id, lf, z);  // energy loss length
     114       65362 :         if (losslen >= std::numeric_limits<double>::max())
     115             :                 return;
     116             : 
     117       64893 :         double step = c->getCurrentStep() / (1 + z); // step size in local frame
     118       64893 :         double loss = step / losslen;  // relative energy loss
     119             : 
     120       64893 :         if (haveElectrons) {
     121           1 :                 double dE = c->current.getEnergy() * loss;  // energy loss
     122           1 :                 int i = round((log10(lf) - 6.05) * 10);  // find closest cdf(Ee|log10(gamma))
     123           1 :                 i = std::min(std::max(i, 0), 69);
     124           1 :                 Random &random = Random::instance();
     125             : 
     126             :                 // draw pairs as long as their energy is smaller than the pair production energy loss
     127        7131 :                 while (dE > 0) {
     128        7131 :                         size_t j = random.randBin(tabSpectrum[i]);
     129        7131 :                         double Ee = pow(10, 6.95 + (j + random.rand()) * 0.1) * eV;
     130        7131 :                         double Epair = 2 * Ee; // NOTE: electron and positron in general don't have same lab frame energy, but averaged over many draws the result is consistent
     131             :                         // if the remaining energy is not sufficient check for random accepting
     132        7131 :                         if (Epair > dE)
     133           1 :                                 if (random.rand() > (dE / Epair))
     134             :                                         break; // not accepted
     135             : 
     136             :                         // create pair and repeat with remaining energy
     137        7130 :                         dE -= Epair;
     138        7130 :                         Vector3d pos = random.randomInterpolatedPosition(c->previous.getPosition(), c->current.getPosition());
     139        7130 :                         c->addSecondary( 11, Ee, pos, 1., interactionTag);
     140        7130 :                         c->addSecondary(-11, Ee, pos, 1., interactionTag);
     141             :                 }
     142             :         }
     143             : 
     144       64893 :         c->current.setLorentzFactor(lf * (1 - loss));
     145       64893 :         c->limitNextStep(limit * losslen);
     146             : }
     147             : 
     148           1 : void ElectronPairProduction::setInteractionTag(std::string tag) {
     149           1 :         interactionTag = tag;
     150           1 : }
     151             : 
     152           2 : std::string ElectronPairProduction::getInteractionTag() const {
     153           2 :         return interactionTag;
     154             : }
     155             : 
     156             : 
     157             : } // namespace crpropa