Skip to content
Snippets Groups Projects

Energy Resolution Fit

Closed Energy Resolution Fit
jihee.kim/reconstruction_benchmarks:Eres_fit into master
Closed Jihee Kim requested to merge jihee.kim/reconstruction_benchmarks:Eres_fit into master
Compare
and
1 file
+ 237
− 0
Compare changes
  • Side-by-side
  • Inline
ecal/scripts/rec_emcal_resolution_analysis.cxx 0 → 100644
+ 237
− 0
 
////////////////////////////////////////
 
// Read reconstruction ROOT output file
 
// Plot energy resolution
 
// 2021/01/20
 
////////////////////////////////////////
 
 
#include "ROOT/RDataFrame.hxx"
 
#include <iostream>
 
 
#include "dd4pod/Geant4ParticleCollection.h"
 
#include "dd4pod/CalorimeterHitCollection.h"
 
#include "dd4pod/TrackerHitCollection.h"
 
#include "eicd/CalorimeterHitCollection.h"
 
#include "eicd/CalorimeterHitData.h"
 
#include "eicd/ClusterCollection.h"
 
#include "eicd/ClusterData.h"
 
 
#include "TCanvas.h"
 
#include "TStyle.h"
 
#include "TMath.h"
 
#include "TH1.h"
 
#include "TH1D.h"
 
 
// If root cannot find a dd4pod header make sure to add the include dir to ROOT_INCLUDE_PATH:
 
// export ROOT_INCLUDE_PATH=$HOME/include:$ROOT_INCLUDE_PATH
 
// export ROOT_INCLUDE_PATH=/home/jihee/stow/development/include:$ROOT_INCLUDE_PATH
 
 
using ROOT::RDataFrame;
 
using namespace ROOT::VecOps;
 
 
void rec_emcal_resolution_analysis(const char* input_fname = "rec_emcal_uniform_electrons.root")
 
{
 
// Setting for graphs
 
gROOT->SetStyle("Plain");
 
gStyle->SetOptFit(1);
 
gStyle->SetLineWidth(2);
 
gStyle->SetPadTickX(1);
 
gStyle->SetPadTickY(1);
 
gStyle->SetPadGridX(1);
 
gStyle->SetPadGridY(1);
 
gStyle->SetPadLeftMargin(0.14);
 
gStyle->SetPadRightMargin(0.14);
 
 
ROOT::EnableImplicitMT();
 
ROOT::RDataFrame d0("events", input_fname);
 
 
// Number of Clusters
 
auto ncluster = [] (const std::vector<eic::ClusterData>& evt) {return (int) evt.size(); };
 
 
// Cluster Energy [GeV]
 
auto clusterE = [] (const std::vector<eic::ClusterData>& evt) {
 
std::vector<double> result;
 
for (const auto& i: evt)
 
result.push_back(i.energy);
 
return result;
 
};
 
 
// Thrown Energy [GeV]
 
auto E_thr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
 
std::vector<double> result;
 
result.push_back(TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz + input[2].mass*input[2].mass));
 
return result;
 
};
 
 
// Energy Resolution = delta(E)/E
 
// delta(E) = (Erec - Ethr)
 
// E = Ethr
 
auto E_res = [] (const std::vector<double>& Erec, const std::vector<double>& Ethr) {
 
std::vector<double> result;
 
for (const auto& E1 : Ethr) {
 
for (const auto& E2 : Erec) {
 
result.push_back((E2-E1)/E1);
 
}
 
}
 
return result;
 
};
 
 
// Define variables
 
auto d1 = d0.Define("ncluster", ncluster, {"EcalClusters"})
 
.Define("clusterE", clusterE, {"EcalClusters"})
 
.Define("E_thr", E_thr, {"mcparticles2"})
 
.Define("E_res", E_res, {"clusterE","E_thr"})
 
;
 
 
// Define Cuts
 
// d2: Select Events with One Cluster
 
// d3: Select Events which determined their center of cluster not being on the edge of detector in addition to d2 cut
 
auto d2 = d1.Filter("ncluster==1");
 
auto d3 = d2.Filter([=] (const std::vector<eic::ClusterData>& evt) {
 
for(const auto& i: evt) {
 
auto pos_x = i.position.x;
 
auto pos_y = i.position.y;
 
auto radial_dist = TMath::Sqrt(pos_x*pos_x + pos_y*pos_y);
 
if (radial_dist > 8.0 && radial_dist < 30.0)
 
return true;
 
}
 
return false;}, {"EcalClusters"});
 
 
// Define histograms
 
auto hEresNC1 = d2.Histo1D({"hEresNC1", "Energy Resolution CE vs Ethr w/ NC1; #DeltaE/E; Events", 100, -0.5, 0.5}, "E_res");
 
auto hEresNC1EGCUT = d3.Histo1D({"hEresNC1EGCUT", "Energy Resolution CE vs Ethr w/ NC1 EGCUT; #DeltaE/E; Events", 100, -0.5, 0.5}, "E_res");
 
 
// Energy Resolution Graphs
 
// Divided into serveral energy bin group
 
// Do gaussian git to extract mean and sigma
 
int nbin = 12; // number of energy bin group
 
double maxEnergy = 30.0; // maximum of energy in dataset
 
double multiple = maxEnergy / (double)nbin; // energy interval between energy bin group
 
double maximumE; // high end in each energy bin group
 
double minimumE; // low end in each energy bin group
 
double mean[nbin]; // mean from gaussian fit
 
double sigma[nbin]; // sigma from gaussian fit
 
double dmean[nbin]; // error on mean
 
double dsigma[nbin]; // error on sigma
 
double gEbin[nbin]; // center of energy in each group
 
double invSqrtgEbin[nbin]; // inverse sqrt center of energy in each group
 
double resolution[nbin]; // resolution = sigma/mean
 
double resolutionError[nbin]; // error on resolution
 
double resolutionErrorSigma; // error on resolution sigma term
 
double resolutionErrorMean; // error on resolution mean term
 
double xerror[nbin]; // set to 0.0 for all
 
// Loop over each energy bin group
 
// Define hitogram and do gaussian fit
 
// Save fit results; mean, sigma, error on mean & sigam
 
// Calculate resolution and resolution error
 
for(int n=1; n <= nbin; n++) {
 
maximumE = n * multiple;
 
minimumE = maximumE - multiple;
 
gEbin[n-1] = (maximumE + minimumE) / 2.0;
 
invSqrtgEbin[n-1] = 1.0 / TMath::Sqrt(gEbin[n-1]);
 
// Define histogram
 
auto h = d3.Filter([=] (const std::vector<eic::ClusterData>& evt) {
 
for (const auto& i: evt) {
 
auto eng = i.energy;
 
if (eng >= minimumE && eng < maximumE)
 
return true;
 
}
 
return false;}, {"EcalClusters"})
 
.Histo1D({"h", "Energy Resolution; #DeltaE/E; Events", 100, -0.5, 0.5}, "E_res");
 
// Gaussian fit
 
h->Fit("gaus","W,E");
 
// Save fit results
 
// how to access how many events in histogram: h->GetEntries();
 
if(h->GetFunction("gaus")) {
 
mean[n-1] = TMath::Abs(h->GetFunction("gaus")->GetParameter(1));
 
dmean[n-1] = h->GetFunction("gaus")->GetParError(1);
 
sigma[n-1] = h->GetFunction("gaus")->GetParameter(2);
 
dsigma[n-1] = h->GetFunction("gaus")->GetParError(2);
 
resolution[n-1] = sigma[n-1] / mean[n-1];
 
resolutionErrorSigma = dsigma[n-1] / mean[n-1];
 
resolutionErrorMean = dmean[n-1] * sigma[n-1];
 
resolutionError[n-1] = TMath::Sqrt(resolutionErrorSigma*resolutionErrorSigma + resolutionErrorMean*resolutionErrorMean);
 
xerror[n-1] = 0.0;
 
}
 
// Draw histogram, but
 
// Just save fit results
 
//TCanvas *c = new TCanvas("c", "c", 500, 500);
 
//h->GetYaxis()->SetTitleOffset(1.4);
 
//h->SetLineWidth(2);
 
//h->SetLineColor(kBlue);
 
//h->Fit("gaus","W,E");
 
//h->DrawClone();
 
//c->SaveAs("results/h.png");
 
}
 
 
// Generate graphs using gaussian fit results
 
// Cluster energy VS Simga
 
TGraphErrors* gCEVsSigmaNC1EGCUT = new TGraphErrors(nbin,gEbin,sigma,xerror,dsigma);
 
TCanvas *c1 = new TCanvas("c1","c1",500,500);
 
gCEVsSigmaNC1EGCUT->SetTitle("Cluster Energy vs #sigma");
 
gCEVsSigmaNC1EGCUT->GetYaxis()->SetTitleOffset(1.4);
 
gCEVsSigmaNC1EGCUT->GetXaxis()->SetTitle("Cluster Energy [GeV]");
 
gCEVsSigmaNC1EGCUT->GetYaxis()->SetTitle("#sigma");
 
gCEVsSigmaNC1EGCUT->GetYaxis()->SetRangeUser(0.0,0.05);
 
gCEVsSigmaNC1EGCUT->SetMarkerStyle(21);
 
gCEVsSigmaNC1EGCUT->DrawClone("AEP");
 
c1->SaveAs("results/emcal_electrons_gCEVsSigmaNC1EGCUT.png");
 
c1->SaveAs("results/emcal_electrons_gCEVsSigmaNC1EGCUT.pdf");
 
// Cluster energy VS Mean
 
TGraphErrors* gCEVsMeanNC1EGCUT = new TGraphErrors(nbin,gEbin,mean,xerror,dmean);
 
TCanvas *c2 = new TCanvas("c2","c2",500,500);
 
gCEVsMeanNC1EGCUT->SetTitle("Cluster Energy vs #mu");
 
gCEVsMeanNC1EGCUT->GetYaxis()->SetTitleOffset(1.4);
 
gCEVsMeanNC1EGCUT->GetXaxis()->SetTitle("Cluster Energy [GeV]");
 
gCEVsMeanNC1EGCUT->GetYaxis()->SetTitle("#mu");
 
gCEVsMeanNC1EGCUT->GetYaxis()->SetRangeUser(0.0,0.1);
 
gCEVsMeanNC1EGCUT->SetMarkerStyle(21);
 
gCEVsMeanNC1EGCUT->DrawClone("AEP");
 
c2->SaveAs("results/emcal_electrons_gCEVsMeanNC1EGCUT.png");
 
c2->SaveAs("results/emcal_electrons_gCEVsMeanNC1EGCUT.pdf");
 
// Inverse of sqrt(cluster energy) VS resolution
 
TGraphErrors* gInvSqCEVsSigmaENC1EGCUT = new TGraphErrors(nbin,invSqrtgEbin,resolution,xerror,resolutionError);
 
TCanvas *c3 = new TCanvas("c3","c3",500,500);
 
gInvSqCEVsSigmaENC1EGCUT->SetTitle("#sigma/E [%] vs 1/sqrt(Cluster Energy)");
 
gInvSqCEVsSigmaENC1EGCUT->GetYaxis()->SetTitleOffset(1.4);
 
gInvSqCEVsSigmaENC1EGCUT->GetXaxis()->SetTitle("1/sqrt(Cluster Energy)");
 
gInvSqCEVsSigmaENC1EGCUT->GetYaxis()->SetTitle("#sigma/E [%]");
 
gInvSqCEVsSigmaENC1EGCUT->GetYaxis()->SetRangeUser(0.0,1.5);
 
gInvSqCEVsSigmaENC1EGCUT->GetXaxis()->SetRangeUser(0.0,1.5);
 
gInvSqCEVsSigmaENC1EGCUT->SetMarkerStyle(21);
 
gInvSqCEVsSigmaENC1EGCUT->DrawClone("AEP");
 
c3->SaveAs("results/emcal_electrons_gInvSqCEVsSimgaENC1EGCUT.png");
 
c3->SaveAs("results/emcal_electrons_gInvSqCEVsSigmaENC1EGCUT.pdf");
 
 
// Fit Function
 
// sigma/E = S/sqrt(E) + N/E + C; added in quadrature
 
// S: Stochastic term
 
// N: Niose term
 
// C: Constant term
 
TF1 *f1 = new TF1("f1","sqrt([0]*[0] + pow([1]/sqrt(x),2) + pow([2]/x,2))", 0.5, 30.0);
 
f1->SetParNames("C", "S", "N"); // param names
 
f1->SetParameters(0.30,2.80,0.12); // initial param values
 
f1->SetLineWidth(2);
 
f1->SetLineColor(kRed);
 
 
// Cluster energy VS resolution
 
TGraphErrors* gCEVsSigmaENC1EGCUT = new TGraphErrors(nbin,gEbin,resolution,xerror,resolutionError);
 
TCanvas *c4 = new TCanvas("c4","c4",500,500);
 
gCEVsSigmaENC1EGCUT->SetTitle("Cluster Energy vs #sigma/E [%]");
 
gCEVsSigmaENC1EGCUT->GetYaxis()->SetTitleOffset(1.4);
 
gCEVsSigmaENC1EGCUT->GetXaxis()->SetTitle("Cluster Energy [GeV]");
 
gCEVsSigmaENC1EGCUT->GetYaxis()->SetTitle("#sigma/E [%]");
 
gCEVsSigmaENC1EGCUT->GetYaxis()->SetRangeUser(0.0,1.5);
 
gCEVsSigmaENC1EGCUT->SetMarkerStyle(21);
 
gCEVsSigmaENC1EGCUT->Fit("f1","W,E");
 
gCEVsSigmaENC1EGCUT->DrawClone("AEP");
 
c4->SaveAs("results/emcal_electrons_gCEVsSimgaENC1EGCUT.png");
 
c4->SaveAs("results/emcal_electrons_gCEVsSigmaENC1EGCUT.pdf");
 
 
// Event Counts
 
auto nevents_thrown = d1.Count();
 
auto nevents_cluster1 = d2.Count();
 
auto nevents_cluster1cut = d3.Count();
 
 
std::cout << "Number of Events: " << (*nevents_cluster1) << " / " << (*nevents_thrown) << " = single cluster events/all \n";
 
std::cout << "Number of Events: " << (*nevents_cluster1cut) << " / " << (*nevents_thrown) << " = single cluster events that passed edge cut/all \n";
 
}