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