drawSystsShiftingNDdata_updatedAna.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 const unsigned int kNumQuantiles = 4;
42 
43 // Put a "NOvA Simulation" tag in the corner
44 void Simulation()
45 {
46  TLatex* prelim = new TLatex(.9, .95, "NOvA Simulation");
47  prelim->SetTextColor(kGray+1);
48  prelim->SetNDC();
49  prelim->SetTextSize(2/30.);
50  prelim->SetTextAlign(32);
51  prelim->Draw();
52 }
53 
54 void Quartile(int q)
55 {
56  TLatex* quart = new TLatex(.2, .8, q < 4? Form("Q%i",q+1) : "All Qs.");
57  quart->SetNDC();
58  quart->SetTextSize(2/30.);
59  quart->Draw();
60 }
61 
62 TH1* getResidualDiff(TH1* h1, TH1* h2)
63 {
64  TH1* h = (TH1*) h1->Clone("h");
65  h->Reset();
66 
67  for(int k = 1; k <= h1->GetNbinsX(); k++)
68  {
69  double val1 = (double) h1->GetBinContent(k);
70  double val2 = (double) h2->GetBinContent(k);
71 
72  if(val1 * val2 > 0)
73  {
74  h->SetBinContent(k, (val1-val2) * 200./(val1 + val2));
75  }
76  }
77  return h;
78 }
79 
81 {
82 
83  std::cout << "Drawing plots..." << std::endl;
84 
85  gStyle->SetStatFormat("6.3g"); // Significant figures. Fix if this is not satisfactory
86  gStyle->SetPaintTextFormat("6.3g");
87  gStyle->SetFitFormat("6.3g");
88 
91 
92  calc.SetTh23(TMath::ASin(sqrt(0.404))); // Published second ana
93  calc.SetDmsq32(2.67e-3); // Same
94 
95  std::vector<const ISyst*> systs = getAllNumuSystsSA();
96 
97  const unsigned int kNumSysts = systs.size();
98 
99  PredictionInterp* predictionND[kNumSysts*2][kNumQuantiles+1];
100  PredictionNoExtrap* predictionFD[kNumSysts*2][kNumQuantiles+1];
101  PredictionNoExtrap* nominalPred[kNumQuantiles+1];
102 
103  TFile* fPred = new TFile( "Predictions/Prediction_period2_ccE_1NDQuantile.root" );
104 
105  std::string predNameND = "ccE_shiftNDdata";
106  std::string predNameFD = "ccE_noExtrapFD";
107 
108  for(unsigned int i = 0; i < kNumQuantiles +1; i++)
109  {
110  nominalPred[i] = LoadFrom<PredictionNoExtrap>(fPred, Form("nominalPred_%i",i)).release();
111 
112  for(unsigned int k = 0; k < kNumSysts; k++)
113  {
114  predictionND[2*k][i] = LoadFrom<PredictionInterp>(fPred, Form( (predNameND+"_Pos_"+systs[k]->ShortName()+"_%i").c_str(), i) ).release();
115 
116  predictionND[2*k+1][i] = LoadFrom<PredictionInterp>(fPred, Form( (predNameND+"_Neg_"+systs[k]->ShortName()+"_%i").c_str(), i) ).release();
117 
118  // Now FD (no extrap)
119 
120  predictionFD[2*k][i] = LoadFrom<PredictionNoExtrap>(fPred, Form( (predNameFD+"_Pos_"+systs[k]->ShortName()+"_%i").c_str(), i) ).release();
121 
122  predictionFD[2*k+1][i] = LoadFrom<PredictionNoExtrap>(fPred, Form( (predNameFD+"_Neg_"+systs[k]->ShortName()+"_%i").c_str(), i) ).release();
123  } // end loop ver systs
124  } // end loop over quantiles
125 
126  std::cout << "Loading went OK. Plotting!" << std::endl;
127 
128  for(unsigned int i = 0; i < kNumQuantiles +1; i++)
129  {
130  TH1* hNominal = nominalPred[i]->Predict(&calc).ToTH1(kAna2017POT);
131  hNominal->Scale(0.1, "width");
132 
133  hNominal->SetLineColor(kBlack);
134  double ymax = hNominal->GetBinContent(hNominal->GetMaximumBin()) * 2.0;
135  hNominal->GetYaxis()->SetRangeUser(0, 15);
136  hNominal->GetYaxis()->SetTitle("Events / (0.1 GeV)");
137  hNominal->GetYaxis()->CenterTitle();
138  hNominal->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
139  hNominal->GetXaxis()->CenterTitle();
140 
141  for (unsigned int k = 0; k < kNumSysts; k++)
142  {
143  std::string fullName = "ccE_"+systs[k]->ShortName();
144 
145  TCanvas* c = new TCanvas("c",(fullName).c_str());
146 
147  c->SetBottomMargin(0.15);
148  c->SetLeftMargin(0.15);
149 
150  TH1* hShiftND_up = predictionND[2*k][i]->Predict(&calc).ToTH1(kAna2017POT);
151  hShiftND_up->Scale(0.1, "width");
152  TH1* hShiftFD_up = predictionFD[2*k][i]->Predict(&calc).ToTH1(kAna2017POT);
153  hShiftFD_up->Scale(0.1, "width");
154 
155  TH1* hShiftND_dn = predictionND[2*k+1][i]->Predict(&calc).ToTH1(kAna2017POT);
156  hShiftND_dn->Scale(0.1, "width");
157  TH1* hShiftFD_dn = predictionFD[2*k+1][i]->Predict(&calc).ToTH1(kAna2017POT);
158  hShiftFD_dn->Scale(0.1, "width");
159 
160  hShiftND_up->SetLineColor(kOrange+1);
161  hShiftND_up->SetLineStyle(2);
162 
163  hShiftND_dn->SetLineColor(kViolet+1);
164  hShiftND_dn->SetLineStyle(2);
165 
166  hShiftFD_up->SetMarkerStyle(20);
167  hShiftFD_up->SetMarkerColor(kOrange+1);
168  hShiftFD_up->SetLineWidth(0);
169  hShiftFD_up->SetLineColor(0);
170  hShiftFD_dn->SetMarkerStyle(21);
171  hShiftFD_dn->SetMarkerColor(kViolet+1);
172  hShiftFD_dn->SetLineWidth(0);
173  hShiftFD_dn->SetLineColor(0);
174 
175  hNominal->Draw("hist");
176  hShiftND_up->Draw("hist same");
177  hShiftND_dn->Draw("hist same");
178  hShiftFD_up->Draw("EX0 same");
179  hShiftFD_dn->Draw("EX0 same");
180 
181  TLegend* leg = new TLegend(0.4,0.6,0.88,0.88);
182  leg->AddEntry(hNominal, "Nominal with FD/ND","l");
183  leg->AddEntry(hShiftND_up, TString("+#sigma in " + systs[k]->LatexName()) + " in ND","l");
184  leg->AddEntry(hShiftFD_up, TString("+#sigma in " + systs[k]->LatexName()) + " in FD","p");
185  leg->AddEntry(hShiftND_dn, TString("-#sigma in " + systs[k]->LatexName()) + " in ND","l");
186  leg->AddEntry(hShiftFD_dn, TString("-#sigma in " + systs[k]->LatexName()) + " in FD","p");
187  leg->SetFillColorAlpha(10,0.);
188  leg->Draw();
189 
190  Simulation();
191  Quartile(i);
192 
193  c->SaveAs(Form( ("plots/"+fullName+"_Q%i.root").c_str(), i));
194  c->SaveAs(Form( ("plots/"+fullName+"_Q%i.pdf").c_str(), i));
195 
196  // Now the ratios
197 
198  c->Clear();
199 
200  hShiftND_up->Divide(hNominal);
201  hShiftND_dn->Divide(hNominal);
202 
203  hShiftFD_up->Divide(hNominal);
204  hShiftFD_dn->Divide(hNominal);
205 
206  hShiftND_up->GetYaxis()->SetRangeUser(0.8,1.2);
207  //if(k == 0 && (kNumSysts <= 4)) hShiftND_up->GetYaxis()->SetRangeUser(0.5, 2.0);
208  hShiftND_up->GetYaxis()->SetTitle("Ratio to nominal MC");
209  hShiftND_up->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
210  hShiftND_up->GetYaxis()->CenterTitle();
211  hShiftND_up->GetXaxis()->CenterTitle();
212 
213  hShiftND_up->SetMarkerColor(kOrange+1);
214  hShiftND_up->SetLineWidth(2);
215  hShiftND_up->SetLineStyle(1);
216  hShiftND_up->SetLineColor(kOrange+1);
217  hShiftND_dn->SetMarkerColor(kOrange+1);
218  hShiftND_dn->SetLineWidth(2);
219  hShiftND_dn->SetLineStyle(1);
220  hShiftND_dn->SetLineColor(kOrange+1);
221 
222  hShiftFD_up->SetMarkerColor(kViolet+1);
223  hShiftFD_up->SetLineWidth(2);
224  hShiftFD_up->SetLineStyle(2);
225  hShiftFD_up->SetLineColor(kViolet+1);
226  hShiftFD_dn->SetMarkerColor(kViolet+1);
227  hShiftFD_dn->SetLineWidth(2);
228  hShiftFD_dn->SetLineStyle(2);
229  hShiftFD_dn->SetLineColor(kViolet+1);
230 
231  hShiftND_up->Draw("hist");
232  hShiftND_dn->Draw("hist same");
233  hShiftFD_up->Draw("hist same");
234  hShiftFD_dn->Draw("hist same");
235 
236  TGraph* grShade = ShadeBetweenHistograms(hShiftND_up, hShiftFD_up);
237  grShade->SetFillColorAlpha(kRed-9, 0.35);
238  grShade->SetLineColor(kRed-9);
239  grShade->Draw("F");
240 
241  TGraph* grShade2 = ShadeBetweenHistograms(hShiftND_dn, hShiftFD_dn);
242  grShade2->SetFillColorAlpha(kRed-9, 0.35);
243  grShade2->SetLineColor(kRed-9);
244  grShade2->Draw("F");
245 
246  TGraph* g1 = new TGraph;
247  g1->SetPoint(0, -1, 1);
248  g1->SetPoint(1, 10, 1);
249 
250  g1->SetLineWidth(1);
251  g1->SetLineColor(kBlack);
252  g1->SetLineStyle(2);
253 
254  g1->Draw("L");
255 
256  leg->Clear();
257  leg->AddEntry(hShiftND_up, TString("#pm#sigma extrap. ND shift in " + systs[k]->LatexName()),"l");
258  leg->AddEntry(hShiftFD_up, TString("#pm#sigma FD shift " + systs[k]->LatexName()),"l");
259  leg->SetFillColorAlpha(10,0.);
260  leg->Draw();
261 
262  Simulation();
263  Quartile(i);
264 
265  c->SaveAs(Form( ("plots/"+fullName+"_ratio_Q%i.root").c_str(), i));
266  c->SaveAs(Form( ("plots/"+fullName+"_ratio_Q%i.pdf").c_str(), i));
267 
268  // And finally the residual difference
269 
270  c->Clear();
271 
272  TH1* hResup = getResidualDiff(hShiftND_up, hShiftFD_up);
273 
274  TH1* hResdn = getResidualDiff(hShiftND_dn, hShiftFD_dn);
275 
276  hResup->GetYaxis()->SetTitle("Residual difference (%)");
277  hResup->GetYaxis()->SetRangeUser( -10, 10 );
278  //if(k == 0 && (kNumSysts <= 4)) hResup->GetYaxis()->SetRangeUser(-50, 50);
279 
280  hResup->SetLineColor(kOrange+1);
281  hResdn->SetLineColor(kViolet+1);
282 
283  hResup->Draw("hist");
284  hResdn->Draw("hist same");
285 
286  leg->Clear();
287  leg->AddEntry(hResup, TString("+ #sigma shift in " + systs[k]->LatexName() + " FD minus ND"),"l");
288  leg->AddEntry(hResdn, TString("- #sigma shift in " + systs[k]->LatexName() + " FD minus ND"),"l");
289  leg->SetFillColorAlpha(10,0.);
290  leg->Draw();
291 
292  TGraph* g0 = new TGraph;
293  g0->SetPoint(0, -1, 0);
294  g0->SetPoint(1, 10, 0);
295 
296  g0->SetLineWidth(1);
297  g0->SetLineColor(kBlack);
298  g0->SetLineStyle(2);
299 
300  g0->Draw("L");
301 
302  Simulation();
303  Quartile(i);
304 
305  c->SaveAs(Form( ("plots/"+fullName+"_residual_Q%i.root").c_str(), i));
306  c->SaveAs(Form( ("plots/"+fullName+"_residual_Q%i.pdf").c_str(), i));
307 
308  c->Close();
309  } // end loop over systematics
310 
311  } // end loop over quantiles
312 
313  std::cout << "All plots finished!" << std::endl;
314 
315  fPred->Close();
316 
317 } // 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
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 kAna2017POT
Definition: Exposures.h:174
TH1 * getResidualDiff(TH1 *h1, TH1 *h2)
osc::OscCalcDumb calc
Double_t ymax
Definition: plot.C:25
void drawSystsShiftingNDdata_updatedAna()
Spectrum Predict(osc::IOscCalc *calc) const override
TGraph * ShadeBetweenHistograms(TH1 *hmin, TH1 *hmax)
Definition: Plots.cxx:1234
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
const unsigned int kNumQuantiles
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
Spectrum Predict(osc::IOscCalc *calc) const override
TH1F * h1
enum BeamMode kViolet
void Quartile(int q)
Prediction that just uses FD MC, with no extrapolation.
Float_t e
Definition: plot.C:35
enum BeamMode string