IndirectHitParticleAssns_tool.cc

Go to the documentation of this file.

#include "larana/T0Finder/AssociationsTools/IHitParticleAssociations.h"

#include "lardataobj/RecoBase/Hit.h"

#include "nusimdata/SimulationBase/MCParticle.h"

#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Utilities/ToolMacros.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

namespace t0 {
  //
  // Class:       IndirectHitParticleAssns
  // Module Type: art tool
  // File:        IndirectHitParticleAssns.h
  //
  //              This provides MC truth information by using output
  //              reco Hit <--> MCParticle associations
  //
  // Configuration parameters:
  //
  // TruncMeanFraction     - the fraction of waveform bins to discard when
  //
  // Created by Tracy Usher (usher@slac.stanford.edu) on November 21, 2017
  //

  class IndirectHitParticleAssns : virtual public IHitParticleAssociations {
  public:
    explicit IndirectHitParticleAssns(fhicl::ParameterSet const& pset);

    // provide for initialization
    void reconfigure(fhicl::ParameterSet const& pset) override;

    void CreateHitParticleAssociations(art::Event&, HitParticleAssociations*) override;

  private:
    std::vector<art::InputTag> fHitModuleLabelVec;
    art::InputTag fMCParticleModuleLabel;
    art::InputTag fHitPartAssnsModuleLabel;
  };

  //----------------------------------------------------------------------------
  IndirectHitParticleAssns::IndirectHitParticleAssns(fhicl::ParameterSet const& pset)
  {
    reconfigure(pset);

    // Report.
    mf::LogInfo("IndirectHitParticleAssns") << "Configured\n";
  }

  //----------------------------------------------------------------------------
  void IndirectHitParticleAssns::reconfigure(fhicl::ParameterSet const& pset)
  {
    fHitPartAssnsModuleLabel = pset.get<art::InputTag>("HitPartAssnsLabel");
    fMCParticleModuleLabel = pset.get<art::InputTag>("MCParticleModuleLabel");
    fHitModuleLabelVec = pset.get<std::vector<art::InputTag>>("HitModuleLabelVec");
  }

  //----------------------------------------------------------------------------
  void IndirectHitParticleAssns::CreateHitParticleAssociations(
    art::Event& evt,
    HitParticleAssociations* hitPartAssns)
  {
    // This function handles the "indirect" creation of hit<-->MCParticle associations by trying to
    // use the previously created Hit<-->MCParticle associations
    // Its pretty much a brute force approach here... but time is short!
    //
    // First step is to recover the preexisting associations

    // Get a handle for the associations...
    auto partHitAssnsHandle = evt.getValidHandle<HitParticleAssociations>(fHitPartAssnsModuleLabel);
    auto mcParticleHandle =
      evt.getValidHandle<std::vector<simb::MCParticle>>(fMCParticleModuleLabel);

    if (!partHitAssnsHandle.isValid() || !mcParticleHandle.isValid()) {
      throw cet::exception("IndirectHitParticleAssns")
        << "===>> NO MCParticle <--> Hit associations found for run/subrun/event: " << evt.run()
        << "/" << evt.subRun() << "/" << evt.id().event();
    }

    // Loop over input hit collections
    for (const auto& inputTag : fHitModuleLabelVec) {
      // Look up the hits we want to process as well, since if they are not there then no point in proceeding
      art::Handle<std::vector<recob::Hit>> hitListHandle;
      evt.getByLabel(inputTag, hitListHandle);

      if (!hitListHandle.isValid()) {
        mf::LogInfo("IndirectHitParticleAssns")
          << "===>> NO Hit collection found to process for run/subrun/event: " << evt.run() << "/"
          << evt.subRun() << "/" << evt.id().event() << "\n";
        continue;
      }

      // Go through the associations and build out our (hopefully sparse) data structure
      using ParticleDataPair = std::pair<size_t, const anab::BackTrackerHitMatchingData*>;
      using MCParticleDataSet = std::set<ParticleDataPair>;
      using TickToPartDataMap = std::unordered_map<raw::TDCtick_t, MCParticleDataSet>;
      using ChannelToTickPartDataMap = std::unordered_map<raw::ChannelID_t, TickToPartDataMap>;

      ChannelToTickPartDataMap chanToTickPartDataMap;

      // Build out the maps between hits/particles
      for (HitParticleAssociations::const_iterator partHitItr = partHitAssnsHandle->begin();
           partHitItr != partHitAssnsHandle->end();
           ++partHitItr) {
        const art::Ptr<simb::MCParticle>& mcParticle = partHitItr->first;
        const art::Ptr<recob::Hit>& recoHit = partHitItr->second;
        const anab::BackTrackerHitMatchingData* data = &partHitAssnsHandle->data(partHitItr);

        TickToPartDataMap& tickToPartDataMap = chanToTickPartDataMap[recoHit->Channel()];

        for (raw::TDCtick_t tick = recoHit->PeakTimeMinusRMS(); tick <= recoHit->PeakTimePlusRMS();
             tick++) {
          tickToPartDataMap[tick].insert(ParticleDataPair(mcParticle.key(), data));
        }
      }

      // Armed with the map, process the hit list
      for (size_t hitIdx = 0; hitIdx < hitListHandle->size(); hitIdx++) {
        art::Ptr<recob::Hit> hit(hitListHandle, hitIdx);

        TickToPartDataMap& tickToPartDataMap = chanToTickPartDataMap[hit->Channel()];

        if (tickToPartDataMap.empty()) {
          mf::LogInfo("IndirectHitParticleAssns")
            << "No channel information found for hit " << hit << "\n";
          continue;
        }

        // Keep track of results
        MCParticleDataSet particleDataSet;

        // Go through the ticks in this hit and recover associations
        for (raw::TDCtick_t tick = hit->PeakTimeMinusRMS(); tick <= hit->PeakTimePlusRMS();
             tick++) {
          TickToPartDataMap::iterator hitInfoItr = tickToPartDataMap.find(tick);

          if (hitInfoItr != tickToPartDataMap.end()) {
            for (const auto& partData : hitInfoItr->second)
              particleDataSet.insert(partData);
          }
        }

        // Now create new associations for the hit in question
        for (const auto& partData : particleDataSet)
          hitPartAssns->addSingle(
            art::Ptr<simb::MCParticle>(mcParticleHandle, partData.first), hit, *partData.second);
      }
    }

    return;
  }

  //----------------------------------------------------------------------------

  DEFINE_ART_CLASS_TOOL(IndirectHitParticleAssns)
}