PredictionNumuFAHadE.cxx
Go to the documentation of this file.
2 
4 
5 #include "TDirectory.h"
6 #include "TH1.h"
7 #include "TMatrixD.h"
8 #include "TObjString.h"
9 #include "TVectorD.h"
10 
11 #include "CAFAna/Core/HistCache.h"
14 
15 #include "Utilities/func/MathUtil.h"
16 
17 namespace ana
18 {
19  /// \brief Cubic interpolation
20  ///
21  /// For each adjacent set of four points we determine coefficients for a
22  /// cubic which will be the curve between the center two. We constrain the
23  /// function to match the two center points and to have the right mean
24  /// gradient at them. This causes this patch to match smoothly with the next
25  /// one along. The resulting function is continuous and first and second
26  /// differentiable.
27  ///
28  /// Data is given at (x1-1,y0), (x1,y1), (x1+1,y2), (x1+2,y3) and the return
29  /// value is y(x). Valid when x1 < x < x1+1
30  double CubicInterp(double x, double x1,
31  double y0, double y1, double y2, double y3)
32  {
33  x -= x1;
34  assert(x >= 0 && x <= 1);
35 
36  // The matrices are simply the inverses of writing out the constraints
37  // expressed in the comments.
38 
39  const double v[4] = {y1, y2, (y2-y0)/2, (y3-y1)/2};
40  const double m[16] = { 2, -2, 1, 1,
41  -3, 3, -2, -1,
42  0, 0, 1, 0,
43  1, 0, 0, 0};
44  const TVectorD res = TMatrixD(4, 4, m) * TVectorD(4, v);
45 
46  return res(0)*util::cube(x) + res(1)*util::sqr(x) + res(2)*x + res(3);
47  }
48 
49  /// \brief Cubic interpolation when there is no y0 point.
50  ///
51  /// See comment on \ref CubicInterp for background.
52  ///
53  /// At the ends of the range we fit a quadratic instead with only one
54  /// constraint on the slope.
55  ///
56  /// Valid when x<x1+1
57  double CubicInterpLeftEdge(double x, double x1,
58  double y1, double y2, double y3)
59  {
60  x -= x1;
61  assert(x <= 1);
62 
63  const double v[3] = {y1, y2, (y3-y1)/2};
64  const double m[9] = { 1, -1, 1,
65  -2, 2, -1,
66  1, 0, 0};
67  const TVectorD res = TMatrixD(3, 3, m) * TVectorD(3, v);
68 
69  return res(0)*util::sqr(x) + res(1)*x + res(2);
70  }
71 
72  /// \brief Cubic interpolation when there is no y3 point
73  ///
74  /// See comments on \ref CubicInterp and \ref CubicInterpLeftEdge
75  ///
76  /// Valid when x>x1
77  double CubicInterpRightEdge(double x, double x1,
78  double y0, double y1, double y2)
79  {
80  x -= x1;
81  assert(x >= 0);
82 
83  const double v[3] = {y1, y2, (y2-y0)/2};
84  const double m[9] = {-1, 1, -1,
85  0, 0, 1,
86  1, 0, 0};
87  const TVectorD res = TMatrixD(3, 3, m) * TVectorD(3, v);
88 
89  return res(0)*util::sqr(x) + res(1)*x + res(2);
90  }
91 
92  //----------------------------------------------------------------------
95  const IPredictionGenerator& predGen,
97  {
98  // Doesn't matter if user passes that syst or not, we always handle it
99  // ourselves
100  std::vector<const ISyst*> cleanedSysts;
101  for(const ISyst* syst: systs){
102  if(syst->ShortName() != "AbsHadEScale") cleanedSysts.push_back(syst);
103  }
104 
105  for(int i = 0; i < 13; ++i){
106  const double hadEshift = -3+i*.5;
107 
108  SystShifts shift(&kHadEnergyScaleSyst, hadEshift);
109  fHadEs[i] = new PredictionInterp(cleanedSysts, osc, predGen, loaders, shift);
110  }
111  }
112 
113  //----------------------------------------------------------------------
115  {
116  for(int i = 0; i < 13; ++i) delete fHadEs[i];
117  }
118 
119  //----------------------------------------------------------------------
121  {
122  // No systs
123  return fHadEs[6]->Predict(calc);
124  }
125 
126  //----------------------------------------------------------------------
128  Flavors::Flavors_t flav,
130  Sign::Sign_t sign) const
131  {
132  // No systs
133  return fHadEs[6]->PredictComponent(calc, flav, curr, sign);
134  }
135 
136  //----------------------------------------------------------------------
138  const SystShifts& shift) const
139  {
140  return PredictComponentSyst(calc, shift,
143  Sign::kBoth);
144  }
145 
146  //----------------------------------------------------------------------
148  const SystShifts& shift,
149  Flavors::Flavors_t flav,
151  Sign::Sign_t sign) const
152  {
153  if(shift.GetShift(&kHadEnergyScaleSyst) == 0){
154  return fHadEs[6]->PredictComponentSyst(calc, shift, flav, curr, sign);
155  }
156 
157  SystShifts cleanedShift = shift;
158  cleanedShift.SetShift(&kHadEnergyScaleSyst, 0);
159 
160  // Convert hadE scale to our indexing, find the closest neighbours
161  const double x = shift.GetShift(&kHadEnergyScaleSyst)*2 + 6;
162 
163  int x1 = int(x);
164  if(x1 < 0) x1 = 0; // ensure x1 is valid
165  if(x1 > 11) x1 = 11; // ensure x1 and x2 are valid
166  const int x2 = x1+1;
167  const int x0 = x1-1; // no guarantees on x0
168  const int x3 = x1+2; // or x3
169 
170  const Spectrum s1 = fHadEs[x1]->PredictComponentSyst(calc, cleanedShift, flav, curr, sign);
171  const Spectrum s2 = fHadEs[x2]->PredictComponentSyst(calc, cleanedShift, flav, curr, sign);
172 
173  std::unique_ptr<TH1D> h0((x0 >= 0) ? fHadEs[x0]->PredictComponentSyst(calc, cleanedShift, flav, curr, sign).ToTH1(s1.POT()) : 0);
174  std::unique_ptr<TH1D> h1(s1.ToTH1(s1.POT()));
175  std::unique_ptr<TH1D> h2(s2.ToTH1(s1.POT()));
176  std::unique_ptr<TH1D> h3((x3 <= 12) ? fHadEs[x3]->PredictComponentSyst(calc, cleanedShift, flav, curr, sign).ToTH1(s1.POT()) : 0);
177 
178  for(int i = 0; i < h1->GetNbinsX()+2; ++i){
179  const double y0 = h0 ? h0->GetBinContent(i) : 0;
180  const double y1 = h1->GetBinContent(i);
181  const double y2 = h2->GetBinContent(i);
182  const double y3 = h3 ? h3->GetBinContent(i) : 0;
183 
184  double y = -1;
185  if(h0 && h3){
186  y = CubicInterp(x, x1, y0, y1, y2, y3);
187  }
188  if(!h0){
189  y = CubicInterpLeftEdge(x, x1, y1, y2, y3);
190  }
191  if(!h3){
192  y = CubicInterpRightEdge(x, x1, y0, y1, y2);
193  }
194  if(!h0 && !h3) assert(0 && "Not reached");
195 
196  // Result is written into h1
197  h1->SetBinContent(i, y);
198  } // end for i
199 
200  TH1D* todel = h2.release();
201  HistCache::Delete(todel);
202  if(h0){
203  todel = h0.release();
204  HistCache::Delete(todel);
205  }
206  if(h3){
207  todel = h3.release();
208  HistCache::Delete(todel);
209  }
210 
211  return Spectrum(std::move(h1), s1.GetLabels(), s1.GetBinnings(), s1.POT(), s1.Livetime());
212  }
213 
214  //----------------------------------------------------------------------
215  void PredictionNumuFAHadE::SaveTo(TDirectory* dir) const
216  {
217  TDirectory* tmp = gDirectory;
218 
219  dir->cd();
220  TObjString("PredictionNumuFAHadE").Write("type");
221 
222  for(int i = 0; i < 13; ++i){
223  fHadEs[i]->SaveTo(dir->mkdir(TString::Format("p%d", i)));
224  }
225 
226  tmp->cd();
227  }
228 
229  //----------------------------------------------------------------------
230  std::unique_ptr<PredictionNumuFAHadE> PredictionNumuFAHadE::LoadFrom(TDirectory* dir)
231  {
232  TObjString* tag = (TObjString*)dir->Get("type");
233  assert(tag);
234  assert(tag->GetString() == "PredictionNumuFAHadE");
235 
236  std::unique_ptr<PredictionNumuFAHadE> ret(new PredictionNumuFAHadE);
237 
238  for(int i = 0; i < 13; ++i){
239  ret->fHadEs[i] = ana::LoadFrom<IPrediction>(dir->GetDirectory(TString::Format("p%d", i))).release();
240  }
241 
242  return ret;
243  }
244 
245 } // namespace
Implements systematic errors by interpolation between shifted templates.
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
virtual Spectrum PredictComponent(osc::IOscCalculator *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
std::vector< SystGroupDef > systs
Definition: syst_header.h:384
General interface to oscillation calculators.
Definition: FwdDeclare.h:15
Float_t y1[n_points_granero]
Definition: compare.C:5
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 void SaveTo(TDirectory *dir) const
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
Float_t x1[n_points_granero]
Definition: compare.C:5
TH1F * h3
Definition: berger.C:36
Collection of SpectrumLoaders for many configurations.
Definition: Loaders.h:17
virtual void SaveTo(TDirectory *dir) const override
Float_t tmp
Definition: plot.C:36
T cube(T x)
More efficient cube function than pow(x,3)
Definition: MathUtil.h:26
virtual Spectrum Predict(osc::IOscCalculator *calc) const override
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
virtual Spectrum PredictComponentSyst(osc::IOscCalculator *calc, const SystShifts &syst, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const
virtual Spectrum PredictComponentSyst(osc::IOscCalculator *calc, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
Interactions of both types.
Definition: IPrediction.h:42
double CubicInterp(double x, double x1, double y0, double y1, double y2, double y3)
Cubic interpolation.
static void Delete(TH1D *&h)
Definition: HistCache.cxx:215
double CubicInterpRightEdge(double x, double x1, double y0, double y1, double y2)
Cubic interpolation when there is no y3 point.
T GetShift(const ISyst *syst) const
virtual Spectrum PredictSyst(osc::IOscCalculator *calc, const SystShifts &shift) const override
WARNING! Experts only. For numu FA final result contours.
double POT() const
Definition: Spectrum.h:263
Oscillation probability calculators.
Definition: Calcs.h:5
def sign(x)
Definition: canMan.py:204
TH1F * h2
Definition: plot.C:45
const HadEnergyScaleSyst kHadEnergyScaleSyst(.05)
Definition: EnergySysts.h:35
double CubicInterpLeftEdge(double x, double x1, double y1, double y2, double y3)
Cubic interpolation when there is no y0 point.
TH1F * h1
std::unique_ptr< IPrediction > LoadFrom< IPrediction >(TDirectory *dir)
Definition: IPrediction.cxx:38
TVectorT< double > TVectorD
Definition: Utilities.h:18
TDirectory * dir
Definition: macro.C:5
std::vector< Loaders * > loaders
Definition: syst_header.h:385
static std::unique_ptr< PredictionNumuFAHadE > LoadFrom(TDirectory *dir)
string syst
Definition: plotSysts.py:176
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
Given loaders and an MC shift, Generate() generates an IPrediction.
std::vector< std::string > GetLabels() const
Definition: Spectrum.h:298
TMatrixT< double > TMatrixD
Definition: Utilities.h:16
All neutrinos, any flavor.
Definition: IPrediction.h:26
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:266
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:70
std::vector< Binning > GetBinnings() const
Definition: Spectrum.h:299