plot_ND_DataMC_energybinning.C
Go to the documentation of this file.
7 #include "CAFAna/Core/Binning.h"
8 #include "CAFAna/Core/Spectrum.h"
9 #include "CAFAna/Core/Var.h"
11 #include "CAFAna/Vars/HistAxes.h"
14 #include "Utilities/rootlogon.C"
15 
16 #include "TCanvas.h"
17 #include "TGraph.h"
18 #include "TGaxis.h"
19 #include "TLegend.h"
20 #include "TLegendEntry.h"
21 #include "TFile.h"
22 #include "TSystem.h"
23 
24 
25 #include <iostream>
26 #include <string>
27 #include <sstream>
28 #include <fstream>
29 #include <iomanip>
30 
31 using namespace ana;
32 
33 TString sNumu_Nue = "nue";
34 std::vector<TString> hcNames = {"fhc", "rhc"}; const unsigned int nHC = hcNames.size();
35 TString inDir = "/nova/ana/users/essmith/ND_DataMC2020/";
36 //TString inDir = "/nova/app/users/essmith/plotlists_dev/3FlavorAna/Ana2020/ND_DataMC/";
37 std::string outDir = "plots_energybinning/";
38 
39 void SetAxis(TH1* h, bool visible, int color, int style = 1, bool scale = true);
40 TString POTstr(double pot);
41 double DerivePowScale(TH1* h);
42 TH1* MakeRatio(TH1* num, TH1* denom, int color=1);
43 void OneLine();
44 void Zero();
45 void Legend(TH1* h_to_clone);
46 void RatioLegend();
47 void NeutrinoMode(TString mode);
48 
50 
51  std::vector<HistDef> defs;
52  AddHistDefRecoND(defs);
53  AddHistDefSlice(defs);
54  AddHistDefNueND(defs);
55  AddHistDefNueShower(defs);
56  AddHistDefNuePID(defs);
58  AddHistDefProngCVN(defs);
59  AddHistDefNueEnergy(defs);
60  unsigned int nVars = defs.size();
61  std::vector<std::string> varNames;
62  for(unsigned int var = 0; var<nVars; ++var){
63  varNames.emplace_back(defs[var].name);
64  }
65 
66  std::vector<std::string> selNames = {"Presel", "CVNeLow", "CVNeHigh"};
67  unsigned int nSels = selNames.size();
68  std::vector<std::string> compNames = {"nue", "numu", "nc", "antinue", "antinumu"};
69  unsigned int nComps = compNames.size();
70 
71  std::unique_ptr<Spectrum> specData;
72  double datapot[nHC][nVars][nSels];
73  TH1* hData[nHC][nVars][nSels];
74  TH1* hMC[nHC][nVars][nSels][nComps];
75  TH1* hDecomp[nHC][nVars][nSels][nComps];
76 
77  unique_ptr<Spectrum> tempspec;
78 
79  for(unsigned int HC = 0; HC<nHC; ++HC){
80  std::cout<<"Loading files for "<<hcNames[HC]<<std::endl;
81  TString nueDataName = inDir+hcNames[HC]+"_nue_data.root";
82  TString nueMCName = inDir+hcNames[HC]+"_nue_mc.root";
83  TString nueMCDecompName = inDir+hcNames[HC]+"_nue_mc_decomp.root";
84 
85  TFile* nueDataFile = TFile::Open(nueDataName);
86  TFile* nueMCFile = TFile::Open(nueMCName);
87  TFile* nueMCDecompFile = TFile::Open(nueMCDecompName);
88 
89  for(unsigned int var = 0; var<nVars; ++var){
90  for(unsigned int sel = 0; sel<nSels; ++sel){
91  TString name = varNames[var]+"_"+selNames[sel];
92  if(!name.Contains("CVN2d_CVNeLow")) continue;
93 
94  specData = Spectrum::LoadFrom(nueDataFile, name);
95  datapot[nHC][nVars][nSels] = specData->POT();
96  hData[HC][var][sel] = specData->ToTH1(datapot[nHC][nVars][nSels]);
97  for(unsigned int comp = 0; comp<nComps; ++comp){
98  tempspec = Spectrum::LoadFrom(nueMCFile, name+"_"+compNames[comp]);
99  hMC[HC][var][sel][comp] = tempspec->ToTH1(datapot[nHC][nVars][nSels]);
100 
101  tempspec = Spectrum::LoadFrom(nueMCDecompFile, name+"_"+compNames[comp]);
102  hDecomp[HC][var][sel][comp] = tempspec->ToTH1(datapot[nHC][nVars][nSels]);
103  }
104  }
105  }
106 
107  // all up the components to get total MC
108  for(unsigned int var = 0; var<nVars; ++var){
109  for(unsigned int sel = 0; sel<nSels; ++sel){
110  TString name = varNames[var]+"_"+selNames[sel];
111  if(!name.Contains("CVN2d_CVNeLow")) continue;
112  TString forfile = "_"+sNumu_Nue+"_"+hcNames[HC];
113 
114  TH1* totalMC = (TH1*)hMC[HC][var][sel][0]->Clone("totalMC");
115  TH1* totalDecomp = (TH1*)hDecomp[HC][var][sel][0]->Clone("totalDecomp");
116 
117  for(unsigned int comp = 1; comp<nComps; ++comp){
118  totalMC->Add(hMC[HC][var][sel][comp]);
119  totalDecomp->Add(hDecomp[HC][var][sel][comp]);
120  }
121 
122  TH1* rMC;
123  TH1* rDecomp;
124  // histograms for ratios
125  rMC = MakeRatio(hData[HC][var][sel], totalMC, 17);
126  rDecomp = MakeRatio(hData[HC][var][sel], totalDecomp);
127 
128  // Lets start drawing
129  // Split plot
130  TCanvas* c = new TCanvas(name+forfile, name+forfile);
131  double ysplit=0.35;
132  TPad* p1 =new TPad("p1", "", 0,0,1,1);
133  TPad* p2 =new TPad("p2","",0,0,1,1);
134  p1->SetBottomMargin(ysplit);
135  p2->SetTopMargin(1-ysplit);
136  p1->SetFillStyle(0);
137  p2->SetFillStyle(0);
138 
139  p1->cd();
140 
141 
142  std::vector <int> colors = {kMagenta,
143  kGreen+2,
144  kAzure,
145  kMagenta-2,
146  kGreen-3,
147  kRed};
148  std::vector<std::string> labels = {"#nu_{e} CC",
149  "#nu_{#mu} CC",
150  "NC",
151  "#bar{#nu}_{e} CC",
152  "#bar{#nu}_{#mu} CC"};
153 
154  totalDecomp->GetYaxis()->SetTitle(POTstr(datapot[nHC][nVars][nSels]));
155  totalDecomp->GetYaxis()->CenterTitle();
156 
157  SetAxis(totalDecomp, true, colors[nComps]);
158  double maxY = totalDecomp->GetMaximum();
159  totalDecomp->GetYaxis()->SetRangeUser(0, 1.30*maxY);
160  totalDecomp->Draw("hist");
161  Legend(totalDecomp);
162 
164  Nue2018ThreeBinLabels(1.42*maxY, 0.06, kGray+3, true);
165 
166 
167  SetAxis(totalMC, false, colors[nComps], 7);
168  totalMC->Draw("hist same");
169  for(unsigned int comp = 0; comp<nComps; comp++){
170  SetAxis(hDecomp[HC][var][sel][comp], false, colors[comp]);
171  hDecomp[HC][var][sel][comp]->Draw("hist same");
172  SetAxis(hMC[HC][var][sel][comp], false, colors[comp], 7);
173  hMC[HC][var][sel][comp]->Draw("hist same");
174  }
175  hData[HC][var][sel]->SetMarkerStyle(kFullCircle);
176  SetAxis(hData[HC][var][sel], false, 1);
177  hData[HC][var][sel]->Draw("ep same");
178 
179  Nue2018ThreeBinAxis(totalDecomp, false, true, true);
180 
181  totalDecomp->GetXaxis()->SetLabelSize(0.);
182  totalDecomp->GetXaxis()->SetTitleSize(0.);
183  totalDecomp->GetYaxis()->SetTitleSize(0.75*totalDecomp->GetYaxis()->GetTitleSize());
184  totalDecomp->GetYaxis()->SetRangeUser(0.001,totalDecomp->GetMaximum());
185 
186  // The empty bins of the histograms still get drawn at y=0. Try to
187  // draw over these colored lines with black
188  Zero();
189 
190  // draw the ratio
191  p2->cd();
192 
193  SetAxis(rDecomp, true, 1, 1, false);
194  rDecomp->Draw("ep");
195 
196  SetAxis(rMC, false, 17, 1, false);
197  rMC->Draw("ep same");
198  OneLine();
200  Nue2018ThreeBinAxis(rDecomp, true, true, true);
201 
202  rDecomp->GetYaxis()->SetTitleSize(0.75*rDecomp->GetYaxis()->GetTitleSize());
203  rDecomp->GetXaxis()->SetLabelSize(0.75*rDecomp->GetYaxis()->GetTitleSize());
204  rMC->GetYaxis()->SetRangeUser(.5,1.5);
205  rDecomp->GetYaxis()->SetRangeUser(.5,1.5);
206  RatioLegend();
207  Preliminary();
208  NeutrinoMode(hcNames[HC]);
209  p1->Draw();
210  p2->Draw("same");
211 
212  gPad->Print(outDir+name+forfile+".pdf");
213  }
214  }
215  nueDataFile->Close();
216  nueMCFile->Close();
217  nueMCDecompFile->Close();
218 
219  }
220 
221  return;
222 }
223 
224 void SetAxis(TH1* h, bool visible, int color, int style, bool scale) {
225  if(scale) h->Scale(.001);
226  h->SetLineStyle(style);
227  h->SetLineColor(color);
228  if(!visible){
229  h->SetTickLength(0.);
230  h->SetLabelSize(0.);
231  }
232  double norm = DerivePowScale(h);
233 }
234 
235 TString POTstr(double pot) {
236  return TString(TString::Format("10^{3} Events / %.3g#times10^{20} POT", pot/1e20));
237 }
238 
239 double DerivePowScale(TH1* h) {
240  // "Next smallest power of 1000 than maximum bin" - from plot_validation_datamc.py
241  return 3*(int(log10(h->GetMaximum())))/3;
242 }
243 
244 TH1* MakeRatio(TH1* num, TH1* denom, int color) {
245  TH1* ret = (TH1*)num->Clone();
246  ret->Divide(denom);
247 
248  ret->SetMarkerStyle(kFullCircle);
249  ret->SetMarkerColor(color);
250  ret->SetLineColor(color);
251  ret->GetYaxis()->SetRangeUser(.5, 1.5);
252  ret->GetYaxis()->SetTitle("Data / MC");
253  ret->GetXaxis()->SetNdivisions(1900, kFALSE);
254  // ret->GetXaxis()->SetTickLength(1.);
255  return ret;
256 }
257 
258 void Zero() {
259  TGraph* g = new TGraph();
260  g->SetPoint(0, -1e10,0.003);
261  g->SetPoint(1, 1e10, 0.003);
262  g->SetLineWidth(2);
263  g->SetLineStyle(1);
264  g->Draw("same");
265 
266 }
267 
268 void OneLine(){
269  TGraph* g = new TGraph();
270  g->SetPoint(0, -1e10,1);
271  g->SetPoint(1, 1e10, 1);
272  g->SetLineWidth(2);
273  g->SetLineStyle(7);
274  g->Draw("same");
275 }
276 
277 void Legend(TH1* h_to_clone){
278  std::vector <std::pair<int, TString>> entries = {{kBlack, "ND data"},
279  {kRed, "Total MC"},
280  {kAzure, "NC"},
281  {kGreen+2, "#nu_{#mu} CC"},
282  {kMagenta, "#nu_{e} CC"},
283  {kGreen-3, "#bar{#nu}_{#mu} CC"},
284  {kMagenta-2, "#bar{#nu}_{e} CC"},
285  {kBlack, "Uncorr. MC"},
286  {17, "Uncorr. MC"},
287  {1, "Decomp MC"}};
288 
289 
290  TLegend* leg = new TLegend(.82, .5, .99, .92);
291 
292  leg->SetFillColor(kWhite);
293  leg->SetBorderSize(1);
294  leg->SetMargin(.22);
295 
296  TH1* h = (TH1F*)h_to_clone->Clone();
297  leg->SetTextSize(0.5*h->GetXaxis()->GetTitleSize());
298  h->SetLineStyle(1);
299  h->SetLineColor(14);
300  h->SetMarkerStyle(kFullCircle);
301 
302  leg->AddEntry(h, entries[0].second, "lep");
303  for(int i = 1; i < (int) entries.size()-2; i++){
304  h->SetLineColor(entries[i].first);
305  if(i == 7) h->SetLineStyle(7); // Uncorr MC
306  else h->SetLineStyle(1);
307  if(i == 8){
308  h->SetMarkerColor(entries[i].first); // Ratio to Uncorr MC
309  h->SetMarkerStyle(kFullCircle);
310  leg->AddEntry(h->Clone(), entries[i].second, "lep");
311  continue;
312  }
313  if(i == 9) {
314  h->SetMarkerColor(entries[i].first);
315  h->SetMarkerStyle(kFullCircle);
316  leg->AddEntry(h->Clone(), entries[i].second, "lep");
317  continue;
318  }
319  leg->AddEntry(h->Clone(), entries[i].second, "l");
320  }
321  leg->SetTextSize(1.3*leg->GetTextSize());
322  leg->Draw();
323 
324 }
325 
326 void RatioLegend(){
327  TLegend* leg = new TLegend(.12, .13, .47, .18);
328  // leg->SetFillStyle(0);
329  leg->SetNColumns(2);
330  leg->SetBorderSize(1);
331  TH1* h = new TH1F("","", 1, 0, 1);
332  h->SetMarkerColor(17);
333  h->SetLineColor(17);
334  h->SetMarkerStyle(kFullCircle);
335  leg->AddEntry(h->Clone(), "Uncorr. MC", "lep");
336  h->SetMarkerColor(1);
337  h->SetLineColor(1);
338  leg->AddEntry(h->Clone(), "Decomp MC", "lep");
339  leg->Draw();
340 }
341 
342 void NeutrinoMode(TString mode){
343  TString txt = mode == "fhc"? "Neutrino Mode" : "Antineutrino Mode";
344  TLatex* lb = new TLatex(.1, .95, txt);
345  lb->SetNDC();
346  lb->SetTextAlign(13);
347  lb->SetTextColor(kGray+3);
348  lb->SetTextSize(2/45.);
349  lb->Draw();
350 
351 }
void Nue2018ThreeBinAxis(THStack *axes, bool drawLabels, bool merged, bool coreOnly)
const XML_Char * name
Definition: expat.h:151
void Legend(TH1 *h_to_clone)
enum BeamMode kRed
TH1 * MakeRatio(TH1 *num, TH1 *denom, int color=1)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void AddHistDefNueShower(std::vector< HistDef > &hd)
void SetAxis(TH1 *h, bool visible, int color, int style=1, bool scale=true)
void plot_ND_DataMC_energybinning()
void AddHistDefSlice(std::vector< HistDef > &hd)
double maxY
Definition: plot_hist.C:10
void AddHistDefNueSelectionExtras(std::vector< HistDef > &hd)
double DerivePowScale(TH1 *h)
std::string outDir
Double_t scale
Definition: plot.C:25
int colors[6]
Definition: tools.h:1
void Nue2018ThreeBinDivisions(bool coreOnly, const int color, const int style)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
void NeutrinoMode(TString mode)
const unsigned int nHC
#define pot
void AddHistDefRecoND(std::vector< HistDef > &hd)
const HistDef defs[kNumVars]
Definition: vars.h:15
void AddHistDefProngCVN(std::vector< HistDef > &hd)
void Preliminary()
void AddHistDefNuePID(std::vector< HistDef > &hd)
OStream cout
Definition: OStream.cxx:6
void RatioLegend()
Float_t norm
T log10(T number)
Definition: d0nt_math.hpp:120
int num
Definition: f2_nu.C:119
void AddHistDefNueND(std::vector< HistDef > &hd)
void Nue2018ThreeBinLabels(const double yNDC, const double textsize, const int color, const bool nd)
TString POTstr(double pot)
TString sNumu_Nue
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
const std::string selNames[kNumSels]
Definition: vars.h:46
void AddHistDefNueEnergy(std::vector< HistDef > &hd)
enum BeamMode kGreen
std::vector< TString > hcNames
enum BeamMode string