PredictionInterpTemplates.cxx
Go to the documentation of this file.
2 
3 #include "CAFAna/Core/Loaders.h"
5 
7 
8 #include "TH2.h"
9 #include "TObjString.h"
10 #include "TVectorD.h"
11 
12 namespace ana
13 {
14  //----------------------------------------------------------------------
16  OscChannel chan) const
17  {
18  TH1D* h1 = s.ToTH1(s.POT(), kPOT);
19  TH2* h2 = ToTH2(s, s.POT(), kPOT,
21 
22  assert(h1->GetNbinsX() == 30);
23  assert(h2->GetNbinsX() == 3);
24  assert(h2->GetNbinsY() == 10);
25 
26  for(int i = 0; i < h1->GetNbinsX(); ++i){
27  const int ipid = i / 10;
28  const int iE = i % 10;
29 
30  // Double-check indexing
31  assert(AlmostEqual(h1->GetBinContent(i+1), h2->GetBinContent(ipid+1, iE+1)));
32 
33  const double w = fSyst->WeightFor(chan, fSigma,
34  h2->GetXaxis()->GetBinCenter(ipid+1),
35  h2->GetYaxis()->GetBinCenter(iE+1));
36 
37  h1->SetBinContent(i+1, w*h1->GetBinContent(i+1));
38  }
39 
40  delete h2;
41  return Spectrum(h1, s.GetLabels(), s.GetBinnings(), s.POT(), s.Livetime());
42  }
43 
44  //----------------------------------------------------------------------
46  Flavors::Flavors_t flav,
48  Sign::Sign_t sign) const
49  {
50  // Just predict something to get the right binning
52  ret.Clear();
53 
54  if(curr & Current::kCC){
55  if(flav & Flavors::kNuEToNuE)
56  ret += Weight(fPred->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, sign), kNuEToNuECC);
57 
58  if(flav & Flavors::kNuMuToNuMu)
59  ret += Weight(fPred->PredictComponent(calc, Flavors::kNuMuToNuMu, Current::kCC, sign), kNuMuToNuMuCC);
60 
61  if(flav & Flavors::kNuEToNuMu)
62  ret += Weight(fPred->PredictComponent(calc, Flavors::kNuEToNuMu, Current::kCC, sign), kNuEToNuMuCC);
63 
64  if(flav & Flavors::kNuEToNuTau)
65  ret += Weight(fPred->PredictComponent(calc, Flavors::kNuEToNuTau, Current::kCC, sign), kNuEToNuTauCC);
66 
67  if(flav & Flavors::kNuMuToNuE)
68  ret += Weight(fPred->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, sign), kNuMuToNuECC);
69 
70  if(flav & Flavors::kNuMuToNuTau)
71  ret += Weight(fPred->PredictComponent(calc, Flavors::kNuMuToNuTau, Current::kCC, sign), kNuMuToNuTauCC);
72  }
73  if(curr & Current::kNC){
74  assert(flav == Flavors::kAll); // Don't know how to calculate anything else
75  assert(sign == Sign::kBoth); // Why would you want to split NCs out by sign?
76 
77  ret += Weight(fPred->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth), kNCToNC);
78  }
79 
80  return ret;
81  }
82 
83  //----------------------------------------------------------------------
84  void PredictionTemplateShift::SaveTo(TDirectory* dir) const
85  {
86  TDirectory* tmp = gDirectory;
87 
88  dir->cd();
89  TObjString("PredictionTemplateShift").Write("type");
90 
91  fPred->SaveTo(dir->mkdir("pred"));
92 
93  TObjString sSyst = fSyst->ShortName().c_str();
94  sSyst.Write("syst");
95 
96  TVectorD vSigma(1);
97  vSigma[0] = fSigma;
98  vSigma.Write("sigma");
99 
100  tmp->cd();
101  }
102 
103  //----------------------------------------------------------------------
104  std::unique_ptr<PredictionTemplateShift> PredictionTemplateShift::LoadFrom(TDirectory* dir)
105  {
106  TObjString* tag = (TObjString*)dir->Get("type");
107  assert(tag);
108  assert(tag->GetString() == "PredictionTemplateShift");
109 
110  IPrediction* pred = ana::LoadFrom<IPrediction>(dir->GetDirectory("pred")).release();
111 
112  TObjString* tosSyst = (TObjString*)dir->Get("syst");
113  assert(tosSyst);
114  const std::string ssyst = tosSyst->GetString().Data();
115 
116  // Ugh, this is pretty dreadful. TODO think about how to serialize
117  // systematics better. ROOT6 to the rescue?
118  const NueSystFromHist* syst = 0;
121  if(s->ShortName() == ssyst){
122  if(syst != 0 && syst != s){
123  std::cout << "Error: multiple systs found matching '"
124  << ssyst << "'" << std::endl;
125  abort();
126  }
127  syst = s;
128  }
129  }
130  }
131  if(!syst){
132  std::cout << "Error: no syst found matching '" << ssyst << "'" << std::endl;
133  abort();
134  }
135 
136  TVectorD* vSigma = (TVectorD*)dir->Get("sigma");
137  assert(vSigma);
138  assert(vSigma->GetNrows() == 1);
139 
140  return std::unique_ptr<PredictionTemplateShift>(new PredictionTemplateShift(pred, syst, (*vSigma)[0]));
141  }
142 
143  //----------------------------------------------------------------------
145  PredictionInterpTemplates(std::vector<const NueSystFromHist*> systs,
147  IPrediction* nom)
148  {
149  fOscOrigin = osc->Copy();
150 
151  for(const NueSystFromHist* syst: systs){
152  ShiftedPreds sp;
153  sp.systName = syst->ShortName();
154 
155  // Genie reweight shifts aren't able to provide +/-3 sigma
156  if(syst->IsGenieReweight())
157  sp.shifts = {-2, -1, 0, +1, +2};
158  else
159  sp.shifts = {-3, -2, -1, 0, +1, +2, +3};
160 
161  for(int sigma: sp.shifts){
162  sp.preds.push_back(new PredictionTemplateShift(nom, syst, sigma));
163  }
164 
165  fPreds.emplace(syst, sp);
166  } // end for syst
167 
168  fPredNom = std::unique_ptr<IPrediction>(nom);
169  }
170 
171  //----------------------------------------------------------------------
173  {
174  }
175 
176 } // namespace
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
std::vector< SystGroupDef > systs
Definition: syst_header.h:384
std::string systName
What systematic we&#39;re interpolating.
General interface to oscillation calculators.
Definition: FwdDeclare.h:15
(&#39; appearance&#39;)
Definition: IPrediction.h:18
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:553
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
virtual const std::string & ShortName() const final
The name printed out to the screen.
Definition: ISyst.h:27
virtual void SaveTo(TDirectory *dir) const
(&#39;beam &#39;)
Definition: IPrediction.h:15
void SaveTo(TDirectory *dir) const override
virtual _IOscCalculator * Copy() const =0
Spectrum Weight(const Spectrum &s, OscChannel chan) const
void Clear()
Definition: Spectrum.cxx:878
Float_t tmp
Definition: plot.C:36
Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
const Binning kNueSAEnergyBinning
The energy part of the SA 2D binning.
Definition: Binning.h:100
const XML_Char * s
Definition: expat.h:262
Charged-current interactions.
Definition: IPrediction.h:39
std::vector< double > shifts
Shift values sampled.
double WeightFor(OscChannel chan, double sigma, double pid, double nueenergy) const
Binning PIDBinning() const
NuESystPID
Definition: NueSysts.h:12
bool AlmostEqual(double a, double b)
Definition: Utilities.cxx:732
double POT() const
Definition: Spectrum.h:263
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
Definition: Utilities.cxx:374
Oscillation probability calculators.
Definition: Calcs.h:5
def sign(x)
Definition: canMan.py:204
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
TH1F * h2
Definition: plot.C:45
PredictionTemplateShift(IPrediction *pred, const NueSystFromHist *syst, double sigma)
TH1F * h1
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir)
Definition: IPrediction.cxx:38
PredictionInterpTemplates(std::vector< const NueSystFromHist * > systs, osc::IOscCalculator *osc, IPrediction *nom)
TDirectory * dir
Definition: macro.C:5
std::vector< const NueSystFromHist * > getAllNueSecondAnaTemplateSysts(NuESystPID pid)
string syst
Definition: plotSysts.py:176
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
static std::unique_ptr< PredictionTemplateShift > LoadFrom(TDirectory *dir)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
std::vector< std::string > GetLabels() const
Definition: Spectrum.h:298
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
Float_t w
Definition: plot.C:20
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:266
std::vector< Binning > GetBinnings() const
Definition: Spectrum.h:299