Bringing back the old digi/rec-producer for testing/debugging purposes

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<TRandom3> random{nullptr};
+  std::unique_ptr<ldmx::NoiseGenerator> noiseGenerator{nullptr};
+  int readout_threshold{2};
+};
+
+}  // namespace hcal
+
+#endif /* HCALSIMPLEDIGIANDRECPRODUCER_H */

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

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<std::string>("input_coll_name");
+  input_pass_name = ps.getParameter<std::string>("input_pass_name");
+  output_coll_name = ps.getParameter<std::string>("output_coll_name");
+  mev_per_mip = ps.getParameter<double>("mev_per_mip");
+  pe_per_mip = ps.getParameter<double>("pe_per_mip");
+  attenuation_length = ps.getParameter<double>("attenuation_length");
+  readout_threshold = ps.getParameter<int>("readout_threshold");
+  mean_noise = ps.getParameter<double>("mean_noise");
+  position_resolution = ps.getParameter<double>("position_resolution");
+}
+
+void HcalSimpleDigiAndRecProducer::SetupRandomNumberGeneration() {
+  if (noiseGenerator == nullptr) {
+    noiseGenerator = std::make_unique<ldmx::NoiseGenerator>(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>(
+            framework::RandomNumberSeedService::CONDITIONS_OBJECT_NAME);
+    noiseGenerator->seedGenerator(
+        rseed.getSeed("HcalSimpleDigiAndRecProducer::NoiseGenerator"));
+  }
+  if (random == nullptr) {
+    const framework::RandomNumberSeedService& rseed =
+        getCondition<framework::RandomNumberSeedService>(
+            framework::RandomNumberSeedService::CONDITIONS_OBJECT_NAME);
+    random = std::make_unique<TRandom3>(
+        rseed.getSeed("HcalSimpleDigiAndRecProducer"));
+  }
+}
+
+void HcalSimpleDigiAndRecProducer::produce(framework::Event& event) {
+  const auto& hcalGeometry = getCondition<ldmx::HcalGeometry>(
+      ldmx::HcalGeometry::CONDITIONS_OBJECT_NAME);
+
+  SetupRandomNumberGeneration();
+
+  std::vector<ldmx::HcalHit> hcalRecHits;
+
+  auto simHits{event.getCollection<ldmx::SimCalorimeterHit>(input_coll_name,
+                                                            input_pass_name)};
+  std::unordered_map<unsigned int, std::vector<const ldmx::SimCalorimeterHit*>>
+      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<const ldmx::SimCalorimeterHit*>{&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<double> 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);