diff --git a/Hcal/include/Hcal/HcalSimpleDigiAndRecProducer.h b/Hcal/include/Hcal/HcalSimpleDigiAndRecProducer.h new file mode 100644 index 000000000..a6ab3a77a --- /dev/null +++ b/Hcal/include/Hcal/HcalSimpleDigiAndRecProducer.h @@ -0,0 +1,42 @@ +#ifndef HCALSIMPLEDIGIANDRECPRODUCER_H +#define HCALSIMPLEDIGIANDRECPRODUCER_H +#include + +#include "DetDescr/DetectorID.h" +#include "DetDescr/HcalGeometry.h" +#include "DetDescr/HcalID.h" +#include "Framework/Configure/Parameters.h" +#include "Framework/EventProcessor.h" +#include "Framework/RandomNumberSeedService.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 onNewRun(const ldmx::RunHeader& runHeader) override; + void configure(framework::config::Parameters& ps) override; + 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_{}; + std::mt19937 rng_{}; + std::unique_ptr> position_resolution_smear_{ + nullptr}; + std::unique_ptr random_{nullptr}; + std::unique_ptr noiseGenerator_{nullptr}; + int readout_threshold_{2}; +}; + +} // namespace hcal + +#endif /* HCALSIMPLEDIGIANDRECPRODUCER_H */ diff --git a/Hcal/python/digi.py b/Hcal/python/digi.py index e4b64f57c..b86e9ee46 100644 --- a/Hcal/python/digi.py +++ b/Hcal/python/digi.py @@ -219,3 +219,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/Hcal/src/Hcal/HcalSimpleDigiAndRecProducer.cxx b/Hcal/src/Hcal/HcalSimpleDigiAndRecProducer.cxx new file mode 100644 index 000000000..8e580f270 --- /dev/null +++ b/Hcal/src/Hcal/HcalSimpleDigiAndRecProducer.cxx @@ -0,0 +1,145 @@ +#include "Hcal/HcalSimpleDigiAndRecProducer.h" + +#include "Hcal/Event/HcalHit.h" +#include "SimCore/Event/SimCalorimeterHit.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_smear_ = + std::make_unique>( + 0.0, ps.getParameter("position_resolution")); +} + +void HcalSimpleDigiAndRecProducer::onNewRun(const ldmx::RunHeader&) { + noiseGenerator_ = std::make_unique(mean_noise_, false); + // hard-code this number, create noise hits for non-zero PEs! + noiseGenerator_->setNoiseThreshold(1); + const framework::RandomNumberSeedService& rseed = + getCondition( + framework::RandomNumberSeedService::CONDITIONS_OBJECT_NAME); + noiseGenerator_->seedGenerator( + rseed.getSeed("HcalSimpleDigiAndRecProducer::NoiseGenerator")); + rng_.seed(rseed.getSeed("HcalSimpleDigiAndRecProducer")); +} + +void HcalSimpleDigiAndRecProducer::produce(framework::Event& event) { + const auto& hcalGeometry = getCondition( + ldmx::HcalGeometry::CONDITIONS_OBJECT_NAME); + + 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 = hcalRecHits.emplace_back(); + 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 += (*position_resolution_smear_)(rng_); + } else { + xpos = stripCenter.x(); + ypos += (*position_resolution_smear_)(rng_); + } + 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{ + std::poisson_distribution(mean_pe_close + mean_noise_)(rng_)}; + int PE_far{ + std::poisson_distribution(mean_pe_far + mean_noise_)(rng_)}; + 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{std::poisson_distribution(mean_pe + mean_noise_)(rng_)}; + 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); + recHit.setOrientation(static_cast(orientation)); + } + event.add(output_coll_name_, hcalRecHits); +} + +} // namespace hcal + +DECLARE_PRODUCER_NS(hcal, HcalSimpleDigiAndRecProducer);