562 std::cout <<
" RecoNeutrinos: " << recoNeutrinoVector.size() << std::endl;
563 std::cout <<
" RecoParticles: " << recoParticleVector.size() << std::endl;
589 if (trueHitsToParticles.empty()) {
591 throw cet::exception(
"LArPandora") <<
" PFParticleMonitoring::analyze - no sim channels "
592 "found, backtracker module must be set in FHiCL "
609 std::cout <<
" TrueParticles: " << particlesToTruth.size() << std::endl;
610 std::cout <<
" TrueEvents: " << truthToParticles.size() << std::endl;
611 std::cout <<
" MatchedParticles: " << trueParticlesToHits.size() << std::endl;
614 if (trueParticlesToHits.empty()) {
633 for (PFParticleVector::const_iterator iter = recoParticleVector.begin(),
634 iterEnd = recoParticleVector.end();
637 const art::Ptr<recob::PFParticle> recoParticle = *iter;
655 recoParticleMap, recoParticlesToHits, recoNeutrinosToHits, recoHitsToNeutrinos);
657 truthToParticles, trueParticlesToHits, trueNeutrinosToHits, trueHitsToNeutrinos);
662 recoNeutrinosToHits, trueHitsToNeutrinos, matchedNeutrinos, matchedNeutrinoHits);
664 for (MCTruthToHits::const_iterator iter = trueNeutrinosToHits.begin(),
665 iterEnd = trueNeutrinosToHits.end();
668 const art::Ptr<simb::MCTruth> trueEvent = iter->first;
669 const HitVector& trueHitVector = iter->second;
671 if (trueHitVector.empty())
continue;
673 if (!trueEvent->NeutrinoSet())
continue;
675 const simb::MCNeutrino trueNeutrino(trueEvent->GetNeutrino());
676 const simb::MCParticle trueParticle(trueNeutrino.Nu());
679 m_mcPdg = trueParticle.PdgCode();
694 m_mcDirX = trueParticle.Px() / trueParticle.P();
695 m_mcDirY = trueParticle.Py() / trueParticle.P();
696 m_mcDirZ = trueParticle.Pz() / trueParticle.P();
746 for (HitVector::const_iterator hIter1 = trueHitVector.begin(),
747 hIterEnd1 = trueHitVector.end();
750 if (recoHitsToNeutrinos.find(*hIter1) == recoHitsToNeutrinos.end())
754 MCTruthToPFParticles::const_iterator pIter1 = matchedNeutrinos.find(trueEvent);
755 if (matchedNeutrinos.end() != pIter1) {
756 const art::Ptr<recob::PFParticle> recoParticle = pIter1->second;
765 std::cout <<
" Warning: Found neutrino with an invalid PDG code " << std::endl;
767 PFParticlesToHits::const_iterator pIter2 = recoNeutrinosToHits.find(recoParticle);
768 if (recoParticlesToHits.end() != pIter2) {
769 const HitVector& recoHitVector = pIter2->second;
771 for (HitVector::const_iterator hIter2 = recoHitVector.begin(),
772 hIterEnd2 = recoHitVector.end();
775 if (trueHitsToNeutrinos.find(*hIter2) == trueHitsToNeutrinos.end())
779 MCTruthToHits::const_iterator pIter3 = matchedNeutrinoHits.find(trueEvent);
780 if (matchedNeutrinoHits.end() != pIter3) {
781 const HitVector& matchedHitVector = pIter3->second;
793 PFParticlesToVertices::const_iterator pIter4 =
794 recoParticlesToVertices.find(recoParticle);
795 if (recoParticlesToVertices.end() != pIter4) {
797 if (!vertexVector.empty()) {
799 std::cout <<
" Warning: Found particle with more than one associated vertex "
802 const art::Ptr<recob::Vertex> recoVertex = *(vertexVector.begin());
803 double xyz[3] = {0.0, 0.0, 0.0};
804 recoVertex->XYZ(xyz);
840 recoParticlesToHits, trueHitsToParticles, matchedParticles, matchedParticleHits);
843 for (MCParticlesToHits::const_iterator iter = trueParticlesToHits.begin(),
844 iterEnd = trueParticlesToHits.end();
847 const art::Ptr<simb::MCParticle> trueParticle = iter->first;
848 const HitVector& trueHitVector = iter->second;
850 if (trueHitVector.empty())
continue;
852 m_mcPdg = trueParticle->PdgCode();
928 m_mcVtxX = trueParticle->Vx(startT);
929 m_mcVtxY = trueParticle->Vy(startT);
930 m_mcVtxZ = trueParticle->Vz(startT);
943 const double Ptot(trueParticle->P(startT));
946 m_mcDirX = trueParticle->Px(startT) / Ptot;
947 m_mcDirY = trueParticle->Py(startT) / Ptot;
948 m_mcDirZ = trueParticle->Pz(startT) / Ptot;
952 catch (cet::exception&
e) {
956 MCParticlesToMCTruth::const_iterator nuIter = particlesToTruth.find(trueParticle);
957 if (particlesToTruth.end() == nuIter)
958 throw cet::exception(
"LArPandora") <<
" PFParticleMonitoring::analyze --- Found a true "
959 "particle without any ancestry information ";
961 const art::Ptr<simb::MCTruth> trueEvent = nuIter->second;
963 if (trueEvent->NeutrinoSet()) {
964 const simb::MCNeutrino neutrino(trueEvent->GetNeutrino());
966 m_mcIsCC = ((
simb::kCC == neutrino.CCNC()) ? 1 : 0);
971 const art::Ptr<simb::MCParticle> parentParticle(
973 const art::Ptr<simb::MCParticle> primaryParticle(
975 m_mcParentPdg = ((parentParticle != trueParticle) ? parentParticle->PdgCode() : 0);
978 m_mcIsDecay = (
"Decay" == trueParticle->Process());
980 catch (cet::exception& e) {
984 bool foundSpacePoints(
false);
986 for (HitVector::const_iterator hIter1 = trueHitVector.begin(),
987 hIterEnd1 = trueHitVector.end();
990 const art::Ptr<recob::Hit>
hit = *hIter1;
992 HitsToSpacePoints::const_iterator hIter2 = hitsToSpacePoints.find(hit);
993 if (hitsToSpacePoints.end() == hIter2)
continue;
995 const art::Ptr<recob::SpacePoint> spacepoint = hIter2->second;
996 const double X(spacepoint->XYZ()[0]);
998 if (!foundSpacePoints) {
1001 foundSpacePoints =
true;
1010 for (HitVector::const_iterator hIter1 = trueHitVector.begin(),
1011 hIterEnd1 = trueHitVector.end();
1012 hIter1 != hIterEnd1;
1014 if (recoHitsToParticles.find(*hIter1) == recoHitsToParticles.end())
1024 MCParticlesToPFParticles::const_iterator pIter1 = matchedParticles.find(trueParticle);
1025 if (matchedParticles.end() != pIter1) {
1026 const art::Ptr<recob::PFParticle> recoParticle = pIter1->second;
1027 m_pfoPdg = recoParticle->PdgCode();
1031 const art::Ptr<recob::PFParticle> parentParticle =
1035 const art::Ptr<recob::PFParticle> primaryParticle =
1039 PFParticlesToHits::const_iterator pIter2 = recoParticlesToHits.find(recoParticle);
1040 if (recoParticlesToHits.end() == pIter2)
1041 throw cet::exception(
"LArPandora")
1042 <<
" PFParticleMonitoring::analyze --- Found a reco particle without any hits ";
1044 const HitVector& recoHitVector = pIter2->second;
1046 for (HitVector::const_iterator hIter2 = recoHitVector.begin(),
1047 hIterEnd2 = recoHitVector.end();
1048 hIter2 != hIterEnd2;
1050 if (trueHitsToParticles.find(*hIter2) == trueHitsToParticles.end())
1054 MCParticlesToHits::const_iterator pIter3 = matchedParticleHits.find(trueParticle);
1055 if (matchedParticleHits.end() == pIter3)
1056 throw cet::exception(
"LArPandora") <<
" PFParticleMonitoring::analyze --- Found a "
1057 "matched true particle without matched hits ";
1059 const HitVector& matchedHitVector = pIter3->second;
1071 PFParticlesToVertices::const_iterator pIter4 = recoParticlesToVertices.find(recoParticle);
1072 if (recoParticlesToVertices.end() != pIter4) {
1074 if (!vertexVector.empty()) {
1076 std::cout <<
" Warning: Found particle with more than one associated vertex "
1079 const art::Ptr<recob::Vertex> recoVertex = *(vertexVector.begin());
1080 double xyz[3] = {0.0, 0.0, 0.0};
1081 recoVertex->XYZ(xyz);
1090 PFParticlesToTracks::const_iterator pIter5 = recoParticlesToTracks.find(recoParticle);
1091 if (recoParticlesToTracks.end() != pIter5) {
1093 if (!trackVector.empty()) {
1095 std::cout <<
" Warning: Found particle with more than one associated track "
1098 const art::Ptr<recob::Track> recoTrack = *(trackVector.begin());
1099 const auto& vtxPosition = recoTrack->Vertex();
1100 const auto& endPosition = recoTrack->End();
1101 const auto& vtxDirection = recoTrack->VertexDirection();
1118 m_pfoIsStitched = (particlesToT0s.end() != particlesToT0s.find(recoParticle));
1130 <<
" trueNu=" << m_mcNuPdg <<
", truePdg=" <<
m_mcPdg
static void BuildPFParticleHitMaps(const PFParticleVector &particleVector, const PFParticlesToSpacePoints &particlesToSpacePoints, const SpacePointsToHits &spacePointsToHits, PFParticlesToHits &particlesToHits, HitsToPFParticles &hitsToParticles, const DaughterMode daughterMode=kUseDaughters)
Build mapping between PFParticles and Hits using PFParticle/SpacePoint/Hit maps.
std::map< art::Ptr< recob::Hit >, art::Ptr< simb::MCParticle > > HitsToMCParticles
int CountHitsByType(const int view, const HitVector &hitVector) const
Count the number of reconstructed hits in a given wire plane.
static void BuildPFParticleMap(const PFParticleVector &particleVector, PFParticleMap &particleMap)
Build particle maps for reconstructed particles.
static int GetParentNeutrino(const PFParticleMap &particleMap, const art::Ptr< recob::PFParticle > daughterParticle)
Return the parent neutrino PDG code (or zero for cosmics) for a given reconstructed particle...
bool m_useDaughterPFParticles
void GetRecoToTrueMatches(const PFParticlesToHits &recoNeutrinosToHits, const HitsToMCTruth &trueHitsToNeutrinos, MCTruthToPFParticles &matchedNeutrinos, MCTruthToHits &matchedNeutrinoHits) const
Perform matching between true and reconstructed neutrino events.
static bool IsNeutrino(const art::Ptr< recob::PFParticle > particle)
Determine whether a particle has been reconstructed as a neutrino.
std::string m_backtrackerLabel
std::map< art::Ptr< simb::MCTruth >, HitVector > MCTruthToHits
std::map< art::Ptr< simb::MCParticle >, art::Ptr< simb::MCTruth > > MCParticlesToMCTruth
std::string m_geantModuleLabel
static art::Ptr< recob::PFParticle > GetFinalStatePFParticle(const PFParticleMap &particleMap, const art::Ptr< recob::PFParticle > daughterParticle)
Return the final-state parent particle by navigating up the chain of parent/daughter associations...
std::map< art::Ptr< simb::MCParticle >, art::Ptr< recob::PFParticle > > MCParticlesToPFParticles
bool m_addDaughterMCParticles
then echo echo For and will not be changed by echo further linking echo echo B echo The symbol is in the uninitialized data multiple common symbols may appear with the echo same name If the symbol is defined the common echo symbols are treated as undefined references For more echo details on common see the discussion of warn common echo in *Note Linker see the discussion of warn common echo in *Note Linker such as a global int variable echo as opposed to a large global array echo echo I echo The symbol is an indirect reference to another symbol This echo is a GNU extension to the a out object file format which is echo rarely used echo echo N echo The symbol is a debugging symbol echo echo R echo The symbol is in a read only data section echo echo S echo The symbol is in an uninitialized data section for small echo objects echo echo T echo The symbol is in the the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo W echo The symbol is a weak symbol that has not been specifically echo tagged as a weak object symbol When a weak defined symbol echo is linked with a normal defined the normal defined echo symbol is used with no error When a weak undefined symbol echo is linked and the symbol is not the value of the echo weak symbol becomes zero with no error echo echo echo The symbol is a stabs symbol in an a out object file In echo this the next values printed are the stabs other echo the stabs desc and the stab type Stabs symbols are echo used to hold debugging information For more echo see *Note or object file format specific echo echo For Mac OS X
static art::Ptr< recob::PFParticle > GetParentPFParticle(const PFParticleMap &particleMap, const art::Ptr< recob::PFParticle > daughterParticle)
Return the top-level parent particle by navigating up the chain of parent/daughter associations...
std::map< art::Ptr< recob::PFParticle >, TrackVector > PFParticlesToTracks
std::map< art::Ptr< recob::PFParticle >, VertexVector > PFParticlesToVertices
std::map< int, art::Ptr< recob::PFParticle > > PFParticleMap
void BuildTrueNeutrinoHitMaps(const MCTruthToMCParticles &truthToParticles, const MCParticlesToHits &trueParticlesToHits, MCTruthToHits &trueNeutrinosToHits, HitsToMCTruth &trueHitsToNeutrinos) const
Build mapping from true neutrinos to hits.
std::map< int, art::Ptr< simb::MCParticle > > MCParticleMap
bool m_useDaughterMCParticles
std::string m_hitfinderLabel
static void CollectSpacePoints(const art::Event &evt, const std::string &label, SpacePointVector &spacePointVector, SpacePointsToHits &spacePointsToHits)
Collect the reconstructed SpacePoints and associated hits from the ART event record.
Charged-current interactions.
double GetLength(const art::Ptr< simb::MCParticle > trueParticle, const int startT, const int endT) const
Find the length of the true particle trajectory through the active region of the detector.
int m_nTrueWithoutRecoHits
True hits which don't belong to any reconstructed particle - "available".
std::map< art::Ptr< recob::PFParticle >, T0Vector > PFParticlesToT0s
std::map< art::Ptr< recob::Hit >, art::Ptr< simb::MCTruth > > HitsToMCTruth
static art::Ptr< simb::MCParticle > GetParentMCParticle(const MCParticleMap &particleMap, const art::Ptr< simb::MCParticle > daughterParticle)
Return the top-level parent particle by navigating up the chain of parent/daughter associations...
void BuildRecoNeutrinoHitMaps(const PFParticleMap &recoParticleMap, const PFParticlesToHits &recoParticlesToHits, PFParticlesToHits &recoNeutrinosToHits, HitsToPFParticles &recoHitsToNeutrinos) const
Build mapping from reconstructed neutrinos to hits.
std::map< art::Ptr< simb::MCParticle >, HitVector > MCParticlesToHits
std::vector< art::Ptr< recob::PFParticle > > PFParticleVector
std::map< art::Ptr< simb::MCTruth >, MCParticleVector > MCTruthToMCParticles
std::map< art::Ptr< simb::MCTruth >, art::Ptr< recob::PFParticle > > MCTruthToPFParticles
std::vector< art::Ptr< recob::SpacePoint > > SpacePointVector
std::map< art::Ptr< recob::PFParticle >, HitVector > PFParticlesToHits
bool m_printDebug
switch for print statements (TODO: use message service!)
static void CollectVertices(const art::Event &evt, const std::string &label, VertexVector &vertexVector, PFParticlesToVertices &particlesToVertices)
Collect the reconstructed PFParticles and associated Vertices from the ART event record.
std::map< art::Ptr< recob::Hit >, art::Ptr< recob::SpacePoint > > HitsToSpacePoints
std::vector< art::Ptr< simb::MCParticle > > MCParticleVector
std::vector< art::Ptr< recob::Track > > TrackVector
static void CollectTracks(const art::Event &evt, const std::string &label, TrackVector &trackVector, PFParticlesToTracks &particlesToTracks)
Collect the reconstructed PFParticles and associated Tracks from the ART event record.
static void CollectPFParticles(const art::Event &evt, const std::string &label, PFParticleVector &particleVector)
Collect the reconstructed PFParticles from the ART event record.
void GetStartAndEndPoints(const art::Ptr< simb::MCParticle > trueParticle, int &startT, int &endT) const
Find the start and end points of the true particle in the active region of detector.
static void CollectHits(const art::Event &evt, const std::string &label, HitVector &hitVector)
Collect the reconstructed Hits from the ART event record.
static void SelectNeutrinoPFParticles(const PFParticleVector &inputParticles, PFParticleVector &outputParticles)
Select reconstructed neutrino particles from a list of all reconstructed particles.
std::string m_particleLabel
std::vector< art::Ptr< recob::Hit > > HitVector
double m_pfoStraightLength
std::map< art::Ptr< recob::Hit >, art::Ptr< recob::PFParticle > > HitsToPFParticles
double m_mcStraightLength
std::vector< art::Ptr< recob::Vertex > > VertexVector
Planes which measure W (third view for Bo, MicroBooNE, etc).
std::vector< art::Ptr< anab::T0 > > T0Vector
static void CollectT0s(const art::Event &evt, const std::string &label, T0Vector &t0Vector, PFParticlesToT0s &particlesToT0s)
Collect a vector of T0s from the ART event record.
bool m_addDaughterPFParticles
static void BuildMCParticleHitMaps(const art::Event &evt, const HitVector &hitVector, const SimChannelVector &simChannelVector, HitsToTrackIDEs &hitsToTrackIDEs)
Collect the links from reconstructed hits to their true energy deposits.
static void CollectMCParticles(const art::Event &evt, const std::string &label, MCParticleVector &particleVector)
Collect a vector of MCParticle objects from the ART event record.
static art::Ptr< simb::MCParticle > GetFinalStateMCParticle(const MCParticleMap &particleMap, const art::Ptr< simb::MCParticle > daughterParticle)
Return the final-state parent particle by navigating up the chain of parent/daughter associations...
int m_nRecoWithoutTrueHits
Reconstructed hits which don't belong to any true particle - "missing".
BEGIN_PROLOG could also be cout
static bool IsFinalState(const PFParticleMap &particleMap, const art::Ptr< recob::PFParticle > daughterParticle)
Determine whether a particle has been reconstructed as a final-state particle.
std::map< art::Ptr< recob::SpacePoint >, art::Ptr< recob::Hit > > SpacePointsToHits
bool m_disableRealDataCheck
Whether to check if the input file contains real data before accessing MC information.
static void BuildMCParticleMap(const MCParticleVector &particleVector, MCParticleMap &particleMap)
Build particle maps for true particles.
1.8.5