From 6bd3562e833329b571a61eba858f29da48320834 Mon Sep 17 00:00:00 2001
From: Jihee Kim <jihee.kim@anl.gov>
Date: Wed, 21 Oct 2020 16:12:32 -0500
Subject: [PATCH] Created unified script to read reconstructed and thrown
information and Added acceptance calcuation
Debugging
Debugging
Debugging
Debugging - Removed mcparticles2 and CrystalEcalHits2 from reading
Debugging - more CrystalEcalHits2 info removed
---
ecal/emcal_electrons.sh | 12 +-
ecal/scripts/rec_emcal_electrons_reader.C | 370 ++++++++++++++++++++++
2 files changed, 377 insertions(+), 5 deletions(-)
create mode 100644 ecal/scripts/rec_emcal_electrons_reader.C
diff --git a/ecal/emcal_electrons.sh b/ecal/emcal_electrons.sh
index b8af5703..271f8e9e 100644
--- a/ecal/emcal_electrons.sh
+++ b/ecal/emcal_electrons.sh
@@ -55,12 +55,14 @@ popd
pwd
mkdir -p results
-root -b -q "ecal/scripts/makeplot.C(${E_start}, ${E_end}, \"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\", \"results/rec_${JUGGLER_FILE_NAME_TAG}.txt\")"
-root -b -q "ecal/scripts/makeplot_input.C(\"${JUGGLER_DETECTOR}/${JUGGLER_SIM_FILE}\", \"results/sim_${JUGGLER_FILE_NAME_TAG}.txt\")"
+#rootls ${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}
+root -b -q "ecal/scripts/rec_emcal_electrons_reader.C(${E_start}, ${E_end}, \"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\")"
+#root -b -q "ecal/scripts/makeplot.C(${E_start}, ${E_end}, \"${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}\", \"results/rec_${JUGGLER_FILE_NAME_TAG}.txt\")"
+#root -b -q "ecal/scripts/makeplot_input.C(\"${JUGGLER_DETECTOR}/${JUGGLER_SIM_FILE}\", \"results/sim_${JUGGLER_FILE_NAME_TAG}.txt\")"
-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\")"
+#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\")"
#mkdir -p sim_output
#cp "${JUGGLER_DETECTOR}/${JUGGLER_REC_FILE}" sim_output/.
diff --git a/ecal/scripts/rec_emcal_electrons_reader.C b/ecal/scripts/rec_emcal_electrons_reader.C
new file mode 100644
index 00000000..38d59ba3
--- /dev/null
+++ b/ecal/scripts/rec_emcal_electrons_reader.C
@@ -0,0 +1,370 @@
+////////////////////////////////////////
+// Read reconstruction ROOT output file
+// Plot variables
+////////////////////////////////////////
+int rec_emcal_electrons_reader(double e_start = 1.0, double e_end = 1.0, const char* input_fname = "sim_output/sim_emcal_electrons_output.root")
+{
+ // Setting figures
+ gROOT->SetStyle("Plain");
+ gStyle->SetLineWidth(3);
+ gStyle->SetOptStat("nem");
+ gStyle->SetPadTickX(1);
+ gStyle->SetPadTickY(1);
+ gStyle->SetPadGridX(1);
+ gStyle->SetPadGridY(1);
+ gStyle->SetPadLeftMargin(0.14);
+ gStyle->SetPadRightMargin(0.14);
+
+ // Input ROOT file
+ TFile *f = new TFile(input_fname,"READ");
+ TTree *t = (TTree *)f->Get("events");
+
+ // Set Branch status and addressed
+ t->SetMakeClass(1);
+ t->SetBranchStatus("*", 0);
+/*
+ Int_t mcparticles2_;
+ t->SetBranchStatus("mcparticles2", 1);
+ t->SetBranchAddress("mcparticles2", &mcparticles2_);
+ Int_t CrystalEcalHits2_;
+ t->SetBranchStatus("CrystalEcalHits2", 1);
+ t->SetBranchAddress("CrystalEcalHits2", &CrystalEcalHits2_);
+*/
+ Int_t RecoEcalHits_;
+ t->SetBranchStatus("RecoEcalHits", 1);
+ t->SetBranchAddress("RecoEcalHits", &RecoEcalHits_);
+ Int_t EcalClusters_;
+ t->SetBranchStatus("EcalClusters", 1);
+ t->SetBranchAddress("EcalClusters", &EcalClusters_);
+/*
+ const Int_t kMaxmcparticles2 = 100000;
+ Double_t px[kMaxmcparticles2];
+ Double_t py[kMaxmcparticles2];
+ Double_t pz[kMaxmcparticles2];
+ Double_t mass[kMaxmcparticles2];
+ t->SetBranchStatus("mcparticles2.psx",1);
+ t->SetBranchStatus("mcparticles2.psy",1);
+ t->SetBranchStatus("mcparticles2.psz",1);
+ t->SetBranchStatus("mcparticles2.mass",1);
+ t->SetBranchAddress("mcparticles2.psx",px);
+ t->SetBranchAddress("mcparticles2.psy",py);
+ t->SetBranchAddress("mcparticles2.psz",pz);
+ t->SetBranchAddress("mcparticles2.mass",mass);
+ const Int_t kMaxCrystalEcalHits2 = 100000;
+ Double_t truth_deposit[kMaxCrystalEcalHits2];
+ Double_t energyDeposit[kMaxCrystalEcalHits2];
+ t->SetBranchStatus("CrystalEcalHits2.truth.deposit",1);
+ t->SetBranchStatus("CrystalEcalHits2.energyDeposit",1);
+ t->SetBranchAddress("CrystalEcalHits2.truth.deposit",truth_deposit);
+ t->SetBranchAddress("CrystalEcalHits2.energyDeposit",energyDeposit);
+*/
+ const Int_t kMaxRecoEcalHits = 100000;
+ Double_t rec_x_pos[kMaxRecoEcalHits];
+ Double_t rec_y_pos[kMaxRecoEcalHits];
+ Double_t rec_energy[kMaxRecoEcalHits];
+ t->SetBranchStatus("RecoEcalHits.position.x",1);
+ t->SetBranchStatus("RecoEcalHits.position.y",1);
+ t->SetBranchStatus("RecoEcalHits.energy",1);
+ t->SetBranchAddress("RecoEcalHits.position.x",rec_x_pos);
+ t->SetBranchAddress("RecoEcalHits.position.y",rec_y_pos);
+ t->SetBranchAddress("RecoEcalHits.energy",rec_energy);
+ const Int_t kMaxEcalClusters = 100000;
+ Double_t cluster_x_pos[kMaxEcalClusters];
+ Double_t cluster_y_pos[kMaxEcalClusters];
+ Double_t cluster_z_pos[kMaxEcalClusters];
+ Float_t cluster_energy[kMaxEcalClusters];
+ t->SetBranchStatus("EcalClusters.position.x",1);
+ t->SetBranchStatus("EcalClusters.position.y",1);
+ t->SetBranchStatus("EcalClusters.position.z",1);
+ t->SetBranchStatus("EcalClusters.energy",1);
+ t->SetBranchAddress("EcalClusters.position.x",cluster_x_pos);
+ t->SetBranchAddress("EcalClusters.position.y",cluster_y_pos);
+ t->SetBranchAddress("EcalClusters.position.z",cluster_z_pos);
+ t->SetBranchAddress("EcalClusters.energy",cluster_energy);
+
+ // Setting for Canvas
+ TCanvas *c1 = new TCanvas("c1", "c1", 500, 500);
+ TCanvas *c2 = new TCanvas("c2", "c2", 500, 500);
+ TCanvas *c3 = new TCanvas("c3", "c3", 500, 500);
+ TCanvas *c4 = new TCanvas("c4", "c4", 500, 500);
+ TCanvas *c5 = new TCanvas("c5", "c5", 500, 500);
+ TCanvas *c6 = new TCanvas("c6", "c6", 500, 500);
+ TCanvas *c7 = new TCanvas("c7", "c7", 500, 500);
+ TCanvas *c8 = new TCanvas("c8", "c8", 500, 500);
+ TCanvas *c9 = new TCanvas("c9", "c9", 500, 500);
+ TCanvas *c10 = new TCanvas("c10", "c10", 700, 500);
+ TCanvas *c11 = new TCanvas("c11", "c11", 700, 500);
+ TCanvas *c12 = new TCanvas("c12", "c12", 700, 500);
+ TCanvas *c13 = new TCanvas("c13", "c13", 700, 500);
+ TCanvas *c14 = new TCanvas("c14", "c14", 700, 500);
+
+ // Declare histograms
+ TH1D *h1 = new TH1D("h1","Scattering Angle(#theta)", 100,130.0,180.0);
+ TH1D *h2 = new TH1D("h2","Pseudo-rapidity(#eta)", 100,-5.0,0.0);
+ TH2D *h3 = new TH2D("h3","Cluster E vs Pseudo-rapidity", 100,e_start-0.5,e_end+0.5,100,-5.0,0.0);
+ TH1D *h4 = new TH1D("h4","Reconstructed energy per event", 100,e_start-0.5,e_end+0.5);
+ TH1D *h5 = new TH1D("h5","Number of Clusters per event", 5,-0.5,4.5);
+ TH1D *h6 = new TH1D("h6","Scattering Angle(#theta) with CUT", 100,130.0,180.0);
+ TH1D *h7 = new TH1D("h7","Pseudo-rapidity(#eta) with CUT", 100,-5.0,0.0);
+ TH2D *h8 = new TH2D("h8","Cluster Hit Position", 62,-62.0,62.0,62,-62.0,62.0);
+ TH2D *h9 = new TH2D("h9","All Hit Position", 62,-62.0,62.0,62,-62.0,62.0);
+ TH1D *h10 = new TH1D("hEnergyRes","Energy Resolution", 100,-0.3,0.3);
+ TH1D *h11 = new TH1D("hEnergyResCUT","Energy Resolution with CUT", 100,-0.3,0.3);
+ TH1D *h12 = new TH1D("h12","Thrown momentum", 61,e_start-0.5,e_end+0.5);
+ TH1D *h13 = new TH1D("h13","Thrown momentum for reconstructed particle", 61,e_start-0.5,e_end+0.5);
+ TH1D *h14 = new TH1D("h14","Ratio p_{rec}/p_{thr}", 61,e_start-0.5,e_end+0.5);
+
+ // Declare ellipse for boundary of crystal calorimeter
+ TEllipse *ell1 = new TEllipse(0.0, 0.0, 60.0, 60.0);
+ ell1->SetFillStyle(4000);
+ TEllipse *ell2 = new TEllipse(0.0, 0.0, 12.0, 12.0);
+ ell2->SetFillStyle(4000);
+
+ // Total number of entries
+ Int_t nentries = t->GetEntries();
+
+ // Variables are used in calculation
+ Double_t r; // Radius [cm]
+ Double_t phi; // Azimuth [degree]
+ Double_t theta; // Inclination [degree]
+ Double_t eta; // Pseudo-rapidity [unitless]
+ Float_t cluster_e; // Cluster energy [GeV]
+ Float_t total_cluster_e; // Add up clusters per event [GeV]
+ Double_t total_thr_e; // Thrown energy [GeV]
+ Double_t total_truth_sim_e; // Add up truth simulated deposit energy [GeV]
+ Double_t total_sim_e; // Add up truth simulated deposit energy [GeV]
+ Double_t momentum2; // momentum squared
+ Double_t momentum; // momentum
+ Double_t mass2; // mass squared
+ Double_t diff_energy; // delta(E) = reconstructed energy - thrown energy
+ Double_t eng_res; // delta(E)/thrown energy
+
+ // Loop over event by event
+ for (int ievent = 0; ievent < nentries; ievent++)
+ {
+ // Read event by event
+ t->GetEntry(ievent);
+ // Read number of hits/clusters
+/*
+ Int_t nmcparticle = 2;
+ Int_t nCrystalEcalHits = CrystalEcalHits2_;
+*/
+ Int_t nreconhits = RecoEcalHits_;
+ Int_t ncluster = EcalClusters_;
+ // Initialize total energy variables
+ total_thr_e = 0.0;
+ total_truth_sim_e = 0.0;
+ total_sim_e = 0.0;
+ total_cluster_e = 0.0;
+ // Thrown energy, momentum, and mass
+/*
+ momentum2 = px[nmcparticle]*px[nmcparticle]+py[nmcparticle]*py[nmcparticle]+pz[nmcparticle]*pz[nmcparticle];
+ momentum = TMath::Sqrt(momentum2);
+ mass2 = mass[nmcparticle]*mass[nmcparticle];
+ total_thr_e = sqrt(momentum2 + mass2)/1.e+3;
+ h12->Fill(momentum, 1.0);
+
+ // Loop over simulated hit by simulated hit
+ for (int isimhit=0; isimhit < nCrystalEcalHits; isimhit++)
+ {
+ total_truth_sim_e += truth_deposit[isimhit];
+ total_sim_e += energyDeposit[isimhit]/1.e+3;
+ }
+*/
+ // Loop over reconstructed hit by reconstructed hit
+ for(int ireconhit=0; ireconhit < nreconhits; ireconhit++)
+ h9->Fill(rec_x_pos[ireconhit],rec_y_pos[ireconhit], 1.0);
+
+ h5->Fill(ncluster, 1.0);
+
+ // Loop over cluster by cluster
+ for (int icluster=0; icluster < ncluster; icluster++)
+ {
+ r = TMath::Sqrt((cluster_x_pos[icluster]*cluster_x_pos[icluster]) +
+ (cluster_y_pos[icluster]*cluster_y_pos[icluster]) +
+ (cluster_z_pos[icluster]*cluster_z_pos[icluster]));
+ phi = TMath::ATan(cluster_y_pos[icluster]/cluster_x_pos[icluster]) * TMath::RadToDeg();
+ theta = TMath::ACos(cluster_z_pos[icluster] / r) * TMath::RadToDeg();
+ eta = -1.0 * TMath::Log(TMath::Tan((theta*TMath::DegToRad())/2.0));
+ cluster_e = cluster_energy[icluster] / 1.e+3;
+ total_cluster_e += cluster_e;
+ }
+
+ // Select events with one cluster
+ if(ncluster == 1)
+ {
+ diff_energy = total_cluster_e - total_thr_e;
+ eng_res = (diff_energy / total_thr_e);
+ h10->Fill(eng_res, 1.0);
+
+ h1->Fill(theta, 1.0);
+ h2->Fill(eta, 1.0);
+ h3->Fill(cluster_e, eta, 1.0);
+ h4->Fill(total_cluster_e, 1.0);
+ h8->Fill(cluster_x_pos[0],cluster_y_pos[0], 1.0);
+ h10->Fill(eng_res, 1.0);
+
+ if(total_cluster_e > 0.2*total_thr_e)
+ {
+ h6->Fill(theta, 1.0);
+ h7->Fill(eta, 1.0);
+ h11->Fill(eng_res, 1.0);
+ }
+/*
+ if(total_cluster_e > 0.9*total_thr_e)
+ h13->Fill(momentum, 1.0);
+*/
+ }
+ }
+ // Drawing and Saving figures
+ c1->cd();
+ h1->SetLineColor(kBlue);
+ h1->SetLineWidth(2);
+ h1->GetYaxis()->SetRangeUser(0.0,h1->GetMaximum()+50.0);
+ h1->GetXaxis()->SetTitle("#theta [degree]");
+ h1->GetYaxis()->SetTitle("events");
+ h1->GetYaxis()->SetTitleOffset(1.4);
+ h1->DrawClone();
+ c1->SaveAs("results/electron_theta_hist_0GeVto30GeV.png");
+ c1->SaveAs("results/electron_theta_hist_0GeVto30GeV.pdf");
+
+ c2->cd();
+ h2->SetLineColor(kBlue);
+ h2->SetLineWidth(2);
+ h2->GetYaxis()->SetRangeUser(0.0,h2->GetMaximum()+50.0);
+ h2->GetXaxis()->SetTitle("#eta");
+ h2->GetYaxis()->SetTitle("events");
+ h2->GetYaxis()->SetTitleOffset(1.4);
+ h2->DrawClone();
+ c2->SaveAs("results/electron_eta_hist_0GeVto30GeV.png");
+ c2->SaveAs("results/electron_eta_hist_0GeVto30GeV.pdf");
+
+ c3->cd();
+ h3->GetXaxis()->SetTitle("Cluster energy [GeV]");
+ h3->GetYaxis()->SetTitle("#eta");
+ h3->GetYaxis()->SetTitleOffset(1.4);
+ h3->DrawClone("COLZ");
+ c3->SaveAs("results/eletron_E_vs_eta_hist_0GeVto30GeV.png");
+ c3->SaveAs("results/eletron_E_vs_eta_hist_0GeVto30GeV.pdf");
+
+ c4->cd();
+ c4->SetLogy(1);
+ h4->SetLineColor(kBlue);
+ h4->SetLineWidth(2);
+ h4->GetXaxis()->SetTitle("reconstructed energy [GeV]");
+ h4->GetYaxis()->SetTitle("events");
+ h4->GetYaxis()->SetTitleOffset(1.4);
+ h4->DrawClone();
+ c4->SaveAs("results/electron_Erec_hist_0GeVto30GeV.png");
+ c4->SaveAs("results/electron_Erec_hist_0GeVto30GeV.pdf");
+
+ c5->cd();
+ c5->SetLogy(1);
+ h5->SetLineColor(kBlue);
+ h5->SetLineWidth(2);
+ h5->GetXaxis()->SetTitle("Number of Clusters");
+ h5->GetYaxis()->SetTitle("events");
+ h5->GetYaxis()->SetTitleOffset(1.4);
+ h5->DrawClone();
+ c5->SaveAs("results/electron_ncluster_hist_0GeVto30GeV.png");
+ c5->SaveAs("results/electron_ncluster_hist_0GeVto30GeV.pdf");
+
+ c6->cd();
+ h6->SetLineColor(kBlue);
+ h6->SetLineWidth(2);
+ h6->GetYaxis()->SetRangeUser(0.0,h1->GetMaximum()+50.0);
+ h6->GetXaxis()->SetTitle("#theta [degree]");
+ h6->GetYaxis()->SetTitle("events");
+ h6->GetYaxis()->SetTitleOffset(1.4);
+ h6->DrawClone();
+ c6->SaveAs("results/electron_theta_hist_CUT_0GeVto30GeV.png");
+ c6->SaveAs("results/electron_theta_hist_CUT_0GeVto30GeV.pdf");
+
+ c7->cd();
+ h7->SetLineColor(kBlue);
+ h7->SetLineWidth(2);
+ h7->GetYaxis()->SetRangeUser(0.0,h2->GetMaximum()+50.0);
+ h7->GetXaxis()->SetTitle("#eta");
+ h7->GetYaxis()->SetTitle("events");
+ h7->GetYaxis()->SetTitleOffset(1.4);
+ h7->DrawClone();
+ c7->SaveAs("results/electron_eta_hist_CUT_0GeVto30GeV.png");
+ c7->SaveAs("results/electron_eta_hist_CUT_0GeVto30GeV.pdf");
+
+ c8->cd();
+ h8->GetXaxis()->SetTitle("Hit position X [cm]");
+ h8->GetYaxis()->SetTitle("Hit position Y [cm]");
+ h8->GetYaxis()->SetTitleOffset(1.4);
+ h8->DrawClone("COLZ");
+ ell1->Draw("same");
+ ell2->Draw("same");
+ c8->SaveAs("results/electron_hit_pos_cluster_0GeVto30GeV.png");
+ c8->SaveAs("results/electron_hit_pos_cluster_0GeVto30GeV.pdf");
+
+ c9->cd();
+ h9->GetXaxis()->SetTitle("Hit position X [cm]");
+ h9->GetYaxis()->SetTitle("Hit position Y [cm]");
+ h9->GetYaxis()->SetTitleOffset(1.4);
+ h9->DrawClone("COLZ");
+ ell1->Draw("same");
+ ell2->Draw("same");
+ c9->SaveAs("results/electron_hit_pos_all_0GeVto30GeV.png");
+ c9->SaveAs("results/electron_hit_pos_all_0GeVto30GeV.pdf");
+
+ c10->cd();
+ h10->SetLineColor(kBlue);
+ h10->SetLineWidth(2);
+ h10->GetXaxis()->SetTitle("#DeltaE/E");
+ h10->GetYaxis()->SetTitle("events");
+ h10->GetYaxis()->SetTitleOffset(1.4);
+ h10->Fit("gaus");
+ h10->DrawClone();
+ c10->SaveAs("results/hEngRes.png");
+ c10->SaveAs("results/hEngRes.pdf");
+
+ c11->cd();
+ h11->SetLineColor(kBlue);
+ h11->SetLineWidth(2);
+ h11->GetXaxis()->SetTitle("#DeltaE/E");
+ h11->GetYaxis()->SetTitle("events");
+ h11->GetYaxis()->SetTitleOffset(1.4);
+ h11->Fit("gaus");
+ h11->DrawClone();
+ c11->SaveAs("results/hEngRes_CUT.png");
+ c11->SaveAs("results/hEngRes_CUT.pdf");
+
+ c12->cd();
+ h12->SetLineColor(kBlue);
+ h12->SetLineWidth(2);
+ h12->GetXaxis()->SetTitle("Reconstructed Momentum [GeV]");
+ h12->GetYaxis()->SetTitle("events");
+ h12->GetYaxis()->SetTitleOffset(1.4);
+ h12->DrawClone();
+ c12->SaveAs("results/hThrMom.png");
+ c12->SaveAs("results/hThrMom.pdf");
+
+ c13->cd();
+ h13->SetLineColor(kBlue);
+ h13->SetLineWidth(2);
+ h13->GetXaxis()->SetTitle("Thrown Momentum [GeV]");
+ h13->GetYaxis()->SetTitle("events");
+ h13->GetYaxis()->SetTitleOffset(1.4);
+ h13->DrawClone();
+ c13->SaveAs("results/hThrMom_Rec.png");
+ c13->SaveAs("results/hThrMom_Rec.pdf");
+
+ c14->cd();
+ h14 = (TH1D*)h12->Clone();
+ h14->Divide(h13);
+ h14->SetTitle("Ratio p_{rec}/p_{thr}");
+ h14->SetLineColor(kBlue);
+ h14->SetLineWidth(2);
+ h14->GetXaxis()->SetTitle("p [GeV]");
+ h14->GetYaxis()->SetTitle("p_{rec}/p_{thr}");
+ h14->GetYaxis()->SetTitleOffset(1.4);
+ h14->DrawClone();
+ c14->SaveAs("results/hAcceptance.png");
+ c14->SaveAs("results/hAcceptance.pdf");
+
+ return 0;
+}
--
GitLab