Select Git revision
emcal_electrons_analysis.cxx

Jihee Kim authored
emcal_electrons_analysis.cxx 6.59 KiB
////////////////////////////////////////
// Read reconstruction ROOT output file
// Plot variables
////////////////////////////////////////
#include "ROOT/RDataFrame.hxx"
#include <iostream>
#include "dd4pod/Geant4ParticleCollection.h"
#include "eicd/ClusterCollection.h"
#include "eicd/ClusterData.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TMath.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 emcal_electrons_analysis(const char* input_fname = "rec_electron_10kEvt.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);
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<float> result;
for (const auto& i: evt)
result.push_back(i.energy);
return result;
};
// Scattering Angle [degree]
auto theta = [] (const std::vector<eic::ClusterData>& evt) {
std::vector<float> result;
for(const auto& i: evt) {
auto r = TMath::Sqrt(i.position.x*i.position.x + i.position.y*i.position.y + i.position.z*i.position.z);
result.push_back(TMath::ACos(i.position.z/r)*TMath::RadToDeg());
}
return result;
};
// Pseudo-rapidity
auto eta = [] (const std::vector<eic::ClusterData>& evt) {
std::vector<float> result;
for(const auto& i: evt) {
auto r = TMath::Sqrt(i.position.x*i.position.x + i.position.y*i.position.y + i.position.z*i.position.z);
auto th = TMath::ACos(i.position.z/r)*TMath::RadToDeg();
result.push_back(-1.0 * TMath::Log(TMath::Tan((th*TMath::DegToRad())/2.0)));
}
return result;
};
// Mass [GeV]
auto mass2 = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
std::vector<float> result;
result.push_back(input[2].mass*input[2].mass);
return result;
};
// Momentum [GeV]
auto p_rec = [](const std::vector<float>& energy_term, const std::vector<float>& mass_term){return TMath::Sqrt((energy_term*energy_term)-(mass_term)); };
//auto p_rec = [] (const std::vector<eic::ClusterData>& evt) {
// std::vector<float> result;
// for (const auto& i: evt)
// result.push_back(TMath::Sqrt(i.energy*i.energy - mass2));
//return result;
//};
// Thrown Momentum [GeV]
auto p_thr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
std::vector<float> result;
result.push_back(TMath::Sqrt(input[2].psx*input[2].psx + input[2].psy*input[2].psy + input[2].psz*input[2].psz));
return result;
};
// Thrown Energy [GeV]
auto E_thr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
std::vector<float> 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;
};
// Add Acceptance and Erec-Ethr for resolution
auto d1 = d0.Define("ncluster", ncluster, {"EcalClusters"})
.Define("clusterE", clusterE, {"EcalClusters"})
.Define("theta", theta, {"EcalClusters"})
.Define("eta", eta, {"EcalClusters"})
.Define("mass2", mass2, {"mcparticles2"})
.Define("p_rec", p_rec, {"clusterE","mass2"})
.Define("p_thr", p_thr, {"mcparticles2"})
.Define("E_thr", E_thr, {"mcparticles2"})
;
// Define Histograms
auto hNCluster = d1.Histo1D({"hNCluster", "Number of Clusters; # of Clusters; Events", 5, -0.5, 4.5}, "ncluster");
auto hClusterE = d1.Histo1D({"hClusterE", "Cluster Energy; Cluster Energy [GeV]; Events", 100, -0.5, 30.5}, "clusterE");
auto hTheta = d1.Histo1D({"hTheta", "Scattering Angle; #theta [degree]; Events", 100, 130.0, 180.0}, "theta");
auto hEta = d1.Histo1D({"hEta", "Pseudo-rapidity; #eta; Events", 100, -5.0, 0.0}, "eta");
auto hPrec = d1.Histo1D({"hPrec", "Reconstructed Momentum; p_{rec}[GeV]; Events", 100, -0.5, 30.5}, "p_rec");
auto hPthr = d1.Histo1D({"hPthr", "Thrown Momentum; p_{thr}[GeV]; Events", 100, -0.5, 30.5}, "p_thr");
auto hEthr = d1.Histo1D({"hEthr", "Thrown Energy; E_{thr}[GeV]; Events", 100, -0.5, 30.5}, "E_thr");
// Draw Histograms
TCanvas *c1 = new TCanvas("c1", "c1", 500, 500);
c1->SetLogy(1);
hNCluster->GetYaxis()->SetTitleOffset(1.6);
hNCluster->SetLineWidth(2);
hNCluster->SetLineColor(kBlue);
hNCluster->DrawClone();
c1->SaveAs("results/emcal_electrons_ncluster.png");
c1->SaveAs("results/emcal_electrons_ncluster.pdf");
TCanvas *c2 = new TCanvas("c2", "c2", 500, 500);
hClusterE->GetYaxis()->SetTitleOffset(1.4);
hClusterE->SetLineWidth(2);
hClusterE->SetLineColor(kBlue);
hClusterE->DrawClone();
c2->SaveAs("results/emcal_electrons_clusterE.png");
c2->SaveAs("results/emcal_electrons_clusterE.pdf");
TCanvas *c3 = new TCanvas("c3", "c3", 500, 500);
hTheta->GetYaxis()->SetTitleOffset(1.4);
hTheta->SetLineWidth(2);
hTheta->SetLineColor(kBlue);
hTheta->DrawClone();
c3->SaveAs("results/emal_electrons_theta.png");
c3->SaveAs("results/emal_electrons_theta.pdf");
TCanvas *c4 = new TCanvas("c4", "c4", 500, 500);
hEta->GetYaxis()->SetTitleOffset(1.4);
hEta->SetLineWidth(2);
hEta->SetLineColor(kBlue);
hEta->DrawClone();
c4->SaveAs("results/emcal_electrons_eta.png");
c4->SaveAs("results/emcal_electrons_eta.pdf");
TCanvas *c5 = new TCanvas("c5", "c5", 500, 500);
hPrec->GetYaxis()->SetTitleOffset(1.4);
hPrec->SetLineWidth(2);
hPrec->SetLineColor(kBlue);
hPrec->DrawClone();
c5->SaveAs("results/emcal_electrons_Prec.png");
c5->SaveAs("results/emcal_electrons_Prec.pdf");
TCanvas *c6 = new TCanvas("c6", "c6", 500, 500);
hPthr->GetYaxis()->SetTitleOffset(1.4);
hPthr->SetLineWidth(2);
hPthr->SetLineColor(kBlue);
hPthr->DrawClone();
c6->SaveAs("results/emcal_electrons_Pthr.png");
c6->SaveAs("results/emcal_electrons_Pthr.pdf");
TCanvas *c7 = new TCanvas("c7", "c7", 500, 500);
hEthr->GetYaxis()->SetTitleOffset(1.4);
hEthr->SetLineWidth(2);
hEthr->SetLineColor(kBlue);
hEthr->DrawClone();
c7->SaveAs("results/emcal_electrons_Ethr.png");
c7->SaveAs("results/emcal_electrons_Ethr.pdf");
//h->Fit("gaus");
}