nhcal_basic_distribution/scripts/basic_distribution_energy_resolution.cxx

0001 #include <cmath>
0002 #include <fstream>
0003 #include <iostream>
0004 #include <string>
0005 #include <vector>
0006
0007 #include <TH1D.h>
0008 #include <TH2D.h>
0009 #include <TFile.h>
0010 #include <TMath.h>
0011 #include <TCanvas.h>
0012 #include <TStyle.h>
0013 #include <TGraphErrors.h>
0014 #include <TF1.h>
0015 #include <TProfile.h>
0016 #include <TLegend.h>
0017
0018 #include "TROOT.h"
0019
0020 using namespace std;
0021 using namespace TMath;
0022
0023 string changeExtension(const string& path, const string& new_ext)
0024 {
0025     size_t pos = path.find_last_of('.');
0026     if (pos != string::npos)
0027         return path.substr(0, pos) + new_ext;
0028     return path + new_ext;
0029 }
0030
0031 int basic_distribution_energy_resolution(const string &filename, string outname_png){
0032
0033     gStyle->SetTitleSize(0.045, "XYZ");
0034     gStyle->SetLabelSize(0.04, "XYZ");
0035     gStyle->SetPadLeftMargin(0.25);
0036     gStyle->SetPadRightMargin(0.15);
0037     gStyle->SetPadBottomMargin(0.15);
0038     gStyle->SetPadTopMargin(0.10);
0039
0040     string outname_pdf = changeExtension(outname_png, ".pdf");
0041
0042     TFile *file = TFile::Open(filename.c_str(), "READ");
0043     if (!file || file->IsZombie())
0044     {
0045         cerr << "Cannot open file: " << filename << endl;
0046         return 1;
0047     }
0048
0049     TH2D *h_energyRes = (TH2D*)file->Get("h_energyRes");
0050     TProfile *p_energyRes = (TProfile*)file->Get("p_energyRes");
0051     TH1D *h_nHCal_hit_contrib_time   = (TH1D*)file->Get("h_nHCal_hit_contrib_time");
0052     TH1D *h_nHCal_hit_contrib_energy = (TH1D*)file->Get("h_nHCal_hit_contrib_energy");
0053     TH2D *h_nHCal_hit_contrib_energy_vs_time       = (TH2D*)file->Get("h_nHCal_hit_contrib_energy_vs_time");
0054     TH2D *h_nHCal_hit_contrib_energy_vs_telap      = (TH2D*)file->Get("h_nHCal_hit_contrib_energy_vs_telap");
0055     TH2D *h_nHCal_hit_contrib_energy_vs_time_total = (TH2D*)file->Get("h_nHCal_hit_contrib_energy_vs_time_total");
0056
0057     if (!h_energyRes || !p_energyRes || !h_nHCal_hit_contrib_time || !h_nHCal_hit_contrib_energy ||
0058     !h_nHCal_hit_contrib_energy_vs_time || !h_nHCal_hit_contrib_energy_vs_telap ||
0059     !h_nHCal_hit_contrib_energy_vs_time_total)
0060     {
0061         cerr << "Cannot retrieve histograms from file" << endl;
0062         file->Close();
0063         return 1;
0064     }
0065
0066     TGraphErrors *g_resolution = new TGraphErrors();
0067
0068     for (int i = 1; i <= p_energyRes->GetNbinsX(); i++) {
0069         double Ekin = p_energyRes->GetBinCenter(i);
0070         double mean = p_energyRes->GetBinContent(i);
0071         double rms = p_energyRes->GetBinError(i);
0072         int entries = p_energyRes->GetBinEntries(i);
0073
0074         if (entries < 10 || mean <= 0) continue;
0075
0076         int n = g_resolution->GetN();
0077         g_resolution->SetPoint(n, Ekin, rms/mean);
0078         g_resolution->SetPointError(n, 0, rms/mean * sqrt(1.0/entries));
0079     }
0080
0081     TF1 *fit = new TF1("fit", "[0]*pow(x, -0.5) + [1] + [2]*pow(x, 1)", 0.1, 6);
0082     fit->SetParameters(0.5, 0.05, 0.001);
0083     fit->SetParNames("a (stochastic)", "b (constant)", "c (noise)");
0084     g_resolution->Fit(fit, "R");
0085
0086     TCanvas *c = new TCanvas("c", "Energy Resolution", 1600, 800);
0087     c->Divide(2, 1);
0088     c->cd(1);
0089
0090     h_energyRes->Draw("COLZ");
0091     p_energyRes->SetLineWidth(3);
0092     p_energyRes->SetLineColor(kRed);
0093     p_energyRes->SetMarkerColor(kRed);
0094     p_energyRes->Draw("SAME");
0095
0096     c->cd(2);
0097     g_resolution->SetTitle("Energy Resolution;E_{kin} [GeV];#sigma/E");
0098     g_resolution->SetMarkerStyle(20);
0099     g_resolution->SetLineWidth(2);
0100     g_resolution->SetLineColor(kBlue);
0101     g_resolution->Draw("APE");
0102     fit->Draw("SAME");
0103
0104     double par_a = fit->GetParameter(0);
0105     double par_b = fit->GetParameter(1);
0106     double par_c = fit->GetParameter(2);
0107
0108     TLegend *leg = new TLegend(0.35, 0.65, 0.80, 0.80);
0109     leg->SetTextSize(0.025);
0110     leg->AddEntry(g_resolution, "Data", "p");
0111     leg->AddEntry(fit, Form("Fit: %.3f/#sqrt{E} + %.3f + %.3f#timesE", par_a, par_b, par_c), "l");
0112     leg->Draw();
0113
0114     c->SaveAs(outname_png.c_str());
0115     c->SaveAs(outname_pdf.c_str());
0116
0117     TCanvas *c_neutronThresholds = new TCanvas("c_neutronThresholds", "c_neutronThresholds", 1600, 800);
0118     c_neutronThresholds->Divide(3, 2);
0119     c_neutronThresholds->cd(1); h_nHCal_hit_contrib_time->Draw();
0120     c_neutronThresholds->cd(2); h_nHCal_hit_contrib_energy->Draw();
0121     c_neutronThresholds->cd(3); h_nHCal_hit_contrib_energy_vs_time->Draw("COLZ");
0122     c_neutronThresholds->cd(4); h_nHCal_hit_contrib_energy_vs_telap->Draw("COLZ");
0123     c_neutronThresholds->cd(5); h_nHCal_hit_contrib_energy_vs_time_total->Draw("COLZ");
0124
0125     auto replaceStr = [](string& s, const string& from, const string& to) {
0126         size_t pos = s.find(from);
0127         if (pos != string::npos) s.replace(pos, from.size(), to);
0128     };
0129
0130     string neutron_png = outname_png;
0131     string neutron_pdf = changeExtension(outname_png, ".pdf");
0132
0133     replaceStr(neutron_png, "energy_resolution_analysis", "neutronThresholds");
0134     replaceStr(neutron_pdf, "energy_resolution_analysis", "neutronThresholds");
0135
0136     c_neutronThresholds->SaveAs(neutron_png.c_str());
0137     c_neutronThresholds->SaveAs(neutron_pdf.c_str());
0138     delete c_neutronThresholds;
0139
0140     file->Close();
0141
0142     return 0;
0143 }