Resolve "Add layer energy-deposit profile to energy resolution benchmark"

Closes #26 (closed)

benchmarks/barrel_ecal/scripts/emcal_barrel_energy_scan_analysis.cxx

@@ -18,6 +18,10 @@
 #include "TH1D.h"
 #include "TGraphErrors.h"
 
+#include "DD4hep/Detector.h"
+#include "DDG4/Geant4Data.h"
+#include "DDRec/CellIDPositionConverter.h"
+
 using ROOT::RDataFrame;
 using namespace ROOT::VecOps;
 
@@ -44,7 +48,21 @@ void save_canvas(TCanvas* c, std::string var_label, std::string E_label, std::st
   save_canvas(c, label_with_E);
 }
 
-std::tuple <double, double, double, double> extract_sampling_fraction_parameters(std::string particle_label, std::string E_label)
+dd4hep::Detector& load_detector(std::string detector_path = "", std::string detector_name = "topside")
+{
+  if(std::getenv("DETECTOR_PATH")) {
+    detector_path = std::getenv("DETECTOR_PATH");
+  }
+  if(std::getenv("JUGGLER_DETECTOR")) {
+    detector_name = std::getenv("JUGGLER_DETECTOR");
+  }
+
+  dd4hep::Detector& detector = dd4hep::Detector::getInstance();
+  detector.fromCompact(fmt::format("{}/{}.xml", detector_path,detector_name));
+  return detector;
+}
+
+std::tuple <double, double, double, double> extract_sampling_fraction_parameters(std::string particle_label, std::string E_label, dd4hep::Detector& detector)
 {
   std::string input_fname = fmt::format("sim_output/energy_scan/{}/sim_emcal_barrel_{}.root", E_label, particle_label);
   ROOT::EnableImplicitMT();
@@ -73,6 +91,25 @@ std::tuple <double, double, double, double> extract_sampling_fraction_parameters
     return sampled / thrown;
   };
 
+  // Energy deposited in layers 
+
+  //dd4hep::rec::CellIDPositionConverter cellid_converter(detector);
+
+  auto decoder = detector.readout("EcalBarrelHits").idSpec().decoder();
+  fmt::print("{}\n", decoder->fieldDescription());
+  auto layer_index = decoder->index("layer");
+  fmt::print(" layer index is {}.\n", layer_index);
+
+  auto Layer = [=](const std::vector<dd4pod::CalorimeterHitData>& evt) {
+    std::vector<double> vec_layer_ids;
+    for (const auto& i: evt) {
+      auto layer_id = decoder->get(i.cellID, layer_index);
+      vec_layer_ids.push_back(layer_id);
+    }
+    RVec<double> layer_ids(vec_layer_ids.data(), vec_layer_ids.size());
+    return layer_ids;
+   };
+
   // Define variables
   auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
                 .Define("nhits", nhits, {"EcalBarrelHits"})
@@ -91,8 +128,35 @@ std::tuple <double, double, double, double> extract_sampling_fraction_parameters
       {"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 500, 0.0, 0.5},
       "Esim");
   auto hfsam = d1.Histo1D(
-      {"hfsam", "Sampling Fraction; Sampling Fraction; Events", 200, 0.0, 0.05},
+      {"hfsam", "Sampling Fraction; Sampling Fraction; Events", 800, 0.0, 0.2},
       "fsam");
+      
+  // Draw energy per layer
+  for(int layer=1; layer<21; layer++) {
+    auto Esim_layer = [=](const std::vector<dd4pod::CalorimeterHitData>& evt) {
+      auto layer_edep = 0.0;
+      for (const auto& i: evt) {
+        if (decoder->get(i.cellID, layer_index) == layer) {
+          layer_edep += i.energyDeposit;
+        }
+      }
+      //std::cout << "layer: " << layer << " layer_edep: " << layer_edep << std::endl;
+      return layer_edep;
+    };
+
+    auto d2 = d0.Define(fmt::format("Esim_layer_{}",layer).c_str(), Esim_layer, {"EcalBarrelHits"});
+    std::cout << "Layer to process: " << layer << std::endl;
+
+    auto hEsim_layer = d2.Histo1D({fmt::format("hEsim_layer_{}",layer).c_str(),fmt::format("Energy Deposit in layer {}; Energy Deposit [Gev]; Events",layer).c_str(), 500, 0.0, 0.01}, fmt::format("Esim_layer_{}",layer));
+
+    TCanvas* clayer = new TCanvas(fmt::format("c1_layer{}",layer).c_str(), fmt::format("c1_layer{}",layer).c_str(), 700, 500);
+    clayer->cd();
+    auto h = hEsim_layer->DrawCopy();
+    double up_fit = h->GetMean() + 5*h->GetStdDev();
+    double down_fit = h->GetMean() - 5*h->GetStdDev();
+    h->GetXaxis()->SetRangeUser(0.,up_fit);
+    save_canvas(clayer, fmt::format("Esim_layer{}",layer), E_label, particle_label);
+  }
 
   // Event Counts
   auto nevents_thrown = d1.Count();
@@ -191,12 +255,13 @@ void emcal_barrel_energy_scan_analysis(std::string particle_label = "electron")
   gStyle->SetPadRightMargin(0.14);
 
     auto scanned_energies = read_scanned_energies(fmt::format("sim_output/emcal_barrel_energy_scan_points_{}.txt", particle_label));
+    auto detector = load_detector();
     
     TGraphErrors gr_fsam(scanned_energies.size()-1);
     TGraphErrors gr_fsam_res(scanned_energies.size()-1);
 
     for (const auto& E_label : scanned_energies) {
-        auto [fsam, fsam_res, fsam_err, fsam_res_err] = extract_sampling_fraction_parameters(particle_label, E_label);
+        auto [fsam, fsam_res, fsam_err, fsam_res_err] = extract_sampling_fraction_parameters(particle_label, E_label, detector);
         auto E = std::stod(E_label);
 
         gr_fsam.SetPoint(gr_fsam.GetN(),E,100*fsam);