HitCreator.cxx

Go to the documentation of this file.

// declaration header
#include "lardata/ArtDataHelper/HitCreator.h"

// C/C++ standard library
#include <utility> // std::move()
#include <algorithm> // std::accumulate(), std::max()
#include <limits> // std::numeric_limits<>
#include <cassert>

// art libraries
#include "canvas/Utilities/Exception.h"
#include "canvas/Persistency/Common/FindOneP.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"

// LArSoft libraries
#include "larcore/Geometry/Geometry.h"
// #include "RawData/RawDigit.h"
#include "lardataobj/RecoBase/Wire.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardata/Utilities/MakeIndex.h"


namespace {
  
  template <typename Left, typename Right, typename Metadata>
  void ClearAssociations(art::Assns<Left, Right, Metadata>& assns) {
    art::Assns<Left, Right, Metadata> empty;
    assns.swap(empty);
  } // ClearAssociations()
  
} // local namespace


namespace recob {
  
  //****************************************************************************
  //***  HitCreator
  //----------------------------------------------------------------------
  HitCreator::HitCreator(
    raw::RawDigit const& digits,
    geo::WireID const&   wireID,
    raw::TDCtick_t       start_tick,
    raw::TDCtick_t       end_tick,
    float                rms,
    float                peak_time,
    float                sigma_peak_time,
    float                peak_amplitude,
    float                sigma_peak_amplitude,
    float                hit_integral,
    float                hit_sigma_integral,
    float                summedADC,
    short int            multiplicity,
    short int            local_index,
    float                goodness_of_fit,
    int                  dof
    ):
    hit(
      digits.Channel(),
      start_tick,
      end_tick,
      peak_time,
      sigma_peak_time,
      rms,
      peak_amplitude,
      sigma_peak_amplitude,
      summedADC,
      hit_integral,
      hit_sigma_integral,
      multiplicity,
      local_index,
      goodness_of_fit,
      dof,
      art::ServiceHandle<geo::Geometry>()->View(digits.Channel()),
      art::ServiceHandle<geo::Geometry>()->SignalType(digits.Channel()),
      wireID
      )
  {} // HitCreator::HitCreator(RawDigit)
  
  
  //----------------------------------------------------------------------
  HitCreator::HitCreator(
    recob::Wire const&   wire,
    geo::WireID const&   wireID,
    raw::TDCtick_t       start_tick,
    raw::TDCtick_t       end_tick,
    float                rms,
    float                peak_time,
    float                sigma_peak_time,
    float                peak_amplitude,
    float                sigma_peak_amplitude,
    float                hit_integral,
    float                hit_sigma_integral,
    float                summedADC,
    short int            multiplicity,
    short int            local_index,
    float                goodness_of_fit,
    int                  dof
    ):
    hit(
      wire.Channel(),
      start_tick,
      end_tick,
      peak_time,
      sigma_peak_time,
      rms,
      peak_amplitude,
      sigma_peak_amplitude,
      summedADC,
      hit_integral,
      hit_sigma_integral,
      multiplicity,
      local_index,
      goodness_of_fit,
      dof,
      wire.View(),
      art::ServiceHandle<geo::Geometry>()->SignalType(wire.Channel()),
      wireID
      )
  {} // HitCreator::HitCreator(Wire)
  
  
  //----------------------------------------------------------------------
  HitCreator::HitCreator(
    recob::Wire const&   wire,
    geo::WireID const&   wireID,
    raw::TDCtick_t       start_tick,
    raw::TDCtick_t       end_tick,
    float                rms,
    float                peak_time,
    float                sigma_peak_time,
    float                peak_amplitude,
    float                sigma_peak_amplitude,
    float                hit_integral,
    float                hit_sigma_integral,
    short int            multiplicity,
    short int            local_index,
    float                goodness_of_fit,
    int                  dof
    ):
    HitCreator(
      wire, wireID, start_tick, end_tick,
      rms, peak_time, sigma_peak_time, peak_amplitude, sigma_peak_amplitude,
      hit_integral, hit_sigma_integral,
      std::accumulate(
        wire.SignalROI().begin() + start_tick,
        wire.SignalROI().begin() + end_tick,
        0.
        ), // sum of ADC counts between start_tick and end_tick
      multiplicity, local_index,
      goodness_of_fit, dof
      )
  {} // HitCreator::HitCreator(Wire; no summed ADC)
  
  
  //----------------------------------------------------------------------
  HitCreator::HitCreator(
    recob::Wire const&        wire,
    geo::WireID const&        wireID,
    float                     rms,
    float                     peak_time,
    float                     sigma_peak_time,
    float                     peak_amplitude,
    float                     sigma_peak_amplitude,
    float                     hit_integral,
    float                     hit_sigma_integral,
    float                     summedADC,
    short int                 multiplicity,
    short int                 local_index,
    float                     goodness_of_fit,
    int                       dof,
    RegionOfInterest_t const& signal
    ):
    HitCreator(
      wire, wireID, signal.begin_index(), signal.end_index(),
      rms, peak_time, sigma_peak_time, peak_amplitude, sigma_peak_amplitude,
      hit_integral, hit_sigma_integral, summedADC, multiplicity, local_index,
      goodness_of_fit, dof
      )
  {} // HitCreator::HitCreator(Wire; RoI)
  
  
  //----------------------------------------------------------------------
  HitCreator::HitCreator(
    recob::Wire const&        wire,
    geo::WireID const&        wireID,
    float                     rms,
    float                     peak_time,
    float                     sigma_peak_time,
    float                     peak_amplitude,
    float                     sigma_peak_amplitude,
    float                     hit_integral,
    float                     hit_sigma_integral,
    float                     summedADC,
    short int                 multiplicity,
    short int                 local_index,
    float                     goodness_of_fit,
    int                       dof,
    size_t                    iSignalRoI
    ):
    HitCreator(
      wire, wireID, rms, peak_time, sigma_peak_time, peak_amplitude, sigma_peak_amplitude,
      hit_integral, hit_sigma_integral, summedADC, multiplicity, local_index,
      goodness_of_fit, dof, wire.SignalROI().range(iSignalRoI)
      )
  {} // HitCreator::HitCreator(Wire; RoI index)
  
  
  HitCreator::HitCreator(recob::Hit const& from): hit(from) {}
  
  
  HitCreator::HitCreator(recob::Hit const& from, geo::WireID const& wireID):
    hit(from)
  {
    hit.fWireID = wireID;
  } // HitCreator::HitCreator(new wire ID)
  
  
  
  //****************************************************************************
  //***  HitAndAssociationsWriterBase
  //----------------------------------------------------------------------
  void HitAndAssociationsWriterBase::put_into() {
    assert(event);
    if (hits) event->put(std::move(hits), prod_instance);
    if (WireAssns) event->put(std::move(WireAssns), prod_instance);
    if (RawDigitAssns) event->put(std::move(RawDigitAssns), prod_instance);
  } // HitAndAssociationsWriterBase::put_into()
  
  
  //****************************************************************************
  //***  HitCollectionCreator
  //----------------------------------------------------------------------
  void HitCollectionCreator::emplace_back(
    recob::Hit&& hit,
    art::Ptr<recob::Wire> const& wire, art::Ptr<raw::RawDigit> const& digits
  ) {
    
    // add the hit to the collection
    hits->emplace_back(std::move(hit));
    
    CreateAssociationsToLastHit(wire, digits);
  } // HitCollectionCreator::emplace_back(Hit&&)
  
  
  //----------------------------------------------------------------------
  void HitCollectionCreator::emplace_back(
    recob::Hit const& hit,
    art::Ptr<recob::Wire> const& wire, art::Ptr<raw::RawDigit> const& digits
  ) {
    
    // add the hit to the collection
    hits->push_back(hit);
    
    CreateAssociationsToLastHit(wire, digits);
  } // HitCollectionCreator::emplace_back(Hit)
  
  
  //----------------------------------------------------------------------
  void HitCollectionCreator::put_into() {
    if (!hits) {
      throw art::Exception(art::errors::LogicError)
        << "HitCollectionCreator is trying to put into the event"
        " a hit collection that was never created!\n";
    }
    HitAndAssociationsWriterBase::put_into();
  } // HitCollectionCreator::put_into()
  
  
  //----------------------------------------------------------------------
  void HitCollectionCreator::CreateAssociationsToLastHit(
    art::Ptr<recob::Wire> const& wire, art::Ptr<raw::RawDigit> const& digits
  ) {
    // if no association is required, we are done
    if (!WireAssns && !RawDigitAssns) return;
    
    // art pointer to the hit we just created
    HitPtr_t hit_ptr(CreatePtrToLastHit());
    
    // association with wires
    if (WireAssns && wire.isNonnull())
      WireAssns->addSingle(wire, hit_ptr); // if it fails, it throws
    
    // association with wires
    if (RawDigitAssns && digits.isNonnull())
      RawDigitAssns->addSingle(digits, hit_ptr); // if it fails, it throws
    
  } // HitCollectionCreator::CreateAssociationsToLastHit()
  
  
  //****************************************************************************
  //***  HitCollectionAssociator
  //----------------------------------------------------------------------
  void HitCollectionAssociator::use_hits
    (std::unique_ptr<std::vector<recob::Hit>>&& srchits)
  {
    hits = std::move(srchits);
  } // HitCollectionAssociator::use_hits()
  
  
  //----------------------------------------------------------------------
  void HitCollectionAssociator::put_into() {
    prepare_associations();
    HitAndAssociationsWriterBase::put_into();
  } // HitCollectionAssociator::put_into()
  
  
  //----------------------------------------------------------------------
  void HitCollectionAssociator::prepare_associations
    (std::vector<recob::Hit> const& srchits)
  {
    if (!RawDigitAssns && !WireAssns) return; // no associations needed
    assert(event);
    
    // we make the associations anew
    if (RawDigitAssns) ClearAssociations(*RawDigitAssns);
    if (WireAssns)     ClearAssociations(*WireAssns);
    
    // the following is true is we want associations with digits
    // but we don't know where digits are; in that case, we try to use wires
    const bool bUseWiresForDigits = RawDigitAssns && (digits_label == "");
    
    if (WireAssns || bUseWiresForDigits) {
      // do we use wires for digit associations too?
      
      // get the wire collection
      art::ValidHandle<std::vector<recob::Wire>> hWires
        = event->getValidHandle<std::vector<recob::Wire>>(wires_label);
      
      // fill a map of wire index vs. channel number
      std::vector<size_t> WireMap
        = util::MakeIndex(*hWires, std::mem_fn(&recob::Wire::Channel));
      
      // use raw rigit - wire association, assuming they have been produced
      // by the same producer as the wire and with the same instance name;
      // we don't check whether the data product is found, but the following
      // code will have FindOneP throw if that was not the case
      // (that's what we would do here anyway, maybe with a better message...)
      std::unique_ptr<art::FindOneP<raw::RawDigit>> WireToDigit;
      if (bUseWiresForDigits) {
        WireToDigit.reset
          (new art::FindOneP<raw::RawDigit>(hWires, *event, wires_label));
      }
      
      // add associations, hit by hit:
      for (size_t iHit = 0; iHit < srchits.size(); ++iHit) {
        
        // find the channel
        size_t iChannel = size_t(srchits[iHit].Channel()); // forcibly converted
        
        // find the wire associated to that channel
        size_t iWire = std::numeric_limits<size_t>::max();
        if (iChannel < WireMap.size()) iWire = WireMap[iChannel];
        if (iWire == std::numeric_limits<size_t>::max()) {
          throw art::Exception(art::errors::LogicError)
            << "No wire associated to channel #" << iChannel << " whence hit #"
            << iHit << " comes!\n";
        } // if no channel
        
        // make the association with wires
        if (WireAssns) {
          art::Ptr<recob::Wire> wire(hWires, iWire);
          WireAssns->addSingle(wire, CreatePtr(iHit));
        }
        
        if (bUseWiresForDigits) {
          // find the digit associated to that channel
          art::Ptr<raw::RawDigit> const& digit = WireToDigit->at(iWire);
          if (digit.isNull()) {
            throw art::Exception(art::errors::LogicError)
              << "No raw digit associated to channel #" << iChannel
              << " whence hit #" << iHit << " comes!\n";
          } // if no channel
          
          // make the association
          RawDigitAssns->addSingle(digit, CreatePtr(iHit));
        } // if create digit associations through wires
      } // for hit
      
    } // if wire associations
    
    if (RawDigitAssns && !bUseWiresForDigits) {
      // get the digit collection
      art::ValidHandle<std::vector<raw::RawDigit>> hDigits
        = event->getValidHandle<std::vector<raw::RawDigit>>(digits_label);
      
      // fill a map of wire index vs. channel number
      std::vector<size_t> DigitMap
        = util::MakeIndex(*hDigits, std::mem_fn(&raw::RawDigit::Channel));
      
      // add associations, hit by hit:
      for (size_t iHit = 0; iHit < srchits.size(); ++iHit) {
        
        // find the channel
        size_t iChannel = size_t(srchits[iHit].Channel()); // forcibly converted
        
        // find the digit associated to that channel
        size_t iDigit = std::numeric_limits<size_t>::max();
        if (iChannel < DigitMap.size()) iDigit = DigitMap[iChannel];
        if (iDigit == std::numeric_limits<size_t>::max()) {
          throw art::Exception(art::errors::LogicError)
            << "No raw digit associated to channel #" << iChannel
            << " whence hit #" << iHit << " comes!\n";
        } // if no channel
        
        // make the association
        art::Ptr<raw::RawDigit> digit(hDigits, iDigit);
        RawDigitAssns->addSingle(digit, CreatePtr(iHit));
        
      } // for hit
    } // if we have rawdigit label
    
  } // HitCollectionAssociator::put_into()
  
  
  //****************************************************************************
  //***  HitRefinerAssociator
  //----------------------------------------------------------------------
  void HitRefinerAssociator::use_hits
    (std::unique_ptr<std::vector<recob::Hit>>&& srchits)
  {
    hits = std::move(srchits);
  } // HitRefinerAssociator::use_hits()
  
  
  //----------------------------------------------------------------------
  void HitRefinerAssociator::put_into() {
    prepare_associations();
    HitAndAssociationsWriterBase::put_into();
  } // HitRefinerAssociator::put_into()
  
  
  //----------------------------------------------------------------------
  void HitRefinerAssociator::prepare_associations
    (std::vector<recob::Hit> const& srchits)
  {
    if (!RawDigitAssns && !WireAssns) return; // no associations needed
    assert(event);
    
    // we make the associations anew
    if (RawDigitAssns) ClearAssociations(*RawDigitAssns);
    
    // read the hits; this is going to hurt performances...
    // no solution to that until there is a way to have a lazy read
    art::ValidHandle<std::vector<recob::Hit>> hHits
      = event->getValidHandle<std::vector<recob::Hit>>(hits_label);
    
    // now get the associations
    if (WireAssns) {
      // we make the associations anew
      ClearAssociations(*WireAssns);
      
      // find the associations between the hits and the wires
      art::FindOneP<recob::Wire> HitToWire(hHits, *event, hits_label);
      if (!HitToWire.isValid()) {
        throw art::Exception(art::errors::ProductNotFound)
          << "Can't find the associations between hits and wires produced by '"
          << hits_label << "'!\n";
      } // if no association
      
      // fill a map of wire vs. channel number
      std::vector<art::Ptr<recob::Wire>> WireMap;
      for (size_t iAssn = 0; iAssn < HitToWire.size(); ++iAssn) {
        art::Ptr<recob::Wire> wire = HitToWire.at(iAssn);
        if (wire.isNull()) continue;
        size_t channelID = (size_t) wire->Channel();
        if (WireMap.size() <= channelID) // expand the map of necessary
          WireMap.resize(std::max(channelID + 1, 2 * WireMap.size()), {});
        WireMap[channelID] = std::move(wire);
      } // for
      
      // now go through all the hits...
      for (size_t iHit = 0; iHit < srchits.size(); ++iHit) {
        recob::Hit const& hit = srchits[iHit];
        size_t channelID = (size_t) hit.Channel();
        
        // no association if there is no wire to associate with
        if ((channelID >= WireMap.size()) || !WireMap[channelID]) continue;
        
        // create an association using the same wire pointer
        WireAssns->addSingle(WireMap[channelID], CreatePtr(iHit));
      } // for hits
    } // if wire associations
    
    // now get the associations
    if (RawDigitAssns) {
      // we make the associations anew
      ClearAssociations(*RawDigitAssns);
      
      // find the associations between the hits and the raw digits
      art::FindOneP<raw::RawDigit> HitToDigits(hHits, *event, hits_label);
      if (!HitToDigits.isValid()) {
        throw art::Exception(art::errors::ProductNotFound)
          << "Can't find the associations between hits and raw digits"
          << " produced by '" << hits_label << "'!\n";
      } // if no association
      
      // fill a map of digits vs. channel number
      std::vector<art::Ptr<raw::RawDigit>> DigitMap;
      for (size_t iAssn = 0; iAssn < HitToDigits.size(); ++iAssn) {
        art::Ptr<raw::RawDigit> digits = HitToDigits.at(iAssn);
        if (digits.isNull()) continue;
        size_t channelID = (size_t) digits->Channel();
        if (DigitMap.size() <= channelID) // expand the map of necessary
          DigitMap.resize(std::max(channelID + 1, 2 * DigitMap.size()), {});
        DigitMap[channelID] = std::move(digits);
      } // for
      
      // now go through all the hits...
      for (size_t iHit = 0; iHit < srchits.size(); ++iHit) {
        recob::Hit const& hit = srchits[iHit];
        size_t channelID = (size_t) hit.Channel();
        
        // no association if there is no digits to associate with
        if ((channelID >= DigitMap.size()) || !DigitMap[channelID]) continue;
        
        // create an association using the same digits pointer
        RawDigitAssns->addSingle(DigitMap[channelID], CreatePtr(iHit));
      } // for hits
    } // if digit associations
    
  } // HitRefinerAssociator::put_into()
  
  
  //----------------------------------------------------------------------
} // namespace recob