NDDataMC.C
Go to the documentation of this file.
1 
2 // Make Data and MC Spectra for validation
3 // makeFile = true saves the spectra for ldrData, ldrMC
4 // makeFile = false reruns formatting for generic validation
5 
6 // Ana01 nus cuts and variables
7 // Same macro should work for different datasets,
8 // as long as all vars and cuts are proper
9 
11 #include "CAFAna/Core/Progress.h"
12 #include "CAFAna/Core/MultiVar.h"
14 #include "CAFAna/Core/Spectrum.h"
15 #include "CAFAna/Cuts/SpillCuts.h"
16 #include "NuXAna/Cuts/NusCuts17.h"
17 #include "CAFAna/Cuts/TruthCuts.h"
18 #include "CAFAna/Vars/Vars.h"
19 #include "NuXAna/Vars/NusVars.h"
21 #include "OscLib/IOscCalc.h"
24 
25 using namespace ana;
26 
27 #include "TFile.h"
28 
29 #include <iostream>
30 #include <vector>
31 
32 void NDDataMC(const std::string predFileName, const bool makeFile=false,
33  const std::string ldrData="", const std::string ldrMC="")
34 {
35  // Variable groups from NusVars
36  std::vector<HistDef> nus_defs = getAllNusAna01NDHistDefs();
37 
38  // Cut flow. Change the names/cuts here only if new spectra are going to be saved
39  struct CutDef
40  {
42  Cut cut;
43  };
44 
45  std::vector<CutDef > sels = {
46  {"0_NoCut", kNoCut},
47  {"1_NusEventQuality", kNus17EventQuality},
48  {"2_NusFiducial", kNus17EventQuality && kNus17NDFiducial},
49  {"3_NusPresel", kNus17NDPresel},
50  {"4_kNusND", kNus17ND},
51  {"5_kNusNDExtra", kNus17NDExtra},
52  {"6_kNusNDPtp", kNus17NDExtra && (kPartPtp <= 0.8)}
53  };
54 
55  MultiVar kCutFlowVar(
56  [sels](const caf::SRProxy* sr)
57  {
58  std::vector<double> passes;
59  for(unsigned int i = 0; i < sels.size(); ++i){
60  if(sels[i].cut(sr)) passes.push_back(i);
61  }
62  return passes;
63  });
64 
65  const Binning flowBins = Binning::Simple(6, 0, 6);
66  const Var kDummy(
67  [](const caf::SRProxy* sr)
68  {
69  return 1;
70  });
71  nus_defs.insert(nus_defs.begin(), {"_dummy", {"Dummy", Binning::Simple(3, -0.5, 2.5), kDummy}});
72 
73  //Save Data and MC spectra. Prediction -> get MC components
74 
75  if(makeFile){
76  SpectrumLoader loaderData(ldrData.c_str());
77  SpectrumLoader loaderMC(ldrMC.c_str());
78 
79  loaderData.SetSpillCut(kStandardSpillCuts);
80  loaderMC.SetSpillCut(kStandardSpillCuts);
81 
82  std::vector<std::vector<Spectrum*> > spects;
83  std::vector<std::vector<Spectrum*> > spectsMC;
84  std::vector<std::vector<CheatDecomp*> > preds;
85 
86  for(int selIdx = 0; selIdx < (int)sels.size(); ++selIdx){
87  std::vector<Spectrum*> temp_spec;
88  std::vector<Spectrum*> temp_specMC;
89  std::vector<CheatDecomp*> temp_pred;
90 
91  for(int varIdx = 0; varIdx < (int)nus_defs.size(); ++varIdx){
92  const HistAxis& axis = nus_defs[varIdx].axis;
93  temp_spec.emplace_back(new Spectrum(loaderData, axis, sels[selIdx].cut, kNoShift));
94  temp_specMC.emplace_back(new Spectrum(loaderMC, axis, sels[selIdx].cut, kNoShift, kXSecCVWgt2017*kPPFXFluxCVWgt));
95  temp_pred.emplace_back(new CheatDecomp(loaderMC,
96  axis,
97  sels[selIdx].cut,
98  kNoShift,
99  kXSecCVWgt2017*kPPFXFluxCVWgt));
100  }
101  spects.push_back(temp_spec);
102  spectsMC.push_back(temp_specMC);
103  preds.push_back(temp_pred);
104  }
105 
106  auto spfl_data = new Spectrum("Cut Flow", flowBins, loaderData, kCutFlowVar, kNoCut);
107  auto spfl_ndmc = new Spectrum("Cut Flow", flowBins, loaderMC, kCutFlowVar, kNoCut);
108 
109  loaderMC.Go();
110  loaderData.Go();
111 
112  TFile* fout = new TFile(predFileName.c_str(), "RECREATE");
113  for(int selIdx = 0; selIdx < (int)sels.size(); ++selIdx){
114  const char* cut = sels[selIdx].name.c_str();
115  for(int varIdx = 0; varIdx < (int)nus_defs.size(); ++varIdx){
116  const char* name = nus_defs[varIdx].name.c_str();
117  spects[selIdx][varIdx]->SaveTo(fout, TString::Format("spect_%s_%s", name, cut));
118  spectsMC[selIdx][varIdx]->SaveTo(fout, TString::Format("spectMC_%s_%s", name, cut));
119  preds[selIdx][varIdx]->SaveTo(fout, TString::Format("pred_%s_%s", name, cut));
120 
121  }
122  }
123 
124  spfl_data->SaveTo(fout, "spect_data_cutflow");
125  spfl_ndmc->SaveTo(fout, "spect_mc_cutflow");
126  delete fout;
127  return;
128  } //end makeFile
129  else{
130  // Play with the format for the generic validation scripts
131  // Here I'm taking the one file and splitting data and MC
132  // Data and MC will be treated as "datasets" = separate files
133  // Alternatively, you could treat them as folders, or categories in a group
134  // Your mileage may vary
135 
136  TFile * predFile = new TFile (predFileName.c_str(), "READ");
137  TString validFileName1 (predFileName); validFileName1.ReplaceAll(".root","_mc_validation.root");
138  TString validFileName2 (predFileName); validFileName2.ReplaceAll(".root","_data_validation.root");
139 
140  TFile * validFile1 = new TFile(validFileName1,"RECREATE");
141  TFile * validFile2 = new TFile(validFileName2,"RECREATE");
142 
143  int kNumSels = sels.size();
144  int kNumVars = nus_defs.size();
145 
146 
147  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
148  const char* cut = sels[selIdx].name.c_str();
149  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
150  const char* name = nus_defs[varIdx].name.c_str();
151 
152  auto spec = Spectrum::LoadFrom(predFile, TString::Format("spect_%s_%s", name, cut));
153  auto specMC = Spectrum::LoadFrom(predFile, TString::Format("spectMC_%s_%s", name, cut));
154  auto pred = CheatDecomp::LoadFrom(predFile, TString::Format("pred_%s_%s", name, cut));
155 
156 
157  spec->SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}", name,cut));
158  specMC->SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}", name,cut));
159 
160  // 1D Spectra only
161  bool components = true;
162  if(components){
163  TDirectory* comp = validFile1->GetDirectory("MCComponents");
164  if(!comp) comp = validFile1->mkdir("MCComponents");
165 
166  pred->NCTotalComponent(). SaveTo(comp, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=0_nc}", name,cut));
167  pred->NueComponent(). SaveTo(comp, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=1_nue}", name,cut));
168  pred->NumuComponent(). SaveTo(comp, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=2_numu}", name,cut));
169  pred->AntiNueComponent(). SaveTo(comp, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=3_antinue}", name,cut));
170  pred->AntiNumuComponent().SaveTo(comp, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=4_antinumu}", name,cut));
171 
172  (*specMC-(pred->TotalMC())).SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Component,cat=5_other}",name,cut));
173  }
174  } // end var
175  auto spfl = Spectrum::LoadFrom(predFile, "spect_data_cutflow");
176  spfl->SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}","_cutflow",cut));
177 
178  spfl = Spectrum::LoadFrom(predFile, "spect_mc_cutflow");
179  spfl->SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}","_cutflow",cut));
180 
181  } // end cut
182  std::cout << "Formatted outputs saved to " << validFileName1 << " and " << validFileName2 << std::endl;
183  delete validFile1;
184  delete validFile2;
185 
186  }
187 }
void NDDataMC(const std::string predFileName, const bool makeFile=false, const std::string ldrData="", const std::string ldrMC="")
Definition: NDDataMC.C:32
const XML_Char * name
Definition: expat.h:151
const Var kPartPtp([](const caf::SRProxy *sr){if(std::isnan(1.*sr->sel.nuecosrej.partptp)|| std::isinf(1.*sr->sel.nuecosrej.partptp)) return-5.f;return float(sr->sel.nuecosrej.partptp);})
Definition: NusVars.h:74
const int kNumVars
Definition: vars.h:14
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< HistDef > getAllNusAna01NDHistDefs()
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
dictionary components
const TString name
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2038
const Cut kNus17NDFiducial([](const caf::SRProxy *sr){assert(sr->vtx.elastic.IsValid &&"Must apply EQ cuts");if(sr->vtx.elastic.vtx.X()< -100.0) return false;if(sr->vtx.elastic.vtx.X() > 100.0) return false;if(sr->vtx.elastic.vtx.Y()< -100.0) return false;if(sr->vtx.elastic.vtx.Y() > 100.0) return false;if(sr->vtx.elastic.vtx.Z()< 150.0) return false;if(sr->vtx.elastic.vtx.Z() > 1000.0) return false;return true;})
Definition: NusCuts17.h:103
void SetSpillCut(const SpillCut &cut)
static std::unique_ptr< CheatDecomp > LoadFrom(TDirectory *dir, const std::string &name)
Definition: CheatDecomp.cxx:70
const Cut kNus17NDPresel
Definition: NusCuts17.h:114
const Cut sels[kNumSels]
Definition: vars.h:44
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:607
const Cut kNus17ND
Full Nus17 ND analysis selection.
Definition: NusCuts17.h:120
Struct to hold cut information.
const int kNumSels
Definition: vars.h:43
std::vector< float > Spectrum
Definition: Constants.h:527
const SystShifts kNoShift
Definition: SystShifts.h:115
OStream cout
Definition: OStream.cxx:6
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
const Cut cut
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Var kDummy([](const caf::SRProxy *sr){return 1;})
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:46
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:38
const Cut kNus17EventQuality([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return false;if(sr->sel.nuecosrej.hitsperplane >=8) return false;if(sr->vtx.elastic.fuzzyk.npng==0) return false;if(sr->slc.ncontplanes<=2) return false;return true;})
Data Quality cuts from docdb 21113.
Definition: NusCuts17.h:21
Just return the ND truth spectra as the decomposition.
Definition: CheatDecomp.h:10
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:117
static constexpr Double_t sr
Definition: Munits.h:164