nue_fd_mc_validation.C
Go to the documentation of this file.
1 // makeFile = true saves the spectra for ldrData, ldrMC
2 // makeFile = false reruns formatting for generic validation
3 
4 // SA nue cuts and variables
5 // Same macro should work for different datasets,
6 // as long as all vars and cuts are proper
7 
10 #include "CAFAna/Core/Progress.h"
11 #include "CAFAna/Core/Loaders.h"
12 #include "CAFAna/Core/MultiVar.h"
14 #include "CAFAna/Core/Spectrum.h"
15 #include "CAFAna/Cuts/SpillCuts.h"
16 #include "CAFAna/Cuts/TruthCuts.h"
17 #include "CAFAna/Cuts/NueCutsSecondAna.h"
18 #include "CAFAna/Vars/Vars.h"
20 #include "OscLib/IOscCalc.h"
21 #include "OscLib/OscCalcDumb.h"
23 
24 using namespace ana;
25 
26 #include "TFile.h"
27 
28 #include <iostream>
29 #include <vector>
30 
31 void nue_fd_mc_validation(const std::string predFileName, const bool makeFile=false,
32  const std::string ldrNonSwap="", const std::string ldrFluxSwap="",
33  const std::string ldrTauSwap="")
34 {
35  // Variable groups from NueVarsExtra
36  // Variables from NueVars.h are used here or in the cuts
37 
38  std::vector <HistDef> defs;
39 
40  AddNueHistDefSecondAnaBasic(defs);
41  AddNueHistDefVertexFD(defs);
42  AddNueHistDefShowerFD(defs);
43  AddNueHistDefShower(defs);
44  AddNueHistDefPIDs(defs);
45  AddNueHistDefForLid(defs);
46  AddNueHistDefForLem(defs);
47  AddNueHistDefForCvn(defs);
50  AddNueHistDefTruth(defs);
51  // Variables that shouldn't be reweighted
52  // Set to -1 for no reweights anywhere
53  // const int kIdUnweight = -1;/
54  const int kIdUnweight = defs.size();
55  AddNueHistDefWeight(defs);
56  const Var kkWeight = kPPFXFluxCVWgt * kXSecCVWgt2017;
57 
58  // Cut flow. Change the names/cuts here only if new spectra are going to be saved
59  struct CutDef
60  {
62  Cut cut;
63  };
64 
65  std::vector <CutDef > sels = {
66 // {"0_NoCut", kNoCut},
67  {"1_NueSecondAnaDQ",kNueSecondAnaDQ},
68  {"2_NueFDSecondAnaDecaf",kNueFDSecondAnaDecafCut},
69  {"3_NueSecondAnaFullPresel", kNueSecondAnaFullPresel},
70  {"4_NueSAFDCVNSsb" ,kNueSAFDCVNSsb},
71  {"5_NueSAFDCVNSb", kNueSecondAnaFullPresel && kNueSecondAnaCVNSb},
72  };
73 
74  //Vars for a quick visual comparison of cutflow
75  MultiVar kCutFlowVar (
76  [sels](const caf::SRProxy* sr){
77  std::vector<double> passes;
78  for (unsigned int i=0;i<sels.size();++i){
79  if(sels[i].cut(sr)) passes.push_back(i);
80  }
81  if(sr->mc.nnu==0)passes.push_back(sels.size()+1);
82  return passes;
83  });
84  const Binning flowBins = Binning::Simple(10,0,10);
85  const Var kDummy ( [] (const caf::SRProxy *sr){return 1;});
86  defs.insert(defs.begin(),{"_dummy", {"Dummy", Binning::Simple(3, -0.5,2.5), kDummy}});
87 
88  if(makeFile){
89 
90  SpectrumLoader loaderNonSwap(ldrNonSwap);
91  SpectrumLoader loaderFlxSwap(ldrFluxSwap);
92 
93  NullLoader loaderTauSwap;
94  std::cerr << "WARNING: not trying to load tau swap" << std::endl;
95 
96  loaderNonSwap.SetSpillCut(kStandardSpillCuts);
97  loaderFlxSwap.SetSpillCut(kStandardSpillCuts);
98  loaderTauSwap.SetSpillCut(kStandardSpillCuts);
99 
100  std::vector < std::vector<Spectrum *> > other_nonswap;
101  std::vector < std::vector<Spectrum *> > other_flxswap;
102  std::vector < std::vector<IPrediction*> > preds;
103 
104  for(int selIdx = 0; selIdx<(int)sels.size(); ++selIdx){
105  std::vector<Spectrum *> temp_nonswap;
106  std::vector<Spectrum *> temp_flxswap;
107  std::vector<IPrediction *> temp_pred;
108 
109  for(int varIdx = 0; varIdx < (int)defs.size(); ++varIdx){
110  const HistAxis& axis = defs[varIdx].axis;
111  Var weight = (varIdx < kIdUnweight) ? kkWeight : kUnweighted;
112 
113  auto sp_nonswap = new Spectrum(loaderNonSwap, axis, sels[selIdx].cut && !kHasNeutrino, kNoShift, weight);
114  auto sp_flxswap = new Spectrum(loaderFlxSwap, axis, sels[selIdx].cut && !kHasNeutrino, kNoShift, weight);
115  temp_nonswap.emplace_back(sp_nonswap);
116  temp_flxswap.emplace_back(sp_flxswap);
117 
118  temp_pred.emplace_back
119  (new PredictionNoExtrap(loaderNonSwap,loaderFlxSwap,loaderTauSwap,
120  axis.GetLabels()[0],axis.GetBinnings()[0],axis.GetVars()[0],
121  sels[selIdx].cut, kNoShift, weight));
122  }
123  other_nonswap.push_back(temp_nonswap);
124  other_flxswap.push_back(temp_flxswap);
125  preds.push_back(temp_pred);
126  }
127 
128  auto spfl_nonswap = new Spectrum ("Cut Flow",flowBins,loaderNonSwap, kCutFlowVar, kNoCut);
129  auto spfl_flxswap = new Spectrum ("Cut Flow",flowBins,loaderFlxSwap, kCutFlowVar, kNoCut);
130 
131  loaderNonSwap.Go();
132  loaderFlxSwap.Go();
133 
134  TFile* fout = new TFile(predFileName.c_str(), "RECREATE");
135  Progress prog("Saving...");
136  int nsteps=0;
137  for(int selIdx = 0; selIdx < (int)sels.size(); ++selIdx){
138  const char* cut = sels[selIdx].name.c_str();
139  for(int varIdx = 0; varIdx < (int)defs.size(); ++varIdx){
140  const char* name = defs[varIdx].name.c_str();
141 
142  other_nonswap[selIdx][varIdx]->SaveTo(fout, TString::Format("other_nonswap_%s_%s", name, cut));
143  other_flxswap[selIdx][varIdx]->SaveTo(fout, TString::Format("other_flxswap_%s_%s", name, cut));
144  preds[selIdx][varIdx] ->SaveTo(fout, TString::Format("pred_%s_%s", name, cut));
145  nsteps ++;
146  prog.SetProgress((double)nsteps/(sels.size()*defs.size()+1));
147  }
148  }
149  spfl_nonswap->SaveTo(fout, "spect_nonswap_cutflow");
150  spfl_flxswap->SaveTo(fout, "spect_flxswap_cutflow");
151  prog.Done();
152  delete fout;
153  std::cout << "Saved prediction to " << predFileName << "\n";
154  return;
155  } //end makeFile
156  else{
157 
158  TFile * predFile = new TFile (predFileName.c_str(), "READ");
159  TString validFileName1 (predFileName);
160  validFileName1.ReplaceAll(".root","_mc_validation.root");
161  TFile * validFile1 = new TFile(validFileName1,"RECREATE");
162 
163  int kNumSels = sels.size();
164  int kNumVars = defs.size();
165  // Don't make all the canvases and .json right now
166 // kNumSels = 3;
167  // kNumVars = 5;
168 
169  bool components = 1;
170  TDirectory* comp = 0;
171  if(components) comp = validFile1->mkdir("MCComponents");
172 
173  Progress prog("Arranging spectra...");
174  int nsteps =0;
175 
177  osc::IOscCalc * dumb = new osc::OscCalcDumb();
179  const char* oscName;
180 
181  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
182  const char* cut = sels[selIdx].name.c_str();
183 
184  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
185  const char* name = defs[varIdx].name.c_str();
186 
188  predFile->GetDirectory(TString::Format("pred_%s_%s", name, cut)));
189  auto other_nonswap= Spectrum::LoadFrom(
190  predFile->GetDirectory(TString::Format("other_nonswap_%s_%s", name, cut)));
191  auto other_flxswap = Spectrum::LoadFrom(
192  predFile->GetDirectory(TString::Format("other_flxswap_%s_%s", name, cut)));
193  auto other = *other_nonswap + *other_flxswap;
194 
195  for(int i=0;i<2;i++){
196  if(i==0){
197  calc = noosc;
198  oscName = "NoOscillations";
199  }
200  if(i==1) {
201  calc = dumb;
202  oscName="DumbOscillations";
203  }
204 
205  (pred->Predict(calc)+other).SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Osc,cat=%s}", name,cut,oscName));
206 
207  if (components){
208 
209  pred->PredictComponent(calc, Flavors::kNuMuToNuE,Current::kCC,Sign::kBoth).
210  SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Osc,cat=%s}{group=Component,cat=0_sig}", name,cut,oscName));
211  pred->PredictComponent(calc, Flavors::kNuEToNuE,Current::kCC,Sign::kBoth).
212  SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Osc,cat=%s}{group=Component,cat=1_beam}", name,cut,oscName));
213  pred->PredictComponent(calc, Flavors::kAll,Current::kNC,Sign::kBoth).
214  SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Osc,cat=%s}{group=Component,cat=2_nc}", name,cut,oscName));
215  pred->PredictComponent(calc, Flavors::kAllNuMu,Current::kCC,Sign::kBoth).
216  SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Osc,cat=%s}{group=Component,cat=3_numu}", name,cut,oscName));
217  pred->PredictComponent(calc, Flavors::kAllNuTau,Current::kCC,Sign::kBoth).
218  SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Osc,cat=%s}{group=Component,cat=4_tau}", name,cut,oscName));
219  other.
220  SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Osc,cat=%s}{group=Component,cat=5_other}", name,cut,oscName));
221  }
222  nsteps++;
223  prog.SetProgress((double)nsteps/(kNumSels*kNumVars*2.));
224  }
225  } //end var
226 
227  for(auto oscName:{"NoOscillations","DumbOscillations"}){
228  auto spfl = Spectrum::LoadFrom(predFile, "spect_nonswap_cutflow");
229  spfl->SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Osc,cat=%s}", "_cutflow_nonswap",cut,oscName));
230  spfl = Spectrum::LoadFrom(predFile, "spect_flxswap_cutflow");
231  spfl->SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Osc,cat=%s}", "_cutflow_flxswap",cut,oscName));
232  }
233  }// end cut
234  prog.Done();
235  delete validFile1;
236  std::cout << "Formatted outputs saved to " << validFileName1 << std::endl;
237 
238 
239  }
240 }
const XML_Char * name
Definition: expat.h:151
void AddNueHistDefWeight(std::vector< HistDef > &hd)
const int kNumVars
Definition: vars.h:14
void AddNueHistDefShowerFD(std::vector< HistDef > &hd)
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
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const std::vector< T > & GetVars() const
Definition: HistAxis.h:92
(&#39; appearance&#39;)
Definition: IPrediction.h:18
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
dictionary components
const Var weight
const TString name
const int nsteps
(&#39;beam &#39;)
Definition: IPrediction.h:15
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
OStream cerr
Definition: OStream.cxx:7
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
void SetSpillCut(const SpillCut &cut)
osc::OscCalcDumb calc
_NoOscillations< double > NoOscillations
Definition: IOscCalc.h:67
Charged-current interactions.
Definition: IPrediction.h:39
void AddNueHistDefForLem(std::vector< HistDef > &)
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
Struct to hold cut information.
Simple oscillation probability calculator that has no solar term or mass hierarchy or delta so it&#39;s s...
Definition: OscCalcDumb.h:16
caf::StandardRecord * sr
virtual void Go() override
Load all the registered spectra.
const HistDef defs[kNumVars]
Definition: vars.h:15
const int kNumSels
Definition: vars.h:43
osc::NoOscillations noosc
Definition: NueSystFuncs.h:44
std::vector< float > Spectrum
Definition: Constants.h:610
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
void AddNueHistDefShower(std::vector< HistDef > &hd)
const Cut kHasNeutrino([](const caf::SRProxy *sr){return(sr->mc.nnu!=0);})
Check if MC slice has neutrino information (useful for in-and-out tests)
Definition: TruthCuts.h:61
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
const std::vector< Binning > & GetBinnings() const
Definition: LabelsAndBins.h:69
void AddNueHistDefVertexFD(std::vector< HistDef > &hd)
void AddNueHistDefTruth(std::vector< HistDef > &hd)
void nue_fd_mc_validation(const std::string predFileName, const bool makeFile=false, const std::string ldrNonSwap="", const std::string ldrFluxSwap="", const std::string ldrTauSwap="")
void SetProgress(double frac)
Update the progress fraction between zero and one.
Definition: Progress.cxx:41
const Cut cut
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
Neutral-current interactions.
Definition: IPrediction.h:40
A simple ascii-art progress bar.
Definition: Progress.h:9
void AddNueHistDefForLid(std::vector< HistDef > &hd)
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
const Var kDummy([](const caf::SRProxy *sr){return 1;})
Dummy loader that doesn&#39;t load any files.
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.
All neutrinos, any flavor.
Definition: IPrediction.h:26
const Var kXSecCVWgt2017
Definition: XsecTunes.h:37
static std::unique_ptr< PredictionNoExtrap > LoadFrom(TDirectory *dir, const std::string &name)
void Done()
Call this when action is completed.
Definition: Progress.cxx:90
void AddNueHistDefConfusion(std::vector< HistDef > &hd)
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:96
void AddNueHistDefPIDs(std::vector< HistDef > &hd)
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
const std::vector< std::string > & GetLabels() const
Definition: LabelsAndBins.h:68
enum BeamMode string