Functions
nue_data_mc_validation.C File Reference
#include "CAFAna/Decomp/CheatDecomp.h"
#include "CAFAna/Core/Progress.h"
#include "CAFAna/Core/MultiVar.h"
#include "CAFAna/Core/SpectrumLoader.h"
#include "CAFAna/Core/Spectrum.h"
#include "CAFAna/Cuts/SpillCuts.h"
#include "CAFAna/Cuts/NueCutsSecondAna.h"
#include "CAFAna/Vars/Vars.h"
#include "3FlavorAna/Vars/NueVarsExtra.h"
#include "OscLib/IOscCalc.h"
#include "CAFAna/Vars/GenieWeights.h"
#include "TFile.h"
#include <iostream>
#include <vector>

Go to the source code of this file.

Functions

void nue_data_mc_validation (const std::string predFileName, const bool makeFile=false, const std::string ldrData="", const std::string ldrMC="")
 

Function Documentation

void nue_data_mc_validation ( const std::string  predFileName,
const bool  makeFile = false,
const std::string  ldrData = "",
const std::string  ldrMC = "" 
)

Definition at line 22 of file nue_data_mc_validation.C.

References ana::AddNueHistDefConfusion(), ana::AddNueHistDefForLem(), ana::AddNueHistDefForLid(), ana::AddNueHistDefPIDs(), ana::AddNueHistDefSelectionExtras(), ana::AddNueHistDefShower(), ana::AddNueHistDefShowerND(), ana::AddNueHistDefTruth(), ana::AddNueHistDefVertexND(), ana::AddNueHistDefWeight(), allInOneTrainingPlots::axis, om::cout, ana::CutDef::cut, ana::defs, ana::Progress::Done(), allTimeWatchdog::endl, genie::utils::style::Format(), submit_syst::fout, MECModelEnuComparisons::i, makeTrainCVSamples::int, ana::kDummy, ana::kNoCut, ana::kNoShift, PandAna.cut.analysis_cuts::kNueNDCVNSsb, ana::kNumSels, ana::kNumVars, ana::kPPFXFluxCVWgt, ana::kStandardDQCuts, ana::kStandardSpillCuts, ana::kUnweighted, ana::kXSecCVWgt2017, ana::CheatDecomp::LoadFrom(), ana::Spectrum::LoadFrom(), ana::CutDef::name, nsteps, plot_validation_datamc::pred, add_attributes::prog, ana::SaveTo(), ana::sels, ana::Progress::SetProgress(), ana::SpectrumLoaderBase::SetSpillCut(), ana::Binning::Simple(), sr, string, and ana::weight.

24 {
25  // Variable groups from NueVarsExtra
26  // Variables from NueVars.h are used here or in the cuts
27 
28  std::vector <HistDef> defs;
29  AddNueHistDefSecondAnaBasic(defs);
30  AddNueHistDefVertexND(defs);
31  AddNueHistDefShowerND(defs);
32  AddNueHistDefShower(defs);
33  AddNueHistDefPIDs(defs);
34  AddNueHistDefForLid(defs);
35  AddNueHistDefForLem(defs);
36  AddNueHistDefForCvn(defs);
39  AddNueHistDefTruth(defs);
40 
41  //Don't apply weights to next / all variables
42  //const int kIdUnweight = -1;//defs.size();
43  const int kIdUnweight = defs.size();
44  AddNueHistDefWeight(defs);
45  const Var kkWeight = kPPFXFluxCVWgt * kXSecCVWgt2017;
46 
47  // Cut flow. Change the names/cuts here only if new spectra are going to be saved
48  struct CutDef
49  {
51  Cut cut;
52  };
53 
54  std::vector <CutDef > sels = {
55  {"0_NoCut", kNoCut},
56  {"1_NueSecondAnaDQ",kNueSecondAnaDQ},
57  {"2_NueNDSecondAnaDecaf",kNueNDSecondAnaDecafCut},
58  {"3_NueNDSecondAnaPresel", kNueNDSecondAnaPresel},
59  {"4_NueNDCVNSsb" ,kNueNDCVNSsb},
60  {"5_NueNDCVNSb",kNueNDSecondAnaPresel && kNueSecondAnaCVNSb},
61  };
62 
63  //Vars for a quick visual comparison of cutflow
64  MultiVar kCutFlowVar ([sels](const caf::SRProxy* sr){
65  std::vector<double> passes;
66  for (unsigned int i=0;i<sels.size();++i){
67  if(sels[i].cut(sr)) passes.push_back(i);
68  }
69  return passes;
70  });
71  const Binning flowBins = Binning::Simple(10,0,10);
72  const Var kDummy ( [] (const caf::SRProxy *sr){return 1;});
73  defs.insert(defs.begin(),{"_dummy", {"Dummy", Binning::Simple(3, -0.5,2.5), kDummy}});
74 
75  //Save Data and MC spectra
76 
77  if(makeFile){
78  SpectrumLoader loaderData(ldrData.c_str());
79  SpectrumLoader loaderMC(ldrMC.c_str());
80 
81  loaderData.SetSpillCut(kStandardSpillCuts);
82  loaderMC.SetSpillCut(kStandardDQCuts);
83 
84  std::vector < std::vector<Spectrum*> > spects;
85  std::vector < std::vector<Spectrum*> > specMC; //useful when truth info is missing
86  std::vector < std::vector<IDecomp*> > preds;
87 
88  for(int selIdx = 0; selIdx<(int)sels.size(); ++selIdx){
89  std::vector<Spectrum *> temp_spec;
90  std::vector<Spectrum *> temp_spMC;
91  std::vector<IDecomp *> temp_pred;
92  for(int varIdx = 0; varIdx < (int)defs.size(); ++varIdx){
93  const HistAxis& axis = defs[varIdx].axis;
94  Var weight = (varIdx < kIdUnweight) ? kkWeight:kUnweighted;
95 
96  temp_spec.emplace_back(new Spectrum (loaderData, axis, sels[selIdx].cut, kNoShift));
97  temp_spMC.emplace_back(new Spectrum (loaderMC, axis, sels[selIdx].cut, kNoShift, weight));
98  temp_pred.emplace_back(new CheatDecomp(loaderMC, axis, sels[selIdx].cut, kNoShift, weight));
99  }
100  spects.push_back(temp_spec);
101  specMC.push_back(temp_spMC);
102  preds.push_back(temp_pred);
103  }
104 
105  auto spfl_data = new Spectrum ("Cut Flow",flowBins,loaderData, kCutFlowVar, kNoCut);
106  auto spfl_ndmc = new Spectrum ("Cut Flow",flowBins,loaderMC, kCutFlowVar, kNoCut);
107 
108  loaderMC.Go();
109  loaderData.Go();
110 
111  TFile* fout = new TFile(predFileName.c_str(), "RECREATE");
112  Progress prog("Saving...");
113  int nsteps=0;
114  for(int selIdx = 0; selIdx < (int)sels.size(); ++selIdx){
115  const char* cut = sels[selIdx].name.c_str();
116  for(int varIdx = 0; varIdx < (int)defs.size(); ++varIdx){
117  const char* name = defs[varIdx].name.c_str();
118  spects[selIdx][varIdx]->SaveTo(fout, TString::Format("spect_%s_%s", name, cut));
119  specMC[selIdx][varIdx]->SaveTo(fout, TString::Format("speMC_%s_%s", name, cut));
120  preds[selIdx][varIdx]->SaveTo(fout, TString::Format("pred_%s_%s", name, cut));
121  nsteps ++;
122  prog.SetProgress((double)nsteps/(sels.size()*defs.size()+1));
123  }
124  }
125  spfl_data->SaveTo(fout, "spect_data_cutflow");
126  spfl_ndmc->SaveTo(fout, "spect_mc_cutflow");
127  prog.Done();
128 
129  delete fout;
130  std::cout << "Saved prediction to " << predFileName << "\n";
131  return;
132  } //end makeFile
133  else{
134  // Play with the format for the generic validation scripts
135  // Here I'm taking the one file and splitting data and MC
136  // Data and MC will be treated as "datasets" = separate files
137  // Alternatively, you could treat them as folders, or categories in a group
138  // Your mileage may vary
139 
140  TFile * predFile = new TFile (predFileName.c_str(), "READ");
141  TString validFileName1 (predFileName); validFileName1.ReplaceAll(".root","_mc_validation.root");
142  TString validFileName2 (predFileName); validFileName2.ReplaceAll(".root","_data_validation.root");
143 
144  TFile * validFile1 = new TFile(validFileName1,"RECREATE");
145  TFile * validFile2 = new TFile(validFileName2,"RECREATE");
146 
147  int kNumSels = sels.size();
148  int kNumVars = defs.size();
149 
150  // Don't make all the canvases right now
151 // kNumSels = 2;
152 // kNumVars = 2;
153  bool components = 0;
154 
155  Progress prog ("Arranging spectra...");
156  int nsteps = 0;
157 
158  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
159  const char* cut = sels[selIdx].name.c_str();
160 
161  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
162  const char* name = defs[varIdx].name.c_str();
163  auto spec = Spectrum::LoadFrom(predFile, TString::Format("spect_%s_%s", name, cut).Data());
164  auto spMC = Spectrum::LoadFrom(predFile, TString::Format("speMC_%s_%s", name, cut).Data());
165  auto pred = CheatDecomp::LoadFrom (predFile, TString::Format("pred_%s_%s", name, cut).Data());
166 
167  // Data and MC for variable "name" will be plotted together
168  // Plots will have selector for cut level {group=Cut,cat=%s}
169  // Optional MCComponents folder will have extra selector; total MC = all_sum will be created as well
170 
171  spec->SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}", name,cut));
172  spMC->SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}", name,cut));
173 
174  if (components){
175  TDirectory* comp = validFile1->GetDirectory("MCComponents");
176  if(!comp) comp = validFile1->mkdir("MCComponents");
177  pred->NueComponent(). SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Component,cat=0_nue}", name,cut));
178  pred->NumuComponent(). SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Component,cat=1_numu}", name,cut));
179  pred->AntiNueComponent(). SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Component,cat=3_antinue}", name,cut));
180  pred->AntiNumuComponent().SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Component,cat=4_antinumu}", name,cut));
181  pred->NCTotalComponent(). SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Component,cat=2_nc}", name,cut));
182  (*spMC-(pred->TotalMC())). SaveTo(comp, TString::Format("%s_{group=Cut,cat=%s}{group=Component,cat=5_other}", name,cut));
183  }
184  prog.SetProgress((double)nsteps/(kNumSels*kNumVars));
185  nsteps++;
186  } //end var
187  auto spfl = Spectrum::LoadFrom(predFile, "spect_data_cutflow");
188  spfl->SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}", "_cutflow",cut));
189 
190  spfl = Spectrum::LoadFrom(predFile, "spect_mc_cutflow");
191  spfl->SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}", "_cutflow",cut));
192 
193  }// end cut
194  prog.Done();
195  std::cout << "Closing files...\n";
196  delete validFile1;
197  delete validFile2;
198  std::cout << "Formatted outputs saved to " << validFileName1 << " and " << validFileName2 << std::endl;
199 
200  }
201 }
const XML_Char * name
Definition: expat.h:151
void AddNueHistDefWeight(std::vector< HistDef > &hd)
const int kNumVars
Definition: vars.h:14
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var weight
const int nsteps
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
std::vector< double > Spectrum
Definition: Constants.h:743
void SetSpillCut(const SpillCut &cut)
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void AddNueHistDefShowerND(std::vector< ana::HistDef > &hd)
void AddNueHistDefForLem(std::vector< HistDef > &)
void AddNueHistDefSelectionExtras(std::vector< HistDef > &hd)
const Cut sels[kNumSels]
Definition: vars.h:44
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.
caf::StandardRecord * sr
const Var kDummy([](const caf::SRProxy *sr){return 1;})
const HistDef defs[kNumVars]
Definition: vars.h:15
const int kNumSels
Definition: vars.h:43
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
void AddNueHistDefShower(std::vector< HistDef > &hd)
void AddNueHistDefTruth(std::vector< HistDef > &hd)
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
A simple ascii-art progress bar.
Definition: Progress.h:9
void AddNueHistDefForLid(std::vector< HistDef > &hd)
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
void AddNueHistDefVertexND(std::vector< HistDef > &hd)
const Var kXSecCVWgt2017
Definition: XsecTunes.h:36
std::string name
void AddNueHistDefConfusion(std::vector< HistDef > &hd)
Just return the ND truth spectra as the decomposition.
Definition: CheatDecomp.h:10
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:96
void AddNueHistDefPIDs(std::vector< HistDef > &hd)
enum BeamMode string