Candidate Splitting module

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)

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"

include/crpropa/module/CandidateSplitting.h

@@ -0,0 +1,80 @@
+#ifndef CRPROPA_CANDIDATEPLITTING_H
+#define CRPROPA_CANDIDATEPLITTING_H
+
+#include <string>
+#include <iostream>
+#include <cmath>
+#include <cstdlib>
+#include <sstream>
+
+#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<double> 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<double>& getEnergyBins() const;
+
+};
+/** @}*/
+
+} // namespace crpropa
+#endif // CRPROPA_PARTICLESPLITTING_H

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<int>;
 %include "crpropa/module/Tools.h"

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<Candidate> 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<double>& 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
+

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 <stdexcept>
+#include <cmath>
+
+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