Line data    Source code

       1             : #include "crpropa/module/EMInverseComptonScattering.h"
       2             : #include "crpropa/Units.h"
       3             : #include "crpropa/Random.h"
       4             : #include "crpropa/Common.h"
       5             : 
       6             : #include <fstream>
       7             : #include <limits>
       8             : #include <stdexcept>
       9             : 
      10             : namespace crpropa {
      11             : 
      12             : static const double mec2 = mass_electron * c_squared;
      13             : 
      14           5 : EMInverseComptonScattering::EMInverseComptonScattering(ref_ptr<PhotonField> photonField, bool havePhotons, double thinning, double limit) {
      15           5 :         setPhotonField(photonField);
      16           5 :         setHavePhotons(havePhotons);
      17           5 :         setLimit(limit);
      18           5 :         setThinning(thinning);
      19           5 : }
      20             : 
      21          17 : void EMInverseComptonScattering::setPhotonField(ref_ptr<PhotonField> photonField) {
      22          17 :         this->photonField = photonField;
      23          17 :         std::string fname = photonField->getFieldName();
      24          17 :         setDescription("EMInverseComptonScattering: " + fname);
      25          51 :         initRate(getDataPath("EMInverseComptonScattering/rate_" + fname + ".txt"));
      26          51 :         initCumulativeRate(getDataPath("EMInverseComptonScattering/cdf_" + fname + ".txt"));
      27          17 : }
      28             : 
      29           6 : void EMInverseComptonScattering::setHavePhotons(bool havePhotons) {
      30           6 :         this->havePhotons = havePhotons;
      31           6 : }
      32             : 
      33           5 : void EMInverseComptonScattering::setLimit(double limit) {
      34           5 :         this->limit = limit;
      35           5 : }
      36             : 
      37           5 : void EMInverseComptonScattering::setThinning(double thinning) {
      38           5 :         this->thinning = thinning;
      39           5 : }
      40             : 
      41          17 : void EMInverseComptonScattering::initRate(std::string filename) {
      42          17 :         std::ifstream infile(filename.c_str());
      43             : 
      44          17 :         if (!infile.good())
      45           0 :                 throw std::runtime_error("EMInverseComptonScattering: could not open file " + filename);
      46             : 
      47             :         // clear previously loaded tables
      48          17 :         tabEnergy.clear();
      49          17 :         tabRate.clear();
      50             : 
      51        4879 :         while (infile.good()) {
      52        4862 :                 if (infile.peek() != '#') {
      53             :                         double a, b;
      54             :                         infile >> a >> b;
      55        4794 :                         if (infile) {
      56        4777 :                                 tabEnergy.push_back(pow(10, a) * eV);
      57        4777 :                                 tabRate.push_back(b / Mpc);
      58             :                         }
      59             :                 }
      60        4862 :                 infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
      61             :         }
      62          17 :         infile.close();
      63          17 : }
      64             : 
      65          17 : void EMInverseComptonScattering::initCumulativeRate(std::string filename) {
      66          17 :         std::ifstream infile(filename.c_str());
      67             : 
      68          17 :         if (!infile.good())
      69           0 :                 throw std::runtime_error("EMInverseComptonScattering: could not open file " + filename);
      70             : 
      71             :         // clear previously loaded tables
      72          17 :         tabE.clear();
      73          17 :         tabs.clear();
      74          17 :         tabCDF.clear();
      75             :         
      76             :         // skip header
      77          85 :         while (infile.peek() == '#')
      78          68 :                 infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
      79             : 
      80             :         // read s values in first line
      81             :         double a;
      82             :         infile >> a; // skip first value
      83        3000 :         while (infile.good() and (infile.peek() != '\n')) {
      84             :                 infile >> a;
      85        2983 :                 tabs.push_back(pow(10, a) * eV * eV);
      86             :         }
      87             : 
      88             :         // read all following lines: E, cdf values
      89        4794 :         while (infile.good()) {
      90             :                 infile >> a;
      91        4794 :                 if (!infile)
      92             :                         break;  // end of file
      93        4777 :                 tabE.push_back(pow(10, a) * eV);
      94             :                 std::vector<double> cdf;
      95      843000 :                 for (int i = 0; i < tabs.size(); i++) {
      96             :                         infile >> a;
      97      838223 :                         cdf.push_back(a / Mpc);
      98             :                 }
      99        4777 :                 tabCDF.push_back(cdf);
     100             :         }
     101          17 :         infile.close();
     102          17 : }
     103             : 
     104             : // Class to calculate the energy distribution of the ICS photon and to sample from it
     105             : class ICSSecondariesEnergyDistribution {
     106             :         private:
     107             :                 std::vector< std::vector<double> > data;
     108             :                 std::vector<double> s_values;
     109             :                 size_t Ns;
     110             :                 size_t Nrer;
     111             :                 double s_min;
     112             :                 double s_max;
     113             :                 double dls;
     114             : 
     115             :         public:
     116             :                 // differential cross-section, see Lee '96 (arXiv:9604098), eq. 23 for x = Ee'/Ee
     117             :                 double dSigmadE(double x, double beta) {
     118     1000000 :                         double q = ((1 - beta) / beta) * (1 - 1./x);
     119     1000000 :                         return ((1 + beta) / beta) * (x + 1./x + 2 * q + q * q);
     120             :                 }
     121             : 
     122             :                 // create the cumulative energy distribution of the up-scattered photon
     123           1 :                 ICSSecondariesEnergyDistribution() {
     124           1 :                         Ns = 1000;
     125           1 :                         Nrer = 1000;
     126           1 :                         s_min = mec2 * mec2;
     127           1 :                         s_max = 2e23 * eV * eV;
     128           1 :                         dls = (log(s_max) - log(s_min)) / Ns;
     129           1 :                         data = std::vector< std::vector<double> >(1000, std::vector<double>(1000));
     130           1 :                         std::vector<double> data_i(1000);
     131             : 
     132             :                         // tabulate s bin borders
     133           1 :                         s_values = std::vector<double>(1001);
     134        1002 :                         for (size_t i = 0; i < Ns + 1; ++i)
     135        1001 :                                 s_values[i] = s_min * exp(i*dls);
     136             : 
     137             : 
     138             :                         // for each s tabulate cumulative differential cross section
     139        1001 :                         for (size_t i = 0; i < Ns; i++) {
     140        1000 :                                 double s = s_min * exp((i+0.5) * dls);
     141        1000 :                                 double beta = (s - s_min) / (s + s_min);
     142        1000 :                                 double x0 = (1 - beta) / (1 + beta);
     143        1000 :                                 double dlx = -log(x0) / Nrer;
     144             : 
     145             :                                 // cumulative midpoint integration
     146        1000 :                                 data_i[0] = dSigmadE(x0, beta) * expm1(dlx);
     147     1000000 :                                 for (size_t j = 1; j < Nrer; j++) {
     148      999000 :                                         double x = x0 * exp((j+0.5) * dlx);
     149      999000 :                                         double dx = exp((j+1) * dlx) - exp(j * dlx);
     150      999000 :                                         data_i[j] = dSigmadE(x, beta) * dx;
     151      999000 :                                         data_i[j] += data_i[j-1];
     152             :                                 }
     153        1000 :                                 data[i] = data_i;
     154             :                         }
     155           1 :                 }
     156             : 
     157             :                 // draw random energy for the up-scattered photon Ep(Ee, s)
     158           1 :                 double sample(double Ee, double s) {
     159           1 :                         size_t idx = std::lower_bound(s_values.begin(), s_values.end(), s) - s_values.begin();
     160           1 :                         std::vector<double> s0 = data[idx];
     161           1 :                         Random &random = Random::instance();
     162           1 :                         size_t j = random.randBin(s0) + 1; // draw random bin (upper bin boundary returned)
     163           1 :                         double beta = (s - s_min) / (s + s_min);
     164           1 :                         double x0 = (1 - beta) / (1 + beta);
     165           1 :                         double dlx = -log(x0) / Nrer;
     166           1 :                         double binWidth = x0 * (exp(j * dlx) - exp((j-1) * dlx));
     167           1 :                         double Ep = (x0 * exp((j-1) * dlx) + binWidth) * Ee;
     168           2 :                         return std::min(Ee, Ep); // prevent Ep > Ee from numerical inaccuracies
     169             :                 }
     170             : };
     171             : 
     172           1 : void EMInverseComptonScattering::performInteraction(Candidate *candidate) const {
     173             :         // scale the particle energy instead of background photons
     174           1 :         double z = candidate->getRedshift();
     175           1 :         double E = candidate->current.getEnergy() * (1 + z);
     176             : 
     177           1 :         if (E < tabE.front() or E > tabE.back())
     178             :                 return;
     179             : 
     180             :         // sample the value of s
     181           1 :         Random &random = Random::instance();
     182           1 :         size_t i = closestIndex(E, tabE);
     183           1 :         size_t j = random.randBin(tabCDF[i]);
     184           1 :         double s_kin = pow(10, log10(tabs[j]) + (random.rand() - 0.5) * 0.1);
     185           1 :         double s = s_kin + mec2 * mec2;
     186             : 
     187             :         // sample electron energy after scattering
     188           1 :         static ICSSecondariesEnergyDistribution distribution;
     189           1 :         double Enew = distribution.sample(E, s);
     190             : 
     191             :         // add up-scattered photon
     192           1 :         double Esecondary = E - Enew;
     193           1 :         double f = Enew / E;
     194           1 :         if (havePhotons) {
     195           1 :                 if (random.rand() < pow(1 - f, thinning)) {
     196           1 :                         double w = 1. / pow(1 - f, thinning);
     197           1 :                         Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
     198           1 :                         candidate->addSecondary(22, Esecondary / (1 + z), pos, w, interactionTag);
     199             :                 }
     200             :         }
     201             : 
     202             :         // update the primary particle energy; do this after adding the secondary to correctly set the secondary's parent
     203           1 :         candidate->current.setEnergy(Enew / (1 + z));
     204             : }
     205             : 
     206       65201 : void EMInverseComptonScattering::process(Candidate *candidate) const {
     207             :         // check if electron / positron
     208       65201 :         int id = candidate->current.getId();
     209       65201 :         if (abs(id) != 11)
     210             :                 return;
     211             : 
     212             :         // scale the particle energy instead of background photons
     213           1 :         double z = candidate->getRedshift();
     214           1 :         double E = candidate->current.getEnergy() * (1 + z);
     215             : 
     216           1 :         if (E < tabEnergy.front() or (E > tabEnergy.back()))
     217             :                 return;
     218             : 
     219             :         // interaction rate
     220           1 :         double rate = interpolate(E, tabEnergy, tabRate);
     221           1 :         rate *= pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z);
     222             : 
     223             :         // run this loop at least once to limit the step size
     224           1 :         double step = candidate->getCurrentStep();
     225           1 :         Random &random = Random::instance();
     226             :         do {
     227           1 :                 double randDistance = -log(random.rand()) / rate;
     228             : 
     229             :                 // check for interaction; if it doesn't ocurr, limit next step
     230           1 :                 if (step < randDistance) {
     231           1 :                         candidate->limitNextStep(limit / rate);
     232           1 :                         return;
     233             :                 }
     234           0 :                 performInteraction(candidate);
     235             : 
     236             :                 // repeat with remaining step
     237           0 :                 step -= randDistance;
     238           0 :         } while (step > 0);
     239             : }
     240             : 
     241           1 : void EMInverseComptonScattering::setInteractionTag(std::string tag) {
     242           1 :         interactionTag = tag;
     243           1 : }
     244             : 
     245           2 : std::string EMInverseComptonScattering::getInteractionTag() const {
     246           2 :         return interactionTag;
     247             : }
     248             : 
     249             : } // namespace crpropa