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