FDDataMC.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 
9 #include "CAFAna/Core/Loaders.h"
10 #include "CAFAna/Core/Spectrum.h"
11 #include "CAFAna/Cuts/SpillCuts.h"
12 #include "CAFAna/Cuts/TimingCuts.h"
13 #include "NuXAna/Cuts/NusCuts17.h"
14 #include "CAFAna/Vars/Vars.h"
15 #include "NuXAna/Vars/NusVars.h"
17 #include "OscLib/IOscCalc.h"
19 #include "CAFAna/Analysis/Calcs.h"
22 
23 using namespace ana;
24 
25 #include "TFile.h"
26 
27 #include <iostream>
28 #include <vector>
29 
30 void FDDataMC(const std::string predFileName,
31  const bool makeFile=false,
32  const std::string ldrFDnon="",
33  const std::string ldrFDswp="",
34  const std::string ldrCosmic="")
35 {
36 
37  std::vector<HistDef> defs = getAllNusAna01FDHistDefs();
38 
39  // Cut flow. Change the names/cuts here only if new spectra are going to be saved
40  struct CutDef
41  {
43  Cut cut;
44  };
45 
46  std::vector<CutDef> sels = {
47  {"0_Veto", kApplySecondAnalysisMask && kNus17Veto},
48  {"1_NusEventQuality", kApplySecondAnalysisMask && kNus17Veto && kNus17EventQuality},
49  {"2_NusPresel", kNus17FDPresel},
50  {"3_kNusBackCut", kNus17FDPresel && kNus17BackwardCut},
51  {"4_kNusCosRej", kNus17FDPresel && kNus17BackwardCut && kNus17NCCosRejAlt},
52  {"5_kNusCVN", kNus17FDSansTimePtpAlt},
53  {"6_kNusTime", kNus17FDSansTimePtpAlt && !(kNus17SlcTimeGap && kNus17SlcDist)},
54  {"7_kNusFD", kNus17FDAlt},
55  {"8_kNusFDNhit", kNus17FDExtraAlt}
56  };
57 
58  MultiVar kCutFlowVar (
59  [sels](const caf::SRProxy* sr)
60  {
61  std::vector<double> passes;
62  for(unsigned int i = 0; i < sels.size(); ++i){
63  if(sels[i].cut(sr)) passes.push_back(i);
64  }
65  return passes;
66  });
67  const Binning flowBins = Binning::Simple(6,0,6);
68  const Var kDummy ([] (const caf::SRProxy *sr){return 1;});
69  defs.insert(defs.begin(),{"_dummy", {"Dummy", Binning::Simple(3, -0.5,2.5), kDummy}});
70  if(makeFile){
71  SpectrumLoader loaderMCnon(ldrFDnon.c_str());
72  SpectrumLoader loaderMCswp(ldrFDswp.c_str());
73  SpectrumLoader loaderCosmic(ldrCosmic.c_str(), DataSource::kCosmic);
74  loaderMCnon.SetSpillCut(kStandardSpillCuts);
75  loaderMCswp.SetSpillCut(kStandardSpillCuts);
76  loaderCosmic.SetSpillCut(kStandardSpillCuts);
77 
78  std::vector<std::vector<Spectrum*> > spects;
79  std::vector<std::vector<Spectrum*> > spectsMC;
80  std::vector<std::vector<IPrediction*> > preds;
81 
82  for(int selIdx = 0; selIdx < (int)sels.size(); ++selIdx){
83  std::vector<Spectrum *> temp_spec_non;
84  std::vector<Spectrum *> temp_specMC;
85  std::vector<IPrediction *> temp_pred;
86 
87  for(int varIdx = 0; varIdx < (int)defs.size(); ++varIdx){
88  const HistAxis& axis = defs[varIdx].axis;
89  temp_spec_non.emplace_back(new Spectrum(loaderCosmic,
90  axis,
91  sels[selIdx].cut&&kInTimingSideband,
92  kNoShift));
93 
94  temp_specMC.emplace_back(new Spectrum(loaderMCnon,
95  axis,
96  sels[selIdx].cut && kInBeamSpill,
97  kNoShift,
99 
100  temp_pred.emplace_back(new PredictionNoExtrap(loaderMCnon,
101  loaderMCswp,
102  axis.GetLabels()[0],
103  axis.GetBinnings()[0],
104  axis.GetVars()[0],
105  sels[selIdx].cut&&kInBeamSpill,
106  kNoShift,
108  }
109  spects.push_back(temp_spec_non);
110  preds.push_back(temp_pred);
111  }
112 
113  loaderCosmic.Go();
114  loaderMCnon.Go();
115  loaderMCswp.Go();
116 
117  TFile* fout = new TFile(predFileName.c_str(), "RECREATE");
118  for(int selIdx = 0; selIdx < (int)sels.size(); ++selIdx){
119  const char* cut = sels[selIdx].name.c_str();
120  for(int varIdx = 0; varIdx < (int)defs.size(); ++varIdx){
121  const char* name = defs[varIdx].name.c_str();
122  spects[selIdx][varIdx]->SaveTo(fout, TString::Format("spect_%s_%s", name, cut));
123  preds[selIdx][varIdx]->SaveTo(fout, TString::Format("pred_%s_%s", name, cut));
124 
125  }
126  }
127  delete fout;
128  return;
129  } //end makeFile
130  else{
131  // Play with the format for the generic validation scripts
132  // Here I'm taking the one file and splitting data and MC
133  // Data and MC will be treated as "datasets" = separate files
134  // Alternatively, you could treat them as folders, or categories in a group
135  // Your mileage may vary
136 
137  TFile * predFile = new TFile (predFileName.c_str(), "READ");
138  TString validFileName1 (predFileName); validFileName1.ReplaceAll(".root","_fhc_mc_validation.root");
139  TString validFileName2 (predFileName); validFileName2.ReplaceAll(".root","_cosmic_validation.root");
140 
141  TFile * validFile1 = new TFile(validFileName1,"RECREATE");
142  TFile * validFile2 = new TFile(validFileName2,"RECREATE");
143 
144  int kNumSels = sels.size();
145  int kNumVars = defs.size();
146 
148 
149 
150  for(int selIdx = 0; selIdx < kNumSels; ++selIdx){
151  const char* cut = sels[selIdx].name.c_str();
152  for(int varIdx = 0; varIdx < kNumVars; ++varIdx){
153  const char* name = defs[varIdx].name.c_str();
154 
155  auto spec = Spectrum::LoadFrom(predFile, TString::Format("spect_%s_%s", name, cut));
156  auto pred = PredictionNoExtrap::LoadFrom(predFile, TString::Format("pred_%s_%s", name, cut));
157 
158  spec->SaveTo(validFile2, TString::Format("%s{group=Cut,cat=%s}", name,cut));
159 
160 
161  pred->PredictComponent(calc, ana::Flavors::kAll, ana::Current::kBoth, ana::Sign::kBoth). SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=6_mc_total}", name,cut));
162  pred->PredictComponent(calc, ana::Flavors::kAll, ana::Current::kNC, ana::Sign::kBoth). SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=0_nc}", name,cut));
163  pred->PredictComponent(calc, ana::Flavors::kNuEToNuE, ana::Current::kCC, ana::Sign::kNu). SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=1_nue_surv}", name,cut));
164  pred->PredictComponent(calc, ana::Flavors::kNuEToNuE, ana::Current::kCC, ana::Sign::kAntiNu). SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=2_antinue_surv}", name,cut));
165  pred->PredictComponent(calc, ana::Flavors::kNuMuToNuMu, ana::Current::kCC, ana::Sign::kNu). SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=3_numu_surv}", name,cut));
166  pred->PredictComponent(calc, ana::Flavors::kNuMuToNuMu, ana::Current::kCC, ana::Sign::kAntiNu).SaveTo(validFile1, TString::Format("%s{group=Cut,cat=%s}{group=Component,cat=4_antinumu_surv}", name,cut));
167 
168 
169  }
170  }
171  std::cout << "Formatted outputs saved to "
172  << validFileName1
173  << " and "
174  << validFileName2
175  << std::endl;
176  delete validFile1;
177  delete validFile2;
178 
179  }
180 }
const XML_Char * name
Definition: expat.h:151
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
Antineutrinos-only.
Definition: IPrediction.h:50
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const std::vector< T > & GetVars() const
Definition: HistAxis.h:92
const Cut kApplySecondAnalysisMask([](const caf::SRProxy *sr){if(sr->hdr.det!=caf::kFARDET) return true; std::bitset< 14 > binary(sr->hdr.dibmask);std::pair< int, int > planesA=calcFirstLastLivePlane(sr->slc.firstplane, binary);std::pair< int, int > planesB=calcFirstLastLivePlane(sr->slc.lastplane, binary);if((planesA.first!=planesB.first)||(planesA.second!=planesB.second)) return false;return((planesA.second-planesA.first+1)/64 >=4);})
Definition: AnalysisMasks.h:18
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const TString name
(&#39;beam &#39;)
Definition: IPrediction.h:15
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
const Cut kNus17BackwardCut
FD backward photon rejection cut, inherited from nue, from docdb 21113.
Definition: NusCuts17.h:45
void SetSpillCut(const SpillCut &cut)
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
const Cut kInTimingSideband([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return(kInTimingSideband_before(sr)|| kInTimingSideband_after(sr));else return(kInTimingSideband_before(sr)|| kInTimingSideband_afterA(sr)|| kInTimingSideband_afterB(sr));}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return(kInTimingSideband_before.Livetime(spill)+ kInTimingSideband_after.Livetime(spill));else return(kInTimingSideband_before.Livetime(spill)+ kInTimingSideband_afterA.Livetime(spill)+ kInTimingSideband_afterB.Livetime(spill));}, [](const caf::SRSpillProxy *spill){return 0;})
Definition: TimingCuts.h:12
std::vector< HistDef > getAllNusAna01FDHistDefs()
Charged-current interactions.
Definition: IPrediction.h:39
Interactions of both types.
Definition: IPrediction.h:42
const Cut sels[kNumSels]
Definition: vars.h:44
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
void FDDataMC(const std::string predFileName, const bool makeFile=false, const std::string ldrFDnon="", const std::string ldrFDswp="", const std::string ldrCosmic="")
Definition: FDDataMC.C:30
const Cut kNus17SlcDist
Definition: NusCuts17.h:55
Struct to hold cut information.
caf::StandardRecord * sr
const Cut kNus17NCCosRejAlt
Until nccosrej branch is in CAFs, must use these cuts.
Definition: NusCuts17.h:70
const HistDef defs[kNumVars]
Definition: vars.h:15
const int kNumSels
Definition: vars.h:43
std::vector< float > Spectrum
Definition: Constants.h:610
Neutrinos-only.
Definition: IPrediction.h:49
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
const std::vector< Binning > & GetBinnings() const
Definition: LabelsAndBins.h:69
const Cut kNus17Veto([](const caf::SRProxy *sr){return sr->sel.veto.keep;})
Definition: NusCuts17.h:25
const Cut cut
Collaborates with Spectrum and OscillatableSpectrum to fill spectra from CAF files.
const Cut kNus17SlcTimeGap
Closest slice cuts based on approximation in space and time, from from docdb 21113.
Definition: NusCuts17.h:54
const Cut kNus17FDAlt
Definition: NusCuts17.h:92
Neutral-current interactions.
Definition: IPrediction.h:40
const Cut kInBeamSpill([](const caf::SRProxy *sr){if(sr->spill.run > util::kLastBadTimingRun) return kInBeamSpill_main(sr);else return kInBeamSpill_main(sr)||kInBeamSpill_shifted(sr);}, [](const caf::SRSpillProxy *spill){if(spill->run > util::kLastBadTimingRun) return kInBeamSpill_main.Livetime(spill);else return kInBeamSpill_main.Livetime(spill)+kInBeamSpill_shifted.Livetime(spill);}, [](const caf::SRSpillProxy *spill) -> double{return spill->spillpot;})
Does the event fall inside the window we call the beam spill?
Definition: TimingCuts.h:8
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
const Var kDummy([](const caf::SRProxy *sr){return 1;})
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)
const Cut kNus17FDPresel
The Nus17 preselection cuts for the Far Detector from docdb 21113.
Definition: NusCuts17.h:38
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
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
const std::vector< std::string > & GetLabels() const
Definition: LabelsAndBins.h:68
enum BeamMode string