ModularExtrapComponent.cxx
Go to the documentation of this file.
5 #include "CAFAna/Core/Ratio.h"
7 #include "CAFAna/Cuts/BeamCuts.h"
8 #include "CAFAna/Vars/Vars.h"
9 
10 #include "TH1.h"
11 #include "TH3.h"
12 #include "TDirectory.h"
13 #include "TObjString.h"
14 
15 #include <iostream>
16 
17 namespace ana
18 {
19  //---------------------------------------------------------------------------
20  // Definition to satisfy declaration in Core/LoadFromFile.h
21  template<> std::unique_ptr<ModularExtrapComponent>
23  {
24  TObjString* ptag = (TObjString*)dir->Get((name+"/type").c_str());
25  assert(ptag);
26 
27  const TString tag = ptag->GetString();
28  delete ptag;
29 
30  if(tag == "NoReweight")
31  return NoReweight::LoadFrom(dir, name);
32  if(tag == "TruthReweight")
33  return TruthReweight::LoadFrom(dir, name);
34  if(tag == "RecoReweight")
35  return RecoReweight::LoadFrom(dir, name);
36 
37  std::cerr << "Unknown Component Extrapolation type '"
38  << tag << "'" << std::endl;
39  abort();
40  }
41 
42  //---------------------------------------------------------------------------
43 
45  {
46  if (fQuiet)
47  return;
48 
49  if (fBins.empty())
50  return;
51 
52  std::cerr << "\nWARNING: There were attempted divisions by empty bins (for which a fallback to no reweighting was used)"
53  << "\n during extrapolation in the following bins (check your MC stats):"
54  << std::endl;
55 
56  for (const auto & tuple : this->fBins)
57  {
58  std::cerr << "\t Channel: " << std::get<0>(tuple)
59  << "\t\t Location: " << std::get<1>(tuple)
60  << "\t\t Bin Index: " << std::get<2>(tuple)
61  << std::endl;
62  }
63 
64  }
65 
66  //---------------------------------------------------------------------------
67 
69  {
70  if (!fEvaluated)
71  {
72  fCache = Eval();
73  fEvaluated = true;
74  }
75  return fCache;
76  }
77 
78  //---------------------------------------------------------------------------
79 
81  const IDecomp& decomp,
82  const DecompResult dr
83  ){
84  switch (dr){
85  case DecompResult::nue : return decomp.NueComponent();
86  case DecompResult::numu : return decomp.NumuComponent();
87  case DecompResult::nuebar : return decomp.AntiNueComponent();
88  case DecompResult::numubar : return decomp.AntiNumuComponent();
89  case DecompResult::NCtot : return decomp.NCTotalComponent();
90  case DecompResult::NCbar : return decomp.NCAntiComponent();
91  case DecompResult::NC : return decomp.NCComponent();
92  }
93  assert( 0 && "Bad DecompResult" );
94  }
95 
96  //---------------------------------------------------------------------------
97 
99  const DecompResult dr
100  ){
101  switch (dr){
102  case DecompResult::nue : return "Nue";
103  case DecompResult::numu : return "Numu";
104  case DecompResult::nuebar : return "NueBar";
105  case DecompResult::numubar : return "NumuBar";
106  case DecompResult::NC : return "NC";
107  case DecompResult::NCbar : return "NCbar";
108  case DecompResult::NCtot : return "NCtot";
109  }
110  assert( 0 && "Bad DecompResult" );
111  }
112 
113  //---------------------------------------------------------------------------
114 
116  const std::string str
117  ){
118  if (str=="Nue") return DecompResult::nue;
119  else if (str=="Numu") return DecompResult::numu;
120  else if (str=="NueBar") return DecompResult::nuebar;
121  else if (str=="NumuBar") return DecompResult::numubar;
122  else if (str=="NC") return DecompResult::NC;
123  else if (str=="NCbar") return DecompResult::NCbar;
124  else if (str=="NCtot") return DecompResult::NCtot;
125  else assert( 0 && "Bad DecompResult String" );
126  }
127 
128  //---------------------------------------------------------------------------
129 
131  const Spectrum& num,
132  const Spectrum& denom,
133  const std::string component,
134  const std::string location,
135  const Spectrum& mult
136  ){
137  const Eigen::ArrayXd numa = num.GetEigen(1e20);
138  const Eigen::ArrayXd denoma = denom.GetEigen(1e20);
139  const Eigen::ArrayXd multa = mult.GetEigen(1e20);
140 
141  Eigen::ArrayXd ratioa = numa;
142  ratioa.setZero();
143 
144  assert( (numa.size() == denoma.size()) && "Bin Mismatch" );
145  assert( (numa.size() == multa.size()) && "Bin Mismatch" );
146 
147  static DivByZeroCounter counter;
148  for (int bin = 0; bin < numa.size(); ++bin)
149  {
150  if ( denoma[bin] != 0 ){
151  ratioa[bin] = numa[bin] / denoma[bin];
152  } else {
153  ratioa[bin] = 1;
154  if ( numa[bin] != 0
155  || multa[bin] != 0 )
156  counter.fBins.insert(std::make_tuple(component, location, bin));
157  }
158  }
159 
160  return Ratio( std::move(ratioa), num.GetLabels(), num.GetBinnings() );
161  }
162 
163  bool ModularExtrapComponent::fQuiet = false;
164 
165  //---------------------------------------------------------------------------
166 
169  const HistAxis& axis,
170  const Cut& fdcut,
171  const SystShifts& shiftMC,
172  const Var& weight,
173  const Cut& flavors,
174  SpectrumLoaderBase& extraloaderswap,
175  SpectrumLoaderBase& extraloadertau )
176  : fRecoFD( loader, axis, fdcut && flavors, shiftMC, weight )
177  {
178  extraloaderswap.AddReweightableSpectrum(
179  fRecoFD, axis.GetVar1D(), kTrueE, fdcut && flavors, shiftMC, weight );
180  extraloadertau.AddReweightableSpectrum(
181  fRecoFD, axis.GetVar1D(), kTrueE, fdcut && flavors, shiftMC, weight );
182  }
183 
185  {return fRecoFD;}
186 
187  void NoReweight::SaveTo(TDirectory* dir, const std::string& name) const
188  {
189  TDirectory* tmp = gDirectory;
190 
191  dir = dir->mkdir(name.c_str()); // switch to subdir
192  dir->cd();
193 
194  TObjString("NoReweight").Write("type");
195  fRecoFD.SaveTo(dir, "RecoFD");
196 
197  dir->Write();
198  delete dir;
199 
200  tmp->cd();
201  }
202 
203  std::unique_ptr<NoReweight>
204  NoReweight::LoadFrom(TDirectory* dir, const std::string& name)
205  {
206  dir = dir->GetDirectory(name.c_str()); // switch to subdir
207  assert(dir);
208 
209  std::unique_ptr<OscillatableSpectrum> s = OscillatableSpectrum::LoadFrom(dir, "RecoFD");
210 
211  delete dir;
212 
213  return std::unique_ptr<NoReweight>(new NoReweight(*s));
214  }
215 
216  //---------------------------------------------------------------------------
217 
219  SpectrumLoaderBase& ndloader,
220  const HistAxis& axisFD,
221  const HistAxis& axisND,
222  const Cut& fdcut,
223  const SystShifts& shiftMC,
224  const Var& weight,
227  const Cut& ndcut,
228  const IDecomp& decomposition,
229  const DecompResult dr,
230  const Cut& ndflavor,
231  SpectrumLoaderBase& fdloader,
232  const Cut& fdflavors
233  )
234  : fRecoToTrueND(ndloader, axisND, ndcut && ndflavor, shiftMC, weight),
235  fTrueToRecoFD(fdloader, axisFD, fdcut && fdflavors, shiftMC, weight),
236  fDecomp(decomposition),
237  fOwnDecomp(false),
238  fDecompRes(dr),
239  fLabel(label),
240  fLatex(latex)
241  {}
242 
244  {
245  if(fOwnDecomp) delete &fDecomp;
246  }
247 
249  {
250 
251  //Copy to local variables because reweighting is in-place
252  OscillatableSpectrum recoToTrueND(fRecoToTrueND);
253  OscillatableSpectrum trueToRecoFD(fTrueToRecoFD);
254 
255  //Get ND data from Decomp
257 
258  //Compute Data/MC Ratio in reco energy bins to get divide-by-zero warnings
260  decompresult, fRecoToTrueND.Unoscillated(),
261  fLabel, "MC ND Reco",
263 
264  //ND Reco->True
265  recoToTrueND.ReweightToRecoSpectrum( decompresult );
266 
267  //Compute Data/MC Ratio in true energy bins
268  Ratio dataMCtrue = FormSmartRatio(
269  recoToTrueND.TrueEnergy(), fRecoToTrueND.TrueEnergy(),
270  fLabel, "MC ND Truth",
272 
273  // Multiply by Data/MC Ratio and add in FD truth information
275  * dataMCtrue );
276 
277  return trueToRecoFD;
278 
279  }
280 
281  void TruthReweight::SaveTo(TDirectory* dir, const std::string& name) const
282  {
283  TDirectory* tmp = gDirectory;
284 
285  dir = dir->mkdir(name.c_str()); // switch to subdir
286  dir->cd();
287 
288  TObjString("TruthReweight").Write("type");
289  fRecoToTrueND.SaveTo(dir, "RecoToTrueND");
290  fTrueToRecoFD.SaveTo(dir, "TrueToRecoFD");
291  fDecomp.SaveTo(dir, "Decomp");
292  TObjString(DRToString(fDecompRes).c_str()).Write("DecompRes");
293  TObjString(fLabel.c_str()).Write("Label");
294  TObjString(fLatex.c_str()).Write("Latex");
295 
296  dir->Write();
297  delete dir;
298 
299  tmp->cd();
300  }
301 
302  std::unique_ptr<TruthReweight>
303  TruthReweight::LoadFrom(TDirectory* dir, const std::string& name)
304  {
305  dir = dir->GetDirectory(name.c_str()); // switch to subdir
306  assert(dir);
307 
308  TObjString* dr = (TObjString*)dir->Get("DecompRes");
309  assert(dr);
310  TObjString* label = (TObjString*)dir->Get("Label");
311  TObjString* latex = (TObjString*)dir->Get("Latex");
312  assert(label);
313  assert(latex);
314 
316  *(OscillatableSpectrum::LoadFrom(dir, "RecoToTrueND")),
317  *(OscillatableSpectrum::LoadFrom(dir, "TrueToRecoFD")),
318  *(ana::LoadFrom<IDecomp>(dir, "Decomp").release()),
319  StringToDR(dr->GetString().Data()),
320  label->GetString().Data(),
321  latex->GetString().Data()
322  );
323  // We know we have the only copy because we just loaded it
324  ret->fOwnDecomp = true;
325 
326  delete dir;
327 
328  delete label;
329  delete latex;
330  delete dr;
331  return std::unique_ptr<TruthReweight>(ret);
332  }
333 
334  //---------------------------------------------------------------------------
335 
337  SpectrumLoaderBase& ndloader,
338  const HistAxis& axis,
339  const Cut& fdcut,
340  const SystShifts& shiftMC,
341  const Var& weight,
344  const Cut& ndcut,
345  const IDecomp& decomposition,
346  const DecompResult dr,
347  const Cut& ndflavor,
348  SpectrumLoaderBase& fdloader,
349  const Cut& fdflavors,
350  SpectrumLoaderBase& extrafdloaderswap,
351  SpectrumLoaderBase& extrafdloadertau
352  )
353  : fRecoND(ndloader, axis, ndcut && ndflavor, shiftMC, weight),
354  fTrueToRecoFD(fdloader, axis, fdcut && fdflavors, shiftMC, weight),
355  fDecomp(&decomposition),
356  fOwnDecomp(false),
357  fDecompRes(dr),
358  fLabel(label),
359  fLatex(latex)
360  {
361  extrafdloaderswap.AddReweightableSpectrum(
362  fTrueToRecoFD, axis.GetVar1D(), kTrueE, fdcut && fdflavors, shiftMC, weight);
363  extrafdloadertau.AddReweightableSpectrum(
364  fTrueToRecoFD, axis.GetVar1D(), kTrueE, fdcut && fdflavors, shiftMC, weight);
365  }
366 
368  {
369  if(fOwnDecomp) delete fDecomp;
370  }
371 
373  {
374 
375  //Copy to local variable because reweighting is in-place
377 
378  //Get ND data from Decomp
379  Spectrum decompresult(GetDecompResult(*fDecomp,fDecompRes));
380 
381  //Compute Data/MC Ratio
382  Ratio dataMC = FormSmartRatio(
383  decompresult, fRecoND,
384  fLabel, "MC ND Reco",
386 
387  // Multiply by Data/MC Ratio and add in FD truth information
389 
390  return result;
391 
392  }
393 
394  void RecoReweight::SaveTo(TDirectory* dir, const std::string& name) const
395  {
396  TDirectory* tmp = gDirectory;
397 
398  dir = dir->mkdir(name.c_str()); // switch to subdir
399  dir->cd();
400 
401  TObjString("RecoReweight").Write("type");
402  fRecoND.SaveTo(dir, "RecoND");
403  fTrueToRecoFD.SaveTo(dir, "TrueToRecoFD");
404  fDecomp->SaveTo(dir, "Decomp");
405  TObjString(DRToString(fDecompRes).c_str()).Write("DecompRes");
406  TObjString(fLabel.c_str()).Write("Label");
407  TObjString(fLatex.c_str()).Write("Latex");
408 
409  dir->Write();
410  delete dir;
411 
412  tmp->cd();
413  }
414 
415  std::unique_ptr<RecoReweight>
416  RecoReweight::LoadFrom(TDirectory* dir, const std::string& name)
417  {
418  dir = dir->GetDirectory(name.c_str()); // switch to subdir
419  assert(dir);
420 
421  assert(dir->GetDirectory("RecoND"));
422  assert(dir->GetDirectory("TrueToRecoFD"));
423  assert(dir->GetDirectory("Decomp"));
424  TObjString* dr = (TObjString*)dir->Get("DecompRes");
425  assert(dr);
426  TObjString* label = (TObjString*)dir->Get("Label");
427  TObjString* latex = (TObjString*)dir->Get("Latex");
428  assert(label);
429  assert(latex);
430 
432  *(Spectrum::LoadFrom(dir, "RecoND")),
433  *(OscillatableSpectrum::LoadFrom(dir, "TrueToRecoFD")),
434  *(ana::LoadFrom<IDecomp>(dir, "Decomp")).release(),
435  StringToDR(dr->GetString().Data()),
436  label->GetString().Data(),
437  latex->GetString().Data()
438  );
439  // We know we have the only copy because we just loaded it
440  ret->fOwnDecomp = true;
441  delete dr;
442  delete label;
443  delete latex;
444  delete dir;
445  return std::unique_ptr<RecoReweight>(ret);
446  }
447 
448 }
static std::unique_ptr< OscillatableSpectrum > LoadFrom(TDirectory *dir, const std::string &name)
const std::vector< Binning > & GetBinnings() const
Definition: Spectrum.h:256
const XML_Char * name
Definition: expat.h:151
void latex(double x, double y, std::string txt, double ang=0, int align=22, double size=0.042)
static std::unique_ptr< TruthReweight > LoadFrom(TDirectory *dir, const std::string &name)
std::unique_ptr< ModularExtrapComponent > LoadFrom< ModularExtrapComponent >(TDirectory *dir, const std::string &label)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
virtual Spectrum AntiNueComponent() const =0
#define location
OscillatableSpectrum fTrueToRecoFD
virtual Spectrum NumuComponent() const =0
const Var weight
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
virtual Spectrum NCAntiComponent() const =0
const Eigen::ArrayXd & GetEigen() const
NB these don&#39;t have POT scaling. For expert high performance ops only!
Definition: Spectrum.h:188
static std::string DRToString(DecompResult)
Helper function to turn a DecompResult into a string (for storage).
OscillatableSpectrum fRecoToTrueND
const OscillatableSpectrum & Return() const
Interface so that result of Eval() is called only once and cached.
OStream cerr
Definition: OStream.cxx:7
static DecompResult StringToDR(std::string)
Helper function to turn a string into a DecompResult (for loading).
virtual Spectrum NCComponent() const =0
const DecompResult fDecompRes
Float_t tmp
Definition: plot.C:36
void ReweightToRecoSpectrum(const Spectrum &target)
Recale bins so that Unweighted will return target.
OscillatableSpectrum Eval() const override
Core extrapolation math.
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const char * label
void SaveTo(TDirectory *dir, const std::string &name) const override
OscillatableSpectrum fTrueToRecoFD
const XML_Char * s
Definition: expat.h:262
NoReweight(SpectrumLoaderBase &loader, const HistAxis &axis, const Cut &fdcut, const SystShifts &shiftMC, const Var &weight, const Cut &flavors, SpectrumLoaderBase &extraloaderswap=kNullLoader, SpectrumLoaderBase &extraloadertau=kNullLoader)
void SaveTo(TDirectory *dir, const std::string &name) const override
const Var kTrueE([](const caf::SRProxy *sr){return(sr->mc.nnu==0)?0.:float(sr->mc.nu[0].E);})
Definition: Vars.cxx:85
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:546
std::vector< std::string > flavors
RecoReweight(SpectrumLoaderBase &ndloader, const HistAxis &axis, const Cut &fdcut, const SystShifts &shiftMC, const Var &weight, std::string label, std::string latex, const Cut &ndcut, const IDecomp &decomposition, const DecompResult dr, const Cut &ndflavor, SpectrumLoaderBase &fdloader, const Cut &fdflavors, SpectrumLoaderBase &extrafdloaderswap=kNullLoader, SpectrumLoaderBase &extrafdloadertau=kNullLoader)
Extrapolates component using truth-over-truth method.
DecompResult
Simple way to remember what to ask the decomposition for.
void SaveTo(TDirectory *dir, const std::string &name) const
Definition: Spectrum.cxx:517
static std::unique_ptr< NoReweight > LoadFrom(TDirectory *dir, const std::string &name)
loader
Definition: demo0.py:10
virtual Spectrum AntiNumuComponent() const =0
float bin[41]
Definition: plottest35.C:14
virtual OscillatableSpectrum Eval() const =0
Core extrapolation math.
Represent the ratio between two spectra.
Definition: Ratio.h:8
virtual void AddReweightableSpectrum(ReweightableSpectrum &spect, const Var &xvar, const Var &yvar, const Cut &cut, const SystShifts &shift, const Var &wei)
For use by the constructors of ReweightableSpectrum subclasses.
Base class for the various types of spectrum loader.
OscillatableSpectrum fRecoFD
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
static Spectrum GetDecompResult(const IDecomp &, DecompResult)
Helper function to pick out single Spectrum from a decomposition.
TDirectory * dir
Definition: macro.C:5
OscillatableSpectrum Eval() const override
Core extrapolation math.
int num
Definition: f2_nu.C:119
virtual Spectrum NCTotalComponent() const
Definition: IDecomp.h:18
void SaveTo(TDirectory *dir, const std::string &name) const override
std::set< std::tuple< std::string, std::string, int > > fBins
Standard interface to all decomposition techniques.
Definition: IDecomp.h:13
assert(nhit_max >=nhit_nbins)
static Ratio FormSmartRatio(const Spectrum &num, const Spectrum &denom, std::string component, std::string location, const Spectrum &mult)
Form Ratio, but be aware of zero division.
virtual void SaveTo(TDirectory *dir, const std::string &name) const =0
void ReweightToTrueSpectrum(const Spectrum &target)
Rescale bins so that WeightingVariable will return target.
TruthReweight(SpectrumLoaderBase &ndloader, const HistAxis &axisFD, const HistAxis &axisND, const Cut &fdcut, const SystShifts &shiftMC, const Var &weight, std::string label, std::string latex, const Cut &ndcut, const IDecomp &decomposition, const DecompResult dr, const Cut &ndflavor, SpectrumLoaderBase &fdloader, const Cut &fdflavors)
const std::vector< std::string > & GetLabels() const
Definition: Spectrum.h:255
std::unique_ptr< IDecomp > LoadFrom< IDecomp >(TDirectory *dir, const std::string &label)
Definition: IDecomp.cxx:53
Spectrum with true energy information, allowing it to be oscillated
const DecompResult fDecompRes
OscillatableSpectrum Eval() const override
Core extrapolation math.
void SaveTo(TDirectory *dir, const std::string &name) const
virtual Spectrum NueComponent() const =0
static std::unique_ptr< RecoReweight > LoadFrom(TDirectory *dir, const std::string &name)
T GetVar1D() const
Definition: HistAxis.cxx:43
Extrapolates using reco-over-reco method.
gm Write()