NusAna2020Systs.cxx
Go to the documentation of this file.
1 #include "CAFAna/Core/Registry.h"
2 #include "CAFAna/Core/Sample.h"
11 
12 #include "TFile.h"
13 #include "TKey.h"
14 
15 #include <cmath>
16 #include <cassert>
17 #include <iostream>
18 
19 using std::cout;
20 using std::endl;
21 using std::string;
22 using std::vector;
23 
24 using ana::covmx::Sample;
25 
26 namespace ana
27 {
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  string InputPath() {
54  string ret;
55  if (kUseCVMFS) {
56  const char* cpath = getenv("NUSDATA_DIR");
57  if (!cpath)
58  throw std::runtime_error("nusdata directory not found! Exiting.");
59  ret = string(cpath) + "/nus20";
60  } else {
61  ret = kNus20Path;
62  }
63  if (kUseXRootD) ret = pnfs2xrootd(ret);
64  return ret;
65  }
66 
67  // Check whether a key syst is one-sided
68  bool IsOneSided(string systKey) {
69  for (string key : kOneSidedSysts) {
70  if (systKey.find(key) != string::npos) return true;
71  }
72  return false;
73  } // function IsOneSided
74 
75  // Check whether a syst should be uncorrelated between detectors
76  bool IsDetUncorrelated(string systKey) {
77  for (string key : kUncorrelatedSysts) {
78  if (systKey.find(key) != string::npos) return true;
79  }
80  return false;
81  } // function IsDetUncorrelated
82 
83  vector<string> GetKeySystNames() {
84  vector<string> ret;
85  for (string name : kDetectorSystNames) {
86  if (IsDetUncorrelated(name)) {
87  ret.push_back("key_neardet_"+name);
88  ret.push_back("key_fardet_"+name);
89  } else {
90  ret.push_back("key_"+name);
91  }
92  }
93  return ret;
94  } // function GetKeySystNames
95 
96  vector<const KeySyst*> GetKeySysts() {
97  vector<const KeySyst*> ret;
98  vector<string> keys = GetKeySystNames();
99  for (string key : keys) {
100  const KeySyst* keySyst = nullptr;
101  const ISyst* syst = Registry<ISyst>::ShortNameToPtr(key, true);
102  if (!syst) {
103  keySyst = new KeySyst(key, key);
104  } else {
105  keySyst = dynamic_cast<const KeySyst*>(syst);
106  }
107  ret.push_back(keySyst);
108  }
109  return ret;
110  } // function GetKeySysts
111 
112  // This function should be in CovMxSysts after the CommonAna migration
113  vector<const ISyst*> LoadSystsFromFile(string filePath, string dirName)
114  {
115 
116  DontAddDirectory guard;
117 
118  // Open up input directory
119  TFile* inFile = TFile::Open(filePath.c_str(), "read");
120  TDirectory* dir = inFile->GetDirectory(dirName.c_str());
121 
122  // Loop over everything in this directory and load it as a systematic
123  std::vector<const ISyst*> systs;
124  TIter next(dir->GetListOfKeys());
125  TKey* key;
126  while ((key = (TKey*)next())) {
127  // Check if this systematic has already been loaded
128  std::string shortName = ((TObjString*)((TDirectory*)key->ReadObj())->Get("name"))
129  ->GetString().Data();
130 
131  const ISyst* syst = Registry<ISyst>::ShortNameToPtr(shortName, true);
132  // If it's already been loaded, just return the existing version
133  if (syst) {
134  systs.push_back(syst);
135  } else {
136  // If it hasn't already been loaded, load it
137  cout << "Syst loading code works, but requires an overhaul." << endl;
138  TObjString* str = (TObjString*)((TDirectory*)key->ReadObj())->Get("type");
139  if (str->GetString() == "NueSyst")
140  systs.push_back(NueSyst::LoadFrom((TDirectory*)key->ReadObj()).release());
141  else if (str->GetString() == "NumuSyst")
142  systs.push_back(NumuSyst::LoadFrom((TDirectory*)key->ReadObj()).release());
143  else if (str->GetString() == "NCSyst")
144  systs.push_back(NCSyst::LoadFrom((TDirectory*)key->ReadObj()).release());
145  else
146  throw std::runtime_error("Syst " + str->GetString() + " not recognised!");
147  //systs.push_back(NuISyst::LoadFrom((TDirectory*)key->ReadObj()).release());
148  }
149  }
150 
151  inFile->Close();
152  return systs;
153 
154  } // function LoadSystsFromFile
155 
156  void AddNusAna2020NormSysts(vector<const ISyst*> &systs, Sample& s)
157  {
158  systs.push_back(&kNusAna2020NormPOT);
159  systs.push_back(&kNusAna2020NormMissingLepton);
160  if (s.detector == covmx::kNearDet) {
161  systs.push_back(&kNusAna2020NormNDMass);
162  if (s.selection == covmx::kCCNumu)
163  systs.push_back(&kNusAna2020NormNDNumuOverlay);
164  if (s.selection == covmx::kNC)
165  systs.push_back(&kNusAna2020NormNDNCOverlay);
166  }
167  if (s.detector == covmx::kFarDet)
168  systs.push_back(&kNusAna2020NormFDMass);
169  }
170 
171  void AddNusAna2020NormSysts(vector<const ISyst*> &systs)
172  {
173  systs.push_back(&kNusAna2020NormPOT);
174  systs.push_back(&kNusAna2020NormMissingLepton);
175  systs.push_back(&kNusAna2020NormNDMass);
176  systs.push_back(&kNusAna2020NormFDMass);
177  systs.push_back(&kNusAna2020NormNDNumuOverlay);
178  systs.push_back(&kNusAna2020NormNDNCOverlay);
179  }
180 
181  void AddNusAna2020FileSysts(vector<const ISyst*> &systs, Sample& s)
182  {
183  // Open systfile and get directory
184  string filePath = InputPath() + "/systs/isysts.root";
185  string dirName = "isysts_" + s.GetTag();
186 
187  for (const ISyst* syst : LoadSystsFromFile(filePath, dirName)) {
188  systs.push_back(syst);
189  }
190  }
191 
192  void AddNusAna2020XSecSysts(vector<const ISyst*> &systs)
193  {
194  for (const ISyst* syst : getAllXsecNuTruthSysts_2020()) {
195  // Skip MEC systs that are tuned to ND data
196  if (syst->ShortName() == "MECShape2020Nu" ||
197  syst->ShortName() == "MECShape2020AntiNu") continue;
198  systs.push_back(syst);
199  }
200  }
201 
202  void AddNusAna2020MECSysts ( vector<const ISyst*> &systs)
203  {
204  systs.push_back(&kNusAna2020MECGauss2DNorm_1);
205  systs.push_back(&kNusAna2020MECGauss2DMeanQ0_1);
206  systs.push_back(&kNusAna2020MECGauss2DMeanQ3_1);
207  systs.push_back(&kNusAna2020MECGauss2DSigmaQ0_1);
208  systs.push_back(&kNusAna2020MECGauss2DSigmaQ3_1);
209  systs.push_back(&kNusAna2020MECGauss2DCorr_1);
210  systs.push_back(&kNusAna2020MECGauss2DNorm_2);
211  systs.push_back(&kNusAna2020MECGauss2DMeanQ0_2);
212  systs.push_back(&kNusAna2020MECGauss2DMeanQ3_2);
213  systs.push_back(&kNusAna2020MECGauss2DSigmaQ0_2);
214  systs.push_back(&kNusAna2020MECGauss2DSigmaQ3_2);
215  systs.push_back(&kNusAna2020MECGauss2DCorr_2);
216  systs.push_back(&kNusAna2020MECBaseline);
217  }
218 
219  void AddNusAna2020BeamSysts (vector<const ISyst*> &systs)
220  {
221  // keeping the first 5 flux PCs for 2020
222  for (int i = 0; i < 5; ++i)
223  {
224  systs.push_back(GetFluxPrincipals2020(i));
225  }
226  }
227 
228  void AddNusAna2020KaonSyst(vector<const ISyst*> &systs)
229  {
230  systs.push_back(&kNusAna2020KaonSyst);
231  }
232 
233  void AddNusAna2020NeutronSyst(vector<const ISyst*> &systs)
234  {
235  systs.push_back(&kNeutronVisEScalePrimariesSyst2018);
236  }
237 
238  void AddNusAna2020TauSyst(vector<const ISyst*> &systs)
239  {
240  systs.push_back(&kNusAna2020TauSyst);
241  }
242 
243  void AddNusAna2020NumuSysts(vector<const ISyst*> &systs)
244  {
245  systs.push_back(&kCorrMuEScaleSyst2020);
246  systs.push_back(&kUnCorrMuCatMuESyst2020);
247  systs.push_back(&kUnCorrNDMuEScaleSyst2020);
248  systs.push_back(&kUnCorrFDMuEScaleSyst2020);
249  systs.push_back(&kPileupMuESyst2020);
250  }
251 
252  //////////////////////////////////////////////////////////////////////
253  vector<const ISyst*> getNusAna2020AllSysts(covmx::Sample& s)
254  {
255 
256  vector<const ISyst*> ret;
257  AddNusAna2020NormSysts(ret, s);
258  AddNusAna2020FileSysts(ret, s);
266 
267  return ret;
268  }
269 
270  //////////////////////////////////////////////////////////////////////
271  vector<const ISyst*> getNusAna2020AllSysts()
272  {
273 
274  vector<const ISyst*> ret;
275 
277  for (const KeySyst* syst : GetKeySysts())
278  ret.push_back(dynamic_cast<const ISyst*>(syst));
286 
287  return ret;
288  }
289 
290  // Set systematic aliases for a given sample
291  void SetSystAlias(Sample& sample) {
292 
293  // Mask off normalisation systs for samples that don't apply
294  if (sample.detector == covmx::kNearDet)
295  sample.SetSystAlias(&kNusAna2020NormFDMass, nullptr);
296  if (sample.detector == covmx::kFarDet)
297  sample.SetSystAlias(&kNusAna2020NormNDMass, nullptr);
298  if (sample.detector == covmx::kFarDet || sample.selection == covmx::kCCNumu)
299  sample.SetSystAlias(&kNusAna2020NormNDNCOverlay, nullptr);
300  if (sample.detector == covmx::kFarDet || sample.selection == covmx::kNC)
301  sample.SetSystAlias(&kNusAna2020NormNDNumuOverlay, nullptr);
302 
303  // Loop over systematic names
304  for (string systName : kDetectorSystNames) {
305  // Get systematic for this sample
306 
307  // If this is detector-uncorrelated, then we need the key for
308  // THIS detector to work, and the key for the OTHER detector
309  // to be a null pointer
310  if (IsDetUncorrelated(systName)) {
311  // ND syst
312  string keyND = "key_neardet_"+systName; // ND key
313  const ISyst* keySyst = Registry<ISyst>::ShortNameToPtr(keyND, true);
314  if (!keySyst) keySyst = new KeySyst(keyND, keyND);
315  const ISyst* childSyst;
316  if (sample.detector == covmx::kNearDet) {
317  childSyst = Registry<ISyst>::ShortNameToPtr(sample.GetTag()+"_"+systName);
318  } else {
319  childSyst = nullptr;
320  }
321  sample.SetSystAlias(keySyst, childSyst);
322 
323  std::string keyFD = "key_fardet_"+systName; // FD key
324  keySyst = Registry<ISyst>::ShortNameToPtr(keyFD, true);
325  if (!keySyst) keySyst = new KeySyst(keyFD, keyFD);
326  if (sample.detector == covmx::kFarDet) {
327  childSyst = Registry<ISyst>::ShortNameToPtr(sample.GetTag()+"_"+systName);
328  } else {
329  childSyst = nullptr;
330  }
331  sample.SetSystAlias(keySyst, childSyst);
332  } else {
333  // Get key systematic
334  const ISyst* childSyst = Registry<ISyst>::ShortNameToPtr(sample.GetTag()+"_"+systName);
335  std::string key = "key_"+systName;
336  const ISyst* keySyst = Registry<ISyst>::ShortNameToPtr(key, true);
337  if (!keySyst) keySyst = new KeySyst(key, key);
338  // Set alias from key to child
339  sample.SetSystAlias(keySyst, childSyst); // This should automatically set nullptr if the key is true
340  }
341  } // for syst name
342 
343  // Turn off numu systs for non-numu samples
344  if (sample.selection != covmx::kCCNumu) {
345  vector<const ISyst*> numuSysts;
346  AddNusAna2020NumuSysts(numuSysts);
347  for (const ISyst* keySyst : numuSysts) sample.SetSystAlias(keySyst, nullptr);
348  }
349 
350  } // function SetSystAlias
351 
352 }
const XML_Char * name
Definition: expat.h:151
const MECDoubleGaussEnhSystNux kNusAna2020MECGauss2DNorm_2(kGauss2DNorm_2Nux,"Norm_2Nux")
keys
Reco plots.
Definition: caf_analysis.py:46
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
const bool kUseCVMFS
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
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")
const std::string kNus20Path
static std::unique_ptr< NumuSyst > LoadFrom(TDirectory *dir)
LoadFrom implementation for NumuSyst.
Definition: CovMxSysts.cxx:77
std::string getenv(std::string const &name)
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()
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 bool kUseXRootD
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