EcalRecProducer.cxx

Go to the documentation of this file.

 
#include "Ecal/EcalRecProducer.h"
 
#include <algorithm>
#include <cmath>
 
#include "DetDescr/EcalGeometry.h"
#include "DetDescr/EcalID.h"
#include "Ecal/EcalReconConditions.h"
#include "Ecal/Event/EcalHit.h"
#include "Framework/RandomNumberSeedService.h"
#include "Recon/Event/HgcrocDigiCollection.h"
#include "SimCore/Event/SimCalorimeterHit.h"
#include "Tools/PulseFitter.h"
 
namespace ecal {
 
EcalRecProducer::EcalRecProducer(const std::string& name,
                                 framework::Process& process)
    : Producer(name, process) {}
 
void EcalRecProducer::configure(framework::config::Parameters& ps) {
  // collection names
  digi_coll_name_ = ps.get<std::string>("digi_coll_name");
  digi_pass_name_ = ps.get<std::string>("digi_pass_name");
  sim_hit_coll_name_ = ps.get<std::string>("sim_hit_coll_name");
  sim_hit_pass_name_ = ps.get<std::string>("sim_hit_pass_name");
  rec_hit_coll_name_ = ps.get<std::string>("rec_hit_coll_name");
 
  layer_weights_ = ps.get<std::vector<double>>("layer_weights");
  second_order_energy_correction_ =
      ps.get<double>("second_order_energy_correction");
 
  mip_si_energy_ = ps.get<double>("mip_si_energy");
  clock_cycle_ = ps.get<double>("clock_cycle");
  charge_per_mip_ = ps.get<double>("charge_per_mip");
 
  amplitude_method_ = ps.get<std::string>("amplitude_method", "soi");
  fit_n_samples_ = ps.get<int>("fit_n_samples", 0);
  fit_fixed_t0_ = ps.get<bool>("fit_fixed_t0", false);
  fit_t0_offset_ = ps.get<double>("fit_t0_offset", 0.);
  phase_resolution_ = ps.get<double>("phase_resolution", 0.);
  phase_name_ = ps.get<std::string>("phase_name", "ecal_digi_phase");
  if (amplitude_method_ == "fit") {
    auto shape = ps.get<std::vector<double>>("pulse_shape");
    if (shape.size() != 5) {
      EXCEPTION_RAISE("InvalidConfig", "pulse_shape needs 5 parameters.");
    }
    pulse_shape_ = std::make_unique<ldmx::HgcrocShape>(
        shape[0], shape[1], shape[2], shape[3], shape[4]);
  } else if (amplitude_method_ != "soi" and amplitude_method_ != "max") {
    EXCEPTION_RAISE("InvalidConfig", "Unknown amplitude_method '" +
                                         amplitude_method_ +
                                         "', use soi, max or fit.");
  }
}
 
void EcalRecProducer::onNewRun(const ldmx::RunHeader&) {
  if (fit_fixed_t0_) {
    const auto& rseed = getCondition<framework::RandomNumberSeedService>(
        framework::RandomNumberSeedService::CONDITIONS_OBJECT_NAME);
    rng_.seed(rseed.getSeed(getName() + "::phase"));
  }
}
 
void EcalRecProducer::produce(framework::Event& event) {
  // Get the Ecal Geometry
  const auto& geometry = getCondition<ldmx::EcalGeometry>(
      ldmx::EcalGeometry::CONDITIONS_OBJECT_NAME);
 
  // Safety check: ensure layer_weights_ covers all geometry layers.
  // A mismatch here typically means an outdated geometry configuration
  // (e.g. v14) is being used with layer weights for a different geometry
  // version (e.g. v15).
  if (geometry.getNumLayers() > static_cast<int>(layer_weights_.size())) {
    EXCEPTION_RAISE(
        "InvalidConfig",
        "The number of layers in the Ecal geometry (" +
            std::to_string(geometry.getNumLayers()) +
            ") exceeds the number of configured layer_weights (" +
            std::to_string(layer_weights_.size()) +
            "). This is likely caused by a mismatch between the geometry "
            "version and the reconstruction configuration. Please ensure "
            "that the correct geometry is being used.");
  }
 
  // Get the reconstruction parameters
  EcalReconConditions the_conditions(
      getCondition<conditions::DoubleTableCondition>(
          EcalReconConditions::CONDITIONS_NAME));
 
  std::vector<ldmx::EcalHit> ecal_rec_hits;
  auto ecal_digis = event.getObject<ldmx::HgcrocDigiCollection>(
      digi_coll_name_, digi_pass_name_);
  const int i_soi = ecal_digis.getSampleOfInterestIndex();
 
  // one measured phase per event, shared by all channels
  double t0_fixed{0.};
  if (fit_fixed_t0_) {
    t0_fixed =
        event.getEventHeader().getFloatParameter(phase_name_) + fit_t0_offset_;
    if (phase_resolution_ > 0.) {
      t0_fixed += std::normal_distribution<double>(0., phase_resolution_)(rng_);
    }
  }
  // loop through digis
  for (auto digi : ecal_digis) {
    // ID from first digi sample
    //  assuming rest of samples have same ID
    ldmx::EcalID id(digi.id());
 
    // ID to real space position
    auto [x_, y_, z_] = geometry.getPosition(id);
 
    // the TOT measurement is in the sample where the pulse fell back below
    // threshold, which is the SOI only for in-time hits (iss #1944)
    int i_tot = digi.totSampleIndex();
 
    // get the estimated charge deposited from digi samples
    double charge(0.);
    // TOA with respect to the 25ns clock window that measured it
    double hit_time(0.);
 
    if (i_tot >= 0) {
      // TOT - number of clock ticks that pulse was over threshold
      //  this is related to the amplitude of the pulse approximately through a
      //  linear drain rate the amplitude of the pulse is related to the energy
      //  deposited
      auto tot_sample = digi.at(i_tot);
 
      // convert the time over threshold into a total energy deposited in the
      // silicon
      //  (time over threshold [ns] - pedestal) * gain
      charge = (tot_sample.tot() - the_conditions.totPedestal(id)) *
               the_conditions.totGain(id);
 
      // TOA is measured from the start of its own sample
      hit_time = (i_tot - i_soi) * clock_cycle_ +
                 tot_sample.toa() * (clock_cycle_ / 1024);
 
      ldmx_log(trace) << "Recon { TOA: " << hit_time << " ns } ";
      ldmx_log(trace) << "TOT Mode (sample " << i_tot << ") -> "
                      << tot_sample.tot() << "TDC -> " << charge << " fC";
    } else {
      hit_time = digi.soi().toa() * (clock_cycle_ / 1024);  // ns
      ldmx_log(trace) << "Recon { TOA: " << hit_time << " ns } ";
 
      // ADC mode of readout
      // ADC - voltage measurement at a specific time of the pulse
      // Pulse Shape:
      //  p[0]/(1.0+exp(p[1](t-p[2]+p[3]-p[4])))/(1.0+exp(p[5]*(t-p[6]+p[3]-p[4])))
      //  p[0] = amplitude to be fit (TBD)
      //  p[1] = -0.345 shape parameter - rate of up slope
      //  p[2] = 70.6547 shape parameter - time of up slope relative to shape
      //  fit p[3] = 77.732 shape parameter - time of peak relative to shape fit
      //  p[4] = peak time to be fit (TBD)
      //  p[5] = 0.140068 shape parameter - rate of down slope
      //  p[6] = 87.7649 shape paramter - time of down slope relative to shape
      //  fit
      // These measurements can be used to fit the pulse shape if TOT is not
      // available. For now, we simply take the measurement of the SOI as the
      // peak amplitude.
 
      const double pedestal = the_conditions.adcPedestal(id);
      if (amplitude_method_ == "soi") {
        charge = (digi.soi().adcT() - pedestal) * the_conditions.adcGain(id);
      } else {
        std::vector<short> adc;
        for (unsigned int i{0}; i < digi.size(); i++)
          adc.push_back(static_cast<short>(digi.at(i).adcT()));
        int i_max = std::max_element(adc.begin(), adc.end()) - adc.begin();
        if (amplitude_method_ == "max") {
          charge = (adc[i_max] - pedestal) * the_conditions.adcGain(id);
          // TOA is measured from the start of its own sample
          hit_time = (i_max - i_soi) * clock_cycle_ +
                     digi.at(i_max).toa() * (clock_cycle_ / 1024);
        } else {
          // window starts one sample before the peak to catch the rise
          const int n_adc = static_cast<int>(adc.size());
          const int n =
              fit_n_samples_ > 0 ? std::min(fit_n_samples_, n_adc) : n_adc;
          int i_peak = i_max;
          if (fit_fixed_t0_) {
            i_peak =
                i_soi + static_cast<int>(std::lround(t0_fixed / clock_cycle_));
          }
          const int first = std::clamp(i_peak - 1, 0, n_adc - n);
          std::vector<short> window(adc.begin() + first,
                                    adc.begin() + first + n);
          // peak within half a clock of the window, else noise tails fake
          // pulses
          ldmx::PulseFitter fitter(
              *pulse_shape_, (first - i_soi) * clock_cycle_, clock_cycle_,
              pedestal, 1., (first - i_soi - 0.5) * clock_cycle_,
              (first + n - 1 - i_soi + 0.5) * clock_cycle_, 0.25);
          auto r = fit_fixed_t0_ ? fitter.fitFixed(window, t0_fixed)
                                 : fitter.fit(window);
          charge = r.converged ? r.amplitude * the_conditions.adcGain(id) : 0.;
          hit_time = r.t0;
        }
      }
 
      ldmx_log(trace) << "ADC Mode -> " << charge << " fC";
    }
 
    if (charge <= 0) continue;
 
    double num_mips_equivalent = charge / charge_per_mip_;
    double energy_deposited_in_si = num_mips_equivalent * mip_si_energy_;
 
    ldmx_log(trace) << " -> " << num_mips_equivalent << " equiv MIPs -> "
                    << energy_deposited_in_si << " MeV";
 
    // incorporate layer_ weights
    double reconstructed_energy =
        (num_mips_equivalent *
             layer_weights_.at(
                 id.layer())       // energy lost in non-sensitive layers
         + energy_deposited_in_si  // energy deposited in Si itself
         ) *
        second_order_energy_correction_;
 
    // copy over information to rec hit structure in new collection
    ldmx::EcalHit rec_hit;
    rec_hit.setID(id.raw());
    rec_hit.setXPos(x_);
    rec_hit.setYPos(y_);
    rec_hit.setZPos(z_);
    rec_hit.setAmplitude(energy_deposited_in_si);
    rec_hit.setEnergy(reconstructed_energy);
    rec_hit.setTime(hit_time);
 
    ecal_rec_hits.push_back(rec_hit);
  }
 
  if (event.exists(sim_hit_coll_name_, sim_hit_pass_name_)) {
    // ecal sim hits_ exist ==> label which hits_ are real and which are pure
    // noise
    auto ecal_sim_hits{event.getCollection<ldmx::SimCalorimeterHit>(
        sim_hit_coll_name_, sim_hit_pass_name_)};
    std::set<int> real_hits;
    for (auto const& sim_hit : ecal_sim_hits) real_hits.insert(sim_hit.getID());
    for (auto& hit : ecal_rec_hits)
      hit.setNoise(real_hits.find(hit.getID()) == real_hits.end());
  }
 
  // add collection to event bus
  event.add(rec_hit_coll_name_, ecal_rec_hits);
}
 
}  // namespace ecal
 
DECLARE_PRODUCER(ecal::EcalRecProducer);