make Histogram Utils and move non struct stuff from Struct headers

manager/CMakeLists.txt

@@ -26,7 +26,7 @@ if(NOT CPU_ONLY)
 endif()
 
 #If compiling with GPU is not enabled MaCh3GPUCompilerOptions will be empty
-target_link_libraries(Manager PUBLIC MaCh3CompilerOptions MaCh3GPUCompilerOptions yaml-cpp spdlog ROOT::Tree ROOT::Hist)
+target_link_libraries(Manager PUBLIC MaCh3CompilerOptions MaCh3GPUCompilerOptions yaml-cpp spdlog NuOscillator ROOT::Tree ROOT::Hist)
 target_link_libraries(Manager PRIVATE MaCh3Warnings)
 
 target_include_directories(Manager PUBLIC

mcmc/MCMCProcessor.h

@@ -38,6 +38,7 @@
 
 // MaCh3 includes
 #include "mcmc/StatisticalUtils.h"
+#include "samplePDF/HistogramUtils.h"
 
 //KS: Joy of forward declaration https://gieseanw.wordpress.com/2018/02/25/the-joys-of-forward-declarations-results-from-the-real-world/
 class TChain;

mcmc/StatisticalUtils.h

@@ -11,6 +11,7 @@
 //MaCh3 includes
 #include "samplePDF/Structs.h"
 #include "manager/manager.h"
+#include "samplePDF/HistogramUtils.h"
 
 /// @file StatisticalUtils.h
 /// @brief Utility functions for statistical interpretations in MaCh3

samplePDF/CMakeLists.txt

@@ -2,13 +2,14 @@ set(HEADERS
     samplePDFBase.h
     samplePDFFDBase.h
     Structs.h
-	FarDetectorCoreInfoStruct.h
+    HistogramUtils.h
+    FarDetectorCoreInfoStruct.h
 )
 
 add_library(SamplePDF SHARED
     samplePDFBase.cpp
     samplePDFFDBase.cpp
-    Structs.cpp
+    HistogramUtils.cpp
 )
 
 target_link_libraries(SamplePDF PUBLIC Splines NuOscillator Covariance)

samplePDF/Structs.cpp → samplePDF/HistogramUtils.cpp

@@ -1,124 +1,7 @@
-#include "samplePDF/Structs.h"
-
-#pragma GCC diagnostic push
-#pragma GCC diagnostic ignored "-Wfloat-conversion"
-#include "Constants/OscillatorConstants.h"
-#pragma GCC diagnostic pop
-
 #include "TList.h"
 #include "TObjArray.h"
 
-namespace MaCh3Utils {
-
-  // *****************************
-  // Get the mass of a particle from the PDG
-  // In GeV, not MeV!
-  double GetMassFromPDG(int PDG) {
-  // *****************************
-
-    switch (abs(PDG)) {  
-    case 11:
-      return 0.511E-3;
-      break;
-    case 13:
-      return 105.658E-3;
-      break;
-    case 15:
-      return 1.77682;
-      break;
-    case 22:
-      return 0.;
-      break;
-    case 211:
-      return 139.57E-3;
-      break;
-    case 111:
-      return 134.98E-3;
-      break;
-    case 2112:
-      return 939.565E-3;
-      break;
-    case 2212:
-      return 938.27E-3;
-      break;
-    //Oxygen nucleus
-    case 1000080160:
-      return 14.89926;
-      break;
-	//eta
-	case 221:
-	  return 547.862E-3;
-	  break;
-	  //K^0 (s or l)
-	case 311:
-	case 130:
-	case 310:
-	  return 497.611E-3;
-	  break;
-	case 321:
-	  return 493.677E-3;
-	  break;
-	// Lamda baryon
-	case 3122:
-	  return 1115.683E-3;
-	  break;
-    case 12:
-    case 14:
-    case 16:
-      return 0.0;
-      break;
-    default:
-      MACH3LOG_ERROR("Haven't got a saved mass for PDG: {}", PDG);
-      MACH3LOG_ERROR("Please implement me!");
-      throw MaCh3Exception(__FILE__, __LINE__);
-    } // End switch
-    MACH3LOG_ERROR("Warning, didn't catch a saved mass");
-    return 0;
-  }
-
-  int PDGToNuOscillatorFlavour(int NuPdg){
-    int NuOscillatorFlavour = _BAD_INT_;
-    switch(std::abs(NuPdg)){
-      case NuPDG::kNue:
-        NuOscillatorFlavour = NuOscillator::kElectron;
-        break;
-      case NuPDG::kNumu:
-        NuOscillatorFlavour = NuOscillator::kMuon;
-        break;
-      case NuPDG::kNutau:
-        NuOscillatorFlavour = NuOscillator::kTau;
-        break;
-      default:
-        MACH3LOG_ERROR("Unknown Nuetrino PDG {}, cannot convert to NuOscillator type", NuPdg);
-        break;
-    }
-
-    //This is very cheeky but if the PDG is negative then multiply the PDG by -1
-    // This is consistent with the treatment that NuOscillator expects as enums only
-    // exist for the generic matter flavour and not the anti-matter version
-    if(NuPdg < 0){NuOscillatorFlavour *= -1;}
-
-    return NuOscillatorFlavour;
-  }
-
-
-  //DB Anything added here must be of the form 2^X, where X is an integer
-  //
-  //DB Used to contain which DetIDs are supported
-  std::unordered_map<int,int>KnownDetIDsMap({
-      {0,1},    //ND
-      {1,8},    //FD
-      {2,16},   //SK1Rmu
-      {3,32},   //Nova
-      {4,64},   //Atm SubGeV e-like
-      {5,128},  //Atm SubGeV mu-like 
-      {6,256},  //Atm MultiGeV e-like
-      {7,512},  //Atm MultiGeV mu-like
-	  });
-
-  int nKnownDetIDs = int(KnownDetIDsMap.size());
-
-}
+#include "samplePDF/HistogramUtils.h"
 
 // **************************************************
 //KS: ROOT changes something with binning when moving from ROOT 5 to ROOT 6. If you open ROOT5 produced file with ROOT6 you will be missing 9 last bins
@@ -448,7 +331,6 @@ double PolyIntegralWidth(TH2Poly *Histogram) {
 //KS: Remove fitted TF1 from hist to make comparison easier
 void RemoveFitter(TH1D* hist, const std::string& name) {
 // *********************
-
   TList *listOfFunctions = hist->GetListOfFunctions();
   TF1 *fitter = dynamic_cast<TF1*>(listOfFunctions->FindObject(name.c_str()));
 

samplePDF/HistogramUtils.h

@@ -0,0 +1,72 @@
+#pragma once
+
+// ROOT include
+#include "TGraphAsymmErrors.h"
+#include "TH2Poly.h"
+#include "TLorentzVector.h"
+#include "TObjString.h"
+
+// MaCh3 inlcudes
+#include "samplePDF/Structs.h"
+
+/// @brief WP: Helper function for calculating unbinned Integral of TH2Poly i.e including overflow
+double OverflowIntegral(TH2Poly* poly);
+
+/// @brief WP: Helper function for calculating binned Integral of TH2Poly i.e not including overflow
+double NoOverflowIntegral(TH2Poly* poly);
+
+/// @brief WP: Poly Projectors
+TH1D* PolyProjectionX(TObject* poly, std::string TempName, const std::vector<double>& xbins, const bool computeErrors = false);
+/// @brief WP: Poly Projectors
+TH1D* PolyProjectionY(TObject* poly, std::string TempName, const std::vector<double>& ybins, const bool computeErrors = false);
+
+/// @brief KS: Convert TH2D to TH2Poly
+TH2D* ConvertTH2PolyToTH2D(TH2Poly *poly, TH2D *TH2Dhist);
+/// @brief KS: Convert TH2Poly to TH2D
+TH2Poly* ConvertTH2DtoTH2Poly(TH2D *TH2Dhist);
+
+/// @brief WP: Helper to Normalise histograms
+TH2Poly* NormalisePoly(TH2Poly* Histogram);
+
+/// @brief WP: Helper to scale th2poly analogous to th2d scale with option "width"
+TH2Poly* PolyScaleWidth(TH2Poly *Histogram, double scale);
+/// @brief WP: Helper to calc integral of th2poly analogous to th2d integra; with option "width"
+double PolyIntegralWidth(TH2Poly *Histogram);
+
+/// @brief KS: ROOT changes something with binning when moving from ROOT 5 to ROOT 6. If you open ROOT5 produced file with ROOT6 you will be missing 9 last bins
+/// @brief However if you use ROOT6 and have ROOT6 file exactly the same code will work. Something have changed with how TH2Poly bins are stored in TFile
+/// @param file ROOT file that we will make version checks
+void CheckTH2PolyFileVersion(TFile *file);
+
+/// @brief KS: Remove fitted TF1 from hist to make comparison easier
+void RemoveFitter(TH1D* hist, const std::string& name);
+
+/// @brief Used by sigma variation, check how 1 sigma changes spectra
+/// @param sigmaArrayLeft sigma var hist at -1 or -3 sigma shift
+/// @param sigmaArrayCentr sigma var hist at prior values
+/// @param sigmaArrayRight sigma var hist at +1 or +3 sigma shift
+/// @param title A tittle for returned object
+/// @return A `TGraphAsymmErrors` object that visualizes the sigma variation of spectra, showing confidence intervals between different sigma shifts.
+TGraphAsymmErrors* MakeAsymGraph(TH1D* sigmaArrayLeft, TH1D* sigmaArrayCentr, TH1D* sigmaArrayRight, const std::string& title);
+
+/// @brief Helper to check if files exist or not
+inline std::string file_exists(std::string filename) {
+  std::ifstream infile(filename.c_str());
+  if (!infile.good()) {
+    MACH3LOG_ERROR("*** ERROR ***");
+    MACH3LOG_ERROR("File {} does not exist", filename);
+    MACH3LOG_ERROR("Please try again");
+    MACH3LOG_ERROR("*************");
+    throw MaCh3Exception(__FILE__ , __LINE__ );
+  }
+  return filename;
+}
+
+/// @brief DB Get the Cherenkov momentum threshold in MeV
+double returnCherenkovThresholdMomentum(int PDG);
+
+/// @brief Recalculate Q^2 after Eb shift. Takes in shifted lepton momentum, lepton angle, and true neutrino energy
+double CalculateQ2(double PLep, double PUpd, double EnuTrue, double InitialQ2 = 0.0);
+
+/// @brief Recalculate Enu after Eb shift. Takes in shifted lepton momentum, lepton angle, and binding energy change, and if nu/anu
+double CalculateEnu(double PLep, double cosTheta, double EB, bool neutrino);