drawSystsShiftingNDdata.C
Go to the documentation of this file.
1 #include "CAFAna/Cuts/Cuts.h"
4 #include "CAFAna/Cuts/NumuCuts.h"
6 #include "CAFAna/Core/Loaders.h"
8 #include "CAFAna/Core/Spectrum.h"
10 #include "CAFAna/Vars/Vars.h"
11 #include "CAFAna/Vars/NumuVars.h"
14 #include "CAFAna/Analysis/Calcs.h"
16 #include "CAFAna/Analysis/Plots.h"
17 
26 #include "CAFAna/Systs/NumuSysts.h"
27 #include "CAFAna/Vars/FitVars.h"
28 
30 
31 #include "TStyle.h"
32 #include "TCanvas.h"
33 #include "TGraph.h"
34 #include "TGraphAsymmErrors.h"
35 #include "TH1.h"
36 #include "TH2.h"
37 #include "TLegend.h"
38 #include "TLatex.h"
39 #include "TArrow.h"
40 
41 using namespace ana;
42 
43 const unsigned int kNumQuantiles = 4;
44 const std::vector<std::string> kOutputExts = {".png", ".eps", ".pdf", ".root"};
45 
46 
47 const std::map<const ISyst*, std::string> kSystRename
48 {
49  { &kMAQEGenieReducedSyst2018, "M_{A}^{QE}" },
50  { &kRPACCQEEnhSyst2018, "RPA CCQE (enh)" },
51  { &kRPACCQESuppSyst2018, "RPA CCQE (supp)" },
52  { &kRPARESSyst2018, "RPA RES" },
53  { &kMECShapeSyst2018Nu, "MEC 4-mom. xfer shape" },
54  // yeah -- use the same text as for nu.
55  // these plots are illustrative only, and we won't make the WS versions.
56  { &kMECShapeSyst2018AntiNu, "MEC 4-mom. xfer shape" },
57  { &kMECEnuShapeSyst2018Nu, "MEC E_{#nu} shape" },
58  { &kMECEnuShapeSyst2018AntiNu, "MEC E_{#nu} shape" }, // etc.
59  { &kMECInitStateNPFracSyst2018Nu, "MEC it'l. state np fraction" },
60  { &kMECInitStateNPFracSyst2018AntiNu, "MEC it'l. state np fraction" },
61 
62  { GetGenieKnobSyst(rwgt::fReweightMaCCRES), "M_{A}^{RES}" },
63  { GetGenieKnobSyst(rwgt::fReweightMvCCRES), "M_{V}^{RES}" },
64 };
65 
66 // Put a "NOvA Simulation" tag in the corner .
67 // this was clashing with the .rootlogon...
68 namespace numu {
69  void Simulation() {
70  TLatex *prelim = new TLatex(.9, .95, "NOvA Simulation");
71  prelim->SetTextColor(kGray + 1);
72  prelim->SetNDC();
73  prelim->SetTextSize(2 / 30.);
74  prelim->SetTextAlign(32);
75  prelim->Draw();
76  }
77 }
78 
80 {
81  auto beamLabel = new TLatex(.15, .92, (beam == Loaders::kFHC ? "Neutrino" : "Antineutrino") + TString(" beam"));
82  beamLabel->SetNDC();
83  beamLabel->SetTextSize(0.04);
84  beamLabel->Draw();
85 }
86 
87 void WriteQuartile(unsigned int q)
88 {
89  if (q >= kNumQuantiles)
90  return;
91 
92  TLatex* quart = new TLatex(.17, .75, Form("Quartile %i",q+1));
93  quart->SetNDC();
94  quart->SetTextSize(0.05);
95  quart->Draw();
96 }
97 
98 void WriteSystName(const ISyst * syst, bool top=false)
99 {
100  std::string name = (kSystRename.find(syst) != kSystRename.end()) ? kSystRename.at(syst) : syst->LatexName();
101  TLatex * text = new TLatex(0.17, top ? 0.83 : 0.2, name.c_str());
102  text->SetNDC();
103  text->SetTextSize(1.5/30.);
104  text->Draw();
105 }
106 
107 TH1* getResidualDiff(TH1* h1, TH1* h2)
108 {
109  TH1* h = (TH1*) h1->Clone("h");
110  h->Reset();
111 
112  for(int k = 1; k <= h1->GetNbinsX(); k++)
113  {
114  double val1 = (double) h1->GetBinContent(k);
115  double val2 = (double) h2->GetBinContent(k);
116 
117  if(val1 * val2 > 0)
118  {
119  h->SetBinContent(k, (val1-val2) * 200./(val1 + val2));
120  }
121  }
122  return h;
123 }
124 
125 TDirectory * GetDir(TFile * f, const TString & dirName)
126 {
127 // std::cout << "Loading directory: " << dirName << std::endl;
128  return f->GetDirectory(dirName);
129 }
130 
131 TH1 * GetQuantilePredictionHist(unsigned int quantNum, IPrediction * preds[], osc::IOscCalculator * calc, double POT)
132 {
133  if (quantNum < kNumQuantiles)
134  return preds[quantNum]->Predict(calc).ToTH1(POT);
135 
136  // rather than taking the prediction that's extrapolated using all 4 quartiles to get the to get the sum,
137  // sum up the predictions from each individually extrapolated quartile to get the total.
138  TH1 * h = preds[0]->Predict(calc).ToTH1(POT);
139  for (unsigned int i = 1; i < kNumQuantiles; i++)
140  h->Add(preds[i]->Predict(calc).ToTH1(POT));
141  return h;
142 }
143 
145 {
146  Loaders::FluxType hornCurr;
147  if (hornCurrStr == "FHC")
148  hornCurr = Loaders::kFHC;
149  else if (hornCurrStr == "RHC")
150  hornCurr = Loaders::kRHC;
151  else
152  {
153  std::cerr << "Unrecognized horn current string: '" << hornCurrStr << "'. Options are FHC, RHC." << std::endl;
154  return;
155  }
156  std::transform(hornCurrStr.begin(), hornCurrStr.end(), hornCurrStr.begin(), ::tolower);
157 
158  std::cout << "Loading predictions...";
159  std::cout << std::flush;
160 
161  gStyle->SetStatFormat("6.3g"); // Significant figures. Fix if this is not satisfactory
162  gStyle->SetPaintTextFormat("6.3g");
163  gStyle->SetFitFormat("6.3g");
164 
165  double POT = hornCurr == Loaders::kFHC ? kAna2018FHCPOT : kAna2018RHCPOT;
166 
169 
170  calc.SetTh23(TMath::ASin(sqrt(0.58))); // 2018 joint fit result
171  calc.SetDmsq32(2.51e-3); // Same
172 
173  std::vector<const ISyst*> systs = getAllDirectNumuSysts2018();
174 
175  const unsigned int kNumSysts = systs.size();
176 
177  IPrediction* predictionND[kNumSysts*2][kNumQuantiles];
178  IPrediction* predictionFD[kNumSysts*2][kNumQuantiles];
179  IPrediction* nominalPred[kNumQuantiles];
180 
181  TFile* fPred = new TFile( predFile.c_str() );
182 
183  std::string predNameND = "ccE_shiftNDdata";
184  std::string predNameFD = "ccE_noExtrapFD";
185 
186  for(unsigned int i = 0; i < kNumQuantiles; i++)
187  {
188  nominalPred[i] = LoadFrom<PredictionNoExtrap>(GetDir(fPred, Form("nominalPred_%i",i))).release();
189 
190  for(unsigned int k = 0; k < kNumSysts; k++)
191  {
192  predictionND[2*k][i] = LoadFrom<PredictionInterp>(GetDir(fPred, Form( (predNameND+"_Pos_"+systs[k]->ShortName()+"_%i").c_str(), i) )).release();
193 
194  predictionND[2*k+1][i] = LoadFrom<PredictionInterp>(GetDir(fPred, Form( (predNameND+"_Neg_"+systs[k]->ShortName()+"_%i").c_str(), i) )).release();
195 
196  // Now FD (no extrap)
197 
198  predictionFD[2*k][i] = LoadFrom<PredictionNoExtrap>(GetDir(fPred, Form( (predNameFD+"_Pos_"+systs[k]->ShortName()+"_%i").c_str(), i) )).release();
199 
200  predictionFD[2*k+1][i] = LoadFrom<PredictionNoExtrap>(GetDir(fPred, Form( (predNameFD+"_Neg_"+systs[k]->ShortName()+"_%i").c_str(), i) )).release();
201  } // end loop ver systs
202  } // end loop over quantiles
203 
204  std::cout << " done." << std::endl;
205 
206  std::cout << "Loading went OK. Plotting!" << std::endl;
207 
208  for(unsigned int i = 0; i < kNumQuantiles +1; i++)
209  {
210  TH1* hNominal = GetQuantilePredictionHist(i, nominalPred, &calc, POT);
211  hNominal->Scale(0.1, "width");
212 
213  hNominal->SetLineColor(kBlack);
214  double ymax = hNominal->GetBinContent(hNominal->GetMaximumBin()) * 2.0;
215  hNominal->GetYaxis()->SetRangeUser(0, 15);
216  hNominal->GetYaxis()->SetTitle("Events / (0.1 GeV)");
217  hNominal->GetYaxis()->CenterTitle();
218  hNominal->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
219  hNominal->GetXaxis()->CenterTitle();
220 
221  for (unsigned int k = 0; k < kNumSysts; k++)
222  {
223  std::string fullName = hornCurrStr+"_ccE_"+systs[k]->ShortName();
224 
225  TCanvas* c = new TCanvas("c",(fullName).c_str());
226 
227  c->SetBottomMargin(0.15);
228  c->SetLeftMargin(0.15);
229 
230  TH1* hShiftND_up = GetQuantilePredictionHist(i, predictionND[2*k], &calc, POT);
231  hShiftND_up->Scale(0.1, "width");
232  TH1* hShiftFD_up = GetQuantilePredictionHist(i, predictionFD[2*k], &calc, POT);
233  hShiftFD_up->Scale(0.1, "width");
234 
235  TH1* hShiftND_dn = GetQuantilePredictionHist(i, predictionND[2*k+1], &calc, POT);
236  hShiftND_dn->Scale(0.1, "width");
237  TH1* hShiftFD_dn = GetQuantilePredictionHist(i, predictionFD[2*k+1], &calc, POT);
238  hShiftFD_dn->Scale(0.1, "width");
239 
240  hShiftND_up->SetLineColor(kOrange+1);
241  hShiftND_up->SetLineStyle(2);
242 
243  hShiftND_dn->SetLineColor(kViolet+1);
244  hShiftND_dn->SetLineStyle(2);
245 
246  hShiftFD_up->SetMarkerStyle(20);
247  hShiftFD_up->SetMarkerColor(kOrange+1);
248  hShiftFD_up->SetLineWidth(0);
249  hShiftFD_up->SetLineColor(0);
250  hShiftFD_dn->SetMarkerStyle(21);
251  hShiftFD_dn->SetMarkerColor(kViolet+1);
252  hShiftFD_dn->SetLineWidth(0);
253  hShiftFD_dn->SetLineColor(0);
254 
255  hNominal->Draw("hist");
256  hShiftND_up->Draw("hist same");
257  hShiftND_dn->Draw("hist same");
258  hShiftFD_up->Draw("EX0 same");
259  hShiftFD_dn->Draw("EX0 same");
260 
261  auto leg = std::make_unique<TLegend>(0.55,0.55,0.95,0.8);
262  leg->AddEntry(hNominal, "Nominal with FD/ND","l");
263  leg->AddEntry(hShiftND_up, "+ 1#sigma in ND","l");
264  leg->AddEntry(hShiftFD_up, "+ 1#sigma in FD","p");
265  leg->AddEntry(hShiftND_dn, "- 1#sigma in ND","l");
266  leg->AddEntry(hShiftFD_dn, "- 1#sigma in FD","p");
267  leg->SetFillStyle(0);
268  leg->Draw();
269 
270  WriteSystName(systs[k], true);
272  WriteQuartile(i);
273  WriteBeam(hornCurr);
274 
275  auto qstr = (i < kNumQuantiles) ? std::to_string(i) : "all";
276  util::SaveObj(c, fullName + "_Q" + qstr, outDir, kOutputExts, true);
277 
278  // Now the ratios
279 
280  c->Clear();
281 
282  hShiftND_up->Divide(hNominal);
283  hShiftND_dn->Divide(hNominal);
284 
285  hShiftFD_up->Divide(hNominal);
286  hShiftFD_dn->Divide(hNominal);
287 
288  hShiftND_up->GetYaxis()->SetRangeUser(0.8,1.2);
289  //if(k == 0 && (kNumSysts <= 4)) hShiftND_up->GetYaxis()->SetRangeUser(0.5, 2.0);
290  hShiftND_up->GetYaxis()->SetTitle("Ratio to nominal MC");
291  hShiftND_up->GetXaxis()->SetTitle("Reconstructed neutrino energy (GeV)");
292  hShiftND_up->GetYaxis()->CenterTitle();
293  hShiftND_up->GetXaxis()->CenterTitle();
294 
295  hShiftND_up->SetMarkerColor(kOrange+1);
296  hShiftND_up->SetLineWidth(2);
297  hShiftND_up->SetLineStyle(1);
298  hShiftND_up->SetLineColor(kOrange+1);
299  hShiftND_dn->SetMarkerColor(kOrange+1);
300  hShiftND_dn->SetLineWidth(2);
301  hShiftND_dn->SetLineStyle(1);
302  hShiftND_dn->SetLineColor(kOrange+1);
303 
304  hShiftFD_up->SetMarkerColor(kViolet+1);
305  hShiftFD_up->SetLineWidth(2);
306  hShiftFD_up->SetLineStyle(2);
307  hShiftFD_up->SetLineColor(kViolet+1);
308  hShiftFD_dn->SetMarkerColor(kViolet+1);
309  hShiftFD_dn->SetLineWidth(2);
310  hShiftFD_dn->SetLineStyle(2);
311  hShiftFD_dn->SetLineColor(kViolet+1);
312 
313  hShiftND_up->Draw("hist");
314  hShiftND_dn->Draw("hist same");
315  hShiftFD_up->Draw("hist same");
316  hShiftFD_dn->Draw("hist same");
317 
318  TGraph* grShade = ShadeBetweenHistograms(hShiftND_up, hShiftFD_up);
319  grShade->SetFillColorAlpha(kRed-9, 0.35);
320  grShade->SetLineColor(kRed-9);
321  grShade->Draw("F");
322 
323  TGraph* grShade2 = ShadeBetweenHistograms(hShiftND_dn, hShiftFD_dn);
324  grShade2->SetFillColorAlpha(kRed-9, 0.35);
325  grShade2->SetLineColor(kRed-9);
326  grShade2->Draw("F");
327 
328  TGraph* g1 = new TGraph;
329  g1->SetPoint(0, -1, 1);
330  g1->SetPoint(1, 10, 1);
331 
332  g1->SetLineWidth(1);
333  g1->SetLineColor(kBlack);
334  g1->SetLineStyle(2);
335 
336  g1->Draw("L");
337 
338  leg = std::make_unique<TLegend>(0.4,0.7,0.9,0.88);
339  leg->SetFillStyle(0);
340  leg->AddEntry(hShiftND_up, TString("#pm 1#sigma extrap. ND shift"),"l");
341  leg->AddEntry(hShiftFD_up, TString("#pm 1#sigma FD shift"),"l");
342  leg->SetFillColorAlpha(10,0.);
343  leg->Draw();
344 
345  WriteSystName(systs[k]);
347  WriteQuartile(i);
348  WriteBeam(hornCurr);
349 
350  util::SaveObj(c, fullName + "_ratio_Q" + qstr, outDir, kOutputExts, true);
351 
352  // And finally the residual difference
353 
354  c->Clear();
355 
356  TH1* hResup = getResidualDiff(hShiftND_up, hShiftFD_up);
357 
358  TH1* hResdn = getResidualDiff(hShiftND_dn, hShiftFD_dn);
359 
360  hResup->GetYaxis()->SetTitle("Residual difference (%)");
361  hResup->GetYaxis()->SetRangeUser( -10, 10 );
362  //if(k == 0 && (kNumSysts <= 4)) hResup->GetYaxis()->SetRangeUser(-50, 50);
363 
364  hResup->SetLineColor(kOrange+1);
365  hResdn->SetLineColor(kViolet+1);
366 
367  hResup->Draw("hist");
368  hResdn->Draw("hist same");
369 
370  leg->Clear();
371  leg->AddEntry(hResup, TString("+1 #sigma shift in FD minus ND"),"l");
372  leg->AddEntry(hResdn, TString("-1 #sigma shift in FD minus ND"),"l");
373  leg->SetFillColorAlpha(10,0.);
374  leg->Draw();
375 
376  TGraph* g0 = new TGraph;
377  g0->SetPoint(0, -1, 0);
378  g0->SetPoint(1, 10, 0);
379 
380  g0->SetLineWidth(1);
381  g0->SetLineColor(kBlack);
382  g0->SetLineStyle(2);
383 
384  g0->Draw("L");
385 
386  WriteSystName(systs[k]);
388  WriteQuartile(i);
389  WriteBeam(hornCurr);
390 
391  util::SaveObj(c, fullName + "_residual_Q" + qstr, outDir, kOutputExts, true);
392 
393  c->Close();
394  } // end loop over systematics
395 
396  } // end loop over quantiles
397 
398  std::cout << "All plots finished!" << std::endl;
399 
400  fPred->Close();
401 
402 } // End of function
const NOvARwgtSyst kMECShapeSyst2018AntiNu("MECShape2018AntiNu","MEC 2018 (q_{0}, |#vec{q}|) response, antineutrinos", novarwgt::kMECQ0Q3RespSystNubar2018)
Definition: MECSysts.h:12
void drawSystsShiftingNDdata(std::string hornCurrStr, std::string predFile, std::string outDir)
const XML_Char * name
Definition: expat.h:151
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
std::vector< SystGroupDef > systs
Definition: syst_header.h:384
const NOvARwgtSyst kRPACCQEEnhSyst2018("RPAShapeenh2018","RPA shape: higher-Q^{2} enhancement (2018)", novarwgt::kRPACCQEEnhSyst2018)
Definition: RPASysts.h:14
General interface to oscillation calculators.
Definition: FwdDeclare.h:15
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:553
Optimized version of OscCalculatorPMNS.
Definition: StanTypedefs.h:28
void WriteQuartile(unsigned int q)
const NOvARwgtSyst kMAQEGenieReducedSyst2018(genie::rew::GSyst::AsString(genie::rew::EGSyst(rwgt::fReweightMaCCQE))++"_reduced_2018", genie::rew::GSyst::AsString(genie::rew::EGSyst(rwgt::fReweightMaCCQE))++"_reduced_2018", novarwgt::kMAQEGenieReducedSyst2018)
2018 &#39;reduced&#39; M_A^QE shift. See documentation in NOvARwgt
Definition: XSecSysts.h:51
T sqrt(T number)
Definition: d0nt_math.hpp:156
void WriteBeam(Loaders::FluxType beam)
const NOvARwgtSyst kMECEnuShapeSyst2018AntiNu("MECEnuShapeAntiNu","MEC E_{#nu} shape, antineutrinos", novarwgt::kMECEnuShapeSyst2017_NubarOnly)
Definition: MECSysts.h:18
const NOvARwgtSyst kMECShapeSyst2018Nu("MECShape2018Nu","MEC 2018 (q_{0}, |#vec{q}|) response, neutrinos", novarwgt::kMECQ0Q3RespSystNu2018)
Definition: MECSysts.h:11
OStream cerr
Definition: OStream.cxx:7
const std::vector< std::string > kOutputExts
const NOvARwgtSyst kMECInitStateNPFracSyst2018AntiNu("MECInitStateNPFracAntiNu","MEC initial state np fraction, antineutrinos", novarwgt::kMECInitStateNPFracSyst2017_NubarOnly)
Definition: MECSysts.h:21
void ResetOscCalcToDefault(osc::IOscCalculatorAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
void SetTh23(const T &th23) override
const NOvARwgtSyst kMECEnuShapeSyst2018Nu("MECEnuShapeNu","MEC E_{#nu} shape, neutrinos", novarwgt::kMECEnuShapeSyst2017_NuOnly)
Definition: MECSysts.h:17
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
void SetDmsq32(const T &dmsq32) override
const double kAna2018RHCPOT
Definition: Exposures.h:208
std::vector< const ISyst * > getAllDirectNumuSysts2018()
Definition: NumuSysts.cxx:123
Double_t ymax
Definition: plot.C:25
TH1 * getResidualDiff(TH1 *h1, TH1 *h2)
void Simulation()
const NOvARwgtSyst kMECInitStateNPFracSyst2018Nu("MECInitStateNPFracNu","MEC initial state np fraction, neutrinos", novarwgt::kMECInitStateNPFracSyst2017_NuOnly)
Definition: MECSysts.h:20
TGraph * ShadeBetweenHistograms(TH1 *hmin, TH1 *hmax)
Definition: Plots.cxx:907
TLatex * prelim
Definition: Xsec_final.C:133
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
TH1F * h2
Definition: plot.C:45
TH1F * h1
virtual const std::string & LatexName() const final
The name used on plots (ROOT&#39;s TLatex syntax)
Definition: ISyst.h:30
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TH1 * GetQuantilePredictionHist(unsigned int quantNum, IPrediction *preds[], osc::IOscCalculator *calc, double POT)
const std::map< const ISyst *, std::string > kSystRename
string syst
Definition: plotSysts.py:176
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const NOvARwgtSyst kRPACCQESuppSyst2018("RPAShapesupp2018","RPA shape: low-Q^{2} suppression (2018)", novarwgt::kRPACCQESuppSyst2018)
Definition: RPASysts.h:15
void SaveObj(const TObject *obj, const std::string &filenameStub, const std::string &dirName="", const std::vector< std::string > exts={".png",".eps",".root"}, bool silent=false)
Definition: ROOTHelpers.h:21
const double kAna2018FHCPOT
Definition: Exposures.h:207
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void WriteSystName(const ISyst *syst, bool top=false)
TDirectory * GetDir(TFile *f, const TString &dirName)
Float_t e
Definition: plot.C:35
const NOvARwgtSyst * GetGenieKnobSyst(rwgt::ReweightLabel_t knobIdx, std::string altName, std::string altLabel)
Convenience function to get a GENIE knob syst. (Allows using the GENIE knob name & description as the...
Definition: XSecSysts.cxx:117
TGeoVolume * top
Definition: make_fe_box.C:9
const unsigned int kNumQuantiles
h
Definition: demo3.py:41
const NOvARwgtSyst kRPARESSyst2018("RPAShapeRES2018","RPA shape: resonance events (2018)", novarwgt::kRPARESSyst2018)
Definition: RPASysts.h:20