drawSystsShiftingNDdata.C
Go to the documentation of this file.
1 #include "CAFAna/Cuts/Cuts.h"
6 #include "CAFAna/Core/Loaders.h"
8 #include "CAFAna/Core/Spectrum.h"
10 #include "CAFAna/Vars/Vars.h"
13 #include "OscLib/OscCalcPMNSOpt.h"
14 #include "CAFAna/Analysis/Calcs.h"
16 #include "CAFAna/Analysis/Plots.h"
17 
27 #include "CAFAna/Vars/FitVars.h"
28 
29 #include "TStyle.h"
30 #include "TCanvas.h"
31 #include "TGraph.h"
32 #include "TGraphAsymmErrors.h"
33 #include "TH1.h"
34 #include "TH2.h"
35 #include "TLegend.h"
36 #include "TLatex.h"
37 #include "TArrow.h"
38 
39 using namespace ana;
40 
41 // Put a "NOvA Simulation" tag in the corner
42 void Simulation()
43 {
44  TLatex* prelim = new TLatex(.9, .95, "NOvA Simulation");
45  prelim->SetTextColor(kGray+1);
46  prelim->SetNDC();
47  prelim->SetTextSize(2/30.);
48  prelim->SetTextAlign(32);
49  prelim->Draw();
50 }
51 
52 TH1* getResidualDiff(TH1* h1, TH1* h2)
53 {
54  TH1* h = (TH1*) h1->Clone("h");
55  h->Reset();
56 
57  for(int k = 1; k <= h1->GetNbinsX(); k++)
58  {
59  double val1 = (double) h1->GetBinContent(k);
60  double val2 = (double) h2->GetBinContent(k);
61 
62  if(val1 * val2 > 0)
63  {
64  h->SetBinContent(k, (val1-val2) * 200./(val1 + val2));
65  }
66  }
67  return h;
68 }
69 
71 {
72 
73  std::cout << "Drawing plots..." << std::endl;
74 
75  gStyle->SetStatFormat("6.3g"); // Significant figures. Fix if this is not satisfactory
76  gStyle->SetPaintTextFormat("6.3g");
77  gStyle->SetFitFormat("6.3g");
78 
81 
82  calc.SetTh23(TMath::ASin(sqrt(0.404))); // Published second ana
83  calc.SetDmsq32(2.67e-3); // Same
84 
85  static const MECq0ShapeSyst2017 kMECq0ShapeSyst;
86  static const MECInitStateNPFracSyst2017 kMEInitStateNPFracSyst;
87  static const MECEnuShapeSyst2017 kMEEnuShapeSyst;
88 
89  std::vector<const ISyst*> systs = getAllNumuSystsSA();
90 
91  systs.push_back(&kMECq0ShapeSyst);
92  systs.push_back(&kMEInitStateNPFracSyst);
93  systs.push_back(&kMEEnuShapeSyst);
94 
95  const unsigned int kNumSysts = systs.size();
96 
97  PredictionCombinePeriods* predictionND[2*kNumSysts];
98 
99  PredictionInterp* predictionND_period1[2*kNumSysts];
100  PredictionInterp* predictionND_period2[2*kNumSysts];
101  PredictionInterp* predictionND_epoch3b[2*kNumSysts];
102  PredictionInterp* predictionND_epoch3c[2*kNumSysts];
103 
104  PredictionCombinePeriods* predictionFD[2*kNumSysts];
105 
106  PredictionNoExtrap* predictionFD_period1[2*kNumSysts];
107  PredictionNoExtrap* predictionFD_period2[2*kNumSysts];
108  PredictionNoExtrap* predictionFD_epoch3b[2*kNumSysts];
109  PredictionNoExtrap* predictionFD_epoch3c[2*kNumSysts];
110 
111  PredictionCombinePeriods* nominalPred;
112 
113  PredictionNoExtrap* nominalPred_period1;
114  PredictionNoExtrap* nominalPred_period2;
115  PredictionNoExtrap* nominalPred_epoch3b;
116  PredictionNoExtrap* nominalPred_epoch3c;
117 
118  TFile* fPred_period1 = new TFile( "Predictions/Prediction_period1_ccE.root" );
119  TFile* fPred_period2 = new TFile( "Predictions/Prediction_period2_ccE.root" );
120  TFile* fPred_epoch3b = new TFile( "Predictions/Prediction_epoch3b_ccE.root" );
121  TFile* fPred_epoch3c = new TFile( "Predictions/Prediction_epoch3c_ccE.root" );
122 
123  std::string predNameND = "ccE_shiftNDdata";
124  std::string predNameFD = "ccE_noExtrapFD";
125 
126  nominalPred_period1 = LoadFrom<PredictionNoExtrap>(fPred_period1, "nominalPred").release();
127  nominalPred_period2 = LoadFrom<PredictionNoExtrap>(fPred_period2, "nominalPred").release();
128  nominalPred_epoch3b = LoadFrom<PredictionNoExtrap>(fPred_epoch3b, "nominalPred").release();
129  nominalPred_epoch3c = LoadFrom<PredictionNoExtrap>(fPred_epoch3c, "nominalPred").release();
130 
131  nominalPred = new PredictionCombinePeriods({{nominalPred_period1, kSecondAnaPeriod1POT},
132  {nominalPred_period2, kSecondAnaPeriod2POT},
133  {nominalPred_epoch3b, kSecondAnaEpoch3bPOT},
134  {nominalPred_epoch3c, kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT}});
135 
136  for(unsigned int k = 0; k < kNumSysts; k++)
137  {
138  predictionND_period1[2*k] = LoadFrom<PredictionInterp>(fPred_period1, (predNameND+"_Pos_"+systs[k]->ShortName()).c_str() ).release();
139 
140  predictionND_period2[2*k] = LoadFrom<PredictionInterp>(fPred_period2, (predNameND+"_Pos_"+systs[k]->ShortName()).c_str() ).release();
141 
142  predictionND_epoch3b[2*k] = LoadFrom<PredictionInterp>(fPred_epoch3b, (predNameND+"_Pos_"+systs[k]->ShortName()).c_str() ).release();
143 
144  predictionND_epoch3c[2*k] = LoadFrom<PredictionInterp>(fPred_epoch3c, (predNameND+"_Pos_"+systs[k]->ShortName()).c_str() ).release();
145 
146  predictionND_period1[2*k+1] = LoadFrom<PredictionInterp>(fPred_period1, (predNameND+"_Neg_"+systs[k]->ShortName()).c_str() ).release();
147 
148  predictionND_period2[2*k+1] = LoadFrom<PredictionInterp>(fPred_period2, (predNameND+"_Neg_"+systs[k]->ShortName()).c_str() ).release();
149 
150  predictionND_epoch3b[2*k+1] = LoadFrom<PredictionInterp>(fPred_epoch3b, (predNameND+"_Neg_"+systs[k]->ShortName()).c_str() ).release();
151 
152  predictionND_epoch3c[2*k+1] = LoadFrom<PredictionInterp>(fPred_epoch3c, (predNameND+"_Neg_"+systs[k]->ShortName()).c_str() ).release();
153 
154  predictionND[2*k] = new PredictionCombinePeriods({ {predictionND_period1[2*k], kSecondAnaPeriod1POT},
155  {predictionND_period2[2*k], kSecondAnaPeriod2POT},
156  {predictionND_epoch3b[2*k], kSecondAnaEpoch3bPOT},
157  {predictionND_epoch3c[2*k], kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT} });
158 
159  predictionND[2*k+1] = new PredictionCombinePeriods({ {predictionND_period1[2*k+1], kSecondAnaPeriod1POT},
160  {predictionND_period2[2*k+1], kSecondAnaPeriod2POT},
161  {predictionND_epoch3b[2*k+1], kSecondAnaEpoch3bPOT},
162  {predictionND_epoch3c[2*k+1], kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT} });
163 
164 
165  // Now FD (no extrap)
166 
167  predictionFD_period1[2*k] = LoadFrom<PredictionNoExtrap>(fPred_period1, (predNameFD+"_Pos_"+systs[k]->ShortName()).c_str() ).release();
168 
169  predictionFD_period2[2*k] = LoadFrom<PredictionNoExtrap>(fPred_period2, (predNameFD+"_Pos_"+systs[k]->ShortName()).c_str() ).release();
170 
171  predictionFD_epoch3b[2*k] = LoadFrom<PredictionNoExtrap>(fPred_epoch3b, (predNameFD+"_Pos_"+systs[k]->ShortName()).c_str() ).release();
172 
173  predictionFD_epoch3c[2*k] = LoadFrom<PredictionNoExtrap>(fPred_epoch3c, (predNameFD+"_Pos_"+systs[k]->ShortName()).c_str() ).release();
174 
175  predictionFD_period1[2*k+1] = LoadFrom<PredictionNoExtrap>(fPred_period1, (predNameFD+"_Neg_"+systs[k]->ShortName()).c_str() ).release();
176 
177  predictionFD_period2[2*k+1] = LoadFrom<PredictionNoExtrap>(fPred_period2, (predNameFD+"_Neg_"+systs[k]->ShortName()).c_str() ).release();
178 
179  predictionFD_epoch3b[2*k+1] = LoadFrom<PredictionNoExtrap>(fPred_epoch3b, (predNameFD+"_Neg_"+systs[k]->ShortName()).c_str() ).release();
180 
181  predictionFD_epoch3c[2*k+1] = LoadFrom<PredictionNoExtrap>(fPred_epoch3c, (predNameFD+"_Neg_"+systs[k]->ShortName()).c_str() ).release();
182 
183  predictionFD[2*k] = new PredictionCombinePeriods({ {predictionFD_period1[2*k], kSecondAnaPeriod1POT},
184  {predictionFD_period2[2*k], kSecondAnaPeriod2POT},
185  {predictionFD_epoch3b[2*k], kSecondAnaEpoch3bPOT},
186  {predictionFD_epoch3c[2*k], kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT} });
187 
188  predictionFD[2*k+1] = new PredictionCombinePeriods({ {predictionFD_period1[2*k+1], kSecondAnaPeriod1POT},
189  {predictionFD_period2[2*k+1], kSecondAnaPeriod2POT},
190  {predictionFD_epoch3b[2*k+1], kSecondAnaEpoch3bPOT},
191  {predictionFD_epoch3c[2*k+1], kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT} });
192  }
193 
194  std::cout << "Loading went OK. Plotting!" << std::endl;
195 
196  TH1* hNominal = nominalPred->Predict(&calc).ToTH1(kSecondAnaPOT);
197 
198  hNominal->SetLineColor(kBlack);
199  hNominal->GetYaxis()->SetRangeUser(0, 30);
200  hNominal->GetYaxis()->SetTitle("Events / (0.25 GeV)");
201  hNominal->GetYaxis()->CenterTitle();
202  hNominal->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
203  hNominal->GetXaxis()->CenterTitle();
204 
205  for (unsigned int k = 0; k < kNumSysts; k++)
206  {
207  std::string fullName = "ccE_"+systs[k]->ShortName();
208 
209  TCanvas* c = new TCanvas("c",(fullName).c_str());
210 
211  c->SetBottomMargin(0.15);
212  c->SetLeftMargin(0.15);
213 
214  TH1* hShiftND_up = predictionND[2*k]->Predict(&calc).ToTH1(kSecondAnaPOT);
215  TH1* hShiftFD_up = predictionFD[2*k]->Predict(&calc).ToTH1(kSecondAnaPOT);
216 
217  TH1* hShiftND_dn = predictionND[2*k+1]->Predict(&calc).ToTH1(kSecondAnaPOT);
218  TH1* hShiftFD_dn = predictionFD[2*k+1]->Predict(&calc).ToTH1(kSecondAnaPOT);
219 
220  hShiftND_up->SetLineColor(kOrange+1);
221  hShiftND_up->SetLineStyle(2);
222 
223  hShiftND_dn->SetLineColor(kViolet+1);
224  hShiftND_dn->SetLineStyle(2);
225 
226  hShiftFD_up->SetMarkerStyle(20);
227  hShiftFD_up->SetMarkerColor(kOrange+1);
228  hShiftFD_up->SetLineWidth(0);
229  hShiftFD_up->SetLineColor(0);
230  hShiftFD_dn->SetMarkerStyle(21);
231  hShiftFD_dn->SetMarkerColor(kViolet+1);
232  hShiftFD_dn->SetLineWidth(0);
233  hShiftFD_dn->SetLineColor(0);
234 
235  hNominal->Draw("hist");
236  hShiftND_up->Draw("hist same");
237  hShiftND_dn->Draw("hist same");
238  hShiftFD_up->Draw("EX0 same");
239  hShiftFD_dn->Draw("EX0 same");
240 
241  TLegend* leg = new TLegend(0.4,0.6,0.88,0.88);
242  leg->AddEntry(hNominal, "Nominal with FD/ND","l");
243  leg->AddEntry(hShiftND_up, TString("+#sigma in " + systs[k]->LatexName()) + " in ND","l");
244  leg->AddEntry(hShiftFD_up, TString("+#sigma in " + systs[k]->LatexName()) + " in FD","p");
245  leg->AddEntry(hShiftND_dn, TString("-#sigma in " + systs[k]->LatexName()) + " in ND","l");
246  leg->AddEntry(hShiftFD_dn, TString("-#sigma in " + systs[k]->LatexName()) + " in FD","p");
247  leg->SetFillColorAlpha(10,0.);
248  leg->Draw();
249 
250  Simulation();
251 
252  c->SaveAs(("plots/"+fullName+".root").c_str());
253  c->SaveAs(("plots/"+fullName+".pdf").c_str());
254 
255  // Now the ratios
256 
257  c->Clear();
258 
259  hShiftND_up->Divide(hNominal);
260  hShiftND_dn->Divide(hNominal);
261 
262  hShiftFD_up->Divide(hNominal);
263  hShiftFD_dn->Divide(hNominal);
264 
265  hShiftND_up->GetYaxis()->SetRangeUser(0.8,1.2);
266  //if(k == 0 && (kNumSysts <= 4)) hShiftND_up->GetYaxis()->SetRangeUser(0.5, 2.0);
267  hShiftND_up->GetYaxis()->SetTitle("Ratio to nominal MC");
268  hShiftND_up->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
269  hShiftND_up->GetYaxis()->CenterTitle();
270  hShiftND_up->GetXaxis()->CenterTitle();
271 
272  hShiftND_up->SetMarkerColor(kOrange+1);
273  hShiftND_up->SetLineWidth(2);
274  hShiftND_up->SetLineStyle(1);
275  hShiftND_up->SetLineColor(kOrange+1);
276  hShiftND_dn->SetMarkerColor(kOrange+1);
277  hShiftND_dn->SetLineWidth(2);
278  hShiftND_dn->SetLineStyle(1);
279  hShiftND_dn->SetLineColor(kOrange+1);
280 
281  hShiftFD_up->SetMarkerColor(kViolet+1);
282  hShiftFD_up->SetLineWidth(2);
283  hShiftFD_up->SetLineStyle(2);
284  hShiftFD_up->SetLineColor(kViolet+1);
285  hShiftFD_dn->SetMarkerColor(kViolet+1);
286  hShiftFD_dn->SetLineWidth(2);
287  hShiftFD_dn->SetLineStyle(2);
288  hShiftFD_dn->SetLineColor(kViolet+1);
289 
290  hShiftND_up->Draw("hist");
291  hShiftND_dn->Draw("hist same");
292  hShiftFD_up->Draw("hist same");
293  hShiftFD_dn->Draw("hist same");
294 
295  TGraph* grShade = ShadeBetweenHistograms(hShiftND_up, hShiftFD_up);
296  grShade->SetFillColorAlpha(kRed-9, 0.35);
297  grShade->SetLineColor(kRed-9);
298  grShade->Draw("F");
299 
300  TGraph* grShade2 = ShadeBetweenHistograms(hShiftND_dn, hShiftFD_dn);
301  grShade2->SetFillColorAlpha(kRed-9, 0.35);
302  grShade2->SetLineColor(kRed-9);
303  grShade2->Draw("F");
304 
305  TGraph* g1 = new TGraph;
306  g1->SetPoint(0, -1, 1);
307  g1->SetPoint(1, 10, 1);
308 
309  g1->SetLineWidth(1);
310  g1->SetLineColor(kBlack);
311  g1->SetLineStyle(2);
312 
313  g1->Draw("L");
314 
315  leg->Clear();
316  leg->AddEntry(hShiftND_up, TString("#pm#sigma extrap. ND shift in " + systs[k]->LatexName()),"l");
317  leg->AddEntry(hShiftFD_up, TString("#pm#sigma FD shift " + systs[k]->LatexName()),"l");
318  leg->SetFillColorAlpha(10,0.);
319  leg->Draw();
320 
321  Simulation();
322 
323  c->SaveAs(("plots/"+fullName+"_ratio.root").c_str());
324  c->SaveAs(("plots/"+fullName+"_ratio.pdf").c_str());
325 
326  // And finally the residual difference
327 
328  c->Clear();
329 
330  TH1* hResup = getResidualDiff(hShiftND_up, hShiftFD_up);
331 
332  TH1* hResdn = getResidualDiff(hShiftND_dn, hShiftFD_dn);
333 
334  hResup->GetYaxis()->SetTitle("Residual difference (%)");
335  hResup->GetYaxis()->SetRangeUser( -10, 10 );
336  //if(k == 0 && (kNumSysts <= 4)) hResup->GetYaxis()->SetRangeUser(-50, 50);
337 
338  hResup->SetLineColor(kOrange+1);
339  hResdn->SetLineColor(kViolet+1);
340 
341  hResup->Draw("hist");
342  hResdn->Draw("hist same");
343 
344  leg->Clear();
345  leg->AddEntry(hResup, TString("+ #sigma shift in " + systs[k]->LatexName() + " FD minus ND"),"l");
346  leg->AddEntry(hResdn, TString("- #sigma shift in " + systs[k]->LatexName() + " FD minus ND"),"l");
347  leg->SetFillColorAlpha(10,0.);
348  leg->Draw();
349 
350  TGraph* g0 = new TGraph;
351  g0->SetPoint(0, -1, 0);
352  g0->SetPoint(1, 10, 0);
353 
354  g0->SetLineWidth(1);
355  g0->SetLineColor(kBlack);
356  g0->SetLineStyle(2);
357 
358  g0->Draw("L");
359 
360  Simulation();
361 
362  c->SaveAs(("plots/"+fullName+"_residual.root").c_str());
363  c->SaveAs(("plots/"+fullName+"_residual.pdf").c_str());
364 
365  c->Close();
366  }
367 
368 
369  fPred_period1->Close();
370  fPred_period2->Close();
371  fPred_epoch3b->Close();
372  fPred_epoch3c->Close();
373 
374  std::cout << "All plots finished!" << std::endl;
375 
376 
377 } // End of function
enum BeamMode kOrange
Implements systematic errors by interpolation between shifted templates.
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const double kSecondAnaPeriod2POT
Definition: Exposures.h:74
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:148
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double kSecondAnaEpoch3dPOT
Definition: Exposures.h:77
void Simulation()
const double kSecondAnaEpoch3bPOT
Definition: Exposures.h:75
const double kSecondAnaPeriod1POT
Definition: Exposures.h:73
osc::OscCalcDumb calc
TH1 * getResidualDiff(TH1 *h1, TH1 *h2)
TGraph * ShadeBetweenHistograms(TH1 *hmin, TH1 *hmax)
Definition: Plots.cxx:1236
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
void SetTh23(const T &th23) override
TLatex * prelim
Definition: Xsec_final.C:133
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
void SetDmsq32(const T &dmsq32) override
TH1F * h1
enum BeamMode kViolet
const double kSecondAnaEpoch3cPOT
Definition: Exposures.h:76
Sum MC predictions from different periods scaled according to data POT targets.
Prediction that just uses FD MC, with no extrapolation.
Float_t e
Definition: plot.C:35
const double kSecondAnaPOT
Definition: Exposures.h:87
void drawSystsShiftingNDdata(std::string hornCurrStr, std::string predFile, std::string outDir)
Spectrum Predict(osc::IOscCalc *calc) const override
enum BeamMode string