MakeSystRatios.C
Go to the documentation of this file.
1 
2 #include "TString.h"
3 #include "TFile.h"
4 
6 
7 #include "CAFAna/Core/ISyst.h"
8 #include "CAFAna/Core/Spectrum.h"
10 #include "NuXAna/Cuts/NusCuts18.h"
11 #include "CAFAna/Vars/HistAxes.h"
12 #include "NuXAna/Vars/HistAxes.h"
13 #include "CAFAna/Vars/XsecTunes.h"
17 
18 using namespace ana;
19 
20 void MakeSystRatios(TString opt) {
21 
22  bool nd = opt.Contains("neardet");
23 
24  // Here get loaders
26  if (nd) loaders.SetND(true);
28 
29  // This should surely be a vector of loaders and cuts, right?
30 
31  auto det = nd ? caf::kNEARDET : caf::kFARDET;
32 
33  std::vector<std::tuple<SpectrumLoaderBase*, const Cut, std::string>> components;
34 
35  SpectrumLoaderBase* nonswap = &loaders.GetLoader(det,
37 
38  // Unoscillated components (for both detectors)
39 
40  components.push_back(std::make_tuple(nonswap, kIsBeamNue && !kIsAntiNu, "cc_nue_to_nue" ));
41  components.push_back(std::make_tuple(nonswap, kIsBeamNue && kIsAntiNu, "cc_nuebar_to_nuebar" ));
42  components.push_back(std::make_tuple(nonswap, kIsNumuCC && !kIsAntiNu, "cc_numu_to_numu" ));
43  components.push_back(std::make_tuple(nonswap, kIsNumuCC && kIsAntiNu, "cc_numubar_to_numubar"));
44  components.push_back(std::make_tuple(nonswap, kIsNC, "nc" ));
45 
46  // CC appearance channels (for far detector only)
47 
48  if (!nd) {
51  SpectrumLoaderBase* tau = &loaders.GetLoader(det,
53  components.push_back(std::make_tuple(swap, kIsSig && !kIsAntiNu, "cc_numu_to_nue" ));
54  components.push_back(std::make_tuple(swap, kIsSig && kIsAntiNu, "cc_numubar_to_nuebar" ));
55  components.push_back(std::make_tuple(swap, kIsNumuApp && !kIsAntiNu, "cc_nue_to_numu" ));
56  components.push_back(std::make_tuple(swap, kIsNumuApp && kIsAntiNu, "cc_nuebar_to_numubar" ));
57  components.push_back(std::make_tuple(tau, kIsTauFromE && !kIsAntiNu, "cc_nue_to_nutau" ));
58  components.push_back(std::make_tuple(tau, kIsTauFromE && kIsAntiNu, "cc_nuebar_to_nutaubar" ));
59  components.push_back(std::make_tuple(tau, kIsTauFromMu && !kIsAntiNu, "cc_numu_to_nutau" ));
60  components.push_back(std::make_tuple(tau, kIsTauFromMu && kIsAntiNu, "cc_numubar_to_nutaubar"));
61  }
62 
63  // PPFX weight vars
64 
65  std::vector<Var> kPPFXFluxUnivWgts;
66  for (size_t i = 0; i < 100; ++i) {
67  const Var tempPPFXWgt([i](const caf::SRProxy* sr){
68  if (!sr->hdr.ismc) return 1.f;
69  if (sr->mc.nnu != 1) return 1.f;
70  if (sr->mc.nu[0].rwgt.ppfx.vuniv[i] <= 0) return 1.f;
71  return (float)sr->mc.nu[0].rwgt.ppfx.vuniv[i];
72  });
73  kPPFXFluxUnivWgts.push_back(tempPPFXWgt);
74  }
75 
76  // GENIE systs
77 
78  std::vector<const ISyst*> systs = getAllXsecSysts_2018();
79 
80  // Analysis cuts and vars
81 
82  const HistAxis axis = nd ? kNCNDAxisE : kNCFDAxisE;
83  const Cut sel = nd ? kNus18ND : kNus18FD;
84  const Var kReweight = kXSecCVWgt2018 * kPPFXFluxCVWgt;
85 
86  // Initialise maps
87 
88  std::vector<Spectrum*> nominal(components.size());
89  std::vector<std::vector<Spectrum*> > ppfx(components.size(),
90  std::vector<Spectrum*>(kPPFXFluxUnivWgts.size(), nullptr));
91  std::vector<std::vector<Spectrum*> > genie_up(components.size(),
92  std::vector<Spectrum*>(systs.size(), nullptr));
93  std::vector<std::vector<Spectrum*> > genie_down(components.size(),
94  std::vector<Spectrum*>(systs.size(), nullptr));
95 
96  // Loop over each component
97 
98  for (size_t c = 0; c < components.size(); ++c) {
99 
100  // Nominal spectra
101 
102  nominal[c] = new Spectrum(*std::get<0>(components[c]),
103  axis, sel && std::get<1>(components[c]), kNoShift, kReweight);
104 
105  for (size_t u = 0; u < 100; ++u)
106  ppfx[c][u] = new Spectrum(*std::get<0>(components[c]),
107  axis, sel && std::get<1>(components[c]), kNoShift, kXSecCVWgt2018 * kPPFXFluxUnivWgts[u]);
108 
109  for (size_t s = 0; s < systs.size(); ++s) {
110 
111  SystShifts shift_up;
112  shift_up.SetShift(systs[s], 1);
113  genie_up[c][s] = new Spectrum(*std::get<0>(components[c]),
114  axis, sel && std::get<1>(components[c]), shift_up, kReweight);
115  SystShifts shift_down;
116  shift_down.SetShift(systs[s], -1);
117  genie_down[c][s] = new Spectrum(*std::get<0>(components[c]),
118  axis, sel && std::get<1>(components[c]), shift_down, kReweight);
119 
120  } // for genie syst
121 
122  } // for component
123 
124  loaders.Go();
125 
126  // Save em all
127 
128  TFile* f = TFile::Open("SystRatios.root", "recreate");
129 
130  TDirectory* dir = f->mkdir("nominal");
131  for (size_t c = 0; c < components.size(); ++c)
132  nominal[c]->SaveTo(dir, std::get<2>(components[c]).c_str());
133 
134  dir = f->mkdir("ppfx");
135  for (size_t u = 0; u < 100; ++u) {
136  TDirectory* udir = dir->mkdir(Form("u%zu",u+1));
137  for (size_t c = 0; c < components.size(); ++c)
138  ppfx[c][u]->SaveTo(udir, std::get<2>(components[c]).c_str());
139  }
140 
141  dir = f->mkdir("genie");
142  for (size_t s = 0; s < systs.size(); ++s) {
143  TDirectory* sdir = dir->mkdir(systs[s]->ShortName().c_str());
144  TDirectory* updir = sdir->mkdir("up");
145  TDirectory* downdir = sdir->mkdir("down");
146  for (size_t c = 0; c < components.size(); ++c) {
147  genie_up[c][s]->SaveTo(updir, std::get<2>(components[c]).c_str());
148  genie_down[c][s]->SaveTo(downdir, std::get<2>(components[c]).c_str());
149  }
150  }
151 
152  // loop and save things here
153 
154  delete f;
155 
156 }
Near Detector underground.
Definition: SREnums.h:10
const Cut kIsNumuCC(CCFlavSel(14, 14))
Select CC .
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< const ISyst * > getAllXsecSysts_2018()
Get master XSec syst list for 2018 analyses.
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
dictionary components
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2137
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
const Cut kIsAntiNu([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return sr->mc.nu[0].pdg< 0;})
Is this truly an antineutrino?
Definition: TruthCuts.h:53
const Cut kIsBeamNue(CCFlavSel(12, 12))
Select CC .
const Cut kIsNumuApp(CCFlavSel(14, 12))
Select CC .
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
const Cut kNus18FD
Definition: NusCuts18.h:100
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
const XML_Char * s
Definition: expat.h:262
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
const Cut kIsTauFromE(CCFlavSel(16, 12))
Select CC .
const HistAxis kNCNDAxisE("Energy Deposited in Scintillator (GeV)", kNCNDBinning, kNus18Energy)
NC covariance hist axes.
Definition: HistAxes.h:20
const Cut kIsTauFromMu(CCFlavSel(16, 14))
Select CC .
if(dump)
const Cut kIsNC([](const caf::SRProxy *sr){if(sr->mc.nnu==0) return false;assert(sr->mc.nnu==1);return!sr->mc.nu[0].iscc;})
Is this a Neutral Current event?
Definition: TruthCuts.h:8
SpectrumLoaderBase & GetLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Retrieve a specific loader.
Definition: Loaders.cxx:129
void MakeSystRatios(TString opt)
const Cut kIsSig(CCFlavSel(12, 14))
Select CC .
caf::StandardRecord * sr
const Cut kNus18ND
Full Nus18 ND analysis selection.
Definition: NusCuts18.h:137
std::vector< float > Spectrum
Definition: Constants.h:610
const SystShifts kNoShift
Definition: SystShifts.cxx:21
Base class for the various types of spectrum loader.
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
const HistAxis kNCFDAxisE("Energy Deposited in Scintillator (GeV)", kNCFDBinning, kNus18Energy)
Definition: HistAxes.h:21
TDirectory * dir
Definition: macro.C:5
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
std::vector< Loaders * > loaders
Definition: syst_header.h:386
void SetND(bool nd)
Definition: Loaders.h:65
caf::Proxy< bool > ismc
Definition: SRProxy.h:242
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod4Loaders.h:96
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80