diff --git a/CMakeLists.txt b/CMakeLists.txt index 6088ba2fc..f66953c2d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -390,6 +390,7 @@ add_library(crpropa SHARED src/module/Acceleration.cpp src/module/Boundary.cpp src/module/BreakCondition.cpp + src/module/CandidateSplitting.cpp src/module/DiffusionSDE.cpp src/module/EMCascade.cpp src/module/EMDoublePairProduction.cpp @@ -656,6 +657,9 @@ if(ENABLE_TESTING) target_link_libraries(testAdiabaticCooling crpropa gtest gtest_main pthread ${COVERAGE_LIBS}) add_test(testAdiabaticCooling testAdiabaticCooling) + add_executable(testCandidateSplitting test/testCandidateSplitting.cpp) + target_link_libraries(testCandidateSplitting crpropa gtest gtest_main pthread ${COVERAGE_LIBS}) + add_test(testCandidateSplitting testCandidateSplitting) if(WITH_GALACTIC_LENSES) add_executable(testGalacticMagneticLens test/testMagneticLens.cpp) diff --git a/include/CRPropa.h b/include/CRPropa.h index 56cb4f291..05743bde1 100644 --- a/include/CRPropa.h +++ b/include/CRPropa.h @@ -29,6 +29,7 @@ #include "crpropa/module/Acceleration.h" #include "crpropa/module/Boundary.h" #include "crpropa/module/BreakCondition.h" +#include "crpropa/module/CandidateSplitting.h" #include "crpropa/module/DiffusionSDE.h" #include "crpropa/module/EMCascade.h" #include "crpropa/module/EMDoublePairProduction.h" diff --git a/include/crpropa/module/CandidateSplitting.h b/include/crpropa/module/CandidateSplitting.h new file mode 100644 index 000000000..514775624 --- /dev/null +++ b/include/crpropa/module/CandidateSplitting.h @@ -0,0 +1,74 @@ +#ifndef CRPROPA_CANDIDATEPLITTING_H +#define CRPROPA_CANDIDATEPLITTING_H + +#include +#include +#include +#include +#include + +#include "crpropa/Vector3.h" +#include "crpropa/Module.h" +#include "crpropa/Units.h" +#include "kiss/logger.h" + + +namespace crpropa { +/** @addtogroup Acceleration +* @{ + */ + +/** +@class CandidateSplitting +@brief Candidates are split into n copies when they gain energy and cross specified energy bins. Weights are set accordingly. + In case of Diffusice Shock Acceleration, splitting number can be adapted to expected spectral index to + compensate for the loss of particles per magnitude in energy +*/ +class CandidateSplitting: public Module { +private: + double nSplit; + double minWeight; + std::vector Ebins; + +public: + + CandidateSplitting(); + + /** Constructor + @param nSplit Number of copies candidates are split + @param Emin Minimal energy for splitting + @param Emax Maximal energy for splitting + @param nBins Number of energy bins + @param minWeight Mimimal Weight + @param log Energy bins in log + */ + CandidateSplitting(int nSplit, double Emin, double Emax, double nBins, double minWeight, bool log = false); + + /** Constructor + @param spectralIndex Expected spectral index determines splitting numbe + @param Emin Minimal energy for splitting + @param nBins Number of bins in energy, with dE(spectralIndex) it determines Emax + */ + CandidateSplitting(double spectralIndex, double Emin, int nBins); + + void process(Candidate *c) const; + + void setEnergyBins(double Emin, double Emax, double nBins, bool log); + + void setEnergyBinsDSA(double Emin, double dE, int n); + + void setNsplit(int n); + + void setMinimalWeight(double w); + + int getNsplit() const; + + double getMinimalWeight() const; + + const std::vector& getEnergyBins() const; + +}; +/** @}*/ + +} // namespace crpropa +#endif // CRPROPA_CANDIDATESPLITTING_H diff --git a/python/2_headers.i b/python/2_headers.i index 6479b7c33..dc659e03b 100644 --- a/python/2_headers.i +++ b/python/2_headers.i @@ -564,6 +564,7 @@ using namespace crpropa; // for usage of namespace in header files, necessary %include "crpropa/module/EMInverseComptonScattering.h" %include "crpropa/module/SynchrotronRadiation.h" %include "crpropa/module/AdiabaticCooling.h" +%include "crpropa/module/CandidateSplitting.h" %template(IntSet) std::set; %include "crpropa/module/Tools.h" diff --git a/src/module/CandidateSplitting.cpp b/src/module/CandidateSplitting.cpp new file mode 100644 index 000000000..6b735d567 --- /dev/null +++ b/src/module/CandidateSplitting.cpp @@ -0,0 +1,118 @@ +#include "crpropa/module/CandidateSplitting.h" + +namespace crpropa { + +CandidateSplitting::CandidateSplitting() { + // no particle splitting if EnergyBins and NSplit are not specified + setNsplit(0); + setMinimalWeight(1.); +} + +CandidateSplitting::CandidateSplitting(int nSplit, double Emin, double Emax, double nBins, double minWeight, bool log) { + setNsplit(nSplit); + setEnergyBins(Emin, Emax, nBins, log); + setMinimalWeight(minWeight); +} + +CandidateSplitting::CandidateSplitting(double spectralIndex, double Emin, int nBins) { + // to use with Diffusive Shock Acceleration + if (spectralIndex > 0){ + throw std::runtime_error( + "CandidateSplitting: spectralIndex > 0 !"); + } + + setNsplit(2); // always split in 2, calculate bins in energy for given spectrum: + double dE = pow(1. / 2, 1. / (spectralIndex + 1)); + setEnergyBinsDSA(Emin, dE, nBins); + setMinimalWeight(1. / pow(2, nBins)); +} + +void CandidateSplitting::process(Candidate *c) const { + double currE = c->current.getEnergy(); + double prevE = c->previous.getEnergy(); + + if (c->getWeight() <= minWeight){ + // minimal weight reached, no splitting + return; + } + if (currE < Ebins[0] || nSplit == 0 ){ + // current energy is smaller than first bin -> no splitting + // or, number of splits = 0 + return; + } + for (size_t i = 0; i < Ebins.size(); ++i){ + + if( prevE < Ebins[i] ){ + // previous energy is in energy bin [i-1, i] + if(currE < Ebins[i]){ + //assuming that dE greater than 0, prevE and E in same energy bin -> no splitting + return; + } + // current energy is in energy bin [i,i+1] or higher -> particle splitting for each crossing + for (size_t j = i; j < Ebins.size(); ++j ){ + + // adapted from Acceleration Module: + c->updateWeight(1. / nSplit); // * 1/n_split + + for (int i = 1; i < nSplit; i++) { + + ref_ptr new_candidate = c->clone(false); + //InteractionTag is PRIM, physically no new particles are created + new_candidate->parent = c; + new_candidate->previous.setEnergy(currE); // so that new candidate is not split again in next step! + } + if (j < Ebins.size()-1 && currE < Ebins[j+1]){ + // candidate is in energy bin [j, j+1] -> no further splitting + return; + } + } + return; + } + } +} + +void CandidateSplitting::setEnergyBins(double Emin, double Emax, double nBins, bool log) { + Ebins.resize(0); + if (Emin > Emax){ + throw std::runtime_error( + "CandidateSplitting: Emin > Emax!"); + } + double dE = (Emax-Emin)/nBins; + for (size_t i = 0; i < nBins; ++i) { + if (log == true) { + Ebins.push_back(Emin * pow(Emax / Emin, i / (nBins - 1.0))); + } else { + Ebins.push_back(Emin + i * dE); + } + } +} + +void CandidateSplitting::setEnergyBinsDSA(double Emin, double dE, int n) { + Ebins.resize(0); + for (size_t i = 1; i < n + 1; ++i) { + Ebins.push_back(Emin * pow(dE, i)); + } +} + +const std::vector& CandidateSplitting::getEnergyBins() const { + return Ebins; +} + +void CandidateSplitting::setNsplit(int n) { + nSplit = n; +} + +void CandidateSplitting::setMinimalWeight(double w) { + minWeight = w; +} + +int CandidateSplitting::getNsplit() const { + return nSplit; +} + +double CandidateSplitting::getMinimalWeight() const { + return minWeight; +} + +} // end namespace crpropa + diff --git a/test/testCandidateSplitting.cpp b/test/testCandidateSplitting.cpp new file mode 100644 index 000000000..dafdd23f6 --- /dev/null +++ b/test/testCandidateSplitting.cpp @@ -0,0 +1,90 @@ +#include "crpropa/Units.h" +#include "crpropa/Common.h" +#include "crpropa/ParticleID.h" +#include "crpropa/module/CandidateSplitting.h" + +#include "gtest/gtest.h" +#include +#include + +namespace crpropa { + +TEST(testCandidateSplitting, SimpleTest) { + int nSplit = 2; + int nBins = 4; + double minWeight = pow(1. / nSplit, 2); + double Emin = 1; // dimensionless for testing + double Emax = 10; + + CandidateSplitting split_lin(nSplit, Emin, Emax, nBins, minWeight); + double dE = (Emax - Emin) / nBins; + EXPECT_DOUBLE_EQ(split_lin.getEnergyBins()[0], Emin); + EXPECT_DOUBLE_EQ(split_lin.getEnergyBins()[1], Emin + dE); + + EXPECT_EQ(split_lin.getNsplit(), nSplit); + EXPECT_DOUBLE_EQ(split_lin.getMinimalWeight(), minWeight); + + CandidateSplitting split_log(nSplit, Emin, Emax, nBins, minWeight, true); + double dE_log = pow(Emax / Emin, 1. / (nBins - 1.0)); + EXPECT_DOUBLE_EQ(split_log.getEnergyBins()[0], Emin); + EXPECT_DOUBLE_EQ(split_log.getEnergyBins()[1], Emin * dE_log); + + double spectralIndex = -2.; + CandidateSplitting split_dsa(spectralIndex, Emin, nBins); + double dE_dsa = pow(1. / 2, 1. / (spectralIndex + 1)); + EXPECT_DOUBLE_EQ(split_dsa.getEnergyBins()[0], Emin * dE_dsa); + EXPECT_DOUBLE_EQ(split_dsa.getEnergyBins()[nBins - 1], Emin * pow(dE_dsa, nBins)); +} + + +TEST(testCandidateSplitting, CheckSplits) { + int nSplit = 2; + int nBins = 3; + double Emin = 1; // dimensionless for testing + double Emax = 10; + double minWeight = pow(1. / nSplit, 4); + + CandidateSplitting splitting(nSplit, Emin, Emax, nBins, minWeight); + Candidate c(nucleusId(1,1),0.5); + double weight = 1.0; + double serial = c.getSerialNumber(); + + splitting.process(&c); // no split + EXPECT_DOUBLE_EQ(c.getWeight(), weight); + EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial); + + c.current.setEnergy(2); + splitting.process(&c); // 1. split + weight = weight/nSplit; + EXPECT_DOUBLE_EQ(c.getWeight(), weight); + EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial + 1); + c.previous.setEnergy(2); + + c.current.setEnergy(6); + splitting.process(&c); // 2. split + weight = weight/nSplit; + EXPECT_DOUBLE_EQ(c.getWeight(), weight); + EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial + 2); + c.previous.setEnergy(6); + + c.current.setEnergy(0.5); + splitting.process(&c); // no split, cooling + EXPECT_DOUBLE_EQ(c.getWeight(), weight); + EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial + 2); + c.previous.setEnergy(0.5); + + c.current.setEnergy(6); + splitting.process(&c); // 3. & 4. split, crosses two boundaries + weight = weight/nSplit/nSplit; + EXPECT_DOUBLE_EQ(c.getWeight(), weight); + EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial + 4); + c.previous.setEnergy(6); + + c.current.setEnergy(8); + splitting.process(&c); // no split, minimal weight reached + EXPECT_DOUBLE_EQ(c.getWeight(), weight); + EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial + 4); + c.previous.setEnergy(8); +} + +} //namespace crpropa