5 #include <EVENT/LCCollection.h>
6 #include <EVENT/SimTrackerHit.h>
7 #include <EVENT/MCParticle.h>
8 #include <EVENT/TrackerHit.h>
9 #include <EVENT/LCRelation.h>
10 #include <UTIL/LCRelationNavigator.h>
21 #include <gear/BField.h>
23 #include <UTIL/BitField64.h>
24 #include <UTIL/ILDConf.h>
27 using namespace lcio ;
28 using namespace marlin ;
40 typedef std::vector< LCRelationNavigator* > Navec;
41 typedef std::vector< SimTrackerHit* > SimTrackerHitVec;
42 typedef std::map< MCParticle* , int > mymap;
43 typedef mymap::iterator myIter;
44 typedef std::multimap< int , MCParticle*> revMCPMap;
45 typedef revMCPMap::iterator revIter;
50 purityMCP GetPurity(TrackerHitVec tvec, Navec navec){
52 int ntrkhits = tvec.size();
53 SimTrackerHitVec simvec;
55 for(
int i = 0; i < ntrkhits; i++){
56 LCObject* trk =
dynamic_cast<LCObject*
>(tvec[i]);
57 const LCObjectVec* lcoVec = NULL;
58 int navsize = navec.size();
59 for(
int i = 0; i < navsize; i++){
61 lcoVec = &navec[i]->getRelatedToObjects(trk);
62 if(lcoVec->size() != 0){
67 int size = lcoVec->size();
73 for(
int jlc = 0; jlc < size; jlc++) simvec.push_back( dynamic_cast<SimTrackerHit*>( lcoVec->at(jlc) ) );
78 for(
int i = 0; i < int(simvec.size()); i++ ){
79 it = map.find(simvec[i]->getMCParticle());
82 if(it != map.end()){ it->second++; }
83 else map.insert(std::make_pair(simvec[i]->getMCParticle(),1));
89 for(myIter itm = map.begin(); itm != map.end(); itm++){
90 revmap.insert(std::make_pair(itm->second,itm->first));
92 MCParticle* dominantmcp = NULL;
94 std::cout <<
"revmap.size() : " << revmap.size() << std::endl;
96 if(revmap.size() != 0){
97 revIter maxit = std::max_element(revmap.begin(),revmap.end());
98 dominantmcp = maxit->second;
100 std::cout <<
"dominantmcp->getPDG() : " << dominantmcp->getPDG() << std::endl;
102 float ndmcp = float(maxit->first);
103 purity = ndmcp/(float(simvec.size()) +
float(bkgcount) );
106 pack.mcp = dominantmcp;
107 pack.purity = purity;
Definition: GetPurity.h:31
Definition: GetPurity.h:37