MakeNus18Systs.C
Go to the documentation of this file.
1 // Make Nus18 systs
2 
4 #include "CAFAna/Core/Loaders.h"
7 #include "CAFAna/Cuts/Cuts.h"
10 #include "CAFAna/Cuts/SpillCuts.h"
17 #include "CAFAna/Systs/BeamSysts.h"
18 #include "CAFAna/Systs/Systs.h"
19 #include "NuXAna/Systs/NusSystsBase.h"
21 #include "CAFAna/Vars/HistAxes.h"
22 #include "NuXAna/Vars/HistAxes.h"
24 
25 
26 using namespace ana;
27 
28 void MakeNus18Systs(TString opt) {
29 
30  TH1::AddDirectory(0);
31 
32  // Get options
33  // TODO add extrap option
34  ECAFType cafType;
35  Loaders::FluxType fluxType;
36  bool extrap = false;
37  if (opt.Contains("caf", TString::kIgnoreCase))
38  cafType = ECAFType::kFullCAF;
39  else if (opt.Contains("concat", TString::kIgnoreCase))
40  cafType = ECAFType::kNusConcat;
41  else {
42  std::cout << "ECAFType required: caf or concat" << std::endl;
43  abort();
44  }
45  if (opt.Contains("fhc", TString::kIgnoreCase))
46  fluxType = Loaders::FluxType::kFHC;
47  else if (opt.Contains("rhc", TString::kIgnoreCase))
48  fluxType = Loaders::FluxType::kRHC;
49  else {
50  std::cout << "FluxType required: fhc or rhc" << std::endl;
51  abort();
52  }
53  if (opt.Contains("extrap", TString::kIgnoreCase))
54  extrap = true;
55 
56  // Set up the loaders
57  Loaders* loaders_nom = new Prod4NomLoaders(cafType, fluxType);
58  Loaders* loaders_up = new Prod4LightLevelLoaders(cafType, fluxType, +1);
59  Loaders* loaders_down = new Prod4LightLevelLoaders(cafType, fluxType, -1);
60  Loaders* loaders_cherenkov = new Prod4CherenkovLoaders(cafType, fluxType);
61  loaders_nom ->SetSpillCut(kStandardSpillCuts);
62  loaders_up ->SetSpillCut(kStandardSpillCuts);
63  loaders_down ->SetSpillCut(kStandardSpillCuts);
64  loaders_cherenkov->SetSpillCut(kStandardSpillCuts);
65 
66  // Set up the new weights
67  const Var kReweight = kXSecCVWgt2018*kPPFXFluxCVWgt;
68  const Var kWeightKaonUp = kReweight * kNus18BeamWeightP;
69  const Var kWeightKaonDown = kReweight * kNus18BeamWeightM;
70 
71  // Set up any syst shifts
72  std::vector<int> sigmas = {-3, -2, -1, 0, 1, 2, 3};
73  std::map<int, SystShifts*> beamShifts;
74  std::map<int, SystShifts*> neutronShifts;
75  for (const auto& sigma : sigmas)
77 
78  // Get axis and cuts
79  HistAxis* nus_axis = new HistAxis(kNus18AxisE);
80  HistAxis* numu_axis = new HistAxis(kNus18BinsNumuCCAxis);
81  Cut* fd_cut = new Cut(kNus18FD);
82  Cut* nd_cut = new Cut(kNus18ND);
83  Cut* numu_cut = new Cut(kNumuCutND2018);
84 
85  // // Make light level syst
86  // NusSystsLoaderShift* lightLevelSyst = new NusSystsLoaderShift("LightLevel");
87  // lightLevelSyst->SetPrediction(nus_axis, numu_axis,
88  // fd_cut, nd_cut, numu_cut,
89  // &kReweight);
90  // lightLevelSyst->AddNominal(loaders_nom);
91  // lightLevelSyst->AddShift("UpLightLevel",
92  // "Light Level Up and Calibration Down",
93  // loaders_up);
94  // lightLevelSyst->AddShift("DownLightLevel",
95  // "Light Level Down and Calibration Up",
96  // loaders_down);
97  // lightLevelSyst->AddShift("Cherenkov",
98  // "Cherenkov Variation",
99  // loaders_cherenkov);
100 
101 // // Make kaon syst
102 // NusSystsWeightShift* kaonSyst = new NusSystsWeightShift("Kaon");
103 // kaonSyst->SetPrediction(nus_axis, numu_axis,
104 // fd_cut, nd_cut, numu_cut,
105 // loaders_nom);
106 // kaonSyst->AddNominal(&kReweight);
107 // kaonSyst->AddShift("KaonUp",
108 // "Kaon Shifted Up",
109 // &kWeightKaonUp);
110 // kaonSyst->AddShift("KaonDown",
111 // "Kaon Shifted Down",
112 // &kWeightKaonDown);
113 
114  // Neutron syst
115  NusSystsWeightShift* neutronSyst = new NusSystsWeightShift("Neutron");
116  neutronSyst->SetPrediction(nus_axis, numu_axis,
117  fd_cut, nd_cut, numu_cut,
118  loaders_nom);
119  neutronSyst->AddNominal(&kReweight);
120  neutronSyst->AddShift("Neutron",
121  "Neutron Primaries",
122  neutronShifts);
123 
124  // Run loaders
125  loaders_nom->Go();
126  loaders_up->Go();
127  loaders_down->Go();
128  loaders_cherenkov->Go();
129 
130  // Run syst makers
131  // lightLevelSyst->Go(6.91e20);
132  kaonSyst->Go(6.91e20);
133  neutronSyst->Go(6.91e20);
134 
135  // Print some shifts
136  // std::vector<std::string> systNames = lightLevelSyst->GetSystNames();
137  // for (const auto& systName : systNames)
138  // std::cout << "Syst " << systName << " has shift +/-" << lightLevelSyst->GetShift(systName, true) << "% (signal), "
139  // << lightLevelSyst->GetShift(systName, false) << "% (background)." << std::endl;
140 // std::vector<std::string> systNames = kaonSyst->GetSystNames();
141 // std::pair<double, double> totalShift = kaonSyst->GetTotalShift();
142 // std::cout << "Kaon: total shift " << totalShift.first*100 << "% (signal), " << totalShift.second*100 << "% (background)" << std::endl;
143 // for (const auto& systName : systNames)
144 // std::cout << "Syst " << systName << " has shift +/-" << kaonSyst->GetShift(systName, true)*100 << "% (signal), "
145 // << kaonSyst->GetShift(systName, false)*100 << "% (background)." << std::endl;
146  std::vector<std::string> systNames = neutronSyst->GetSystNames();
147  std::pair<double, double> totalShift = neutronSyst->GetTotalShift();
148  std::cout << "Neutron: total shift " << totalShift.first*100 << "% (signal), " << totalShift.second*100 << "% (background)" << std::endl;
149  for (const auto& systName : systNames)
150  std::cout << "Syst " << systName << " has shift +/-" << neutronSyst->GetShift(systName, true)*100 << "% (signal), "
151  << neutronSyst->GetShift(systName, false)*100 << "% (background)." << std::endl;
152 
153  return;
154 
155 }
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var kPPFXFluxCVWgt
Definition: PPFXWeights.h:16
const Var kNus18BeamWeightM
Definition: NusVars.h:90
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const Cut kNumuCutND2018
Definition: NumuCuts2018.h:41
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:23
Loaders for Cherenkov paths/definitions.
Definition: Prod4Loaders.h:172
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
const Cut kNus18FD
Definition: NusCuts18.h:100
const Var kNus18BeamWeightP
Definition: NusVars.h:87
const HistAxis kNus18AxisE("Energy Deposited in Scintillator (GeV)", kNus18EnergyBinning, kNus18Energy)
Axes used in Nus18 analysis by nus group.
Definition: HistAxes.h:16
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
Loaders for light level paths/definitions.
Definition: Prod4Loaders.h:145
void MakeNus18Systs(TString opt)
_Cut< caf::SRProxy > Cut
Representation of a cut (selection) to be applied to a caf::StandardRecord object.
Definition: Cut.h:96
const Cut kNus18ND
Full Nus18 ND analysis selection.
Definition: NusCuts18.h:137
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
std::map< std::string, std::string > systNames
Definition: PlotUnfolding.C:32
const Var kXSecCVWgt2018
Definition: XsecTunes.h:49
const HistAxis kNus18BinsNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNus18EnergyBinning, kCCE)
Definition: HistAxes.h:17
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
ECAFType
Definition: Loaders.h:19