Feature/dbarrow257/osc prob dune old core

configs/Variations_Atmospherics.yaml

@@ -0,0 +1,120 @@
+General:
+  OutputFile: "parameter_study_sterile/test_2d_var/test_var_sterile_all param_thetav1_AllMC_trueE_trueCos.root"
+
+  DUNESamples: ["configs/Samples/AtmSample_AllMC.yaml"]
+
+  #Nu-FIT
+  #OscillationParameters: [0.310, 0.582, 0.224, 7.39E-5, 2.5254E-3, -2.498, 25]
+
+  #Nu-FIT 5.2 w. SK atm
+  #OscillationParameters: [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233]
+  OscillationParameters: [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, 0., 0., 0., 0.001, 0., 0.]
+  #OscillationParameters: [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, -0.5, 0.2, 0.4, 0.3, 0.5, 0.25, 0., 0., 0., 0.2, 0.4, 0.4]
+
+  #T2K-like best-fit
+  #OscillationParameters: [0.307, 0.528, 0.0218, 7.53e-5, 2.509e-3, -1.601, 25]
+
+  Variations:
+    - Parameter: 
+        Name: sin2th_14
+        #OscParDefault: [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, 0., 0., 0., 0.001, 0., 0.]
+        VarValues: [0., 0.001, 0.005, 0.01, 0.05, 0.1, 0.5]
+
+    - Parameter: 
+        Name: sin2th_24
+        #OscParDefault: [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, 0., 0., 0., 0.001, 0., 0.]
+        VarValues: [0., 0.001, 0.005, 0.01, 0.05, 0.1, 0.5]
+
+    - Parameter: 
+        Name: sin2th_34
+        #OscParDefault: [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, 0., 0., 0., 0.001, 0., 0.]
+        VarValues: [0., 0.001, 0.005, 0.01, 0.05, 0.1, 0.5]
+
+    - Parameter:
+        Name: delm2_14
+        OscParDefault: [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, 0.03, 0.06, 0.05, 0.001, 0., 0.]
+        VarValues: [0., 0.00001, 0.0001, 0.001, 0.01, 0.1, 1.]
+
+    - Parameter:
+        Name: delta_14
+        OscParDefault: [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, 0.03, 0.06, 0.05, 0.001, 0., 0.]
+        VarValues: [-2.35619449019, -1.57079632679, -0.78539816339, 0., 0.78539816339, 1.57079632679, 2.35619449019, 3.14159265]
+
+    - Parameter:
+        Name: delta_24
+        OscParDefault: [0.303, 0.451, 0.02225, 7.41E-5, 2.45E-3, -2.233, 0.03, 0.06, 0.05, 0.001, 0., 0.]
+        VarValues: [-2.35619449019, -1.57079632679, -0.78539816339, 0., 0.78539816339, 1.57079632679, 2.35619449019, 3.14159265]
+
+  Systematics:
+    XsecCovFile: ["configs/CovObjs/xsec_covariance_DUNE_systs_2022a_FD_v3.yaml"]
+    XsecCovName: "xsec_cov"
+    XsecStepScale: 0.1
+    XsecAtGen: false 
+    OscCovFile: ["configs/CovObjs/OscProb_test.yaml"]
+    OscCovName: "osc_cov"
+
+  Fitter:
+    FitTestLikelihood: false
+
+  MCMC:
+    NSteps: 2000
+    AutoSave: 10000
+
+  Output:
+    FileName: "TestEventRates.root"
+    OUTPUTNAME: "TestLLH.root"
+
+  ProcessMCMC: No
+  Seed: 0
+  Debug: No
+
+"Projections": [
+  {
+    "Name": "RecoNuEnergy", 
+    "VarString": "RecoNeutrinoEnergy",
+    "VarBins": [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4],
+    "KinematicCuts": [
+      {
+        "Range": [0.1,0.8],
+        "Name": "TrueNuEnergy", 
+        "VarString": "TrueNeutrinoEnergy"
+      }, 
+      {
+        "Range": [0.4,1.0],
+        "Name": "TrueCosZ", 
+        "VarString": "TrueCosineZ"
+      }
+    ], 
+    "CategoryCuts": [
+      {
+        "Name": "OscillationChannel_Single", 
+        "VarString": "OscChannel",
+        "Breakdown": [[0.0],  [1.0],  [2.0],  [3.0],  [4.0],  [5.0],  [6.0],  [7.0],  [8.0],  [9.0], [10.0], [11.0]],
+	"Names":     ["nue_x_nue","nue_x_numu","nue_x_nutau","numu_x_nue","numu_x_numu","numu_x_nutau","nuebar_x_nuebar","nuebar_x_numubar","nuebar_x_nutaubar","numubar_x_nuebar","numubar_x_numubar","numubar_x_nutaubar"],
+      },
+      {
+        "Name": "OscillationChannel_Group", 
+        "VarString": "OscChannel",
+        "Breakdown": [[0.0, 1.0, 2.0, 3.0, 4.0, 5.0],  [6.0, 7.0, 8.0, 9.0, 10.0, 11.0]],
+	"Names":     ["Nu","Nubar"],
+      },
+      {
+        "Name": "Mode_Single", 
+        "VarString": "Mode",
+        "Breakdown": [[0.0] , [1.0], [2.0], [3.0], [4.0], [5.0], [6.0], [7.0], [8.0], [9.0], [10.0], [11.0], [12.0], [13.0], [14.0], [15.0], [16.0], [17.0], [18.0], [19.0], [20.0], [21.0], [22.0], [23.0], [24.0], [25.0], [26.0]],
+      },
+      {
+        "Name": "Mode_Group", 
+        "VarString": "Mode",
+        "Breakdown": [ [0.0],   [2.0],   [3.0],   [9.0],  [15.0],  [16.0], [1.0, 4.0, 5.0, 6.0, 7.0, 8.0, 10.0, 11.0, 12.0, 13.0], [14.0, 17.0, 18.0, 19.0, 20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0]],
+	"Names":     ["CCQE", "CCDIS", "CCRES", "CCMEC", "NCDIS", "NCRES",                                                "CCOth",                                                            "NCOth"],
+      }
+    ],
+  },
+
+  {
+    "Name": "TrueNuEnergy",
+    "VarString": "TrueNeutrinoEnergy",
+    "VarBins": [0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3.0,3.5,4.0,4.5,5.0,6.0,7.0,8.0,9.0,10.0],
+  }
+]

src/CMakeLists.txt

@@ -6,6 +6,7 @@ foreach(app
 	Projections
 	SigmaVariation
 	LikelihoodScan
+  Variations
       )
 
     add_executable(${app} ${app}.cpp)

src/LikelihoodScan.cpp

@@ -45,7 +45,7 @@ int main(int argc, char * argv[]) {
 
     MACH3LOG_INFO("Event rate for {} : {:<5.2f}", Sample->GetName(), Sample->get1DHist()->Integral());
     if (Sample->GetNDim() == 1)
-      Sample->addData((TH2D*)DUNEHists.back());
+      Sample->addData((TH1D*)DUNEHists.back());
     else if (Sample->GetNDim() == 2)
       Sample->addData((TH2D*)DUNEHists.back());
   }

src/Variations.cpp

@@ -0,0 +1,130 @@
+#include <iostream>
+#include <chrono>
+#include <iomanip>
+#include <vector>
+
+#include <TH1D.h>
+#include <TH2D.h>
+#include <THStack.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TRint.h>
+#include <TLegend.h>
+#include <TColor.h>
+#include <TMath.h>
+
+#include "mcmc/mcmc.h"
+#include "samplePDFDUNE/MaCh3DUNEFactory.h"
+#include "samplePDFDUNE/StructsDUNE.h"
+
+int main(int argc, char * argv[]) {
+  if(argc == 1){
+    MACH3LOG_ERROR("Usage: bin/EventRatesDUNEBeam config.cfg");
+    return 1;
+  }
+  manager* FitManager = new manager(argv[1]);
+
+
+  //###############################################################################################################################
+  //Create samplePDFFD objects
+
+  covarianceXsec* xsec = nullptr;
+  covarianceOsc* osc = nullptr;
+
+  std::vector<samplePDFFDBase*> DUNEPdfs;
+  MakeMaCh3DuneInstance(FitManager, DUNEPdfs, xsec, osc);
+
+  std::vector<double> oscpars = FitManager->raw()["General"]["OscillationParameters"].as<std::vector<double>>();
+
+  //###############################################################################################################################
+  //Perform reweight and print total integral
+
+  MACH3LOG_INFO("=======================================================");
+  for(samplePDFFDBase* Sample: DUNEPdfs){
+    Sample->reweight();
+    MACH3LOG_INFO("Event rate for {} : {:<5.2f}", Sample->GetName(), Sample->get1DHist()->Integral());
+  }
+
+  //###############################################################################################################################
+  //DB Can't use the core sigma variations as it's entirely set up around the concept of multiple selections per samplePDF object
+  //   Thats not the case in the FD code, which has one selection per samplePDF object
+  //   Consequently have to write out own code
+
+  std::vector<covarianceBase*> CovObjs;
+  CovObjs.emplace_back(xsec);
+  CovObjs.emplace_back(osc);
+
+  MACH3LOG_INFO("=======================================================");
+
+  std::string OutputFileName = FitManager->raw()["General"]["OutputFile"].as<std::string>();
+  TFile* File = TFile::Open(OutputFileName.c_str(),"RECREATE");
+
+
+  for (covarianceBase* CovObj: CovObjs) {
+    MACH3LOG_INFO("Starting Variations for covarianceBase Object: {}",CovObj->getName());
+
+    int nPars = CovObj->getNpars();
+
+    for (int iPar=0;iPar<nPars;iPar++) {
+      std::string ParName = CovObj->GetParName(iPar);
+
+      for (auto const &param : FitManager->raw()["General"]["Variations"]) {
+
+        std::string VarName = (param["Parameter"]["Name"].as<std::string>());
+
+        if(ParName == VarName) {
+
+          MACH3LOG_INFO("\tParameter : {:<30}",ParName);
+
+          if(!param["Parameter"]["OscParDefault"]){ //if specific default values not specified for the parameter then use global default ones
+            CovObj->setParameters(oscpars);
+          }
+          else {
+            CovObj->setParameters((param["Parameter"]["OscParDefault"].as<std::vector<double>>()));
+          }
+
+          std::vector<double> valVariations = (param["Parameter"]["VarValues"].as<std::vector<double>>());
+
+          File->cd();
+          File->mkdir(ParName.c_str());
+          File->cd(ParName.c_str());
+
+          for (size_t iSigVar=0;iSigVar<valVariations.size();iSigVar++) {
+	          double VarVal = valVariations[iSigVar];
+
+	          MACH3LOG_INFO("\t\tParameter Value : {:<10.7f}",VarVal);
+            CovObj->setParProp(iPar,VarVal);
+
+	          for (size_t iSample=0;iSample<DUNEPdfs.size();iSample++) {
+	            std::string SampleName = DUNEPdfs[iSample]->GetName();
+
+	            File->cd(ParName.c_str());
+	            if (iSigVar == 0) {
+	              File->mkdir((ParName+"/"+SampleName).c_str());
+	            }
+	            File->cd((ParName+"/"+SampleName).c_str());
+
+	            DUNEPdfs[iSample]->reweight();
+              TH1* Hist;
+              if(DUNEPdfs[iSample]->GetNDim() == 1)
+	              Hist = DUNEPdfs[iSample]->get1DHist();
+              else if(DUNEPdfs[iSample]->GetNDim() == 2)
+                Hist = DUNEPdfs[iSample]->get2DHist();
+
+	            MACH3LOG_INFO("\t\t\tSample : {:<30} - Integral : {:<10}",SampleName,Hist->Integral());
+
+	            Hist->Write(Form("Variation_%.2e",VarVal));
+	          }
+          }
+
+        }
+      }
+    }
+
+    MACH3LOG_INFO("=======================================================");
+  }
+
+
+  File->Close();
+
+}