fix log integration

Harvard-Neutrino · Nov 6, 2024 · 3095d66 · 3095d66
1 parent 80a2e0b
commit 3095d66
Show file tree

Hide file tree

Showing 2 changed files with 4 additions and 3 deletions.
diff --git 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<double(double, double)> 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<double> & 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
@@ -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);
        }