CRHitRemoval_module.cc

Go to the documentation of this file.

//
// Class:       CRHitRemoval
// Module Type: producer
// File:        CRHitRemoval_module.cc
//
// This module produces RecoBase/Hit objects after removing those
// deemed to be due to CR muons.
//
// Configuration parameters:
//
// CosmicProducerLabels    - a list of cosmic ray producers which should be or'ed
// HitProducerLabel        - the producer of the recob::Hit objects
// PFParticleProducerLabel - the producer of the recob::PFParticles to consider
// TrackProducerLabel      - the producer of the recob::Track objects
// CosmicTagThresholds     - a vector of thresholds to apply to label as cosmic
// EndTickPadding          - # ticks to "pad" the end tick to account for possible
//                           uncertainty in drift velocity
//
// Created by Tracy Usher (usher@slac.stanford.edu) on September 18, 2014
//

#include <cmath>
#include <algorithm>
#include <vector>

#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "larcore/Geometry/Geometry.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Cluster.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/PFParticle.h"
#include "lardata/ArtDataHelper/HitCreator.h"
#include "lardataobj/AnalysisBase/CosmicTag.h"
#include "lardata/Utilities/AssociationUtil.h"
#include "lardata/DetectorInfoServices/DetectorClocksService.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "lardata/DetectorInfoServices/LArPropertiesService.h"

class CRHitRemoval : public art::EDProducer
{
public:
    
    // Copnstructors, destructor.
    explicit CRHitRemoval(fhicl::ParameterSet const & pset);
    virtual ~CRHitRemoval();
    
    // Overrides.
    virtual void reconfigure(fhicl::ParameterSet const & pset);
    virtual void produce(art::Event & e);
    virtual void beginJob();
    virtual void endJob();
    
private:
    // define vector for hits to make sure of uniform use
    using HitPtrVector = std::vector<art::Ptr<recob::Hit>>;
    
    // Methods
    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<raw::RawDigit>&,
                     art::FindOneP<recob::Wire>&,
                     recob::HitCollectionCreator&);
    
    void copyInTimeHits(std::vector< art::Ptr<recob::Hit>>& ,
                        art::FindOneP<raw::RawDigit>&,
                        art::FindOneP<recob::Wire>&,
                        recob::HitCollectionCreator&);
    
    void FilterHits(HitPtrVector& hits, HitPtrVector& used_hits);
    
    // Fcl parameters.
    std::vector<std::string> fCosmicProducerLabels;     
    std::string              fHitProducerLabel;         
    std::string              fPFParticleProducerLabel;  
    std::vector<std::string> fTrackProducerLabels;      
    std::vector<std::string> fAssnProducerLabels;       
    
    std::vector<double>      fCosmicTagThresholds;      
    
    int                      fEndTickPadding;           
    
    int                      fDetectorWidthTicks;       
    int                      fMinTickDrift;             
    int                      fMaxTickDrift;             
    int                      fMaxOutOfTime;             
    
    // Statistics.
    int fNumEvent;        
    int fNumCRRejects;    
};

DEFINE_ART_MODULE(CRHitRemoval)

//----------------------------------------------------------------------------
CRHitRemoval::CRHitRemoval(fhicl::ParameterSet const & pset) :
fNumEvent(0),
fNumCRRejects(0)
{
    reconfigure(pset);
    
    // 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(*this);
    
    // Report.
    mf::LogInfo("CRHitRemoval") << "CRHitRemoval configured\n";
}

//----------------------------------------------------------------------------
CRHitRemoval::~CRHitRemoval()
{}

//----------------------------------------------------------------------------
void CRHitRemoval::reconfigure(fhicl::ParameterSet const & pset)
{
    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);
}

//----------------------------------------------------------------------------
void CRHitRemoval::beginJob()
{
    auto const* detp = lar::providerFrom<detinfo::DetectorPropertiesService>();
    auto const* geo  = lar::providerFrom<geo::Geometry>();
    auto const* ts   = lar::providerFrom<detinfo::DetectorClocksService>();
    
    float samplingRate  = detp->SamplingRate();
    float driftVelocity = detp->DriftVelocity( detp->Efield(), detp->Temperature() ); // cm/us
    
    fDetectorWidthTicks = 2*geo->DetHalfWidth()/(driftVelocity*samplingRate/1000);
    fMinTickDrift       = ts->Time2Tick(ts->TriggerTime()); // this is the hardware trigger time
    fMaxTickDrift       = fMinTickDrift + fDetectorWidthTicks + fEndTickPadding;
}

//----------------------------------------------------------------------------
void CRHitRemoval::produce(art::Event & evt)
{
    ++fNumEvent;
    
    // Start by looking up the original hits
    art::Handle< std::vector<recob::Hit> > hitHandle;
    evt.getByLabel(fHitProducerLabel, hitHandle);
    
    // also get the associated wires and raw digits;
    // we assume they have been created by the same module as the hits
    art::FindOneP<raw::RawDigit> ChannelHitRawDigits(hitHandle, evt, fHitProducerLabel);
    art::FindOneP<recob::Wire>   ChannelHitWires(hitHandle, evt, fHitProducerLabel);
    
    // If there are no hits then there should be no output
    if (!hitHandle.isValid()) return;
    
    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(*this,
                                     evt,
                                     /* doWireAssns */ ChannelHitWires.isValid(),
                                     /* doRawDigitAssns */ ChannelHitRawDigits.isValid()
                                     );
    
    // 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, ChannelHitRawDigits, 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, ChannelHitRawDigits, 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, ChannelHitRawDigits, 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, ChannelHitRawDigits, ChannelHitWires, hcol);
    
    // put the hit collection and associations into the event
    hcol.put_into(evt);
    
}

//----------------------------------------------------------------------------
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)
{
    // 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<raw::RawDigit>&       rawDigitAssns,
                               art::FindOneP<recob::Wire>&         wireAssns,
                               recob::HitCollectionCreator&        newHitCollection)
{
    for(const auto& hitPtr : inputHits)
    {
        art::Ptr<recob::Wire>   wire      = wireAssns.at(hitPtr.key());
        art::Ptr<raw::RawDigit> rawdigits = rawDigitAssns.at(hitPtr.key());
        
        // just copy it
        newHitCollection.emplace_back(*hitPtr, wire, rawdigits);
    }
    
    return;
}

void CRHitRemoval::copyInTimeHits(std::vector< art::Ptr<recob::Hit>>& inputHits,
                                  art::FindOneP<raw::RawDigit>&       rawDigitAssns,
                                  art::FindOneP<recob::Wire>&         wireAssns,
                                  recob::HitCollectionCreator&        newHitCollection)
{
    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());
        art::Ptr<raw::RawDigit> rawdigits = rawDigitAssns.at(hitPtr.key());
        
        // just copy it
        newHitCollection.emplace_back(*hitPtr, wire, rawdigits);
    }
    
    return;
}

//----------------------------------------------------------------------------
// Filter a collection of hits (set difference).
// This function is copied from Track3DKalmanHit_module.cc
//
// Arguments:
//
// hits      - Hit collection from which hits should be removed.
// used_hits - Hits to remove.
//
void CRHitRemoval::FilterHits(HitPtrVector& hits, HitPtrVector& used_hits)
{
    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::endJob()
{
    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";
}