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