fix issues with 2D tabular integral calculation, sampling errors with…

… very short decay lengths, DarkNews errors, and kinematic erros in Dipole portal HNL DIS

projects/detector/private/DetectorModel.cxx

@@ -44,6 +44,24 @@ namespace {
 void string_to_lower(std::string & data) {
     std::transform(data.begin(), data.end(), data.begin(), [](unsigned char c){ return std::tolower(c); });
 }
+double estimateDotError(const std::array<double, 3>& p0, const std::array<double, 3>& p1, const std::array<double, 3>& int_dir) {
+    constexpr double epsilon = std::numeric_limits<double>::epsilon();
+
+    double dot_product = 0.0;
+    double error_sum = 0.0;
+
+    for (int i = 0; i < 3; ++i) {
+        double diff = p1[i] - p0[i];
+        double prod = int_dir[i] * diff;
+        dot_product += prod;
+        error_sum += std::fabs(int_dir[i] * diff);
+    }
+
+    // Estimated error of the expression
+    double error_estimate = epsilon * error_sum;
+
+    return error_estimate;
+}
 
 template <class InIt>
 typename std::iterator_traits<InIt>::value_type accumulate(InIt begin, InIt end) {
@@ -696,6 +714,11 @@ double DetectorModel::GetInteractionDensity(Geometry::IntersectionList const & i
     } else {
         direction.normalize();
     }
+    // If we have only decays, avoid the sector loop
+    // total_decay_length is in m
+    if(targets.empty()) {
+        return 1.0 / total_decay_length;
+    }
     double dot = direction * intersections.direction;
     assert(std::abs(1.0 - std::abs(dot)) < 1e-6);
     double offset = (intersections.position - p0) * direction;
@@ -706,11 +729,6 @@ double DetectorModel::GetInteractionDensity(Geometry::IntersectionList const & i
         dot = 1;
     }
 
-    // If we have only decays, avoid the sector loop
-    // total_decay_length is in m
-    if(targets.empty()) {
-        return 1.0 / total_decay_length;
-    }
 
     double interaction_density = std::numeric_limits<double>::quiet_NaN();
 
@@ -991,11 +1009,22 @@ double DetectorModel::GetInteractionDepthInCGS(Geometry::IntersectionList const
     }
     Vector3D direction = p1 - p0;
     double distance = direction.magnitude();
+    // this dot error turned out to be much smaller than 1e-6 for event that failed the assertion below
+    // double dot_error = estimateDotError(std::array<double, 3>(p0.get()),
+    //                                     std::array<double, 3>(p1.get()),
+    //                                     std::array<double, 3>(intersections.direction));
+
+    // If we have only decays, avoid the sector loop
+    if(targets.empty()) {
+      return distance/total_decay_length; // m / m --> dimensionless
+    }
     if(distance == 0.0) {
         return 0.0;
     }
     direction.normalize();
 
+    // TODO: a better numerical precision check when the traversed distance is very small
+    // this functionally only happens for decays right now, so we just check for decays at the top
     double dot = intersections.direction * direction;
     assert(std::abs(1.0 - std::abs(dot)) < 1e-6);
     double offset = (intersections.position - p0) * direction;
@@ -1006,11 +1035,6 @@ double DetectorModel::GetInteractionDepthInCGS(Geometry::IntersectionList const
         dot = 1;
     }
 
-    // If we have only decays, avoid the sector loop
-    if(targets.empty()) {
-      return distance/total_decay_length; // m / m --> dimensionless
-    }
-
     std::vector<double> interaction_depths(targets.size(), 0.0);
 
     std::function<bool(std::vector<Geometry::Intersection>::const_iterator, std::vector<Geometry::Intersection>::const_iterator, double)> callback =
@@ -1339,6 +1363,11 @@ double DetectorModel::DistanceForInteractionDepthFromPoint(Geometry::Intersectio
         interaction_depth *= -1;
         direction = -direction;
     }
+    // If we have only decays, avoid the sector loop
+    // decay length in m
+    if(targets.empty()) {
+      return interaction_depth * total_decay_length;
+    }
 
     double dot = intersections.direction * direction;
     assert(std::abs(1.0 - std::abs(dot)) < 1e-6);
@@ -1350,11 +1379,6 @@ double DetectorModel::DistanceForInteractionDepthFromPoint(Geometry::Intersectio
         dot = 1;
     }
 
-    // If we have only decays, avoid the sector loop
-    // decay length in m
-    if(targets.empty()) {
-      return interaction_depth * total_decay_length;
-    }
 
     // Recast decay length to cm for density integral
     double total_decay_length_cm = total_decay_length / siren::utilities::Constants::cm;

projects/distributions/private/primary/energy_direction/Tabulated2DFluxDistribution.cxx

@@ -31,11 +31,22 @@ namespace {
 Tabulated2DFluxDistribution::Tabulated2DFluxDistribution() {}
 
 void Tabulated2DFluxDistribution::ComputeIntegral() {
+    // Check if the table is log in x (energy). If so, compute the integral in log space
+    // assuing we are never log in y (zenith)
+    double eMin = energyMin;
+    double eMax = energyMax;
     std::function<double(double, double)> integrand = [&] (double x, double y) -> double {
-        //std::cout << "x " << x << " y " << y << " z " << pow(10,x)*unnormed_pdf(pow(10,x),y) << std::endl;
-        return pow(10,x)*unnormed_pdf(pow(10,x),y)/log(10);
+       return unnormed_pdf(x,y);
     };
-    integral = siren::utilities::simpsonIntegrate2D(integrand, log10(energyMin), log10(energyMax), zenithMin, zenithMax);
+    if (fluxTable.IsLogX()) {
+        eMin = log10(energyMin);
+        eMax = log10(energyMax);
+        integrand =  [&] (double x, double y) -> double {
+            return log(10)*pow(10,x)*unnormed_pdf(pow(10,x),y);
+        };
+    }
+
+    integral = siren::utilities::simpsonIntegrate2D(integrand, eMin, eMax, zenithMin, zenithMax);
 }
 
 void Tabulated2DFluxDistribution::LoadFluxTable() {

projects/interactions/private/DISFromSpline.cxx

@@ -538,9 +538,7 @@ void DISFromSpline::SampleFinalState(dataclasses::CrossSectionDistributionRecord
     double Q2 = 2 * E1_lab * E2_lab * pow(10.0, kin_vars[1] + kin_vars[2]);
     double p1x_lab = std::sqrt(p1_lab.px() * p1_lab.px() + p1_lab.py() * p1_lab.py() + p1_lab.pz() * p1_lab.pz());
     double pqx_lab = (m1*m1 + m3*m3 + 2 * p1x_lab * p1x_lab + Q2 + 2 * E1_lab * E1_lab * (final_y - 1)) / (2.0 * p1x_lab);
-    //double pqx_lab_prime = (Q2 + m3*m3 + 2*E1_lab*E1_lab*final_y)/(2*E1_lab);
     double momq_lab = std::sqrt(m1*m1 + p1x_lab*p1x_lab + Q2 + E1_lab * E1_lab * (final_y * final_y - 1));
-    //double momq_lab_prime = std::sqrt(Q2 + pow(E1_lab*final_y,2));
     // rounding error check
     double pqy_lab;
     if (pqx_lab>momq_lab){
@@ -549,7 +547,6 @@ void DISFromSpline::SampleFinalState(dataclasses::CrossSectionDistributionRecord
     }
     else pqy_lab = std::sqrt(momq_lab*momq_lab - pqx_lab *pqx_lab);
     double Eq_lab = E1_lab * final_y;
-    //double Eq_lab_prime = std::sqrt(momq_lab*momq_lab - Q2);
 
     geom3::UnitVector3 x_dir = geom3::UnitVector3::xAxis();
     geom3::Vector3 p1_mom = p1_lab.momentum();

projects/interactions/private/HNLDISFromSpline.cxx

@@ -346,7 +346,7 @@ double HNLDISFromSpline::DifferentialCrossSection(siren::dataclasses::Particle::
 
 double HNLDISFromSpline::InteractionThreshold(dataclasses::InteractionRecord const & interaction) const {
     // Using center of mass frame
-    // See
+    // require E_cm > m_HNL
     double M_iso = siren::utilities::Constants::isoscalarMass;
     return (std::pow(hnl_mass_ + M_iso,2) - M_iso*M_iso)/(2*M_iso);
 }

projects/interactions/private/HNLDipoleDISFromSpline.cxx

@@ -28,6 +28,7 @@ namespace interactions {
 namespace {
 ///Check whether a given point in phase space is physically realizable.
 ///Based on equation 5 of https://arxiv.org/pdf/1707.08573.pdf
+///Also eq 6-8 of https://arxiv.org/pdf/hep-ph/0208187
 ///\param x Bjorken x of the interaction
 ///\param y Bjorken y of the interaction
 ///\param E Incoming neutrino in energy in the lab frame ($E_\nu$)
@@ -38,7 +39,18 @@ bool kinematicallyAllowed(double x, double y, double E, double M, double m) {
     double W2 = M*M + Q2/x * (1-x);
     double Er = E*y;
     double term = m*m - W2 - 2*x*E*M - x*x*M*M + 2*Er*(x*M + E);
-    return Er*Er - W2 - term*term/(4*E*E) > 0; // equation 5
+    if (Er*Er - W2 - term*term/(4*E*E) <= 0) return false; // equation 5
+    double x_low = m*m / (2*M*(E-m));
+    if (x < x_low || x > 1) return false; // equation 6 of 0208187
+    double a_num = 1 - m*m*(1/(2*M*E*x) + 1/(2*E*E));
+    double a_denom = 2*(1 + M*x/(2*E));
+    double b_num = std::sqrt(pow(1 - m*m/(2*M*E*x),2) - m*m/(E*E));
+    double b_denom = 2*(1 + M*x/(2*E));
+    double a = a_num/a_denom;
+    double b = b_num/b_denom;
+    if (y < a-b || y > a+b) return false; // equation 8 of 0208187
+    return true;
+
 }
 }
 
@@ -326,8 +338,10 @@ double HNLDipoleDISFromSpline::DifferentialCrossSection(siren::dataclasses::Part
 }
 
 double HNLDipoleDISFromSpline::InteractionThreshold(dataclasses::InteractionRecord const & interaction) const {
-    // Consider implementing thershold at some point
-    return 0;
+    // Using center of mass frame
+    // require E_cm > m_HNL
+    double M_iso = siren::utilities::Constants::isoscalarMass;
+    return (std::pow(hnl_mass_ + M_iso,2) - M_iso*M_iso)/(2*M_iso);
 }
 
 void HNLDipoleDISFromSpline::SampleFinalState(dataclasses::CrossSectionDistributionRecord& interaction, std::shared_ptr<siren::utilities::SIREN_random> random) const {
@@ -485,10 +499,13 @@ void HNLDipoleDISFromSpline::SampleFinalState(dataclasses::CrossSectionDistribut
     double Q2 = 2 * E1_lab * E2_lab * pow(10.0, kin_vars[1] + kin_vars[2]);
     double p1x_lab = std::sqrt(p1_lab.px() * p1_lab.px() + p1_lab.py() * p1_lab.py() + p1_lab.pz() * p1_lab.pz());
     double pqx_lab = (m1*m1 + m3*m3 + 2 * p1x_lab * p1x_lab + Q2 + 2 * E1_lab * E1_lab * (final_y - 1)) / (2.0 * p1x_lab);
-    double momq_lab = std::sqrt(m1*m1 + p1x_lab*p1x_lab + Q2 + E1_lab * E1_lab * (final_y * final_y - 1));
+    //double momq_lab = std::sqrt(m1*m1 + p1x_lab*p1x_lab + Q2 + E1_lab * E1_lab * (final_y * final_y - 1));
+    double momq_lab = std::sqrt(Q2 + E1_lab*E1_lab*final_y*final_y);
     double pqy_lab;
     if (pqx_lab>momq_lab){
-        assert(((pqx_lab-momq_lab)/momq_lab)<1e-3);
+        std::cout << "WARNING: DIS sampling resulted in an x component of the momentum transfer larger than the total momentum transfer by ";
+        std::cout << ((pqx_lab-momq_lab)/momq_lab)*100 << "%" << std::endl;
+        std::cout << "Setting y component to zero" << std::endl;
         pqy_lab = 0;
     }
     else pqy_lab = std::sqrt(momq_lab*momq_lab - pqx_lab *pqx_lab);

projects/utilities/public/SIREN/utilities/Integration.h

@@ -150,23 +150,6 @@ double rombergIntegrate(const FuncType& func, double a, double b, double tol=1e-
   throw(std::runtime_error("Integral failed to converge"));
 }
 
-/**
-* @brief Performs a 2D trapezoidal integration of a function
-*
-*
-* @param func the function to integrate
-* @param a the lower bound of integration for the first dimension
-* @param b the upper bound of integration for the first dimension
-* @param c the lower bound of integration for the second dimension
-* @param d the upper bound of integration for the second dimension
-* @param tol the (absolute) tolerance on the error of the integral per dimension
-*/
-template<typename FuncType>
-double trapezoidIntegrate2D(const FuncType& func, double a1, double b1, double a2, double b2, unsigned int order=5){
-
-  return 0;
-}
-
 /**
 * @brief Performs a 2D simpson integration of a function
 *
@@ -179,11 +162,9 @@ double trapezoidIntegrate2D(const FuncType& func, double a1, double b1, double a
 * @param tol the (absolute) tolerance on the error of the integral per dimension
 */
 template<typename FuncType>
-double simpsonIntegrate2D(const FuncType& func, double a1, double b1, double a2, double b2, double tol=1e-6){
+double simpsonIntegrate2D(const FuncType& func, double a1, double b1, double a2, double b2,
+                          double tol=1e-3, const unsigned int N1 = 10, const unsigned int N2=10, const unsigned int maxIter=15){
 
-  const unsigned int N1=10;
-  const unsigned int N2=10;
-  const unsigned int maxIter=10;
   if(tol<0)
      throw(std::runtime_error("Integration tolerance must be positive"));
 
@@ -215,7 +196,6 @@ double simpsonIntegrate2D(const FuncType& func, double a1, double b1, double a2,
        }
      }
      estimate *= h1i*h2i/9;
-     std::cout << estimate << " " << std::abs((estimate-prev_estimate)/estimate) << std::endl;
      if(std::abs((estimate-prev_estimate)/estimate) < tol) {
       return estimate;
      }

python/SIREN_Controller.py

@@ -497,14 +497,14 @@ def Initialize(self, injection_filenames=None, weighter_filename=None):
         self.InitializeWeighter(filename=weighter_filename)
 
     # Generate events using the self.injector object
-    def GenerateEvents(self, N=None, fill_tables_at_exit=True):
+    def GenerateEvents(self, N=None, fill_tables_at_exit=True, verbose=True):
         if N is None:
             N = self.events_to_inject
         count = 0
         self.gen_times,self.global_times = [],[]
         prev_time = time.time()
         while (self.injector.InjectedEvents() < self.events_to_inject) and (count < N):
-            print("Injecting Event %d/%d  " % (count, N), end="\r")
+            if verbose: print("Injecting Event %d/%d  " % (count, N), end="\r")
             event = self.injector.GenerateEvent()
             self.events.append(event)
             t = time.time()

python/SIREN_DarkNews.py

@@ -669,6 +669,7 @@ def __init__(self, dec_case, table_dir=None):
 
         # Some variables for storing the decay phase space integrator
         self.decay_integrator = None
+        self.decay_dict = None
         self.decay_norm = None
         self.PS_samples = None
         self.PS_weights = None
@@ -702,6 +703,7 @@ def __init__(self, dec_case, table_dir=None):
     # serialization method
     def get_representation(self):
         return {"decay_integrator":self.decay_integrator,
+                "decay_dict":self.decay_dict,
                 "decay_norm":self.decay_norm,
                 "dec_case":self.dec_case,
                 "PS_samples":self.PS_samples,
@@ -716,7 +718,7 @@ def SetIntegratorAndNorm(self):
         int_file = os.path.join(self.table_dir, "decay_integrator.pkl")
         if os.path.isfile(int_file):
             with open(int_file, "rb") as ifile:
-                _, self.decay_integrator = pickle.load(ifile)
+                self.decay_dict, self.decay_integrator = pickle.load(ifile)
         # Try to find the normalization information
         norm_file = os.path.join(self.table_dir, "decay_norm.json")
         if os.path.isfile(norm_file):
@@ -832,7 +834,7 @@ def TotalDecayWidth(self, arg1):
                     self.SetIntegratorAndNorm()
                 else:
                     self.total_width = (
-                        self.decay_integrator["diff_decay_rate_0"].mean
+                        self.decay_dict["diff_decay_rate_0"].mean
                         * self.decay_norm["diff_decay_rate_0"]
                     )
             else: