Skip to content
Snippets Groups Projects

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

Compare and
4 files
+ 68
5
Compare changes
  • Side-by-side
  • Inline
Files
4
@@ -18,6 +18,10 @@
@@ -18,6 +18,10 @@
#include "TH1D.h"
#include "TH1D.h"
#include "TGraphErrors.h"
#include "TGraphErrors.h"
 
#include "DD4hep/Detector.h"
 
#include "DDG4/Geant4Data.h"
 
#include "DDRec/CellIDPositionConverter.h"
 
using ROOT::RDataFrame;
using ROOT::RDataFrame;
using namespace ROOT::VecOps;
using namespace ROOT::VecOps;
@@ -73,6 +77,35 @@ std::tuple <double, double, double, double> extract_sampling_fraction_parameters
@@ -73,6 +77,35 @@ std::tuple <double, double, double, double> extract_sampling_fraction_parameters
return sampled / thrown;
return sampled / thrown;
};
};
 
// Energy deposited in layers
 
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));
 
//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
// Define variables
auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
auto d1 = d0.Define("Ethr", Ethr, {"mcparticles"})
.Define("nhits", nhits, {"EcalBarrelHits"})
.Define("nhits", nhits, {"EcalBarrelHits"})
@@ -91,8 +124,35 @@ std::tuple <double, double, double, double> extract_sampling_fraction_parameters
@@ -91,8 +124,35 @@ std::tuple <double, double, double, double> extract_sampling_fraction_parameters
{"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 500, 0.0, 0.5},
{"hEsim", "Energy Deposit; Energy Deposit [GeV]; Events", 500, 0.0, 0.5},
"Esim");
"Esim");
auto hfsam = d1.Histo1D(
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");
"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
// Event Counts
auto nevents_thrown = d1.Count();
auto nevents_thrown = d1.Count();
Loading