makeResolutionPlots.C
Go to the documentation of this file.
1 // I found that the energy estimator used by the nuebar oscillation analysis
2 // was not sufficient for this analysis due to training on a sample of
3 // neutrinos with energy less than 4.5 GeV.
4 //
5 // This script marks the first piece of training a new estimator.
6 // Following docdb-26696, we need to normalize the true energy distribution
7 // to avoid biasing the energy estimation toward our peak beam energy.
8 //
9 // Plot the energy resolution of the various energy estimators
10 // and correct for the shift
11 
12 // dddoyle@colostate.edu
13 
14 #include "CAFAna/Core/Spectrum.h"
15 
16 #include "TH1D.h"
17 #include "TCanvas.h"
18 #include "TF1.h"
19 #include "TStyle.h"
20 #include "TGraph.h"
21 
25 #include "NDAna/nuebarcc_inc/EnergyEstimator/NuebarEnergyEstimatorHelper.h"
26 
27 using namespace ana;
28 
29 void PrintResolutionPlot(TH1D* h, std::string name, TFile* fout, std::string fit);
30 void PrintFitResults(TH1D* h, std::string name, std::string fit);
31 void PrintMeanRMSScanPlot(std::vector<ResolutionScan*> res_scan, std::string dir, TFile* fout);
32 
33 std::vector<std::string> fits = {"poly2", "poly2x",
34  "poly3", "poly3x",
35  "nue2018"};
36 
37 
39 {
40  std::string directory = "/nova/ana/users/ddoyle/NuebarEnergyEstimatorFiles/";
41  //TFile* fRes = TFile::Open((directory + "energy_resolution_spectra_training_sample.root").c_str());
42  TFile* fRes = TFile::Open((directory + "energy_resolution_spectra.root").c_str());
43  TFile* fout = new TFile((directory+"Plots/energy_resolution_plots.root").c_str(), "recreate");
44  fRes->cd();
45  std::vector<bool> isShift = {true, false};
46 
47  /////////////////////////////////////////////////////////////////////////
48  // Inclusive Resolution
49  /////////////////////////////////////////////////////////////////////////
50  std::vector<Spectrum*> sRes_shift;
51  std::vector<Spectrum*> sRes_noshift;
52  std::vector<Spectrum*> sRes_shift_shape;
53  std::vector<Spectrum*> sRes_noshift_shape;
54  std::vector<TH1D*> hRes_shift;
55  std::vector<TH1D*> hRes_noshift;
56  std::vector<TH1D*> hRes_shift_shape;
57  std::vector<TH1D*> hRes_noshift_shape;
58  for(std::string fit : fits) {
59  sRes_noshift.push_back(Spectrum::LoadFrom(fRes, ("res_" + fit).c_str())).release();
60  sRes_shift.push_back(Spectrum::LoadFrom(fRes, ("res_" + fit + "_shifted").c_str())).release();
61  sRes_noshift_shape.push_back(Spectrum::LoadFrom(fRes, ("res_nocut_" + fit + "_shape").c_str())).release();
62  sRes_shift_shape.push_back(Spectrum::LoadFrom(fRes, ("res_nocut_" + fit + "_shifted_shape").c_str())).release();
63  }
64  double POT = sRes_shift[0]->POT();
65  for(int i = 0; i < (int) fits.size(); i++) {
66  std::cout << "Fitting " << fits[i] << "flat noshift" << std::endl;
67  hRes_noshift.push_back(sRes_noshift[i]->ToTH1(POT));
68  hRes_noshift[i]->SetTitle(fits[i].c_str());
69  hRes_noshift[i]->Fit("gaus");
70  hRes_noshift[i]->Fit("gaus", "", "", -0.2, 0.2);
71  std::cout << "Fitting " << fits[i] << "flat shift" << std::endl;
72  hRes_shift.push_back(sRes_shift[i]->ToTH1(POT));
73  hRes_shift[i]->SetTitle(fits[i].c_str());
74  hRes_shift[i]->Fit("gaus");
75  hRes_shift[i]->Fit("gaus", "", "", -0.2, 0.2);
76 
77  // shifted
78  std::cout << "Fitting " << fits[i] << "shape noshift" << std::endl;
79  hRes_noshift_shape.push_back(sRes_noshift_shape[i]->ToTH1(POT));
80  hRes_noshift_shape[i]->SetTitle(fits[i].c_str());
81  hRes_noshift_shape[i]->Fit("gaus");
82  hRes_noshift_shape[i]->Fit("gaus", "", "", -0.2, 0.2);
83  std::cout << "Fitting " << fits[i] << "shape shift" << std::endl;
84  hRes_shift_shape.push_back(sRes_shift_shape[i]->ToTH1(POT));
85  hRes_shift_shape[i]->SetTitle(fits[i].c_str());
86  hRes_shift_shape[i]->Fit("gaus");
87  hRes_shift_shape[i]->Fit("gaus", "", "", -0.2, 0.2);
88  }
89 
91  for(int i = 0; i < (int) fits.size(); i ++) {
92  PrintResolutionPlot(hRes_noshift[i],
93  TString::Format("%sPlots/resolution_%s_noshift.pdf",
94  directory.c_str(),
95  fits[i].c_str()).Data(),
96  fout,
97  "");
98  PrintResolutionPlot(hRes_shift[i],
99  TString::Format("%sPlots/resolution_%s_shift.pdf",
100  directory.c_str(),
101  fits[i].c_str()).Data(),
102  fout,
103  "");
104  PrintResolutionPlot(hRes_noshift_shape[i],
105  TString::Format("%sPlots/resolution_%s_noshift_shape.pdf",
106  directory.c_str(),
107  fits[i].c_str()).Data(),
108  fout,
109  "");
110  PrintResolutionPlot(hRes_shift_shape[i],
111  TString::Format("%sPlots/resolution_%s_shift_shape.pdf",
112  directory.c_str(),
113  fits[i].c_str()).Data(),
114  fout,
115  fits[i].c_str());
116  }
117 
118  /////////////////////////////////////////////////////////////////////////
119  // Resolution by interaction
120  /////////////////////////////////////////////////////////////////////////
121  std::vector<InteractionSpectra*> int_specs;
122  for(std::string fit : fits) {
123  for(bool is_shifted : {true, false}) {
124  for(bool is_flat : {true, false}) {
125  int_specs.push_back(new InteractionSpectra(fRes, is_shifted, is_flat, fit));
126  }
127  }
128  }
129  for(InteractionSpectra* spec : int_specs ) {
130  spec->MakePlot(directory + "Plots/", fout);
131  }
132 
133  /////////////////////////////////////////////////////////////////////////
134  // Resolution scan
135  /////////////////////////////////////////////////////////////////////////
136  std::vector<ResolutionScan*> res_scan_specs;
137  for(std::string fit : fits) {
138  for(bool is_shifted : {true, false}) {
139  for(bool is_flat : {true, false}) {
140  // if(fit=="nue2018" && is_shifted) continue;
141  for(std::string xvar : {"trueE", "EME", "HADE"}) {
142  res_scan_specs.push_back(new ResolutionScan(fRes, is_shifted,
143  is_flat,
144  fit,
145  xvar));
146  }
147  }
148  }
149  }
150  for(ResolutionScan* scan : res_scan_specs ) {
151  scan->MakePlot(directory + "Plots/", fout);
152  }
153 
154  /////////////////////////////////////////////////////////////////////////
155  // Mean and RMS of Resolution scans by fit
156  /////////////////////////////////////////////////////////////////////////
157  std::vector<ResolutionScan*> res_mean_scan_trueE_shift;
158  std::vector<ResolutionScan*> res_mean_scan_trueE_noshift;
159  std::vector<ResolutionScan*> res_mean_scan_EME_shift;
160  std::vector<ResolutionScan*> res_mean_scan_EME_noshift;
161  std::vector<ResolutionScan*> res_mean_scan_HADE_shift;
162  std::vector<ResolutionScan*> res_mean_scan_HADE_noshift;
163  for(std::string fit : fits) {
164  res_mean_scan_trueE_shift .push_back(new ResolutionScan(fRes, true, true, fit, "trueE"));
165  res_mean_scan_EME_shift .push_back(new ResolutionScan(fRes, true, true, fit, "EME"));
166  res_mean_scan_HADE_shift .push_back(new ResolutionScan(fRes, true, true, fit, "HADE"));
167  res_mean_scan_trueE_noshift.push_back(new ResolutionScan(fRes, false, true, fit, "trueE"));
168  res_mean_scan_EME_noshift .push_back(new ResolutionScan(fRes, false, true, fit, "EME"));
169  res_mean_scan_HADE_noshift .push_back(new ResolutionScan(fRes, false, true, fit, "HADE"));
170  }
171  PrintMeanRMSScanPlot(res_mean_scan_trueE_shift, directory + "Plots/", fout);
172  PrintMeanRMSScanPlot(res_mean_scan_trueE_noshift, directory + "Plots/", fout);
173  PrintMeanRMSScanPlot(res_mean_scan_EME_shift, directory + "Plots/", fout);
174  PrintMeanRMSScanPlot(res_mean_scan_EME_noshift, directory + "Plots/", fout);
175  PrintMeanRMSScanPlot(res_mean_scan_HADE_shift, directory + "Plots/", fout);
176  PrintMeanRMSScanPlot(res_mean_scan_HADE_noshift, directory + "Plots/", fout);
177 
178 }
179 
180 void PrintMeanRMSScanPlot(std::vector<ResolutionScan*> res_scan, std::string dir, TFile* fout)
181 {
182  TDirectory* tmp = gDirectory;
183 
184  std::vector<Int_t> colors = {kViolet+1, kCyan, kMagenta, kGreen, kOrange+10};
185  std::string xvar = res_scan[0]->GetXVar();
186  bool is_shifted = res_scan[0]->IsShifted();
187  std::string plot_name_mean = "mean_res_by_fit_" + xvar;
188  std::string plot_name_rms = "rms_by_fit_" + xvar;
189  plot_name_mean += is_shifted ? "_shift.pdf" : "_noshift.pdf";
190  plot_name_rms += is_shifted ? "_shift.pdf" : "_noshift.pdf";
191  std::vector<TH1D*> means;
192  std::vector<TH1D*> rmss;
193  for(ResolutionScan* res : res_scan) {
194  means.push_back(res->MeanScan());
195  rmss.push_back(res->RMSScan());
196  }
197 
198 
199  TCanvas* c = new TCanvas();
200  c->cd();
201  // TLegend* leg = new TLegend(0.117, 0.58, 0.47, 0.73);
202  TLegend* leg_mean = new TLegend();
203  TLegend* leg_rms = new TLegend();
204  for(int i = 0; i < (int) res_scan.size(); i++) {
205  means[i]->SetLineColor(colors[i]);
206  means[i]->SetLineStyle(kSolid);
207  if(i==0) {
208  means[i]->GetXaxis()->SetTitle(xvar.c_str());
209  means[i]->GetYaxis()->SetTitle("Mean of resolution");
210  means[i]->Draw("hist el");
211  }
212  means[i]->Draw("same el");
213  leg_mean->AddEntry(means[i], fits[i].c_str(), "l");
214  leg_rms->AddEntry(rmss[i], fits[i].c_str(), "l");
215  }
216  leg_mean->Draw("same");
217  means[0]->GetYaxis()->SetRangeUser(-0.5, 0.5);
218  TGraph* g = new TGraph();
219  g->SetPoint(0, -1, 0);
220  g->SetPoint(1, 11, 0);
221  g->SetLineColor(kRed);
222  g->SetLineStyle(7);
223  g->Draw("same");
224  fout->cd();
225  c->SaveAs((dir + plot_name_mean).c_str());
226  c->Write((dir + plot_name_mean).c_str());
227 
228  c->cd();
229  for(int i = 0; i < (int) res_scan.size(); i++) {
230  rmss[i]->SetLineColor(colors[i]);
231  if(i==0) {
232  rmss[i]->GetXaxis()->SetTitle(xvar.c_str());
233  rmss[i]->GetYaxis()->SetTitle("RMS of resolution");
234  rmss[i]->GetYaxis()->SetRangeUser(0, 0.5);
235  rmss[i]->Draw("hist e");
236  }
237  rmss[i]->Draw("same hist e");
238  }
239  leg_rms->Draw("same");
240  fout->cd();
241  c->SaveAs((dir + plot_name_rms).c_str());
242  c->Write((dir + plot_name_rms).c_str());
243 
244  delete c;
245  tmp->cd();
246 }
247 
249 {
250  std::ofstream fit_results;
251  std::string fit_results_name = fit + "_fit_results_shifted.txt";
252  if(fit == "") fit_results_name = "trash.txt";
253  fit_results.open("FitResults/" + fit_results_name, std::ios::app);
254  TF1* func = h->GetFunction("gaus");
255  double avg = h->GetMean();
256  double RMS = h->GetRMS();
257  double fit_mean = func->GetParameter(1);
258  double fit_std = func->GetParameter(2);
259  std::cout << "-----------------------------------------\n";
260  std::cout << "Gaussian fit Results for " << name << std::endl;
261  std::cout << "Histo MEAN: " << avg << std::endl;
262  std::cout << "Histo RMS: " << RMS << std::endl;
263  std::cout << "Func MEAN: " << fit_mean << std::endl;
264  std::cout << "Func STD: " << fit_std << std::endl;
265  std::cout << "-----------------------------------------\n";
266  fit_results << "mean=" << avg << std::endl;
267 }
268 
270 {
271  TDirectory* tmp = gDirectory;
272  TGraph* g = new TGraph();
273  g->SetPoint(0, 0, -10);
274  g->SetPoint(1, 0, 1e6);
275  g->SetLineColor(kRed);
276  g->SetLineStyle(2);
277 
278  TF1* func = h->GetFunction("gaus");
279  double avg = h->GetMean();
280  double RMS = h->GetRMS();
281  double fit_mean = func->GetParameter(1);
282  double fit_std = func->GetParameter(2);
283  func->SetLineColor(kRed);
284  h->SetStats(kTRUE);
285  gStyle->SetOptFit(0101);
286  gStyle->SetOptStat(1001100);
287  TCanvas* c = new TCanvas();
288  c->cd();
289 
290  h->Draw("hist");
291  func->Draw("same");
292  g->Draw("same");
293  fout->cd();
294  c->SaveAs(name.c_str());
295  c->Write(name.c_str());
296  delete c;
297 
298  PrintFitResults(h, name, fit);
299 
300  tmp->cd();
301 }
const XML_Char * name
Definition: expat.h:151
enum BeamMode kOrange
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void PrintFitResults(TH1D *h, std::string name, std::string fit)
string directory
projection from multiple chains
Float_t tmp
Definition: plot.C:36
int colors[6]
Definition: tools.h:1
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
void makeResolutionPlots()
double func(double x, double y)
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
TDirectory * dir
Definition: macro.C:5
enum BeamMode kViolet
app
Definition: demo.py:189
std::vector< std::string > fits
void PrintResolutionPlot(TH1D *h, std::string name, TFile *fout, std::string fit)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
enum BeamMode kGreen
void PrintMeanRMSScanPlot(std::vector< ResolutionScan * > res_scan, std::string dir, TFile *fout)
enum BeamMode string