Merge pull request #442 from sophieaerdker/CandSplitting

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,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)

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,74 @@
+#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 {
+/** @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<double> 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<double>& getEnergyBins() const;
+
+};
+/** @}*/
+
+} // namespace crpropa
+#endif // CRPROPA_CANDIDATESPLITTING_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,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<Candidate> 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<double>& 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
+

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