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