Line data    Source code

       1             : #include "crpropa/module/ElasticScattering.h"
       2             : #include "crpropa/Units.h"
       3             : #include "crpropa/ParticleID.h"
       4             : #include "crpropa/ParticleMass.h"
       5             : #include "crpropa/Random.h"
       6             : 
       7             : #include <cmath>
       8             : #include <limits>
       9             : #include <sstream>
      10             : #include <fstream>
      11             : #include <stdexcept>
      12             : 
      13             : namespace crpropa {
      14             : 
      15             : const double ElasticScattering::lgmin = 6.;  // minimum log10(Lorentz-factor)
      16             : const double ElasticScattering::lgmax = 14.; // maximum log10(Lorentz-factor)
      17             : const size_t ElasticScattering::nlg = 201;   // number of Lorentz-factor steps
      18             : const double ElasticScattering::epsmin = log10(2 * eV) + 3;    // log10 minimum photon background energy in nucleus rest frame for elastic scattering
      19             : const double ElasticScattering::epsmax = log10(2 * eV) + 8.12; // log10 maximum photon background energy in nucleus rest frame for elastic scattering
      20             : const size_t ElasticScattering::neps = 513; // number of photon background energies in nucleus rest frame
      21             : 
      22           2 : ElasticScattering::ElasticScattering(ref_ptr<PhotonField> f) {
      23           2 :         setPhotonField(f);
      24           2 : }
      25             : 
      26           4 : void ElasticScattering::setPhotonField(ref_ptr<PhotonField> photonField) {
      27           4 :         this->photonField = photonField;
      28           4 :         std::string fname = photonField->getFieldName();
      29           4 :         setDescription("ElasticScattering: " + fname);
      30          12 :         initRate(getDataPath("ElasticScattering/rate_" + fname.substr(0,3) + ".txt"));
      31          12 :         initCDF(getDataPath("ElasticScattering/cdf_" + fname.substr(0,3) + ".txt"));
      32           4 : }
      33             : 
      34           4 : void ElasticScattering::initRate(std::string filename) {
      35           4 :         std::ifstream infile(filename.c_str());
      36           4 :         if (not infile.good())
      37           0 :                 throw std::runtime_error("ElasticScattering: could not open file " + filename);
      38             : 
      39           4 :         tabRate.clear();
      40             : 
      41         824 :         while (infile.good()) {
      42         824 :                 if (infile.peek() == '#') {
      43          16 :                         infile.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
      44          16 :                         continue;
      45             :                 }
      46             :                 double r;
      47             :                 infile >> r;
      48         808 :                 if (!infile)
      49             :                         break;
      50         804 :                 tabRate.push_back(r / Mpc);
      51             :         }
      52             : 
      53           4 :         infile.close();
      54           4 : }
      55             : 
      56           4 : void ElasticScattering::initCDF(std::string filename) {
      57           4 :         std::ifstream infile(filename.c_str());
      58           4 :         if (not infile.good())
      59           0 :                 throw std::runtime_error("ElasticScattering: could not open file " + filename);
      60             : 
      61           4 :         tabCDF.clear();
      62             :         std::string line;
      63             :         double a;
      64         820 :         while (std::getline(infile, line)) {
      65         816 :                 if (line[0] == '#')
      66          12 :                         continue;
      67             : 
      68         804 :                 std::stringstream lineStream(line);
      69             :                 lineStream >> a;
      70             : 
      71         804 :                 std::vector<double> cdf(neps);
      72      413256 :                 for (size_t i = 0; i < neps; i++) {
      73             :                         lineStream >> a;
      74      412452 :                         cdf[i] = a;
      75             :                 }
      76         804 :                 tabCDF.push_back(cdf);
      77         804 :         }
      78             : 
      79           4 :         infile.close();
      80           4 : }
      81             : 
      82           1 : void ElasticScattering::process(Candidate *candidate) const {
      83           1 :         int id = candidate->current.getId();
      84           1 :         double z = candidate->getRedshift();
      85             : 
      86           1 :         if (not isNucleus(id))
      87             :                 return;
      88             : 
      89           1 :         double lg = log10(candidate->current.getLorentzFactor() * (1 + z));
      90           1 :         if ((lg < lgmin) or (lg > lgmax))
      91             :                 return;
      92             : 
      93           1 :         int A = massNumber(id);
      94           1 :         int Z = chargeNumber(id);
      95           1 :         int N = A - Z;
      96             : 
      97           1 :         double step = candidate->getCurrentStep();
      98           7 :         while (step > 0) {
      99             : 
     100           7 :                 double rate = interpolateEquidistant(lg, lgmin, lgmax, tabRate);
     101           7 :                 rate *= Z * N / double(A);  // TRK scaling
     102           7 :                 rate *= pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z);  // cosmological scaling
     103             : 
     104             :                 // check for interaction
     105           7 :                 Random &random = Random::instance();
     106           7 :                 double randDist = -log(random.rand()) / rate;
     107           7 :                 if (step < randDist)
     108           1 :                         return;
     109             : 
     110             :                 // draw random background photon energy from CDF
     111           6 :                 size_t i = floor((lg - lgmin) / (lgmax - lgmin) * (nlg - 1)); // index of closest gamma tabulation point
     112           6 :                 size_t j = random.randBin(tabCDF[i]) - 1; // index of next lower tabulated eps value
     113             :                 double binWidth = (epsmax - epsmin) / (neps - 1); // logarithmic bin width
     114           6 :                 double eps = pow(10, epsmin + (j + random.rand()) * binWidth);
     115             : 
     116             :                 // boost to lab frame
     117           6 :                 double cosTheta = 2 * random.rand() - 1;
     118           6 :                 double E = eps * candidate->current.getLorentzFactor() * (1. - cosTheta);
     119             : 
     120           6 :                 Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
     121           6 :                 candidate->addSecondary(22, E, pos, 1., interactionTag);
     122             : 
     123             :                 // repeat with remaining step
     124           6 :                 step -= randDist;
     125             :         }
     126             : }
     127             : 
     128           0 : void ElasticScattering::setInteractionTag(std::string tag) {
     129           0 :         this -> interactionTag = tag;
     130           0 : }
     131             : 
     132           0 : std::string ElasticScattering::getInteractionTag() const {
     133           0 :         return interactionTag;
     134             : }
     135             : 
     136             : } // namespace crpropa