From f784d31eb9af5ef10b5a1449f4dfeb7222de0ea1 Mon Sep 17 00:00:00 2001 From: Sophie Aerdker Date: Wed, 4 Oct 2023 13:42:53 +0200 Subject: [PATCH 1/5] Candidate Splitting module --- CMakeLists.txt | 5 + include/CRPropa.h | 1 + include/crpropa/module/CandidateSplitting.h | 80 ++++++++++ python/2_headers.i | 1 + src/module/CandidateSplitting.cpp | 154 ++++++++++++++++++++ test/testCandidateSplitting.cpp | 90 ++++++++++++ 6 files changed, 331 insertions(+) create mode 100644 include/crpropa/module/CandidateSplitting.h create mode 100644 src/module/CandidateSplitting.cpp create mode 100644 test/testCandidateSplitting.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 6088ba2fc..327ef54ea 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,10 @@ 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..7d4479568 --- /dev/null +++ b/include/crpropa/module/CandidateSplitting.h @@ -0,0 +1,80 @@ +#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 { + +/** +@class CandidateSplitting +@brief Candidates are split into n copies when they 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 n_split; + double minWeight; + std::vector Ebins; + +public: + + CandidateSplitting(); + + CandidateSplitting(int n_split, double Emin, double Emax, double n_bins, double minWeight); + /** Constructor + @param n_split Number of copies candidates are split + @param Emin Minimal energy for splitting + @param Emax Maximal energy for splitting + @param minWeight Mimimal Weight + @param n_bins Number of energy bins + */ + CandidateSplitting(int n_split, double Emin, double Emax, double n_bins, double minWeight, bool log); + /** Constructor + @param n_split Number of copies candidates are split + @param Emin Minimal energy for splitting + @param Emax Maximal energy for splitting + @param n_bins Number of energy bins + @param minWeight Mimimal Weight + @param log Energy bins in log + */ + + CandidateSplitting(double SpectralIndex, double Emin, int nBins); + /** Constructor + @param SpectralIndex Absolute value of expected spectral index determines splitting number + @param Emin Minimal energy for splitting + @param nBins Number of bins in energy, with dE(spectralIndex) it determines Emax + */ + + void process(Candidate *c) const; + + void setEnergyBins(double Emin, double Emax, double n_bins, 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_PARTICLESPLITTING_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..700e1eb32 --- /dev/null +++ b/src/module/CandidateSplitting.cpp @@ -0,0 +1,154 @@ +#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 n_split, double Emin, double Emax, double n_bins, double minWeight) { + // linear energy bins: + setNsplit(n_split); + setEnergyBins(Emin, Emax, n_bins, false); + setMinimalWeight(minWeight); +} + +CandidateSplitting::CandidateSplitting(int n_split, double Emin, double Emax, double n_bins, double minWeight, bool log) { + setNsplit(n_split); + setEnergyBins(Emin, Emax, n_bins, 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] || n_split == 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. / n_split); // * 1/n_split + + for (int i = 1; i < n_split; i++) { + + ref_ptr new_candidate = c->clone(false); + new_candidate->parent = c; + uint64_t snr = Candidate::getNextSerialNumber(); + new_candidate->setSerialNumber(snr); + new_candidate->previous.setEnergy(currE); // so that new candidate is not split again in next step! + //InteractionTag is PRIM, physically no new particles are created + c->addSecondary(new_candidate); + Candidate::setNextSerialNumber(snr + 1); + } + + + 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 n_bins, bool log) { + + Ebins.resize(0); + + if (Emin > Emax){ + throw std::runtime_error( + "CandidateSplitting: Emin > Emax!"); + } + + double dE = (Emax-Emin)/n_bins; + + for (size_t i = 0; i < n_bins; ++i) { + if (log == true) { + Ebins.push_back(Emin * pow(Emax / Emin, i / (n_bins - 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) { + + n_split = n; +} + +void CandidateSplitting::setMinimalWeight(double w) { + + minWeight = w; +} + +int CandidateSplitting::getNsplit() const { + + return n_split; +} + +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..ab5f24115 --- /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 n_split = 2; + int n_bins = 4; + double minWeight = pow(1./n_split, 2); + double Emin = 1*GeV; + double Emax = 10*GeV; + + CandidateSplitting split_lin(n_split, Emin, Emax, n_bins, minWeight); + double dE = (Emax-Emin)/n_bins; + EXPECT_DOUBLE_EQ(split_lin.getEnergyBins()[0], Emin); + EXPECT_DOUBLE_EQ(split_lin.getEnergyBins()[1], Emin+dE); + + EXPECT_EQ(split_lin.getNsplit(), n_split); + EXPECT_DOUBLE_EQ(split_lin.getMinimalWeight(), minWeight); + + CandidateSplitting split_log(n_split, Emin, Emax, n_bins, minWeight, true); + double dE_log = pow(Emax / Emin, 1./(n_bins - 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, n_bins); + 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()[n_bins-1], Emin * pow(dE_dsa, n_bins)); + + +} + + +TEST(testCandidateSplitting, CheckSplits) { + + int n_split = 2; + int n_bins = 3; + double Emin = 1*GeV; + double Emax = 10*GeV; + double minWeight = pow(1./n_split, 4); + + CandidateSplitting splitting(n_split, Emin, Emax, n_bins, minWeight); + Candidate c(nucleusId(1,1),0.5*GeV); + double weight = 1.0; + + splitting.process(&c); // no split + EXPECT_DOUBLE_EQ(c.getWeight(), weight); + + c.current.setEnergy(2*GeV); + splitting.process(&c); // 1. split + weight = weight/n_split; + EXPECT_DOUBLE_EQ(c.getWeight(), weight); + c.previous.setEnergy(2*GeV); + + c.current.setEnergy(6*GeV); + splitting.process(&c); // 2. split + weight = weight/n_split; + EXPECT_DOUBLE_EQ(c.getWeight(), weight); + c.previous.setEnergy(6*GeV); + + c.current.setEnergy(0.5*GeV); + splitting.process(&c); // no split, cooling + EXPECT_DOUBLE_EQ(c.getWeight(), weight); + c.previous.setEnergy(0.5*GeV); + + c.current.setEnergy(6*GeV); + splitting.process(&c); // 3. & 4. split, crosses two boundaries + weight = weight/n_split/n_split; + EXPECT_DOUBLE_EQ(c.getWeight(), weight); + c.previous.setEnergy(6*GeV); + + c.current.setEnergy(8*GeV); + splitting.process(&c); // no split, minimal weight reached + EXPECT_DOUBLE_EQ(c.getWeight(), weight); + c.previous.setEnergy(8*GeV); + +} + + +} //namespace crpropa From f93baa0f0f21b2ed135bba991776499c613d72e5 Mon Sep 17 00:00:00 2001 From: Sophie Aerdker Date: Fri, 6 Oct 2023 10:27:56 +0200 Subject: [PATCH 2/5] addressing review requests by JD --- CMakeLists.txt | 1 - include/crpropa/module/CandidateSplitting.h | 28 ++++----- src/module/CandidateSplitting.cpp | 63 +++++++-------------- test/testCandidateSplitting.cpp | 43 ++++++-------- 4 files changed, 51 insertions(+), 84 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 327ef54ea..f66953c2d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -661,7 +661,6 @@ if(ENABLE_TESTING) 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) target_link_libraries(testGalacticMagneticLens crpropa gtest gtest_main pthread ${COVERAGE_LIBS}) diff --git a/include/crpropa/module/CandidateSplitting.h b/include/crpropa/module/CandidateSplitting.h index 7d4479568..70ec97baf 100644 --- a/include/crpropa/module/CandidateSplitting.h +++ b/include/crpropa/module/CandidateSplitting.h @@ -17,14 +17,13 @@ namespace crpropa { /** @class CandidateSplitting -@brief Candidates are split into n copies when they cross specified energy bins. Weights are set accordingly. +@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 n_split; + double nSplit; double minWeight; std::vector Ebins; @@ -32,34 +31,35 @@ class CandidateSplitting: public Module { CandidateSplitting(); - CandidateSplitting(int n_split, double Emin, double Emax, double n_bins, double minWeight); /** Constructor - @param n_split Number of copies candidates are split + @param nSplit Number of copies candidates are split @param Emin Minimal energy for splitting @param Emax Maximal energy for splitting @param minWeight Mimimal Weight - @param n_bins Number of energy bins + @param nBins Number of energy bins */ - CandidateSplitting(int n_split, double Emin, double Emax, double n_bins, double minWeight, bool log); + CandidateSplitting(int nSplit, double Emin, double Emax, double nBins, double minWeight); + /** Constructor - @param n_split Number of copies candidates are split + @param nSplit Number of copies candidates are split @param Emin Minimal energy for splitting @param Emax Maximal energy for splitting - @param n_bins Number of energy bins + @param nBins Number of energy bins @param minWeight Mimimal Weight @param log Energy bins in log */ - - CandidateSplitting(double SpectralIndex, double Emin, int nBins); + CandidateSplitting(int nSplit, double Emin, double Emax, double nBins, double minWeight, bool log); + /** Constructor - @param SpectralIndex Absolute value of expected spectral index determines splitting number + @param spectralIndex Absolute value of expected spectral index determines splitting number @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 n_bins, bool log); + void setEnergyBins(double Emin, double Emax, double nBins, bool log); void setEnergyBinsDSA(double Emin, double dE, int n); @@ -77,4 +77,4 @@ class CandidateSplitting: public Module { /** @}*/ } // namespace crpropa -#endif // CRPROPA_PARTICLESPLITTING_H +#endif // CRPROPA_CANDIDATESPLITTING_H diff --git a/src/module/CandidateSplitting.cpp b/src/module/CandidateSplitting.cpp index 700e1eb32..8e89ed838 100644 --- a/src/module/CandidateSplitting.cpp +++ b/src/module/CandidateSplitting.cpp @@ -8,37 +8,33 @@ CandidateSplitting::CandidateSplitting() { setMinimalWeight(1.); } -CandidateSplitting::CandidateSplitting(int n_split, double Emin, double Emax, double n_bins, double minWeight) { +CandidateSplitting::CandidateSplitting(int nSplit, double Emin, double Emax, double nBins, double minWeight) { // linear energy bins: - setNsplit(n_split); - setEnergyBins(Emin, Emax, n_bins, false); + setNsplit(nSplit); + setEnergyBins(Emin, Emax, nBins, false); setMinimalWeight(minWeight); } -CandidateSplitting::CandidateSplitting(int n_split, double Emin, double Emax, double n_bins, double minWeight, bool log) { - setNsplit(n_split); - setEnergyBins(Emin, Emax, n_bins, log); +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) { +CandidateSplitting::CandidateSplitting(double spectralIndex, double Emin, int nBins) { // to use with Diffusive Shock Acceleration - - if (SpectralIndex <= 0){ + 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)); + double dE = pow(1. / 2, 1. / (-spectralIndex + 1)); setEnergyBinsDSA(Emin, dE, nBins); - setMinimalWeight(1./pow(2, nBins)); - + setMinimalWeight(1. / pow(2, nBins)); } - void CandidateSplitting::process(Candidate *c) const { - double currE = c->current.getEnergy(); double prevE = c->previous.getEnergy(); @@ -46,7 +42,7 @@ void CandidateSplitting::process(Candidate *c) const { // minimal weight reached, no splitting return; } - if (currE < Ebins[0] || n_split == 0 ){ + if (currE < Ebins[0] || nSplit == 0 ){ // current energy is smaller than first bin -> no splitting // or, number of splits = 0 return; @@ -63,9 +59,9 @@ void CandidateSplitting::process(Candidate *c) const { for (size_t j = i; j < Ebins.size(); ++j ){ // adapted from Acceleration Module: - c->updateWeight(1. / n_split); // * 1/n_split + c->updateWeight(1. / nSplit); // * 1/n_split - for (int i = 1; i < n_split; i++) { + for (int i = 1; i < nSplit; i++) { ref_ptr new_candidate = c->clone(false); new_candidate->parent = c; @@ -76,79 +72,58 @@ void CandidateSplitting::process(Candidate *c) const { c->addSecondary(new_candidate); Candidate::setNextSerialNumber(snr + 1); } - - 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 n_bins, bool log) { +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)/n_bins; - - for (size_t i = 0; i < n_bins; ++i) { + double dE = (Emax-Emin)/nBins; + for (size_t i = 0; i < nBins; ++i) { if (log == true) { - Ebins.push_back(Emin * pow(Emax / Emin, i / (n_bins - 1.0))); + 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) { - - n_split = n; + nSplit = n; } void CandidateSplitting::setMinimalWeight(double w) { - minWeight = w; } int CandidateSplitting::getNsplit() const { - - return n_split; + return nSplit; } double CandidateSplitting::getMinimalWeight() const { - return minWeight; } - } // end namespace crpropa diff --git a/test/testCandidateSplitting.cpp b/test/testCandidateSplitting.cpp index ab5f24115..2949dec4b 100644 --- a/test/testCandidateSplitting.cpp +++ b/test/testCandidateSplitting.cpp @@ -10,46 +10,41 @@ namespace crpropa { TEST(testCandidateSplitting, SimpleTest) { - - - int n_split = 2; - int n_bins = 4; - double minWeight = pow(1./n_split, 2); + int nSplit = 2; + int nBins = 4; + double minWeight = pow(1. / nSplit, 2); double Emin = 1*GeV; double Emax = 10*GeV; - CandidateSplitting split_lin(n_split, Emin, Emax, n_bins, minWeight); - double dE = (Emax-Emin)/n_bins; + 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(), n_split); + EXPECT_EQ(split_lin.getNsplit(), nSplit); EXPECT_DOUBLE_EQ(split_lin.getMinimalWeight(), minWeight); - CandidateSplitting split_log(n_split, Emin, Emax, n_bins, minWeight, true); - double dE_log = pow(Emax / Emin, 1./(n_bins - 1.0)); + 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, n_bins); - double dE_dsa = pow(1./2, 1./(-spectralIndex + 1)); + 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()[n_bins-1], Emin * pow(dE_dsa, n_bins)); - - + EXPECT_DOUBLE_EQ(split_dsa.getEnergyBins()[nBins-1], Emin * pow(dE_dsa, nBins)); } TEST(testCandidateSplitting, CheckSplits) { - - int n_split = 2; - int n_bins = 3; + int nSplit = 2; + int nBins = 3; double Emin = 1*GeV; double Emax = 10*GeV; - double minWeight = pow(1./n_split, 4); + double minWeight = pow(1. / nSplit, 4); - CandidateSplitting splitting(n_split, Emin, Emax, n_bins, minWeight); + CandidateSplitting splitting(nSplit, Emin, Emax, nBins, minWeight); Candidate c(nucleusId(1,1),0.5*GeV); double weight = 1.0; @@ -58,13 +53,13 @@ TEST(testCandidateSplitting, CheckSplits) { c.current.setEnergy(2*GeV); splitting.process(&c); // 1. split - weight = weight/n_split; + weight = weight/nSplit; EXPECT_DOUBLE_EQ(c.getWeight(), weight); c.previous.setEnergy(2*GeV); c.current.setEnergy(6*GeV); splitting.process(&c); // 2. split - weight = weight/n_split; + weight = weight/nSplit; EXPECT_DOUBLE_EQ(c.getWeight(), weight); c.previous.setEnergy(6*GeV); @@ -75,7 +70,7 @@ TEST(testCandidateSplitting, CheckSplits) { c.current.setEnergy(6*GeV); splitting.process(&c); // 3. & 4. split, crosses two boundaries - weight = weight/n_split/n_split; + weight = weight/nSplit/nSplit; EXPECT_DOUBLE_EQ(c.getWeight(), weight); c.previous.setEnergy(6*GeV); @@ -83,8 +78,6 @@ TEST(testCandidateSplitting, CheckSplits) { splitting.process(&c); // no split, minimal weight reached EXPECT_DOUBLE_EQ(c.getWeight(), weight); c.previous.setEnergy(8*GeV); - } - } //namespace crpropa From 234c04ed0fa8b0b3481be835d55fb22dc251343e Mon Sep 17 00:00:00 2001 From: Sophie Aerdker Date: Fri, 6 Oct 2023 15:30:53 +0200 Subject: [PATCH 3/5] removed manual update of serial number and included serial number in splitting test --- src/module/CandidateSplitting.cpp | 6 +----- test/testCandidateSplitting.cpp | 33 +++++++++++++++++++------------ 2 files changed, 21 insertions(+), 18 deletions(-) diff --git a/src/module/CandidateSplitting.cpp b/src/module/CandidateSplitting.cpp index 8e89ed838..c60719555 100644 --- a/src/module/CandidateSplitting.cpp +++ b/src/module/CandidateSplitting.cpp @@ -64,13 +64,9 @@ void CandidateSplitting::process(Candidate *c) const { 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; - uint64_t snr = Candidate::getNextSerialNumber(); - new_candidate->setSerialNumber(snr); new_candidate->previous.setEnergy(currE); // so that new candidate is not split again in next step! - //InteractionTag is PRIM, physically no new particles are created - c->addSecondary(new_candidate); - Candidate::setNextSerialNumber(snr + 1); } if (j < Ebins.size()-1 && currE < Ebins[j+1]){ // candidate is in energy bin [j, j+1] -> no further splitting diff --git a/test/testCandidateSplitting.cpp b/test/testCandidateSplitting.cpp index 2949dec4b..8eca2cf1f 100644 --- a/test/testCandidateSplitting.cpp +++ b/test/testCandidateSplitting.cpp @@ -40,44 +40,51 @@ TEST(testCandidateSplitting, SimpleTest) { TEST(testCandidateSplitting, CheckSplits) { int nSplit = 2; int nBins = 3; - double Emin = 1*GeV; - double Emax = 10*GeV; + 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*GeV); + 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*GeV); + c.current.setEnergy(2); splitting.process(&c); // 1. split weight = weight/nSplit; EXPECT_DOUBLE_EQ(c.getWeight(), weight); - c.previous.setEnergy(2*GeV); + EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial + 1); + c.previous.setEnergy(2); - c.current.setEnergy(6*GeV); + c.current.setEnergy(6); splitting.process(&c); // 2. split weight = weight/nSplit; EXPECT_DOUBLE_EQ(c.getWeight(), weight); - c.previous.setEnergy(6*GeV); + EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial + 2); + c.previous.setEnergy(6); - c.current.setEnergy(0.5*GeV); + c.current.setEnergy(0.5); splitting.process(&c); // no split, cooling EXPECT_DOUBLE_EQ(c.getWeight(), weight); - c.previous.setEnergy(0.5*GeV); + EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial + 2); + c.previous.setEnergy(0.5); - c.current.setEnergy(6*GeV); + c.current.setEnergy(6); splitting.process(&c); // 3. & 4. split, crosses two boundaries weight = weight/nSplit/nSplit; EXPECT_DOUBLE_EQ(c.getWeight(), weight); - c.previous.setEnergy(6*GeV); + EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial + 4); + c.previous.setEnergy(6); - c.current.setEnergy(8*GeV); + c.current.setEnergy(8); splitting.process(&c); // no split, minimal weight reached EXPECT_DOUBLE_EQ(c.getWeight(), weight); - c.previous.setEnergy(8*GeV); + EXPECT_DOUBLE_EQ(c.getNextSerialNumber(), serial + 4); + c.previous.setEnergy(8); } } //namespace crpropa From ffdc93774a35b7f141a65dd0dff7d551928db87d Mon Sep 17 00:00:00 2001 From: Sophie Aerdker Date: Mon, 9 Oct 2023 09:38:44 +0200 Subject: [PATCH 4/5] negative power-law index, default arguments in constructor --- include/crpropa/module/CandidateSplitting.h | 11 +---------- src/module/CandidateSplitting.cpp | 13 +++---------- test/testCandidateSplitting.cpp | 18 +++++++++--------- 3 files changed, 13 insertions(+), 29 deletions(-) diff --git a/include/crpropa/module/CandidateSplitting.h b/include/crpropa/module/CandidateSplitting.h index 70ec97baf..bedf091db 100644 --- a/include/crpropa/module/CandidateSplitting.h +++ b/include/crpropa/module/CandidateSplitting.h @@ -31,15 +31,6 @@ class CandidateSplitting: public Module { CandidateSplitting(); - /** Constructor - @param nSplit Number of copies candidates are split - @param Emin Minimal energy for splitting - @param Emax Maximal energy for splitting - @param minWeight Mimimal Weight - @param nBins Number of energy bins - */ - CandidateSplitting(int nSplit, double Emin, double Emax, double nBins, double minWeight); - /** Constructor @param nSplit Number of copies candidates are split @param Emin Minimal energy for splitting @@ -48,7 +39,7 @@ class CandidateSplitting: public Module { @param minWeight Mimimal Weight @param log Energy bins in log */ - CandidateSplitting(int nSplit, double Emin, double Emax, double nBins, double minWeight, bool log); + CandidateSplitting(int nSplit, double Emin, double Emax, double nBins, double minWeight, bool log = false); /** Constructor @param spectralIndex Absolute value of expected spectral index determines splitting number diff --git a/src/module/CandidateSplitting.cpp b/src/module/CandidateSplitting.cpp index c60719555..6b735d567 100644 --- a/src/module/CandidateSplitting.cpp +++ b/src/module/CandidateSplitting.cpp @@ -8,13 +8,6 @@ CandidateSplitting::CandidateSplitting() { setMinimalWeight(1.); } -CandidateSplitting::CandidateSplitting(int nSplit, double Emin, double Emax, double nBins, double minWeight) { - // linear energy bins: - setNsplit(nSplit); - setEnergyBins(Emin, Emax, nBins, false); - setMinimalWeight(minWeight); -} - CandidateSplitting::CandidateSplitting(int nSplit, double Emin, double Emax, double nBins, double minWeight, bool log) { setNsplit(nSplit); setEnergyBins(Emin, Emax, nBins, log); @@ -23,13 +16,13 @@ CandidateSplitting::CandidateSplitting(int nSplit, double Emin, double Emax, do CandidateSplitting::CandidateSplitting(double spectralIndex, double Emin, int nBins) { // to use with Diffusive Shock Acceleration - if (spectralIndex <= 0){ + if (spectralIndex > 0){ throw std::runtime_error( - "CandidateSplitting: spectralIndex <= 0 !"); + "CandidateSplitting: spectralIndex > 0 !"); } setNsplit(2); // always split in 2, calculate bins in energy for given spectrum: - double dE = pow(1. / 2, 1. / (-spectralIndex + 1)); + double dE = pow(1. / 2, 1. / (spectralIndex + 1)); setEnergyBinsDSA(Emin, dE, nBins); setMinimalWeight(1. / pow(2, nBins)); } diff --git a/test/testCandidateSplitting.cpp b/test/testCandidateSplitting.cpp index 8eca2cf1f..dafdd23f6 100644 --- a/test/testCandidateSplitting.cpp +++ b/test/testCandidateSplitting.cpp @@ -13,13 +13,13 @@ TEST(testCandidateSplitting, SimpleTest) { int nSplit = 2; int nBins = 4; double minWeight = pow(1. / nSplit, 2); - double Emin = 1*GeV; - double Emax = 10*GeV; + double Emin = 1; // dimensionless for testing + double Emax = 10; CandidateSplitting split_lin(nSplit, Emin, Emax, nBins, minWeight); - double dE = (Emax-Emin)/nBins; + double dE = (Emax - Emin) / nBins; EXPECT_DOUBLE_EQ(split_lin.getEnergyBins()[0], Emin); - EXPECT_DOUBLE_EQ(split_lin.getEnergyBins()[1], Emin+dE); + EXPECT_DOUBLE_EQ(split_lin.getEnergyBins()[1], Emin + dE); EXPECT_EQ(split_lin.getNsplit(), nSplit); EXPECT_DOUBLE_EQ(split_lin.getMinimalWeight(), minWeight); @@ -27,13 +27,13 @@ TEST(testCandidateSplitting, SimpleTest) { 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); + EXPECT_DOUBLE_EQ(split_log.getEnergyBins()[1], Emin * dE_log); - double spectralIndex = 2.; + 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)); + 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)); } From 612a1b6233f463663b0c576f75637ff2bebc2817 Mon Sep 17 00:00:00 2001 From: Sophie Aerdker Date: Mon, 9 Oct 2023 09:50:26 +0200 Subject: [PATCH 5/5] improve documentation --- include/crpropa/module/CandidateSplitting.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/include/crpropa/module/CandidateSplitting.h b/include/crpropa/module/CandidateSplitting.h index bedf091db..514775624 100644 --- a/include/crpropa/module/CandidateSplitting.h +++ b/include/crpropa/module/CandidateSplitting.h @@ -14,6 +14,9 @@ namespace crpropa { +/** @addtogroup Acceleration +* @{ + */ /** @class CandidateSplitting @@ -42,7 +45,7 @@ class CandidateSplitting: public Module { CandidateSplitting(int nSplit, double Emin, double Emax, double nBins, double minWeight, bool log = false); /** Constructor - @param spectralIndex Absolute value of expected spectral index determines splitting number + @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 */