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            8 :         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              :         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              :         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