Draft: Picker24/feature/globalbin vectors kineparamtidy

CMakeLists.txt

@@ -45,10 +45,6 @@ if(NOT TARGET ROOT::ROOT)
   cmessage(FATAL_ERROR "MaCh3 Expected dependency target: ROOT::ROOT")
 endif()
 
-if(DEFINED ROOT_CXX_STANDARD AND ROOT_CXX_STANDARD GREATER CMAKE_CXX_STANDARD)
-  set(CMAKE_CXX_STANDARD ${ROOT_CXX_STANDARD})
-endif()
-
 ############################  DUNEAnaObj  ####################################
 
 find_package(duneanaobj)
@@ -74,7 +70,7 @@ find_package(MaCh3)
 if(NOT MaCh3_FOUND)
   CPMFindPackage(
       NAME MaCh3
-      GIT_TAG "develop"
+      GIT_TAG "picker24/feature/GenericBinningAndReturnKinematicParameterTidy"
       GITHUB_REPOSITORY mach3-software/MaCh3
   )
 else()
@@ -89,7 +85,7 @@ endif()
 
 ############################  C++ Compiler  ####################################
 if (NOT DEFINED CMAKE_CXX_STANDARD OR "${CMAKE_CXX_STANDARD} " STREQUAL " ")
-  SET(CMAKE_CXX_STANDARD 11)
+  SET(CMAKE_CXX_STANDARD 17)
 endif()
 
 if(DEFINED ROOT_CXX_STANDARD AND ROOT_CXX_STANDARD GREATER CMAKE_CXX_STANDARD)

configs/EventRates_Beam.yaml

@@ -1,11 +1,17 @@
 General:
   OutputFile: "DuneEventRates.root"
-  DUNESamples: ["configs/Samples/SamplePDFDune_FHC_numuselec.yaml", "configs/Samples/SamplePDFDune_FHC_nueselec.yaml", "configs/Samples/SamplePDFDune_RHC_numuselec.yaml", "configs/Samples/SamplePDFDune_RHC_nueselec.yaml"]	
+  DUNESamples:
+    - "configs/Samples/SamplePDFDune_FHC_numuselec.yaml" 
+    # - "configs/Samples/SamplePDFDune_FHC_nueselec.yaml"
+    # - "configs/Samples/SamplePDFDune_RHC_numuselec.yaml"
+    # - "configs/Samples/SamplePDFDune_RHC_nueselec.yaml"
+
   #Nu-FIT
   #OscillationParameters: [0.310, 0.582, 0.224, 7.39E-5, 2.5254E-3, -2.498]
   # T2K-like best-fit
   OscillationParameters: [0.307, 0.528, 0.0218, 7.53e-5, 2.509e-3, -1.601, 1284.9, 2.848]
   OscillatorConfigName: "configs/OscillatorObj.yaml"
+
   Systematics:
     XsecCovFile: ["configs/CovObjs/xsec_covariance_DUNE_systs_2022a_FD_v3_xsec.yaml"]
     XsecCovName: "xsec_cov"
@@ -22,6 +28,7 @@ General:
   Output:
     FileName: "TestEventRates.root"
     OUTPUTNAME: "TestLLH.root"
+
   ProcessMCMC: No
   Seed: 0
   Debug: No

configs/Samples/SamplePDFDune_FHC_numuselec.yaml

@@ -10,101 +10,124 @@ SelectionCuts:
     Bounds: [ 50.0, 1244.0 ]
   - KinematicStr: "CVNNumu"
     Bounds: [ 0.5, 999 ]
+
 Binning:
-  XVarStr : "RecoNeutrinoEnergy" 
-  XVarBins: [0.,  0.5,  1.,  1.25, 1.5, 1.75, 2., 2.25, 2.5, 2.75, 3., 3.25, 3.5, 3.75, 4., 5., 6., 10.]
+  Axes:
+# Appropriate values for VarStr
+# TrueNeutrinoEnergy
+# RecoNeutrinoEnergy
+# CVNNumu
+# CVNNue
+# q0
+# q3
+# ERecQE
+# ELepRec
+# EHadRec
+    - VarStr: "ELepRec" 
+      Uniform: [25,0,5]
+      Title: "E_{lep.}^{Rec} [GeV]"
+    - VarStr: "EHadRec" 
+      Uniform: [25,0,5]
+      Title: "E_{had.}^{Vis} [GeV]"
+    - VarStr: "RecoNeutrinoEnergy" 
+      Uniform: [10,0.5,5.5]
+      Title: "E_{#nu}^{Rec} [GeV]"
+  # if you want to use the generic binning with a single 
+  # axis, useful for calculated VarStr that cannot be referenced
+  ForceGeneric: True
+
 DUNESampleBools:
   iselike: no
 InputFiles:
   mtupleprefix: "inputs/DUNE_CAF_files/FD_FHC_ger_"
   mtuplesuffix: "_numuselec.root"
-  splineprefix: "inputs/DUNE_spline_files/FD_FHC_ger_"
-  splinesuffix: "_numuselec_splines.root"
+  # splineprefix: "inputs/DUNE_spline_files/FD_FHC_ger_"
+  # splinesuffix: "_numuselec_splines.root"
 DetID: 24
 NuOsc:
   NuOscConfigFile: "configs/NuOsc/CUDAProb3Linear.yaml"
 NSubSamples: 12
 SubSamples:
   - name: "FHC_numu_x_numu"
     mtuplefile: "numu_x_numu"
-    splinefile: "numu_x_numu"
+    # splinefile: "numu_x_numu"
     samplevecno: 0
     nutype: 14
     oscnutype: 14
     signal: false
   - name: "FHC_nue_x_nue"
     mtuplefile: "nue_x_nue"
-    splinefile: "nue_x_nue"
+    # splinefile: "nue_x_nue"
     samplevecno: 1
     nutype: 12
     oscnutype: 12
     signal: false
   - name: "FHC_numubar_x_numubar"
     mtuplefile: "numubar_x_numubar"
-    splinefile: "numubar_x_numubar"
+    # splinefile: "numubar_x_numubar"
     samplevecno: 2
     nutype: -14
     oscnutype: -14
     signal: false
   - name: "FHC_nuebar_x_nuebar"
     mtuplefile: "nuebar_x_nuebar"
-    splinefile: "nuebar_x_nuebar"
+    # splinefile: "nuebar_x_nuebar"
     samplevecno: 3
     nutype: -12
     oscnutype: -12
     signal: false
   - name: "FHC_numu_x_nue"
     mtuplefile: "numu_x_nue"
-    splinefile: "numu_x_nue"
+    # splinefile: "numu_x_nue"
     samplevecno: 4
     nutype: 14
     oscnutype: 12
     signal: true
   - name: "FHC_numubar_x_nuebar"
     mtuplefile: "numubar_x_nuebar"
-    splinefile: "numubar_x_nuebar"
+    # splinefile: "numubar_x_nuebar"
     samplevecno: 5
     nutype: -14
     oscnutype: -12 
     signal: true
   - name: "FHC_nue_x_numu"
     mtuplefile: "nue_x_numu"
-    splinefile: "nue_x_numu"
+    # splinefile: "nue_x_numu"
     samplevecno: 6
     nutype: 12
     oscnutype: 14
     signal: true
   - name: "FHC_nuebar_x_numubar"
     mtuplefile: "nuebar_x_numubar"
-    splinefile: "nuebar_x_numubar"
+    # splinefile: "nuebar_x_numubar"
     samplevecno: 7
     nutype: -12
     oscnutype: -14
     signal: true
   - name: "FHC_numu_x_nutau"
     mtuplefile: "numu_x_nutau"
-    splinefile: "numu_x_nutau"
+    # splinefile: "numu_x_nutau"
     samplevecno: 8
     nutype: 14
     oscnutype: 16
     signal: true
   - name: "FHC_nue_x_nutau"
     mtuplefile: "nue_x_nutau"
-    splinefile: "nue_x_nutau"
+    # splinefile: "nue_x_nutau"
     samplevecno: 9
     nutype: 12
     oscnutype: 16
     signal: true
   - name: "FHC_numubar_x_nutaubar"
     mtuplefile: "numubar_x_nutaubar"
-    splinefile: "numubar_x_nutaubar"
+    # splinefile: "numubar_x_nutaubar"
     samplevecno: 10
     nutype: -14
     oscnutype: -16
     signal: true
   - name: "FHC_nuebar_x_nutaubar"
     mtuplefile: "nuebar_x_nutaubar"
-    splinefile: "nuebar_x_nutaubar"
+    # splinefile: "nuebar_x_nutaubar"
     samplevecno: 11
     nutype: -12
     oscnutype: -16

samplePDFDUNE/StructsDUNE.h

@@ -6,85 +6,91 @@ struct dunemc_base {
   int oscnutype;
   bool signal; // true if signue
   int nEvents; // how many MC events are there
-  int *Target; //Target the interaction was on
+  std::vector<int> Target; //Target the interaction was on
 
-  double *rw_erec;
-  double *rw_erec_shifted;
-  double *rw_erec_had;
-  double *rw_erec_lep;
-  double *rw_yrec;
-  double *rw_eRecoP;
-  double *rw_eRecoPip;
-  double *rw_eRecoPim;
-  double *rw_eRecoPi0;
-  double *rw_eRecoN;
-
-  double *rw_LepE;
-  double *rw_eP;
-  double *rw_ePip;
-  double *rw_ePim;
-  double *rw_ePi0;
-  double *rw_eN;
-
-  double *rw_etru;
-  double *rw_mom;
-  double *rw_theta;
-  double *rw_Q2;
-
-  double *rw_cvnnumu;
-  double *rw_cvnnue;
-  double *rw_cvnnumu_shifted;
-  double *rw_cvnnue_shifted;
-  int *rw_reco_nue;
-  int *rw_reco_numu;
-  double *rw_berpaacvwgt;
-  int    *rw_isCC;
-  int    *rw_nuPDGunosc;
-  int    *rw_nuPDG;
-  int    *rw_run;
-  bool    *rw_isFHC;
-  double *rw_vtx_x;
-  double *rw_vtx_y;
-  double *rw_vtx_z;
+  std::vector<double> rw_erec;
+  std::vector<double> rw_erec_shifted;
+  std::vector<double> rw_erec_had;
+  std::vector<double> rw_erec_lep;
+  std::vector<double> rw_yrec;
+  std::vector<double> rw_eRecoP;
+  std::vector<double> rw_eRecoPip;
+  std::vector<double> rw_eRecoPim;
+  std::vector<double> rw_eRecoPi0;
+  std::vector<double> rw_eRecoN;
+
+  std::vector<double> true_q0;
+  std::vector<double> true_q3;
+
+  std::vector<double> rw_LepE;
+  std::vector<double> rw_eP;
+  std::vector<double> rw_ePip;
+  std::vector<double> rw_ePim;
+  std::vector<double> rw_ePi0;
+  std::vector<double> rw_eN;
+
+  std::vector<double> rw_etru;
+  std::vector<double> rw_mom;
+  std::vector<double> rw_theta;
+  std::vector<double> rw_Q2;
+
+  std::vector<double> rw_cvnnumu;
+  std::vector<double> rw_cvnnue;
+  std::vector<double> rw_cvnnumu_shifted;
+  std::vector<double> rw_cvnnue_shifted;
+  std::vector<int> rw_reco_nue;
+  std::vector<int> rw_reco_numu;
+  std::vector<double> rw_berpaacvwgt;
+  std::vector<int>    rw_isCC;
+  std::vector<int>    rw_nuPDGunosc;
+  std::vector<int>    rw_nuPDG;
+  std::vector<int>    rw_run;
+  bool *              rw_isFHC;
+  std::vector<double> rw_vtx_x;
+  std::vector<double> rw_vtx_y;
+  std::vector<double> rw_vtx_z;
   double dummy_y;
-  double *rw_reco_q;
-  double *reco_numu;
+  std::vector<double> rw_reco_q;
+  std::vector<double> reco_numu;
 
   double pot_s;
   double norm_s;
 
-  double *beam_w;
-  double *flux_w;
+  std::vector<double> beam_w;
+  std::vector<double> flux_w;
 
-  int *mode;
-  int *isbound;
+  std::vector<int> mode;
+  std::vector<int> isbound;
 
-  double *rw_truecz;
+  std::vector<double> rw_truecz;
 
-  int *nproton; ///< number of (post-FSI) primary protons
-  int *nneutron; ///< number of (post-FSI) primary neutrons
-  int *npip; ///< number of (post-FSI) primary pi+
-  int *npim; ///< number of (post-FSI) primary pi-
-  int *npi0; ///< number of (post-FSI) primary pi0
+  std::vector<int> nproton; ///< number of (post-FSI) primary protons
+  std::vector<int> nneutron; ///< number of (post-FSI) primary neutrons
+  std::vector<int> npip; ///< number of (post-FSI) primary pi+
+  std::vector<int> npim; ///< number of (post-FSI) primary pi-
+  std::vector<int> npi0; ///< number of (post-FSI) primary pi0
 
-  int *ntruemuon; //number of true muons
-  int *ntruemuonprim; //number of true primary muons
-  int *nrecomuon; //number of reconstructed muons
-  double *nmuonsratio; //number of reco muons divided by number of true muons
+  std::vector<int> ntruemuon; //number of true muons
+  std::vector<int> ntruemuonprim; //number of true primary muons
+  std::vector<int> nrecomuon; //number of reconstructed muons
+  std::vector<double> nmuonsratio; //number of reco muons divided by number of true muons
 
-  double *rw_lep_pT;  //transverse lepton momentum
-  double *rw_lep_pZ; //parallel lepton momentum
-  double *rw_reco_vtx_x;
-  double *rw_reco_vtx_y;
-  double *rw_reco_vtx_z;
-  double *rw_reco_rad;
-  double *rw_rad;
+  std::vector<double> rw_lep_pT;  //transverse lepton momentum
+  std::vector<double> rw_lep_pZ; //parallel lepton momentum
+  std::vector<double> rw_reco_vtx_x;
+  std::vector<double> rw_reco_vtx_y;
+  std::vector<double> rw_reco_vtx_z;
+  std::vector<double> rw_reco_rad;
+  std::vector<double> rw_rad;
 
-  double *rw_elep_reco;
-  double *rw_elep_true;
+  std::vector<double> rw_elep_reco;
+  std::vector<double> rw_elep_true;
 
-  int *nrecoparticles;
-  bool *in_fdv;
+  std::vector<int> nrecoparticles;
+  //unfortunately std::vector<bool> is a nightmare, leave this as a pointer
+  bool * in_fdv;
+
+  std::vector<double> global_bin_number; //for using generic binning, this value holds the global bin number for each event
 };
 
 // ********************************