Merge pull request #62 from mach3-software/feature_Binned_spline_CI

Feat: Binned spline CI

CIValidations/CovarianceValidations.cpp

@@ -17,7 +17,7 @@ int main(int argc, char *argv[])
   std::vector<std::string> xsecCovMatrixFile = {"Inputs/SystematicModel.yaml"};
   auto xsec = std::make_unique<covarianceXsec>(xsecCovMatrixFile, "xsec_cov");
 
-  std::vector<double> ParProp = {1.05, 0.90, 1.10, 1.05, 0.25, 1.70, 3.20, -1.10, -1.70};
+  std::vector<double> ParProp = {1.05, 0.90, 1.10, 1.05, 1.05, 1.05, 1.05, 1.05, 1.70, 3.20, -1.10, -1.70};
   xsec->setParameters(ParProp);
   xsec->printNominalCurrProp();
 
@@ -148,7 +148,7 @@ std::string yamlContent = R"(
   //KS: Let's make Doctor Wallace proud
   Adapt->initialiseAdaption(AdaptSetting);
 
-  std::vector<double> ParAdapt = {1.05, 0.90, 1.10, 1.05, 0.25, 1.70, 3.20, -1.10, -1.70, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 };
+  std::vector<double> ParAdapt = {1.05, 0.90, 1.10, 1.05, 1.05, 1.05, 1.05, 1.05, 1.70, 3.20, -1.10, -1.70, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0 };
   Adapt->setParameters(ParAdapt);
   bool increase = true;
   for(int i = 0; i < 50000; ++i ) {

CIValidations/LLHValidation.cpp

@@ -20,7 +20,7 @@ int main(int argc, char *argv[])
   }
 
   TFile* file = new TFile("LLH_Test.root", "UPDATE");
-  std::vector<std::string> Names = {"Norm_Param_0", "Norm_Param_1", "Norm_Param_2"};
+  std::vector<std::string> Names = {"Norm_Param_0", "Norm_Param_1", "Norm_Param_2", "BinnedSplineParam1", "BinnedSplineParam2", "BinnedSplineParam3", "BinnedSplineParam4", "BinnedSplineParam5"};
   std::ofstream outFile("LLH.txt");
   for (size_t i = 0; i < Names.size(); ++i) {
     std::string histPath = "Sample_LLH/" + Names[i] + "_sam";

CIValidations/NuOscillatorInterfaceValidations.cpp

@@ -1,3 +1,4 @@
+#include <Constants/OscillatorConstants.h>
 #include "Utils/Comparison.h"
 
 #include "Oscillator/OscillatorFactory.h"
@@ -25,7 +26,7 @@ int main() {
   auto osc = OscillatorFactory().CreateOscillator(
       "NuOscillatorInterfaceValidations.yml");
 
-  std::vector<double> Energies;
+  std::vector<FLOAT_T> Energies;
   double start = -3;
   double end = 3;
   size_t nbins = 100;
@@ -37,7 +38,7 @@ int main() {
   osc->SetEnergyArrayInCalcer(Energies);
   osc->Setup();
 
-  std::vector<double> OscParProp = {0.3,      0.5, 0.020, 7.53e-5,
+  std::vector<FLOAT_T> OscParProp = {0.3,      0.5, 0.020, 7.53e-5,
                                     2.494e-3, 0.0, 1300,  2.6, 0.5};
   osc->CalculateProbabilities(OscParProp);
 

CIValidations/Scripts/CreateBinnedSplineFile.cpp

@@ -0,0 +1,79 @@
+#include "TRandom3.h"
+#include "TSpline.h"
+#include "TFile.h"
+#include "TGraph.h"
+#include "TString.h"
+#include "TH3F.h"
+
+#include <iostream>
+#include <algorithm>
+
+void CreateBinnedSplineFile(){
+
+  int nTrueEnergyBins = 5;
+  int nXBins = 10;
+  //To make splines in Etrue, x and y change this to be a larger number
+  int nYBins = 1;
+  double TrueEnergyBins[] = {0., 1.0, 2.0, 3.0, 4.0, 5.0};
+  double XBins[] = {0., 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0};
+  double YBins[] = {0.,10.0};
+
+  //Let's try out 5 systematics
+  //This gives us a total of 5*10*10*5 = 2500 splines
+  std::vector<int> SystematicKnots = {5, 7, 7, 2, 3};//, 7, 2, 7, 7, 7};
+  std::vector<std::vector<std::string>> SystematicModeNames = {{"ccqe"}, {"cc1pipm"}, {"2p2h"}, {"other"}, {"ccqe"}};
+  //Values to centre the gaussian response on
+  std::vector<double> SystematicMeanResponse = {1.0, 2.0, 1.0, 1.5, 1.0};
+  std::vector<int> SystematicNominalKnot = {3, 4, 4, 1, 2};
+  std::vector<std::string> SystematicNames = {"mysyst1", "mysyst2", "mysyst3","mysyst4","mysyst5"};
+
+  if(SystematicKnots.size() != SystematicNames.size() && SystematicNames.size() != SystematicNominalKnot.size() && SystematicKnots.size() != SystematicMeanResponse.size()) {
+    std::cerr << "Systematic vectors are not of equal length this might cause problems, please fix!" << std::endl;
+    throw;
+  }
+
+  TRandom3 Random = TRandom3(0);
+
+  auto OutputFile = std::unique_ptr<TFile>(TFile::Open("BinnedSplinesTutorialInputs.root", "RECREATE"));
+
+  //Firstly let's write the dev.tmp.0.0 file which sets the binning
+  TH3F* BinningHist = new TH3F("dev_tmp.0.0", "dev_tmp.0.0", nTrueEnergyBins, TrueEnergyBins, nXBins, XBins, nYBins, YBins);
+
+  for(auto iSyst = 0 ; iSyst < SystematicNames.size() ; ++iSyst){
+    for(auto SystematicModeName : SystematicModeNames[iSyst]){
+      for(auto SystematicKnotNumber : SystematicKnots) {
+        for(auto TrueEnergyBin_i = 0 ; TrueEnergyBin_i < nTrueEnergyBins ; ++TrueEnergyBin_i){
+          for(auto XBin_i = 0 ; XBin_i < nXBins ; ++XBin_i){
+            for(auto YBin_i = 0 ; YBin_i < nYBins ; ++YBin_i){
+              double knot_w = 1;
+              TGraph *graph = new TGraph(SystematicKnotNumber);
+
+              for(auto iKnot = 0 ; iKnot < SystematicKnotNumber ; ++iKnot) { 
+                //Check on if you are the nominal knot
+                if(iKnot != SystematicNominalKnot[iSyst]) {
+                  //Do a random throw from a gaussian 
+                  knot_w = Random.Gaus(SystematicMeanResponse[iSyst], 1.0);
+                  knot_w = std::max(0., knot_w);
+                  // std::cout << "Knot_w is " << knot_w << std::endl;
+                  //point number, x-val, y-val
+                }
+                graph->SetPoint(iKnot, iKnot, knot_w);
+              }
+              TSpline3 *Spline = new TSpline3(Form("dev.%s.%s.sp.%i.%i.%i", SystematicNames[iSyst].c_str(), SystematicModeName.c_str(), TrueEnergyBin_i, XBin_i, YBin_i),graph);
+              Spline->SetName(Form("dev.%s.%s.sp.%i.%i.%i", SystematicNames[iSyst].c_str(), SystematicModeName.c_str(), TrueEnergyBin_i, XBin_i, YBin_i));
+              //This makes things slow but removes many backup-cycles being saved to the file
+              Spline->Write(Spline->GetName(), TDirectoryFile::kOverwrite);
+              delete graph;
+              delete Spline;
+            }
+          }
+        }
+      }
+    }
+  }
+
+  OutputFile->Write();
+  OutputFile->Close();
+
+  return;
+}