Merge pull request #491 from LeanderSchlegel/empair

Small fix for EMPairProduction when haveElectrons=False and small Optimizations
CRPropa · Jun 4, 2024 · dc3299b · dc3299b
2 parents a6dbe04 + 14c7835
commit dc3299b
Show file tree

Hide file tree

Showing 3 changed files with 9 additions and 10 deletions.
diff --git a/src/module/EMDoublePairProduction.cpp b/src/module/EMDoublePairProduction.cpp
@@ -78,7 +78,6 @@ void EMDoublePairProduction::performInteraction(Candidate *candidate) const {
 
 	double f = Ee / E;
 
-	if (haveElectrons) {
 		if (random.rand() < pow(1 - f, thinning)) {
 			double w = 1. / pow(1 - f, thinning);
 			candidate->addSecondary( 11, Ee / (1 + z), pos, w, interactionTag);
@@ -87,7 +86,6 @@ void EMDoublePairProduction::performInteraction(Candidate *candidate) const {
 			double w = 1. / pow(f, thinning);
 			candidate->addSecondary(-11, Ee / (1 + z), pos, w, interactionTag);
 		}
-	}
 }
 
 void EMDoublePairProduction::process(Candidate *candidate) const {

diff --git a/src/module/EMInverseComptonScattering.cpp b/src/module/EMInverseComptonScattering.cpp
@@ -189,9 +189,9 @@ void EMInverseComptonScattering::performInteraction(Candidate *candidate) const
 	double Enew = distribution.sample(E, s);
 
 	// add up-scattered photon
-	double Esecondary = E - Enew;
-	double f = Enew / E;
 	if (havePhotons) {
+		double Esecondary = E - Enew;
+		double f = Enew / E;
 		if (random.rand() < pow(1 - f, thinning)) {
 			double w = 1. / pow(1 - f, thinning);
 			Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());

diff --git a/src/module/EMPairProduction.cpp b/src/module/EMPairProduction.cpp
@@ -173,13 +173,17 @@ class PPSecondariesEnergyDistribution {
 };
 
 void EMPairProduction::performInteraction(Candidate *candidate) const {
-	// scale particle energy instead of background photon energy
-	double z = candidate->getRedshift();
-	double E = candidate->current.getEnergy() * (1 + z);
+
+	// photon is lost after interacting
+	candidate->setActive(false);
 
 	// check if secondary electron pair needs to be produced
 	if (not haveElectrons)
 		return;
+
+	// scale particle energy instead of background photon energy
+	double z = candidate->getRedshift();
+	double E = candidate->current.getEnergy() * (1 + z);
 
 	// check if in tabulated energy range
 	if (E < tabE.front() or (E > tabE.back()))
@@ -203,9 +207,6 @@ void EMPairProduction::performInteraction(Candidate *candidate) const {
 	if (not std::isfinite(Ee) || not std::isfinite(Ep))
 		return;
 
-	// photon is lost after interacting
-	candidate->setActive(false);
-
 	// sample random position along current step
 	Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
 	// apply sampling