Merge pull request #315 from rafaelab/urb

add missing URB photon fields

CHANGELOG.md

@@ -8,6 +8,8 @@
 
 ### New features:
 
+* New model for the radio background added: Nitu et al. 2021 measurements.
+* New model for the radio background added: Fixsen et al. 2011 (ARCADE-2) measurements.
 * Weighted sampling thinning of electromagnetic processes  (EMPairProduction,
   EMInverseComptonScattering, EMDoublePairProduction, EMTripletPairProduction).
 * Planck JF12b variant of the JF12Field. See arXiv:1601.00546. Thanks to
@@ -21,7 +23,7 @@
 * New class-based interface for turbulent fields introduced
 * New turbulence modules implemented:
   - GridTurbulence (with the bendover scale) which should in general be used
-    instead of initTurbulence (before, it was implmeneted as
+    instead of initTurbulence (before, it was implemented as
     initTurbulenceWithBendover);
   - SimpleGridTurbulence which provides the exact field as initTurbulence;
   - HelicalGridTurbulence which provides the exact field as initHelicalTurbulence

CMakeLists.txt

@@ -318,10 +318,12 @@ endif()
 if(APPLE)
   add_definitions(-arch x86_64)
 endif(APPLE)
+
+# ----------------------------------------------------------------------------
 # Download data files (interaction data, masses, decay data ...)
 # ----------------------------------------------------------------------------
-OPTION(DOWNLOAD_DATA "Download CRProap Data files" ON)
-set(CRPROPA_DATAFILE_VER "2020-05-25")
+OPTION(DOWNLOAD_DATA "Download CRPropa data files" ON)
+set(CRPROPA_DATAFILE_VER "2020-11-28")
 if(DOWNLOAD_DATA)
   message("-- Downloading data file from crpropa.desy.de ~ 50 MB")
   file(DOWNLOAD

include/crpropa/PhotonBackground.h

@@ -186,6 +186,47 @@ class IRB_Stecker16_lower: public TabularPhotonField {
 	IRB_Stecker16_lower() : TabularPhotonField("IRB_Stecker16_lower", true) {}
 };
 
+/**
+ @class URB
+ @brief Extragalactic background light model from Protheroe & Biermann 1996
+
+ Source info:
+ DOI:10.1016/S0927-6505(96)00041-2
+ https://www.sciencedirect.com/science/article/abs/pii/S0927650596000412
+ */
+class URB_Protheroe96: public TabularPhotonField {
+public:
+	URB_Protheroe96() : TabularPhotonField("URB_Protheroe96", false) {}
+};
+
+/**
+ @class URB
+ @brief Extragalactic background light model based on ARCADE2 observations, by Fixsen et al.
+ Note that this model does not cover the same energy range as other URB models. Here, only ~10 MHz - 10 GHz is considered. 
+ Therefore, it only makes sense to use this model in very specific studies.
+  
+ Source info:
+ DOI:10.1088/0004-637X/734/1/5 
+ https://iopscience.iop.org/article/10.1088/0004-637X/734/1/5
+ */
+class URB_Fixsen11: public TabularPhotonField {
+public:
+	URB_Fixsen11() : TabularPhotonField("URB_Fixsen11", false) {}
+};
+
+/**
+ @class URB
+ @brief Extragalactic background light model by Nitu et al. 
+  
+ Source info:
+ DOI:10.1016/j.astropartphys.2020.102532
+ https://www.sciencedirect.com/science/article/pii/S0927650520301043?
+ */
+class URB_Nitu21: public TabularPhotonField {
+public:
+	URB_Nitu21() : TabularPhotonField("URB_Nitu21", false) {}
+};
+
 /**
  @class BlackbodyPhotonField
  @brief Photon field decorator for black body photon fields.

include/crpropa/module/EMDoublePairProduction.h

@@ -1,9 +1,11 @@
 #ifndef CRPROPA_EMDOUBLEPAIRPRODUCTION_H
 #define CRPROPA_EMDOUBLEPAIRPRODUCTION_H
 
+#include <fstream>
+#include <cmath>
+
 #include "crpropa/Module.h"
 #include "crpropa/PhotonBackground.h"
-#include <fstream>
 
 namespace crpropa {
 

include/crpropa/module/EMInverseComptonScattering.h

@@ -1,9 +1,11 @@
 #ifndef CRPROPA_EMINVERSECOMPTONSCATTERING_H
 #define CRPROPA_EMINVERSECOMPTONSCATTERING_H
 
+#include <fstream>
+#include <cmath>
+
 #include "crpropa/Module.h"
 #include "crpropa/PhotonBackground.h"
-#include <fstream>
 
 namespace crpropa {
 

include/crpropa/module/EMPairProduction.h

@@ -1,9 +1,12 @@
 #ifndef CRPROPA_EMPAIRPRODUCTION_H
 #define CRPROPA_EMPAIRPRODUCTION_H
 
+#include <fstream>
+#include <cmath>
+
 #include "crpropa/Module.h"
 #include "crpropa/PhotonBackground.h"
-#include <fstream>
+
 
 namespace crpropa {
 

include/crpropa/module/EMTripletPairProduction.h

@@ -1,9 +1,11 @@
 #ifndef CRPROPA_EMTRIPLETPAIRPRODUCTION_H
 #define CRPROPA_EMTRIPLETPAIRPRODUCTION_H
 
+#include <fstream>
+#include <cmath>
+
 #include "crpropa/Module.h"
 #include "crpropa/PhotonBackground.h"
-#include <fstream>
 
 namespace crpropa {
 /**

src/module/EMPairProduction.cpp

@@ -157,9 +157,9 @@ class PPSecondariesEnergyDistribution {
 			Random &random = Random::instance();
 			size_t j = random.randBin(s0) + 1;
 
-			double s_min = 4 * mec2 * mec2;
-			double beta = sqrt(1 - s_min / s);
-			double x0 = (1 - beta) / 2.;
+			double s_min = 4. * mec2 * mec2;
+			double beta = sqrtl(1. - s_min / s);
+			double x0 = (1. - beta) / 2.;
 			double dx = log((1 + beta) / (1 - beta)) / N;
 			double binWidth = x0 * (exp(j*dx) - exp((j-1)*dx));
 			if (random.rand() < 0.5)
@@ -199,6 +199,10 @@ void EMPairProduction::performInteraction(Candidate *candidate) const {
 	double Ep = E - Ee;
 	double f = Ep / E;
 
+	// for some backgrounds Ee=nan due to precision limitations.
+	if (not std::isfinite(Ee) || not std::isfinite(Ep))
+		return;
+
 	// sample random position along current step
 	Vector3d pos = random.randomInterpolatedPosition(candidate->previous.getPosition(), candidate->current.getPosition());
 	double w0 = candidate->getWeight();
@@ -230,16 +234,21 @@ void EMPairProduction::process(Candidate *candidate) const {
 	double rate = interpolate(E, tabEnergy, tabRate);
 	rate *= pow_integer<2>(1 + z) * photonField->getRedshiftScaling(z);
 
-	Random &random = Random::instance();
-	double randDistance = -log(random.rand()) / rate;
+	// run this loop at least once to limit the step size 
 	double step = candidate->getCurrentStep();
-	if (step < randDistance) {
-		candidate->limitNextStep(limit / rate);
-		return;
-	} else { // perform interaction and stop tracking particle 
-		performInteraction(candidate);
-		return;
-	}
+	Random &random = Random::instance();
+	do {
+		double randDistance = -log(random.rand()) / rate;
+		// check for interaction; if it doesn't ocurr, limit next step
+		if (step < randDistance) { 
+			candidate->limitNextStep(limit / rate);
+		} else {
+			performInteraction(candidate);
+			return;
+		}
+		step -= randDistance; 
+	} while (step > 0.);
+
 }
 
 } // namespace crpropa