MakeSysts.C
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////
2 // MakeSysts.C
3 // Jeremy Hewes (jhewes15@fnal.gov)
4 //
5 // Make systematics for 2020 sterile covariance matrix analysis
6 //////////////////////////////////////////////////////////////
7 
13 
17 
18 #include "NuXAna/Cuts/NusCuts20.h"
21 
22 using std::pair;
23 using std::string;
24 using std::vector;
25 
26 using ana::covmx::Sample;
27 
28 using namespace ana;
29 
30 // -------------------------------------------------------------------------------------
31 void MakeSysts(TString opt) {
32 
33  DontAddDirectory guard;
34 
35  // Define exposures
36  double potND = kNux20FHCNDPOT;
37  double potFD = kNux20FHCFDPOT;;
38 
39  Sample sample = GetSampleFromOptString(opt);
40 
42 
44  if (sample.selection == covmx::kCCNue) CAFType = ECAFType::kNueConcat;
45  else if (sample.selection == covmx::kCCNumu) CAFType = ECAFType::kNumuConcat;
46  else if (sample.selection == covmx::kNC) CAFType = ECAFType::kNusConcat;
47 
48  // Load quantiles from file if we have any numu
49 
50  vector<Cut> HadEFracQuantCuts;
51  if (sample.selection == covmx::kCCNumu && sample.quantile != covmx::kNoQ) {
52  string quantDir = std::getenv("NUMUDATA_DIR");
53  string infile_quant = quantDir + "/lib/ana2020/Quantiles/quantiles_" + sample.GetPolarity() + "_full_numu2020.root";
54  TFile* infile = TFile::Open( pnfs2xrootd(infile_quant).c_str() );
55  TH2* spec2DFD = (TH2*)infile->FindObjectAny("FDSpec2D");
56  HadEFracQuantCuts = QuantileCutsFromTH2(spec2DFD, kNumuCCOptimisedAxis2020, kHadEFracAxis, 4);
57  infile->Close();
58  } // if numu selection found
59 
60  Prod5NomLoaders* loadersNom = new Prod5NomLoaders (CAFType, fluxType);
61  Prod5LightLevelLoaders* loadersLightUp = new Prod5LightLevelLoaders(CAFType, fluxType, +1);
62  Prod5LightLevelLoaders* loadersLightDown = new Prod5LightLevelLoaders(CAFType, fluxType, -1);
63  Prod5CherenkovLoaders* loadersCherenkov = new Prod5CherenkovLoaders (CAFType, fluxType);
64  Prod5CalibDriftLoaders* loadersCalibDrift = new Prod5CalibDriftLoaders(CAFType, fluxType);
65  Prod5AbsCalibLoaders* loadersCalibUp = new Prod5AbsCalibLoaders (CAFType, fluxType, +1);
66  Prod5AbsCalibLoaders* loadersCalibDown = new Prod5AbsCalibLoaders (CAFType, fluxType, -1);
67  Prod5CalibShapeLoaders* loadersCalibShape = new Prod5CalibShapeLoaders(CAFType, fluxType);
68 
69  loadersNom ->SetSpillCut(kStandardSpillCuts);
70  loadersLightUp ->SetSpillCut(kStandardSpillCuts);
71  loadersLightDown ->SetSpillCut(kStandardSpillCuts);
72  loadersCherenkov ->SetSpillCut(kStandardSpillCuts);
73  loadersCalibDrift->SetSpillCut(kStandardSpillCuts);
74  loadersCalibUp ->SetSpillCut(kStandardSpillCuts);
75  loadersCalibDown ->SetSpillCut(kStandardSpillCuts);
76  loadersCalibShape->SetSpillCut(kStandardSpillCuts);
77 
78  if (sample.detector == covmx::kNearDet) {
79  loadersNom ->SetND(true);
80  loadersLightUp ->SetND(true);
81  loadersLightDown ->SetND(true);
82  loadersCherenkov ->SetND(true);
83  loadersCalibDrift->SetND(true);
84  loadersCalibUp ->SetND(true);
85  loadersCalibDown ->SetND(true);
86  loadersCalibShape->SetND(true);
87  }
88 
89  // Set up the weights
90  Var kReweight = kPPFXFluxCVWgt;
91 
92  // Get POT
93  double pot = (sample.detector == covmx::kNearDet) ? potND : potFD;
94 
95  // Sigma variations for any syst shifts
96  vector<int> sigmas = {-3, -2, -1, 0, 1, 2, 3};
97 
98  const HistAxis* axis = nullptr;
99  const HistAxis* ncNDnumuAxis = nullptr;
100  const Cut* ndCut = nullptr;
101  const Cut* fdCut = nullptr;
102  const Cut* ncNDnumuCut = nullptr;
103 
104  // Neutral current axis and cuts
105  if (sample.selection == covmx::kNC) {
106  if (sample.detector == covmx::kNearDet) axis = &kNC20NDAxisE;
107  else axis = &kNC20FDAxisE;
108  ndCut = &kNus20NDCuts;
109  fdCut = &kNus20FDCuts_ML;
110  }
111 
112  // Numu axis and cuts
113  else if (sample.selection == covmx::kCCNumu) {
114  axis = &kNumuCCOptimisedAxis2020;
115  if (sample.quantile == covmx::kNoQ) {
116  ndCut = &kNumu2020ND;
117  fdCut = &kNumu2020FD_ML;
118  } else {
119  ndCut = new const Cut(kNumu2020ND && HadEFracQuantCuts[sample.quantile]);
120  fdCut = new const Cut(kNumu2020FD_ML && HadEFracQuantCuts[sample.quantile]);
121  }
122  }
123 
124  // Nue axis and cuts
125  else if (sample.selection == covmx::kCCNue) {
126  axis = &kNue2020Axis; /// kNue2020AxisMergedPeripheral if Peripheral is collapsed, e.g. for FD only predictions
127  ndCut = &kNue2020ND;
128  fdCut = &kNue2020FDAllSamples_ML;
129  }
130 
131  // Get prediction type
132  // NB predType is currently passed to NusSystsMaker, so each NusSystsMaker
133  // can only handle one prediction type.
134  // It could instead be passed to SystMaker::AddSystematic, to be used when
135  // making NusSystsShiftMaker (this is actually all that happens under the hood
136  // right now), if we wanted the maker to be able to handle systematics from
137  // multiple sources.
139  if (sample.detector == covmx::kFarDet) predType = SystPredType::kFD;
140  else predType = SystPredType::kND;
141 
142  // Set up systematic maker
143  SystematicsMaker* systsMaker = new SystematicsMaker(sample.GetTag().c_str(), predType);
144 
145  systsMaker->SetAxes(axis, ncNDnumuAxis);
146  systsMaker->SetCuts(fdCut, ndCut, ncNDnumuCut);
147  systsMaker->SetNominalLoaders(loadersNom, loadersNom);
148  systsMaker->SetNominalWeight(&kReweight);
149 
150  // Calibration systematics
151  SystMaker* absCalibSyst = new SystMaker("AbsCalib", "Absolute Calibration", sample.GetTag());
152  SystMakerLoaderShift* absCalib = new SystMakerLoaderShift("AbsCalib", "Absolute Calibration");
153  absCalib->AddSigma(loadersCalibUp, +1);
154  absCalib->AddSigma(loadersCalibDown, -1);
155  absCalibSyst->AddShift(absCalib);
156  systsMaker->AddSystematic(absCalibSyst);
157 
158  SystMaker* shapeCalibSyst = new SystMaker("ShapeCalib", "Shape Calibration", sample.GetTag());
159  SystMakerLoaderShift* shapeCalib = new SystMakerLoaderShift("ShapeCalib", "Shape Calibration");
160  shapeCalib->AddSigma(loadersCalibShape, +1);
161  shapeCalibSyst->AddShift(shapeCalib);
162  systsMaker->AddSystematic(shapeCalibSyst);
163 
164  // Light level systematic
165  SystMaker* lightLevelSyst = new SystMaker("LightLevel", "Light Level", sample.GetTag());
166  SystMakerLoaderShift* lightLevelShift = new SystMakerLoaderShift("LightLevel", "Light Level");
167  lightLevelShift->AddSigma(loadersLightUp, +1);
168  lightLevelShift->AddSigma(loadersLightDown, -1);
169  lightLevelSyst->AddShift(lightLevelShift);
170  systsMaker->AddSystematic(lightLevelSyst);
171 
172  // Cherenkov systematic
173  SystMaker* cherenkovSyst = new SystMaker("Cherenkov", "Cherenkov", sample.GetTag());
174  SystMakerLoaderShift* cherenkovShift = new SystMakerLoaderShift("Cherenkov", "Cherenkov");
175  cherenkovShift->AddSigma(loadersCherenkov, +1);
176  cherenkovSyst->AddShift(cherenkovShift);
177  systsMaker->AddSystematic(cherenkovSyst);
178 
179  // Drift systematic
180  SystMaker* calibDriftSyst = new SystMaker("CalibDrift", "CalibDrift", sample.GetTag());
181  SystMakerLoaderShift* calibDriftShift = new SystMakerLoaderShift("CalibDrift", "CalibDrift");
182  calibDriftShift->AddSigma(loadersCalibDrift, +1);
183  calibDriftSyst->AddShift(calibDriftShift);
184  systsMaker->AddSystematic(calibDriftSyst);
185 
186  cout << "Running Go() on all loaders." << endl;
187 
188  // Run loaders
189  loadersNom ->Go();
190  loadersLightUp ->Go();
191  loadersLightDown ->Go();
192  loadersCherenkov ->Go();
193  loadersCalibDrift->Go();
194  loadersCalibUp ->Go();
195  loadersCalibDown ->Go();
196  loadersCalibShape->Go();
197 
198  // Save SystematicsMaker objects
199  TFile* outfile = new TFile("systmakers.root", "recreate");
200  systsMaker->SaveTo(outfile, sample.GetTag());
201  delete outfile; // clean up
202 
203 }
204 
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
const HistAxis kHadEFracAxis("E_{had.} / E_{#nu}", Binning::Simple(200, 0, 1), kHadEFrac)
HistAxis that implements Hadronic Energy fraction binning used by L.Vinton to derive Hadronic Energy ...
Definition: HistAxes.h:30
const Cut kNumu2020ND
Definition: NumuCuts2020.h:57
std::string GetTag() const
Definition: Sample.cxx:94
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
const HistAxis kNumuCCOptimisedAxis2020("Reconstructed Neutrino Energy (GeV)", kNumuCCEOptimisedBinning, kNumuE2020)
Definition: HistAxes.h:25
covmx::Sample GetSampleFromOptString(TString optString)
Function to take an option TString and return a single associated covmx::Sample.
Definition: Utilities.h:250
const Cut kNue2020ND
Definition: NueCuts2020.h:178
const double kNux20FHCFDPOT
Definition: Utilities.h:41
const Cut kNus20NDCuts
Definition: NusCuts20.h:102
void AddShift(SystMakerShift *shift)
Add SystMakerShift to this syst maker.
Definition: SystMaker.cxx:1015
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
const HistAxis kNC20FDAxisE("Energy Deposited in Scintillator (GeV)", kNCFDBinning, kNus20Energy)
Definition: HistAxes.h:21
Loaders for absolute calibration paths/definitions.
Definition: Prod5Loaders.h:122
Loaders for Cherenkov paths/definitions.
Definition: Prod5Loaders.h:173
void SetAxes(const HistAxis *axis, const HistAxis *numuAxis=nullptr)
Definition: SystMaker.cxx:1611
std::string GetPolarity() const
Definition: Sample.cxx:52
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
string infile
void SetNominalLoaders(Loaders *loader, Loaders *loader_light=nullptr)
Definition: SystMaker.cxx:1625
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod5Loaders.h:101
std::string getenv(std::string const &name)
Selection selection
Definition: Sample.h:95
const HistAxis kNue2020Axis("NuE Energy / Analysis Bin", kNue2020Binning, kNue2020AnaBin)
Use this Axis for Ana2020, official Axis.
Definition: NueCuts2020.h:195
void SetCuts(const Cut *fd_cut, const Cut *nd_cut, const Cut *numu_cut=nullptr)
Definition: SystMaker.cxx:1617
Detector detector
Definition: Sample.h:97
#define pot
void AddSystematic(SystMaker *syst)
Definition: SystMaker.cxx:1573
Loaders for calibration drift (detector aging) paths/definitions.
Definition: Prod5Loaders.h:218
OStream cout
Definition: OStream.cxx:6
Loaders for light level paths/definitions.
Definition: Prod5Loaders.h:147
const HistAxis kNC20NDAxisE("Energy Deposited in Scintillator (GeV)", kNCNDBinning, kNus20Energy)
NC covariance hist axes.
Definition: HistAxes.h:20
Quantile quantile
Definition: Sample.h:98
_Cut< caf::SRProxy, caf::SRSpillProxy > Cut
Definition: Cut.h:9
void SetND(bool nd)
Definition: Loaders.h:65
std::vector< Cut > QuantileCutsFromTH2(TH2 *quantileHist, const HistAxis &independentAxis, const HistAxis &quantileAxis, const unsigned int &numQuantiles, const bool verbose)
: Do the same as the QuantileCuts function but taking in the TH2 instead of making it...
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
void SetNominalWeight(const Var *weight)
Definition: SystMaker.cxx:1631
const Cut kNue2020FDAllSamples_ML
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
const Cut kNus20FDCuts_ML
Polarity polarity
Definition: Sample.h:96
Loaders for calibration shape paths/definitions.
Definition: Prod5Loaders.h:196
const Cut kNumu2020FD_ML
const double kNux20FHCNDPOT
Definition: Utilities.h:40
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: SystMaker.cxx:1648
FILE * outfile
Definition: dump_event.C:13
SystPredType
Definition: SystMaker.h:38
ECAFType
Definition: Loaders.h:19
enum BeamMode string