nue_data_mc_validation.C
Go to the documentation of this file.
1 // Make Data and MC Spectra for validation
2 // makeFile = true saves the spectra for ldrData, ldrMC
3 // makeFile = false reruns formatting for generic validation
4 
5 // SA nue cuts and variables
6 // Same macro should work for different datasets,
7 // as long as all vars and cuts are proper
8 
9 // TODO: add logic for FD
10 // TODO: add logic to split data and MC
11 
12 #include "CAFAna/Core/Utilities.h"
15 #include "CAFAna/Core/Spectrum.h"
16 #include "CAFAna/Cuts/SpillCuts.h"
17 #include "CAFAna/Cuts/NueCutsSecondAna.h"
20 #include "CAFAna/Vars/Vars.h"
26 #include "OscLib/IOscCalc.h"
32 
33 using namespace ana;
34 
35 #include "TFile.h"
36 
37 #include <iostream>
38 #include <vector>
39 
41  const std::string predFileName,
42  const bool makeFile=false)
43 {
44  // Variable groups from NueVarsExtra
45  // Variables from NueVars.h are used here or in the cuts
46  assert(beamMode == "rhc" || beamMode == "fhc");
47 
48  std::string ldrData;
49  std::string ldrMC;
50  if(beamMode == "rhc") {
51  ldrData = "prod_caf_R17-09-05-prod4recopreview.f_nd_numi_rhc_full_v1_addShortSimpleCVN_goodruns";
52  ldrMC = "prod_caf_R17-11-14-prod4reco.e_nd_genie_nonswap_rhc_nova_v08_full_v1";
53  }
54  else {
55  ldrData = "prod_caf_R17-09-05-prod4recopreview.f_nd_numi_fhc_full_v1_addShortSimpleCVN_goodruns";
56  ldrMC = "prod_caf_R17-11-14-prod4reco.d_nd_genie_nonswap_fhc_nova_v08_full_v1";
57  }
58 
59  std::vector <HistDef> defs;
60 
61  AddNueHistDefBasic(defs);
62  AddNueHistDefVertexND(defs);
64  AddNueHistDefShower(defs);
65  AddNueHistDefPIDs(defs);
66  AddNueHistDefForCvn(defs);
68  AddNueHistDefProngCVN(defs);
69 
70  // Cut flow. Change the names/cuts here only if new spectra are going to be saved
71  struct CutDef
72  {
74  Cut cut;
75  };
76 
77  std::vector <CutDef> sels = {
78  {"0_Presel",kNue2018NDPresel},
79  {"1_CVNeLow",kNue2018NDCVNSsb},
80  };
81  if(beamMode == "rhc")
82  sels.push_back({"2_CVNeHigh",kNue2018NDPresel && kCVNSSe >= .98});
83  else
84  sels.push_back({"2_CVNeHigh",kNue2018NDPresel && kCVNSSe >= .96});
85 
86  std::string decompFName = "nue2018_decomp_ratios_rhc_prop.root";
87  std::string decompShortName = "prop";
88  if(beamMode == "fhc") {
89  decompFName = "nue2018_decomp_ratios_fhc.root";
90  decompShortName = "combo";
91  }
92 
93  const Var kDecompWeight(DecompWeightFunc(FindCAFAnaDir()+
94  "/data/decomp/" + decompFName,
95  decompShortName,
97 
98  if(makeFile){
99  SpectrumLoader loaderData(ldrData.c_str());
100  SpectrumLoader loaderMC(ldrMC.c_str());
101 
102  loaderData.SetSpillCut(kStandardSpillCuts);
103  loaderMC.SetSpillCut(kStandardDQCuts);
104 
105  std::vector < std::vector<Spectrum*> > spects;
106  std::vector < std::vector<IPrediction*> > preds_mc;
107  std::vector < std::vector<IPrediction*> > preds_mc_decomp;
108 
109  for(int selIdx = 0; selIdx<(int)sels.size(); ++selIdx){
110  std::vector<Spectrum *> temp_spec;
111  std::vector<IPrediction *> temp_pred_mc;
112  std::vector<IPrediction *> temp_pred_mc_decomp;
113 
114  for(int varIdx = 0; varIdx < (int)defs.size(); ++varIdx){
115  const HistAxis& axis = defs[varIdx].axis;
116  temp_spec.emplace_back( new Spectrum(loaderData, axis, sels[selIdx].cut, kNoShift));
117  temp_pred_mc.emplace_back(new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
118  axis.GetLabels()[0], axis.GetBinnings()[0], axis.GetVars()[0],
119  sels[selIdx].cut, kNoShift, kXSecCVWgt2018*kPPFXFluxCVWgt)); //TODO Broken w new cafs
120  temp_pred_mc_decomp.emplace_back(new PredictionNoExtrap(loaderMC, kNullLoader, kNullLoader,
121  axis.GetLabels()[0], axis.GetBinnings()[0],
122  axis.GetVars()[0], sels[selIdx].cut,
123  kNoShift,
124  kXSecCVWgt2018*kPPFXFluxCVWgt*kDecompWeight));
125  }
126  spects.push_back(temp_spec);
127  preds_mc.push_back(temp_pred_mc);
128  preds_mc_decomp.push_back(temp_pred_mc_decomp);
129  }
130 
131  loaderMC.Go();
132  loaderData.Go();
133 
134  TFile* fout = new TFile(predFileName.c_str(), "RECREATE");
135  for(int selIdx = 0; selIdx < (int)sels.size(); ++selIdx){
136  const char* cut = sels[selIdx].name.c_str();
137  std::cout << cut << std::endl;
138  for(int varIdx = 0; varIdx < (int)defs.size(); ++varIdx){
139  const char* name = defs[varIdx].name.c_str();
140  std::string tempname = defs[varIdx].name;
141  spects[selIdx][varIdx]->SaveTo(fout, TString::Format("spect_%s_%s", name, cut));
142  preds_mc[selIdx][varIdx]->SaveTo(fout, TString::Format("pred_mc_%s_%s", name, cut));
143  preds_mc_decomp[selIdx][varIdx]->SaveTo(fout, TString::Format("pred_mc_decomp_%s_%s", name, cut));
144 
145  }
146  }
147  delete fout;
148  return;
149  } //end makeFile
150  else{
151  // Play with the format for the generic validation scripts
152  // Here I'm taking the one file and splitting data and MC
153  // Data and MC will be treated as "datasets" = separate files
154  // Alternatively, you could treat them as folders, or categories in a group
155  // Your mileage may vary
156 
157  TFile * predFile = new TFile (predFileName.c_str(), "READ");
158  TString validFileName1 (predFileName); validFileName1.ReplaceAll(".root","_mc_validation.root");
159  TString validFileName3 (predFileName); validFileName3.ReplaceAll(".root","_data_validation.root");
160  TString validFileName2 (predFileName); validFileName2.ReplaceAll(".root","_mc_decomp_validation.root");
161 
162  TFile * validFile1 = new TFile(validFileName1,"RECREATE");
163  TFile * validFile3 = new TFile(validFileName3,"RECREATE");
164  TFile * validFile2 = new TFile(validFileName2,"RECREATE");
165 
166  int kNumSels = sels.size();
167  int kNumVars = defs.size();
168 
169  // Don't make all the canvases and .json right now
170  // kNumSels = 4;
171  // kNumVars = 2;
172 
173  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
174  const char* cut = sels[selIdx].name.c_str();
175  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
176  const char* name = defs[varIdx].name.c_str();
177 
178  auto spec = Spectrum::LoadFrom(predFile, TString::Format("spect_%s_%s", name, cut));
179  auto pred_mc = PredictionNoExtrap::LoadFrom (predFile,TString::Format("pred_mc_%s_%s", name, cut).Data());
180  auto pred_mc_decomp = PredictionNoExtrap::LoadFrom(predFile, TString::Format("pred_mc_decomp_%s_%s", name, cut).Data());
181 
182 
183  // Subfolders -> new pages (not shown here)
184  // Data and MC for variable "name" will be plotted together
185  // Plots will have selector for cut level {group=Cut,cat=%s}
186  // Plots will have selector for MC component component; total MC = all_sum will be created as well
187  // Add functions to operate on variable, cut or axis labels
188 
189  spec->SaveTo(validFile3, TString::Format("%s{group=Cut,cat=%s}", name,cut));
190 
191  pred_mc->ComponentCC(12,12). Unoscillated().SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=0_nue}", name,cut));
192  pred_mc->ComponentCC(14,14). Unoscillated().SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=1_numu}", name,cut));
193  pred_mc->ComponentCC(-12,-12).Unoscillated().SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=3_antinue}", name,cut));
194  pred_mc->ComponentCC(-14,-14).Unoscillated().SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=4_antinumu}", name,cut));
195  pred_mc->ComponentNCTotal(). SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=2_nc}", name,cut));
196 
197 
198  pred_mc_decomp->ComponentCC(12,12). Unoscillated().SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=0_nue}", name,cut));
199  pred_mc_decomp->ComponentCC(14,14). Unoscillated().SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=1_numu}", name,cut));
200  pred_mc_decomp->ComponentCC(-12,-12).Unoscillated().SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=3_antinue}", name,cut));
201  pred_mc_decomp->ComponentCC(-14,-14).Unoscillated().SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=4_antinumu}", name,cut));
202  pred_mc_decomp->ComponentNCTotal(). SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=2_nc}", name,cut));
203 
204  }
205  }
206  std::cout << "Formatted outputs saved to " << validFileName1 << ", " << validFileName2 << " and " << validFileName3 << std::endl;
207  delete validFile1;
208  delete validFile3;
209 
210  }
211 }
212 
const XML_Char * name
Definition: expat.h:151
const int kNumVars
Definition: vars.h:14
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const std::vector< T > & GetVars() const
Definition: HistAxis.h:92
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const TString name
void AddNueHistDefBasic(std::vector< HistDef > &hd)
const Cut kNue2018NDCVNSsb
Definition: NueCuts2018.h:153
void SetSpillCut(const SpillCut &cut)
std::string FindCAFAnaDir()
Definition: Utilities.cxx:204
const Cut kNue2018NDPresel
Definition: NueCuts2018.h:145
const Var kNue2018AnaBin([](const caf::SRProxy *sr){int selBin=kNue2018SelectionBin(sr);float nuE=kNueEnergy2018(sr);int nuEBin=nuE/0.5;assert(nuEBin<=8 &&"An event with nuE > 4.5 should never happen");int anaBin=9 *selBin+nuEBin;return anaBin;})
Use this Analysis Binning for Ana2018, official Binning.
Definition: NueCuts2018.h:171
void AddNueHistDefShowerND(std::vector< ana::HistDef > &hd)
void AddNueHistDefSelectionExtras(std::vector< HistDef > &hd)
const Cut sels[kNumSels]
Definition: vars.h:44
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
const SpillCut kStandardDQCuts([](const caf::SRSpillProxy *spill){if(spill->dcmedgematchfrac==0 &&spill->fracdcm3hits==0 &&spill->nmissingdcmslg==0) return bool(spill->isgoodspill); if(spill->det==caf::kNEARDET && (spill->fracdcm3hits > 0.45|| spill->nmissingdcms > 0)) return false; if(spill->eventincomplete) return false; if(spill->det==caf::kFARDET && spill->nmissingdcmslg > 0) return false; if(spill->det==caf::kFARDET && !spill->ismc && spill->dcmedgematchfrac<=0.2) return false;return true;})
Cut out events with a noisy detector or with parts missing.
Definition: SpillCuts.h:16
Struct to hold cut information.
const HistDef defs[kNumVars]
Definition: vars.h:15
const int kNumSels
Definition: vars.h:43
std::vector< float > Spectrum
Definition: Constants.h:610
const Var kCVNSSe([](const caf::SRProxy *sr){throw std::runtime_error("kCVNSSe is no longer available. Fix your macro so you don't use it.");return-5.;})
2018 nue PID
Definition: Vars.h:52
void AddNueHistDefProngCVN(std::vector< HistDef > &hd)
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
void AddNueHistDefShower(std::vector< HistDef > &hd)
const std::vector< Binning > & GetBinnings() const
Definition: LabelsAndBins.h:69
const Cut cut
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
assert(nhit_max >=nhit_nbins)
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
Prediction that just uses FD MC, with no extrapolation.
void AddNueHistDefVertexND(std::vector< HistDef > &hd)
static std::unique_ptr< PredictionNoExtrap > LoadFrom(TDirectory *dir, const std::string &name)
void AddNueHistDefPIDs(std::vector< HistDef > &hd)
void nue_data_mc_validation(const std::string beamMode, const std::string predFileName, const bool makeFile=false)
static NullLoader kNullLoader
Dummy loader that doesn&#39;t load any files.
const std::vector< std::string > & GetLabels() const
Definition: LabelsAndBins.h:68
enum BeamMode string