Add spatial plots

benchmarks/imaging_ecal/scripts/emcal_barrel_energy_scan_analysis.cxx

@@ -15,8 +15,12 @@
 #include "eicd/ClusterData.h"
 #include "eicd/RawCalorimeterHitCollection.h"
 #include "eicd/RawCalorimeterHitData.h"
+#include "eicd/ClusterLayerData.h"
+#include "eicd/Cluster3DInfoData.h"
+#include "eicd/Cluster2DInfoData.h"
 
 #include "TCanvas.h"
+#include "TVector3.h"
 #include "TStyle.h"
 #include "TMath.h"
 #include "TH1.h"
@@ -73,6 +77,27 @@ std::tuple <double, double, double, double, double, double, double, double> extr
     return energy;
   };
 
+  // Phi
+  auto Phithr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+    auto p = input[2];
+    TVector3 mom(p.ps.x, p.ps.y, p.ps.z);
+    return mom.Phi();
+  };
+
+  // Theta
+  auto Thetathr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+    auto p = input[2];
+    TVector3 mom(p.ps.x, p.ps.y, p.ps.z);
+    return mom.Theta();
+  };
+
+  // Eta
+  auto Etathr = [](std::vector<dd4pod::Geant4ParticleData> const& input) {
+    auto p = input[2];
+    TVector3 mom(p.ps.x, p.ps.y, p.ps.z);
+    return mom.Eta();
+  };
+
   // Reconstructed Energy [GeV]
   auto Erec = [](const std::vector<eic::CalorimeterHitData>& evt) {
     auto total_edep = 0.0;
@@ -100,25 +125,82 @@ std::tuple <double, double, double, double, double, double, double, double> extr
     return total_edep;
   };
 
+  // Avarage theta
+  auto EClusInfoImgTheta = [](const std::vector<eic::Cluster3DInfoData>& evt) {
+    double THETA = 0;
+    int count = 0;
+    for (const auto& i: evt){
+     THETA += i.polar.theta;
+     count++;
+    }
+    return THETA / count;
+  };
+
+  // Avarage Phi
+  auto EClusInfoImgPhi = [](const std::vector<eic::Cluster3DInfoData>& evt) {
+    double PHI = 0;
+    int count = 0;
+    for (const auto& i: evt){
+     PHI += i.polar.phi;
+     count++;
+    }
+    return PHI / count;
+  };
+
+  // Avarage Eta
+  auto EClusInfoImgEta = [](const std::vector<eic::Cluster3DInfoData>& evt) {
+    double ETA = 0;
+    int count = 0;
+    for (const auto& i: evt){
+     ETA += i.eta;
+     count++;
+    }
+    return ETA / count;
+  };
+
   // Sampling fraction = Esampling / Ethrown
   auto fsam = [](const double sampled, const double thrown) {
     return sampled / thrown;
   };
 
+  // Delta Theta
+  auto deltaTheta = [](const double Thetathr, const double Thetarec) {
+    return (Thetathr-Thetarec)/Thetathr;
+  };
+
+  // Delta Phi
+  auto deltaTheta = [](const double Phithr, const double Phirec) {
+    return (Phithr-Phirec)/Phithr;
+  };
+
+  // Delta Eta
+  auto deltaTheta = [](const double Etathr, const double Etarec) {
+    return (Etathr-Etarec)/Etathr;
+  };
+
   // Define variables
   auto d1 = d0.Define("Ethr", Ethr, {"mcparticles2"})
-                .Define("ErecImg", Erec, {"RecoEcalBarrelImagingHits"})
-                .Define("ErecScFi", Erec, {"EcalBarrelScFiHitsReco"})
-                .Define("EdigiImg", Edigi, {"DigiEcalBarrelImagingHits"})
-                .Define("EdigiScFi", Edigi, {"EcalBarrelScFiHitsDigi"})
-                .Define("NClusterImg", ncluster, {"EcalBarrelImagingClusters"})
-                .Define("NClusterScFi", ncluster, {"EcalBarrelScFiClusters"})
-                .Define("EClusterImg", Ecluster, {"EcalBarrelImagingClusters"})
-                .Define("EClusterScFi", Ecluster, {"EcalBarrelScFiClusters"})
-                .Define("fsamRecImg", fsam, {"ErecImg", "Ethr"})
-                .Define("fsamRecScFi", fsam, {"ErecScFi", "Ethr"})
-                .Define("fsamClusterImg", fsam, {"EClusterImg", "Ethr"})
-                .Define("fsamClusterScFi", fsam, {"EClusterScFi", "Ethr"});
+              .Define("Thetathr", Thetathr, {"mcparticles2"})
+              .Define("Phithr", Phithr, {"mcparticles2"})
+              .Define("Etathr", Etathr, {"mcparticles2"})
+              .Define("ErecImg", Erec, {"RecoEcalBarrelImagingHits"})
+              .Define("ErecScFi", Erec, {"EcalBarrelScFiHitsReco"})
+              .Define("EdigiImg", Edigi, {"DigiEcalBarrelImagingHits"})
+              .Define("EdigiScFi", Edigi, {"EcalBarrelScFiHitsDigi"})
+              .Define("NClusterImg", ncluster, {"EcalBarrelImagingClusters"})
+              .Define("NClusterScFi", ncluster, {"EcalBarrelScFiClusters"})
+              .Define("EClusterImg", Ecluster, {"EcalBarrelImagingClusters"})
+              .Define("EClusterScFi", Ecluster, {"EcalBarrelScFiClusters"})
+              .Define("fsamRecImg", fsam, {"ErecImg", "Ethr"})
+              .Define("fsamRecScFi", fsam, {"ErecScFi", "Ethr"})
+              .Define("fsamClusterImg", fsam, {"EClusterImg", "Ethr"})
+              .Define("fsamClusterScFi", fsam, {"EClusterScFi", "Ethr"})
+              .Define("EClusInfoImgTheta", EClusInfoImgTheta, {"EcalBarrelImagingClustersInfo"})
+              .Define("EClusInfoImgPhi", EClusInfoImgPhi, {"EcalBarrelImagingClustersInfo"})
+              .Define("EClusInfoImgEta", EClusInfoImgEta, {"EcalBarrelImagingClustersInfo"})
+              .Define("DeltaTheta",deltaTheta,{"Thetathr","EClusInfoImgTheta"})
+              .Define("DeltaPhi",deltaTheta,{"Phithr","EClusInfoImgPhi"})
+              .Define("DeltaEta",deltaTheta,{"Etathr","EClusInfoImgEta"});
 
   // Define Histograms
   auto hEthr = d1.Histo1D({"hEthr", "Thrown Energy; Thrown Energy [GeV]; Events", 100, 0.0, 25.0},"Ethr");
@@ -136,8 +218,14 @@ std::tuple <double, double, double, double, double, double, double, double> extr
   auto hfsamScFi = d1.Histo1D({"hfsamScFi", "Cluster Energy/E true; Cluster Energy/E true; Events", 100, 0.8, 1.2},"fsamClusterScFi");
   auto hfsamRecScFi = d1.Histo1D({"hfsamRecScFi", "Reco Hits Energy/E true; Reco Hits Energy/E true; Events", 50, 0.0, 0.25},"fsamRecScFi");
 
+  auto d2 = d0.Filter("NClusterImg==1");
+  auto hDeltaPhi = d2.Histo1D({"hDeltaPhi", "Phi resolution; (Phi true - Phi reco)/Phi true; Events", 700, -3.5, 3.5},"DeltaPhi");
+  auto hDeltaTheta = d2.Histo1D({"hDeltaTheta", "Theta resolution; (Theta true - Theta reco)/Theta true; Events", 700, -3.5, 3.5},"DeltaTheta");
+  auto hDeltaEta = d2.Histo1D({"hDeltaEta", "Eta resolution; (Eta true - Eta reco)/Eta true; Events", 700, -3.5, 3.5},"DeltaEta");
+  
   // Event Counts
   auto nevents_thrown = d1.Count();
+
   std::cout << "Number of Thrown Events: " << (*nevents_thrown) << "\n";
 
   // Draw Histograms
@@ -187,6 +275,34 @@ std::tuple <double, double, double, double, double, double, double, double> extr
   save_canvas(c2, "E_digi_rec", E_label, particle_label);
 
 
+  TCanvas* c2R = new TCanvas("c2R", "c2R", 1400, 1000);
+  c2R->Divide(2,2);
+  c2R->cd(1);
+  hDeltaPhi->GetYaxis()->SetTitleOffset(1.4);
+  hDeltaPhi->SetLineWidth(2);
+  hDeltaPhi->SetLineColor(kBlue);
+  hDeltaPhi->DrawClone();
+  auto up_fit = hDeltaPhi->GetMean() + 5*hDeltaPhi->GetStdDev();
+  auto down_fit = hDeltaPhi->GetMean() - 5*hDeltaPhi->GetStdDev();
+  hDeltaPhi->GetXaxis()->SetRangeUser(down_fit,up_fit);
+  c2R->cd(2);
+  hDeltaTheta->GetYaxis()->SetTitleOffset(1.4);
+  hDeltaTheta->SetLineWidth(2);
+  hDeltaTheta->SetLineColor(kBlue);
+  up_fit = hDeltaTheta->GetMean() + 5*hDeltaTheta->GetStdDev();
+  down_fit = hDeltaTheta->GetMean() - 5*hDeltaTheta->GetStdDev();
+  hDeltaTheta->GetXaxis()->SetRangeUser(down_fit,up_fit);
+  hDeltaTheta->DrawClone();
+  c2R->cd(3);
+  hDeltaEta->GetYaxis()->SetTitleOffset(1.4);
+  hDeltaEta->SetLineWidth(2);
+  hDeltaEta->SetLineColor(kBlue);
+  hDeltaEta->DrawClone();
+  up_fit = hDeltaEta->GetMean() + 5*hDeltaEta->GetStdDev();
+  down_fit = hDeltaEta->GetMean() - 5*hDeltaEta->GetStdDev();
+  hDeltaEta->GetXaxis()->SetRangeUser(down_fit,up_fit);
+  save_canvas(c2R, "spatial_res", E_label, particle_label);
+
 {
   TCanvas* c3 = new TCanvas("c3", "c3", 1400, 500);
   c3->Divide(2,1);

benchmarks/imaging_ecal/scripts/emcal_barrel_energy_scan_analysis_pi_e.cxx

@@ -204,7 +204,7 @@ std::tuple <double, double> extract_sampling_fraction_parameters(std::string E_l
       }
       std::cout << "efficiency cut bin: " << efficiency_cut_bin << " efficiency: " << 1. - (helectron->Integral(1,efficiency_cut_bin))/(*ngen_helectron) << std::endl;
       auto ngen_hpion = dpions.Count();
-      auto pion_suppresion = hpion->Integral(efficiency_cut_bin, nbins)/(*ngen_hpion);
+      auto pion_suppresion = (*ngen_hpion)/hpion->Integral(efficiency_cut_bin, nbins);
       std::cout << "pion supression: " << pion_suppresion << std::endl;
     }
 
@@ -283,17 +283,18 @@ std::tuple <double, double> extract_sampling_fraction_parameters(std::string E_l
       helectron->SetLineWidth(2);
       helectron->SetLineColor(kBlue);
       helectron->DrawClone();
+      auto up_range = helectron->GetMean() + 3*helectron->GetStdDev();
+      auto down_range = helectron->GetMean() - 3*helectron->GetStdDev();
+      // h->SetLineWidth(2);
+      // h->SetLineColor(kBlue);
+    
+      helectron->GetXaxis()->SetRange(helectron->GetXaxis()->GetBinUpEdge(1), up_range); // skip 0th bin 
       hpion->SetLineColor(kRed);
       hpion->DrawClone("same");
 
       // auto h = hEsim_layer->DrawCopy();
 
-      // auto up_range = h->GetMean() + 3*h->GetStdDev();
-      // auto down_range = h->GetMean() - 3*h->GetStdDev();
-      // h->SetLineWidth(2);
-      // h->SetLineColor(kBlue);
-    
-      // h->GetXaxis()->SetRange(h->GetXaxis()->GetBinUpEdge(1), up_range); // skip 0th bin 
+
       // auto mean_layer = h->GetMean();
       // auto rms_layer = h->GetStdDev();
       // h->GetXaxis()->SetRange(); // reset the range