Compute Q2 more robustly before checking bounds

Harvard-Neutrino · Jan 19, 2024 · f5a6020 · f5a6020
1 parent e43bc10
commit f5a6020
Show file tree

Hide file tree

Showing 3 changed files with 8 additions and 5 deletions.
diff --git a/projects/interactions/private/DISFromSpline.cxx b/projects/interactions/private/DISFromSpline.cxx
@@ -322,10 +322,10 @@ double DISFromSpline::DifferentialCrossSection(dataclasses::InteractionRecord co
     double x = Q2 / (2.0 * p2.dot(q));
     double lepton_mass = particleMass(interaction.signature.secondary_types[lepton_index]);
 
-    return DifferentialCrossSection(primary_energy, x, y, lepton_mass);
+    return DifferentialCrossSection(primary_energy, x, y, lepton_mass, Q2);
 }
 
-double DISFromSpline::DifferentialCrossSection(double energy, double x, double y, double secondary_lepton_mass) const {
+double DISFromSpline::DifferentialCrossSection(double energy, double x, double y, double secondary_lepton_mass, double Q2) const {
     double log_energy = log10(energy);
     // check preconditions
     if(log_energy < differential_cross_section_.lower_extent(0)
@@ -339,7 +339,9 @@ double DISFromSpline::DifferentialCrossSection(double energy, double x, double y
     // we assume that:
     // the target is stationary so its energy is just its mass
     // the incoming neutrino is massless, so its kinetic energy is its total energy
-    double Q2 = 2.0 * energy * target_mass_ * x * y;
+    if(std::isnan(Q2)) {
+        Q2 = 2.0 * energy * target_mass_ * x * y;
+    }
     if(Q2 < minimum_Q2_) // cross section not calculated, assumed to be zero
         return 0;
 
@@ -402,6 +404,7 @@ void DISFromSpline::SampleFinalState(dataclasses::InteractionRecord& interaction
     double E1_lab = p1_lab.e();
     double E2_lab = p2_lab.e();
 
+
     // The out-going particle always gets at least enough energy for its rest mass
     double yMax = 1 - m / primary_energy;
     double logYMax = log10(yMax);

diff --git a/projects/interactions/private/pybindings/DISFromSpline.h b/projects/interactions/private/pybindings/DISFromSpline.h
@@ -33,7 +33,7 @@ void register_DISFromSpline(pybind11::module_ & m) {
         .def("TotalCrossSection",overload_cast<LI::dataclasses::Particle::ParticleType, double>(&DISFromSpline::TotalCrossSection, const_))
         .def("TotalCrossSection",overload_cast<LI::dataclasses::Particle::ParticleType, double, LI::dataclasses::Particle::ParticleType>(&DISFromSpline::TotalCrossSection, const_))
         .def("DifferentialCrossSection",overload_cast<LI::dataclasses::InteractionRecord const &>(&DISFromSpline::DifferentialCrossSection, const_))
-        .def("DifferentialCrossSection",overload_cast<double, double, double, double>(&DISFromSpline::DifferentialCrossSection, const_))
+        .def("DifferentialCrossSection",overload_cast<double, double, double, double, double>(&DISFromSpline::DifferentialCrossSection, const_))
         .def("InteractionThreshold",&DISFromSpline::InteractionThreshold)
         .def("GetPossibleTargets",&DISFromSpline::GetPossibleTargets)
         .def("GetPossibleTargetsFromPrimary",&DISFromSpline::GetPossibleTargetsFromPrimary)

diff --git a/projects/interactions/public/LeptonInjector/interactions/DISFromSpline.h b/projects/interactions/public/LeptonInjector/interactions/DISFromSpline.h
@@ -69,7 +69,7 @@ friend cereal::access;
     double TotalCrossSection(LI::dataclasses::Particle::ParticleType primary, double energy) const;
     double TotalCrossSection(LI::dataclasses::Particle::ParticleType primary, double energy, LI::dataclasses::Particle::ParticleType target) const override;
     double DifferentialCrossSection(dataclasses::InteractionRecord const &) const override;
-    double DifferentialCrossSection(double energy, double x, double y, double secondary_lepton_mass) const;
+    double DifferentialCrossSection(double energy, double x, double y, double secondary_lepton_mass, double Q2=std::numeric_limits<double>::quiet_NaN()) const;
     double InteractionThreshold(dataclasses::InteractionRecord const &) const override;
     void SampleFinalState(dataclasses::InteractionRecord &, std::shared_ptr<LI::utilities::LI_random> random) const override;