diff --git a/Trigger/TrigTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEvent/TrigIDHitStats.h b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEvent/TrigIDHitStats.h new file mode 100644 index 0000000000000000000000000000000000000000..95e054ee88396f482eec53151bb51acf8c28eb5e --- /dev/null +++ b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEvent/TrigIDHitStats.h @@ -0,0 +1,73 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +/* +* TrigInDetTrackTruth Helper class to keep numbers of hits in ID sub-detectors +* +* Copied from Tracking/TrkEvent/TrkTruthData/TrkTruthData/SubDetHitStatistics.h +*/ + +#ifndef TRIGIDHITSTATS_TRUTH_H +#define TRIGIDHITSTATS_TRUTH_H + +#include <cstring> +#include "CLIDSvc/CLASS_DEF.h" + + +class TrigIDHitStats{ + + public: + + enum IDSubDetType { PIX=0, SCT, TRT, NUM_SUBDETECTORS }; + + /** Constructors: POOL needs default constructor */ + TrigIDHitStats() { + memset(numHits, 0, NUM_SUBDETECTORS); + } + + virtual ~TrigIDHitStats() { } + + unsigned char& operator[](IDSubDetType i) { return numHits[i]; } + + const unsigned char& operator[](IDSubDetType i) const { return numHits[i]; } + + unsigned int total() const { + unsigned char tot=0; + for(unsigned i=0; i<NUM_SUBDETECTORS; i++) { + tot += numHits[i]; + } + return (unsigned int)tot; + } + + unsigned int pixhits() const { return (unsigned int)numHits[PIX];} + unsigned int scthits() const { return (unsigned int)numHits[SCT];} + unsigned int trthits() const { return (unsigned int)numHits[TRT];} + + + TrigIDHitStats& operator+=(const TrigIDHitStats& b) { + for(unsigned i=0; i<NUM_SUBDETECTORS; i++) { + numHits[i] += b.numHits[i]; + } + return *this; + } + +private: + + // For InDet, the largest typical number of measurements per track is 36 in + // the TRT. One byte is enough to keep any of the numbers. + unsigned char numHits[NUM_SUBDETECTORS]; + +}; + +template<class Ostream> Ostream& operator<<(Ostream& os, const TrigIDHitStats& m) { + os<<"IDSubDetStat("; + for(unsigned i=0; i<TrigIDHitStats::NUM_SUBDETECTORS; i++) { + os<<unsigned(m[TrigIDHitStats::IDSubDetType(i)]); + if(i!=TrigIDHitStats::NUM_SUBDETECTORS) os<<","; + } + return os<<")"; +} + +CLASS_DEF( TrigIDHitStats , 212921831 , 1 ) +#endif diff --git a/Trigger/TrigTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEvent/TrigInDetTrackTruth.h b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEvent/TrigInDetTrackTruth.h new file mode 100644 index 0000000000000000000000000000000000000000..39ee6d583ad7c1792c89f3b36fae0f2c20452b34 --- /dev/null +++ b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEvent/TrigInDetTrackTruth.h @@ -0,0 +1,110 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + + +/************************************************************************** + ** + ** File: TrigInDetTrackTruth.h + ** + ** Description: - Stores a vector of pointers to GenParticles which match + ** a TrigInDetTrack and associated matching quality + ** quantities (nr of common hits only for now) + ** + ** Author: R.Goncalo + ** + ** Created: Sat Jan 18 19:55:56 GMT 2006 + ** Modified: + ** + ** + ** + **************************************************************************/ +#ifndef TRIGINDETTRACK_TRUTH_H +#define TRIGINDETTRACK_TRUTH_H + +#include "CLIDSvc/CLASS_DEF.h" + +#include "TrigInDetEvent/TrigInDetTrack.h" +#include "TrigInDetTruthEvent/TrigIDHitStats.h" + +#include "HepMC/GenParticle.h" +#include "HepMC/GenVertex.h" +#include "GeneratorObjects/HepMcParticleLink.h" + +#include "TrigSteeringEvent/MessageSvcProvider.h" +#include <iostream> + +using namespace std; + +class TrigInDetTrackTruth { + + public: + + /** Constructors: POOL needs default constructor */ + TrigInDetTrackTruth() : + best_match_hits(-1), best_Si_match_hits(-1), best_TRT_match_hits(-1), m_true_part_vec(0), m_nr_common_hits(0), m_family_tree() + { }; + + /** initialized constructor: easier way to construct an instance + if there is just one true particle associated with a track */ + TrigInDetTrackTruth(const HepMcParticleLink p_tru_part, TrigIDHitStats hits) + { + // add this particle to the vector + m_true_part_vec.push_back(p_tru_part); + + // add number of hits + m_nr_common_hits.push_back( hits ); + + // set the best match index (only 1 match so far) + best_match_hits = 0; + best_Si_match_hits = 0; + best_TRT_match_hits = 0; + }; + + // Destructor + virtual ~TrigInDetTrackTruth() { }; + + // accessors: to add links to GenParticles to object and get + // index if already filled; addMatch returns index of new entry + int addMatch(HepMcParticleLink p_tru_part, TrigIDHitStats hits); + int index(HepMcParticleLink&) const; + + // to get quantities + const HepMcParticleLink* bestMatch() const; + const HepMcParticleLink* bestSiMatch() const; + const HepMcParticleLink* bestTRTMatch() const; + const HepMcParticleLink* truthMatch(unsigned int i) const; + unsigned int nrMatches() const; + unsigned int nrCommonHits(unsigned int i) const; + unsigned int nrCommonSiHits(unsigned int i) const; + unsigned int nrCommonTRTHits(unsigned int i) const; + unsigned int nrCommonHitsBestSi() const; + unsigned int nrCommonHitsBestTRT() const; + + // to handle mother-daughter relationships within matching particles + int updateFamilyTree(); + std::vector< std::pair<unsigned int, unsigned int> > getFamilyTree() const; + bool motherInChain(unsigned int) const; + int motherIndexInChain(unsigned int) const; + bool daughtersInChain(unsigned int) const; + std::vector<unsigned int> daughterIndicesInChain(unsigned int) const; + + private: + + // reference best match quantities + int best_match_hits; + int best_Si_match_hits; + int best_TRT_match_hits; + + // vector of HepMcParticleLink pointers and matching quantities + std::vector<HepMcParticleLink> m_true_part_vec; + std::vector<TrigIDHitStats> m_nr_common_hits; + // vector<pair<int,int>> to act as bidirectional map "mother index <-> daughter index" + // where the indices refer to the GenPaticle positions in m_true_part_vec vector + std::vector< pair<unsigned int, unsigned int> > m_family_tree; + + +}; + +CLASS_DEF( TrigInDetTrackTruth , 243330893 , 1 ) +#endif diff --git a/Trigger/TrigTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEvent/TrigInDetTrackTruthMap.h b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEvent/TrigInDetTrackTruthMap.h new file mode 100644 index 0000000000000000000000000000000000000000..271bbf9047b6952da349fd0d88cf1ab3d0fa69ca --- /dev/null +++ b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEvent/TrigInDetTrackTruthMap.h @@ -0,0 +1,93 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + + +/************************************************************************** + ** + ** File: TrigInDetTrackTruthMap.h + ** + ** Description: - Stores a vector of pointers to GenParticles which match + ** a TrigInDetTrack and associated matching quality + ** quantities (nr of common hits only, for now) + ** + ** Author: R.Goncalo + ** + ** Created: Sat Jan 18 19:55:56 GMT 2006 + ** Modified: RG - 19 Jun 06 - changed the class to deal with persistency + ** it now behaves a bit like a: + ** map<ElementLink to TrigInDetTrack, TrigInDetTrackTruth> + ** + **************************************************************************/ +#ifndef TRIGINDETTRACK_TRUTH_MAP_H +#define TRIGINDETTRACK_TRUTH_MAP_H + +#include "CLIDSvc/CLASS_DEF.h" + +using std::size_t; +#include "DataModel/DataVector.h" +#include "DataModel/ElementLinkVector.h" + +#include "TrigInDetEvent/TrigInDetTrack.h" +#include "TrigInDetEvent/TrigInDetTrackCollection.h" +#include "TrigInDetTruthEvent/TrigInDetTrackTruth.h" +#include "HepMC/GenParticle.h" + +#include <iostream> + +using namespace std; + +class TrigInDetTrackTruthMap { + + + public: + + /** Constructors: POOL needs default constructor */ + TrigInDetTrackTruthMap() : + m_elink_vec(0), + m_truth_vec(0) + { }; + + // Destructor + virtual ~TrigInDetTrackTruthMap() { }; + + + /** accessors to fill map */ + void addMatch(const TrigInDetTrackCollection* trkColl, + unsigned int trk_indx, + TrigInDetTrackTruth& p_trk_tru); + + /** methods to get truth-match objects */ + // returns true if truth association exists for track + bool hasTruth(const TrigInDetTrack* p_trig_trk) const; + + // returns the track truth association object + const TrigInDetTrackTruth* truth(const TrigInDetTrack* p_trig_trk) const; + + // to make the map more useful: return the link to GenParticle which better + // matches this track according to number of common hits + const HepMcParticleLink* bestMatchSi(const TrigInDetTrack* p_trig_trk) const; + const HepMcParticleLink* bestMatchTRT(const TrigInDetTrack* p_trig_trk) const; + + // return number of hits for best match + int bestMatchSiHits(const TrigInDetTrack* p_trig_trk) const; + int bestMatchTRTHits(const TrigInDetTrack* p_trig_trk) const; + + void print() const; + + size_t size() const; + const TrigInDetTrackTruth* truthi (size_t i) const; + const ElementLink<TrigInDetTrackCollection> trackiLink (size_t i) const; + const TrigInDetTrack* tracki (size_t i) const; + + private: + + // used to be map for fast lookup, but changed to use ElemLink as key + // std::vector< ElementLink< DataVector<TrigInDetTrack> > > m_elink_vec; + //std::vector< ElementLink< TrigInDetTrackCollection > > m_elink_vec; + ElementLinkVector< TrigInDetTrackCollection > m_elink_vec; + std::vector< TrigInDetTrackTruth > m_truth_vec; +}; + +CLASS_DEF( TrigInDetTrackTruthMap , 78130186 , 1 ) +#endif diff --git a/Trigger/TrigTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEventDict.h b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEventDict.h new file mode 100644 index 0000000000000000000000000000000000000000..b185264012771ddfc7c91dee95b98698c85ee84c --- /dev/null +++ b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEventDict.h @@ -0,0 +1,13 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + +#ifndef TRIGINDETTRACKTRUTHDICT_H +#define TRIGINDETTRACKTRUTHDICT_H + + +#include "TrigInDetTruthEvent/TrigIDHitStats.h" +#include "TrigInDetTruthEvent/TrigInDetTrackTruth.h" +#include "TrigInDetTruthEvent/TrigInDetTrackTruthMap.h" + +#endif diff --git a/Trigger/TrigTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEvent/selection.xml b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEvent/selection.xml new file mode 100644 index 0000000000000000000000000000000000000000..b9f38424dbd5e8b36834a216671a16725a32cdc7 --- /dev/null +++ b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/TrigInDetTruthEvent/selection.xml @@ -0,0 +1,11 @@ +<lcgdict> + <class name="TrigInDetTrackTruthMap" id="41581666-F06D-44AE-93B9-D7E912A27AA1" /> + <class name="std::vector< TrigInDetTrackTruth >" /> + <class name="TrigInDetTrackTruth" > + <field name="m_true_part_vec" transient="false" /> + <field name="m_nr_common_hits" transient="false" /> + <field name="m_family_tree" transient="false" /> + </class> + <class name="std::vector<TrigIDHitStats>" /> + <class name="TrigIDHitStats" /> +</lcgdict> diff --git a/Trigger/TrigTruthEvent/TrigInDetTruthEvent/cmt/requirements b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/cmt/requirements new file mode 100644 index 0000000000000000000000000000000000000000..cf23894aab563fa5ceb6edcf346032db19aa7e97 --- /dev/null +++ b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/cmt/requirements @@ -0,0 +1,31 @@ +package TrigInDetTruthEvent + +author Ricardo Goncalo <r.goncalo@rhul.ac.uk> + + +public +use AtlasPolicy AtlasPolicy-* +use AtlasHepMC AtlasHepMC-* External +use CLIDSvc CLIDSvc-* Control +use DataModel DataModel-* Control +use GeneratorObjects GeneratorObjects-* Generators +use TrigInDetEvent TrigInDetEvent-* Trigger/TrigEvent +use TrigSteeringEvent TrigSteeringEvent-* Trigger/TrigEvent +private +use AthenaKernel AthenaKernel-* Control +use AtlasReflex AtlasReflex-* External -no_auto_imports +use GaudiInterface GaudiInterface-* External +end_private + +# library: +library TrigInDetTruthEvent *.cxx +apply_pattern installed_library + +# generate dictionary fillers (pool converters in TrigEventAthenaPool) +private + +apply_pattern lcgdict dict=TrigInDetTruthEvent selectionfile=selection.xml \ + headerfiles="../TrigInDetTruthEvent/TrigInDetTruthEventDict.h" + +end_private + diff --git a/Trigger/TrigTruthEvent/TrigInDetTruthEvent/src/TrigInDetTrackTruth.cxx b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/src/TrigInDetTrackTruth.cxx new file mode 100644 index 0000000000000000000000000000000000000000..6ba64f96fbeba9c0eca33e3287f7685b6ab99803 --- /dev/null +++ b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/src/TrigInDetTrackTruth.cxx @@ -0,0 +1,293 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + + +/************************************************************************** + ** + ** File: TrigInDetTrackTruth.cxx + ** + ** Description: - Stores a vector of pointers to GenParticles which match + ** a TrigInDetTrack and associated matching quality + ** quantities (nr of common hits only for now) + ** + ** Author: R.Goncalo + ** + ** Created: Sat Jan 18 19:55:56 GMT 2006 + ** Modified: RG 1/9/06 - added method for object to update own family tree + ** RG 2/9/06 - was going to change set of two vectors acting + ** as a map<HepMcPart.Link,TrigIDHitStats> into a GaudiUtils::VectorMap to + ** improve readout speed but decided against this since it complicates + ** matters and creates dependencies, and because only a few (~1-5) GenParticles + ** can be expected to match any TrigInDetTrack, so the overhead is small. + ** Changed m_family_tree from map<unsigned int,unsigned int> into + ** vector< pair<unsigned int,unsigned int>> to avoid persistency problems and + ** also to make it easier to use either mother or daughter indices as keys + **************************************************************************/ + +#include "TrigInDetTruthEvent/TrigInDetTrackTruth.h" +#include "AthenaKernel/getMessageSvc.h" +#include "GaudiKernel/IMessageSvc.h" + +/** accessor to fill object: returns index of new entry in vectors */ +int TrigInDetTrackTruth::addMatch(HepMcParticleLink p_tru_part,TrigIDHitStats hits) +{ + std::string thisName("TrigInDetTrackTruth::addMatch"); + MsgStream log(Athena::getMessageSvc(), thisName); + + log << MSG::DEBUG<< "Inserting HepMcParticleLink and TrigIDHitStats to TrigInDetTrackTruth map" << endreq; + + int indx = index(p_tru_part); + if ( indx >= 0 ) { + // HepMcParticleLink already exists: replace in vectors + log << MSG::DEBUG<< "HepMcParticleLink already in map: replacing" + << endreq; + m_true_part_vec[indx] = p_tru_part; + m_nr_common_hits[indx] = hits; + } else { + // push back into vectors to use HepMcParticleLink as key + m_true_part_vec.push_back(p_tru_part); + m_nr_common_hits.push_back( hits ); + indx = m_true_part_vec.size() - 1; + } + + // to know the best match in total nr.hits + if ( best_match_hits == -1 || + (unsigned int)hits.total() > nrCommonHits(best_match_hits) ) + best_match_hits = m_true_part_vec.size() - 1; + + if ( best_Si_match_hits == -1 || + (unsigned int)(hits.pixhits() + hits.scthits())> nrCommonSiHits(best_Si_match_hits) ) + best_Si_match_hits = m_true_part_vec.size() - 1; + + if ( best_TRT_match_hits == -1 || + (unsigned int)hits.trthits() > nrCommonTRTHits(best_TRT_match_hits) ) + best_TRT_match_hits = m_true_part_vec.size() - 1; + + return indx; +} + +/** method to find if a given HepMcParticleLink already exists in "map" and, if + so, what is its index; if not found retruns -1; relies on + HepMcParticleLink overloaded == method (i.e. compares only barcode) */ +int TrigInDetTrackTruth::index(HepMcParticleLink& hep_link) const { + + std::vector<HepMcParticleLink>::const_iterator it1, end=m_true_part_vec.end(); + int indx = 0; + for (it1=m_true_part_vec.begin(); it1 != end; ++it1) { + if (hep_link == (*it1)) { + // HepMcParticleLink found in vector: return index + return indx; + } + ++indx; + } + // HepMcParticleLink not found: return default -1 + return -1; +} + +/** accessor to fill family tree: for each HepMcParticleLink in the internal + vector of HepMC::GenParticles, this method searches for its mother in the + same vector. If the mother is in the vector, the method adds the relation + as a a mother-daughter pair to the family tree vector + note: this method can be used as many times as necessary as long as the + association object is not constant */ +int TrigInDetTrackTruth::updateFamilyTree() +{ + + std::string thisName("TrigInDetTrackTruth::updateFamilyTree"); + MsgStream log(Athena::getMessageSvc(), thisName); + + log << MSG::DEBUG<< "In TrigInDetTrackTruth::updateFamilyTree()" << endreq; + + int nr_mothers_found=0; + int child=-1; + std::vector<HepMcParticleLink>::iterator it1, end=m_true_part_vec.end(); + + for (it1=m_true_part_vec.begin(); it1 != end; ++it1) + { + child++; + /* get production vertex GenParticle pointed to by this link */ + log << MSG::DEBUG<< "Looking for mother of matching particle nr "<<child<<endreq; + + // first get GenParticle pointer + if ( !it1->isValid() ) continue; + + const HepMC::GenParticle* p_child = (*it1); + log << MSG::DEBUG << "GenParticle " << child << " (" << p_child << "); PDG id=" + << p_child->pdg_id() << "; status=" << p_child->status() + << "; pT=" << p_child->momentum().perp() + << "; searches mother..." + << endreq; + + // then get production vertex (check against null) + HepMC::GenVertex* p_child_vtx = p_child->production_vertex(); + if ( p_child_vtx == NULL ) + { + log << MSG::DEBUG<<"GenVertex pointer null: jump to next particle"<<endreq; + continue; + } + log << MSG::DEBUG<< "GenParticle "<< child << " comes from vertex with pointer " + << p_child_vtx << endreq; + + /* find mother: there should be only one for final state particles + (particles which can leave energy deposits in detectors) */ + HepMC::GenVertex::particles_in_const_iterator p_mum = p_child_vtx->particles_in_const_begin(); + + // check a mother was found + if ( p_mum == p_child_vtx->particles_in_const_end() ) + { + log << MSG::DEBUG<< "Mother not found: go to next particle" <<endreq; + continue; + } + log << MSG::DEBUG<< "Mother GenParticle (" << *p_mum << ") found; PDG id=" + << (*p_mum)->pdg_id() << "; status=" << (*p_mum)->status() + << "; pT=" << (*p_mum)->momentum().perp() + << "; does it match track?" + << endreq; + // mother is (*p_mum); still have to see if it is a match to this track + std::vector<HepMcParticleLink>::iterator it2=m_true_part_vec.begin(); + + bool mum_found=false; + for (unsigned int mum=0; it2 != end; ++it2, ++mum) + { + log << MSG::DEBUG << "* Trying daughter index=" << child + << " and mother index=" << mum << endreq; + if ( *p_mum == *it2 ) + { // mother also matches track + m_family_tree.push_back( std::pair<unsigned int, unsigned int>(mum,child) ); + mum_found=true; + nr_mothers_found++; + + log << MSG::DEBUG << "* Mother also matches track! " + << nr_mothers_found + << " mother-daughter relations found so far" << endreq; + log << MSG::DEBUG << "Daughter "<< child <<" (PDG id=" + << p_child->pdg_id() << "; pT=" << p_child->momentum().perp() + << ") comes from mother " << mum << " (PDG id=" + << (*p_mum)->status() << "; pT=" << p_child->momentum().perp() + << ")" << endreq; + } + } + if (!mum_found) log << MSG::DEBUG << "* Mother doesn't match track" + << endreq; + } + return nr_mothers_found; +} + +/** returns best match according to the number of hits */ +const HepMcParticleLink* TrigInDetTrackTruth::bestMatch() const +{ + return &(m_true_part_vec[best_match_hits]); +} + +/** returns best match according to the number of hits */ +const HepMcParticleLink* TrigInDetTrackTruth::bestSiMatch() const +{ + return &(m_true_part_vec[best_Si_match_hits]); +} + +/** returns best match according to the number of hits */ +const HepMcParticleLink* TrigInDetTrackTruth::bestTRTMatch() const +{ + return &(m_true_part_vec[best_TRT_match_hits]); +} + +/** returns matching true particle number i */ +const HepMcParticleLink* TrigInDetTrackTruth::truthMatch(unsigned int i) const +{ + if (i < m_true_part_vec.size()) return &(m_true_part_vec[i]); + else return NULL; +} + + +/** returns number of common hits from true particle i and TrigInDetTrack */ +unsigned int TrigInDetTrackTruth::nrCommonHits(unsigned int i) const +{ + if (i < m_true_part_vec.size()) return m_nr_common_hits[i].total(); + else return 0; +} + +/** returns number of common hits from true particle i and TrigInDetTrack */ +unsigned int TrigInDetTrackTruth::nrCommonSiHits(unsigned int i) const +{ + if (i < m_true_part_vec.size()) return (m_nr_common_hits[i].pixhits() + m_nr_common_hits[i].scthits()); + else return 0; +} + +/** returns number of common hits from true particle i and TrigInDetTrack */ +unsigned int TrigInDetTrackTruth::nrCommonTRTHits(unsigned int i) const +{ + if (i < m_true_part_vec.size()) return m_nr_common_hits[i].trthits(); + else return 0; +} + +/** returns total number of common hits from best match true particle and TrigInDetTrack */ +unsigned int TrigInDetTrackTruth::nrCommonHitsBestSi() const +{ + return (m_nr_common_hits[best_Si_match_hits].scthits() + m_nr_common_hits[best_Si_match_hits].pixhits()); +} + +/** returns total number of common hits from best match true particle and TrigInDetTrack */ +unsigned int TrigInDetTrackTruth::nrCommonHitsBestTRT() const +{ + return m_nr_common_hits[best_TRT_match_hits].trthits(); +} + + +/** returns number of matching particles */ +unsigned int TrigInDetTrackTruth::nrMatches() const +{ + return m_true_part_vec.size(); +} + +/** returns copy of family tree "map" */ +std::vector< std::pair<unsigned int, unsigned int> > TrigInDetTrackTruth::getFamilyTree() const +{ + return m_family_tree; +} + +/** given index of a GenParticle which matches the track returns index of + its mother, if it also matches the track, or -1 if not */ +int TrigInDetTrackTruth::motherIndexInChain(unsigned int daughter) const { + + if ( m_family_tree.empty() ) return -1; + + std::vector< pair<unsigned int, unsigned int> >::const_iterator it,it_end = m_family_tree.end(); + + for (it = m_family_tree.begin(); it != it_end; ++it) { + if (daughter == (*it).second) return (int)(*it).first; + } + return -1; +} + +/** given index of a GenParticle which matches the track returns true + if its mother also matches the track and false if not */ +bool TrigInDetTrackTruth::motherInChain(unsigned int daughter) const { + return ( motherIndexInChain(daughter) >=0 ); +} + +/** given index of a GenParticle which matches the track returns vector with indices + of its daughters, if they also matches the track, or an empty vector if not */ +std::vector<unsigned int> TrigInDetTrackTruth::daughterIndicesInChain(unsigned int mother) const { + + std::vector<unsigned int> v_indx; + v_indx.clear(); + if ( m_family_tree.empty() ) return v_indx; + + std::vector< pair<unsigned int, unsigned int> >::const_iterator it,it_end = m_family_tree.end(); + unsigned int i=0; + for (it = m_family_tree.begin(); it != it_end; ++it) { + if (mother == (*it).first) { + v_indx.push_back(i); + } + ++i; + } + return v_indx; +} + +/** given index of a GenParticle which matches the track returns true + if it has stable a daughter which also matches the track and false if not */ +bool TrigInDetTrackTruth::daughtersInChain(unsigned int mother) const { + std::vector<unsigned int> v_indx = daughterIndicesInChain(mother); + return ( !(v_indx.empty()) ); +} diff --git a/Trigger/TrigTruthEvent/TrigInDetTruthEvent/src/TrigInDetTrackTruthMap.cxx b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/src/TrigInDetTrackTruthMap.cxx new file mode 100644 index 0000000000000000000000000000000000000000..b8fa32d4a03f6047ccbcb3e5627f44495f21147f --- /dev/null +++ b/Trigger/TrigTruthEvent/TrigInDetTruthEvent/src/TrigInDetTrackTruthMap.cxx @@ -0,0 +1,262 @@ +/* + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration +*/ + + +/************************************************************************** + ** + ** File: TrigInDetTrackTruthMap.cxx + ** + ** Description: - Stores a vector of pointers to GenParticles which match + ** a TrigInDetTrack and associated matching quality + ** quantities (nr of common hits only for now) + ** + ** Author: R.Goncalo + ** + ** Created: Sat Jan 18 19:55:56 GMT 2006 + ** Modified: + ** + ** + ** + **************************************************************************/ + +#include "TrigInDetTruthEvent/TrigInDetTrackTruthMap.h" +#include "AthenaKernel/getMessageSvc.h" +#include "GaudiKernel/IMessageSvc.h" + +#include <iomanip> +#include <iostream> +#include <cassert> + +/** accessors to fill map */ +void TrigInDetTrackTruthMap::addMatch(const TrigInDetTrackCollection* trkColl, + unsigned int trk_indx, + TrigInDetTrackTruth& p_trk_tru) +{ + std::string thisName("TrigInDetTrackTruthMap::addMatch"); + MsgStream log(Athena::getMessageSvc(), thisName); + + log << MSG::DEBUG << "Adding TrigInDetTrack to map" << endreq; + + // check entry for this track doesn't exist (otherwise we have a multimap) + const TrigInDetTrack* p_trig_trk((*trkColl)[trk_indx]); + + if ( !hasTruth( p_trig_trk ) ) { + ElementLink< TrigInDetTrackCollection > p_trk_lnk(*trkColl,trk_indx); + m_elink_vec.push_back(p_trk_lnk); + bool status = m_elink_vec.back().toPersistent(); + if (!status) + log << MSG::DEBUG << "ERROR: could not set ElementLink persistent" << endreq; + m_truth_vec.push_back(p_trk_tru); + } else { + log << MSG::DEBUG << "TrigInDetTrack already in map!" << endreq; + } +} + +/** methods to get truth-match objects */ +// returns true if truth association exists for track +bool TrigInDetTrackTruthMap::hasTruth(const TrigInDetTrack* p_trig_trk) const +{ + // must loop over map because ElementLink is a unidirectional link + bool found=false; + for (unsigned int i=0; i < m_elink_vec.size(); ++i) { + // if ( p_trig_trk == *(m_elink_vec[i]) ) { + if(!m_elink_vec[i].isValid()) + continue; + if ((*(m_elink_vec[i]))->param()) { + if ((*(m_elink_vec[i]))->algorithmId() == p_trig_trk->algorithmId() && + (*(m_elink_vec[i]))->param()->pT() == p_trig_trk->param()->pT() && + (*(m_elink_vec[i]))->param()->eta() == p_trig_trk->param()->eta() && + (*(m_elink_vec[i]))->param()->phi0() == p_trig_trk->param()->phi0()) { + found=true; + break; + } + } + } + // true if this track pointer exists at least once in the vector + return (found); +} + + +// returns the track truth association object +const TrigInDetTrackTruth* TrigInDetTrackTruthMap::truth(const TrigInDetTrack* p_trig_trk) const +{ + std::string thisName("TrigInDetTrackTruthMap::truth"); + MsgStream log(Athena::getMessageSvc(), thisName); + + log << MSG::DEBUG<<"Searching truth for track at ptr="<<p_trig_trk<<endreq; + + // must loop over map because ElementLink is a unidirectional link + for (unsigned int i=0; i < m_elink_vec.size(); ++i) { + // if ( p_trig_trk == *(m_elink_vec[i]) ) { + if(!m_elink_vec[i].isValid()) + continue; + if ((*(m_elink_vec[i]))->param()) { + if ((*(m_elink_vec[i]))->algorithmId() == p_trig_trk->algorithmId() && + (*(m_elink_vec[i]))->param()->pT() == p_trig_trk->param()->pT() && + (*(m_elink_vec[i]))->param()->eta() == p_trig_trk->param()->eta() && + (*(m_elink_vec[i]))->param()->phi0() == p_trig_trk->param()->phi0()) { + // found position in vector corresponding to this track pointer + log << MSG::DEBUG << "Truth match for track at ptr=" << p_trig_trk + << " found in map at index " << i <<endreq; + return &m_truth_vec[i]; + } + } + } + // didn't find it: return null pointer + log << MSG::DEBUG <<"Truth match for track at ptr="<<p_trig_trk <<" not in map"<<endreq; + return NULL; +} + + +// to make the map more useful: return the link to GenParticle which better +// matches this track according to number of common hits +const HepMcParticleLink* TrigInDetTrackTruthMap::bestMatchSi(const TrigInDetTrack* p_trig_trk) const +{ + if ( hasTruth(p_trig_trk) ) { + return ( truth(p_trig_trk)->bestSiMatch() ); + } else { + return NULL; + } +} + +int TrigInDetTrackTruthMap::bestMatchSiHits(const TrigInDetTrack* p_trig_trk) const +{ + if ( hasTruth(p_trig_trk) ) { + return ( truth(p_trig_trk)->nrCommonHitsBestSi() ); + } else { + return 0; + } +} + +// to make the map more useful: return the link to GenParticle which better +// matches this track according to number of common hits +const HepMcParticleLink* TrigInDetTrackTruthMap::bestMatchTRT(const TrigInDetTrack* p_trig_trk) const +{ + if ( hasTruth(p_trig_trk) ) { + return ( truth(p_trig_trk)->bestTRTMatch() ); + } else { + return NULL; + } +} + +int TrigInDetTrackTruthMap::bestMatchTRTHits(const TrigInDetTrack* p_trig_trk) const +{ + if ( hasTruth(p_trig_trk) ) { + return ( truth(p_trig_trk)->nrCommonHitsBestTRT() ); + } else { + return 0; + } +} + + +void TrigInDetTrackTruthMap::print() const +{ + std::string thisName("TrigInDetTrackTruthMap::print"); + MsgStream log(Athena::getMessageSvc(), thisName); + + std::ostringstream oss; + oss << "TrigInDetTruthMap: " << m_elink_vec.size() + << " track-truth associations" << std::endl; + + oss << "---------------------------------------------------------------------------------------------------------------------------------" << std::endl; + oss << "#track|algo| pointer | pT | eta | phi |#match|mother|Sihits|TRThits|ev.index| barcode | pdg id | pT | eta | phi |"<< std::endl; + for (unsigned int i=0; i < m_elink_vec.size(); ++i) { + oss << std::setiosflags(std::ios::dec) << std::setw(6) << i << "|" + << std::setiosflags(std::ios::dec) << std::setw(4) << (*(m_elink_vec[i]))->algorithmId() << "|" + << std::setw(11) << *(m_elink_vec[i]) << "|"; + + if(!m_elink_vec[i].isValid()) { + oss << std::setiosflags(std::ios::dec) << "Invalid TrigInDetTrack link !" + << std::setiosflags(std::ios::dec) << std::setw(8) << "|"; + } + else { + // print reconstructed track parameters + if ((*(m_elink_vec[i]))->param()) { + oss << std::setiosflags(std::ios::dec) << std::setw(14) << (*(m_elink_vec[i]))->param()->pT() << "|" + << std::setiosflags(std::ios::dec) << std::setw(10) << (*(m_elink_vec[i]))->param()->eta() << "|" + << std::setiosflags(std::ios::dec) << std::setw(10) << (*(m_elink_vec[i]))->param()->phi0()<< "|"; + } else { + oss << std::setiosflags(std::ios::dec) << std::setw(15) << "|" + << std::setiosflags(std::ios::dec) << std::setw(11) << "|" + << std::setiosflags(std::ios::dec) << std::setw(11) << "|"; + } + } + + if (m_truth_vec[i].nrMatches() == 0) oss << std::endl; + // print parameters of truth particles which have hits in common with track + for (unsigned int j=0; j < m_truth_vec[i].nrMatches(); ++j) { + + // find if this matching true particle has a mother which also matches the track + HepMcParticleLink p_link( *((m_truth_vec[i]).truthMatch(j)) ); + int child_indx = (m_truth_vec[i]).index(p_link); + int mother_indx = -1; + if (child_indx >= 0) mother_indx = (m_truth_vec[i]).motherIndexInChain(child_indx); + + if (j>0) { + + oss << std::setiosflags(std::ios::dec) << std::setw(7) << "|" + << std::setiosflags(std::ios::dec) << std::setw(5) << "|"; + oss << std::setiosflags(std::ios::dec) << std::setw(15) << "|" + << std::setiosflags(std::ios::dec) << std::setw(11) << "|" + << std::setiosflags(std::ios::dec) << std::setw(11) << "|"; + } + + oss << std::setiosflags(std::ios::dec) << std::setw(6) << j+1 << "|"; + + if (mother_indx >= 0) { + oss << std::setiosflags(std::ios::dec) << std::setw(6) << mother_indx << "|"; + } else { + oss << " -- |"; + } + oss << std::setiosflags(std::ios::dec) << std::setw(6) << (m_truth_vec[i]).nrCommonSiHits(j) << "|" + << std::setiosflags(std::ios::dec) << std::setw(7) << (m_truth_vec[i]).nrCommonTRTHits(j) << "|" + << std::setiosflags(std::ios::dec) << std::setw(8) << (m_truth_vec[i]).truthMatch(j)->eventIndex() << "|" + << std::setiosflags(std::ios::dec) << std::setw(9) << (m_truth_vec[i]).truthMatch(j)->barcode() << "|"; + + if ((m_truth_vec[i]).truthMatch(j)->cptr()) { + + oss << std::setiosflags(std::ios::dec) << std::setw(10) << (m_truth_vec[i]).truthMatch(j)->cptr()->pdg_id() << "|" + << std::setiosflags(std::ios::dec) << std::setw(14) << (m_truth_vec[i]).truthMatch(j)->cptr()->momentum().perp()<< "|" + << std::setiosflags(std::ios::dec) << std::setw(10) << (m_truth_vec[i]).truthMatch(j)->cptr()->momentum().eta() << "|" + << std::setiosflags(std::ios::dec) << std::setw(10) << (m_truth_vec[i]).truthMatch(j)->cptr()->momentum().phi() << "|"; + } else { + oss << std::setiosflags(std::ios::dec) << std::setw(11) << "|" + << std::setiosflags(std::ios::dec) << std::setw(15) << "|" + << std::setiosflags(std::ios::dec) << std::setw(11) << "|" + << std::setiosflags(std::ios::dec) << std::setw(11) << "|"; + } + oss << std::endl; + } + } + oss << "---------------------------------------------------------------------------------------------------------------------------------" << std::endl; + + log << MSG::DEBUG << oss.str() << endreq; +} + +size_t TrigInDetTrackTruthMap::size() const +{ + assert (m_elink_vec.size() == m_truth_vec.size()); + return m_elink_vec.size(); +} + +const TrigInDetTrackTruth* +TrigInDetTrackTruthMap::truthi (size_t i) const +{ + assert (i < m_truth_vec.size()); + return &m_truth_vec[i]; +} + +const TrigInDetTrack* +TrigInDetTrackTruthMap::tracki (size_t i) const +{ + assert (i < m_elink_vec.size()); + return *m_elink_vec[i]; +} + +const ElementLink<TrigInDetTrackCollection> +TrigInDetTrackTruthMap::trackiLink (size_t i) const +{ + assert (i < m_elink_vec.size()); + return m_elink_vec[i]; +}