3FlavorSystHelper.cxx
Go to the documentation of this file.
2 
5 #include "CAFAna/Systs/Systs.h"
10 #include "CAFAna/Systs/BeamSysts.h"
16 #include "3FlavorAna/Systs/GeniePCASyst.h" // GetGeniePrincipals2020
18 #include "3FlavorAna/Systs/MichelTaggingSyst.h"//MichelSystematic
21 
23 
24 #include <iostream>
25 
26 namespace ana
27 {
28  std::vector <const ISyst*> GetSystematics2020(const TString option)
29  {
30  std::vector <const ISyst*> ret;
31 
32  // same as in GetSystematics2018+CalibDrift
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("CalibDrift") || option.Contains("All"))
40  ret.push_back(&kAnaCalibDriftSyst);
41  if (option.Contains("LightLevel") || option.Contains("All"))
42  ret.push_back(&kAnaLightlevelSyst);
43  if (option.Contains("Light_Level_FD") || option.Contains("All"))
44  ret.push_back(&kAnaLightlevelFDSyst);
45  if (option.Contains("Light_Level_ND") || option.Contains("All"))
46  ret.push_back(&kAnaLightlevelNDSyst);
47  if (option.Contains("Cherenkov") || option.Contains("All"))
48  ret.push_back(&kAnaCherenkovSyst);
49 
50  if(option.Contains("XSecSyst") || option.Contains("All"))
51  {
52  std::vector<const ISyst*> geniesysts = getAllXsecSysts_2020();
53  // The tau systematic description is in the Nue executive summary,
54  // p. 25 of doc-26699-v9
55  geniesysts.push_back(&kTauScaleSyst);
56  int minsyst=0, maxsyst=geniesysts.size();
57  if(option.Contains("XSecSyst_Frac"))
58  {
59  minsyst = 0;
60  maxsyst = 3;
61  std::cout << "\n\nOnly looking at a small range: " << std::endl;
62  }
63 
64  if(option.Contains("XSecSyst_1of6")){ minsyst= 0; maxsyst=11;}
65  if(option.Contains("XSecSyst_2of6")){ minsyst=11; maxsyst=22;}
66  if(option.Contains("XSecSyst_3of6")){ minsyst=22; maxsyst=33;}
67  if(option.Contains("XSecSyst_4of6")){ minsyst=33; maxsyst=44;}
68  if(option.Contains("XSecSyst_5of6")){ minsyst=44; maxsyst=55;}
69  if(option.Contains("XSecSyst_6of6")){ minsyst=55; }
70  if(option.Contains("All")){ minsyst=0; }
71 
72  std::cout << minsyst << " -- " << maxsyst << "\n\n";
73 
74  for(int i = minsyst; i < maxsyst; ++i)
75  {
76  if (i >= (int)geniesysts.size()) break;
77  ret.push_back(geniesysts[i]);
78  }
79  }
80 
81  // Flux systematic PCs
82  if(option.Contains("FluxSyst") || option.Contains("All"))
83  {
84  // Selecting 8 principal components to include in the big PDF
85  // In the final fit we expect this number to be smaller ~4-5
86  // TODO: Add functionality to specify desired number of components at
87  // runtime.
88  for (int i = 0; i < 8; i++)
89  ret.push_back(GetFluxPrincipals2020(i));
90  }
91 
92  // Muon energy scale systematic
93  if (option.Contains("MuEnergySyst") || (option.Contains("All")))
94  {
95  ret.push_back(&kUnCorrNDMuEScaleSyst2020);
96  ret.push_back(&kUnCorrMuCatMuESyst2020);
97  ret.push_back(&kPileupMuESyst2020);
98  ret.push_back(&kCorrMuEScaleSyst2020);
99 
100  // should be only numu?
101  ret.push_back(&kUnCorrFDMuEScaleSyst2020);
102  }
103 
104  if ( option.Contains( "LeptonAngleSyst" ) || option.Contains("All") )
105  {
106  ret.push_back( &kLeptonAngleSystNDXZ2020 );
107  ret.push_back( &kLeptonAngleSystNDYZ2020 );
108  ret.push_back( &kLeptonAngleSystFDXZ2020 );
109  ret.push_back( &kLeptonAngleSystFDYZ2020 );
110  }
111 
112  // Summed small xsec systs
113  if (option.Contains("SummedSmall"))
114  {
116  ret.push_back(small);
117  }
118 
119  // Neutron systematic calculation
120  if(option.Contains("NeutronSyst")||(option.Contains("All")))
121  {
122  ret.push_back(&kNeutronVisEScalePrimariesSyst2018);
123  }
124 
125  // Nue signal acceptance systematic
126  if (option.Contains("AcceptanceSyst")|| (option.Contains("All") && option.Contains("nue")))
127  {
128  if(option.Contains("PtExtrap")){
129  ret.push_back(&kNueAcceptSystSignalKinPtExtrap2020FHC); // nue signal acceptance, FHC
130  ret.push_back(&kNueAcceptSystSignalKinPtExtrap2020RHC); // nue signal acceptance, RHC
131  }
132  else{
133  ret.push_back(&kNueAcceptSystSignalKin2020FHC); // nue signal acceptance, FHC
134  ret.push_back(&kNueAcceptSystSignalKin2020RHC); // nue signal acceptance, RHC
135  }
136  }
137 
138  // Add MichelTagging sytematics
139  if (option.Contains("MichelTaggingSyst")|| (option.Contains("All") && option.Contains("nue") && option.Contains("FHC")))
140  {
141  ret.push_back(&kMichelTaggingSyst2020);
142  }
143 
144  if(option.Contains("GENIEPCASyst") || option.Contains("All"))
145  {
146  // Selecting 5 principal components to include in the big PDF
147  // This matches the number we think we will include in the final fit
148  // TODO: Add functionality to specify desired number of components at
149  // runtime.
150  for (int i = 0; i < 12; i++)
151  ret.push_back(GetGeniePrincipals2020Small(i));
152  }
153 
154  /*
155  if (option.Contains("NuMuSmall"))
156  {
157  const ISyst* small = getAna2018SummedSmallXsecSysts(kNumuAna2018);
158  ret.push_back(small);
159  }
160 
161 
162  // Wrong sign scale systematic
163  if(option.Contains("WrongSignScale"))
164  {
165  ret.push_back(&kWrongSignScale);
166  ret.push_back(&kWrongSignScale100);
167  }
168 
169  // Wrong sign scale by energy bin
170  if(option.Contains("WrongSignEnergyBin"))
171  {
172  ret.push_back(&kWrongSignEnergyBin);
173  ret.push_back(&kWrongSignEnergyBinRHC);
174  }
175 
176  */
177  return ret;
178  } // std::vector <const ISyst*> GetSystematics2020(const TString option)
179 
180  std::vector <ShiDef> GetShifts2020(const TString option)
181  {
182  std::vector <ShiDef> ret;
183  TString topstr, lowstr;
184 
185  // Add nominal
186  if (option.Contains("All") || option.Contains("Nominal"))
187  {
188  topstr = "Nominal"; lowstr= "noShift";
189  ret.push_back({topstr, lowstr, kNoShift});
190  }
191 
192  // Two-sided non-file based systematics
193  if(option.Contains("FluxSyst") ||
194  option.Contains("GENIEPCASyst") ||
195  option.Contains("XSecSyst") ||
196  option.Contains("MuEnergySyst") ||
197  option.Contains("LeptonAngleSyst") ||
198  option.Contains("AllSmall") ||
199  option.Contains("NuMuSmall") ||
200  option.Contains("AcceptanceSyst") ||
201  option.Contains("NeutronSyst") ||
202  option.Contains("MichelTaggingSyst") ||
203  //option.Contains("WrongSignScale") ||
204  //option.Contains("WrongSignEnergyBin") ||
205  option.Contains("All"))
206  {
207  auto systs = GetSystematics2020(option);
208 
209  for (auto syst:systs){
210  for (auto sigma:{+2,+1,-1,-2})
211  {
212  TString sig_str = (sigma>0?"plus":"minus");
213  if(abs(sigma)==1) sig_str+="One";
214  if(abs(sigma)==2) sig_str+="Two";
215  ret.push_back({syst->ShortName(),sig_str,SystShifts(syst,sigma)});
216  std::cout << (TString) "Added " + syst->ShortName() + "_" + sig_str
217  << std::endl;
218  }//end sigmas
219  std::cout << std::endl;
220  }//end systs
221  }//if FluxSyst || GENIEPCASyst || XSecSyst || MuEnergySyst || AllSmall || AcceptanceSyst || NeutronSyst || MichelTaggingSyst
222 
223  // File-based systematics
224  if (option.Contains("CalibrationUp") || option.Contains("All"))
225  {
226  topstr = "Calibration"; lowstr ="plusOne";
227  ret.push_back({topstr, lowstr, kNoShift});
228  }
229  if(option.Contains("CalibrationDown") || option.Contains("All"))
230  {
231  topstr = "Calibration"; lowstr ="minusOne";
232  ret.push_back({topstr, lowstr, kNoShift});
233  }
234  if(option.Contains("RelativeCalibUp") || option.Contains("All"))
235  {
236  topstr = "RelativeCalib"; lowstr ="plusTwo";
237  ret.push_back({topstr, lowstr, kNoShift});
238  }
239  if(option.Contains("RelativeCalibDown") || option.Contains("All"))
240  {
241  topstr = "RelativeCalib"; lowstr ="minusTwo";
242  ret.push_back({topstr, lowstr, kNoShift});
243  }
244  if(option.Contains("CalibShape") || option.Contains("All"))
245  {
246  topstr = "CalibShape"; lowstr = "plusOne";
247  ret.push_back({topstr, lowstr, kNoShift});
248  }
249  if(option.Contains("CalibDrift") || option.Contains("All"))
250  {
251  topstr = "CalibDrift"; lowstr = "plusOne";
252  ret.push_back({topstr, lowstr, kNoShift});
253  }
254  if(option.Contains("Cherenkov") || option.Contains("All"))
255  {
256  topstr = "Cherenkov"; lowstr = "plusOne";
257  ret.push_back({topstr, lowstr, kNoShift});
258  }
259  if(option.Contains("LightLevelUp") || option.Contains("All"))
260  {
261  topstr = "Lightlevel"; lowstr = "plusOne";
262  ret.push_back({topstr, lowstr, kNoShift});
263  }
264  if(option.Contains("LightLevelDown") || option.Contains("All"))
265  {
266  topstr = "Lightlevel"; lowstr = "minusOne";
267  ret.push_back({topstr, lowstr, kNoShift});
268  }
269  if(option.Contains("Light_Level_FDUp") || option.Contains("All"))
270  {
271  topstr = "Light_Level_FD"; lowstr = "plusOne";
272  ret.push_back({topstr, lowstr, kNoShift});
273  }
274  if(option.Contains("Light_Level_FDDown") || option.Contains("All"))
275  {
276  topstr = "Light_Level_FD"; lowstr = "minusOne";
277  ret.push_back({topstr, lowstr, kNoShift});
278  }
279  if(option.Contains("Light_Level_NDUp") || option.Contains("All"))
280  {
281  topstr = "Light_Level_ND"; lowstr = "plusOne";
282  ret.push_back({topstr, lowstr, kNoShift});
283  }
284  if(option.Contains("Light_Level_NDDown") || option.Contains("All"))
285  {
286  topstr = "Light_Level_ND"; lowstr = "minusOne";
287  ret.push_back({topstr, lowstr, kNoShift});
288  }
289 
290  return ret;
291  } // std::vector <ShiDef> GetShifts2020(const TString option)
292 
293  std::vector<SNameDef> GetShiftNames2020(const TString option)
294  {
295  std::vector <SNameDef> names;
296 
297  if(option.Contains("All"))
298  names.push_back({"SumSmallXSec3Flavor2020",{"plusOne","minusOne","plusTwo","minusTwo"}});
299  if((option == "Calibration") || option.Contains("All"))
300  names.push_back({"Calibration", {"plusOne","minusOne"}});
301  if((option == "RelativeCalib") || option.Contains("All"))
302  names.push_back({"RelativeCalib",{"plusTwo","minusTwo"}});
303  if((option == "CalibShape") || option.Contains("All"))
304  names.push_back({"CalibShape", {"plusOne"}});
305  if((option == "CalibDrift") || option.Contains("All"))
306  names.push_back({"CalibDrift", {"plusOne"}});
307  if((option == "Cherenkov") || option.Contains("All"))
308  names.push_back({"Cherenkov", {"plusOne"}});
309  if((option == "Lightlevel") || option.Contains("All"))
310  names.push_back({"Lightlevel", {"plusOne","minusOne"}});
311  if((option == "Light_Level_FD") || option.Contains("All"))
312  names.push_back({"Light_Level_FD", {"plusOne","minusOne"}});
313  if((option == "Light_Level_ND") || option.Contains("All"))
314  names.push_back({"Light_Level_ND", {"plusOne","minusOne"}});
315 
316 
317  if (option.Contains("All"))
318  {
319  for(auto c:{"FluxSyst", "MuEnergySyst", "LeptonAngleSyst", "SummedSmall", "NeutronSyst", "AcceptanceSyst", "MichelTaggingSyst", "XSecSyst"})
320  { //nue wants XSec to be last
321  auto systs = GetSystematics2020(c);
322  for (auto const & syst:systs)
323  names.push_back(
324  {syst->ShortName(),{"plusOne","minusOne","plusTwo","minusTwo"}});
325  }//end beam or xsec
326  }
327  else
328  {
329  auto systs = GetSystematics2020(option);
330  for (auto const & syst:systs)
331  names.push_back(
332  {syst->ShortName(), {"plusOne","minusOne","plusTwo","minusTwo"}});
333  }
334 
335  return names;
336  } // std::vector<SNameDef> GetShiftNames2020(const TString option)
337 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const DummyAnaSyst kAnaLightlevelFDSyst("Light_Level_FD","Light_Level_FD")
const NuTruthSystComponentScale kTauScaleSyst("NuTauScale","#nu_{#tau} Scale", kIsTau_NT &&!kIsNC_NT, 0.6, NuTruthSystComponentScale::kLinear)
100% uncertainty scale on taus
Definition: Systs.h:176
void abs(TH1 *hist)
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
const NueAcceptSystSignalKin2020FHC kNueAcceptSystSignalKin2020FHC
const DummyAnaSyst kAnaRelativeCalibSyst("RelativeCalib","RelCalib")
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const LeptonAngleSyst kLeptonAngleSystNDYZ2020("LeptonAngleSystNDYZ2020","Lepton Angle Syst ND YZ", caf::kNEARDET, kAngleShiftYZ, 0.010)
Definition: AngleSysts.h:35
const LeptonAngleSyst kLeptonAngleSystNDXZ2020("LeptonAngleSystNDXZ2020","Lepton Angle Syst ND XZ", caf::kNEARDET, kAngleShiftXZ, 0.010)
Definition: AngleSysts.h:34
std::vector< ShiDef > GetShifts2020(const TString option)
const DummyAnaSyst kAnaLightlevelSyst("Lightlevel","Lightlevel")
const UnCorrFDMuEScaleSyst2020 kUnCorrFDMuEScaleSyst2020(0.0015)
std::vector< const ISyst * > getAllXsecSysts_2020()
const LeptonAngleSyst kLeptonAngleSystFDXZ2020("LeptonAngleSystFDXZ2020","Lepton Angle Syst FD XZ", caf::kFARDET, kAngleShiftXZ, 0.010)
Definition: AngleSysts.h:36
const DummyAnaSyst kAnaCalibrationSyst("Calibration","AbsCalib")
const DummyAnaSyst kAnaCherenkovSyst("Cherenkov","Cherenkov")
const NueAcceptSystSignalKin2020FHC kNueAcceptSystSignalKinPtExtrap2020FHC(ana::kExtrapPt,"accept_signalkin_pTextrap_FHC_2020")
BeamSyst * GetFluxPrincipals2020(int PCIdx)
Definition: BeamSysts.cxx:99
double sigma(TH1F *hist, double percentile)
const MichelTaggingSyst2020 kMichelTaggingSyst2020
const SystShifts kNoShift
Definition: SystShifts.cxx:21
OStream cout
Definition: OStream.cxx:6
std::vector< SNameDef > GetShiftNames2020(const TString option)
std::vector< const ISyst * > GetSystematics2020(const TString option)
const PileupMuESyst2020 kPileupMuESyst2020(0.46, 1.3)
const SummedSyst * get3FlavorAna2020SummedSmallXsecSysts(const EAnaType2020 ana)
const CorrMuEScaleSyst2020 kCorrMuEScaleSyst2020(0.0074, 0.0074, 0.0013)
const UnCorrNDMuEScaleSyst2020 kUnCorrNDMuEScaleSyst2020(0.0013)
const DummyAnaSyst kAnaLightlevelNDSyst("Light_Level_ND","Light_Level_ND")
const UnCorrMuCatMuESyst2020 kUnCorrMuCatMuESyst2020(0.0048)
const DummyAnaSyst kAnaCalibDriftSyst("CalibDrift","CalibDrift")
const LeptonAngleSyst kLeptonAngleSystFDYZ2020("LeptonAngleSystFDYZ2020","Lepton Angle Syst FD YZ", caf::kFARDET, kAngleShiftYZ, 0.010)
Definition: AngleSysts.h:37
const NueAcceptSystSignalKin2020RHC kNueAcceptSystSignalKinPtExtrap2020RHC(ana::kExtrapPt,"accept_signalkin_pTextrap_RHC_2020")
const DummyAnaSyst kAnaCalibShapeSyst("CalibShape","CalibShape")
GeniePCASyst * GetGeniePrincipals2020Small(int PCIdx)
const NueAcceptSystSignalKin2020RHC kNueAcceptSystSignalKin2020RHC