NueSystHelper.cxx
Go to the documentation of this file.
2 
4 #include "CAFAna/Systs/Systs.h"
14 #include "3FlavorAna/Systs/GeniePCASyst.h" // GetGeniePrincipals2018
16 #include "3FlavorAna/Systs/MichelTaggingSyst.h"//MicheclSystematic
18 
20 
21 #include <iostream>
22 
23 namespace ana
24 {
25 
26 // Functions compied from syst_variations
27 
28  std::vector <const ISyst*> GetSystematics2018(const TString option)
29  {
30  std::vector <const ISyst*> ret;
31  std::vector<const ISyst*> tmp;
32 
33  if (option.Contains("RelativeCalib") || option.Contains("All"))
34  ret.push_back(&kAnaRelativeCalibSyst);
35  if (option.Contains("Calibration") || option.Contains("All"))
36  ret.push_back(&kAnaCalibrationSyst);
37  if (option.Contains("CalibShape") || option.Contains("All"))
38  ret.push_back(&kAnaCalibShapeSyst);
39  if (option.Contains("LightLevel") || option.Contains("All"))
40  ret.push_back(&kAnaLightlevelSyst);
41  if (option.Contains("Cherenkov") || option.Contains("All"))
42  ret.push_back(&kAnaCherenkovSyst);
43 
44  if(option.Contains("BeamTransport"))
45  {
46  std::cerr << "Option 'BeamTransport' is deprecated." << std::endl;
47  std::cerr << "'FluxSyst' includes both 'BeamTransport' and 'PPFXSyst'"
48  << std::endl;
49  abort(); // outdated option
50  }
51  if(option.Contains("BeamAllTransport"))
52  {
53  std::cerr << "Option 'BeamAllTransport' is deprecated." << std::endl;
54  std::cerr << "'FluxSyst' includes both 'BeamTransport' and 'PPFXSyst'"
55  << std::endl;
56  abort(); // outdated option
57  }
58  if(option.Contains("PPFXSyst"))
59  {
60  std::cerr << "Option 'PPFXSyst' is deprecated." << std::endl;
61  std::cerr << "'FluxSyst' includes both 'BeamTransport' and 'PPFXSyst'"
62  << std::endl;
63  abort(); // outdated option
64  }
65 
66  if(option.Contains("FluxSyst") || option.Contains("All"))
67  {
68  // Selecting 8 principal components to include in the big PDF
69  // In the final fit we expect this number to be smaller ~4-5
70  // TODO: Add functionality to specify desired number of components at
71  // runtime.
72  for (int i = 0; i < 8; i++)
73  ret.push_back(GetFluxPrincipals2018(i));
74  }
75  if(option.Contains("GENIEPCASyst") || option.Contains("All"))
76  {
77  // Selecting 5 principal components to include in the big PDF
78  // This matches the number we think we will include in the final fit
79  // TODO: Add functionality to specify desired number of components at
80  // runtime.
81  for (int i = 0; i < 5; i++)
82  ret.push_back(GetGeniePrincipals2018Small(i));
83  }
84  if(option.Contains("XSecSyst") || option.Contains("All"))
85  {
86  std::vector<const ISyst*> geniesysts = getAllXsecSysts_2018();
87  // After 2018 analysis some new systematics were added. This returns the
88  // 2018 systematics list plus 5 new systematics:
89  // - kRPACCQEEnhSyst2019
90  // - kRPACCQESuppSyst2019
91  // - kRPARESSyst2019
92  // - kMECShapeSyst2018RPAFixNu
93  // - kMECShapeSyst2018RPAFixAntiNu
94  // ONLY use if you are also applying the RPAFix weights
95  if (option.Contains("RPAFix")){
96  std::cout<<"Replacing 2018RPA systs by RPAFix..."<<std::endl;
97  geniesysts.clear();
98  geniesysts = getAllXsecSysts_2018_RPAFix();
99  }
100  // The tau systematic description is in the Nue executive summary,
101  // p. 25 of doc-26699-v9
102  geniesysts.push_back(&kTauScaleSyst);
103  int minsyst=0, maxsyst=geniesysts.size();
104  if(option.Contains("XSecSyst_Frac"))
105  {
106  minsyst = 0;
107  maxsyst = 3;
108  std::cout << "\n\nOnly looking at a small range: " << std::endl;
109  }
110  if(option.Contains("XSecSyst_RPATest"))
111  {
112  std::cout << "\n\nWe are not using kRPARESSyst2018Test" << std::endl;
113  }
114  if(option.Contains("XSecSyst_COH"))
115  {
116  minsyst = 53;
117  maxsyst = 55;
118  std::cout << "\n\nOnly looking at COHCC and COHNC (which replaced MaCOHpi and R0COHpi): " << std::endl;
119  }
120  if(option.Contains("XSecSyst_1of6")){ minsyst= 0; maxsyst=14;}
121  if(option.Contains("XSecSyst_2of6")){ minsyst=14; maxsyst=28;}
122  if(option.Contains("XSecSyst_3of6")){ minsyst=28; maxsyst=42;}
123  if(option.Contains("XSecSyst_4of6")){ minsyst=42; maxsyst=56;}
124  if(option.Contains("XSecSyst_5of6")){ minsyst=56; maxsyst=70;}
125  if(option.Contains("XSecSyst_6of6")){ minsyst=70; maxsyst=87;}
126  if(option.Contains("All")){ minsyst=0; maxsyst=87;}
127 
128  std::cout << minsyst << " -- " << maxsyst << "\n\n";
129 
130  for(int i = minsyst; i < maxsyst; ++i)
131  {
132  if (i >= (int)geniesysts.size()) break;
133  ret.push_back(geniesysts[i]);
134  }
135  }
136  if (option.Contains("MuEnergySyst") || (option.Contains("All")))
137  {
138  // 1.0% abs (FD + ND) muon E scale applied to muon track length
139  ret.push_back(&kMuEScaleSyst2017); // Nue and Numu
140  if(option.Contains("numu", TString::kIgnoreCase)) // Numu only!
141  {
142  // .76% rel (FD only) muon E scale applied to muon track length
143  ret.push_back(&kRelMuEScaleSyst2017);
144  // 5% hadronic energy scale shift.
145  ret.push_back(&kDirectHadEScaleSyst2017);
146  // The relative hadronic energy scale shift
147  ret.push_back(&kDirectRelHadEScaleSyst2017);
148  }
149  }
150 
151  //Neutron systematic calculation
152  if(option.Contains("NeutronSyst")||(option.Contains("All")))
153  {
154  ret.push_back(&kNeutronVisEScalePrimariesSyst2018);
155  }
156 
157 
158  if(option.Contains ("NueExtrapSyst") ||
159  (option.Contains("All") && option.Contains("nue")))
160  {
161  tmp = {
162  //&kNueExtrapSystSignalKin2017,
163  //&kNueExtrapSystBkg2017
164  };
165  for (auto entry:tmp)
166  ret.push_back(entry);
167  }
168  if (option.Contains("SummedSmall"))
169  {
171  ret.push_back(small);
172  }
173  if (option.Contains("NuMuSmall"))
174  {
176  ret.push_back(small);
177  }
178  // Add acceptance sytematics
179  if (option.Contains("AcceptanceSyst")|| (option.Contains("All") && option.Contains("nue")))
180  {
181  tmp = {
182  &kNueAcceptSystBkg2018FHC, // background, FHC
183  &kNueAcceptSystSignalKin2018FHC, // signal, FHC
184  &kNueAcceptSystBkg2018RHC, // background, RHC
185  &kNueAcceptSystSignalKin2018RHC, // signal, RHC
186  };
187  for (auto entry:tmp)
188  ret.push_back(entry);
189  }
190 
191  // Add MichelTagging sytematics
192  if (option.Contains("MichelTaggingSyst")|| (option.Contains("All") && option.Contains("nue")))
193  {
194  ret.push_back(&kMichelTaggingSyst2018);
195  }
196 
197  // Wrong sign scale systematic
198  if(option.Contains("WrongSignScale"))
199  {
200  ret.push_back(&kWrongSignScale);
201  ret.push_back(&kWrongSignScale100);
202  }
203 
204  // Wrong sign scale by energy bin
205  if(option.Contains("WrongSignEnergyBin"))
206  {
207  ret.push_back(&kWrongSignEnergyBin);
208  ret.push_back(&kWrongSignEnergyBinRHC);
209  }
210 
211 
212  return ret;
213  }// std::vector <const ISyst*> GetSystematics2018
214 
215  std::vector <ShiDef> GetShifts2018(const TString option)
216  {
217  std::vector <ShiDef> ret;
218  TString topstr, lowstr;
219 
220  // Add nominal
221  if (option.Contains("All") || option.Contains("Nominal"))
222  {
223  topstr = "Nominal"; lowstr= "noShift";
224  ret.push_back({topstr, lowstr, kNoShift});
225  }
226  if(option.Contains("FluxSyst") ||
227  option.Contains("GENIEPCASyst") ||
228  option.Contains("XSecSyst") ||
229  option.Contains("MuEnergySyst") ||
230  option.Contains("NueExtrapSyst") ||
231  option.Contains("AllSmall") ||
232  option.Contains("NuMuSmall") ||
233  option.Contains("AcceptanceSyst") ||
234  option.Contains("NeutronSyst") ||
235  option.Contains("MichelTaggingSyst") ||
236  option.Contains("WrongSignScale") ||
237  option.Contains("WrongSignEnergyBin") ||
238  option.Contains("All"))
239  {
240  auto systs = GetSystematics2018(option);
241 
242  for (auto syst:systs){
243  for (auto sigma:{+2,+1,-1,-2})
244  {
245  TString sig_str = (sigma>0?"plus":"minus");
246  if(abs(sigma)==1) sig_str+="One";
247  if(abs(sigma)==2) sig_str+="Two";
248  ret.push_back({syst->ShortName(),sig_str,SystShifts(syst,sigma)});
249  std::cout << (TString) "Added " + syst->ShortName() + "_" + sig_str
250  << std::endl;
251  }//end sigmas
252  std::cout << std::endl;
253  }//end systs
254  }//if FluxSyst || GENIEPCASyst || XSecSyst || NumuEnergySyst ||
255  //NueExtrapSyst || AllSmall || NuMuSmall || WrongSign || DoubleWS
256 
257  if (option.Contains("CalibrationUp") || option.Contains("All"))
258  {
259  topstr = "Calibration"; lowstr ="plusOne";
260  ret.push_back({topstr, lowstr, kNoShift});
261  }
262  if(option.Contains("CalibrationDown") || option.Contains("All"))
263  {
264  topstr = "Calibration"; lowstr ="minusOne";
265  ret.push_back({topstr, lowstr, kNoShift});
266  }
267  if(option.Contains("RelativeCalibUp") || option.Contains("All"))
268  {
269  topstr = "RelativeCalib"; lowstr ="plusTwo";
270  ret.push_back({topstr, lowstr, kNoShift});
271  }
272  if(option.Contains("RelativeCalibDown") || option.Contains("All"))
273  {
274  topstr = "RelativeCalib"; lowstr ="minusTwo";
275  ret.push_back({topstr, lowstr, kNoShift});
276  }
277  if(option.Contains("CalibShape") || option.Contains("All"))
278  {
279  topstr = "CalibShape"; lowstr = "plusOne";
280  ret.push_back({topstr, lowstr, kNoShift});
281  }
282  if(option.Contains("Cherenkov") || option.Contains("All"))
283  {
284  topstr = "Cherenkov"; lowstr = "plusOne";
285  ret.push_back({topstr, lowstr, kNoShift});
286  }
287  if(option.Contains("LightLevelUp") || option.Contains("All"))
288  {
289  topstr = "Lightlevel"; lowstr = "plusOne";
290  ret.push_back({topstr, lowstr, kNoShift});
291  }
292  if(option.Contains("LightLevelDown") || option.Contains("All"))
293  {
294  topstr = "Lightlevel"; lowstr = "minusOne";
295  ret.push_back({topstr, lowstr, kNoShift});
296  }
297  if(option.Contains("LightLevelNom") || option.Contains("All"))
298  {
299  topstr = "Lightlevel"; lowstr = "noShift";
300  ret.push_back({topstr, lowstr, kNoShift});
301  }
302  if(option.Contains("CalibrationUp_2sigmaNoTau") ||
303  option.Contains("All"))
304  {
305  topstr = "Calibration"; lowstr = "plusTwo";
306  ret.push_back({topstr, lowstr, kNoShift});
307  }
308  if(option.Contains("CalibrationDown_2sigmaNoTau") ||
309  option.Contains("All"))
310  {
311  topstr = "Calibration"; lowstr = "minusTwo";
312  ret.push_back({topstr, lowstr, kNoShift});
313  }
314  return ret;
315  } //std::vector <ShiDef> GetShifts2018
316 
317  std::vector<SNameDef> GetShiftNames2018(const TString option)
318  {
319  std::vector <SNameDef> names;
320 
321  if (option.Contains("All"))
322  {
323  names.push_back({"SumSmallXSecJoint2018",
324  {"plusOne","minusOne","plusTwo","minusTwo"}});
325  names.push_back({"extrap_signalkin",
326  {"plusOne","minusOne","plusTwo","minusTwo"}});
327  names.push_back({"extrap_bkg",
328  {"plusOne","minusOne","plusTwo","minusTwo"}});
329  }
330  if((option == "Calibration") || option.Contains("All"))
331  names.push_back({"Calibration", {"plusOne","minusOne"}});
332  if((option == "RelativeCalib") || option.Contains("All"))
333  names.push_back({"RelativeCalib",{"plusTwo","minusTwo"}});
334  if((option == "CalibShape") || option.Contains("All"))
335  names.push_back({"CalibShape", {"plusOne"}});
336  if((option == "Cherenkov") || option.Contains("All"))
337  names.push_back({"Cherenkov", {"plusOne"}});
338  if((option == "Lightlevel") || option.Contains("All"))
339  names.push_back({"Lightlevel", {"plusOne","minusOne"}});
340 
341  if(option.Contains("All") && option.Contains("Reduced") )
342  {
343  std::vector <const ISyst*> systs;
344  if(!option.Contains("Super") )
345  {
346  systs = {&kMECq0ShapeSyst2017,
348  GetGenieKnobSyst(rwgt::fReweightMaCCRES),
349  GetGenieKnobSyst(rwgt::fReweightMaNCRES),
350  GetGenieKnobSyst(rwgt::fReweightMvCCRES),
354  &kRadCorrNue,
358  for (int i =0; i<5; i++)
359  {
360  systs.push_back(GetPPFXPrincipals(i));
361  }
362  }
363  // Some nue specific ones...
364  if(option.Contains("Nue") )
365  {
366  systs.push_back(&kNueExtrapSystSignalKin2017);
367  systs.push_back(&kNueExtrapSystBkg2017);
368  }
369  if(option.Contains("Numu") )
370  {
371  systs.push_back(&kMuEScaleSyst2017);
372  systs.push_back(&kRelMuEScaleSyst2017);
373  systs.push_back(&kDirectHadEScaleSyst2017);
374  systs.push_back(&kDirectRelHadEScaleSyst2017);
376  }
377  // Want kBeamAll
378  systs.push_back(&kBeamAllTransport);
379  // Fill names with the systs we've just added..
380  for (auto const & syst:systs)
381  names.push_back(
382  {syst->ShortName(), {"plusOne","minusOne","plusTwo","minusTwo"}});
383  }// if "All" && "Reduced"
384  if (option.Contains("All"))
385  {
386  for(auto c:{"FluxSyst", "NumuEnergySyst",
387  "NueExtrapSyst","AllSmall","XSecSyst"})
388  { //nue wants XSec to be last
389  auto systs = GetSystematics2018(c);
390  for (auto const & syst:systs)
391  names.push_back(
392  {syst->ShortName(),{"plusOne","minusOne","plusTwo","minusTwo"}});
393  }//end beam or xsec
394  }
395  else
396  {
397  auto systs = GetSystematics2018(option);
398  for (auto const & syst:systs)
399  names.push_back(
400  {syst->ShortName(), {"plusOne","minusOne","plusTwo","minusTwo"}});
401  }
402  return names;
403  }
404 
405 
406 }
BeamSyst * GetPPFXPrincipals(int PCIdx)
Definition: BeamSysts.cxx:63
std::vector< ShiDef > GetShifts2018(const TString option)
const WrongSignScale kWrongSignScale
const NOvARwgtSyst k2ndClassCurrs("2ndclasscurr","Second class currents", novarwgt::kSimpleSecondClassCurrentsSystKnob)
Second-class current syst. See documentation in NOvARwgt.
Definition: XSecSysts.h:71
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const NueAcceptSystSignalKin2018RHC kNueAcceptSystSignalKin2018RHC
std::vector< const ISyst * > getAllXsecSysts_2018()
Get master XSec syst list for 2018 analyses.
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const SummedSyst * getAna2018SummedSmallXsecSysts(const EAnaType2018 ana)
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const WrongSignEnergyBinRHC kWrongSignEnergyBinRHC
const NuTruthSystComponentScale kTauScaleSyst("NuTauScale","#nu_{#tau} Scale", kIsTau_NT &&!kIsNC_NT, 0.6, NuTruthSystComponentScale::kLinear)
100% uncertainty scale on taus
Definition: Systs.h:176
const NOvARwgtSyst kMAQEGenieReducedSyst2017(genie::rew::GSyst::AsString(genie::rew::EGSyst(rwgt::fReweightMaCCQE))++"_reduced", genie::rew::GSyst::AsString(genie::rew::EGSyst(rwgt::fReweightMaCCQE))++"_reduced", novarwgt::kMAQEGenieReducedSyst2017)
2017 &#39;reduced&#39; M_A^QE shift. See documentation in NOvARwgt.
Definition: XSecSysts.h:47
OStream cerr
Definition: OStream.cxx:7
void abs(TH1 *hist)
Float_t tmp
Definition: plot.C:36
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
const DummyAnaSyst kAnaRelativeCalibSyst("RelativeCalib","RelCalib")
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
BeamSyst * GetFluxPrincipals2018(int PCIdx)
Definition: BeamSysts.cxx:80
std::vector< const ISyst * > GetSystematics2018(const TString option)
const NOvARwgtSyst kDISvnCC1pi_2017("DISvnCC1pi_2017","DIS vnCC1pi", novarwgt::kDIS_CC_1pi_nu_n_SystKnob_2017)
Definition: DISSysts.h:59
const NOvARwgtSyst kRadCorrNuebar("radcorrnuebar","Radiative corrections for #bar{#nu}_{e}", novarwgt::kSimpleRadiativeCorrNuebarXsecSystKnob)
Radiative corrections syst (nuebars). See documentation in NOvARwgt.
Definition: XSecSysts.h:67
const MuEScaleSyst2017 kMuEScaleSyst2017(0.0074, 0.0012)
const DummyAnaSyst kAnaLightlevelSyst("Lightlevel","Lightlevel")
const NOvARwgtSyst kRPACCQEEnhSyst2017("RPAShapeenh2017","RPA shape: higher-Q^{2} enhancement (2017)", novarwgt::kRPACCQEEnhSyst2017)
Definition: RPASysts.h:11
const DummyAnaSyst kAnaCalibrationSyst("Calibration","AbsCalib")
const DummyAnaSyst kAnaCherenkovSyst("Cherenkov","Cherenkov")
const RelMuEScaleSyst2017 kRelMuEScaleSyst2017(0.0045, 10.5)
const NueAcceptSystBkg2018RHC kNueAcceptSystBkg2018RHC
RHC.
const BeamSyst kBeamAllTransport((FindCAFAnaDir()+"/data/beam/TABeamSyst_2018Dec11.root"),"totErr","beamTransportComb","Combined Beam Transport Systematics")
All Beam Transport systematics combined in quadratures.
Definition: BeamSysts.h:125
const NueAcceptSystSignalKin2018FHC kNueAcceptSystSignalKin2018FHC
const DirectHadEScaleSyst2017 kDirectHadEScaleSyst2017(0.05)
const NueExtrapSystSignalKin2017 kNueExtrapSystSignalKin2017
const NueAcceptSystBkg2018FHC kNueAcceptSystBkg2018FHC
FHC.
const NueExtrapSystBkg2017 kNueExtrapSystBkg2017
double sigma(TH1F *hist, double percentile)
const WrongSignScale100 kWrongSignScale100
GeniePCASyst * GetGeniePrincipals2018Small(int PCIdx)
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
const DirectRelHadEScaleSyst2017 kDirectRelHadEScaleSyst2017(0.05)
std::vector< const ISyst * > getAllXsecSysts_2018_RPAFix()
const WrongSignEnergyBin kWrongSignEnergyBin
const NOvARwgtSyst kMECq0ShapeSyst2017("MECq0Shape","MEC q_{0} shape", novarwgt::kMECq0ShapeSyst2017)
Definition: MECSysts.h:41
const MichelTaggingSyst2018 kMichelTaggingSyst2018
const NOvARwgtSyst kRPARESSyst2017("RPAShapeRES2017","RPA shape: resonance events", novarwgt::kRPARESSyst2017)
Definition: RPASysts.h:20
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
def entry(str)
Definition: HTMLTools.py:26
const NOvARwgtSyst kRadCorrNue("radcorrnue","Radiative corrections for #nu_{e}", novarwgt::kSimpleRadiativeCorrNueXsecSystKnob)
Radiative corrections syst (nues). See documentation in NOvARwgt.
Definition: XSecSysts.h:64
std::vector< SNameDef > GetShiftNames2018(const TString option)
const DummyAnaSyst kAnaCalibShapeSyst("CalibShape","CalibShape")