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