Line data    Source code

       1             : #include "crpropa/module/CandidateSplitting.h"
       2             : 
       3             : namespace crpropa {
       4             : 
       5           0 : CandidateSplitting::CandidateSplitting() {
       6             :         // no particle splitting if EnergyBins and NSplit are not specified
       7           0 :         setNsplit(0);
       8           0 :         setMinimalWeight(1.);
       9           0 : }
      10             : 
      11           3 : CandidateSplitting::CandidateSplitting(int nSplit, double Emin, double Emax,  double nBins, double minWeight, bool log) {
      12           3 :         setNsplit(nSplit);
      13           3 :         setEnergyBins(Emin, Emax, nBins, log);
      14           3 :         setMinimalWeight(minWeight);
      15           3 : }
      16             : 
      17           1 : CandidateSplitting::CandidateSplitting(double spectralIndex, double Emin, int nBins)  {
      18             :         // to use with Diffusive Shock Acceleration
      19           1 :         if (spectralIndex > 0){
      20           0 :                 throw std::runtime_error(
      21           0 :                                 "CandidateSplitting: spectralIndex > 0 !"); 
      22             :         }
      23             : 
      24           1 :         setNsplit(2); // always split in 2, calculate bins in energy for given spectrum:
      25           1 :         double dE = pow(1. / 2, 1. / (spectralIndex + 1)); 
      26           1 :         setEnergyBinsDSA(Emin, dE, nBins);
      27           1 :         setMinimalWeight(1. / pow(2, nBins));
      28           1 : }
      29             : 
      30           6 : void CandidateSplitting::process(Candidate *c) const {
      31           6 :         double currE = c->current.getEnergy(); 
      32           6 :         double prevE = c->previous.getEnergy();
      33             : 
      34           6 :         if (c->getWeight() <= minWeight){
      35             :                 // minimal weight reached, no splitting
      36             :                 return;
      37             :         }
      38           5 :         if (currE < Ebins[0] || nSplit == 0 ){
      39             :                 // current energy is smaller than first bin -> no splitting
      40             :                 // or, number of splits = 0
      41             :                 return;
      42             :         }
      43           4 :         for (size_t i = 0; i < Ebins.size(); ++i){
      44             :                 
      45           4 :                 if( prevE < Ebins[i] ){
      46             :                         // previous energy is in energy bin [i-1, i]
      47           3 :                         if(currE < Ebins[i]){
      48             :                                 //assuming that dE greater than 0, prevE and E in same energy bin -> no splitting
      49             :                                 return;
      50             :                         }
      51             :                         // current energy is in energy bin [i,i+1] or higher -> particle splitting for each crossing
      52           4 :                         for (size_t j = i; j < Ebins.size(); ++j ){
      53             : 
      54             :                                 // adapted from Acceleration Module:
      55           4 :                                 c->updateWeight(1. / nSplit); // * 1/n_split
      56             : 
      57           8 :                                 for (int i = 1; i < nSplit; i++) {
      58             :                                 
      59           4 :                                         ref_ptr<Candidate> new_candidate = c->clone(false);
      60           4 :                                         new_candidate->parent = c;
      61           4 :                                         new_candidate->previous.setEnergy(currE); // so that new candidate is not split again in next step!
      62             :                                         c->addSecondary(new_candidate);
      63             :                                 }
      64           4 :                                 if (j < Ebins.size()-1 && currE < Ebins[j+1]){
      65             :                                         // candidate is in energy bin [j, j+1] -> no further splitting
      66             :                                         return;
      67             :                                 }
      68             :                         }
      69             :                         return;
      70             :                 }
      71             :         }
      72             : }
      73             : 
      74           3 : void CandidateSplitting::setEnergyBins(double Emin, double Emax, double nBins, bool log) {
      75           3 :         Ebins.resize(0);
      76           3 :         if (Emin > Emax){
      77           0 :                 throw std::runtime_error(
      78           0 :                                 "CandidateSplitting: Emin > Emax!");
      79             :         }
      80           3 :         double dE = (Emax-Emin)/nBins;
      81          14 :         for (size_t i = 0; i < nBins; ++i) {
      82          11 :                 if (log == true) {
      83           4 :                         Ebins.push_back(Emin * pow(Emax / Emin, i / (nBins - 1.0)));
      84             :                 } else {
      85           7 :                         Ebins.push_back(Emin + i * dE);
      86             :                 }
      87             :         }
      88           3 : }
      89             : 
      90           1 : void CandidateSplitting::setEnergyBinsDSA(double Emin, double dE, int n) {
      91           1 :         Ebins.resize(0);
      92           5 :         for (size_t i = 1; i < n + 1; ++i) {
      93           4 :                 Ebins.push_back(Emin * pow(dE, i));
      94             :         }
      95           1 : }
      96             : 
      97           6 : const std::vector<double>& CandidateSplitting::getEnergyBins() const {
      98           6 :         return Ebins;
      99             : }
     100             : 
     101           4 : void CandidateSplitting::setNsplit(int n) {
     102           4 :         nSplit = n;
     103           4 : }
     104             : 
     105           4 : void CandidateSplitting::setMinimalWeight(double w) {
     106           4 :         minWeight = w;
     107           4 : }
     108             : 
     109           1 : int CandidateSplitting::getNsplit() const {
     110           1 :         return nSplit;
     111             : }
     112             : 
     113           1 : double CandidateSplitting::getMinimalWeight() const {
     114           1 :         return minWeight;
     115             : }
     116             : 
     117             : } // end namespace crpropa
     118             :