NusAna2020Systs.cxx
Go to the documentation of this file.
1 #include "CAFAna/Core/Registry.h"
2 #include "CAFAna/Core/Sample.h"
12 
13 #include "TFile.h"
14 #include "TKey.h"
15 
16 #include <cmath>
17 #include <cassert>
18 #include <iostream>
19 
20 using std::cout;
21 using std::endl;
22 using std::string;
23 using std::vector;
24 
25 using ana::covmx::Sample;
26 
27 namespace ana
28 {
29  const covmx::NormSyst kNusAna2020NormPOT( "POT", "POT uncertainty", 0.0055 );
30  const covmx::NormSyst kNusAna2020NormNDMass( "NDMass", "Near detector mass", 0.003 );
31  const covmx::NormSyst kNusAna2020NormFDMass( "FDMass", "Far detector mass", 0.003 );
32  const covmx::NormSyst kNusAna2020NormMissingLepton( "MissingLepton", "Missing lepton bug", 0.008 );
33  const covmx::NormSyst kNusAna2020NormNDNumuOverlay( "NDNumuOverlay", "ND numu selection overlay", 0.0021 );
34  const covmx::NormSyst kNusAna2020NormNDNCOverlay( "NDNCOverlay", "ND NC selection overlay", 0.01 );
35 
38 
52 
53  // Check whether a key syst is one-sided
54  bool IsOneSided(string systKey) {
55  for (string key : kOneSidedSysts) {
56  if (systKey.find(key) != string::npos) return true;
57  }
58  return false;
59  } // function IsOneSided
60 
61  // Check whether a syst should be uncorrelated between detectors
62  bool IsDetUncorrelated(string systKey) {
63  for (string key : kUncorrelatedSysts) {
64  if (systKey.find(key) != string::npos) return true;
65  }
66  return false;
67  } // function IsDetUncorrelated
68 
69  vector<string> GetKeySystNames() {
70  vector<string> ret;
71  for (string name : kDetectorSystNames) {
72  if (IsDetUncorrelated(name)) {
73  ret.push_back("key_neardet_"+name);
74  ret.push_back("key_fardet_"+name);
75  } else {
76  ret.push_back("key_"+name);
77  }
78  }
79  return ret;
80  } // function GetKeySystNames
81 
82  vector<const KeySyst*> GetKeySysts() {
83  vector<const KeySyst*> ret;
84  vector<string> keys = GetKeySystNames();
85  for (string key : keys) {
86  const KeySyst* keySyst = nullptr;
87  const ISyst* syst = Registry<ISyst>::ShortNameToPtr(key, true);
88  if (!syst) {
89  keySyst = new KeySyst(key, key);
90  } else {
91  keySyst = dynamic_cast<const KeySyst*>(syst);
92  }
93  ret.push_back(keySyst);
94  }
95  return ret;
96  } // function GetKeySysts
97 
98  // This function should be in CovMxSysts after the CommonAna migration
99  vector<const ISyst*> LoadSystsFromFile(string filePath, string dirName)
100  {
101 
102  DontAddDirectory guard;
103 
104  // Open up input directory
105  TFile* inFile = TFile::Open(filePath.c_str(), "read");
106  TDirectory* dir = inFile->GetDirectory(dirName.c_str());
107 
108  // Loop over everything in this directory and load it as a systematic
109  std::vector<const ISyst*> systs;
110  TIter next(dir->GetListOfKeys());
111  TKey* key;
112  while ((key = (TKey*)next())) {
113  // Check if this systematic has already been loaded
114  std::string shortName = ((TObjString*)((TDirectory*)key->ReadObj())->Get("name"))
115  ->GetString().Data();
116 
117  const ISyst* syst = Registry<ISyst>::ShortNameToPtr(shortName, true);
118  // If it's already been loaded, just return the existing version
119  if (syst) {
120  systs.push_back(syst);
121  } else {
122  // If it hasn't already been loaded, load it
123  cout << "Syst loading code works, but requires an overhaul." << endl;
124  TObjString* str = (TObjString*)((TDirectory*)key->ReadObj())->Get("type");
125  if (str->GetString() == "NueSyst")
126  systs.push_back(NueSyst::LoadFrom((TDirectory*)key->ReadObj()).release());
127  else if (str->GetString() == "NumuSyst")
128  systs.push_back(NumuSyst::LoadFrom((TDirectory*)key->ReadObj()).release());
129  else if (str->GetString() == "NCSyst")
130  systs.push_back(NCSyst::LoadFrom((TDirectory*)key->ReadObj()).release());
131  else
132  throw std::runtime_error("Syst " + str->GetString() + " not recognised!");
133  //systs.push_back(NuISyst::LoadFrom((TDirectory*)key->ReadObj()).release());
134  }
135  }
136 
137  inFile->Close();
138  return systs;
139 
140  } // function LoadSystsFromFile
141 
142  void AddNusAna2020NormSysts(vector<const ISyst*> &systs, Sample& s)
143  {
144  systs.push_back(&kNusAna2020NormPOT);
145  systs.push_back(&kNusAna2020NormMissingLepton);
146  if (s.detector == covmx::kNearDet) {
147  systs.push_back(&kNusAna2020NormNDMass);
148  if (s.selection == covmx::kCCNumu)
149  systs.push_back(&kNusAna2020NormNDNumuOverlay);
150  if (s.selection == covmx::kNC)
151  systs.push_back(&kNusAna2020NormNDNCOverlay);
152  }
153  if (s.detector == covmx::kFarDet)
154  systs.push_back(&kNusAna2020NormFDMass);
155  }
156 
157  void AddNusAna2020NormSysts(vector<const ISyst*> &systs)
158  {
159  systs.push_back(&kNusAna2020NormPOT);
160  systs.push_back(&kNusAna2020NormMissingLepton);
161  systs.push_back(&kNusAna2020NormNDMass);
162  systs.push_back(&kNusAna2020NormFDMass);
163  systs.push_back(&kNusAna2020NormNDNumuOverlay);
164  systs.push_back(&kNusAna2020NormNDNCOverlay);
165  }
166 
167  void AddNusAna2020FileSysts(vector<const ISyst*> &systs, Sample& s)
168  {
169  // Open systfile and get directory
170  string filePath = InputPath() + "/systs/isysts.root";
171  string dirName = "isysts_" + s.GetTag();
172 
173  for (const ISyst* syst : LoadSystsFromFile(filePath, dirName)) {
174  systs.push_back(syst);
175  }
176  }
177 
178  void AddNusAna2020XSecSysts(vector<const ISyst*> &systs)
179  {
180  for (const ISyst* syst : getAllXsecNuTruthSysts_2020()) {
181  // Skip MEC systs that are tuned to ND data
182  if (syst->ShortName() == "MECShape2020Nu" ||
183  syst->ShortName() == "MECShape2020AntiNu") continue;
184  systs.push_back(syst);
185  }
186  }
187 
188  void AddNusAna2020MECSysts ( vector<const ISyst*> &systs)
189  {
190  systs.push_back(&kNusAna2020MECGauss2DNorm_1);
191  systs.push_back(&kNusAna2020MECGauss2DMeanQ0_1);
192  systs.push_back(&kNusAna2020MECGauss2DMeanQ3_1);
193  systs.push_back(&kNusAna2020MECGauss2DSigmaQ0_1);
194  systs.push_back(&kNusAna2020MECGauss2DSigmaQ3_1);
195  systs.push_back(&kNusAna2020MECGauss2DCorr_1);
196  systs.push_back(&kNusAna2020MECGauss2DNorm_2);
197  systs.push_back(&kNusAna2020MECGauss2DMeanQ0_2);
198  systs.push_back(&kNusAna2020MECGauss2DMeanQ3_2);
199  systs.push_back(&kNusAna2020MECGauss2DSigmaQ0_2);
200  systs.push_back(&kNusAna2020MECGauss2DSigmaQ3_2);
201  systs.push_back(&kNusAna2020MECGauss2DCorr_2);
202  systs.push_back(&kNusAna2020MECBaseline);
203  }
204 
205  void AddNusAna2020BeamSysts (vector<const ISyst*> &systs)
206  {
207  // keeping the first 5 flux PCs for 2020
208  for (int i = 0; i < 5; ++i)
209  {
210  systs.push_back(GetFluxPrincipals2020(i));
211  }
212  }
213 
214  void AddNusAna2020KaonSyst(vector<const ISyst*> &systs)
215  {
216  systs.push_back(&kNusAna2020KaonSyst);
217  }
218 
219  void AddNusAna2020NeutronSyst(vector<const ISyst*> &systs)
220  {
221  systs.push_back(&kNeutronVisEScalePrimariesSyst2018);
222  }
223 
224  void AddNusAna2020TauSyst(vector<const ISyst*> &systs)
225  {
226  systs.push_back(&kNusAna2020TauSyst);
227  }
228 
229  void AddNusAna2020NumuSysts(vector<const ISyst*> &systs)
230  {
231  systs.push_back(&kCorrMuEScaleSyst2020);
232  systs.push_back(&kUnCorrMuCatMuESyst2020);
233  systs.push_back(&kUnCorrNDMuEScaleSyst2020);
234  systs.push_back(&kUnCorrFDMuEScaleSyst2020);
235  systs.push_back(&kPileupMuESyst2020);
236  }
237 
238  //////////////////////////////////////////////////////////////////////
239  vector<const ISyst*> getNusAna2020AllSysts(covmx::Sample& s)
240  {
241 
242  vector<const ISyst*> ret;
243  AddNusAna2020NormSysts(ret, s);
244  AddNusAna2020FileSysts(ret, s);
252 
253  return ret;
254  }
255 
256  //////////////////////////////////////////////////////////////////////
257  vector<const ISyst*> getNusAna2020AllSysts()
258  {
259 
260  vector<const ISyst*> ret;
261 
263  for (const KeySyst* syst : GetKeySysts())
264  ret.push_back(dynamic_cast<const ISyst*>(syst));
272 
273  return ret;
274  }
275 
276  // Set systematic aliases for a given sample
277  void SetSystAlias(Sample& sample) {
278 
279  // Mask off normalisation systs for samples that don't apply
280  if (sample.detector == covmx::kNearDet)
281  sample.SetSystAlias(&kNusAna2020NormFDMass, nullptr);
282  if (sample.detector == covmx::kFarDet)
283  sample.SetSystAlias(&kNusAna2020NormNDMass, nullptr);
284  if (sample.detector == covmx::kFarDet || sample.selection == covmx::kCCNumu)
285  sample.SetSystAlias(&kNusAna2020NormNDNCOverlay, nullptr);
286  if (sample.detector == covmx::kFarDet || sample.selection == covmx::kNC)
287  sample.SetSystAlias(&kNusAna2020NormNDNumuOverlay, nullptr);
288 
289  // Loop over systematic names
290  for (string systName : kDetectorSystNames) {
291  // Get systematic for this sample
292 
293  // If this is detector-uncorrelated, then we need the key for
294  // THIS detector to work, and the key for the OTHER detector
295  // to be a null pointer
296  if (IsDetUncorrelated(systName)) {
297  // ND syst
298  string keyND = "key_neardet_"+systName; // ND key
299  const ISyst* keySyst = Registry<ISyst>::ShortNameToPtr(keyND, true);
300  if (!keySyst) keySyst = new KeySyst(keyND, keyND);
301  const ISyst* childSyst;
302  if (sample.detector == covmx::kNearDet) {
303  childSyst = Registry<ISyst>::ShortNameToPtr(sample.GetTag()+"_"+systName);
304  } else {
305  childSyst = nullptr;
306  }
307  sample.SetSystAlias(keySyst, childSyst);
308 
309  std::string keyFD = "key_fardet_"+systName; // FD key
310  keySyst = Registry<ISyst>::ShortNameToPtr(keyFD, true);
311  if (!keySyst) keySyst = new KeySyst(keyFD, keyFD);
312  if (sample.detector == covmx::kFarDet) {
313  childSyst = Registry<ISyst>::ShortNameToPtr(sample.GetTag()+"_"+systName);
314  } else {
315  childSyst = nullptr;
316  }
317  sample.SetSystAlias(keySyst, childSyst);
318  } else {
319  // Get key systematic
320  const ISyst* childSyst = Registry<ISyst>::ShortNameToPtr(sample.GetTag()+"_"+systName);
321  std::string key = "key_"+systName;
322  const ISyst* keySyst = Registry<ISyst>::ShortNameToPtr(key, true);
323  if (!keySyst) keySyst = new KeySyst(key, key);
324  // Set alias from key to child
325  sample.SetSystAlias(keySyst, childSyst); // This should automatically set nullptr if the key is true
326  }
327  } // for syst name
328 
329  // Turn off numu systs for non-numu samples
330  if (sample.selection != covmx::kCCNumu) {
331  vector<const ISyst*> numuSysts;
332  AddNusAna2020NumuSysts(numuSysts);
333  for (const ISyst* keySyst : numuSysts) sample.SetSystAlias(keySyst, nullptr);
334  }
335 
336  } // function SetSystAlias
337 
338 }
const XML_Char * name
Definition: expat.h:151
const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DNorm_2(kGauss2DNorm_2Nux,"Norm_2Nux")
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void AddNusAna2020FileSysts(vector< const ISyst * > &systs, Sample &s)
void AddNusAna2020XSecSysts(vector< const ISyst * > &systs)
bool IsOneSided(string systKey)
const std::vector< std::string > kOneSidedSysts
vector< const ISyst * > getNusAna2020AllSysts(covmx::Sample &s)
void SetSystAlias(Sample &sample)
const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DMeanQ0_1(kGauss2DMeanQ0_1Nux,"MeanQ0_1Nux")
std::vector< const ISyst * > getAllXsecNuTruthSysts_2020()
Get master XSec syst list for 2020 analyses.
const NusAna2020TauSyst kNusAna2020TauSyst
const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DMeanQ3_2(kGauss2DMeanQ3_2Nux,"MeanQ3_2Nux")
const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DCorr_2(kGauss2DCorr_2Nux,"Corr_2Nux")
void AddNusAna2020NeutronSyst(vector< const ISyst * > &systs)
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
vector< string > GetKeySystNames()
void AddNusAna2020NormSysts(vector< const ISyst * > &systs, Sample &s)
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
const NusAna2020KaonSyst kNusAna2020KaonSyst
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const covmx::NormSyst kNusAna2020NormPOT("POT","POT uncertainty", 0.0055)
const covmx::NormSyst kNusAna2020NormMissingLepton("MissingLepton","Missing lepton bug", 0.008)
const UnCorrFDMuEScaleSyst2020 kUnCorrFDMuEScaleSyst2020(0.0015)
const MECDoubleGaussEnhSystNux kNusAna2020MECBaseline(kBaselineNux,"Baseline_2Nux")
ifstream inFile
Definition: AnaPlotMaker.h:34
const XML_Char * s
Definition: expat.h:262
const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DMeanQ3_1(kGauss2DMeanQ3_1Nux,"MeanQ3_1Nux")
static std::unique_ptr< NumuSyst > LoadFrom(TDirectory *dir)
LoadFrom implementation for NumuSyst.
Definition: CovMxSysts.cxx:77
Selection selection
Definition: Sample.h:95
const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DMeanQ0_2(kGauss2DMeanQ0_2Nux,"MeanQ0_2Nux")
BeamSyst * GetFluxPrincipals2020(int PCIdx)
Definition: BeamSysts.cxx:106
vector< const KeySyst * > GetKeySysts()
void AddNusAna2020TauSyst(vector< const ISyst * > &systs)
OStream cout
Definition: OStream.cxx:6
string InputPath()
Definition: Utilities.cxx:12
const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DSigmaQ0_2(kGauss2DSigmaQ0_2Nux,"SigmaQ0_2Nux")
static std::unique_ptr< NCSyst > LoadFrom(TDirectory *dir)
LoadFrom implementation for NCSyst.
Definition: CovMxSysts.cxx:67
void AddNusAna2020KaonSyst(vector< const ISyst * > &systs)
const covmx::NormSyst kNusAna2020NormNDNCOverlay("NDNCOverlay","ND NC selection overlay", 0.01)
const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DSigmaQ3_2(kGauss2DSigmaQ3_2Nux,"SigmaQ3_2Nux")
std::string dirName
Definition: PlotSpectra.h:47
string GetString(xmlDocPtr xml_doc, string node_path)
TDirectory * dir
Definition: macro.C:5
const PileupMuESyst2020 kPileupMuESyst2020(0.46, 1.3)
vector< const ISyst * > LoadSystsFromFile(string filePath, string dirName)
const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DCorr_1(kGauss2DCorr_1Nux,"Corr_1Nux")
const covmx::NormSyst kNusAna2020NormNDMass("NDMass","Near detector mass", 0.003)
const std::vector< std::string > kDetectorSystNames
void AddNusAna2020MECSysts(vector< const ISyst * > &systs)
const CorrMuEScaleSyst2020 kCorrMuEScaleSyst2020(0.0074, 0.0074, 0.0013)
const UnCorrNDMuEScaleSyst2020 kUnCorrNDMuEScaleSyst2020(0.0013)
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
static std::unique_ptr< NueSyst > LoadFrom(TDirectory *dir)
LoadFrom implementation for NueSyst.
Definition: CovMxSysts.cxx:171
const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DNorm_1(kGauss2DNorm_1Nux,"Norm_1Nux")
const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DSigmaQ0_1(kGauss2DSigmaQ0_1Nux,"SigmaQ0_1Nux")
const std::vector< std::string > kUncorrelatedSysts
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
const UnCorrMuCatMuESyst2020 kUnCorrMuCatMuESyst2020(0.0048)
void AddNusAna2020NumuSysts(vector< const ISyst * > &systs)
void next()
Definition: show_event.C:84
const covmx::NormSyst kNusAna2020NormFDMass("FDMass","Far detector mass", 0.003)
static const T * ShortNameToPtr(const std::string &s, bool allowFail=false)
Definition: Registry.cxx:60
bool IsDetUncorrelated(string systKey)
const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DSigmaQ3_1(kGauss2DSigmaQ3_1Nux,"SigmaQ3_1Nux")
void AddNusAna2020BeamSysts(vector< const ISyst * > &systs)
const covmx::NormSyst kNusAna2020NormNDNumuOverlay("NDNumuOverlay","ND numu selection overlay", 0.0021)
enum BeamMode string