From a989e4bc432f1bbd25a75843d0d50bf9af2d5190 Mon Sep 17 00:00:00 2001
From: Jihee Kim <jihee.kim@anl.gov>
Date: Fri, 5 Mar 2021 19:51:32 +0000
Subject: [PATCH] Update pi0 analysis script
---
benchmarks/ecal/emcal_pi0s.sh | 2 +-
benchmarks/ecal/scripts/emcal_pi0.cxx | 4 +-
.../ecal/scripts/emcal_pion_anlaysis.cxx | 251 ++++++++++++++++++
3 files changed, 254 insertions(+), 3 deletions(-)
create mode 100644 benchmarks/ecal/scripts/emcal_pion_anlaysis.cxx
diff --git a/benchmarks/ecal/emcal_pi0s.sh b/benchmarks/ecal/emcal_pi0s.sh
index cdb93bbd..6dac4a0d 100644
--- a/benchmarks/ecal/emcal_pi0s.sh
+++ b/benchmarks/ecal/emcal_pi0s.sh
@@ -68,7 +68,7 @@ fi
mkdir -p results
-root -b -q "benchmarks/ecal/scripts/makeplot_pi0.C(\"${JUGGLER_REC_FILE}\")"
+root -b -q "benchmarks/ecal/scripts/emcal_pion_anlaysis.cxx(\"${JUGGLER_REC_FILE}\")"
if [[ "$?" -ne "0" ]] ; then
echo "ERROR running analysis script"
exit 1
diff --git a/benchmarks/ecal/scripts/emcal_pi0.cxx b/benchmarks/ecal/scripts/emcal_pi0.cxx
index 69018e78..3d6fd0c5 100644
--- a/benchmarks/ecal/scripts/emcal_pi0.cxx
+++ b/benchmarks/ecal/scripts/emcal_pi0.cxx
@@ -28,8 +28,8 @@ void emcal_pi0(int n_events = 1e6, const char* out_fname = "./data/emcal_pi0_0Ge
// Constraining the solid angle, but larger than that subtended by the detector
// https://indico.bnl.gov/event/7449/contributions/35966/attachments/27177/41430/EIC-DWG-Calo-03192020.pdf
// See a figure on slide 26
- double cos_theta_min = std::cos(M_PI * (120.0 / 180.0));
- double cos_theta_max = std::cos(M_PI);
+ double cos_theta_min = std::cos(M_PI * (160.0 / 180.0));
+ double cos_theta_max = std::cos(M_PI * (175.0 / 180.0));
for (events_parsed = 0; events_parsed < n_events; events_parsed++) {
// FourVector(px,py,pz,e,pdgid,status)
diff --git a/benchmarks/ecal/scripts/emcal_pion_anlaysis.cxx b/benchmarks/ecal/scripts/emcal_pion_anlaysis.cxx
new file mode 100644
index 00000000..ca1c7fa2
--- /dev/null
+++ b/benchmarks/ecal/scripts/emcal_pion_anlaysis.cxx
@@ -0,0 +1,251 @@
+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+
+#include "ROOT/RDataFrame.hxx"
+#include <iostream>
+
+#include "dd4pod/Geant4ParticleCollection.h"
+#include "dd4pod/CalorimeterHitCollection.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 emcal_pion_anlaysis(const char* input_fname = "rec_emcal_uniform_pi0s.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) {
+ //std::cout << (int) evt.size() << endl;
+ return (int) evt.size(); };
+
+ // Angle between two photons [degree]
+ auto theta_two_photons = [] (const std::vector<eic::ClusterData>& evt) {
+ std::vector<double> result;
+ double cluster_posx[2];
+ double cluster_posy[2];
+ double cluster_posz[2];
+ double cluster_energy[2];
+
+ if((int) evt.size() == 2) {
+ std::cout << "Found it" << endl;
+ auto n = 0;
+ for(const auto& i: evt) {
+ //std::cout << n << ": " << i.energy << endl;
+ cluster_posx[n] = i.position.x;
+ cluster_posy[n] = i.position.y;
+ cluster_posz[n] = i.position.z;
+ cluster_energy[n] = i.energy;
+ n++;
+ }
+ }
+ // Calculate invariant mass
+ auto dot_product_pos_clusters = cluster_posx[0]*cluster_posx[1] + cluster_posy[0]*cluster_posy[1] + cluster_posz[0]*cluster_posz[1];
+ auto mag_pos2_cluster_1 = (cluster_posx[0]*cluster_posx[0]) + (cluster_posy[0]*cluster_posy[0]) + (cluster_posz[0]*cluster_posz[0]);
+ auto mag_pos2_cluster_2 = (cluster_posx[1]*cluster_posx[1]) + (cluster_posy[1]*cluster_posy[1]) + (cluster_posz[1]*cluster_posz[1]);
+ auto cosine_clusters = (dot_product_pos_clusters/TMath::Sqrt(mag_pos2_cluster_1*mag_pos2_cluster_2));
+ auto theta_photons = TMath::ACos(cosine_clusters)*TMath::RadToDeg();
+ result.push_back(theta_photons);
+ return result;
+ };
+
+ // Invariant Mass [MeV]
+ auto inv_mass = [] (const std::vector<eic::ClusterData>& evt) {
+ std::vector<double> result;
+ double cluster_posx[2];
+ double cluster_posy[2];
+ double cluster_posz[2];
+ double cluster_energy[2];
+
+ if((int) evt.size() == 2) {
+ //std::cout << "Found it" << endl;
+ auto n = 0;
+ for(const auto& i: evt) {
+ //std::cout << n << ": " << i.energy << endl;
+ cluster_posx[n] = i.position.x;
+ cluster_posy[n] = i.position.y;
+ cluster_posz[n] = i.position.z;
+ cluster_energy[n] = i.energy;
+ n++;
+ }
+ }
+ // Calculate invariant mass
+ auto dot_product_pos_clusters = cluster_posx[0]*cluster_posx[1] + cluster_posy[0]*cluster_posy[1] + cluster_posz[0]*cluster_posz[1];
+ auto mag_pos2_cluster_1 = (cluster_posx[0]*cluster_posx[0]) + (cluster_posy[0]*cluster_posy[0]) + (cluster_posz[0]*cluster_posz[0]);
+ auto mag_pos2_cluster_2 = (cluster_posx[1]*cluster_posx[1]) + (cluster_posy[1]*cluster_posy[1]) + (cluster_posz[1]*cluster_posz[1]);
+ auto cosine_clusters = (dot_product_pos_clusters/TMath::Sqrt(mag_pos2_cluster_1*mag_pos2_cluster_2));
+ auto theta_photons = TMath::ACos(cosine_clusters)*TMath::RadToDeg();
+ auto invariant_mass = TMath::Sqrt(2.0*cluster_energy[0]*cluster_energy[1]*(1.0 - cosine_clusters))*1.e+3;
+ result.push_back(invariant_mass);
+ 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;
+ };
+
+ // Cluster Energy [GeV]
+ auto clusterE = [] (const std::vector<eic::ClusterData>& evt) {
+ std::vector<double> result;
+ auto total_eng = 0.0;
+ for (const auto& i: evt) {
+ total_eng += i.energy;
+ }
+ result.push_back(total_eng);
+ return result;
+ };
+
+ // Energy Resolution
+ 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("theta_two_photons", theta_two_photons, {"EcalClusters"})
+ .Define("inv_mass", inv_mass, {"EcalClusters"})
+ .Define("E_thr", E_thr, {"mcparticles2"})
+ .Define("clusterE", clusterE, {"EcalClusters"})
+ .Define("E_res", E_res, {"clusterE","E_thr"})
+ ;
+
+ // Define Histograms
+ auto hNCluster = d1.Histo1D({"hNCluster", "Number of Clusters; # of Clusters; Events", 5, -0.5, 4.5}, "ncluster");
+ auto hInvMass = d1.Histo1D({"hInvMass", "Invariant Mass; mass [MeV]; Events", 60,0.0,300.0}, "inv_mass");
+
+ // Select Events with One Cluster
+ auto d2 = d1.Filter("ncluster==2");
+ auto hInvMass_NC2 = d2.Histo1D({"hInvMass_NC2", "Invariant Mass; mass [MeV]; Events", 60,0.0,300.0}, "inv_mass");
+ auto hThetaTwoPhotons_NC2 = d2.Histo1D({"hThetaTwoPhotons_NC2", "Angle between Two Photons; angle [degree]; Events", 100,0.0,30.0}, "theta_two_photons");
+ auto hEres_NC2 = d2.Histo1D({"hEres_NC2", "Energy Resolution; #DeltaE/E; Events", 100, -0.5, 0.5}, "E_res");
+
+ // Select Events with Edge CUT
+
+ auto d3 = d2.Filter([=] (const std::vector<eic::ClusterData>& evt) {
+ double cluster_posx[2];
+ double cluster_posy[2];
+
+ auto n = 0;
+ for(const auto& i: evt) {
+ cluster_posx[n] = i.position.x;
+ cluster_posy[n] = i.position.y;
+ n++;
+ }
+ auto radial_dist_c1 = TMath::Sqrt(cluster_posx[0]*cluster_posx[0] + cluster_posy[0]*cluster_posy[0]);
+ auto radial_dist_c2 = TMath::Sqrt(cluster_posx[1]*cluster_posx[1] + cluster_posy[1]*cluster_posy[1]);
+
+ if (radial_dist_c1 > 18.0 && radial_dist_c1 < 30.0 && radial_dist_c2 > 18.0 && radial_dist_c2 < 30.0)
+ return true;
+ return false;}, {"EcalClusters"});
+ auto hEres_NC2_CUT = d3.Histo1D({"hEres_NC2_CUT", "Energy Resolution; #DeltaE/E; Events", 100, -0.5, 0.5}, "E_res");
+ auto hInvMass_NC2_CUT = d3.Histo1D({"hInvMass_NC2_CUT", "Invariant Mass; mass [MeV]; Events", 60,0.0,300.0}, "inv_mass");
+
+ // 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";
+
+ // 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_pi0s_ncluster.png");
+ c1->SaveAs("results/emcal_pi0s_ncluster.pdf");
+
+ TCanvas *c2 = new TCanvas("c2", "c2", 500, 500);
+ hInvMass->GetYaxis()->SetTitleOffset(1.4);
+ hInvMass->SetLineWidth(2);
+ hInvMass->SetLineColor(kBlue);
+ hInvMass->Fit("gaus");
+ hInvMass->DrawClone();
+ c2->SaveAs("results/emcal_pi0s_invariant_mass.png");
+ c2->SaveAs("results/emcal_pi0s_invariant_mass.pdf");
+
+ TCanvas *c3 = new TCanvas("c3", "c3", 500, 500);
+ hInvMass_NC2->GetYaxis()->SetTitleOffset(1.4);
+ hInvMass_NC2->SetLineWidth(2);
+ hInvMass_NC2->SetLineColor(kBlue);
+ hInvMass_NC2->Fit("gaus");
+ hInvMass_NC2->DrawClone();
+ c3->SaveAs("results/emcal_pi0s_invariant_mass_nc2.png");
+ c3->SaveAs("results/emcal_pi0s_invariant_mass_nc2.pdf");
+
+ TCanvas *c4 = new TCanvas("c4", "c4", 500, 500);
+ hThetaTwoPhotons_NC2->GetYaxis()->SetTitleOffset(1.4);
+ hThetaTwoPhotons_NC2->SetLineWidth(2);
+ hThetaTwoPhotons_NC2->SetLineColor(kBlue);
+ hThetaTwoPhotons_NC2->DrawClone();
+ c4->SaveAs("results/emcal_pi0s_angle_two_photons_nc2.png");
+ c4->SaveAs("results/emcal_pi0s_angle_two_photons_nc2.pdf");
+
+ TCanvas *c5 = new TCanvas("c5", "c5", 500, 500);
+ hEres_NC2->GetYaxis()->SetTitleOffset(1.4);
+ hEres_NC2->SetLineWidth(2);
+ hEres_NC2->SetLineColor(kBlue);
+ hEres_NC2->Fit("gaus");
+ hEres_NC2->DrawClone();
+ c5->SaveAs("results/emcal_pi0s_Eres_nc2.png");
+ c5->SaveAs("results/emcal_pi0s_Eres_nc2.pdf");
+
+ TCanvas *c6 = new TCanvas("c6", "c6", 500, 500);
+ hEres_NC2_CUT->GetYaxis()->SetTitleOffset(1.4);
+ hEres_NC2_CUT->SetLineWidth(2);
+ hEres_NC2_CUT->SetLineColor(kBlue);
+ hEres_NC2_CUT->Fit("gaus");
+ hEres_NC2_CUT->DrawClone();
+ c6->SaveAs("results/emcal_pi0s_Eres_nc2_cut.png");
+ c6->SaveAs("results/emcal_pi0s_Eres_nc2_cut.pdf");
+
+ TCanvas *c7 = new TCanvas("c7", "c7", 500, 500);
+ hInvMass_NC2_CUT->GetYaxis()->SetTitleOffset(1.4);
+ hInvMass_NC2_CUT->SetLineWidth(2);
+ hInvMass_NC2_CUT->SetLineColor(kBlue);
+ hInvMass_NC2_CUT->Fit("gaus");
+ hInvMass_NC2_CUT->DrawClone();
+ c7->SaveAs("results/emcal_pi0s_invariant_mass_nc2_cut.png");
+ c7->SaveAs("results/emcal_pi0s_invariant_mass_nc2_cut.pdf");
+
+}
--
GitLab