Update pi0 analysis script

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

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)

benchmarks/ecal/scripts/emcal_pion_anlaysis.cxx

+////////////////////////////////////////
+// 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");
+ 
+}