EcalRecProducer.cxx

Go to the documentation of this file.

 
#include "Ecal/EcalRecProducer.h"
 
#include "DetDescr/EcalGeometry.h"
#include "Ecal/EcalReconConditions.h"
#include "Ecal/Event/EcalHit.h"
#include "Recon/Event/HgcrocDigiCollection.h"
#include "SimCore/Event/SimCalorimeterHit.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");
}
 
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();
  // 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.
 
      charge = (digi.soi().adcT() - the_conditions.adcPedestal(id)) *
               the_conditions.adcGain(id);
 
      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);