From 3095d66766cc653359dde66bfc4fbd1a96d74b34 Mon Sep 17 00:00:00 2001 From: Nicholas Kamp Date: Wed, 6 Nov 2024 18:39:31 -0500 Subject: [PATCH] fix log integration --- .../primary/energy_direction/Tabulated2DFluxDistribution.cxx | 5 +++-- projects/utilities/public/SIREN/utilities/Integration.h | 2 +- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/projects/distributions/private/primary/energy_direction/Tabulated2DFluxDistribution.cxx b/projects/distributions/private/primary/energy_direction/Tabulated2DFluxDistribution.cxx index 5acd448c..06081150 100644 --- a/projects/distributions/private/primary/energy_direction/Tabulated2DFluxDistribution.cxx +++ b/projects/distributions/private/primary/energy_direction/Tabulated2DFluxDistribution.cxx @@ -32,7 +32,8 @@ Tabulated2DFluxDistribution::Tabulated2DFluxDistribution() {} void Tabulated2DFluxDistribution::ComputeIntegral() { std::function integrand = [&] (double x, double y) -> double { - return x*unnormed_pdf(pow(10,x),y); + //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); }; integral = siren::utilities::simpsonIntegrate2D(integrand, log10(energyMin), log10(energyMax), zenithMin, zenithMax); } @@ -92,7 +93,7 @@ void Tabulated2DFluxDistribution::LoadFluxTable(std::vector & energies, table_data.y = zeniths; table_data.f = flux; energy_nodes = energies; - zenith_nodes = energies; + zenith_nodes = zeniths; // If no physical are manually set, use first/last entry of table if(not energy_bounds_set) { diff --git a/projects/utilities/public/SIREN/utilities/Integration.h b/projects/utilities/public/SIREN/utilities/Integration.h index df3b7204..e87b6325 100644 --- a/projects/utilities/public/SIREN/utilities/Integration.h +++ b/projects/utilities/public/SIREN/utilities/Integration.h @@ -205,11 +205,11 @@ double simpsonIntegrate2D(const FuncType& func, double a1, double b1, double a2, if (j==0 || j==N1i) c1 = 1; else if (j%2==0) c1 = 2; else if (j%2==1) c1 = 4; + double x = a1 + j*h1i; for(unsigned int k = 0; k < N2i+1; k++) { if (k==0 || k==N2i) c2 = 1; else if (k%2==0) c2 = 2; else if (k%2==1) c2 = 4; - double x = a1 + j*h1i; double y = a2 + k*h2i; estimate += c1*c2*func(x, y); }