22 int NY = hist->GetNbinsY();
23 double y_0 = hist->GetYaxis()->GetBinLowEdge(1);
24 double y_1 = hist->GetYaxis()->GetBinLowEdge(NY+1);
35 for(
auto j = 1;
j <= rebin->GetNbinsY();
j++) {
36 double lower_bound = 1;
37 for(
auto irebin = 1; irebin <= rebin->GetNbinsX(); irebin++) {
38 double rebin_edge = rebin->GetXaxis()->GetBinLowEdge(irebin+1);
39 for(
auto i = 1;
i <= hist->GetNbinsX();
i++) {
40 int gid = rebin->FindBin(hist->GetXaxis()->GetBinCenter(
i),
41 hist->GetYaxis()->GetBinCenter(
j));
43 rebin->GetBinXYZ(gid, x, y, z);
44 rebin->SetBinContent(x, y, z, rebin->GetBinContent(x, y, z) + hist->GetBinContent(
i,
j));
55 double maxval =
std::max(reco->GetMaximum(), truth->GetMaximum());
57 TCanvas *
c =
new TCanvas();
58 TLegend *
leg =
new TLegend();
59 reco->SetLineColor(kBlack);
60 truth->SetLineColor(
kRed);
61 leg->AddEntry(reco,
"Reco",
"l");
62 leg->AddEntry(truth,
"Truth",
"l");
63 reco->GetXaxis()->SetTitle(
"Double Differential Analysis Bins");
64 reco->GetYaxis()->SetRangeUser(0, maxval*1.2);
66 truth->Draw(
"hist same");
68 c->SaveAs(name.c_str());
74 double maxval =
std::max(reco->GetMaximum(), truth->GetMaximum());
75 reco->GetXaxis()->SetRangeUser(xmin, xmax);
76 TCanvas *
c =
new TCanvas();
77 TLegend *
leg =
new TLegend();
78 reco->SetLineColor(kBlack);
79 truth->SetLineColor(
kRed);
80 leg->AddEntry(reco,
"Reco",
"l");
81 leg->AddEntry(truth,
"Truth",
"l");
82 reco->GetYaxis()->SetRangeUser(0, maxval*1.2);
84 truth->Draw(
"hist same");
86 c->SaveAs(name.c_str());
93 hist->GetYaxis()->CenterTitle();
95 int NX = hist->GetNbinsX();
97 double maxx = hist->GetXaxis()->GetBinLowEdge(NX+1);
98 double range = 0.5 * maxx;
101 auto rms =
new TH1D(
UniqueName().c_str(),
"", NX, hist->GetXaxis()->GetXbins()->GetArray());
102 rms->GetXaxis()->SetTitle(hist->GetXaxis()->GetTitle());
103 rms->GetXaxis()->CenterTitle();
104 for(
int i = 1;
i <= hist->GetNbinsX(); ++
i){
105 auto projection = (TH1F*) hist->ProjectionY(
"",
i,
i);
106 int lowbin = projection->FindBin(-range);
107 int highbin = projection->FindBin(range);
108 if(projection->Integral(lowbin, highbin)) {
109 projection->Fit(
"gaus",
"+Q",
"", -range, range);
110 double irms = projection->GetFunction(
"gaus")->GetParameter(2);
111 rms->SetBinContent(
i,irms);
112 rms->SetBinError(
i, 0);
116 auto bin_boxes =
new TH1D(
UniqueName().c_str(),
"", edges.size() -1, &edges[0]);
117 bin_boxes->SetLineColor(
kRed);
118 rms->SetLineColor(kBlack);
119 for(
auto i = 1;
i <= NX;
i++) {
120 bin_boxes->SetBinContent(
i, bin_boxes->GetBinWidth(
i));
123 double ymax =
rms->GetMaximum() * 1.2;
124 rms->GetYaxis()->SetRangeUser(0, ymax);
125 rms->GetXaxis()->SetRangeUser(xmin, xmax);
126 rms->GetYaxis()->SetTitle(
"Resolution");
127 rms->GetYaxis()->CenterTitle();
128 rms->GetYaxis()->SetTitleOffset(0.65);
131 TCanvas *
c =
new TCanvas();
133 bin_boxes->Draw(
"hist same");
136 TGaxis * right_axis =
new TGaxis(xmax,
141 rms->GetNdivisions(),
143 right_axis->SetLabelSize(
rms->GetYaxis()->GetLabelSize());
144 right_axis->SetLabelFont(
rms->GetYaxis()->GetLabelFont());
145 right_axis->SetLabelColor(
kRed);
146 right_axis->SetTitleSize(
rms->GetYaxis()->GetTitleSize());
147 right_axis->SetTitleFont(
rms->GetYaxis()->GetTitleFont());
148 right_axis->SetLineColor(
kRed);
149 right_axis->SetTitleColor(
kRed);
150 right_axis->SetTitleOffset(0.75);
151 right_axis->SetTitle(
"Bin Width");
152 right_axis->CenterTitle();
156 c->SaveAs(name.c_str());
163 gStyle->SetOptStat(1111);
164 gStyle->SetOptFit(0111);
165 gStyle->SetStatY(0.9);
166 gStyle->SetStatW(0.4);
167 gStyle->SetStatH(0.4);
168 TCanvas *
c =
new TCanvas();
169 if(
logy) c->SetLogy();
170 hist->GetXaxis()->SetTitle(
"Absolute Bias (Reco - True) [GeV]");
172 TF1 *
func = hist->GetFunction(
"gaus");
174 func->SetLineColor(
kRed);
177 c->SaveAs(name.c_str());
179 gStyle->SetOptStat(0);
180 gStyle->SetOptFit(0);
187 hist->GetYaxis()->CenterTitle();
188 int NX = hist->GetNbinsX();
189 auto markers =
new TH1D(
UniqueName().c_str(),
"", NX, hist->GetXaxis()->GetXbins()->GetArray());
190 if(hist->GetXaxis()->GetBinLowEdge(NX+1) !=
markers->GetXaxis()->GetBinLowEdge(NX+1)) {
191 double x0 = hist->GetXaxis()->GetBinLowEdge(1);
192 double x1 = hist->GetXaxis()->GetBinLowEdge(NX+1);
195 double range = 0.5 * hist->GetXaxis()->GetBinLowEdge(NX+1);
196 for(
int i = 1;
i <= hist->GetNbinsX(); ++
i){
197 auto projection = (TH1F*) hist->ProjectionY(
"",
i,
i);
198 int lowbin = projection->FindBin(-range);
199 int highbin = projection->FindBin(range);
200 if(projection->Integral(lowbin, highbin)) {
201 projection->Fit(
"gaus",
"+Q",
"", -range, range);
204 double mean = projection->GetFunction(
"gaus")->GetParameter(1);
205 double rms = projection->GetFunction(
"gaus")->GetParameter(2);
210 if(plot_projections) {
212 hist->GetXaxis()->GetBinLowEdge(
i),
213 hist->GetXaxis()->GetBinLowEdge(
i+1)));
214 projection->GetXaxis()->SetTitle(hist->GetYaxis()->GetTitle());
215 projection->GetYaxis()->SetTitle(
"Events");
216 projection->GetXaxis()->CenterTitle();
217 projection->GetYaxis()->CenterTitle();
218 projection->SetStats(
true);
219 size_t pos = name.find(
".pdf");
221 proj_name.replace(pos, 4,
"_projection_" +
std::to_string(
i) +
".pdf");
226 TCanvas *
c1 =
new TCanvas();
227 c1->SetRightMargin(0.15);
229 if(logz) c1->SetLogz();
232 hist->GetXaxis()->SetRangeUser(xmin, xmax);
233 if(tight) hist->GetYaxis()->SetRangeUser(-1, 1);
237 c1->SaveAs(name.c_str());
243 std::string plot_dump=
"/nova/ana/users/ddoyle/NuebarResolutionStudy/plots")
245 TFile *file_resolution =
new TFile(input_file_name.c_str(),
"READ");
247 if( !(file_resolution) ){
248 std::cout <<
"ResolutionSpectra.root does not exist" 249 <<
" run resolutionscript.C first" <<
std::endl;
253 size_t pos = input_file_name.find(
".root");
255 output_file_name.replace(pos, 5,
"_hists.root");
268 auto standard_axis_reco_angular_abs =*
Spectrum::LoadFrom(file_resolution,
"standard_axis_reco_angular_abs" );
269 auto standard_axis_truth_angular_abs =*
Spectrum::LoadFrom(file_resolution,
"standard_axis_truth_angular_abs");
270 auto standard_axis_reco_energy_abs =*
Spectrum::LoadFrom(file_resolution,
"standard_axis_reco_energy_abs" );
271 auto standard_axis_truth_energy_abs =*
Spectrum::LoadFrom(file_resolution,
"standard_axis_truth_energy_abs" );
273 auto signal_reco_energy_vs_cos =*
Spectrum::LoadFrom(file_resolution,
"signal_reco_energy_vs_cos" );
274 auto signal_truth_energy_vs_cos =*
Spectrum::LoadFrom(file_resolution,
"signal_truth_energy_vs_cos");
276 double MCPOT = reco_angular_frac.POT();
279 TH2 * hreco_angular_frac = reco_angular_frac .ToTH2(MCPOT);
280 TH2 * htruth_angular_frac = truth_angular_frac.ToTH2(MCPOT);
281 TH2 * hreco_energy_frac = reco_energy_frac .ToTH2(MCPOT);
282 TH2 * htruth_energy_frac = truth_energy_frac .ToTH2(MCPOT);
284 TH2 * hreco_angular_abs = reco_angular_abs .ToTH2(MCPOT);
285 TH2 * htruth_angular_abs = truth_angular_abs .ToTH2(MCPOT);
286 TH2 * hreco_energy_abs = reco_energy_abs .ToTH2(MCPOT);
287 TH2 * htruth_energy_abs = truth_energy_abs .ToTH2(MCPOT);
289 TH2 * hstandard_axis_reco_angular_abs = standard_axis_reco_angular_abs .ToTH2(MCPOT);
290 TH2 * hstandard_axis_truth_angular_abs = standard_axis_truth_angular_abs .ToTH2(MCPOT);
291 TH2 * hstandard_axis_reco_energy_abs = standard_axis_reco_energy_abs .ToTH2(MCPOT);
292 TH2 * hstandard_axis_truth_energy_abs = standard_axis_truth_energy_abs .ToTH2(MCPOT);
307 TH1 * hsignal_reco_energy_vs_cos = signal_reco_energy_vs_cos .ToTH1(DATAPOT);
308 TH1 * hsignal_truth_energy_vs_cos = signal_truth_energy_vs_cos.ToTH1(DATAPOT);
312 hreco_angular_frac ->GetYaxis()->SetTitle(
"Fractional Bias (Reco - True)/True");
313 htruth_angular_frac ->GetYaxis()->SetTitle(
"Fractional Bias (Reco - True)/True");
314 hreco_energy_frac ->GetYaxis()->SetTitle(
"Fractional Bias (Reco - True)/True");
315 htruth_energy_frac ->GetYaxis()->SetTitle(
"Fractional Bias (Reco - True)/True");
316 ResolutionPlot(hreco_angular_frac , 0, 1 , plot_dump +
"/angular_reco_res_frac.pdf",
true,
true);
317 ResolutionPlot(htruth_angular_frac , 0, 1 , plot_dump +
"/angular_true_res_frac.pdf",
true,
true);
318 ResolutionPlot(hreco_energy_frac , 0, 10, plot_dump +
"/energy_reco_res_frac.pdf" ,
true,
true);
319 ResolutionPlot(htruth_energy_frac , 0, 10, plot_dump +
"/energy_truth_res_frac.pdf",
true,
true);
321 hreco_angular_abs ->GetYaxis()->SetTitle(
"Absolute Bias (Reco - True) [GeV]");
322 htruth_angular_abs ->GetYaxis()->SetTitle(
"Absolute Bias (Reco - True) [GeV]");
323 hreco_energy_abs ->GetYaxis()->SetTitle(
"Absolute Bias (Reco - True) [GeV]");
324 htruth_energy_abs ->GetYaxis()->SetTitle(
"Absolute Bias (Reco - True) [GeV]");
325 ResolutionPlot(hreco_angular_abs , 0, 1 , plot_dump +
"/angular_reco_res_abs.pdf",
false,
true);
326 ResolutionPlot(htruth_angular_abs , 0, 1 , plot_dump +
"/angular_true_res_abs.pdf",
false,
true);
327 ResolutionPlot(hreco_energy_abs , 0, 10, plot_dump +
"/energy_reco_res_abs.pdf" ,
false,
true);
328 ResolutionPlot(htruth_energy_abs , 0, 10, plot_dump +
"/energy_truth_res_abs.pdf",
false,
true);
330 hstandard_axis_reco_angular_abs ->GetYaxis()->SetTitle(
"Absolute Bias (Reco - True) [GeV]");
331 hstandard_axis_truth_angular_abs ->GetYaxis()->SetTitle(
"Absolute Bias (Reco - True) [GeV]");
332 hstandard_axis_reco_energy_abs ->GetYaxis()->SetTitle(
"Absolute Bias (Reco - True) [GeV]");
333 hstandard_axis_truth_energy_abs ->GetYaxis()->SetTitle(
"Absolute Bias (Reco - True) [GeV]");
334 ResolutionPlot(hstandard_axis_reco_angular_abs , 0, 1 , plot_dump +
"/standard_axis_angular_reco_res_abs.pdf",
false,
true,
true);
335 ResolutionPlot(hstandard_axis_truth_angular_abs , 0, 1 , plot_dump +
"/standard_axis_angular_true_res_abs.pdf",
false,
true,
true);
336 ResolutionPlot(hstandard_axis_reco_energy_abs , 0, 10, plot_dump +
"/standard_axis_energy_reco_res_abs.pdf" ,
false,
true);
337 ResolutionPlot(hstandard_axis_truth_energy_abs , 0, 10, plot_dump +
"/standard_axis_energy_true_res_abs.pdf" ,
false,
true);
349 SignalCountPlot(hsignal_reco_energy_vs_cos, hsignal_truth_energy_vs_cos , plot_dump +
"/signal_energy_vs_cos.pdf" );
351 Spectrum * signal_reco_energy_vs_cos_2d =
new Spectrum(hsignal_reco_energy_vs_cos,
357 Spectrum * signal_truth_energy_vs_cos_2d =
new Spectrum(hsignal_truth_energy_vs_cos,
364 TH2D* hsignal_reco_energy_vs_cos_2d = (TH2D*) signal_reco_energy_vs_cos_2d ->
ToTH2(DATAPOT);
365 TH2D* hsignal_truth_energy_vs_cos_2d = (TH2D*) signal_truth_energy_vs_cos_2d->ToTH2(DATAPOT);
367 TCanvas *
c =
new TCanvas();
368 hsignal_truth_energy_vs_cos_2d->SetTitle(
"");
369 hsignal_truth_energy_vs_cos_2d->SetMarkerColor(kWhite);
370 hsignal_truth_energy_vs_cos_2d->GetYaxis()->SetRangeUser(0, 10);
371 hsignal_truth_energy_vs_cos_2d->GetXaxis()->SetRangeUser(.6, 1);
372 hsignal_truth_energy_vs_cos_2d->Draw(
"colz text");
373 c->SaveAs((plot_dump +
"/signal_truth_energy_vs_cos_2d.pdf").c_str());
375 TH1D * hsignal_reco_energy = hsignal_reco_energy_vs_cos_2d->ProjectionY(
UniqueName().c_str());
376 TH1D * hsignal_reco_cos = hsignal_reco_energy_vs_cos_2d->ProjectionX(
UniqueName().c_str());
377 hsignal_reco_energy ->SetTitle(
"");
378 hsignal_reco_cos ->SetTitle(
"");
380 TH1D * hsignal_truth_energy = hsignal_truth_energy_vs_cos_2d->ProjectionY(
UniqueName().c_str());
381 TH1D * hsignal_truth_cos = hsignal_truth_energy_vs_cos_2d->ProjectionX(
UniqueName().c_str());
382 hsignal_truth_energy ->SetTitle(
"");
383 hsignal_truth_cos ->SetTitle(
"");
385 SignalCountPlot(hsignal_reco_energy, hsignal_truth_energy, 0, 10, plot_dump +
"/signal_energy.pdf");
386 SignalCountPlot(hsignal_reco_cos , hsignal_truth_cos , 0, 10, plot_dump +
"/signal_cos.pdf" );
388 TFile *
output =
new TFile(output_file_name.c_str(),
"recreate");
389 hsignal_reco_energy_vs_cos_2d ->Write(
"hsignal_reco_energy_vs_cos_2d");
390 hsignal_truth_energy_vs_cos_2d->Write(
"hsignal_truth_energy_vs_cos_2d");
391 hreco_angular_frac ->Write(
"hreco_angular_frac");
392 htruth_angular_frac ->Write(
"htruth_angular_frac");
393 hreco_energy_frac ->Write(
"hreco_energy_frac");
394 htruth_energy_frac ->Write(
"htruth_energy_frac");
396 hreco_angular_abs ->Write(
"hreco_angular_abs");
397 htruth_angular_abs ->Write(
"htruth_angular_abs");
398 hreco_energy_abs ->Write(
"hreco_energy_abs");
399 htruth_energy_abs ->Write(
"htruth_energy_abs");
400 hstandard_axis_reco_angular_abs ->Write(
"hstandard_axis_reco_angular_abs");
401 hstandard_axis_truth_angular_abs ->Write(
"hstandard_axis_truth_angular_abs");
402 hstandard_axis_reco_energy_abs ->Write(
"hstandard_axis_reco_energy_abs");
403 hstandard_axis_truth_energy_abs ->Write(
"hstandard_axis_truth_energy_abs");
404 hsignal_reco_energy_vs_cos ->Write(
"hsignal_reco_energy_vs_cos");
405 hsignal_truth_energy_vs_cos ->Write(
"hsignal_truth_energy_vs_cos");
406 hsignal_reco_energy ->Write(
"hsignal_reco_energy");
407 hsignal_truth_energy ->Write(
"hsignal_truth_energy");
408 hsignal_reco_cos ->Write(
"hsignal_reco_cos");
409 hsignal_truth_cos ->Write(
"hsignal_truth_cos");
const std::vector< double > angedges
T max(const caf::Proxy< T > &a, T b)
const std::vector< double > eelecedges
Cuts and Vars for the 2020 FD DiF Study.
std::map< std::string, double > xmax
const ana::HistAxis kTrueElectronEStandardAxis
Float_t x1[n_points_granero]
TH2 * Rebin(TH2 *hist, std::vector< double > edges)
void ResolutionPlot(TH2 *hist, double xmin, double xmax, std::string name, bool tight=false, bool logz=false, bool plot_projections=false)
const ana::HistAxis kRecoElectronCosThetaStandardAxis("Reco cos #theta_{e}", costhetabins, kRecoElectronCosTheta)
Representation of a spectrum in any variable, with associated POT.
const ana::HistAxis kTrueElectronCosThetaStandardAxis
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
double func(double x, double y)
std::vector< float > Spectrum
const ana::HistAxis kRecoElectronEStandardAxis("Reco E_{e} (GeV)", eelecbins, kElectronE)
const double kAna2019RHCPOT
const std::vector< Binning > & GetBinnings() const
void PlotResolution(std::string input_file_name="/nova/ana/users/ddoyle/NuebarResolutionStudy/resolution_plots.root", std::string plot_dump="/nova/ana/users/ddoyle/NuebarResolutionStudy/plots")
void RMSPlot(TH2 *hist, std::vector< double > edges, double xmin, double xmax, std::string name, bool logy=false)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
std::string to_string(ModuleType mt)
void PlotProjection(TH1 *hist, std::string name, bool logy=false)
void SignalCountPlot(TH1 *reco, TH1 *truth, std::string name)
std::string UniqueName()
Return a different string each time, for creating histograms.
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
For use with Var2D.
const std::vector< std::string > & GetLabels() const