MakeNus18Systs.C
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////
2 // MakeNus18Systs.C
3 // Mike Wallbank (wallbank@fnal.gov)
4 //
5 // Make all systematics for Nus18 analyses.
6 //////////////////////////////////////////////////////////////
7 
9 #include "NuXAna/Analysis/NusSystsMaker.h"
11 #include "CAFAna/Core/Loaders.h"
12 #include "CAFAna/Core/SystShifts.h"
13 #include "CAFAna/Core/Utilities.h"
14 #include "CAFAna/Core/Var.h"
15 #include "CAFAna/Cuts/Cuts.h"
16 #include "NuXAna/Cuts/NusCuts18.h"
18 #include "CAFAna/Cuts/SpillCuts.h"
25 #include "CAFAna/Systs/BeamSysts.h"
26 #include "CAFAna/Systs/Systs.h"
28 #include "CAFAna/Systs/XSecSysts.h"
30 #include "CAFAna/Vars/HistAxes.h"
31 #include "NuXAna/Vars/HistAxes.h"
33 
35 
36 using namespace ana;
37 
38 // -------------------------------------------------------------------------------------
39 void MakeNus18Systs(TString opt) {
40 
41  TH1::AddDirectory(0);
42 
43  // Get options
44  // TODO add extrap option
45  bool fhc = opt.Contains("FHC", TString::kIgnoreCase);
46  bool rhc = opt.Contains("RHC", TString::kIgnoreCase);
47  assert(fhc xor rhc);
48 
49  bool caf = opt.Contains("caf", TString::kIgnoreCase);
50  bool concat = opt.Contains("concat", TString::kIgnoreCase);
51  assert(caf xor concat);
52 
53  bool extrap = opt.Contains("extrap", TString::kIgnoreCase);
54  bool fd = opt.Contains("fd", TString::kIgnoreCase);
55  bool nd = opt.Contains("nd", TString::kIgnoreCase);
56  assert(extrap xor fd xor nd);
57 
58  // Set up the loaders
61  Loaders* loaders_nom = new Prod4NomLoaders(cafType, fluxType);
62  Loaders* loaders_light_nom = new Prod4LightLevelLoaders(cafType, fluxType, 0);
63  Loaders* loaders_light_up = new Prod4LightLevelLoaders(cafType, fluxType, +1);
64  Loaders* loaders_light_down = new Prod4LightLevelLoaders(cafType, fluxType, -1);
65  Loaders* loaders_cherenkov = new Prod4CherenkovLoaders(cafType, fluxType);
66  Loaders* loaders_calib_up = new Prod4AbsCalibLoaders(cafType, fluxType, +1);
67  Loaders* loaders_calib_down = new Prod4AbsCalibLoaders(cafType, fluxType, -1);
68  Loaders* loaders_calib_shape = new Prod4CalibShapeLoaders(cafType, fluxType);
69  loaders_nom ->SetSpillCut(kStandardSpillCuts);
70  loaders_light_nom ->SetSpillCut(kStandardSpillCuts);
71  loaders_light_up ->SetSpillCut(kStandardSpillCuts);
72  loaders_light_down ->SetSpillCut(kStandardSpillCuts);
73  loaders_cherenkov ->SetSpillCut(kStandardSpillCuts);
74  loaders_calib_up ->SetSpillCut(kStandardSpillCuts);
75  loaders_calib_down ->SetSpillCut(kStandardSpillCuts);
76  loaders_calib_shape->SetSpillCut(kStandardSpillCuts);
77 
78  // Set up the weights
79  const Var kReweight = kXSecCVWgt2018 * kPPFXFluxCVWgt;
80 
81  // Sigma variations for any syst shifts
82  std::vector<int> sigmas = {-3, -2, -1, 0, 1, 2, 3};
83 
84  // Get axis and cuts
85  HistAxis* nus_axis = new HistAxis(kNus18AxisE);
86  HistAxis* numu_axis = new HistAxis(kNus18BinsNumuCCAxis);
87  Cut* fd_cut = new Cut(kNus18FD);
88  Cut* nd_cut = new Cut(kNus18ND);
89  Cut* numu_cut = new Cut(kNumuCutND2018);
90 
91  // Get prediction type
92  // NB pred_type is currently passed to NusSystsMaker, so each NusSystsMaker
93  // can only handle one prediction type.
94  // It could instead be passed to NusSystMaker::AddSystematic, to be used when
95  // making NusSystsShiftMaker (this is actually all that happens under the hood
96  // right now), if we wanted the maker to be able to handle systematics from
97  // multiple sources.
98  NusSystsPredType pred_type = NusSystsPredType::kUnknown;
99  if (extrap) pred_type = NusSystsPredType::kExtrap;
100  if (fd) pred_type = NusSystsPredType::kFD;
101  if (nd) pred_type = NusSystsPredType::kND;
102 
103  // Set up systematic maker
104  NusSystematicsMaker* nus_systs_maker = new NusSystematicsMaker("Nus18", pred_type);
105  nus_systs_maker->SetAxes(nus_axis, numu_axis);
106  nus_systs_maker->SetCuts(fd_cut, nd_cut, numu_cut);
107  nus_systs_maker->SetNominalLoaders(loaders_nom, loaders_light_nom);
108  nus_systs_maker->SetNominalWeight(&kReweight);
109 
110  // Beam transport systematic
111  std::map<std::string, std::vector<std::string> > beamShiftLabels = GetBeamTranspShiftLabels();
112  std::string beamShiftFile = FindCAFAnaDir() + "/data/beam/TABeamSyst.root";
113  NusSystMaker* beamTransSyst = new NusSystMaker("BeamTrans", "Beam Transport");
114  for (const auto& shiftLabel : beamShiftLabels) {
115  NusSystSystShift* beamTransShift = new NusSystSystShift(shiftLabel.first, shiftLabel.second[1]);
116  for (const auto& sigma : sigmas)
117  beamTransShift->AddSigma(new SystShifts(new BeamSyst(beamShiftFile, shiftLabel.second[0], shiftLabel.second[1]), sigma),
118  sigma);
119  beamTransSyst->AddShift(beamTransShift);
120  }
121  nus_systs_maker->AddSystematic(beamTransSyst);
122 
123  // PPFX systematic
124  std::vector<Var> ppfxFluxUnivWgts = GetPPFXFluxWeights();
125  Var* universeWgt;
126  NusSystMaker* ppfxSyst = new NusSystMaker("PPFX", "PPFX");
127  int universe = 0;
128  for (const auto& ppfxFluxWgt : ppfxFluxUnivWgts) {
129  NusSystWeightShift* ppfxShift = new NusSystWeightShift(Form("Universe%dShift", universe),
130  Form("Universe%dShift", universe));
131  universeWgt = new Var(ppfxFluxWgt * kXSecCVWgt2018);
132  ppfxShift->AddOnOff(universeWgt);
133  ppfxSyst->AddShift(ppfxShift);
134  ++universe;
135  }
136  nus_systs_maker->AddSystematic(ppfxSyst);
137 
138  // GENIE systematics -- general
139  NusSystMaker* genieGeneralSyst = new NusSystMaker("GenieGeneral", "GENIE General");
140  std::map<std::string, const ISyst*> xsecLabels = GetXSecShiftLabels();
141  for (const auto& xsecLabel : xsecLabels) {
142  NusSystSystShift* genieGeneralShift = new NusSystSystShift(xsecLabel.first, xsecLabel.first);
143  for (const auto& sigma : sigmas)
144  genieGeneralShift->AddSigma(new SystShifts(xsecLabel.second, sigma),
145  sigma);
146  genieGeneralSyst->AddShift(genieGeneralShift);
147  }
148  nus_systs_maker->AddSystematic(genieGeneralSyst);
149 
150  // GENIE systematics -- nova
151  NusSystMaker* genieNovaSyst = new NusSystMaker("GenieNova", "GENIE NOvA");
152  std::map<std::string, rwgt::ReweightLabel_t> genieLabels = GetGENIEShiftLabels();
153  for (const auto& genieLabel : genieLabels) {
154  NusSystSystShift* genieNovaShift = new NusSystSystShift(genieLabel.first, genieLabel.first);
155  for (const auto& sigma : sigmas)
156  genieNovaShift->AddSigma(new SystShifts(GetGenieKnobSyst(genieLabel.second), sigma),
157  sigma);
158  genieNovaSyst->AddShift(genieNovaShift);
159  }
160  nus_systs_maker->AddSystematic(genieNovaSyst);
161 
162  // Absolute calibration systs
163  NusSystMaker* absCalibFlatSyst = new NusSystMaker("AbsCalibFlat",
164  "Absolute Calibration Flat");
165  NusSystLoaderShift* absCalibFlat = new NusSystLoaderShift("AbsCalibFlat",
166  "Absolute Calibration Flat");
167  absCalibFlat->AddSigma(loaders_calib_up,
168  +1);
169  absCalibFlat->AddSigma(loaders_calib_down,
170  -1);
171  absCalibFlatSyst->AddShift(absCalibFlat);
172  nus_systs_maker->AddSystematic(absCalibFlatSyst);
173 
174  NusSystMaker* absCalibShapeSyst = new NusSystMaker("AbsCalibShape",
175  "Absolute Calibration Shape");
176  NusSystLoaderShift* absCalibShape = new NusSystLoaderShift("AbsCalibShape",
177  "Absolute Calibration Shape");
178  absCalibShape->AddOnOff(loaders_calib_shape);
179  absCalibShapeSyst->AddShift(absCalibShape);
180  nus_systs_maker->AddSystematic(absCalibShapeSyst);
181 
182  // Relative calibration systs
183  NusSystMaker* relCalibNDFlatSyst = new NusSystMaker("RelCalibNDFlatUp",
184  "Relative Calibration ND Flat");
185  NusSystLoaderShift* relCalibNDFlat = new NusSystLoaderShift("RelCalibNDFlatUp",
186  "Relative Calibration ND Flat");
187  relCalibNDFlat->AddSigma(loaders_calib_up,
188  loaders_nom,
189  +1);
190  relCalibNDFlat->AddSigma(loaders_calib_down,
191  loaders_nom,
192  -1);
193  relCalibNDFlatSyst->AddShift(relCalibNDFlat);
194  nus_systs_maker->AddSystematic(relCalibNDFlatSyst);
195 
196  NusSystMaker* relCalibNDShapeSyst = new NusSystMaker("RelCalibNDShape",
197  "Relative Calibration ND Flat Down");
198  NusSystLoaderShift* relCalibNDShape = new NusSystLoaderShift("RelCalibNDShape",
199  "Relative Calibration ND Shape");
200  relCalibNDShape->AddOnOff(loaders_calib_shape,
201  loaders_nom);
202  relCalibNDShapeSyst->AddShift(relCalibNDShape);
203  nus_systs_maker->AddSystematic(relCalibNDShapeSyst);
204 
205  NusSystMaker* relCalibFDFlatSyst = new NusSystMaker("RelCalibFDFlat",
206  "Relative Calibration FD Flat");
207  NusSystLoaderShift* relCalibFDFlat = new NusSystLoaderShift("RelCalibFDFlat",
208  "Relative Calibration FD Flat");
209  relCalibFDFlat->AddSigma(loaders_calib_up,
210  loaders_nom,
211  +1);
212  relCalibFDFlat->AddSigma(loaders_calib_down,
213  loaders_nom,
214  -1);
215  relCalibFDFlatSyst->AddShift(relCalibFDFlat);
216  nus_systs_maker->AddSystematic(relCalibFDFlatSyst);
217 
218  NusSystMaker* relCalibFDShapeSyst = new NusSystMaker("RelCalibFDShape",
219  "Relative Calibration FD Flat Down");
220  NusSystLoaderShift* relCalibFDShape = new NusSystLoaderShift("RelCalibFDShape",
221  "Relative Calibration FD Shape");
222  relCalibFDShape->AddOnOff(loaders_calib_shape,
223  loaders_nom);
224  relCalibFDShapeSyst->AddShift(relCalibFDShape);
225  nus_systs_maker->AddSystematic(relCalibFDShapeSyst);
226 
227  // Light level systematic
228  NusSystMaker* lightLevelSyst = new NusSystMaker("LightLevel", "Light level");
229  NusSystLoaderShift* lightLevelShift = new NusSystLoaderShift("LightLevel", "Light level");
230  lightLevelShift->AddSigma(loaders_light_up,
231  +1);
232  lightLevelShift->AddSigma(loaders_light_down,
233  -1);
234  lightLevelSyst->AddShift(lightLevelShift);
235  nus_systs_maker->AddSystematic(lightLevelSyst);
236 
237  // Cherenkov systematic
238  NusSystMaker* cherenkovSyst = new NusSystMaker("Cherenkov", "Cherenkov");
239  NusSystLoaderShift* cherenkovShift = new NusSystLoaderShift("Cherenkov", "Cherenkov");
240  cherenkovShift->AddOnOff(loaders_cherenkov);
241  cherenkovSyst->AddShift(cherenkovShift);
242  nus_systs_maker->AddSystematic(cherenkovSyst);
243 
244  // Kaon systematic
245  NusSystMaker* kaonSyst = new NusSystMaker("Kaon", "Kaon");
246  const Var kWeightKaonUp = kReweight * kNus18BeamWeightP;
247  const Var kWeightKaonDown = kReweight * kNus18BeamWeightM;
248  NusSystWeightShift* kaonShift = new NusSystWeightShift("Kaon", "Kaon");
249  kaonShift->AddSigma(&kWeightKaonUp,
250  +1);
251  kaonShift->AddSigma(&kWeightKaonDown,
252  -1);
253  kaonSyst->AddShift(kaonShift);
254  nus_systs_maker->AddSystematic(kaonSyst);
255 
256  // // Neutron systematic
257  // NusSystMaker* neutronSyst = new NusSystMaker("Neutron", "Neutron");
258  // NusSystSystShift* neutronShift = new NusSystSystShift("Neutron", "Neutron Primaries");
259  // for (const auto& sigma : sigmas)
260  // neutronShift->AddSigma(new SystShifts(&kNeutronVisEScalePrimariesSyst2018, sigma),
261  // sigma);
262  // neutronSyst->AddShift(neutronShift);
263  // nus_systs_maker->AddSystematic(neutronSyst);
264 
265  // Tau systematic
266  NusSystMaker* tauSyst = new NusSystMaker("Tau", "Tau");
267  const Var kWeightTauUp = kReweight * kFDTauWgtP;
268  const Var kWeightTauDown = kReweight * kFDTauWgtM;
269  NusSystWeightShift* tauShift = new NusSystWeightShift("Tau", "Tau Oscillations");
270  tauShift->AddSigma(&kWeightTauUp,
271  +1);
272  tauShift->AddSigma(&kWeightTauDown,
273  -1);
274  tauSyst->AddShift(tauShift);
275  nus_systs_maker->AddSystematic(tauSyst);
276 
277  // X-sec tune systematic
278  NusSystMaker* xsecTuneSyst = new NusSystMaker("XSecTune", "Cross-Section Tune");
279  const Var kNoXSecTuneWeight = kPPFXFluxCVWgt;
280  NusSystWeightShift* xsecTuneShift = new NusSystWeightShift("XSecTune", "Cross-Section Tune On/Off");
281  xsecTuneShift->AddOnOff(&kNoXSecTuneWeight);
282  xsecTuneSyst->AddShift(xsecTuneShift);
283  nus_systs_maker->AddSystematic(xsecTuneSyst);
284 
285  // Run loaders
286  loaders_nom ->Go();
287  loaders_light_nom ->Go();
288  loaders_light_up ->Go();
289  loaders_light_down ->Go();
290  loaders_cherenkov ->Go();
291  loaders_calib_up ->Go();
292  loaders_calib_down ->Go();
293  loaders_calib_shape->Go();
294 
295  // Save to file
296  TFile* outFile = new TFile("Nus18SystMaker.root", "RECREATE");
297  nus_systs_maker->SaveTo(outFile, "Nus18Systs"), true;
298  outFile->Close();
299  delete outFile;
300 
301  return;
302 
303 }
_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
Beam systematics:
Definition: BeamSysts.h:38
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
std::map< std::string, std::pair< GENIEWeightLabel, rwgt::ReweightLabel_t > > GetXSecShiftLabels()
Get cross-section systematics labels.
std::string FindCAFAnaDir()
Definition: Utilities.cxx:203
const Cut kNus18FD
Definition: NusCuts18.h:100
std::map< std::string, std::pair< GENIEWeightLabel, rwgt::ReweightLabel_t > > GetGENIEShiftLabels()
Get GENIE systematics labels.
const Var kNus18BeamWeightP
Definition: NusVars.h:87
_Var< caf::SRProxy > Var
Representation of a variable to be retrieved from a caf::StandardRecord object.
Definition: Var.h:74
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
Loaders for calibration shape paths/definitions.
Definition: Prod4Loaders.h:196
std::vector< Var > GetPPFXFluxWeights()
Loaders for absolute calibration paths/definitions.
Definition: Prod4Loaders.h:119
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
Loaders for light level paths/definitions.
Definition: Prod4Loaders.h:145
TFile * outFile
Definition: PlotXSec.C:135
std::map< std::string, std::vector< std::string > > GetBeamTranspShiftLabels()
Get Beam Transport systematics labels.
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)
const Var kXSecCVWgt2018
Definition: XsecTunes.h:48
const Var kFDTauWgtP
Definition: NusVars.h:99
assert(nhit_max >=nhit_nbins)
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
const Var kFDTauWgtM
Definition: NusVars.h:96
This module creates Common Analysis Files.
Definition: FileReducer.h:10
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod4Loaders.h:96
const NOvARwgtSyst * GetGenieKnobSyst(rwgt::ReweightLabel_t knobIdx, std::string altName, std::string altLabel)
Convenience function to get a GENIE knob syst. (Allows using the GENIE knob name & description as the...
Definition: XSecSysts.cxx:119
ECAFType
Definition: Loaders.h:19
void NusSystMaker()
Definition: NusSystMaker.C:81
enum BeamMode string