Inheritance diagram for cluster::ClusterTrackAna:

Public Types
using	ModuleType = EDAnalyzer

Public Member Functions
	ClusterTrackAna (fhicl::ParameterSet const &p)

	ClusterTrackAna (ClusterTrackAna const &)=delete

	ClusterTrackAna (ClusterTrackAna &&)=delete

ClusterTrackAna &	operator= (ClusterTrackAna const &)=delete

ClusterTrackAna &	operator= (ClusterTrackAna &&)=delete

void	analyze (art::Event const &e) override

void	doBeginJob (SharedResources const &resources)

void	doEndJob ()

void	doRespondToOpenInputFile (FileBlock const &fb)

void	doRespondToCloseInputFile (FileBlock const &fb)

void	doRespondToOpenOutputFiles (FileBlock const &fb)

void	doRespondToCloseOutputFiles (FileBlock const &fb)

bool	doBeginRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doEndRun (RunPrincipal &rp, ModuleContext const &mc)

bool	doBeginSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEndSubRun (SubRunPrincipal &srp, ModuleContext const &mc)

bool	doEvent (EventPrincipal &ep, ModuleContext const &mc, std::atomic< std::size_t > &counts_run, std::atomic< std::size_t > &counts_passed, std::atomic< std::size_t > &counts_failed)

ModuleDescription const &	moduleDescription () const

void	setModuleDescription (ModuleDescription const &)

std::array< std::vector< ProductInfo >, NumBranchTypes > const &	getConsumables () const

void	sortConsumables (std::string const &current_process_name)

std::unique_ptr< Worker >	makeWorker (WorkerParams const &wp)

template<typename T , BranchType BT>
ViewToken< T >	consumesView (InputTag const &tag)

template<typename T , BranchType BT>
ViewToken< T >	mayConsumeView (InputTag const &tag)

Protected Member Functions
std::string const &	processName () const

bool	wantAllEvents () const noexcept

bool	wantEvent (ScheduleID id, Event const &e) const

Handle< TriggerResults >	getTriggerResults (Event const &e) const

ConsumesCollector &	consumesCollector ()

template<typename T , BranchType = InEvent>
ProductToken< T >	consumes (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	consumesView (InputTag const &)

template<typename T , BranchType = InEvent>
void	consumesMany ()

template<typename T , BranchType = InEvent>
ProductToken< T >	mayConsume (InputTag const &)

template<typename Element , BranchType = InEvent>
ViewToken< Element >	mayConsumeView (InputTag const &)

template<typename T , BranchType = InEvent>
void	mayConsumeMany ()

Private Member Functions
void	endJob () override

void	FirstLastHitInPlane (unsigned int tpc, unsigned int plane, unsigned int mcpi, unsigned int &firstHitIndex, unsigned int &lastHitIndex)

std::string	HitLocation (unsigned int iht)
	packs hit WireID and PeakTime into a compact format More...

Private Attributes
art::InputTag	fHitModuleLabel

art::InputTag	fClusterModuleLabel

art::InputTag	fTrackModuleLabel

simb::Origin_t	fTruthOrigin

std::vector< int >	fSkipPDGCodes

short	fPrintLevel

unsigned int	fInTPC

float	fBadEP

unsigned int	fEventCnt {0}

unsigned int	fNBadEP {0}

std::array< float, 5 >	Cnts {{0}}

std::array< float, 5 >	EPSums {{0}}

std::array< float, 5 >	ESums {{0}}

std::array< float, 5 >	PSums {{0}}

art::Handle< std::vector< recob::Hit > >	fHitHandle

std::vector< unsigned int >	fHitMCPIndex

bool	fCompareProductIDs
	compare Hit and Cluster-> Hit art product IDs on the first event More...

bool	fFirstPrint {true}

unsigned int	fCurrentRun {0}

Detailed Description

Definition at line 46 of file ClusterTrackAna_module.cc.

Member Typedef Documentation

using art::EDAnalyzer::ModuleType = EDAnalyzer

inherited

Definition at line 22 of file EDAnalyzer.h.

Constructor & Destructor Documentation

cluster::ClusterTrackAna::ClusterTrackAna ( fhicl::ParameterSet const & p )

explicit

Definition at line 99 of file ClusterTrackAna_module.cc.

References Cnts, EPSums, ESums, fBadEP, fClusterModuleLabel, fHitModuleLabel, fInTPC, fPrintLevel, fSkipPDGCodes, fTrackModuleLabel, fTruthOrigin, PSums, and tmp.

                                                                      : EDAnalyzer{pset}
 {
   fHitModuleLabel = pset.get<art::InputTag>("HitModuleLabel");
   bool validModuleLabel = false;
   fClusterModuleLabel = "NA";
   fTrackModuleLabel = "NA";
   fClusterModuleLabel = pset.get<art::InputTag>("ClusterModuleLabel", "NA");
   if (fClusterModuleLabel != "NA") validModuleLabel = true;
   fTrackModuleLabel = pset.get<art::InputTag>("TrackModuleLabel", "NA");
   if (validModuleLabel && fTrackModuleLabel != "NA")
     throw cet::exception("ClusterTrackAna")
       << "You must specify either a ClusterModuleLabel OR a TrackModuleLabel\n";
   if (!validModuleLabel && fTrackModuleLabel != "NA") validModuleLabel = true;
   // origin = 0 (anything), 1(nu), 2(cosmics), 3(SN nu), 4(SingleParticle)
   int tmp = pset.get<int>("TruthOrigin", 0);
   fTruthOrigin = (simb::Origin_t)tmp;
   fPrintLevel = pset.get<short>("PrintLevel", 0);
   if (pset.has_key("SkipPDGCodes")) fSkipPDGCodes = pset.get<std::vector<int>>("SkipPDGCodes");
   fBadEP = pset.get<float>("BadEP", 0.);
   fInTPC = UINT_MAX;
   int intpc = pset.get<int>("InTPC", -1);
   if (intpc >= 0) fInTPC = intpc;
   // do some initialization
   Cnts.fill(0.);
   EPSums.fill(0.);
   ESums.fill(0.);
   PSums.fill(0.);
 } // ClusterTrackAna constructor

cluster::ClusterTrackAna::ClusterTrackAna ( ClusterTrackAna const & )

delete

cluster::ClusterTrackAna::ClusterTrackAna ( ClusterTrackAna && )

delete

Member Function Documentation

void cluster::ClusterTrackAna::analyze ( art::Event const & e )

overridevirtual

Implements art::EDAnalyzer.

Definition at line 129 of file ClusterTrackAna_module.cc.

References util::abs(), sim::ParticleList::begin(), Cnts, simb::MCParticle::E(), util::empty(), sim::ParticleList::end(), EPSums, ESums, art::Event::event(), fBadEP, fClusterModuleLabel, fCompareProductIDs, fCurrentRun, fEventCnt, fFirstPrint, fHitHandle, fHitMCPIndex, fHitModuleLabel, fInTPC, FirstLastHitInPlane(), fNBadEP, fPrintLevel, fSkipPDGCodes, fTrackModuleLabel, fTruthOrigin, Get, art::ProductRetriever::getByLabel(), HitLocation(), cheat::BackTrackerService::HitToTrackIDEs(), art::Handle< T >::id(), ipart, art::Handle< T >::isValid(), simb::kUnknown, art::InputTag::label(), simb::MCParticle::Mass(), simb::MCTruth::Origin(), part, cheat::ParticleInventoryService::ParticleList(), simb::MCParticle::PdgCode(), PSums, art::Event::run(), util::size(), util::to_string(), simb::MCParticle::TrackId(), and cheat::ParticleInventoryService::TrackIdToMCTruth_P().

 {
   // Match hits to MCParticles, then consider reconstructed hits in each TPC and plane
   // to calculate Efficiency, Purity and Efficiency * Purity (aka EP).
 
   ++fEventCnt;
   auto const* geom = lar::providerFrom<geo::Geometry>();
   auto const& wireReadoutGeom = art::ServiceHandle<geo::WireReadout>()->Get();
   //  auto hitsHandle = art::Handle<std::vector<recob::Hit>>();
   if (!evt.getByLabel(fHitModuleLabel, fHitHandle))
     throw cet::exception("ClusterTrackAna")
       << "Failed to get a handle to hit collection '" << fHitModuleLabel.label() << "'\n";
   // get a reference to the MCParticles
   auto mcpHandle = art::Handle<std::vector<simb::MCParticle>>();
   if (!evt.getByLabel("largeant", mcpHandle))
     throw cet::exception("ClusterTrackAna")
       << "Failed to get a handle to MCParticles using largeant\n";
 
   // decide whether to consider cluster -> hit -> MCParticle for any MCParticle origin or for
   // a specific user-specified origin
   bool anySource = (fTruthOrigin == simb::kUnknown);
 
   if (fPrintLevel > 0 && evt.run() != fCurrentRun) {
     mf::LogVerbatim("ClusterTrackAna") << "Run: " << evt.run();
     fCurrentRun = evt.run();
   }
 
   art::ServiceHandle<cheat::BackTrackerService const> bt_serv;
   art::ServiceHandle<cheat::ParticleInventoryService const> pi_serv;
   sim::ParticleList const& plist = pi_serv->ParticleList();
   auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
 
   // make a list of Hit -> MCParticle assns in all TPCs. The first step is
   // to make a list of Geant TrackIDs whose origin was requested by the user.
   std::vector<int> trkIDs;
   // and a vector of MCParticle indices into the mcpHandle vector indexed by trkIDs
   std::vector<unsigned int> MCPIs;
   for (sim::ParticleList::const_iterator ipart = plist.begin(); ipart != plist.end(); ++ipart) {
     const simb::MCParticle* part = (*ipart).second;
     int trackID = part->TrackId();
     art::Ptr<simb::MCTruth> theTruth = pi_serv->TrackIdToMCTruth_P(trackID);
     if (!anySource && theTruth->Origin() != fTruthOrigin) continue;
     int pdg = abs(part->PdgCode());
     bool isCharged = (pdg == 11) || (pdg == 13) || (pdg == 211) || (pdg == 321) || (pdg == 2212);
     if (!isCharged) continue;
     if (std::find(fSkipPDGCodes.begin(), fSkipPDGCodes.end(), pdg) != fSkipPDGCodes.end()) continue;
     float TMeV = 1000 * (part->E() - part->Mass());
     // ignore low energy particles
     if (TMeV < 1) continue;
     // find the MCParticle index and stash it in trkMCP
     unsigned int mcpi = UINT_MAX;
     for (mcpi = 0; mcpi < (*mcpHandle).size(); ++mcpi)
       if ((*mcpHandle)[mcpi].TrackId() == trackID) break;
     if (mcpi == UINT_MAX) {
       std::cout << "Failed to find a MCParticle from ParticleList\n";
       return;
     }
     trkIDs.push_back(trackID);
     MCPIs.push_back(mcpi);
   } // ipart
   if (trkIDs.empty()) return;
 
   // next construct a companion vector of MCParticle indices indexed to the full hit collection
   fHitMCPIndex.resize((*fHitHandle).size(), UINT_MAX);
   // Make a list of the <first, last> MC-matched hit in each TPC. This will be used
   // to only iterate through the range of hits that are interesting
   std::vector<std::pair<unsigned int, unsigned int>> hitRange(geom->NTPC() + 1);
   size_t noMatch = 0;
   size_t nMatch = 0;
   size_t nInTPC = 0;
   for (auto& hr : hitRange)
     hr = std::make_pair(UINT_MAX, UINT_MAX);
   for (size_t iht = 0; iht < (*fHitHandle).size(); ++iht) {
     auto& hit = (*fHitHandle)[iht];
     // only consider hits in a single TPC?
     unsigned int tpc = hit.WireID().TPC;
     if (fInTPC != UINT_MAX && tpc != fInTPC) continue;
     ++nInTPC;
     auto tides = bt_serv->HitToTrackIDEs(clockData, hit);
     if (tides.empty()) {
       ++noMatch;
       continue;
     }
     for (auto& tide : tides) {
       // declare a match to a MCParticle if > 50% of energy is from it
       if (tide.energyFrac < 0.5) continue;
       int trackID = tide.trackID;
       // find the MCParticle index and define the Hit -> MCParticle assn
       for (size_t indx = 0; indx < trkIDs.size(); ++indx) {
         if (trkIDs[indx] == trackID) {
           fHitMCPIndex[iht] = MCPIs[indx];
           break;
         }
       } // indx
       break;
     } // tide
     if (fHitMCPIndex[iht] == UINT_MAX) continue;
     ++nMatch;
     // populate hitRange
     if (hitRange[tpc].first == UINT_MAX) hitRange[tpc].first = iht;
     hitRange[tpc].second = iht;
   } // iht
 
   if (nMatch == 0) {
     fHitMCPIndex.resize(0);
     return;
   }
 
   if (fPrintLevel > 3) {
     // Print the gory details
     mf::LogVerbatim myprt("ClusterTrackAna");
     myprt << "Checking " << trkIDs.size() << " selected MCParticles with the requested TruthOrigin";
     if (fInTPC == UINT_MAX) { myprt << " in all TPCs\n"; }
     else {
       myprt << " in TPC " << fInTPC;
     }
     myprt << " in Run " << evt.run() << " Event " << evt.event();
     myprt << "\n";
     myprt << "Found " << nMatch << " MC-matched hits with the requested TruthOrigin";
     myprt << " out of " << nInTPC << " total hits.\n";
     myprt << "Found " << noMatch << " hits not matched to ANY MCParticle.\n";
     // count the number of hits matched to each MCParticle in the list
     //                     mcpi          count
     std::vector<std::pair<unsigned int, unsigned int>> mcpCnt;
     for (auto mcpi : fHitMCPIndex) {
       if (mcpi == UINT_MAX) continue;
       unsigned short mIndx = 0;
       for (mIndx = 0; mIndx < mcpCnt.size(); ++mIndx)
         if (mcpCnt[mIndx].first == mcpi) break;
       if (mIndx == mcpCnt.size()) mcpCnt.push_back(std::make_pair(mcpi, 0));
       // increment the count of MC-matched hits
       ++mcpCnt[mIndx].second;
     } // mcpi
     myprt << " MCPI TrackID  PDGCode  T(MeV)   nHits    Process";
     for (auto mcpcnt : mcpCnt) {
       myprt << "\n";
       unsigned int mcpi = mcpcnt.first;
       auto& mcp = (*mcpHandle)[mcpi];
       myprt << std::setw(5) << mcpi;
       myprt << std::setw(8) << mcp.TrackId();
       myprt << std::setw(8) << abs(mcp.PdgCode());
       int TMeV = 1000 * (mcp.E() - mcp.Mass());
       myprt << std::setw(9) << TMeV;
       myprt << std::setw(8) << mcpcnt.second;
       myprt << "    " << mcp.Process();
     } // mcpcnt
   }   // fPrintLevel > 3
 
   // fill a vector of hits indices from all clusters or all tracks
   std::vector<std::vector<unsigned int>> recoHits;
   // and a companion vector of indices into the cluster or track collection
   std::vector<unsigned int> recoIndex;
   // handles to these collections to enable printing the cluster or track ID. The
   // handle will be valid for either clusters or tracks
   art::Handle<std::vector<recob::Cluster>> inputClusters;
   art::Handle<std::vector<recob::Track>> inputTracks;
   if (fClusterModuleLabel.label() != "NA") {
     // get MC-matched hits from clusters
     if (!evt.getByLabel(fClusterModuleLabel, inputClusters))
       throw cet::exception("ClusterTrackAna")
         << "Failed to get a handle to the cluster collection '" << fClusterModuleLabel.label()
         << "'\n";
     art::FindManyP<recob::Hit> hitsFromCls(inputClusters, evt, fClusterModuleLabel);
     if (!hitsFromCls.isValid())
       throw cet::exception("ClusterTrackAna") << "Failed to get a handle to Cluster -> Hit assns\n";
     for (unsigned int icl = 0; icl < (*inputClusters).size(); ++icl) {
       std::vector<art::Ptr<recob::Hit>> cluhits = hitsFromCls.at(icl);
       if (cluhits.empty()) continue;
       if (fCompareProductIDs) {
         if (cluhits[0].id() != fHitHandle.id())
           throw cet::exception("ClusterTrackAna")
             << "Hits associated with ClusterModuleLabel are in a different collection than "
                "HitModuleLabel.\n";
         fCompareProductIDs = false;
       } // fCompareProductIDs
       // only load MC-matched hits. Hits that are not MC-matched were either matched to a
       // MCParticle that was not selected or were mis-reconstructed by the hit finder. Neither
       // of these conditions warrant penalizing the reconstruction module performance
       std::vector<unsigned int> hitsIndex;
       for (auto& cluhit : cluhits) {
         if (fHitMCPIndex[cluhit.key()] != UINT_MAX) hitsIndex.push_back(cluhit.key());
       }
       if (hitsIndex.empty()) continue;
       recoIndex.push_back(icl);
       recoHits.push_back(hitsIndex);
     } // icl
   }
   else {
     // get MC-matched hits from tracks
     if (!evt.getByLabel(fTrackModuleLabel, inputTracks))
       throw cet::exception("ClusterTrackAna")
         << "Failed to get a handle to the track collection '" << fTrackModuleLabel.label() << "'\n";
     art::FindManyP<recob::Hit> hitsFromTrk(inputTracks, evt, fTrackModuleLabel);
     if (!hitsFromTrk.isValid())
       throw cet::exception("ClusterTrackAna") << "Failed to get a handle to Track -> Hit assns\n";
     for (unsigned int itk = 0; itk < (*inputTracks).size(); ++itk) {
       std::vector<art::Ptr<recob::Hit>> trkhits = hitsFromTrk.at(itk);
       if (trkhits.empty()) continue;
       if (fCompareProductIDs) {
         if (trkhits[0].id() != fHitHandle.id())
           throw cet::exception("ClusterTrackAna")
             << "Hits associated with TrackModuleLabel are in a different collection than "
                "HitModuleLabel.\n";
         fCompareProductIDs = false;
       } // fCompareProductIDs
       std::vector<unsigned int> hitsIndex;
       for (auto& trkhit : trkhits) {
         if (fHitMCPIndex[trkhit.key()] != UINT_MAX) hitsIndex.push_back(trkhit.key());
       }
       if (hitsIndex.empty()) continue;
       recoIndex.push_back(itk);
       recoHits.push_back(hitsIndex);
     } // itk
   }   // get hits from tracks
 
   for (const auto& tpcid : geom->Iterate<geo::TPCID>()) {
     unsigned int tpc = tpcid.TPC;
     if (hitRange[tpc].first == UINT_MAX) continue;
     // iterate over planes
     for (unsigned short plane = 0; plane < wireReadoutGeom.Nplanes(); ++plane) {
       unsigned int tpcMatHitCnt = 0;
       unsigned int tpcTotHitCnt = 0;
       // create a list of (MCParticle index, matched hit count> pairs
       //  mcParticle  plane
       std::vector<std::pair<unsigned int, float>> mcpCnt;
       // count MCParticle, Cluster/Track hit counts - size matched to mcpCnt
       std::vector<std::vector<std::pair<unsigned int, float>>> mcpClsCnt;
       for (unsigned int iht = hitRange[tpc].first; iht <= hitRange[tpc].second; ++iht) {
         auto& hit = (*fHitHandle)[iht];
         if (hit.WireID().TPC != tpc) continue;
         if (hit.WireID().Plane != plane) continue;
         ++tpcTotHitCnt;
         // ignore hits that were either unmatched to the selected set of PDGCodes
         // or have a not-requested origin
         if (fHitMCPIndex[iht] == UINT_MAX) continue;
         ++tpcMatHitCnt;
         unsigned int mcpi = fHitMCPIndex[iht];
         unsigned short mIndx = 0;
         for (mIndx = 0; mIndx < mcpCnt.size(); ++mIndx)
           if (mcpCnt[mIndx].first == mcpi) break;
         if (mIndx == mcpCnt.size()) {
           // found a MCParticle not in the list yet, so add it
           mcpCnt.push_back(std::make_pair(mcpi, 0));
           mcpClsCnt.resize(mcpCnt.size());
         }
         // increment the number of MC-matched hits
         ++mcpCnt[mIndx].second;
       } // iht
       // ignore TPCs/planes with few MC-matched hits
       if (fPrintLevel > 2) {
         mf::LogVerbatim myprt("ClusterTrackAna");
         myprt << "TPC:" << tpc << " Plane:" << plane << " has " << tpcMatHitCnt << "/"
               << tpcTotHitCnt;
         myprt << " (MC-matched hits) / (all hits)";
       }
       if (tpcMatHitCnt < 3) continue;
       // next iterate over all clusters/tracks and count mc-matched hits that are in this TPC/plane
       for (unsigned int ii = 0; ii < recoHits.size(); ++ii) {
         float nRecoHitsInPlane = 0;
         float nRecoMatHitsInPlane = 0;
         for (auto iht : recoHits[ii]) {
           auto& hit = (*fHitHandle)[iht];
           if (hit.WireID().TPC != tpc) continue;
           if (hit.WireID().Plane != plane) continue;
           ++nRecoHitsInPlane;
           if (fHitMCPIndex[iht] == UINT_MAX) continue;
           ++nRecoMatHitsInPlane;
           unsigned int mcpi = fHitMCPIndex[iht];
           // find the MCParticle index in mcpCnt and use it to count the match entry
           // in mcpClsCnt
           unsigned short mIndx = 0;
           for (mIndx = 0; mIndx < mcpCnt.size(); ++mIndx)
             if (mcpCnt[mIndx].first == mcpi) break;
           if (mIndx == mcpCnt.size()) {
             std::cout << "Logic error: fHitMCPIndex = " << fHitMCPIndex[iht]
                       << " is valid but isn't in the list of MCParticles to use. Please send email "
                          "to baller@fnal.gov.\n";
             continue;
           }
           unsigned short cIndx = 0;
           for (cIndx = 0; cIndx < mcpClsCnt[mIndx].size(); ++cIndx)
             if (mcpClsCnt[mIndx][cIndx].first == ii) break;
           if (cIndx == mcpClsCnt[mIndx].size()) mcpClsCnt[mIndx].push_back(std::make_pair(ii, 0));
           ++mcpClsCnt[mIndx][cIndx].second;
         } // cluhit
         if (nRecoMatHitsInPlane == 0) continue;
       } // ii
       // find the cluster that has the most hits matched to each MC Particle
       for (unsigned short mIndx = 0; mIndx < mcpCnt.size(); ++mIndx) {
         // require at least 3 MC-matched hits
         if (mcpCnt[mIndx].second < 3) continue;
         unsigned int mcpi = mcpCnt[mIndx].first;
         auto& mcp = (*mcpHandle)[mcpi];
         unsigned int pdgCode = abs(mcp.PdgCode());
         unsigned short pIndx = USHRT_MAX;
         if (pdgCode == 11) pIndx = 0;
         if (pdgCode == 13) pIndx = 1;
         if (pdgCode == 211) pIndx = 2;
         if (pdgCode == 321) pIndx = 3;
         if (pdgCode == 2212) pIndx = 4;
         if (mcpClsCnt[mIndx].empty()) {
           // Un-reconstructed MCParticle hits
           ++Cnts[pIndx];
           ++fNBadEP;
           if (fPrintLevel > 0) {
             mf::LogVerbatim myprt("ClusterTrackAna");
             myprt << " MCPI " << mcpi;
             int TMeV = 1000 * (mcp.E() - mcp.Mass());
             myprt << " " << TMeV << " MeV";
             std::string partName = std::to_string(pdgCode);
             if (pdgCode == 11) partName = "El";
             if (pdgCode == 13) partName = "Mu";
             if (pdgCode == 211) partName = "Pi";
             if (pdgCode == 311) partName = "Ka";
             if (pdgCode == 2212) partName = "Pr";
             myprt << std::setw(5) << partName;
             // print out the range of truth-matched hits
             unsigned int firstHitIndex = UINT_MAX;
             unsigned int lastHitIndex = UINT_MAX;
             FirstLastHitInPlane(tpc, plane, mcpi, firstHitIndex, lastHitIndex);
             myprt << " Failed to reconstruct. Truth-matched hit range from ";
             myprt << HitLocation(firstHitIndex);
             myprt << " to ";
             myprt << HitLocation(lastHitIndex);
             myprt << " <- EP = 0!";
           } // fPrintLevel > 0
           continue;
         } // (mcpClsCnt[mIndx].empty()
         std::pair<unsigned int, float> big = std::make_pair(UINT_MAX, 0);
         for (unsigned short cIndx = 0; cIndx < mcpClsCnt[mIndx].size(); ++cIndx) {
           auto& mcc = mcpClsCnt[mIndx][cIndx];
           if (mcc.second > big.second) big = mcc;
         } // cIndx
         if (big.first == UINT_MAX) {
           if (fPrintLevel > 2) {
             mf::LogVerbatim myprt("ClusterTrackAna");
             unsigned int mcpi = mcpCnt[mIndx].first;
             auto& mcp = (*mcpHandle)[mcpi];
             myprt << "Match failed: mcpi " << mcpi << " pdg " << mcp.PdgCode();
           }
           std::cout << "match failed mIndx " << mIndx << "\n";
           continue;
         } // big.first == UINT_MAX
         unsigned int ii = big.first;
         float eff = big.second / mcpCnt[mIndx].second;
         float nRecoHitsInPlane = 0;
         // define some variables to print the range of the cluster/track
         unsigned int firstRecoHitIndex = UINT_MAX;
         unsigned int lastRecoHitIndex = UINT_MAX;
         for (auto iht : recoHits[ii]) {
           auto& hit = (*fHitHandle)[iht];
           if (hit.WireID().TPC != tpc) continue;
           if (hit.WireID().Plane != plane) continue;
           ++nRecoHitsInPlane;
           if (firstRecoHitIndex == UINT_MAX) {
             firstRecoHitIndex = iht;
             lastRecoHitIndex = iht;
           }
           unsigned int wire = (*fHitHandle)[iht].WireID().Wire;
           if (wire < (*fHitHandle)[firstRecoHitIndex].WireID().Wire) firstRecoHitIndex = iht;
           if (wire > (*fHitHandle)[lastRecoHitIndex].WireID().Wire) lastRecoHitIndex = iht;
         } // iht
         float pur = big.second / nRecoHitsInPlane;
         ++Cnts[pIndx];
         float ep = eff * pur;
         EPSums[pIndx] += ep;
         ESums[pIndx] += eff;
         PSums[pIndx] += pur;
         bool hasBadEP = (ep < fBadEP);
         if (hasBadEP) ++fNBadEP;
         bool prt = fPrintLevel > 1 || (fPrintLevel == 1 && hasBadEP);
         if (prt) {
           mf::LogVerbatim myprt("ClusterTrackAna");
           if (fFirstPrint) {
             myprt << "Hit location format is TPC:Plane:Wire:Tick\n";
             myprt << " MCPI     Ptcl T(MeV)  nMCPHits  ____From____ _____To_____";
             if (inputClusters.isValid()) { myprt << "  clsID"; }
             else {
               myprt << "  trkID";
             }
             myprt
               << "  ____From____  _____To_____   nRecoHits nMCPRecoHits    Eff      Pur      EP\n";
             fFirstPrint = false;
           } // fFirstPrint
           myprt << std::setw(5) << mcpi;
           // convert the PDG code into nicer format
           std::string partName = std::to_string(pdgCode);
           if (pdgCode == 11) partName = " Electron";
           if (pdgCode == 13) partName = "     Muon";
           if (pdgCode == 211) partName = "     Pion";
           if (pdgCode == 311) partName = "     Kaon";
           if (pdgCode == 2212) partName = "   Proton";
           myprt << partName;
           int TMeV = 1000 * (mcp.E() - mcp.Mass());
           myprt << std::setw(7) << TMeV;
           // nMCPHits
           myprt << std::setw(10) << mcpCnt[mIndx].second;
           unsigned int firstTruHitIndex = UINT_MAX;
           unsigned int lastTruHitIndex = UINT_MAX;
           FirstLastHitInPlane(tpc, plane, mcpi, firstTruHitIndex, lastTruHitIndex);
           myprt << std::setw(14) << HitLocation(firstTruHitIndex);
           myprt << std::setw(14) << HitLocation(lastTruHitIndex);
           int id = -1;
           if (inputClusters.isValid()) {
             // print cluster info
             auto& cls = (*inputClusters)[recoIndex[ii]];
             id = cls.ID();
           }
           else if (inputTracks.isValid()) {
             // print track info
             auto& trk = (*inputTracks)[recoIndex[ii]];
             id = trk.ID();
           }
           myprt << std::setw(6) << id;
           myprt << std::setw(14) << HitLocation(firstRecoHitIndex);
           myprt << std::setw(14) << HitLocation(lastRecoHitIndex);
           myprt << std::setw(12) << (int)nRecoHitsInPlane;
           myprt << std::setw(13) << (int)big.second;
           myprt << std::fixed << std::setprecision(2);
           myprt << std::setw(8) << eff;
           myprt << std::setw(8) << pur;
           myprt << std::setw(8) << ep;
           if (hasBadEP) myprt << " <- BadEP";
           myprt << " Evt: " << evt.event();
           myprt << " Evt Cnt " << fEventCnt;
         } // prt
       }   // mIndx
     }     // plane
   }       // tpcid
 
   fHitMCPIndex.resize(0);
 
 } // analyze

template<typename T , BranchType BT>

ProductToken< T > art::ModuleBase::consumes ( InputTag const & tag )

protectedinherited

Definition at line 61 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumes().

   {
     return collector_.consumes<T, BT>(tag);
   }

ConsumesCollector & art::ModuleBase::consumesCollector ( )

protectedinherited

Definition at line 57 of file ModuleBase.cc.

References art::ModuleBase::collector_.

   {
     return collector_;
   }

template<typename T , BranchType BT>

void art::ModuleBase::consumesMany ( )

protectedinherited

Definition at line 75 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesMany().

   {
     collector_.consumesMany<T, BT>();
   }

template<typename Element , BranchType = InEvent>

ViewToken<Element> art::ModuleBase::consumesView ( InputTag const & )

protectedinherited

template<typename T , BranchType BT>

ViewToken<T> art::ModuleBase::consumesView ( InputTag const & tag )

inherited

Definition at line 68 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::consumesView().

   {
     return collector_.consumesView<T, BT>(tag);
   }

void art::detail::Analyzer::doBeginJob ( SharedResources const & resources )

inherited

Definition at line 25 of file Analyzer.cc.

Referenced by art::detail::Analyzer::Analyzer().

   {
     setupQueues(resources);
     ProcessingFrame const frame{ScheduleID{}};
     beginJobWithFrame(frame);
   }

bool art::detail::Analyzer::doBeginRun	(	RunPrincipal &	rp,
		ModuleContext const &	mc
	)

inherited

Definition at line 68 of file Analyzer.cc.

References art::ModuleContext::scheduleID().

Referenced by art::detail::Analyzer::Analyzer().

   {
     ProcessingFrame const frame{mc.scheduleID()};
     beginRunWithFrame(std::as_const(rp).makeRun(mc), frame);
     return true;
   }

bool art::detail::Analyzer::doBeginSubRun	(	SubRunPrincipal &	srp,
		ModuleContext const &	mc
	)

inherited

Definition at line 84 of file Analyzer.cc.

References art::ModuleContext::scheduleID().

Referenced by art::detail::Analyzer::Analyzer().

   {
     ProcessingFrame const frame{mc.scheduleID()};
     beginSubRunWithFrame(std::as_const(srp).makeSubRun(mc), frame);
     return true;
   }

void art::detail::Analyzer::doEndJob ( )

inherited

Definition at line 33 of file Analyzer.cc.

Referenced by art::detail::Analyzer::Analyzer().

   {
     ProcessingFrame const frame{ScheduleID{}};
     endJobWithFrame(frame);
   }

bool art::detail::Analyzer::doEndRun	(	RunPrincipal &	rp,
		ModuleContext const &	mc
	)

inherited

Definition at line 76 of file Analyzer.cc.

References art::ModuleContext::scheduleID().

Referenced by art::detail::Analyzer::Analyzer().

   {
     ProcessingFrame const frame{mc.scheduleID()};
     endRunWithFrame(std::as_const(rp).makeRun(mc), frame);
     return true;
   }

bool art::detail::Analyzer::doEndSubRun	(	SubRunPrincipal &	srp,
		ModuleContext const &	mc
	)

inherited

Definition at line 92 of file Analyzer.cc.

References art::ModuleContext::scheduleID().

Referenced by art::detail::Analyzer::Analyzer().

   {
     ProcessingFrame const frame{mc.scheduleID()};
     endSubRunWithFrame(std::as_const(srp).makeSubRun(mc), frame);
     return true;
   }

bool art::detail::Analyzer::doEvent	(	EventPrincipal &	ep,
		ModuleContext const &	mc,
		std::atomic< std::size_t > &	counts_run,
		std::atomic< std::size_t > &	counts_passed,
		std::atomic< std::size_t > &	counts_failed
	)

inherited

Definition at line 100 of file Analyzer.cc.

References e, and art::ModuleContext::scheduleID().

Referenced by art::detail::Analyzer::Analyzer().

   {
     auto const e = std::as_const(ep).makeEvent(mc);
     if (wantEvent(mc.scheduleID(), e)) {
       ++counts_run;
       ProcessingFrame const frame{mc.scheduleID()};
       analyzeWithFrame(e, frame);
       ++counts_passed;
     }
     return true;
   }

void art::detail::Analyzer::doRespondToCloseInputFile ( FileBlock const & fb )

inherited

Definition at line 47 of file Analyzer.cc.

Referenced by art::detail::Analyzer::Analyzer().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToCloseInputFileWithFrame(fb, frame);
   }

void art::detail::Analyzer::doRespondToCloseOutputFiles ( FileBlock const & fb )

inherited

Definition at line 61 of file Analyzer.cc.

Referenced by art::detail::Analyzer::Analyzer().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToCloseOutputFilesWithFrame(fb, frame);
   }

void art::detail::Analyzer::doRespondToOpenInputFile ( FileBlock const & fb )

inherited

Definition at line 40 of file Analyzer.cc.

Referenced by art::detail::Analyzer::Analyzer().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToOpenInputFileWithFrame(fb, frame);
   }

void art::detail::Analyzer::doRespondToOpenOutputFiles ( FileBlock const & fb )

inherited

Definition at line 54 of file Analyzer.cc.

Referenced by art::detail::Analyzer::Analyzer().

   {
     ProcessingFrame const frame{ScheduleID{}};
     respondToOpenOutputFilesWithFrame(fb, frame);
   }

void cluster::ClusterTrackAna::endJob ( )

overrideprivatevirtual

Reimplemented from art::EDAnalyzer.

Definition at line 600 of file ClusterTrackAna_module.cc.

References Cnts, DEFINE_ART_MODULE, EPSums, ESums, fBadEP, fClusterModuleLabel, fEventCnt, fInTPC, fNBadEP, fTrackModuleLabel, fTruthOrigin, art::InputTag::label(), and PSums.

 {
   // output results
   mf::LogVerbatim myprt("ClusterTrackAna");
   myprt << "ClusterTrackAna summary results for " << fEventCnt;
   if (fClusterModuleLabel.label() != "NA") {
     myprt << " events using ClusterModuleLabel: " << fClusterModuleLabel.label();
   }
   else {
     myprt << " events using TrackModuleLabel: " << fTrackModuleLabel.label();
   }
   myprt << " Origin: " << fTruthOrigin;
   if (fInTPC != UINT_MAX) { myprt << " in TPC " << fInTPC; }
   else {
     myprt << " in all TPCs";
   }
   myprt << "\n";
   float cnts = 0;
   for (unsigned short pIndx = 0; pIndx < 5; ++pIndx)
     cnts += Cnts[pIndx];
   if (cnts == 0) {
     myprt << "No ClusterTrackAna results";
     return;
   }
   float sumEP = 0;
   float sumE = 0;
   float sumP = 0;
   myprt << "Efficiency (Eff), Purity (Pur) and Eff * Pur (EP) by selected truth particle types\n";
   std::array<std::string, 5> pName = {{"El", "Mu", "Pi", "K ", "P "}};
   myprt << "particle     Eff     Pur     EP\n";
   myprt << "-------------------------------";
   for (unsigned short pIndx = 0; pIndx < 5; ++pIndx) {
     if (Cnts[pIndx] == 0) continue;
     float ave;
     myprt << "\n   " << pName[pIndx] << "   ";
     myprt << std::fixed << std::setprecision(3);
     ave = ESums[pIndx] / Cnts[pIndx];
     myprt << std::setw(8) << ave;
     ave = PSums[pIndx] / Cnts[pIndx];
     myprt << std::setw(8) << ave;
     ave = EPSums[pIndx] / Cnts[pIndx];
     myprt << std::setw(8) << ave;
     if (pIndx == 0) continue;
     sumEP += EPSums[pIndx];
     sumE += ESums[pIndx];
     sumP += PSums[pIndx];
   } // pIndx
   if (cnts == 0) return;
   myprt << "\n";
   myprt << "Averages for all selected truth particles\n";
   myprt << "Ave Eff " << sumE / cnts;
   myprt << " Ave Pur " << sumP / cnts;
   myprt << " Ave EP " << sumEP / cnts;
   myprt << " nBadEP " << fNBadEP;
   myprt << " (EP < " << std::fixed << std::setprecision(2) << fBadEP << ")";
   myprt << "\n";
   myprt << "MCParticle counts in all TPCs and Planes:";
   for (unsigned short pIndx = 0; pIndx < 5; ++pIndx) {
     if (Cnts[pIndx] == 0) continue;
     myprt << " " << pName[pIndx] << " " << (int)Cnts[pIndx];
   } // pIndx
 } // endJob

void cluster::ClusterTrackAna::FirstLastHitInPlane	(	unsigned int	tpc,
		unsigned int	plane,
		unsigned int	mcpi,
		unsigned int &	firstHitIndex,
		unsigned int &	lastHitIndex
	)

private

Definition at line 574 of file ClusterTrackAna_module.cc.

References fHitHandle, and fHitMCPIndex.

Referenced by analyze().

 {
   // Returns the index of the first hit (lowest wire number) and last hit (highest wire number)
   // matched to the MCParticle indexed by mcpi in the requested tpc, plane
   firstHitIndex = UINT_MAX;
   lastHitIndex = UINT_MAX;
   for (unsigned int iht = 0; iht < (*fHitHandle).size(); ++iht) {
     if (fHitMCPIndex[iht] != mcpi) continue;
     auto& hit = (*fHitHandle)[iht];
     if (hit.WireID().TPC != tpc) continue;
     if (hit.WireID().Plane != plane) continue;
     if (firstHitIndex == UINT_MAX) {
       firstHitIndex = iht;
       lastHitIndex = iht;
     }
     unsigned int wire = (*fHitHandle)[iht].WireID().Wire;
     if (wire < (*fHitHandle)[firstHitIndex].WireID().Wire) firstHitIndex = iht;
     if (wire > (*fHitHandle)[lastHitIndex].WireID().Wire) lastHitIndex = iht;
   } // iht
 } // FirstLastHitInPlane

std::array< std::vector< ProductInfo >, NumBranchTypes > const & art::ModuleBase::getConsumables ( ) const

inherited

Definition at line 43 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::getConsumables().

   {
     return collector_.getConsumables();
   }

Handle< TriggerResults > art::Observer::getTriggerResults ( Event const & e ) const

protectedinherited

Definition at line 75 of file Observer.cc.

References art::ProductRetriever::get(), and art::Observer::selectors_.

Referenced by art::OutputModule::doWriteEvent(), and art::Observer::wantAllEvents().

   {
     if (selectors_) {
       return selectors_->getOneTriggerResults(e);
     }
 
     // The following applies for cases where no SelectEvents entries
     // exist.
     Handle<TriggerResults> h;
     if (e.get(empty_process_name, h)) {
       return h;
     }
     return Handle<TriggerResults>{};
   }

std::string cluster::ClusterTrackAna::HitLocation ( unsigned int iht )

private

packs hit WireID and PeakTime into a compact format

Definition at line 564 of file ClusterTrackAna_module.cc.

References util::to_string().

Referenced by analyze().

 {
   // Put the hit location into a compact human-readable format
   if (iht >= (*fHitHandle).size()) return "NA";
   auto& hit = (*fHitHandle)[iht];
   return std::to_string(hit.WireID().TPC) + ":" + std::to_string(hit.WireID().Plane) + ":" +
          std::to_string(hit.WireID().Wire) + ":" + std::to_string((int)hit.PeakTime());
 } // HitLocation

std::unique_ptr< Worker > art::ModuleBase::makeWorker ( WorkerParams const & wp )

inherited

Definition at line 37 of file ModuleBase.cc.

References art::ModuleBase::doMakeWorker(), and art::NumBranchTypes.

   {
     return doMakeWorker(wp);
   }

template<typename T , BranchType BT>

ProductToken< T > art::ModuleBase::mayConsume ( InputTag const & tag )

protectedinherited

Definition at line 82 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsume().

   {
     return collector_.mayConsume<T, BT>(tag);
   }

template<typename T , BranchType BT>

void art::ModuleBase::mayConsumeMany ( )

protectedinherited

Definition at line 96 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeMany().

   {
     collector_.mayConsumeMany<T, BT>();
   }

template<typename Element , BranchType = InEvent>

ViewToken<Element> art::ModuleBase::mayConsumeView ( InputTag const & )

protectedinherited

template<typename T , BranchType BT>

ViewToken<T> art::ModuleBase::mayConsumeView ( InputTag const & tag )

inherited

Definition at line 89 of file ModuleBase.h.

References art::ModuleBase::collector_, and art::ConsumesCollector::mayConsumeView().

   {
     return collector_.mayConsumeView<T, BT>(tag);
   }

ModuleDescription const & art::ModuleBase::moduleDescription ( ) const

inherited

Definition at line 13 of file ModuleBase.cc.

References art::errors::LogicError.

Referenced by art::OutputModule::doRespondToOpenInputFile(), art::OutputModule::doWriteEvent(), art::Modifier::fillProductDescriptions(), art::OutputModule::makePlugins_(), art::OutputWorker::OutputWorker(), reco::shower::LArPandoraModularShowerCreation::produce(), art::Modifier::registerProducts(), and art::OutputModule::registerProducts().

   {
     if (md_.has_value()) {
       return *md_;
     }
 
     throw Exception{errors::LogicError,
                     "There was an error while calling moduleDescription().\n"}
       << "The moduleDescription() base-class member function cannot be called\n"
          "during module construction.  To determine which module is "
          "responsible\n"
          "for calling it, find the '<module type>:<module "
          "label>@Construction'\n"
          "tag in the message prefix above.  Please contact artists@fnal.gov\n"
          "for guidance.\n";
   }

ClusterTrackAna& cluster::ClusterTrackAna::operator= ( ClusterTrackAna const & )

delete

ClusterTrackAna& cluster::ClusterTrackAna::operator= ( ClusterTrackAna && )

delete

string const & art::Observer::processName ( ) const

protectedinherited

Definition at line 57 of file Observer.cc.

References art::Observer::process_name_.

Referenced by art::FileDumperOutput::printPrincipal().

   {
     return process_name_;
   }

void art::ModuleBase::setModuleDescription ( ModuleDescription const & md )

inherited

Definition at line 31 of file ModuleBase.cc.

References art::ModuleBase::md_.

   {
     md_ = md;
   }

void art::ModuleBase::sortConsumables ( std::string const & current_process_name )

inherited

Definition at line 49 of file ModuleBase.cc.

References art::ModuleBase::collector_, and art::ConsumesCollector::sortConsumables().

   {
     // Now that we know we have seen all the consumes declarations,
     // sort the results for fast lookup later.
     collector_.sortConsumables(current_process_name);
   }

bool art::Observer::wantAllEvents ( ) const

inlineprotectednoexceptinherited

Definition at line 31 of file Observer.h.

References e, art::Observer::getTriggerResults(), art::Observer::wantAllEvents_, and art::Observer::wantEvent().

     {
       return wantAllEvents_;
     }

bool art::Observer::wantEvent	(	ScheduleID	id,
		Event const &	e
	)		const

protectedinherited

Definition at line 63 of file Observer.cc.

References art::Observer::rejectors_, art::Observer::selectors_, and art::Observer::wantAllEvents_.

Referenced by art::OutputModule::doEvent(), art::OutputModule::doWriteEvent(), and art::Observer::wantAllEvents().

   {
     if (wantAllEvents_) {
       return true;
     }
     bool const select_event = selectors_ ? selectors_->matchEvent(id, e) : true;
     bool const reject_event =
       rejectors_ ? rejectors_->matchEvent(id, e) : false;
     return select_event and not reject_event;
   }