MCTruthTestAna Class Reference

Inheritance diagram for MCTruthTestAna:

Public Member Functions
	MCTruthTestAna (fhicl::ParameterSet const &pset)
virtual	~MCTruthTestAna ()
void	beginJob ()
void	endJob ()
void	beginRun (const art::Run &run)
void	reconfigure (fhicl::ParameterSet const &pset)
void	analyze (const art::Event &evt)

Private Attributes
art::InputTag	fHitProducerLabel
art::InputTag	fTrackProducerLabel
std::unique_ptr < truth::IMCTruthMatching >	fMCTruthMatching
TH1D *	fNumIDEHist
TH1D *	fDeltaIDEHist
TH1D *	fNegEnergyHist
TH1D *	fNegNElecHist
TH1D *	fIDEMisMatchHist
TH1D *	fBTEfficiencyHist
TH1D *	fAssocEfficiencyHist
TH2D *	fEfficiencyCompHist
TH1D *	fBTPurityHist
TH1D *	fAssocPurityHist
TH2D *	fPurityCompHist
const geo::GeometryCore *	fGeometry

Detailed Description

Class: MCTruthTestAna Module Type: analyzer File: MCTruthTestAna_module.cc

Author: Tracy Usher E-mail address: usher.nosp@m.@sla.nosp@m.c.sta.nosp@m.nfor.nosp@m.d.edu

This is a test module which will compare the output of the parallel "backtracker" based on Hit <–> MCParticle assocations to the output of the actual BackTracker.

It assumes that the event contains both the simchannel information needed by the BackTracker and the Hit <–> MCParticle associations needed by the alternate backtracker.

Definition at line 79 of file MCTruthTestAna_module.cc.

Constructor & Destructor Documentation

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

explicit

Definition at line 139 of file MCTruthTestAna_module.cc.

    : EDAnalyzer(parameterSet)

{
    fGeometry           = lar::providerFrom<geo::Geometry>();

    // Read in the parameters from the .fcl file.
    this->reconfigure(parameterSet);
}

MCTruthTestAna::~MCTruthTestAna ( )

virtual

Definition at line 151 of file MCTruthTestAna_module.cc.

{}

Member Function Documentation

void MCTruthTestAna::analyze

(

const art::Event &

evt

)

Definition at line 199 of file MCTruthTestAna_module.cc.

{
    // Recover the backtracker
    art::ServiceHandle<cheat::BackTrackerService>       backTracker;
    art::ServiceHandle<cheat::ParticleInventoryService> partInventory;
   
    // "Rebuild" the maps used by the parallel backtracker
    fMCTruthMatching->Rebuild(event);
    
    // Begin the comparisons by simply checking that the number of MCParticles agree between the two
    const sim::ParticleList& particleList = fMCTruthMatching->ParticleList();
    const sim::ParticleList& btPartList   = partInventory->ParticleList();
    
    if (particleList.size() != btPartList.size())
    {
        mf::LogInfo("MCTruthTestAna") << "*** ParticleList size mismatch, BackTracker says: " << btPartList.size() << ", alternate says: " << particleList.size() << std::endl;
    }

    // Flush with the success of passing that test, let's recover the hits and then loop through them
    // and do hit-by-hit comparisons...
    art::Handle<std::vector<recob::Hit> > hitHandle;
    event.getByLabel(fHitProducerLabel, hitHandle);

    auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(event);
    if (hitHandle.isValid())
    {
        int nTotalMisMatches(0);
        int nNegTrackIDs(0);
        
        for(const auto& hit : *hitHandle)
        {
            int nIDEMisMatches(0);
            
            // Check the claimed parentage of the current hit
            std::vector<sim::TrackIDE> trackIDEVec = fMCTruthMatching->HitToTrackID(clockData, hit);
            std::vector<sim::TrackIDE> btTrkIDEVec = backTracker->HitToTrackIDEs(clockData, hit);
            
            int deltaIDs = int(btTrkIDEVec.size()) - int(trackIDEVec.size());
            
            fNumIDEHist->Fill(trackIDEVec.size(),1.);
            fDeltaIDEHist->Fill(deltaIDs, 1.);

            // It can be that the BackTracker is returning IDE's with negative Track ID's and these don't
            // appear in the associations, so bypass them here...
            if (deltaIDs != 0)
            {
                int numNeg = std::count_if(btTrkIDEVec.begin(),btTrkIDEVec.end(),[](const auto& trackIDE){return trackIDE.trackID < 0;});
                
                deltaIDs -= numNeg;
                
                std::sort(btTrkIDEVec.begin(),btTrkIDEVec.end(),[](const auto& left, const auto& right){return left.trackID > right.trackID;});
                
                while(btTrkIDEVec.back().trackID < 0)
                {
                    fNegEnergyHist->Fill(std::min(9.98,double(btTrkIDEVec.back().energy)), 1.);
                    fNegNElecHist->Fill(std::min(99.9,double(btTrkIDEVec.back().numElectrons)), 1.);
                    
                    btTrkIDEVec.pop_back();
                    nNegTrackIDs++;
                }
            }
            
            // Currently assume BackTracker is "correct" so loop over its output and look for match
            for(const auto& btTrkID : btTrkIDEVec)
            {
                std::vector<sim::TrackIDE>::const_iterator trackIDItr = std::find_if(trackIDEVec.begin(),trackIDEVec.end(),[btTrkID](const auto& id){return btTrkID.trackID == id.trackID && btTrkID.energy == id.energy && btTrkID.numElectrons == id.numElectrons;});
                
                if (trackIDItr == trackIDEVec.end()) nIDEMisMatches++;
            }
            
            fIDEMisMatchHist->Fill(nIDEMisMatches, 1.);
            nTotalMisMatches += nIDEMisMatches;
        }
        
        mf::LogInfo("MCTruthTestAna") << "==> Found " << nTotalMisMatches << " between BackTracker and Associations, BT reported " << nNegTrackIDs << " negative track ids \n";
    }
    
    // We're on a roll now!
    // So let's see if we can recover information about tracks
    art::Handle<std::vector<recob::Track>> trackHandle;
    event.getByLabel(fTrackProducerLabel, trackHandle);
    
    if (trackHandle.isValid())
    {
        // We are going to want a vector of art pointers to hits, so build here
        std::vector<art::Ptr<recob::Hit>> hitPtrVector;
        
        art::fill_ptr_vector(hitPtrVector, hitHandle);

        // Recover the associations between tracks and hits
        art::FindManyP<recob::Hit> hitTrackAssns(trackHandle, event, fTrackProducerLabel);
        
        int nBadEffMatches(0);
        int nBadPurMatches(0);

        for(size_t trackIdx = 0; trackIdx < trackHandle->size(); trackIdx++)
        {
            art::Ptr<recob::Track> track(trackHandle,trackIdx);
            
            std::vector<art::Ptr<recob::Hit>> trackHitVec = hitTrackAssns.at(track.key());
            
            std::set<int> trackIDSet = fMCTruthMatching->GetSetOfTrackIDs(clockData, trackHitVec);
            std::set<int> btTrkIDSet = backTracker->GetSetOfTrackIds(clockData, trackHitVec);
            
            // Check for case were we might have negative track IDs from the BackTracker
            if (btTrkIDSet.size() > trackIDSet.size())
            {
                // Note that the container here is a std::set which will be ordered smallest to largest.
                // In this case we want to remove elements from the front until they are positive
                for(std::set<int>::iterator btTrkIDSetItr = btTrkIDSet.begin(); btTrkIDSetItr != btTrkIDSet.end();)
                {
                    if (!(*btTrkIDSetItr < 0)) break;
                    
                    btTrkIDSetItr = btTrkIDSet.erase(btTrkIDSetItr);
                }
            }
            
            if (btTrkIDSet.size() != trackIDSet.size())
                mf::LogDebug("MCTruthTestAna") << "Mismatch in associated track ids, backtracker: " << btTrkIDSet.size() << ", associations: " << trackIDSet.size() << ", track id: " << track.key() << "\n";

            // Ok, turn the "set" of track ID's into a vector of track ID's
            std::vector<int> btTrkIDVec;
            
            for(auto& trackID : btTrkIDSet) btTrkIDVec.emplace_back(trackID);
            
            // And then use this to recover the vectors of hits associated to each MC track
            std::vector<std::vector<art::Ptr<recob::Hit>>> trkHitVecVec;

            for(const auto& tkID : btTrkIDVec)
            {
                std::vector<art::Ptr<recob::Hit>> hitVec = backTracker->TrackIdToHits_Ps(clockData, tkID, hitPtrVector);
                trkHitVecVec.push_back(hitVec);
            }
            // Apply majority logic - we declare the MCParticle with the most hits to be the "winner"
            std::vector<std::vector<art::Ptr<recob::Hit>>>::iterator bestTrkHitVecItr = std::max_element(trkHitVecVec.begin(),trkHitVecVec.end(),[](const auto& a, const auto& b){return a.size() < b.size();});
            
            if (bestTrkHitVecItr == trkHitVecVec.end())
            {
                // I have been watching way too much Star Trek (the original!)
                mf::LogDebug("MCTruthTestAna") << ">>>>>>> ERROR! >>>>>>> ERROR! >>>>>> MUST PURIFY! >>>>>> ERROR!" << "\n";
                continue;
            }
            
            int indexToBestMC = std::distance(trkHitVecVec.begin(),bestTrkHitVecItr);
            int bestMCTrackID = btTrkIDVec.at(indexToBestMC);
            
            std::vector<art::Ptr<recob::Hit>>& bestMCTrackHitVec = *bestTrkHitVecItr;
            
            std::set<int> mcTrackIdxSet = {bestMCTrackID};
            
            double btTrkEffic = backTracker->HitCollectionEfficiency(clockData, mcTrackIdxSet, trackHitVec, bestMCTrackHitVec, geo::k3D);
            double trackEffic = fMCTruthMatching->HitCollectionEfficiency(clockData, mcTrackIdxSet, trackHitVec, bestMCTrackHitVec, geo::k3D);
            
            if (btTrkEffic != trackEffic)
            {
                mf::LogDebug("MCTruthTestAna") << "Efficiency mismatch, track ID: " << bestMCTrackID << ", track: " << track.key() << ", # hits: " << trackHitVec.size() << ", # MC hits: " << bestMCTrackHitVec.size() << ", btTrkEff: " << btTrkEffic << ", trackEff: " << trackEffic << "\n|";
                nBadEffMatches++;
            }
            
            fBTEfficiencyHist->Fill(std::min(1.01,btTrkEffic), 1.);
            fAssocEfficiencyHist->Fill(std::min(1.01,trackEffic), 1.);
            fEfficiencyCompHist->Fill(std::min(1.01,btTrkEffic), std::min(1.01,trackEffic), 1.);
            
            double btTrkPurity = backTracker->HitCollectionPurity(clockData, mcTrackIdxSet, trackHitVec);
            double trackPurity = fMCTruthMatching->HitCollectionPurity(clockData, mcTrackIdxSet, trackHitVec);
            
            if (btTrkPurity != trackPurity)
            {
                mf::LogDebug("MCTruthTestAna") << "Purity mismatch, track ID: " << bestMCTrackID << ", track: " << track.key() << ", # hits: " << trackHitVec.size() << ", # MC hits: " << bestMCTrackHitVec.size() << ", btTrkPurity: " << btTrkPurity << ", trackPurity: " << trackPurity << "\n|";
                nBadPurMatches++;
            }
            
            fBTPurityHist->Fill(std::min(1.01,btTrkPurity), 1.);
            fAssocPurityHist->Fill(std::min(1.01,trackPurity), 1.);
            fPurityCompHist->Fill(std::min(1.01,btTrkPurity),std::min(1.01,trackPurity), 1.);
        }
        
        mf::LogDebug("MCTruthTestAna") << "Event with " << trackHandle->size() << " reconstructed tracks, found " << nBadEffMatches << " efficiency mismatches, " << nBadPurMatches << " purity mismatchs \n";
    }

    return;
}

void MCTruthTestAna::beginJob

(

)

Definition at line 154 of file MCTruthTestAna_module.cc.

{
    // Access ART's TFileService, which will handle creating and writing
    // histograms and n-tuples for us.
    art::ServiceHandle<art::TFileService> tfs;
    
    // The arguments to 'make<whatever>' are the same as those passed
    // to the 'whatever' constructor, provided 'whatever' is a ROOT
    // class that TFileService recognizes.
    
    // Define the histograms. Putting semi-colons around the title
    // causes it to be displayed as the x-axis label if the histogram
    // is drawn.
    fNumIDEHist          = tfs->make<TH1D>("NumIDEs",      ";# ides",                   20,   0.,  20.);
    fDeltaIDEHist        = tfs->make<TH1D>("DeltaIDE",     ";delta",                    20, -10.,  10.);
    fNegEnergyHist       = tfs->make<TH1D>("NegEnergy",    ";energy",                  500,   0.,  10.);
    fNegNElecHist        = tfs->make<TH1D>("NegNumElec",   ";# electrons",             200,   0., 100.);
    fIDEMisMatchHist     = tfs->make<TH1D>("IDEMisMatch",  ";# mismatches",             20,  10.,  10.);
    fBTEfficiencyHist    = tfs->make<TH1D>("BTEfficiency", ";efficiency",              102,   0.,   1.02);
    fAssocEfficiencyHist = tfs->make<TH1D>("AssocEffic",   ";efficiency",              102,   0.,   1.02);
    fEfficiencyCompHist  = tfs->make<TH2D>("Efficiency",   ";BackTrack;Associations",   52,   0.,   1.04, 52, 0., 1.04);
    fBTPurityHist        = tfs->make<TH1D>("BTPurity",     ";Purity",                  102,   0.,   1.02);
    fAssocPurityHist     = tfs->make<TH1D>("AssocPurity",  ";Purity",                  102,   0.,   1.02);
    fPurityCompHist      = tfs->make<TH2D>("Purity",       ";BackTrack;Associations",   52,   0.,   1.04, 52, 0., 1.04);

}

void MCTruthTestAna::beginRun

(

const art::Run &

run

)

Definition at line 182 of file MCTruthTestAna_module.cc.

{}

void MCTruthTestAna::endJob

(

)

Definition at line 382 of file MCTruthTestAna_module.cc.

{
    // Make a call to normalize histograms if so desired

    return;
}

void MCTruthTestAna::reconfigure

(

fhicl::ParameterSet const &

pset

)

Definition at line 185 of file MCTruthTestAna_module.cc.

{
    // Read parameters from the .fcl file. The names in the arguments
    // to p.get<TYPE> must match names in the .fcl file.
    fHitProducerLabel              = pset.get<art::InputTag>("HitModuleLabel");
    fTrackProducerLabel            = pset.get<art::InputTag>("TrackProducerLabel");

    // Get the tool for MC Truth matching
    fMCTruthMatching = art::make_tool<truth::IMCTruthMatching>(pset.get<fhicl::ParameterSet>("MCTruthMatching"));

    return;
}

Member Data Documentation

TH1D* MCTruthTestAna::fAssocEfficiencyHist

private

Definition at line 122 of file MCTruthTestAna_module.cc.

TH1D* MCTruthTestAna::fAssocPurityHist

private

Definition at line 125 of file MCTruthTestAna_module.cc.

TH1D* MCTruthTestAna::fBTEfficiencyHist

private

Definition at line 121 of file MCTruthTestAna_module.cc.

TH1D* MCTruthTestAna::fBTPurityHist

private

Definition at line 124 of file MCTruthTestAna_module.cc.

TH1D* MCTruthTestAna::fDeltaIDEHist

private

Definition at line 117 of file MCTruthTestAna_module.cc.

TH2D* MCTruthTestAna::fEfficiencyCompHist

private

Definition at line 123 of file MCTruthTestAna_module.cc.

const geo::GeometryCore* MCTruthTestAna::fGeometry

private

Definition at line 129 of file MCTruthTestAna_module.cc.

art::InputTag MCTruthTestAna::fHitProducerLabel

private

Definition at line 109 of file MCTruthTestAna_module.cc.

TH1D* MCTruthTestAna::fIDEMisMatchHist

private

Definition at line 120 of file MCTruthTestAna_module.cc.

std::unique_ptr<truth::IMCTruthMatching> MCTruthTestAna::fMCTruthMatching

private

Definition at line 113 of file MCTruthTestAna_module.cc.

TH1D* MCTruthTestAna::fNegEnergyHist

private

Definition at line 118 of file MCTruthTestAna_module.cc.

TH1D* MCTruthTestAna::fNegNElecHist

private

Definition at line 119 of file MCTruthTestAna_module.cc.

TH1D* MCTruthTestAna::fNumIDEHist

private

Definition at line 116 of file MCTruthTestAna_module.cc.

TH2D* MCTruthTestAna::fPurityCompHist

private

Definition at line 126 of file MCTruthTestAna_module.cc.

art::InputTag MCTruthTestAna::fTrackProducerLabel

private

Definition at line 110 of file MCTruthTestAna_module.cc.

---

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

• MCTruthTestAna_module.cc