MakeSysts.C
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////
2 // MakeSysts.C
3 // Jeremy Hewes (jhewes15@fnal.gov)
4 //
5 // Make all systematics for 2019 NC covariance matrix analysis
6 //////////////////////////////////////////////////////////////
7 
10 #include "CAFAna/Core/Loaders.h"
11 #include "CAFAna/Core/SystShifts.h"
12 #include "CAFAna/Core/Utilities.h"
13 #include "CAFAna/Core/Var.h"
14 #include "CAFAna/Cuts/Cuts.h"
15 #include "NuXAna/Cuts/NusCuts18.h"
17 #include "CAFAna/Cuts/SpillCuts.h"
24 #include "CAFAna/Systs/BeamSysts.h"
25 #include "CAFAna/Systs/Systs.h"
28 #include "CAFAna/Systs/XSecSysts.h"
30 #include "CAFAna/Vars/HistAxes.h"
31 #include "NuXAna/Vars/HistAxes.h"
33 #include "NuXAna/Analysis/NusSystsMaker.h"
34 
36 
37 #include "TRandom.h"
38 
39 using namespace ana;
40 
41 // -------------------------------------------------------------------------------------
42 void MakeSysts(TString opt) {
43 
44  DontAddDirectory guard;
45 
46  bool fd = opt.Contains("fardet", TString::kIgnoreCase);
47  bool nd = opt.Contains("neardet", TString::kIgnoreCase);
48  assert(fd xor nd);
49 
50  std::string name = "nus_fhc_";
51 
52  if (nd) name = name + "neardet";
53  else name = name + "fardet";
54 
55  // Get POT
56  double pot = 12e20;
57 
59 
60  Loaders* loaders_nom = new Prod4NomLoaders (ECAFType::kFullCAF, fluxType);
61  Loaders* loaders_light_nom = new Prod4LightLevelLoaders(ECAFType::kFullCAF, fluxType, 0);
62  Loaders* loaders_light_up = new Prod4LightLevelLoaders(ECAFType::kFullCAF, fluxType, +1);
63  Loaders* loaders_light_down = new Prod4LightLevelLoaders(ECAFType::kFullCAF, fluxType, -1);
64  Loaders* loaders_cherenkov = new Prod4CherenkovLoaders (ECAFType::kFullCAF, fluxType);
65  Loaders* loaders_calib_up = new Prod4AbsCalibLoaders (ECAFType::kFullCAF, fluxType, +1);
66  Loaders* loaders_calib_down = new Prod4AbsCalibLoaders (ECAFType::kFullCAF, fluxType, -1);
67  Loaders* loaders_calib_shape = new Prod4CalibShapeLoaders(ECAFType::kFullCAF, fluxType);
68 
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  if (nd) {
79  loaders_nom ->SetND(true);
80  loaders_light_nom ->SetND(true);
81  loaders_light_up ->SetND(true);
82  loaders_light_down ->SetND(true);
83  loaders_cherenkov ->SetND(true);
84  loaders_calib_up ->SetND(true);
85  loaders_calib_down ->SetND(true);
86  loaders_calib_shape->SetND(true);
87 
90 
93 
96 
99 
102 
105  }
106 
107  // Set up the weights
108  const Var kReweight = kXSecCVWgt2018 * kPPFXFluxCVWgt;
109 
110  // Sigma variations for any syst shifts
111  std::vector<int> sigmas = {-3, -2, -1, 0, 1, 2, 3};
112 
113  // Get axis and cuts
114  HistAxis* nus_axis = nullptr;
115  if (nd) nus_axis = new HistAxis(kNCNDAxisE);
116  else if (fd) nus_axis = new HistAxis(kNCFDAxisE);
117  HistAxis* numu_axis = new HistAxis(kNus18BinsNumuCCAxis);
118  Cut* fd_cut = new Cut(kNus18FD);
119  Cut* nd_cut = new Cut(kNus18ND);
120  Cut* numu_cut = new Cut(kNumuCutND2018);
121 
122  // Get prediction type
123  // NB pred_type is currently passed to NusSystsMaker, so each NusSystsMaker
124  // can only handle one prediction type.
125  // It could instead be passed to NusSystMaker::AddSystematic, to be used when
126  // making NusSystsShiftMaker (this is actually all that happens under the hood
127  // right now), if we wanted the maker to be able to handle systematics from
128  // multiple sources.
129  NusSystsPredType pred_type = NusSystsPredType::kUnknown;
130  if (fd) pred_type = NusSystsPredType::kFD;
131  if (nd) pred_type = NusSystsPredType::kND;
132 
133  // Set up systematic maker
134  NusSystematicsMaker* nus_systs_maker = new NusSystematicsMaker("nc_fhc_2019", pred_type);
135  nus_systs_maker->SetAxes(nus_axis, numu_axis);
136  nus_systs_maker->SetCuts(fd_cut, nd_cut, numu_cut);
137  if (nd) nus_systs_maker->SetNominalLoaders(loaders_nom, loaders_nom);
138  else nus_systs_maker->SetNominalLoaders(loaders_nom, loaders_light_nom);
139  nus_systs_maker->SetNominalWeight(&kReweight);
140 
141  // Beam transport systematic
142  std::map<std::string, std::vector<std::string> > beamShiftLabels = GetBeamTranspShiftLabels();
143  std::string beamShiftFile = FindCAFAnaDir() + "/data/beam/TABeamSyst.root";
144 
145  // Let's also do all the beam transport systematics separately
146  for (const auto& shiftLabel : beamShiftLabels) {
147  NusSystMaker* beamTransSyst = new NusSystMaker(shiftLabel.first, shiftLabel.second[1], name);
148  NusSystSystShift* beamTransShift = new NusSystSystShift(shiftLabel.first, shiftLabel.second[1]);
149  for (const auto& sigma : sigmas)
150  beamTransShift->AddSigma(new SystShifts(new BeamSyst(beamShiftFile, shiftLabel.second[0], shiftLabel.second[1]), sigma),
151  sigma);
152  beamTransSyst->AddShift(beamTransShift);
153  nus_systs_maker->AddSystematic(beamTransSyst);
154  }
155 
156  // Calibration systematics
157 
158  NusSystMaker* absCalibSyst = new NusSystMaker("AbsCalib", "Absolute Calibration", name);
159  NusSystLoaderShift* absCalib = new NusSystLoaderShift("AbsCalib", "Absolute Calibration");
160  absCalib->AddSigma(loaders_calib_up, +1);
161  absCalib->AddSigma(loaders_calib_down, -1);
162  absCalibSyst->AddShift(absCalib);
163  nus_systs_maker->AddSystematic(absCalibSyst);
164 
165  NusSystMaker* absCalibShapeSyst = new NusSystMaker("AbsCalibShape", "Absolute Calibration Shape", name);
166  NusSystLoaderShift* absCalibShape = new NusSystLoaderShift("AbsCalibShape", "Absolute Calibration Shape");
167  absCalibShape->AddSigma(loaders_calib_shape, +1);
168  absCalibShapeSyst->AddShift(absCalibShape);
169  nus_systs_maker->AddSystematic(absCalibShapeSyst);
170 
171  // Light level systematics
172 
173  NusSystMaker* lightLevelSyst = new NusSystMaker("LightLevel", "Light Level", name);
174  NusSystLoaderShift* lightLevelShift = new NusSystLoaderShift("LightLevel", "Light Level");
175  lightLevelShift->AddSigma(loaders_light_up, +1);
176  lightLevelShift->AddSigma(loaders_light_down, -1);
177  lightLevelSyst->AddShift(lightLevelShift);
178  nus_systs_maker->AddSystematic(lightLevelSyst);
179 
180  // Cherenkov systematic
181 
182  NusSystMaker* cherenkovSyst = new NusSystMaker("Cherenkov", "Cherenkov", name);
183  NusSystLoaderShift* cherenkovShift = new NusSystLoaderShift("Cherenkov", "Cherenkov");
184  cherenkovShift->AddSigma(loaders_cherenkov, +1);
185  cherenkovSyst->AddShift(cherenkovShift);
186  nus_systs_maker->AddSystematic(cherenkovSyst);
187 
188  // Kaon systematic
189 
190  NusSystMaker* kaonSyst = new NusSystMaker("Kaon", "Kaon", name);
191  const Var kWeightKaonUp = kReweight * kNus18BeamWeightP;
192  const Var kWeightKaonDown = kReweight * kNus18BeamWeightM;
193  NusSystWeightShift* kaonShift = new NusSystWeightShift("Kaon", "Kaon");
194  kaonShift->AddSigma(&kWeightKaonUp, +1);
195  kaonShift->AddSigma(&kWeightKaonDown, -1);
196  kaonSyst->AddShift(kaonShift);
197  nus_systs_maker->AddSystematic(kaonSyst);
198 
199  // Neutron systematic
200 
201  NusSystMaker* neutronSyst = new NusSystMaker("Neutron", "Neutron", name);
202  NusSystSystShift* neutronShift = new NusSystSystShift("Neutron", "Neutron Primaries");
203  for (const auto& sigma : sigmas)
204  neutronShift->AddSigma(new SystShifts(&kNeutronVisEScalePrimariesSyst2018, sigma),
205  sigma);
206  neutronSyst->AddShift(neutronShift);
207  nus_systs_maker->AddSystematic(neutronSyst);
208 
209  // X-sec tune systematic
210 
211  NusSystMaker* xsecTuneSyst = new NusSystMaker("XSecTune", "Cross-Section Tune", name);
212  const Var kNoXSecTuneWeight = kPPFXFluxCVWgt;
213  NusSystWeightShift* xsecTuneShift = new NusSystWeightShift("XSecTune", "Cross-Section Tune On/Off");
214  xsecTuneShift->AddSigma(&kNoXSecTuneWeight, +1);
215  xsecTuneSyst->AddShift(xsecTuneShift);
216  nus_systs_maker->AddSystematic(xsecTuneSyst);
217 
218  // Run loaders
219  loaders_nom ->Go();
220  if (!nd) loaders_light_nom ->Go();
221  loaders_light_up ->Go();
222  loaders_light_down ->Go();
223  loaders_cherenkov ->Go();
224  loaders_calib_up ->Go();
225  loaders_calib_down ->Go();
226  loaders_calib_shape->Go();
227 
228  // Save to file
229  TFile* outfile = new TFile("systsmaker_nus_fhc.root", "RECREATE");
230  nus_systs_maker->SaveTo(outfile, "systsmaker_nus_fhc"), true;
231  if (outfile) outfile->Close();
232 
233 }
Near Detector underground.
Definition: SREnums.h:10
const XML_Char * name
Definition: expat.h:151
_HistAxis< Var > HistAxis
Definition: HistAxis.h:103
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void MakeSysts(TString opt)
Definition: MakeSysts.C:42
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
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
std::string FindCAFAnaDir()
Definition: Utilities.cxx:204
const Cut kNus18FD
Definition: NusCuts18.h:100
const Var kNus18BeamWeightP
Definition: NusVars.h:87
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
Loaders for calibration shape paths/definitions.
Definition: Prod4Loaders.h:196
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
std::map< std::string, std::vector< std::string > > GetBeamTranspShiftLabels()
Get Beam Transport systematics labels.
SpectrumLoaderBase & GetLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Retrieve a specific loader.
Definition: Loaders.cxx:129
#define pot
_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
void SetND(bool nd)
Definition: Loaders.h:65
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
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod4Loaders.h:96
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
FILE * outfile
Definition: dump_event.C:13
void NusSystMaker()
Definition: NusSystMaker.C:81
enum BeamMode string