demoSysts.C
Go to the documentation of this file.
1 
2 #include "TCanvas.h"
3 #include "TLegend.h"
4 #include "TPad.h"
5 #include "TH1D.h"
6 
9 #include "CAFAna/Core/ISyst.h"
10 #include "CAFAna/Core/Ratio.h"
11 #include "CAFAna/Core/Spectrum.h"
13 #include "CAFAna/Cuts/SpillCuts.h"
16 #include "CAFAna/Vars/HistAxes.h"
18 #include "OscLib/OscCalc.h"
20 #include "Utilities/func/MathUtil.h"
21 
22 using namespace ana;
23 
24 namespace demo
25 {
26  // Systs are classes derived from ana::ISyst.
27  // Instances of each Syst class must be unique:
28  // if you instantiate more than one, you'll get a complaint
29  // from the SystRegistry that they're stored in.
30  // (More on that later.)
31 
32  // Let's take a simple first example:
33  // imagine you thought that there was a beam-related uncertainty
34  // that depends on Enu: +1sigma = fixed 10% for Enu < 5 GeV,
35  // and 10% + a function that blows up smoothly from 10% to 50%
36  // from 5 to 15 GeV:
37  // y(E) = 0.1 + 0.4*Heaviside(E-5)*(1-exp(-(E-5)**2/10))
38  class DemoSyst1 : public ISyst
39  {
40  public:
42  // they have to supply a "short" name (used for lookup internally)
43  // and a "Latex" name (which is used on plots).
44  // The "short" name MUST be unique or you'll get runtime errors...
45  : ISyst("DemoSyst1", "Demo syst ##1")
46  {}
47 
48  // the magic happens in TruthShift() (when truth-only) or Shift() (when reco involved).
49  // here, we're just reweighting events according to the function
50  // described above, so we just adjust the weight
51  void TruthShift(double sigma, caf::SRNeutrinoProxy* sr, double& weight) const override
52  {
53 
54  // first calculate the function.
55  double additionalWgt = 0.1;
56  double Enu = sr->E;
57  if (Enu > 5)
58  additionalWgt += 0.4 * (1 - exp(-util::sqr(Enu-5)/10.));
59 
60  // we just constructed the +1sigma function.
61  // now compute the actual shift requested.
62  additionalWgt *= sigma;
63 
64  // obviously the nominal (0 sigma) means leave the event alone,
65  // so that's weight = 1.
66  // increasing by 10% is weight = 1.1, etc.
67  // so we add the change to 1 to get the weight to use.
68  additionalWgt += 1;
69 
70  // we don't want events with negative weights...
71  if (additionalWgt < 0)
72  additionalWgt = 0;
73 
74  // remember to multiply it into the weight that was already there...
75  weight *= additionalWgt;
76  }
77  };
78 
79  // In CAFAna/Syst, you'll find that the instances of the ISyst classes
80  // in .h files are all declared 'extern', which tells the compiler
81  // not to instantiate them based on the header declaration.
82  // They get instantiated by the corresponding lines in the .cxx files.
83  // This way, even if a .h is #included in multiple places, there
84  // will actually only ever be a single instance of the Syst object
85  // in the library that's ultimately made -- the one coming from the .cxx.
86  // This demo is standalone, so the implicit second copy that gets made
87  // by ACLiC (don't worry about it) and you see warnings about in the output
88  // won't hurt anything.
90 
91  // ---------------------------------------
92 
93  // Now imagine a systematic that needs to change the values
94  // of some part of the reconstructed event
95  // (for instance, an energy scale shift.)
96  // That (usually) can't be implemented just by reweighting.
97  // In this example we'll make one that adds a fixed amount of visible energy
98  // for each neutron in the event.
99  class DemoSyst2 : public ISyst
100  {
101  public:
103  : ISyst("DemoSyst2", "Demo syst ##2")
104  {}
105 
106  // we'll be modifying the SRProxy this time.
107  // (that's why it's passed non-const.)
108  void Shift(double sigma, caf::SRProxy* sr, double& weight) const override
109  {
110  // if no truth info, we can't do anything
111  if (sr->mc.nnu != 1)
112  return;
113 
114  unsigned int nNeutron = 0;
115  for (const auto & particle : sr->mc.nu[0].prim)
116  {
117  if (particle.pdg != 2112)
118  continue;
119 
120  // only increase if they were less than 50 MeV to begin with
121  if (particle.visEinslcBirks + particle.daughterVisEinslcBirks >= 0.050)
122  continue;
123 
124  nNeutron++;
125  }
126 
127  // now adjust the numu energy.
128  // we're going to add 25 MeV of visible energy for each neutron
129  // for the +1sigma shift (and subtract it for -1sigma, etc.)
130  double addedE = sigma * 0.025 * nNeutron;
131  if (sr->energy.numu.hadcalE + addedE < 0)
132  addedE = -sr->energy.numu.hadcalE; // don't let the shift make a negative energy deposit
133 
134  sr->energy.numu.hadcalE += addedE; // this goes into the numu energy estimator function
135 
136  }
137  };
138 
139  // now instantiate
141 
142  // ---------------------------------------
143  // makes drawing easier later
146  const ISyst* syst)
147  {
148  SystShifts shifts;
149 
150  // get the spectra from the predictions
151  auto spec_nom = pred->PredictUnoscillated();
152  shifts.SetShift(syst, 1.0);
153  auto spec_shift_up = pred->PredictSyst(calc, shifts);
154  Ratio rat_shift_up(spec_shift_up, spec_nom);
155  shifts.SetShift(syst, -1.0);
156  auto spec_shift_down = pred->PredictSyst(calc, shifts);
157  Ratio rat_shift_down(spec_shift_down, spec_nom);
158 
159  double pot = spec_nom.POT();
160 
161  // prep the drawing area
162  auto c = new TCanvas; // just leak this in the demo.
163  TPad * top;
164  TPad * bottom;
165  SplitCanvas(0.3, top, bottom);
166  auto leg = new TLegend(0.6, 0.6, 0.9, 0.9); // also leaked
167  leg->SetBorderSize(0);
168  leg->SetFillStyle(0);
169 
170  // compare the event rate spectra in the top panel
171  top->cd();
172  auto hup = spec_shift_up.ToTH1(pot, kRed, kDashed);
173  hup->GetXaxis()->SetTitleSize(0); // since it will also be drawn on the bottom panel
174  hup->GetXaxis()->SetLabelSize(0);
175  hup->Draw("hist");
176  auto hnom = spec_nom.ToTH1(pot);
177  hnom->Draw("hist same");
178  auto hdown = spec_shift_down.ToTH1(pot, kBlue, kDashed);
179  hdown->Draw("hist same");
180  leg->AddEntry(hnom, "Nominal", "l");
181  leg->AddEntry(hup, "+1#sigma shift");
182  leg->AddEntry(hdown, "-1#sigma shift");
183  leg->Draw();
184 
185  // compare the ratios to nominal in the bottom
186  bottom->cd();
187  auto hbottom = rat_shift_up.ToTH1(kRed);
188  hbottom->GetYaxis()->SetRangeUser(0.45, 1.55);
189  hbottom->Draw("hist");
190  rat_shift_down.ToTH1(kBlue)->Draw("hist same");
191 
192  return c;
193  }
194 
195 }
196 
197 // ---------------------------------------
198 
199 void demoSysts()
200 {
201  // let's use the systs to create some shifted distributions to prove they work.
202 
205 
206  // we'll use ND predictions for simplicity. (no oscillations to confuse us)
207  osc::NoOscillations noOscCalc;
208  auto & NDLoader = loaders.GetLoader(caf::kNEARDET, Loaders::kMC);
209 
210  // ------
211  // predictions for syst #1
212 
213  // we just want to see the effect on the true spectrum, so don't bother with the reco cuts here
214  std::vector<const ISyst*> systs1 { &demo::kDemoSyst1 };
215  SystShifts shifts1(&demo::kDemoSyst1, 1.0);
216  HistAxis trueEaxis("True neutrino energy (GeV)", Binning::Simple(60, 0, 15), kTrueE);
217  NoOscPredictionGenerator predGen_shift1(NDLoader, trueEaxis, kNoCut, kUnweighted);
218  PredictionInterp pred_trueE_syst1(systs1, &noOscCalc, predGen_shift1, loaders, shifts1);
219 
220  // ------
221  // predictions for syst #2
222 
223  // for this one we want to see the effect on the events with neutrons separately
224  const Cut kHasNeutron([](const caf::SRProxy * sr){
225  if (sr->mc.nnu != 1)
226  return false;
227 
228  for (const auto & particle : sr->mc.nu[0].prim)
229  {
230  if (particle.pdg == 2112)
231  return true;
232  }
233 
234  return false;
235  });
236 
237  // here we need the reco cuts and all.
238  std::vector<const ISyst*> systs2 { &demo::kDemoSyst2 };
239  SystShifts shifts2(&demo::kDemoSyst2, 1.0);
240  NoOscPredictionGenerator predGen_numuCC_ND_withNeutron(NDLoader, kNumuCCAxis, kNumuND && kHasNeutron, kUnweighted);
241  PredictionInterp pred_recoE_withNeutron_syst2(systs2, &noOscCalc, predGen_numuCC_ND_withNeutron, loaders, shifts2);
242  NoOscPredictionGenerator predGen_numuCC_ND_noNeutron(NDLoader, kNumuCCAxis, kNumuND && !kHasNeutron, kUnweighted);
243  PredictionInterp pred_recoE_noNeutron_syst2(systs2, &noOscCalc, predGen_numuCC_ND_noNeutron, loaders, shifts2);
244 
245  // ------
246  // make predictions now
247  loaders.Go();
248 
249 
250  // ------
251  // draw comparison plots
252  demo::DrawUpDownRatioCanvas(&pred_trueE_syst1, &noOscCalc, &demo::kDemoSyst1)->Print("test_syst1.png");
253  demo::DrawUpDownRatioCanvas(&pred_recoE_withNeutron_syst2, &noOscCalc, &demo::kDemoSyst2)->Print("test_syst2_neutron.png");
254  demo::DrawUpDownRatioCanvas(&pred_recoE_noNeutron_syst2, &noOscCalc, &demo::kDemoSyst2)->Print("test_syst2_noneutron.png");
255 };
Near Detector underground.
Definition: SREnums.h:10
Pass neutrinos through unchanged.
Definition: IOscCalc.h:45
Implements systematic errors by interpolation between shifted templates.
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TH1D * ToTH1(Color_t col=kBlack, Style_t style=kSolid) const
Definition: Ratio.cxx:67
void TruthShift(double sigma, caf::SRNeutrinoProxy *sr, double &weight) const override
Definition: demoSysts.C:51
const Var weight
Proxy for caf::SRNeutrino.
Definition: SRProxy.h:510
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
caf::Proxy< caf::SRNumuEnergy > numu
Definition: SRProxy.h:214
Proxy for caf::StandardRecord.
Definition: SRProxy.h:2126
General interface to oscillation calculators.
Definition: StanTypedefs.h:23
caf::Proxy< std::vector< caf::SRNeutrino > > nu
Definition: SRProxy.h:618
Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &shift) const override
void Shift(double sigma, caf::SRProxy *sr, double &weight) const override
Perform the systematic shift.
Definition: demoSysts.C:108
virtual Spectrum PredictUnoscillated() const
Definition: IPrediction.cxx:33
caf::Proxy< short int > nnu
Definition: SRProxy.h:617
caf::Proxy< float > E
Definition: SRProxy.h:520
void SplitCanvas(double ysplit, TPad *&p1, TPad *&p2)
Split the current pad into two vertically stacked pieces, suitable for ratio plots and so on...
Definition: Plots.cxx:1419
caf::Proxy< caf::SREnergyBranch > energy
Definition: SRProxy.h:2136
osc::OscCalcDumb calc
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
The simplest job module possible to illustrate the essential features of a job module.
Definition: demoSysts.C:24
const DemoSyst1 kDemoSyst1
Definition: demoSysts.C:89
const DemoSyst2 kDemoSyst2
Definition: demoSysts.C:140
void Go()
Call Go() on all the loaders.
Definition: Loaders.cxx:162
void SetSpillCut(const SpillCut &cut)
Definition: Loaders.cxx:121
const Var kTrueE([](const caf::SRProxy *sr){assert(sr->mc.nnu==1);return sr->mc.nu[0].E;})
Definition: Vars.cxx:85
const HistAxis kNumuCCAxis("Reconstructed Neutrino Energy (GeV)", kNumuEnergyBinning, kCCE)
Definition: HistAxes.h:8
SpectrumLoaderBase & GetLoader(caf::Det_t det, DataMC datamc, DataSource src=kBeam, SwappingConfig swap=kNonSwap)
Retrieve a specific loader.
Definition: Loaders.cxx:129
const Cut kNumuND
Definition: NumuCuts.h:55
#define pot
caf::StandardRecord * sr
caf::Proxy< float > hadcalE
Definition: SRProxy.h:169
double sigma(TH1F *hist, double percentile)
Represent the ratio between two spectra.
Definition: Ratio.h:8
caf::Proxy< caf::SRTruthBranch > mc
Definition: SRProxy.h:2138
std::vector< Loaders * > loaders
Definition: syst_header.h:386
const SpillCut kStandardSpillCuts
Apply this unless you&#39;re doing something special.
Definition: SpillCuts.h:49
For nominal spectra and reweighting systs (xsec/flux)
Definition: Prod4Loaders.h:96
void demoSysts()
Definition: demoSysts.C:199
#define for
Definition: msvc_pragmas.h:3
enum BeamMode kBlue
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:107
const Var kUnweighted
The simplest possible Var, always 1. Used as a default weight.
Definition: Var.h:96
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:81
const Cut kNoCut
The simplest possible cut: pass everything, used as a default.
Definition: Cut.h:109
TGeoVolume * top
Definition: make_fe_box.C:9
TCanvas * DrawUpDownRatioCanvas(const PredictionInterp *pred, osc::IOscCalc *calc, const ISyst *syst)
Definition: demoSysts.C:144