Skip to content
Snippets Groups Projects
Commit 14155fff authored by Jihee Kim's avatar Jihee Kim
Browse files

Eres fit script into CI pipeline

parent 0b24bac9
No related branches found
No related tags found
1 merge request!53Eres fit script into CI pipeline
Hide whitespace changes
Inline Side-by-side
ecal/emcal_electrons.sh
+ 7
− 0
View file @ 14155fff
......@@ -104,6 +104,13 @@ if [[ "$?" -ne "0" ]] ; then
exit 1
fi
# Script to generate energy resolution plots
root -b -q "ecal/scripts/rec_emcal_resolution_analysis.cxx(\"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\")"
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running analysis script"
exit 1
fi
#paste results/sim_${JUGGLER_FILE_NAME_TAG}.txt results/rec_${JUGGLER_FILE_NAME_TAG}.txt > results/eng_${JUGGLER_FILE_NAME_TAG}.txt
#root -b -q "ecal/scripts/read_eng.C(\"results/eng_${JUGGLER_FILE_NAME_TAG}.root\", \"results/eng_${JUGGLER_FILE_NAME_TAG}.txt\")"
#root -b -q "ecal/scripts/cal_eng_res.C(\"results/eng_${JUGGLER_FILE_NAME_TAG}.root\")"
......
ecal/scripts/rec_emcal_resolution_analysis.cxx
+ 108
− 97
View file @ 14155fff
......@@ -7,8 +7,8 @@
#include "ROOT/RDataFrame.hxx"
#include <iostream>
#include "dd4pod/Geant4ParticleCollection.h"
#include "dd4pod/CalorimeterHitCollection.h"
#include "dd4pod/Geant4ParticleCollection.h"
#include "dd4pod/TrackerHitCollection.h"
#include "eicd/CalorimeterHitCollection.h"
#include "eicd/CalorimeterHitData.h"
......@@ -16,10 +16,10 @@
#include "eicd/ClusterData.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TMath.h"
#include "TH1.h"
#include "TH1D.h"
#include "TMath.h"
#include "TStyle.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
......@@ -45,64 +45,69 @@ void rec_emcal_resolution_analysis(const char* input_fname = "rec_emcal_uniform_
ROOT::RDataFrame d0("events", input_fname);
// Number of Clusters
auto ncluster = [] (const std::vector<eic::ClusterData>& evt) {return (int) evt.size(); };
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;
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;
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
// 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;
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("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"});
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");
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
// Do gaussian fit 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
......@@ -123,106 +128,110 @@ void rec_emcal_resolution_analysis(const char* input_fname = "rec_emcal_uniform_
// 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");
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","",h->GetMean() - 2.0 * h->GetRMS(),h->GetMean() + 2.0 * h->GetRMS());
// 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);
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->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);
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->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);
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->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
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);
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->GetYaxis()->SetRangeUser(0.0, 1.5);
gCEVsSigmaENC1EGCUT->SetMarkerStyle(21);
gCEVsSigmaENC1EGCUT->Fit("f1","W,E");
gCEVsSigmaENC1EGCUT->Fit("f1", "W");
gCEVsSigmaENC1EGCUT->DrawClone("AEP");
c4->SaveAs("results/emcal_electrons_gCEVsSimgaENC1EGCUT.png");
c4->SaveAs("results/emcal_electrons_gCEVsSigmaENC1EGCUT.pdf");
......@@ -232,6 +241,8 @@ void rec_emcal_resolution_analysis(const char* input_fname = "rec_emcal_uniform_
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";
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";
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment