neutKEsyst_plot.C
Go to the documentation of this file.
1 #include "TGaxis.h"
2 #include "TCanvas.h"
3 #include "TFile.h"
4 #include "TH1D.h"
5 #include "TLatex.h"
6 #include "TLegend.h"
7 #include "TStyle.h"
8 #include "TSystem.h"
9 #include "TF1.h"
10 #include "TH2.h"
11 #include "TGraph.h"
12 
14 #include "CAFAna/Analysis/Plots.h"
15 #include "CAFAna/Core/Spectrum.h"
17 
18 #include <iostream>
19 #include <sstream>
20 #include <fstream>
21 #include <vector>
22 #include <string>
23 #include <TROOT.h>
24 #include <TStyle.h>
25 
26 using namespace ana;
27 
28 // ****************************** Some very broad things ******************************
29 
30 // --- Define the number of GENIE interaction types I'm looking at....
31 const std::vector<std::string> GENIEStr = { "All" }; const unsigned int nGENIE = 1;
32 //const std::vector<std::string> GENIEStr = { "All", "MEC", "DIS", "Resonance", "Coherent", "QuasiElactic" }; const unsigned int nGENIE = 6;
33 
34 // --- Define my broad cuts
35 std::vector<std::string> ND_CutNames = {"ND_2018full"}; const unsigned int nND_CutNames = 1;
36 unsigned int FCut_ND = nND_CutNames - 1;
37 std::vector<std::string> FD_CutNames = {"FD_2018full"}; const unsigned int nFD_CutNames = 1;
38 unsigned int FCut_FD = nFD_CutNames - 1;
39 
40 // --- Am I running quantile boundaries?
41 std::vector<std::string> QuantNames = { "AllQuant", "Quant1", "Quant2", "Quant3", "Quant4" }; const unsigned int nQuantNames = 5;
42 
43 // --- FHC or RHC?
44 std::vector<std::string> HCNames = {"fhc", "rhc"}; const unsigned int nHC = 2;
45 std::vector<double> POT = {kAna2018FHCPOT, kAna2018RHCPOT};
46 
47 // --- ND or FD? (for labeling)
48 std::vector<std::string> DetNames = {"ND", "FD"}; const unsigned int nDet = 2;
49 
50 // ****************************** Define a struct ******************************
51 struct MySpectra {
52  TH1D* NuE_nom_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
53  TH1D* HadE_nom_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
54  TH1D* NuE_scale_p_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
55  TH1D* HadE_scale_p_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
56  TH1D* NuE_scale_m_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
57  TH1D* HadE_scale_m_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
58  TH1D* NNeu_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
59  TH1D* VisE_sum_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
60  TH1D* VisE_ND [nGENIE][nND_CutNames][nQuantNames][nHC];
61 
62  TH1D* NuE_nom_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
63  TH1D* HadE_nom_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
64  TH1D* NuE_scale_p_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
65  TH1D* HadE_scale_p_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
66  TH1D* NuE_scale_m_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
67  TH1D* HadE_scale_m_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
68  TH1D* NNeu_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
69  TH1D* VisE_sum_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
70  TH1D* VisE_FD [nGENIE][nFD_CutNames][nQuantNames][nHC];
71 
72 };
73 
74 // ****************************** Define functions ******************************
75 TH1D* SpecToHist( TFile* MyF, std::string SpecName, double POT, std::string AxisTitle );
76 void NormByBinWidth(TH1D* hist);
77 void Print1D(TH1D* hist, std::string plotName);
78 void MakeComparisonPlots(std::vector<TH1D*> TheHists, std::string plotName, std::vector<std::string> Labels, bool var);
79 void MakeComparisonPlotsWithRatio(std::vector<TH1D*> TheHists, std::string plotName, std::vector<std::string> Labels, bool var);
80 TH1D* MakeRatio(TH1D* num, TH1D* denom, int Col, std::string FType, bool var);
81 void FitGauss(TH1D* h);
82 // ****************************** Main code block ******************************
84  // --- I don't want to see the print messages ---
85  gErrorIgnoreLevel = kError;
86 
87  // --- NOvA style ---
88  gStyle->SetOptStat(0);
89  gROOT->SetStyle("novaStyle");
90  TGaxis::SetMaxDigits(4);
91 
92  MySpectra sp;
93 
94  std::string BaseCAFFileDir="/nova/ana/users/essmith/neutronKE/syst/";
95  std::string FDName, NDName;
96 
97  // ----- Spectra -----
98  // pulling spectra out of files and converting to histograms
99  // also want to normalize some of these plots by bin width (but not others)
100  // best to do it here instead of in comparison plot code or I'll normalize
101  // the first one a bunch of times by accident
102 
103  // fhc/rhc
104  for (unsigned int HC = 0; HC<nHC; ++HC){
105 
106  FDName = BaseCAFFileDir+"nKEsystBirks_FD"+HCNames[HC]+".root";
107  TFile *FDFile = TFile::Open(FDName.c_str());
108 
109  NDName = BaseCAFFileDir+"nKEsystBirks_ND"+HCNames[HC]+".root";
110  TFile *NDFile = TFile::Open(NDName.c_str());
111 
112  // interaction types
113  for (unsigned int gen=0; gen<nGENIE; ++gen) {
114 
115  // quantiles
116  for (unsigned int quant=0; quant<nQuantNames; ++quant) {
117 
118  // analysis cuts, ND
119  for(unsigned int det=0; det<nND_CutNames; ++det) {
120 
121  std::string Comb = GENIEStr[gen]+"_"+ND_CutNames[det]+"_"+QuantNames[quant];
122 
123  sp.NuE_nom_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NuE_nom_"+Comb,
124  POT[HC], ";Reco #nu E (GeV);Events/0.1GeV");
125  sp.HadE_nom_ND [gen][det][quant][HC] = SpecToHist(NDFile, "HadE_nom_"+Comb,
126  POT[HC], ";Reco Had E (GeV);Events");
127  sp.NuE_scale_p_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NuE_scale_p_"+Comb,
128  POT[HC], ";Reco #nu E (GeV);Events/0.1GeV");
129  sp.HadE_scale_p_ND [gen][det][quant][HC] = SpecToHist(NDFile, "HadE_scale_p_"+Comb,
130  POT[HC], ";Reco Had E (GeV);Events");
131  sp.NuE_scale_m_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NuE_scale_m_"+Comb,
132  POT[HC], ";Reco #nu E (GeV);Events/0.1GeV");
133  sp.HadE_scale_m_ND [gen][det][quant][HC] = SpecToHist(NDFile, "HadE_scale_m_"+Comb,
134  POT[HC], ";Reco Had E (GeV);Events");
135  sp.NNeu_ND [gen][det][quant][HC] = SpecToHist(NDFile, "NNeu_"+Comb,
136  POT[HC], ";Multiplicity");
137  sp.VisE_sum_ND [gen][det][quant][HC] = SpecToHist(NDFile, "VisE_sum_"+Comb,
138  POT[HC], ";Summed Vis E (GeV)");
139  sp.VisE_ND [gen][det][quant][HC] = SpecToHist(NDFile, "VisE_"+Comb,
140  POT[HC], ";Vis E (GeV)");
141 
142  // normalize the plots with variable binning by bin width
143  NormByBinWidth(sp.NuE_nom_ND[gen][det][quant][HC]);
144  NormByBinWidth(sp.NuE_scale_p_ND[gen][det][quant][HC]);
145  NormByBinWidth(sp.NuE_scale_m_ND[gen][det][quant][HC]);
146  }
147 
148  // analysis cuts, FD
149  for(unsigned int det=0; det<nFD_CutNames; ++det) {
150 
151  std::string Comb = GENIEStr[gen]+"_"+FD_CutNames[det]+"_"+QuantNames[quant];
152 
153  sp.NuE_nom_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NuE_nom_"+Comb,
154  POT[HC], ";Reco #nu E (GeV);Events/0.1GeV");
155  sp.HadE_nom_FD [gen][det][quant][HC] = SpecToHist(FDFile, "HadE_nom_"+Comb,
156  POT[HC], ";Reco Had E (GeV);Events");
157  sp.NuE_scale_p_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NuE_scale_p_"+Comb,
158  POT[HC], ";Reco #nu E (GeV);Events/0.1GeV");
159  sp.HadE_scale_p_FD [gen][det][quant][HC] = SpecToHist(FDFile, "HadE_scale_p_"+Comb,
160  POT[HC], ";Reco Had E (GeV);Events");
161  sp.NuE_scale_m_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NuE_scale_m_"+Comb,
162  POT[HC], ";Reco #nu E (GeV);Events/0.1GeV");
163  sp.HadE_scale_m_FD [gen][det][quant][HC] = SpecToHist(FDFile, "HadE_scale_m_"+Comb,
164  POT[HC], ";Reco Had E (GeV);Events");
165  sp.NNeu_FD [gen][det][quant][HC] = SpecToHist(FDFile, "NNeu_"+Comb,
166  POT[HC], ";Multiplicity");
167  sp.VisE_sum_FD [gen][det][quant][HC] = SpecToHist(FDFile, "VisE_sum_"+Comb,
168  POT[HC], ";Summed Vis E (GeV)");
169  sp.VisE_FD [gen][det][quant][HC] = SpecToHist(FDFile, "VisE_"+Comb,
170  POT[HC], ";Vis E (GeV)");
171 
172  // normalize the plots with variable binning by bin width
173  NormByBinWidth(sp.NuE_nom_FD[gen][det][quant][HC]);
174  NormByBinWidth(sp.NuE_scale_p_FD[gen][det][quant][HC]);
175  NormByBinWidth(sp.NuE_scale_m_FD[gen][det][quant][HC]);
176  }
177  }
178  }
179  }
180 
181  // ----- comparison plots -----
182  //
183 
184  //final cut plots
185  for (unsigned int HC=0; HC<nHC; ++HC){
186  for (unsigned int gen=0; gen<nGENIE; ++gen) {
187  for (unsigned int quant=0; quant<nQuantNames; ++quant) {
188 
189  // --- Now I can make the combination string...
190  std::string Comb = GENIEStr[gen]+"_"+QuantNames[quant]+"_"+HCNames[HC];
191  MakeComparisonPlotsWithRatio({sp.NuE_nom_ND[gen][FCut_ND][quant][HC], sp.NuE_scale_p_ND[gen][FCut_ND][quant][HC]},
192  "NuE_scale_p_ND_"+Comb, {"no scale", "scale up"}, true);
193  MakeComparisonPlotsWithRatio({sp.NuE_nom_ND[gen][FCut_ND][quant][HC], sp.NuE_scale_m_ND[gen][FCut_ND][quant][HC]},
194  "NuE_scale_m_ND_"+Comb, {"no scale", "scale down"}, true);
195  MakeComparisonPlotsWithRatio({sp.HadE_nom_ND[gen][FCut_ND][quant][HC], sp.HadE_scale_p_ND[gen][FCut_ND][quant][HC]},
196  "HadE_scale_p_ND_"+Comb, {"no scale", "scale up"}, false);
197  MakeComparisonPlotsWithRatio({sp.HadE_nom_ND[gen][FCut_ND][quant][HC], sp.HadE_scale_m_ND[gen][FCut_ND][quant][HC]},
198  "HadE_scale_m_ND_"+Comb, {"no scale", "scale down"}, false);
199  Print1D(sp.NNeu_ND [gen][FCut_ND][quant][HC], "NNeu_ND_"+Comb);
200  Print1D(sp.VisE_sum_ND[gen][FCut_ND][quant][HC], "VisE_sum_ND_"+Comb);
201  Print1D(sp.VisE_ND [gen][FCut_ND][quant][HC], "VisE_ND_"+Comb);
202 
203  MakeComparisonPlotsWithRatio({sp.NuE_nom_FD[gen][FCut_FD][quant][HC], sp.NuE_scale_p_FD[gen][FCut_FD][quant][HC]},
204  "NuE_scale_p_FD_"+Comb, {"no scale", "scale up"}, true);
205  MakeComparisonPlotsWithRatio({sp.NuE_nom_FD[gen][FCut_FD][quant][HC], sp.NuE_scale_m_FD[gen][FCut_FD][quant][HC]},
206  "NuE_scale_m_FD_"+Comb, {"no scale", "scale down"}, true);
207  MakeComparisonPlotsWithRatio({sp.HadE_nom_FD[gen][FCut_FD][quant][HC], sp.HadE_scale_p_FD[gen][FCut_FD][quant][HC]},
208  "HadE_scale_p_FD_"+Comb, {"no scale", "scale up"}, false);
209  MakeComparisonPlotsWithRatio({sp.HadE_nom_FD[gen][FCut_FD][quant][HC], sp.HadE_scale_m_FD[gen][FCut_FD][quant][HC]},
210  "HadE_scale_m_FD_"+Comb, {"no scale", "scale down"}, false);
211 
212  Print1D(sp.NNeu_FD [gen][FCut_FD][quant][HC], "NNeu_FD_"+Comb);
213  Print1D(sp.VisE_sum_FD[gen][FCut_FD][quant][HC], "VisE_sum_FD_"+Comb);
214  Print1D(sp.VisE_FD [gen][FCut_FD][quant][HC], "VisE_FD_"+Comb);
215 
216  } // quants
217  } // int. types
218  } // hc
219 
220  return;
221 }
222 
223 //......................................................
224 TH1D* SpecToHist( TFile* MyF, std::string SpecName, double POT, std::string AxisTitle ) {
225  std::unique_ptr<Spectrum> TempSpec = Spectrum::LoadFrom(MyF, SpecName.c_str());
226  TH1D* MyHist = TempSpec->ToTH1(POT);
227  ana::CenterTitles(MyHist);
228  MyHist->SetMinimum(0);
229  MyHist->SetMaximum(MyHist->GetMaximum()*1.1);
230  MyHist->SetTitle(AxisTitle.c_str());
231 
232  //area norm
233 // MyHist->Scale(1./MyHist->Integral());
234 
235  // comment this out to keep the 0bin
236 // MyHist->GetXaxis()->SetRange(2,MyHist->GetNbinsX());
237 
238  return MyHist;
239 }
240 
241 //......................................................
242 void NormByBinWidth(TH1D* hist){
243  hist->Scale(0.1, "width");
244  hist->SetMaximum(hist->GetMaximum()*1.1);
245 }
246 
247 //......................................................
249  TCanvas* hCan = new TCanvas(plotName.c_str(), plotName.c_str());
250  hCan->SetLeftMargin(0.1);
251  hist->DrawCopy("hist");
252  hCan -> SaveAs(TString("Plots/")+TString(hCan->GetName())+TString(".pdf"));
253  hCan -> SaveAs(TString("Plots/")+TString(hCan->GetName())+TString(".png"));
254 
255 }
256 //......................................................
257 void MakeComparisonPlots( std::vector<TH1D*> TheHists, std::string plotName, std::vector<std::string> Labels, bool var ) {
258  TCanvas* hCan = new TCanvas(plotName.c_str(), plotName.c_str());
259  hCan->SetLeftMargin(0.1);
260  int NBins = TheHists[0]->GetNbinsX();
261  double max = 0.0;
262 
263  // --- Legend ---
264  double HistLeg_X1 = 0.7, HistLeg_Y1 = 0.65, HistLeg_X2 = 0.9, HistLeg_Y2 = 0.88;
265  TLegend* Leg = new TLegend(HistLeg_X1, HistLeg_Y1, HistLeg_X2, HistLeg_Y2 );
266  Leg->SetFillStyle(0);
267  Leg->SetBorderSize(0);
268  double size = 0.9*TheHists[0]->GetXaxis()->GetTitleSize();
269  Leg->SetTextSize( size );
270 
271  for (size_t hist=0; hist<TheHists.size(); ++hist) {
272  TheHists[hist]->SetLineColor( 1+hist );
273  if (hist==0) {
274  TheHists[hist]->GetYaxis()->SetTitleOffset(0.7);
275  TheHists[hist]->DrawCopy("hist");
276  if(TString(plotName).Contains("NuE")){
277  std::cout<<"PlotName: "<<plotName[hist]<<" Label: "<<Labels[hist]<<std::endl;
278  FitGauss(TheHists[hist]);
279  }
280  }
281  else{
282  TheHists[hist] -> DrawCopy("hist same");
283  if(TString(plotName).Contains("NuE")){
284  std::cout<<"PlotName: "<<plotName[hist]<<" Label: "<<Labels[hist]<<std::endl;
285  FitGauss(TheHists[hist]);
286  }
287  }
288  Leg->AddEntry(TheHists[hist], Labels[hist].c_str());
289  if(TheHists[hist]->GetMaximum() > max) max = TheHists[hist]->GetMaximum();
290  }
291  // set the new maximum
292  TheHists[0]->SetMaximum(max);
293  Leg->Draw();
294 
295  hCan -> SaveAs(TString("Plots/")+TString(hCan->GetName())+TString(".pdf"));
296  hCan -> SaveAs(TString("Plots/")+TString(hCan->GetName())+TString(".png"));
297 }
298 //......................................................
299 void MakeComparisonPlotsWithRatio( std::vector<TH1D*> TheHists, std::string plotName, std::vector<std::string> Labels, bool var ) {
300  // making two canvases -- one with just the comparison, and one with a ratio
301 
302  TCanvas* hCan = new TCanvas(plotName.c_str(), plotName.c_str());
303  hCan->SetLeftMargin(0.1);
304  int NBins = TheHists[0]->GetNbinsX();
305  double max = 0.0;
306 
307  // --- Legend ---
308  double HistLeg_X1 = 0.7, HistLeg_Y1 = 0.65, HistLeg_X2 = 0.9, HistLeg_Y2 = 0.88;
309  TLegend* Leg = new TLegend(HistLeg_X1, HistLeg_Y1, HistLeg_X2, HistLeg_Y2 );
310  Leg->SetFillStyle(0);
311  Leg->SetBorderSize(0);
312  double size = 0.9*TheHists[0]->GetXaxis()->GetTitleSize();
313  Leg->SetTextSize( size );
314 
315  for (size_t hist=0; hist<TheHists.size(); ++hist) {
316  TheHists[hist]->SetLineColor( 1+hist );
317  if (hist==0) {
318  TheHists[hist]->GetYaxis()->SetTitleOffset(0.7);
319  TheHists[hist]->DrawCopy("hist");
320  if(TString(plotName).Contains("NuE")){
321  std::cout<<"PlotName: "<<plotName<<" Label: "<<Labels[hist]<<std::endl;
322  FitGauss(TheHists[hist]);
323  }
324  }
325  else{
326  TheHists[hist] -> DrawCopy("hist same");
327  if(TString(plotName).Contains("NuE")){
328  std::cout<<"PlotName: "<<plotName<<" Label: "<<Labels[hist]<<std::endl;
329  FitGauss(TheHists[hist]);
330  }
331  }
332  Leg->AddEntry(TheHists[hist], Labels[hist].c_str());
333  if(TheHists[hist]->GetMaximum() > max) max = TheHists[hist]->GetMaximum();
334  }
335  // set the new maximum
336  TheHists[0]->SetMaximum(max);
337  Leg->Draw();
338 
339  // now we're going to make the canvas with the plot + ratio plot on it
340  std::string rCname = "r"+plotName;
341 
342  // --- Get ready to make the canvas with ratios.
343  TCanvas* rC = new TCanvas(rCname.c_str(), rCname.c_str());
344 
345  // Sort out the canvas.
346  rC -> SetBottomMargin(0.);
347  double Spl = 0.3;
348  TPad* P1 = new TPad( "Temp_1", "", 0.0, Spl, 1.0, 1.0, 0 );
349  TPad* P2 = new TPad( "Temp_2", "", 0.0, 0.0, 1.0, Spl, 0 );
350  P2 -> SetRightMargin (.03);
351  P2 -> SetTopMargin (.00);
352  P2 -> SetBottomMargin(.22);
353  P2 -> SetLeftMargin (.1);
354  P2 -> Draw();
355  P1 -> SetRightMargin (.03);
356  P1 -> SetLeftMargin (.1);
357  P1 -> SetTopMargin (.10);
358  P1 -> SetBottomMargin(.00);
359  P1 -> Draw();
360  // Set some label sizes.
361  double Lb1 = 0.07;
362  double Lb2 = 0.10;
363  // --- First, draw the ratios so cd onto Pad2
364  P2 -> cd();
365 
366  TH1D* Rat = MakeRatio(TheHists[0], TheHists[1], kBlack, rCname, var); // Make the ratio.
367  // Set axis ranges etc.
368  Rat->GetYaxis()->SetTitleSize( Lb2 );
369  Rat->GetYaxis()->SetTitleOffset(0.4);
370  Rat->GetYaxis()->SetLabelSize( Lb2 );
371  Rat->GetXaxis()->SetTitleSize( Lb2 );
372  Rat->GetXaxis()->SetLabelSize( Lb2 );
373  Rat->Draw();
374  // Add a reference line
375  double x[2] = {TheHists[0]->GetXaxis()->GetXmin(), TheHists[0]->GetXaxis()->GetXmax()};
376  double y[2] = {1, 1};
377  TGraph *line = new TGraph(2, x, y);
378  line->SetLineColor(kGray);
379  line->SetLineWidth(3);
380  line->SetLineStyle(2);
381  line->Draw("l same");
382 
383  Rat->Draw("axis same");
384  P1->cd();
385 
386  // hist goes on P1
387  double Leg_X1 = 0.75, Leg_Y1 = 0.65, Leg_X2 = 0.95, Leg_Y2 = 0.88;
388  TLegend* Leg2 = new TLegend(Leg_X1, Leg_Y1, Leg_X2, Leg_Y2 );
389  Leg2->SetFillStyle (0);
390  Leg2->SetBorderSize(0);
391  Leg2->SetTextSize( 1.2*size );
392 
393  max = 0.0;
394  for (size_t hist=0; hist<TheHists.size(); ++hist) {
395  TheHists[hist]->SetLineColor( 1+hist );
396  if (hist==0) {
397  TheHists[hist]->GetYaxis()->SetTitleSize( Lb1 );
398  TheHists[hist]->GetYaxis()->SetLabelSize( Lb1 );
399  TheHists[hist]->GetYaxis()->SetTitleOffset( 0.7 );
400  // Remove the x axis labels
401  TheHists[hist]->GetXaxis()->SetLabelSize (0 );
402  TheHists[hist]->GetXaxis()->SetLabelOffset(99);
403  TheHists[hist]->DrawCopy("hist");
404  }
405  else{
406  TheHists[hist]->DrawCopy("hist same");
407  }
408  Leg2->AddEntry(TheHists[hist], Labels[hist].c_str());
409  if(TheHists[hist]->GetMaximum() > max) max = TheHists[hist]->GetMaximum();
410  }
411  // set the new maximum
412  TheHists[0]->SetMaximum(max);
413  Leg2->Draw();
414 
415  hCan -> SaveAs( TString("Plots/")+TString(hCan->GetName())+TString(".pdf") );
416  hCan -> SaveAs( TString("Plots/")+TString(hCan->GetName())+TString(".png") );
417  rC -> SaveAs( TString("Plots/")+TString(rC->GetName())+TString(".pdf") );
418  rC -> SaveAs( TString("Plots/")+TString(rC->GetName())+TString(".png") );
419 
420  return;
421 }
422 
423 TH1D* MakeRatio( TH1D* num, TH1D* denom, int Col, std::string FType, bool var ) {
424  // I know this function is awful but it's 1 am so this is what we get
425  if(var){
426  Double_t bins[] = { 0, 0.75, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.25, 2.5, 2.75, 3, 3.5, 4, 5 };
427  Int_t binnum = sizeof(bins)/sizeof(Double_t) - 1;
428  TH1D* ratHist = new TH1D( TString(FType)+TString("_ratio")+TString(std::to_string(Col))+TString(num->GetName()),
429  TString(";")+TString(num->GetXaxis()->GetTitle())+TString(";Ratio"),
430  binnum,
431  bins
432  );
433  // Don't want any exponents?
434  ratHist->GetYaxis()->SetNoExponent();
435  ratHist->GetXaxis()->SetNoExponent();
436  // Divide the two input histograms...
437  ratHist->Divide( num , denom );
438  // Set Line and Marker colors / styles.
439  ratHist->SetLineColor ( Col );
440  ratHist->SetMarkerColor( Col );
441  ratHist->SetMarkerStyle( 20 );
442  ratHist->SetMarkerSize ( 0 );
443  // Figure out what the range should be...
444  ratHist->GetYaxis()->SetRangeUser( 0, 2.0 );
445  // Return my shiny new ratio plot!
446  return ratHist;
447 
448  }
449 
450  else{
451  TH1D* ratHist = new TH1D( TString(FType)+TString("_ratio")+TString(std::to_string(Col))+TString(num->GetName()),
452  TString(";")+TString(num->GetXaxis()->GetTitle())+TString(";Ratio"),
453  num->GetNbinsX(),
454  num->GetXaxis()->GetXmin(),
455  num->GetXaxis()->GetXmax()
456  );
457  // Don't want any exponents?
458  ratHist->GetYaxis()->SetNoExponent();
459  ratHist->GetXaxis()->SetNoExponent();
460  // Divide the two input histograms...
461  ratHist->Divide( num , denom );
462  // Set Line and Marker colors / styles.
463  ratHist->SetLineColor ( Col );
464  ratHist->SetMarkerColor( Col );
465  ratHist->SetMarkerStyle( 20 );
466  ratHist->SetMarkerSize ( 0 );
467  // Figure out what the range should be...
468  ratHist->GetYaxis()->SetRangeUser( 0, 2.0 );
469  // Return my shiny new ratio plot!
470  return ratHist;
471 
472  }
473 }
474 
475 void FitGauss(TH1D* h)
476 {
477  double mean, bin1, bin2, fwhm;
478  mean = h->GetMean();
479  bin1 = h->FindFirstBinAbove(h->GetMaximum()/2);
480  bin2 = h->FindLastBinAbove(h->GetMaximum()/2);
481  fwhm = h->GetBinCenter(bin2 ) - h->GetBinCenter(bin1 );
482 
483  int ilo = 1;
484  int ihi = h->GetNbinsX();
485 
486  int imx = h->GetMaximumBin();
487  double mx = h->GetBinContent(imx);
488 
489  int ileft;
490  for (ileft=imx; ileft>=ilo; --ileft){
491  if (h->GetBinContent(ileft)<0.5*mx) {
492  ++ileft;
493  break;
494  }
495  }
496 
497  int iright;
498  for (iright=imx; iright<=ihi; ++iright) {
499  if (h->GetBinContent(iright)<0.5*mx) {
500  break;
501  }
502  }
503 
504  TF1 *myFit = new TF1("myFit","gaus", h->GetBinLowEdge(ileft), h->GetBinLowEdge(iright));
505  h->Fit(myFit,"R q 0");
506  myFit->SetLineColor(h->GetLineColor());
507  std::cout<<"Mean: " << myFit->GetParameter(1) << "+/-" << myFit->GetParError(1) <<std::endl;
508  return;
509 }
TH1D * VisE_sum_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
tree Draw("slc.nhit")
TH1D * MakeRatio(TH1D *num, TH1D *denom, int Col, std::string FType, bool var)
std::vector< std::string > FD_CutNames
void neutKEsyst_plot()
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
TH1D * NuE_nom_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
TH1D * HadE_nom_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
TH1D * NNeu_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
void Print1D(TH1D *hist, std::string plotName)
TH1D * NuE_scale_m_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
std::vector< std::string > QuantNames
unsigned int FCut_FD
TH1D * VisE_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
TH1D * NuE_nom_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
TH1D * NuE_scale_p_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
TH1D * HadE_scale_m_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
TH1D * HadE_scale_p_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
void NormByBinWidth(TH1D *hist)
const unsigned int nDet
void FitGauss(TH1D *h)
std::vector< std::string > DetNames
void MakeComparisonPlotsWithRatio(std::vector< TH1D * > TheHists, std::string plotName, std::vector< std::string > Labels, bool var)
const double kAna2018RHCPOT
Definition: Exposures.h:208
const unsigned int nFD_CutNames
TH1D * VisE_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
TH1D * HadE_scale_m_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
TH1D * NuE_scale_m_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
const unsigned int nND_CutNames
const std::vector< std::string > GENIEStr
const unsigned int nQuantNames
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
const Binning bins
Definition: NumuCC_CPiBin.h:8
std::vector< std::string > HCNames
const unsigned int nGENIE
std::vector< std::string > ND_CutNames
TH1D * VisE_sum_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
const unsigned int nHC
TH1D * NNeu_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
int num
Definition: f2_nu.C:119
TH1D * SpecToHist(TFile *MyF, std::string SpecName, double POT, std::string AxisTitle)
TH1D * HadE_nom_FD[nGENIE][nFD_CutNames][nQuantNames][nHC]
unsigned int FCut_ND
void MakeComparisonPlots(std::vector< TH1D * > TheHists, std::string plotName, std::vector< std::string > Labels, bool var)
TH1D * HadE_scale_p_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
const double kAna2018FHCPOT
Definition: Exposures.h:207
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
cosmicTree SaveAs("cosmicTree.root")
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
TH1D * NuE_scale_p_ND[nGENIE][nND_CutNames][nQuantNames][nHC]
c cd(1)
enum BeamMode string