plot_nd_spectra_2018.C
Go to the documentation of this file.
1 #pragma once
2 
8 #include "CAFAna/Core/Spectrum.h"
11 #include "Utilities/rootlogon.C"
12 
13 #include "TFile.h"
14 #include "TH1.h"
15 #include "THStack.h"
16 #include "TCanvas.h"
17 #include "TGaxis.h"
18 #include "TGraph.h"
19 #include "TLatex.h"
20 #include "TLegend.h"
21 #include "TLegendEntry.h"
22 #include "TLine.h"
23 #include "TSystem.h"
24 #include "TPaveLabel.h"
25 #include <fstream>
26 #include <iostream>
27 #include <iomanip>
28 #include <string>
29 #include <cstddef>
30 
31 using namespace ana;
32 
33 void SetAxis(TH1* h, bool visible, int color, int style = 1, bool scale = true);
34 TString POTstr(double pot);
35 double DerivePowScale(TH1* h);
36 TH1* MakeRatio(TH1* num, TH1* denom, int color=1);
37 void One();
38 void Zero();
39 void Legend(TH1* h_to_clone);
40 void RatioLegend();
42 
43 
45 
46  std::string fdir = "/nova/ana/nu_e_ana/Ana2018/DataMC/";
47 
48  TFile* fdata = new TFile((fdir+"20180331_Ana2018_"+mode+"_datamc_decomp_data_validation.root").c_str(), "READ");
49  TFile *fmc = new TFile((fdir+ "20180331_Ana2018_"+mode+"_datamc_decomp_mc_validation.root").c_str(), "READ");
50  TFile* fdecomp = new TFile((fdir + "20180331_Ana2018_"+mode+"_datamc_decomp_mc_decomp_validation.root").c_str(), "READ");
51 
52  if(fmc->IsZombie()) {std::cout << "MC File not found "; return; }
53  if(fdata->IsZombie()) {std::cout << "Data File not found "; return; }
54  if(fdecomp->IsZombie()) {std::cout << "MC Decomp File not found"; return; }
55 
56  // Pulling Data
57  auto spec = Spectrum::LoadFrom(fdata, (TString)"Energy_Ana2018_Bin{group=Cut,cat=1_CVNeLow}");
58  double datapot = spec->POT();
59  TH1* hdata = spec->ToTH1(datapot);
60 
61  std::vector <TH1*> mc;
62  std::vector <TH1*> decomp;
63  std::vector <std::string> comps = {"0_nue", "1_numu", "2_nc", "3_antinue", "4_antinumu"};
64  TH1* tempmc;
65  TH1* tempdecomp;
66  unique_ptr<Spectrum> tempspec;
67  // Load spectra from uncorrected mc and decomp corrected mc files
68  for(std::string comp : comps ) {
69  tempspec = Spectrum::LoadFrom(fmc, (TString::Format("Energy_Ana2018_Bin{group=Cut,cat=1_CVNeLow}{group=Component,cat=%s}",comp.c_str())));
70  tempmc = tempspec->ToTH1(datapot);
71 
72  std::cout << "Loaded MC Spectrum Energy_Ana2018_Bin{group=Cut,cat=1_CVNeLow}{group=Component,cat="+comp+"}" << endl;;
73 
74  tempspec = Spectrum::LoadFrom(fdecomp, (TString::Format("Energy_Ana2018_Bin{group=Cut,cat=1_CVNeLow}{group=Component,cat=%s}",comp.c_str())));
75  tempdecomp = tempspec->ToTH1(datapot);
76 
77  std::cout << "Loaded Decomp Spectrum Energy_Ana2018_Bin{group=Cut,cat=1_CVNeLow}{group=Component,cat="+comp+"}" << endl;
78 
79  mc.push_back(tempmc);
80  decomp.push_back(tempdecomp);
81  }
82 
83  // Add up all of the components to get total MC
84  TH1* totalmc = (TH1*)mc[0]->Clone("totalmc");
85  TH1* totaldecomp = (TH1*)decomp[0]->Clone("totaldecomp");
86  for(int i = 1; i < (int)mc.size(); i++) {
87  totalmc->Add( mc[i] );
88  totaldecomp->Add( decomp[i] );
89  std::cout << "Adding component " + std::to_string(i) << endl;
90  }
91  mc.push_back(totalmc);
92  decomp.push_back(totaldecomp);
93  std::cout << "Added total MC and Decomp" << endl;
94 
95  int tot = mc.size()-1;
96  TH1* rdecomp;
97  TH1* rmc;
98  // make the histograms for the ratios
99  rdecomp = MakeRatio(hdata, decomp[tot]);
100  rmc = MakeRatio(hdata, mc[tot], 17);
101 
102 
103  // Lets start drawing
104  // Split plot
105  TCanvas* c = new TCanvas("c");
106  double ysplit=0.35;
107  TPad* p1 =new TPad("p1", "", 0,0,1,1);
108  TPad* p2 =new TPad("p2","",0,0,1,1);
109  p1->SetBottomMargin(ysplit);
110  p2->SetTopMargin(1-ysplit);
111  p1->SetFillStyle(0);
112  p2->SetFillStyle(0);
113 
114  p1->cd();
115 
116 
117  std::vector <int> colors = {kMagenta,
118  kGreen+2,
119  kAzure,
120  kMagenta-2,
121  kGreen-3,
122  kRed};
123  std::vector<std::string> labels = {"#nu_{e} CC",
124  "#nu_{#mu} CC",
125  "NC",
126  "#bar{#nu}_{e} CC",
127  "#bar{#nu}_{#mu} CC"};
128 
129 
130 
131  decomp[tot]->GetYaxis()->SetTitle(POTstr(datapot));
132  decomp[tot]->GetYaxis()->CenterTitle();
133 
134  SetAxis(decomp[tot], true, colors[tot]);
135  double maxY = decomp[tot]->GetMaximum();
136  decomp[tot]->GetYaxis()->SetRangeUser(0, 1.30*maxY);
137  decomp[tot]->Draw("hist");
138  Legend(decomp[tot]);
139 
140  // Makes it look fancy
142  Nue2018ThreeBinLabels(1.42*maxY, 0.06, kGray+3, true);
143 
144  SetAxis(mc[tot], false, colors[tot], 7);
145  mc[tot]->Draw("hist same");
146  for(int j = 0; j < (int)mc.size()-1 ; j++) {
147  SetAxis(decomp[j], false, colors[j]);
148  decomp[j]->Draw("hist same");
149  SetAxis(mc[j], false, colors[j], 7);
150  mc[j]->Draw("hist same");
151  }
152 
153  hdata->SetMarkerStyle(kFullCircle);
154  SetAxis(hdata, false, 1);
155  hdata->Draw("ep same");
156 
157  // Add fancy axes
158  Nue2018ThreeBinAxis(decomp[tot], false, true, true);
159  decomp[tot]->GetXaxis()->SetLabelSize(0.);
160  decomp[tot]->GetXaxis()->SetTitleSize(0.);
161  decomp[tot]->GetYaxis()->SetTitleSize(0.75*decomp[tot]->GetYaxis()->GetTitleSize());
162  decomp[tot]->GetYaxis()->SetRangeUser(0.001,
163  decomp[tot]->GetMaximum());
164 
165  // The empty bins of the histograms still get drawn at y=0. Try to
166  // draw over these colored lines with black
167  Zero();
168 
169  // Draw the ratio
170  p2->cd();
171 
172  SetAxis(rdecomp, true, 1, 1, false);
173  rdecomp->Draw("ep");
174 
175  SetAxis(rmc, false, 17, 1, false);
176  rmc->Draw("ep same");
177  One();
179  Nue2018ThreeBinAxis(rdecomp, true, true, true);
180 
181  rdecomp->GetYaxis()->SetTitleSize(0.75*rdecomp->GetYaxis()->GetTitleSize());
182  rmc->GetYaxis()->SetRangeUser(.5,1.5);
183  rdecomp->GetYaxis()->SetRangeUser(.5,1.5);
184  RatioLegend();
185  Preliminary();
187 
188  p1->Draw();
189  p2->Draw("same");
190 
191  gPad->Print("Energy_Ana2018_Bin_CVNeLow_split.pdf");
192 }
193 
194 void SetAxis(TH1* h, bool visible, int color, int style, bool scale) {
195  if(scale) h->Scale(.001);
196  h->SetLineStyle(style);
197  h->SetLineColor(color);
198  if(!visible){
199  h->SetTickLength(0.);
200  h->SetLabelSize(0.);
201  }
202  double norm = DerivePowScale(h);
203 }
204 
205 TString POTstr(double pot) {
206  return TString(TString::Format("10^{3} Events / %.3g#times10^{20} POT", pot/1e20));
207 }
208 
209 double DerivePowScale(TH1* h) {
210  // "Next smallest power of 1000 than maximum bin" - from plot_validation_datamc.py
211  return 3*(int(log10(h->GetMaximum())))/3;
212 }
213 
214 TH1* MakeRatio(TH1* num, TH1* denom, int color) {
215  TH1* ret = (TH1*)num->Clone();
216  ret->Divide(denom);
217 
218  ret->SetMarkerStyle(kFullCircle);
219  ret->SetMarkerColor(color);
220  ret->SetLineColor(color);
221  ret->GetYaxis()->SetRangeUser(.5, 1.5);
222  ret->GetYaxis()->SetTitle("Data / MC");
223  ret->GetXaxis()->SetNdivisions(1900, kFALSE);
224  // ret->GetXaxis()->SetTickLength(1.);
225  return ret;
226 }
227 
228 void Zero() {
229  TGraph* g = new TGraph();
230  g->SetPoint(0, -1e10,0.003);
231  g->SetPoint(1, 1e10, 0.003);
232  g->SetLineWidth(2);
233  g->SetLineStyle(1);
234  g->Draw("same");
235 
236 }
237 
238 void One(){
239  TGraph* g = new TGraph();
240  g->SetPoint(0, -1e10,1);
241  g->SetPoint(1, 1e10, 1);
242  g->SetLineWidth(2);
243  g->SetLineStyle(7);
244  g->Draw("same");
245 }
246 
247 void Legend(TH1* h_to_clone){
248  std::vector <std::pair<int, TString>> entries = {{kBlack, "ND data"},
249  {kRed, "Total MC"},
250  {kAzure, "NC"},
251  {kGreen+2, "#nu_{#mu} CC"},
252  {kMagenta, "#nu_{e} CC"},
253  {kGreen-3, "#bar{#nu}_{#mu} CC"},
254  {kMagenta-2, "#bar{#nu}_{e} CC"},
255  {kBlack, "Uncorr. MC"},
256  {17, "Uncorr. MC"},
257  {1, "Decomp MC"}};
258 
259 
260  TLegend* leg = new TLegend(.78, .4, .95, .89);
261 
262  leg->SetFillColor(kWhite);
263  leg->SetBorderSize(1);
264  leg->SetMargin(.22);
265 
266  TH1* h = (TH1F*)h_to_clone->Clone();
267  leg->SetTextSize(0.55*h->GetXaxis()->GetTitleSize());
268  h->SetLineStyle(1);
269  h->SetLineColor(14);
270  h->SetMarkerStyle(kFullCircle);
271 
272  leg->AddEntry(h, entries[0].second, "lep");
273  for(int i = 1; i < (int) entries.size()-2; i++){
274  h->SetLineColor(entries[i].first);
275  if(i == 7) h->SetLineStyle(7); // Uncorr MC
276  else h->SetLineStyle(1);
277  if(i == 8){
278  h->SetMarkerColor(entries[i].first); // Ratio to Uncorr MC
279  h->SetMarkerStyle(kFullCircle);
280  leg->AddEntry(h->Clone(), entries[i].second, "lep");
281  continue;
282  }
283  if(i == 9) {
284  h->SetMarkerColor(entries[i].first);
285  h->SetMarkerStyle(kFullCircle);
286  leg->AddEntry(h->Clone(), entries[i].second, "lep");
287  continue;
288  }
289  leg->AddEntry(h->Clone(), entries[i].second, "l");
290  }
291  leg->SetTextSize(1.3*leg->GetTextSize());
292  leg->Draw();
293 
294 }
295 
296 void RatioLegend(){
297  TLegend* leg = new TLegend(.12, .13, .47, .18);
298  // leg->SetFillStyle(0);
299  leg->SetNColumns(2);
300  leg->SetBorderSize(1);
301  TH1* h = new TH1F("","", 1, 0, 1);
302  h->SetMarkerColor(17);
303  h->SetLineColor(17);
304  h->SetMarkerStyle(kFullCircle);
305  leg->AddEntry(h->Clone(), "Uncorr. MC", "lep");
306  h->SetMarkerColor(1);
307  h->SetLineColor(1);
308  leg->AddEntry(h->Clone(), "Decomp MC", "lep");
309  leg->Draw();
310 }
311 
313  std::string txt = mode == "fhc"? "Neutrino Mode" : "Antineutrino Mode";
314  TLatex* lb = new TLatex(.1, .95, txt.c_str());
315  lb->SetNDC();
316  lb->SetTextAlign(13);
317  lb->SetTextColor(kGray+3);
318  lb->SetTextSize(2/45.);
319  lb->Draw();
320 
321 }
void Nue2018ThreeBinAxis(THStack *axes, bool drawLabels, bool merged, bool coreOnly)
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void Legend(TH1 *h_to_clone)
correl_xv GetYaxis() -> SetDecimals()
double maxY
Definition: plot_hist.C:10
TFile * fmc
Definition: bdt_com.C:14
Double_t scale
Definition: plot.C:25
int colors[6]
Definition: tools.h:1
void Nue2018ThreeBinDivisions(bool coreOnly, const int color, const int style)
void RatioLegend()
void NeutrinoMode(std::string mode)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
void plot_nd_spectra_2018(std::string mode="fhc")
void SetAxis(TH1 *h, bool visible, int color, int style=1, bool scale=true)
#define pot
const double j
Definition: BetheBloch.cxx:29
std::vector< std::string > comps
void Preliminary()
TString POTstr(double pot)
OStream cout
Definition: OStream.cxx:6
Float_t norm
T log10(T number)
Definition: d0nt_math.hpp:120
TH1 * MakeRatio(TH1 *num, TH1 *denom, int color=1)
void Zero()
int num
Definition: f2_nu.C:119
void One()
void Nue2018ThreeBinLabels(const double yNDC, const double textsize, const int color, const bool nd)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
double DerivePowScale(TH1 *h)
enum BeamMode kGreen
enum BeamMode string