CovMxSysts.cxx
Go to the documentation of this file.
2 
4 
5 #include "TKey.h"
6 #include "TVectorD.h"
7 
8 #include <iostream>
9 
10 using std::cout;
11 using std::endl;
12 using std::make_pair;
13 using std::make_unique;
14 using std::map;
15 using std::pair;
16 using std::runtime_error;
17 using std::string;
18 using std::unique_ptr;
19 
20 namespace ana {
21 
22  // -------------------------------------------------------------------------------------
23  /// NumuSyst basic constructor
24  NumuSyst::NumuSyst(string name, string label, string sampleName,
25  map<int, pair<TH1*, TH1*> > shifts)
26  : NuISyst(name, label, sampleName, shifts)
27  {}
28 
29  // -------------------------------------------------------------------------------------
30  /// Implementation of standard ISyst::Shift function for NumuSyst
31  void NumuSyst::Shift(double sigma, caf::SRProxy* sr, double& weight) const {
32 
33  // Find the energy estimate for this slice
34  double calE = kNumuE2020(sr);
35 
36  // Find the signal or background oscillation channel
38  int sign;
39  if (sr->spill.isFHC) sign = 1;
40  else if (sr->spill.isRHC) sign = -1;
41  else throw runtime_error("Event is neither FHC nor RHC!");
42 
43  if (sr->mc.nnu == 0) {
44  chan = SystChannel::kBkg;
45  }
46  else if (sr->mc.nu[0].iscc && sr->mc.nu[0].pdg == sign*14 && sr->mc.nu[0].pdgorig == sign*14) {
47  chan = SystChannel::kSig;
48  }
49  else if (!sr->mc.nu[0].iscc || sr->mc.nu[0].pdg != sign*14 || sr->mc.nu[0].pdgorig != sign*14) {
50  chan = SystChannel::kBkg;
51  }
52  if (chan == SystChannel::kUnknown) {
53  cout << "Error: Unknown oscillation channel." << endl;
54  abort();
55  }
56 
57  weight *= WeightFor(chan, sigma, calE);
58 
59  } // function NumuSyst::Shift
60 
61  // -------------------------------------------------------------------------------------
62  /// SaveTo implementation for NumuSyst
63  void NumuSyst::SaveTo(TDirectory* dir) const {
64 
65  TDirectory* tmp = gDirectory;
66  dir->cd();
67  TObjString("NumuSyst").Write("type");
68  tmp->cd();
69  NuISyst::SaveTo(dir);
70 
71  } // function NumuSyst::SaveTo
72 
73  // -------------------------------------------------------------------------------------
74  /// LoadFrom implementation for NumuSyst
75  unique_ptr<NumuSyst> NumuSyst::LoadFrom(TDirectory* dir) {
76 
77  DontAddDirectory guard;
78 
79  // Check this is a NumuSyst
80  TObjString* tag = (TObjString*)dir->Get("type");
81  assert(tag);
82  string sTag = tag->GetString().Data();
83  assert(sTag == "NumuSyst");
84 
85  // Load names
86  TObjString* name = (TObjString*)dir->Get("name");
87  assert(name);
88  string sName = name->GetString().Data();
89  TObjString* label = (TObjString*)dir->Get("label");
90  assert(label);
91  string sLabel = label->GetString().Data();
92  TObjString* sampleName = (TObjString*)dir->Get("samplename");
93  assert(sampleName);
94  string sSampleName = sampleName->GetString().Data();
95 
96  map<int, pair<TH1*, TH1*>> shifts;
97 
98  // Load all shifts
99  TIter next(dir->GetDirectory("Sigmas")->GetListOfKeys());
100  TKey* key;
101  while ((key = (TKey*)next())) {
102  TDirectory* sigmaDir = (TDirectory*)key->ReadObj();
103  const TVectorD vSigma = *(TVectorD*)sigmaDir->Get("sigma");
104  TH1* hSignal = (TH1*)sigmaDir->Get("sig");
105  TH1* hBackground = (TH1*)sigmaDir->Get("bkg");
106  shifts[vSigma[0]] = make_pair(hSignal, hBackground);
107  }
108 
109  return make_unique<NumuSyst>(sName, sLabel, sSampleName, shifts);
110 
111  } // function NumuSyst::SaveTo
112 
113  // -------------------------------------------------------------------------------------
114  /// Implementation of standard ISyst::Shift function for NueSyst
115  NueSyst::NueSyst(string name, string label, string sampleName,
116  map<int, pair<TH1*, TH1*> > shifts)
117  : NuISyst(name, label, sampleName, shifts)
118  {}
119 
120  // -------------------------------------------------------------------------------------
121  /// Implementation of standard ISyst::Shift function for NueSyst
122  void NueSyst::Shift(double sigma, caf::SRProxy* sr, double& weight) const {
123 
124  if (sr->hdr.det == caf::kNEARDET && !kNue2020ND(sr)) return;
125  if (sr->hdr.det == caf::kFARDET && !kNue2020FDAllSamples_ML(sr)) return;
126 
127  // Find the energy estimate for this slice
128  double calE = kNue2020AnaBin(sr);
129 
130  // Find the signal or background oscillation channel
132  int sign;
133  if (sr->spill.isFHC) sign = 1;
134  else if (sr->spill.isRHC) sign = -1;
135  else throw runtime_error("Event is neither FHC nor RHC!");
136 
137  if (sr->mc.nnu == 0) {
138  chan = SystChannel::kBkg;
139  }
140  else if (sr->mc.nu[0].iscc && sr->mc.nu[0].pdg == sign*12 && sr->mc.nu[0].pdgorig == sign*14) {
141  chan = SystChannel::kSig;
142  }
143  else if (!sr->mc.nu[0].iscc || sr->mc.nu[0].pdg != sign*12 || sr->mc.nu[0].pdgorig != sign*14) {
144  chan = SystChannel::kBkg;
145  }
146  if (chan == SystChannel::kUnknown) {
147  cout << "Error: Unknown oscillation channel." << endl;
148  abort();
149  }
150 
151  weight *= WeightFor(chan, sigma, calE);
152 
153  } // function NueSyst::Shift
154 
155  // -------------------------------------------------------------------------------------
156  /// SaveTo implementation for NueSyst
157  void NueSyst::SaveTo(TDirectory* dir) const {
158 
159  TDirectory* tmp = gDirectory;
160  dir->cd();
161  TObjString("NueSyst").Write("type");
162  tmp->cd();
163  NuISyst::SaveTo(dir);
164 
165  } // function NueSyst::SaveTo
166 
167  // -------------------------------------------------------------------------------------
168  /// LoadFrom implementation for NueSyst
169  unique_ptr<NueSyst> NueSyst::LoadFrom(TDirectory* dir) {
170 
171  DontAddDirectory guard;
172 
173  // Check this is a NueSyst
174  TObjString* tag = (TObjString*)dir->Get("type");
175  assert(tag);
176  string sTag = tag->GetString().Data();
177  assert(sTag == "NueSyst");
178 
179  // Load names
180  TObjString* name = (TObjString*)dir->Get("name");
181  assert(name);
182  string sName = name->GetString().Data();
183  TObjString* label = (TObjString*)dir->Get("label");
184  assert(label);
185  string sLabel = label->GetString().Data();
186  TObjString* sampleName = (TObjString*)dir->Get("samplename");
187  assert(sampleName);
188  string sSampleName = sampleName->GetString().Data();
189 
190  map<int, pair<TH1*, TH1*>> shifts;
191 
192  // Load all shifts
193  TIter next(dir->GetDirectory("Sigmas")->GetListOfKeys());
194  TKey* key;
195  while ((key = (TKey*)next())) {
196  TDirectory* sigmaDir = (TDirectory*)key->ReadObj();
197  const TVectorD vSigma = *(TVectorD*)sigmaDir->Get("sigma");
198  TH1* hSignal = (TH1*)sigmaDir->Get("sig");
199  TH1* hBackground = (TH1*)sigmaDir->Get("bkg");
200  shifts[vSigma[0]] = make_pair(hSignal, hBackground);
201  }
202 
203  return make_unique<NueSyst>(sName, sLabel, sSampleName, shifts);
204 
205  } // function NumuSyst::SaveTo
206 
207 } // namespace ana
caf::Proxy< caf::SRSpill > spill
Definition: SRProxy.h:2125
Near Detector underground.
Definition: SREnums.h:10
const XML_Char * name
Definition: expat.h:151
Far Detector at Ash River.
Definition: SREnums.h:11
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const Var weight
caf::Proxy< caf::SRHeader > hdr
Definition: SRProxy.h:2119
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2108
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:597
virtual double WeightFor(SystChannel channel, double sigma, double calE) const
Get weight for a given energy, channel and sigma shift.
Definition: CovMxSysts.cxx:80
NueSyst(std::string name, std::string label, std::string sample_name="", std::map< int, std::pair< TH1 *, TH1 * > > shifts={})
Implementation of standard ISyst::Shift function for NueSyst.
Definition: CovMxSysts.cxx:115
caf::Proxy< short int > nnu
Definition: SRProxy.h:596
Float_t tmp
Definition: plot.C:36
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
const Cut kNue2020ND
Definition: NueCuts2020.h:178
const char * label
SystChannel
Definition: CovMxSysts.h:12
const Var kNumuE2020
Definition: NumuEFxs.h:217
caf::Proxy< bool > isFHC
Definition: SRProxy.h:1354
virtual void SaveTo(TDirectory *dir) const override
SaveTo implementation for NumuSyst.
Definition: CovMxSysts.cxx:63
static std::unique_ptr< NumuSyst > LoadFrom(TDirectory *dir)
LoadFrom implementation for NumuSyst.
Definition: CovMxSysts.cxx:75
virtual void SaveTo(TDirectory *dir) const override
SaveTo implementation for NueSyst.
Definition: CovMxSysts.cxx:157
caf::StandardRecord * sr
virtual void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Implementation of standard ISyst::Shift function for NueSyst.
Definition: CovMxSysts.cxx:122
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2120
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TDirectory * dir
Definition: macro.C:5
const Var kNue2020AnaBin([](const caf::SRProxy *sr){int selBin=kNue2020SelectionBin(sr);float nuE=kNueEnergy2020(sr);int nuEBin=nuE/0.5;assert(nuEBin<=8 &&"An event with nuE > 4.5 should never happen");int anaBin=9 *selBin+nuEBin;return anaBin;})
Use this Analysis Binning for Ana2020, official Binning.
Definition: NueCuts2020.h:191
assert(nhit_max >=nhit_nbins)
static std::unique_ptr< NueSyst > LoadFrom(TDirectory *dir)
LoadFrom implementation for NueSyst.
Definition: CovMxSysts.cxx:169
virtual void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Implementation of standard ISyst::Shift function for NumuSyst.
Definition: CovMxSysts.cxx:31
caf::Proxy< bool > isRHC
Definition: SRProxy.h:1355
const Cut kNue2020FDAllSamples_ML
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
virtual void SaveTo(TDirectory *dir) const
SaveTo implementation for NuISyst.
Definition: CovMxSysts.cxx:51
NumuSyst(std::string name, std::string label, std::string sample_name="", std::map< int, std::pair< TH1 *, TH1 * > > shifts={})
NumuSyst basic constructor.
Definition: CovMxSysts.cxx:24
void next()
Definition: show_event.C:84
def sign(x)
Definition: canMan.py:197
caf::Proxy< caf::Det_t > det
Definition: SRProxy.h:231