demoPlotSystBands.C
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Taken largely from macros by Luke, and Diana
4 // I've trimmed it down a bit for simplicity
5 // if you want further examples like how to include cosmics or plotting the components (CC, NC, WS)
6 // they can be found here:
7 // CAFAna/numu/Analysis2018/BoxOpening/plot_predictions.C
8 // CAFAna/numu/Analysis2018/BoxOpening/plot_datapredictions.C
9 //
10 ////////////////////////////////////////////////////////////////////////////////
11 
12 #include "CAFAna/numu/Analysis2018/includes.h"
13 #include "CAFAna/numu/Analysis2018/settings.h"
14 #include "CAFAna/numu/Analysis2018/numu_tools.h"
16 
17 using namespace ana;
18 
20  bool maxmix = false)
21 {
22  /////////////////////////////////////////////////////////////
23  //
24  // all the setup stuff
25  //
26  /////////////////////////////////////////////////////////////
27 
28  double pot = kAna2018FHCPOT;
30 
31  // some fit points
32  const double dCP = 1.5;
33  const double ssth23 = 0.585;
34  const double dmsq32 = 2.501e-3;
35 
36 
37  // make the calculator
39  calc->SetTh23(asin(sqrt(ssth23)));
40  calc->SetDmsq32(dmsq32);
41 
42 
43  // which systs (choose wisely)
44  std::vector< const ISyst * > systs = {};
45  systs.push_back(&kMuEScaleSyst2017);
46  systs.push_back(&kRelMuEScaleSyst2017);
47  systs.push_back(&kNeutronVisEScalePrimariesSyst2018);
48 
49 
50 
51  /////////////////////////////////////////////////////////////
52  //
53  // Loading in predictions and making experiments
54  // you could also load your real data here along with cosmics
55  //
56  /////////////////////////////////////////////////////////////
57  std::vector< PredictionSystJoint2018* > predictions;
58  std::vector< Spectrum > s_fakedata;
59  predictions.push_back(new PredictionSystJoint2018(kNuMu, calc, "fhc", -1, systs, "/nova/ana/nu_mu_ana/Ana2018/Predictions/pred_numuconcat_fhc__numu2018.root"));
60  s_fakedata.push_back(predictions.back()->Predict(calc).FakeData(pot));
61  for(int quant = 1; quant < 5; ++quant){
62  predictions.push_back(new PredictionSystJoint2018(kNuMu, calc, "fhc", quant, systs, "/nova/ana/nu_mu_ana/Ana2018/Predictions/pred_numuconcat_fhc__numu2018.root"));
63  s_fakedata.push_back(predictions.back()->Predict(calc).FakeData(pot));
64  }
65 
66 
67  /////////////////////////////////////////////////////////////
68  //
69  // doing some real things
70  //
71  /////////////////////////////////////////////////////////////
72 
73  // Since here we want to scale the spectra to Events / 0.1 GeV
74  // we have to make histograms from the predictions and data
75  // otherwise you could use the PlotWithSystErrorBand_Quant constructor that
76  // takes the whole ISyst vector and all the work would be done right now
77  std::vector<TH1*> hDataUnscaled;
78  std::vector<TH1*> hData;
79  std::vector<TH1*> hPred;
80  for(unsigned int quant = 0; quant < 5; quant++){
81  hPred.push_back(predictions[quant]->Predict(calc).ToTH1(pot));
82  hData.push_back(s_fakedata[quant].ToTH1(pot));
83  hDataUnscaled.push_back(s_fakedata[quant].ToTH1(pot));
84  }
85 
86  // scale that shiz
87  for(unsigned int quant = 0; quant < 5; quant++){
88  hData[quant]->Scale(0.1, "width");
89  hPred[quant]->Scale(0.1, "width");
90  }
91 
92  // data needs asymm graph and scaled predictions need to be put in a spectrum
93  std::vector<TGraph*> grData;
94  std::vector<Spectrum> s_predictions_scaled;
95  for(unsigned int quant = 0; quant < 5; quant++){
96  grData.push_back(graphAsymmErrorScaled(hData[quant], hDataUnscaled[quant]));
97  s_predictions_scaled.push_back(Spectrum(hPred[quant], pot, livetime));
98  }
99 
100  // Keep track of up and dn shifts like this because of the scaling
101  std::vector<Spectrum> ups[5], dns[5];
102  for(unsigned int quant = 0; quant < 5; quant++){
103  for(const ISyst* syst: systs){
104  SystShifts shifts;
105 
106  shifts.SetShift(syst, +1);
107  TH1* htempup = predictions[quant]->PredictSyst(calc, shifts).ToTH1(pot);
108  htempup->Scale(0.1, "width");
109  ups[quant].push_back(Spectrum(htempup, pot, livetime));
110 
111  shifts.SetShift(syst, -1);
112  TH1* htempdn = predictions[quant]->PredictSyst(calc, shifts).ToTH1(pot);
113  htempdn->Scale(0.1, "width");
114  dns[quant].push_back(Spectrum(htempdn, pot, livetime));
115  }// systs loop
116  }// quantiles loop
117 
118 
119  /////////////////////////////////////////////////////////////
120  //
121  // plot systs as a group
122  //
123  /////////////////////////////////////////////////////////////
124  // axis labels for later
125  TLatex* xTitle = new TLatex(0.5, 0.03, "Reconstructed Energy (GeV)");
126  xTitle->SetNDC();
127  xTitle->SetTextAlign(22);
128  TLatex* yTitle = new TLatex(0.02, 0.5, "Events / 0.1 GeV");
129  yTitle->SetNDC();
130  yTitle->SetTextAlign(22);
131  yTitle->SetTextAngle(90.);
132 
133  std::string n_allquants = "all_quants";
134  std::string n_spltquant = "splt_quant";
135 
136  // first do all quants
137  TCanvas *c_allquants = new TCanvas("c_allquants", n_allquants.c_str(), 600, 500);
138  PlotWithSystErrorBand_Quant(0, s_predictions_scaled[0], ups[0], dns[0], pot, kViolet+1, kViolet-9, 11, true);
139  grData[0]->Draw("p0 same");
140  drawQuantLabel(0);
141  xTitle->Draw();
142  yTitle->Draw();
143  TLegend *leg_allquants = new TLegend(0.52,0.70,0.98,0.85,"","NDC");
144  leg_allquants->AddEntry(hPred[0], "Prediction","l");
145  TLegendEntry *entry_allquants=leg_allquants->AddEntry("error","1-#sigma syst.","bf");
146  entry_allquants->SetFillStyle(1001);
147  entry_allquants->SetFillColor(color_band);
148  entry_allquants->SetLineColor(color_band);
149  leg_allquants->SetTextSize(0.05); //no border for legend
150  leg_allquants->SetBorderSize(0); //no border for legend
151  leg_allquants->SetFillStyle(0); //fill colour is white
152  leg_allquants->Draw();
153 
154  std::string label_name = "All Systs";
155  CornerLabel(label_name);
156  c_allquants->SaveAs(("preds_syst_bands_" + n_allquants + ".png").c_str());
157 
158  // now for the split quantiles
159  TCanvas *c_spltquant;
160  TPad *pad1, *pad2, *pad3, *pad4;
161  SplitCanvasQuant(c_spltquant, pad1, pad2, pad3, pad4);
162  TLegend *leg_spltquant = new TLegend(0.09,0.78,0.27,0.88,"","NDC");
163  leg_spltquant->AddEntry(hPred[1], "Prediction","l");
164  TLegendEntry *entry_spltquant=leg_spltquant->AddEntry("error","1-#sigma syst.","bf");
165  entry_spltquant->SetFillStyle(1001);
166  entry_spltquant->SetFillColor(color_band);
167  entry_spltquant->SetLineColor(color_band);
168  leg_spltquant->SetTextSize(0.025); //no border for legend
169  leg_spltquant->SetBorderSize(0); //no border for legend
170  leg_spltquant->SetFillStyle(0); //fill colour is white
171 
172  for(unsigned int quant = 1; quant < 5; quant++){
173  if(quant == 1) pad1->cd();
174  if(quant == 2) pad2->cd();
175  if(quant == 3) pad3->cd();
176  if(quant == 4) pad4->cd();
177  PlotWithSystErrorBand_Quant(quant, s_predictions_scaled[quant], ups[quant], dns[quant], pot, kViolet+1, kViolet-9, 3.5, true);
178  grData[quant]->Draw("p0 same");
179  drawQuantLabel(quant);
180  }
181 
182  for(auto & p:{pad1,pad2,pad3,pad4}){
183  p->RedrawAxis();
184  }
185  pad1->cd();
186  leg_spltquant->Draw();
187  CornerLabel(label_name);
188  c_spltquant->SetTitle(n_spltquant.c_str());
189  c_spltquant->Update();
190  c_spltquant->cd();
191  xTitle->Draw();
192  yTitle->Draw();
193  c_spltquant->SaveAs(("preds_syst_bands_" + n_spltquant + ".png").c_str());
194 
195 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
double ssth23
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
Loads shifted spectra from files.
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
void CornerLabel(std::string &s)
Definition: numu_tools.h:145
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
virtual void SetDmsq32(const T &dmsq32)=0
const MuEScaleSyst2017 kMuEScaleSyst2017(0.0074, 0.0012)
const RelMuEScaleSyst2017 kRelMuEScaleSyst2017(0.0045, 10.5)
TGraph * graphAsymmErrorScaled(TH1 *histScaled, TH1 *hist, double overallScale)
Definition: Plots.cxx:980
double dCP
#define pot
std::vector< float > Spectrum
Definition: Constants.h:728
double dmsq32
double livetime
Definition: saveFDMCHists.C:21
TGraphAsymmErrors * PlotWithSystErrorBand_Quant(const int quant, IPrediction *pred, const std::vector< const ISyst * > &systs, osc::IOscCalc *calc, double pot, int col, int errCol, float maxy, bool newaxis)
Definition: Plots.cxx:543
enum BeamMode kViolet
TPaveText * drawQuantLabel(int quantId=0)
Definition: numu_tools.h:4
virtual void SetTh23(const T &th23)=0
TPad * pad3
Definition: analysis.C:13
void SplitCanvasQuant(TCanvas *&canvas, TPad *&pad1, TPad *&pad2, TPad *&pad3, TPad *&pad4)
Definition: Plots.cxx:1439
const double kAna2018FHCPOT
Definition: Exposures.h:207
TPad * pad2
Definition: analysis.C:13
void demoPlotSystBands(std::string beam="fhc", bool maxmix=false)
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:81
Color_t color_band
Definition: saveFDMCHists.C:46
const double kAna2018FHCLivetime
Definition: Exposures.h:213
TPad * pad1
Definition: analysis.C:13
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string