diff --git a/projects/detector/private/DetectorModel.cxx b/projects/detector/private/DetectorModel.cxx index 359ff59e..55cf0330 100644 --- a/projects/detector/private/DetectorModel.cxx +++ b/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& p0, const std::array& p1, const std::array& int_dir) { + constexpr double epsilon = std::numeric_limits::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 typename std::iterator_traits::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::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(p0.get()), + // std::array(p1.get()), + // std::array(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 interaction_depths(targets.size(), 0.0); std::function::const_iterator, std::vector::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; diff --git a/projects/distributions/private/primary/energy_direction/Tabulated2DFluxDistribution.cxx b/projects/distributions/private/primary/energy_direction/Tabulated2DFluxDistribution.cxx index 06081150..62abe5a8 100644 --- a/projects/distributions/private/primary/energy_direction/Tabulated2DFluxDistribution.cxx +++ b/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 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() { diff --git a/projects/interactions/private/DISFromSpline.cxx b/projects/interactions/private/DISFromSpline.cxx index fa1caf03..ecd74e6b 100644 --- a/projects/interactions/private/DISFromSpline.cxx +++ b/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(); diff --git a/projects/interactions/private/HNLDISFromSpline.cxx b/projects/interactions/private/HNLDISFromSpline.cxx index 10213e70..97267ccd 100644 --- a/projects/interactions/private/HNLDISFromSpline.cxx +++ b/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); } diff --git a/projects/interactions/private/HNLDipoleDISFromSpline.cxx b/projects/interactions/private/HNLDipoleDISFromSpline.cxx index ed6447be..d31885ea 100644 --- a/projects/interactions/private/HNLDipoleDISFromSpline.cxx +++ b/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 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); diff --git a/projects/utilities/public/SIREN/utilities/Integration.h b/projects/utilities/public/SIREN/utilities/Integration.h index e87b6325..9662eb9d 100644 --- a/projects/utilities/public/SIREN/utilities/Integration.h +++ b/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 -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 -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; } diff --git a/python/SIREN_Controller.py b/python/SIREN_Controller.py index dd7dc35c..f466d89d 100644 --- a/python/SIREN_Controller.py +++ b/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() diff --git a/python/SIREN_DarkNews.py b/python/SIREN_DarkNews.py index d63cb0ad..695bbaad 100644 --- a/python/SIREN_DarkNews.py +++ b/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: