CRHitRemoval Class Reference

Inheritance diagram for CRHitRemoval:

Public Member Functions
	CRHitRemoval (fhicl::ParameterSet const &pset)
virtual void	produce (art::Event &e)
virtual void	beginJob ()
	Begin job method. More...
virtual void	endJob ()
	End job method. More...

Private Types
using	HitPtrVector = std::vector< art::Ptr< recob::Hit >>

Private Member Functions
void	collectPFParticleHits (const recob::PFParticle *pfParticle, const art::Handle< std::vector< recob::PFParticle >> &pfParticleHandle, const art::FindManyP< recob::Cluster > &partToClusAssns, const art::FindManyP< recob::Hit > &clusToHitAssns, HitPtrVector &hitVec)
void	copyAllHits (std::vector< art::Ptr< recob::Hit >> &, art::FindOneP< recob::Wire > &, recob::HitCollectionCreator &)
void	copyInTimeHits (std::vector< art::Ptr< recob::Hit >> &, art::FindOneP< recob::Wire > &, recob::HitCollectionCreator &)
void	FilterHits (HitPtrVector &hits, HitPtrVector &used_hits)

Private Attributes
std::vector< std::string >	fCosmicProducerLabels
	List of cosmic tagger producers. More...
std::string	fHitProducerLabel
	The full collection of hits. More...
std::string	fPFParticleProducerLabel
	PFParticle producer. More...
std::vector< std::string >	fTrackProducerLabels
	Track producer. More...
std::vector< std::string >	fAssnProducerLabels
	Track to PFParticle assns producer. More...
std::vector< double >	fCosmicTagThresholds
	Thresholds for tagging. More...
int	fEndTickPadding
	Padding the end tick. More...
int	fDetectorWidthTicks
	Effective drift time in ticks. More...
int	fMinTickDrift
	Starting tick. More...
int	fMaxTickDrift
	Ending tick. More...
int	fMaxOutOfTime
	Max hits that can be out of time before rejecting. More...
int	fNumEvent
	Number of events seen. More...
int	fNumCRRejects
	Number of tracks produced. More...

Detailed Description

Definition at line 48 of file CRHitRemoval_module.cc.

Member Typedef Documentation

using CRHitRemoval::HitPtrVector = std::vector<art::Ptr<recob::Hit>>

private

Definition at line 60 of file CRHitRemoval_module.cc.

Constructor & Destructor Documentation

CRHitRemoval::CRHitRemoval ( fhicl::ParameterSet const &  pset )

explicit

Constructor.

Arguments:

pset - Fcl parameters.

Definition at line 109 of file CRHitRemoval_module.cc.

  : EDProducer{pset}, fNumEvent(0), fNumCRRejects(0)
{
  fCosmicProducerLabels = pset.get<std::vector<std::string>>("CosmicProducerLabels");
  fHitProducerLabel = pset.get<std::string>("HitProducerLabel");
  fPFParticleProducerLabel = pset.get<std::string>("PFParticleProducerLabel");
  fTrackProducerLabels = pset.get<std::vector<std::string>>("TrackProducerLabels");
  fAssnProducerLabels = pset.get<std::vector<std::string>>("AssnProducerLabels");
  fCosmicTagThresholds = pset.get<std::vector<double>>("CosmicTagThresholds");
  fEndTickPadding = pset.get<int>("EndTickPadding", 50);
  fMaxOutOfTime = pset.get<int>("MaxOutOfTime", 4);

  // let HitCollectionCreator declare that we are going to produce
  // hits and associations with wires and raw digits
  // (with no particular product label)
  recob::HitCollectionCreator::declare_products(producesCollector(), "", true, false);

  // Report.
  mf::LogInfo("CRHitRemoval") << "CRHitRemoval configured\n";
}

Member Function Documentation

void CRHitRemoval::beginJob ( )

virtual

Begin job method.

Definition at line 133 of file CRHitRemoval_module.cc.

{
  auto const* geo = lar::providerFrom<geo::Geometry>();
  auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();
  auto const detp =
    art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataForJob(clock_data);

  float const samplingRate = sampling_rate(clock_data);
  float const driftVelocity = detp.DriftVelocity(detp.Efield(), detp.Temperature()); // cm/us

  fDetectorWidthTicks = 2 * geo->DetHalfWidth() / (driftVelocity * samplingRate / 1000);
  fMinTickDrift =
    clock_data.Time2Tick(clock_data.TriggerTime()); // this is the hardware trigger time
  fMaxTickDrift = fMinTickDrift + fDetectorWidthTicks + fEndTickPadding;
}

void CRHitRemoval::collectPFParticleHits ( const recob::PFParticle *  pfParticle, const art::Handle< std::vector< recob::PFParticle >> &  pfParticleHandle, const art::FindManyP< recob::Cluster > &  partToClusAssns, const art::FindManyP< recob::Hit > &  clusToHitAssns, HitPtrVector &  hitVec  )

private

Find all hits in PFParticle hierarchy

Arguments:

pfParticle - the top level PFParticle to have hits removed pfParticleHandle - handle to the PFParticle objects partToClusAssns - list of PFParticle to Cluster associations clusToHitAssns - list of Cluster to Hit associations hitVec - the current list of hits

This recursively called method will remove all hits associated to an input PFParticle and, in addition, will call itself for all daughters of the input PFParticle

Definition at line 419 of file CRHitRemoval_module.cc.

{
  // Recover the clusters associated to the input PFParticle
  std::vector<art::Ptr<recob::Cluster>> clusterVec = partToClusAssns.at(pfParticle->Self());

  int minTick(fMinTickDrift);
  int maxTick(fMaxTickDrift);

  // Loop over the clusters and grab the associated hits
  for (const auto& cluster : clusterVec) {
    std::vector<art::Ptr<recob::Hit>> clusHitVec = clusToHitAssns.at(cluster.key());
    hitVec.insert(hitVec.end(), clusHitVec.begin(), clusHitVec.end());

    int minHitTick = (*std::min_element(clusHitVec.begin(),
                                        clusHitVec.end(),
                                        [](const auto& hit, const auto& min) {
                                          return hit->PeakTimeMinusRMS() < min->PeakTimeMinusRMS();
                                        }))
                       ->PeakTimeMinusRMS();
    int maxHitTick = (*std::max_element(clusHitVec.begin(),
                                        clusHitVec.end(),
                                        [](const auto& hit, const auto& max) {
                                          return hit->PeakTimePlusRMS() > max->PeakTimePlusRMS();
                                        }))
                       ->PeakTimePlusRMS();

    if (minHitTick < minTick) minTick = minHitTick;
    if (maxHitTick < maxTick) maxTick = maxHitTick;
  }

  // Loop over the daughters of this particle and remove their hits as well
  for (const auto& daughterId : pfParticle->Daughters()) {
    art::Ptr<recob::PFParticle> daughter(pfParticleHandle, daughterId);

    collectPFParticleHits(
      daughter.get(), pfParticleHandle, partToClusAssns, clusToHitAssns, hitVec);
  }

  return;
}

void CRHitRemoval::copyAllHits ( std::vector< art::Ptr< recob::Hit >> &  inputHits, art::FindOneP< recob::Wire > &  wireAssns, recob::HitCollectionCreator &  newHitCollection  )

private

Definition at line 466 of file CRHitRemoval_module.cc.

{
  for (const auto& hitPtr : inputHits) {
    art::Ptr<recob::Wire> wire = wireAssns.at(hitPtr.key());

    // just copy it
    newHitCollection.emplace_back(*hitPtr, wire);
  }

  return;
}

void CRHitRemoval::copyInTimeHits ( std::vector< art::Ptr< recob::Hit >> &  inputHits, art::FindOneP< recob::Wire > &  wireAssns, recob::HitCollectionCreator &  newHitCollection  )

private

Definition at line 481 of file CRHitRemoval_module.cc.

{
  for (const auto& hitPtr : inputHits) {
    // Check on out of time hits
    if (hitPtr->PeakTimeMinusRMS() < fMinTickDrift || hitPtr->PeakTimePlusRMS() > fMaxTickDrift)
      continue;

    art::Ptr<recob::Wire> wire = wireAssns.at(hitPtr.key());

    // just copy it
    newHitCollection.emplace_back(*hitPtr, wire);
  }

  return;
}

void CRHitRemoval::endJob ( )

virtual

End job method.

Definition at line 528 of file CRHitRemoval_module.cc.

{
  double aveCRPerEvent = fNumEvent > 0 ? double(fNumCRRejects) / double(fNumEvent) : 0.;

  mf::LogInfo("CRHitRemoval") << "CRHitRemoval statistics:\n"
                              << "  Number of events = " << fNumEvent << "\n"
                              << "  Number of Cosmic Rays found = " << fNumCRRejects << ", "
                              << aveCRPerEvent << " average/event";
}

void CRHitRemoval::FilterHits ( HitPtrVector &  hits, HitPtrVector &  used_hits  )

private

Definition at line 509 of file CRHitRemoval_module.cc.

{
  if (used_hits.size() > 0) {
    // Make sure both hit collections are sorted.
    std::stable_sort(hits.begin(), hits.end());
    std::stable_sort(used_hits.begin(), used_hits.end());

    // Do set difference operation.
    std::vector<art::Ptr<recob::Hit>>::iterator it = std::set_difference(
      hits.begin(), hits.end(), used_hits.begin(), used_hits.end(), hits.begin());

    // Truncate hit collection.
    hits.erase(it, hits.end());
  }
}

void CRHitRemoval::produce ( art::Event &  evt )

virtual

Produce method.

Arguments:

evt - Art event.

This is the primary method. The goal is to produce a list of recob::Hit objects which are a "clean" subset of all hits and which are believed to be due to a neutrino interaction. It does this by considering input CosmicTag objects, relating them to PFParticles/Tracks and removing the hits associated to those objects which are believed to be Cosmic Rays.

Definition at line 163 of file CRHitRemoval_module.cc.

{
  ++fNumEvent;

  // Start by looking up the original hits
  art::Handle<std::vector<recob::Hit>> hitHandle;
  evt.getByLabel(fHitProducerLabel, hitHandle);

  // If there are no hits then there should be no output
  if (!hitHandle.isValid()) return;

  // also get the associated wires
  // Note that we no longer assume RawDigit associations and will not preserve them if they exist
  // we assume they have been created by the same module as the hits
  art::FindOneP<recob::Wire> ChannelHitWires(hitHandle, evt, fHitProducerLabel);

  HitPtrVector ChHits;
  art::fill_ptr_vector(ChHits, hitHandle);

  // this object contains the hit collection
  // and its associations to wires and raw digits:
  recob::HitCollectionCreator hcol(evt, ChannelHitWires.isValid(), false);

  // Now recover thre remaining collections of objects in the event store that we need
  // Recover the PFParticles that are ultimately associated to tracks
  // This will be the key to recovering the hits
  art::Handle<std::vector<recob::PFParticle>> pfParticleHandle;
  evt.getByLabel(fPFParticleProducerLabel, pfParticleHandle);

  // Without a valid collection of PFParticles we can't do the hit removal
  if (!pfParticleHandle.isValid()) {
    copyAllHits(ChHits, ChannelHitWires, hcol);

    // put the hit collection and associations into the event
    hcol.put_into(evt);

    return;
  }

  // Recover the clusters so we can do associations to the hits
  // In theory the clusters come from the same producer as the PFParticles
  art::Handle<std::vector<recob::Cluster>> clusterHandle;
  evt.getByLabel(fPFParticleProducerLabel, clusterHandle);

  // If there are no clusters then something is really wrong
  if (!clusterHandle.isValid()) {
    copyAllHits(ChHits, ChannelHitWires, hcol);

    // put the hit collection and associations into the event
    hcol.put_into(evt);

    return;
  }

  // Build the list of cosmic tag producers, associations and thresholds
  std::vector<art::Handle<std::vector<anab::CosmicTag>>> cosmicHandleVec;
  std::vector<std::map<int, std::vector<art::Ptr<recob::Track>>>> cosmicToTrackVecMapVec;
  std::vector<std::map<int, std::vector<art::Ptr<recob::PFParticle>>>> trackToPFParticleVecMapVec;
  std::vector<double> thresholdVec;

  for (size_t idx = 0; idx != fCosmicProducerLabels.size(); idx++) {
    std::string& handleLabel(fCosmicProducerLabels[idx]);

    art::Handle<std::vector<anab::CosmicTag>> cosmicHandle;
    evt.getByLabel(handleLabel, cosmicHandle);

    if (cosmicHandle.isValid()) {
      // Look up the associations to tracks
      art::FindManyP<recob::Track> cosmicTrackAssns(cosmicHandle, evt, handleLabel);

      // Get the handle to the tracks
      art::Handle<std::vector<recob::Track>> trackHandle;
      evt.getByLabel(fTrackProducerLabels[idx], trackHandle);

      // This should be the case
      if (trackHandle.isValid()) {
        std::map<int, std::vector<art::Ptr<recob::Track>>> cosmicToTrackVecMap;
        std::map<int, std::vector<art::Ptr<recob::PFParticle>>> trackToPFParticleVecMap;

        // Loop through the cosmic tags
        for (size_t cosmicIdx = 0; cosmicIdx < cosmicHandle->size(); cosmicIdx++) {
          art::Ptr<anab::CosmicTag> cosmicTag(cosmicHandle, cosmicIdx);

          std::vector<art::Ptr<recob::Track>> cosmicToTrackVec =
            cosmicTrackAssns.at(cosmicTag.key());

          cosmicToTrackVecMap[cosmicTag.key()] = cosmicToTrackVec;

          art::FindManyP<recob::PFParticle> trackPFParticleAssns(
            trackHandle, evt, fAssnProducerLabels[idx]);

          for (auto& track : cosmicToTrackVec) {
            std::vector<art::Ptr<recob::PFParticle>> trackToPFParticleVec =
              trackPFParticleAssns.at(track.key());

            for (auto& pfParticle : trackToPFParticleVec)
              trackToPFParticleVecMap[track.key()].push_back(pfParticle);
          }
        }

        // Store these collections
        cosmicHandleVec.emplace_back(cosmicHandle);
        cosmicToTrackVecMapVec.emplace_back(cosmicToTrackVecMap);
        trackToPFParticleVecMapVec.emplace_back(trackToPFParticleVecMap);
        thresholdVec.push_back(fCosmicTagThresholds[idx]);
      }
    }
  }

  // No cosmic tags then nothing to do here
  if (cosmicHandleVec.empty()) {
    copyAllHits(ChHits, ChannelHitWires, hcol);

    // put the hit collection and associations into the event
    hcol.put_into(evt);

    return;
  }

  // Now we walk up the remaining list of associations needed to go from CR tags to hits
  // We'll need to go from tracks to PFParticles
  // From PFParticles we go to clusters
  art::FindManyP<recob::Cluster> clusterAssns(pfParticleHandle, evt, fPFParticleProducerLabel);

  // Likewise, recover the collection of associations to hits
  art::FindManyP<recob::Hit> clusterHitAssns(clusterHandle, evt, fPFParticleProducerLabel);

  // No point double counting hits
  std::set<const recob::PFParticle*> taggedSet;

  // Start the identification of hits to remove. The outer loop is over the various producers of
  // the CosmicTag objects we're examininig
  for (size_t idx = 0; idx != cosmicHandleVec.size(); idx++) {
    // Obviously, dereference the handle and associations
    const art::Handle<std::vector<anab::CosmicTag>>& cosmicHandle(cosmicHandleVec[idx]);
    const std::map<int, std::vector<art::Ptr<recob::Track>>>& crTagTrackVec(
      cosmicToTrackVecMapVec[idx]);
    const std::map<int, std::vector<art::Ptr<recob::PFParticle>>>& trackToPFParticleVecMap(
      trackToPFParticleVecMapVec[idx]);

    for (size_t crIdx = 0; crIdx != cosmicHandle->size(); crIdx++) {
      art::Ptr<anab::CosmicTag> cosmicTag(cosmicHandle, crIdx);

      // If this was tagged as a CR muon then we have work to do!
      if (cosmicTag->CosmicScore() > thresholdVec[idx]) {
        // Recover the associated track
        std::vector<art::Ptr<recob::Track>> trackVec(crTagTrackVec.at(cosmicTag.key()));

        // Loop over the tracks (almost always only 1)
        for (const auto& track : trackVec) {
          std::map<int, std::vector<art::Ptr<recob::PFParticle>>>::const_iterator
            trackToPFParticleVecItr = trackToPFParticleVecMap.find(track.key());

          if (trackToPFParticleVecItr != trackToPFParticleVecMap.end()) {
            // Loop through them
            for (const auto& pfParticlePtr : trackToPFParticleVecItr->second) {
              // Get bare pointer
              const recob::PFParticle* pfParticle = pfParticlePtr.get();

              // A cosmic ray must be a primary (by fiat)
              while (!pfParticle->IsPrimary())
                pfParticle =
                  art::Ptr<recob::PFParticle>(pfParticleHandle, pfParticle->Parent()).get();

              // Add to our list of tagged PFParticles
              taggedSet.insert(pfParticle);
            }
          }
        }
      }
    }
  }

  // If no PFParticles have been tagged then nothing to do
  if (!taggedSet.empty()) {
    // This may all seem backwards... but what you want to do is remove all the hits which are associated to tagged
    // cosmic ray tracks/PFParticles and what you want to leave is all other hits. This includes hits on other PFParticles
    // as well as unassociated (as yet) hits.
    // One also needs to deal with a slight complication with hits in 2D which are shared between untagged PFParticles and
    // tagged ones... we should leave these in.

    // All that is left is to go through the PFParticles which have not been tagged and collect up all their hits
    // to output to our new collection
    // Note that this SHOULD take care of the case of shared 2D hits automagically since if the PFParticle has not
    // been tagged and it shares hits we'll pick those up here.
    HitPtrVector taggedHits;
    HitPtrVector untaggedHits;

    // Loop through the PFParticles and build out the list of hits on untagged PFParticle trees
    for (const auto& pfParticle : *pfParticleHandle) {
      // Start with only primaries
      if (!pfParticle.IsPrimary()) continue;

      // Temporary container for these hits
      HitPtrVector tempHits;

      // Find the hits associated to this untagged PFParticle
      collectPFParticleHits(&pfParticle, pfParticleHandle, clusterAssns, clusterHitAssns, tempHits);

      // One more possible chance at identifying tagged hits...
      // Check these hits to see if any lie outside time window
      bool goodHits(true);

      if (taggedSet.find(&pfParticle) != taggedSet.end()) goodHits = false;

      if (goodHits) {
        int nOutOfTime(0);

        for (const auto& hit : tempHits) {
          // Check on out of time hits
          if (hit->PeakTimeMinusRMS() < fMinTickDrift || hit->PeakTimePlusRMS() > fMaxTickDrift)
            nOutOfTime++;

          if (nOutOfTime > fMaxOutOfTime) {
            goodHits = false;
            break;
          }
        }
      }

      if (goodHits)
        std::copy(tempHits.begin(), tempHits.end(), std::back_inserter(untaggedHits));
      else
        std::copy(tempHits.begin(), tempHits.end(), std::back_inserter(taggedHits));
    }

    // First task - remove hits from the tagged hit collection that are inthe untagged hits (shared hits)
    FilterHits(taggedHits, untaggedHits);

    // Now filter the tagged hits from the total hit collection
    FilterHits(ChHits, taggedHits);
  }

  // Copy our new hit collection to the output
  copyInTimeHits(ChHits, ChannelHitWires, hcol);

  // put the hit collection and associations into the event
  hcol.put_into(evt);
}

Member Data Documentation

std::vector<std::string> CRHitRemoval::fAssnProducerLabels

private

Track to PFParticle assns producer.

Definition at line 84 of file CRHitRemoval_module.cc.

std::vector<std::string> CRHitRemoval::fCosmicProducerLabels

private

List of cosmic tagger producers.

Definition at line 80 of file CRHitRemoval_module.cc.

std::vector<double> CRHitRemoval::fCosmicTagThresholds

private

Thresholds for tagging.

Definition at line 86 of file CRHitRemoval_module.cc.

int CRHitRemoval::fDetectorWidthTicks

private

Effective drift time in ticks.

Definition at line 90 of file CRHitRemoval_module.cc.

int CRHitRemoval::fEndTickPadding

private

Padding the end tick.

Definition at line 88 of file CRHitRemoval_module.cc.

std::string CRHitRemoval::fHitProducerLabel

private

The full collection of hits.

Definition at line 81 of file CRHitRemoval_module.cc.

int CRHitRemoval::fMaxOutOfTime

private

Max hits that can be out of time before rejecting.

Definition at line 93 of file CRHitRemoval_module.cc.

int CRHitRemoval::fMaxTickDrift

private

Ending tick.

Definition at line 92 of file CRHitRemoval_module.cc.

int CRHitRemoval::fMinTickDrift

private

Starting tick.

Definition at line 91 of file CRHitRemoval_module.cc.

int CRHitRemoval::fNumCRRejects

private

Number of tracks produced.

Definition at line 97 of file CRHitRemoval_module.cc.

int CRHitRemoval::fNumEvent

private

Number of events seen.

Definition at line 96 of file CRHitRemoval_module.cc.

std::string CRHitRemoval::fPFParticleProducerLabel

private

PFParticle producer.

Definition at line 82 of file CRHitRemoval_module.cc.

std::vector<std::string> CRHitRemoval::fTrackProducerLabels

private

Track producer.

Definition at line 83 of file CRHitRemoval_module.cc.

---

The documentation for this class was generated from the following file:

• CRHitRemoval_module.cc