Line data    Source code

       1             : #include "crpropa/module/SynchrotronRadiation.h"
       2             : #include "crpropa/Units.h"
       3             : #include "crpropa/Random.h"
       4             : 
       5             : #include <fstream>
       6             : #include <limits>
       7             : #include <stdexcept>
       8             : 
       9             : namespace crpropa {
      10             : 
      11           0 : SynchrotronRadiation::SynchrotronRadiation(ref_ptr<MagneticField> field, bool havePhotons, double thinning, int nSamples, double limit) {
      12           0 :         setField(field);
      13           0 :         setBrms(0);
      14           0 :         initSpectrum();
      15           0 :         setHavePhotons(havePhotons);
      16           0 :         setLimit(limit);
      17           0 :         setSecondaryThreshold(1e6 * eV);
      18           0 :         setMaximumSamples(nSamples);
      19           0 : }
      20             : 
      21           1 : SynchrotronRadiation::SynchrotronRadiation(double Brms, bool havePhotons, double thinning, int nSamples, double limit) {
      22           1 :         setBrms(Brms);
      23           1 :         initSpectrum();
      24           1 :         setHavePhotons(havePhotons);
      25           1 :         setLimit(limit);
      26           1 :         setSecondaryThreshold(1e6 * eV);
      27           1 :         setMaximumSamples(nSamples);
      28           1 : }
      29             : 
      30           0 : void SynchrotronRadiation::setField(ref_ptr<MagneticField> f) {
      31           0 :         this->field = f;
      32           0 : }
      33             : 
      34           0 : ref_ptr<MagneticField> SynchrotronRadiation::getField() {
      35           0 :         return field;
      36             : }
      37             : 
      38           1 : void SynchrotronRadiation::setBrms(double Brms) {
      39           1 :         this->Brms = Brms;
      40           1 : }
      41             : 
      42           0 : double SynchrotronRadiation::getBrms() {
      43           0 :         return Brms;
      44             : }
      45             : 
      46           1 : void SynchrotronRadiation::setHavePhotons(bool havePhotons) {
      47           1 :         this->havePhotons = havePhotons;
      48           1 : }
      49             : 
      50           0 : bool SynchrotronRadiation::getHavePhotons() {
      51           0 :         return havePhotons;
      52             : }
      53             : 
      54           0 : void SynchrotronRadiation::setThinning(double thinning) {
      55           0 :         this->thinning = thinning;
      56           0 : }
      57             : 
      58           0 : double SynchrotronRadiation::getThinning() {
      59           0 :         return thinning;
      60             : }
      61             : 
      62           1 : void SynchrotronRadiation::setLimit(double limit) {
      63           1 :         this->limit = limit;
      64           1 : }
      65             : 
      66           0 : double SynchrotronRadiation::getLimit() {
      67           0 :         return limit;
      68             : }
      69             : 
      70           1 : void SynchrotronRadiation::setMaximumSamples(int nmax) {
      71           1 :         maximumSamples = nmax;
      72           1 : }
      73             : 
      74           0 : int SynchrotronRadiation::getMaximumSamples() {
      75           0 :         return maximumSamples;
      76             : }
      77             : 
      78           1 : void SynchrotronRadiation::setSecondaryThreshold(double threshold) {
      79           1 :         secondaryThreshold = threshold;
      80           1 : }
      81             : 
      82           0 : double SynchrotronRadiation::getSecondaryThreshold() const {
      83           0 :         return secondaryThreshold;
      84             : }
      85             : 
      86           1 : void SynchrotronRadiation::initSpectrum() {
      87           2 :         std::string filename = getDataPath("Synchrotron/spectrum.txt");
      88           1 :         std::ifstream infile(filename.c_str());
      89             : 
      90           1 :         if (!infile.good())
      91           0 :                 throw std::runtime_error("SynchrotronRadiation: could not open file " + filename);
      92             : 
      93             :         // clear previously loaded interaction rates
      94           1 :         tabx.clear();
      95           1 :         tabCDF.clear();
      96             : 
      97        1407 :         while (infile.good()) {
      98        1406 :                 if (infile.peek() != '#') {
      99             :                         double a, b;
     100             :                         infile >> a >> b;
     101        1402 :                         if (infile) {
     102        1401 :                                 tabx.push_back(pow(10, a));
     103        1401 :                                 tabCDF.push_back(b);
     104             :                         }
     105             :                 }
     106        1406 :                 infile.ignore(std::numeric_limits < std::streamsize > ::max(), '\n');
     107             :         }
     108           1 :         infile.close();
     109           2 : }
     110             : 
     111           1 : void SynchrotronRadiation::process(Candidate *candidate) const {
     112           1 :         double charge = fabs(candidate->current.getCharge());
     113           1 :         if (charge == 0)
     114           0 :                 return; // only charged particles
     115             : 
     116             :         // calculate gyroradius, evaluated at the current position
     117           1 :         double z = candidate->getRedshift();
     118             :         double B;
     119           1 :         if (field.valid()) {
     120           0 :                 Vector3d Bvec = field->getField(candidate->current.getPosition(), z);
     121           0 :                 B = Bvec.cross(candidate->current.getDirection()).getR();
     122             :         } else {
     123           1 :                 B = sqrt(2. / 3) * Brms; // average perpendicular field component
     124             :         }
     125           1 :         B *= pow(1 + z, 2); // cosmological scaling
     126           1 :         double Rg = candidate->current.getMomentum().getR() / charge / B;
     127             : 
     128             :         // calculate energy loss
     129           1 :         double lf = candidate->current.getLorentzFactor();
     130           1 :         double dEdx = 1. / 6 / M_PI / epsilon0 * pow(lf * lf - 1, 2) * pow(charge / Rg, 2); // Jackson p. 770 (14.31)
     131           1 :         double step = candidate->getCurrentStep() / (1 + z); // step size in local frame
     132           1 :         double dE = step * dEdx;
     133             : 
     134             :         // apply energy loss and limit next step
     135           1 :         double E = candidate->current.getEnergy();
     136           1 :         candidate->current.setEnergy(E - dE);
     137           1 :         candidate->limitNextStep(limit * E / dEdx);
     138             : 
     139             :         // optionally add secondary photons
     140           1 :         if (not(havePhotons))
     141             :                 return;
     142             : 
     143             :         // check if photons with energies > 14 * Ecrit are possible
     144           1 :         double Ecrit = 3. / 4 * h_planck / M_PI * c_light * pow(lf, 3) / Rg;
     145           1 :         if (14 * Ecrit < secondaryThreshold)
     146             :                 return;
     147             : 
     148             :         // draw photons up to the total energy loss
     149             :         // if maximumSamples is reached before that, compensate the total energy afterwards
     150           1 :         Random &random = Random::instance();
     151             :         double dE0 = dE;
     152             :         std::vector<double> energies;
     153             :         int counter = 0;
     154     3623058 :         while (dE > 0) {
     155             :                 // draw random value between 0 and maximum of corresponding cdf
     156             :                 // choose bin of s where cdf(x) = cdf_rand -> x_rand
     157     3623058 :                 size_t i = random.randBin(tabCDF); // draw random bin (upper bin boundary returned)
     158     3623058 :                 double binWidth = (tabx[i] - tabx[i-1]);
     159     3623058 :                 double x = tabx[i-1] + random.rand() * binWidth; // draw random x uniformly distributed in bin
     160     3623058 :                 double Ephoton = x * Ecrit;
     161             : 
     162             :                 // if the remaining energy is not sufficient check for random accepting
     163     3623058 :                 if (Ephoton > dE) {
     164           1 :                         if (random.rand() > (dE / Ephoton))
     165             :                                 break; // not accepted
     166             :                 }
     167             : 
     168             :                 // only activate the "per-step" sampling if maximumSamples is explicitly set.
     169     3623057 :                 if (maximumSamples > 0) {
     170           0 :                         if (counter >= maximumSamples) 
     171             :                                 break;                  
     172             :                 }
     173             : 
     174             :                 // store energies in array
     175     3623057 :                 energies.push_back(Ephoton);
     176             : 
     177             :                 // energy loss
     178     3623057 :                 dE -= Ephoton;
     179             : 
     180             :                 // counter for sampling break condition;
     181     3623057 :                 counter++;
     182             :         }
     183             : 
     184             :         // while loop before gave total energy which is just a fraction of the required
     185             :         double w1 = 1;
     186           1 :         if (maximumSamples > 0 && dE > 0)
     187           0 :                 w1 = 1. / (1. - dE / dE0); 
     188             : 
     189             :         // loop over sampled photons and attribute weights accordingly
     190     3623058 :         for (int i = 0; i < energies.size(); i++) {
     191     3623057 :                 double Ephoton = energies[i];
     192     3623057 :                 double f = Ephoton / (E - dE0);
     193     3623057 :                 double w = w1 / pow(f, thinning);
     194             : 
     195             :                 // thinning procedure: accepts only a few random secondaries
     196     3623057 :                 if (random.rand() < pow(f, thinning)) {
     197     3623057 :                         Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
     198     3623057 :                         if (Ephoton > secondaryThreshold) // create only photons with energies above threshold
     199     3311363 :                                 candidate->addSecondary(22, Ephoton, pos, w, interactionTag);
     200             :                 }
     201             :         }
     202             : }
     203             : 
     204           0 : std::string SynchrotronRadiation::getDescription() const {
     205           0 :         std::stringstream s;
     206           0 :         s << "Synchrotron radiation";
     207           0 :         if (field.valid())
     208           0 :                 s << " for specified magnetic field";
     209             :         else
     210           0 :                 s << " for Brms = " << Brms / nG << " nG";
     211           0 :         if (havePhotons)
     212           0 :                 s << ", synchrotron photons E > " << secondaryThreshold / eV << " eV";
     213             :         else
     214           0 :                 s << ", no synchrotron photons";
     215           0 :         if (maximumSamples > 0)
     216           0 :                 s << "maximum number of photon samples: " << maximumSamples;
     217           0 :         if (thinning > 0)
     218           0 :                 s << "thinning parameter: " << thinning; 
     219           0 :         return s.str();
     220           0 : }
     221             : 
     222           1 : void SynchrotronRadiation::setInteractionTag(std::string tag) {
     223           1 :         interactionTag = tag;
     224           1 : }
     225             : 
     226           2 : std::string SynchrotronRadiation::getInteractionTag() const {
     227           2 :         return interactionTag;
     228             : }
     229             : 
     230             : } // namespace crpropa