From a6a4dab991af7f9decc7247dc6d5236983260f33 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Einar=20El=C3=A9n?= Date: Sun, 7 Apr 2024 14:51:44 +0200 Subject: [PATCH 1/3] Quick and dirty implementation --- include/Hcal/HcalSimpleDigiAndRecProducer.h | 40 ++++++ python/digi.py | 30 ++++ src/Hcal/HcalSimpleDigiAndRecProducer.cxx | 150 ++++++++++++++++++++ 3 files changed, 220 insertions(+) create mode 100644 include/Hcal/HcalSimpleDigiAndRecProducer.h create mode 100644 src/Hcal/HcalSimpleDigiAndRecProducer.cxx diff --git a/include/Hcal/HcalSimpleDigiAndRecProducer.h b/include/Hcal/HcalSimpleDigiAndRecProducer.h new file mode 100644 index 0000000..79da032 --- /dev/null +++ b/include/Hcal/HcalSimpleDigiAndRecProducer.h @@ -0,0 +1,40 @@ +#ifndef HCALSIMPLEDIGIANDRECPRODUCER_H +#define HCALSIMPLEDIGIANDRECPRODUCER_H +#include "DetDescr/DetectorID.h" +#include "DetDescr/HcalGeometry.h" +#include "DetDescr/HcalID.h" +#include "Framework/Configure/Parameters.h" +#include "Framework/EventDef.h" +#include "Framework/EventProcessor.h" +#include "Framework/RandomNumberSeedService.h" +#include "TRandom3.h" +#include "Tools/NoiseGenerator.h" + +namespace hcal { +class HcalSimpleDigiAndRecProducer : public framework::Producer { + public: + HcalSimpleDigiAndRecProducer(const std::string& name, + framework::Process& process) + : framework::Producer{name, process} {} + ~HcalSimpleDigiAndRecProducer() override = default; + void configure(framework::config::Parameters& ps) override; + void SetupRandomNumberGeneration(); + void produce(framework::Event& event) override; + + private: + std::string input_coll_name{}; + std::string input_pass_name{}; + std::string output_coll_name{}; + double mev_per_mip{}; + double pe_per_mip{}; + double attenuation_length{}; + double mean_noise{}; + double position_resolution{}; + std::unique_ptr random{nullptr}; + std::unique_ptr noiseGenerator{nullptr}; + int readout_threshold{2}; +}; + +} // namespace hcal + +#endif /* HCALSIMPLEDIGIANDRECPRODUCER_H */ diff --git a/python/digi.py b/python/digi.py index 46a1d7b..aa8ed27 100644 --- a/python/digi.py +++ b/python/digi.py @@ -215,3 +215,33 @@ def __init__(self, instance_name = 'hcalDoubleRecon', pass_name = '', coll_name self.pass_name = pass_name self.rec_coll_name = rec_coll_name self.rec_pass_name = rec_pass_name + +class HcalSimpleDigiAndRecProducer(Producer) : + """Configuration for Digitization producer in the HCal + Sets all parameters to reasonable defaults. + Examples + -------- + from LDMX.EventProc.hcal import HcalDigiProducer + p.sequence.append( HcalDigiProducer() ) + """ + + def __init__(self,name = 'hcalSimpleDigiAndRec') : + super().__init__(name,'hcal::HcalSimpleDigiAndRecProducer','Hcal') + self.input_coll_name = 'HcalSimHits' + self.input_pass_name = '' + self.output_coll_name = 'HcalRecHits' + + self.mean_noise = 0.02 + self.readout_threshold= 1 + self.strips_side_lr_per_layer = 12 + self.num_side_lr_hcal_layers = 26 + self.strips_side_tb_per_layer = 12 + self.num_side_tb_hcal_layers = 28 + self.strips_back_per_layer = 60 # n strips correspond to 5 cm wide bars + self.num_back_hcal_layers = 96 + self.super_strip_size = 1 # 1 = 5 cm readout, 2 = 10 cm readout, ... + self.mev_per_mip = 4.66 # measured 1.4 MeV for a 6mm thick tile, so for 20mm bar = 1.4*20/6 + self.pe_per_mip = 68. # PEs per MIP at 1m (assume 80% attentuation of 1m) + self.attenuation_length = 5. # this is in m + self.position_resolution = 150. # this is in mm + self.sim_hit_pass_name = '' #use any pass available diff --git a/src/Hcal/HcalSimpleDigiAndRecProducer.cxx b/src/Hcal/HcalSimpleDigiAndRecProducer.cxx new file mode 100644 index 0000000..6699afd --- /dev/null +++ b/src/Hcal/HcalSimpleDigiAndRecProducer.cxx @@ -0,0 +1,150 @@ +#include "Hcal/HcalSimpleDigiAndRecProducer.h" + +namespace hcal { + +void HcalSimpleDigiAndRecProducer::configure( + framework::config::Parameters& ps) { + input_coll_name = ps.getParameter("input_coll_name"); + input_pass_name = ps.getParameter("input_pass_name"); + output_coll_name = ps.getParameter("output_coll_name"); + mev_per_mip = ps.getParameter("mev_per_mip"); + pe_per_mip = ps.getParameter("pe_per_mip"); + attenuation_length = ps.getParameter("attenuation_length"); + readout_threshold = ps.getParameter("readout_threshold"); + mean_noise = ps.getParameter("mean_noise"); + position_resolution = ps.getParameter("position_resolution"); +} + +void HcalSimpleDigiAndRecProducer::SetupRandomNumberGeneration() { + if (noiseGenerator == nullptr) { + noiseGenerator = std::make_unique(mean_noise, false); + noiseGenerator->setNoiseThreshold( + 1); // hard-code this number, create noise hits for non-zero PEs! + } + if (!noiseGenerator->hasSeed()) { + const framework::RandomNumberSeedService& rseed = + getCondition( + framework::RandomNumberSeedService::CONDITIONS_OBJECT_NAME); + noiseGenerator->seedGenerator( + rseed.getSeed("HcalSimpleDigiAndRecProducer::NoiseGenerator")); + } + if (random == nullptr) { + const framework::RandomNumberSeedService& rseed = + getCondition( + framework::RandomNumberSeedService::CONDITIONS_OBJECT_NAME); + random = std::make_unique( + rseed.getSeed("HcalSimpleDigiAndRecProducer")); + } +} + +void HcalSimpleDigiAndRecProducer::produce(framework::Event& event) { + const auto& hcalGeometry = getCondition( + ldmx::HcalGeometry::CONDITIONS_OBJECT_NAME); + + SetupRandomNumberGeneration(); + + std::vector hcalRecHits; + + auto simHits{event.getCollection(input_coll_name, + input_pass_name)}; + std::unordered_map> + hits_by_id{}; + // Important, has to be a reference so that we don't take the address of a + // variable that goes out of scope! + for (const auto& hit : simHits) { + auto id{hit.getID()}; + auto found{hits_by_id.find(id)}; + if (found == hits_by_id.end()) { + hits_by_id[id] = std::vector{&hit}; + } else { + hits_by_id[id].push_back(&hit); + } + } + for (const auto& [barID, simhits_in_bar] : hits_by_id) { + ldmx::HcalHit recHit{}; + double edep{}; + double time{}; + std::vector pos{0, 0, 0}; + for (auto hit : simhits_in_bar) { + edep += hit->getEdep(); + time += hit->getTime() * edep; + auto hitPos{hit->getPosition()}; + pos[0] += hitPos[0] * edep; + pos[1] += hitPos[1] * edep; + pos[2] += hitPos[2] * edep; + } + ldmx::HcalID hitID{barID}; + + // Position smearing + double mean_pe{(edep / mev_per_mip) * pe_per_mip}; + double xpos{pos[0] / edep}; + double ypos{pos[1] / edep}; + double zpos{pos[2] / edep}; + time /= edep; + + auto orientation{hcalGeometry.getScintillatorOrientation(barID)}; + double half_total_width{ + hcalGeometry.getHalfTotalWidth(hitID.section(), hitID.layer())}; + double scint_bar_length{hcalGeometry.getScintillatorLength(hitID)}; + + auto stripCenter{hcalGeometry.getStripCenterPosition(hitID)}; + if (hitID.section() == ldmx::HcalID::HcalSection::BACK) { + double distance_along_bar = + (orientation == + ldmx::HcalGeometry::ScintillatorOrientation::horizontal) + ? xpos + : ypos; + if (orientation == + ldmx::HcalGeometry::ScintillatorOrientation::horizontal) { + ypos = stripCenter.y(); + xpos += random->Gaus(0, position_resolution); + } else { + xpos = stripCenter.x(); + ypos += random->Gaus(0, position_resolution); + } + zpos = stripCenter.z(); + // Attenuation + mean_pe *= exp(1. / attenuation_length); + double mean_pe_close = + mean_pe * exp(-1. * + ((half_total_width - distance_along_bar) / + (scint_bar_length * 0.5)) / + attenuation_length); + double mean_pe_far = + mean_pe * exp(-1. * + ((half_total_width + distance_along_bar) / + (scint_bar_length * 0.5)) / + attenuation_length); + int PE_close{random->Poisson(mean_pe_close + mean_noise)}; + int PE_far{random->Poisson(mean_pe_far + mean_noise)}; + recHit.setPE(PE_close + PE_far); + recHit.setMinPE(std::min(PE_close, PE_far)); + } else { + // Side HCAL, no attenuation business since single ended readout + int PE{random->Poisson(mean_pe + mean_noise)}; + recHit.setPE(PE); + recHit.setMinPE(PE); + // TODO: Look into this + xpos = stripCenter.x(); + ypos = stripCenter.y(); + zpos = stripCenter.z(); + } + + recHit.setID(hitID.raw()); + recHit.setXPos(xpos); + recHit.setNoise(false); + recHit.setYPos(ypos); + recHit.setZPos(zpos); + recHit.setTime(time); + recHit.setSection(hitID.section()); + recHit.setStrip(hitID.strip()); + recHit.setLayer(hitID.layer()); + recHit.setEnergy(edep); + hcalRecHits.push_back(recHit); + } + event.add(output_coll_name, hcalRecHits); +} + +} // namespace hcal + +DECLARE_PRODUCER_NS(hcal, HcalSimpleDigiAndRecProducer); From 72af7ae6bec189d08e5350ddd04b7b8eba426075 Mon Sep 17 00:00:00 2001 From: Cristina Mantilla Suarez Date: Fri, 3 May 2024 12:30:26 -0500 Subject: [PATCH 2/3] Update src/Hcal/HcalSimpleDigiAndRecProducer.cxx Co-authored-by: Tom Eichlersmith <31970302+tomeichlersmith@users.noreply.github.com> --- src/Hcal/HcalSimpleDigiAndRecProducer.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Hcal/HcalSimpleDigiAndRecProducer.cxx b/src/Hcal/HcalSimpleDigiAndRecProducer.cxx index 6699afd..961591d 100644 --- a/src/Hcal/HcalSimpleDigiAndRecProducer.cxx +++ b/src/Hcal/HcalSimpleDigiAndRecProducer.cxx @@ -61,7 +61,7 @@ void HcalSimpleDigiAndRecProducer::produce(framework::Event& event) { } } for (const auto& [barID, simhits_in_bar] : hits_by_id) { - ldmx::HcalHit recHit{}; + ldmx::HcalHit& recHit = hcalRecHits.emplace_back(); double edep{}; double time{}; std::vector pos{0, 0, 0}; From 31924bb69087a3d80708b528383f27eafed108d5 Mon Sep 17 00:00:00 2001 From: Cristina Mantilla Suarez Date: Fri, 3 May 2024 12:30:36 -0500 Subject: [PATCH 3/3] Update src/Hcal/HcalSimpleDigiAndRecProducer.cxx Co-authored-by: Tom Eichlersmith <31970302+tomeichlersmith@users.noreply.github.com> --- src/Hcal/HcalSimpleDigiAndRecProducer.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/src/Hcal/HcalSimpleDigiAndRecProducer.cxx b/src/Hcal/HcalSimpleDigiAndRecProducer.cxx index 961591d..9189e19 100644 --- a/src/Hcal/HcalSimpleDigiAndRecProducer.cxx +++ b/src/Hcal/HcalSimpleDigiAndRecProducer.cxx @@ -140,7 +140,6 @@ void HcalSimpleDigiAndRecProducer::produce(framework::Event& event) { recHit.setStrip(hitID.strip()); recHit.setLayer(hitID.layer()); recHit.setEnergy(edep); - hcalRecHits.push_back(recHit); } event.add(output_coll_name, hcalRecHits); }