MarlinTrkProcessors  2.4.1
GetPurity.h
1 #ifndef GETPURITY_h
2 #define GETPURITY_h 1
3 
4 
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>
11 
12 #include <iostream>
13 #include <algorithm>
14 #include <functional>
15 #include <map>
16 #include <vector>
17 #include <cmath>
18 #include <climits>
19 
20 
21 #include <gear/BField.h>
22 
23 #include <UTIL/BitField64.h>
24 #include <UTIL/ILDConf.h>
25 
26 
27 using namespace lcio ;
28 using namespace marlin ;
29 
30 
31 class purityMCP{
32  public:
33  float purity;
34  MCParticle* mcp;
35 };
36 
38  public:
39  GetPurityUtil(){;}
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;
46  //navecはsimtrackerhit,trackerhitの奴を使いそうなやつだけpush_backしてからGetPurityの引数に入れる。
47 
48 
49 
50  purityMCP GetPurity(TrackerHitVec tvec, Navec navec){
51  //pMCPにはdominantMCPのアドレスを指した状態で返ってくる。
52  int ntrkhits = tvec.size();
53  SimTrackerHitVec simvec;
54  int bkgcount = 0;
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++){
60  if(navec[i] != NULL){
61  lcoVec = &navec[i]->getRelatedToObjects(trk);//見つからない場合はNULLを返すわけじゃないので、sizeが0かどうかで判定する。
62  if(lcoVec->size() != 0){
63  break;
64  }
65  }
66  }
67  int size = lcoVec->size();
68  simvec.reserve(size);
69  if(size == 0){
70  bkgcount++;
71  }
72  else{
73  for(int jlc = 0; jlc < size; jlc++) simvec.push_back( dynamic_cast<SimTrackerHit*>( lcoVec->at(jlc) ) );
74  }
75  }
76  mymap map;
77  myIter it;
78  for(int i = 0; i < int(simvec.size()); i++ ){
79  it = map.find(simvec[i]->getMCParticle());
80  //( it != map.end() ) ? ( it->second++ ) : ( map.insert(std::make_pair(simvec[i]->getMCParticle(),1)) );
81  //上の表記は何故かできない模様。
82  if(it != map.end()){ it->second++; }
83  else map.insert(std::make_pair(simvec[i]->getMCParticle(),1));
84  }
85 
86 
87  float purity = 0 ;
88  revMCPMap revmap;
89  for(myIter itm = map.begin(); itm != map.end(); itm++){
90  revmap.insert(std::make_pair(itm->second,itm->first));
91  }
92  MCParticle* dominantmcp = NULL;
93 #if 0
94  std::cout << "revmap.size() : " << revmap.size() << std::endl;
95 #endif
96  if(revmap.size() != 0){
97  revIter maxit = std::max_element(revmap.begin(),revmap.end());
98  dominantmcp = maxit->second;
99 #if 0
100  std::cout << "dominantmcp->getPDG() : " << dominantmcp->getPDG() << std::endl;
101 #endif
102  float ndmcp = float(maxit->first);
103  purity = ndmcp/(float(simvec.size()) + float(bkgcount) );
104  }
105  purityMCP pack;
106  pack.mcp = dominantmcp;
107  pack.purity = purity;
108  return pack;
109  }
110 };
111 
112 
113 
114 
115 
116 
117 
118 
119 #endif
Definition: GetPurity.h:31
Definition: GetPurity.h:37