CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
ConversionTrackPairFinder.cc
Go to the documentation of this file.
1 #include <iostream>
2 #include <vector>
3 #include <memory>
4 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionTrackPairFinder.h"
5 // Framework
6 #include "FWCore/MessageLogger/interface/MessageLogger.h"
7 //
8 
9 //
10 #include <vector>
11 #include <map>
12 
13 
14 //using namespace std;
15 
16 ConversionTrackPairFinder::ConversionTrackPairFinder( ){
17 
18  LogDebug("ConversionTrackPairFinder") << " CTOR " << "\n";
19 
20 }
21 
22 ConversionTrackPairFinder::~ConversionTrackPairFinder() {
23 
24  LogDebug("ConversionTrackPairFinder") << " DTOR " << "\n";
25 
26 }
27 
28 
29 
30 
31  std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors> ConversionTrackPairFinder::run(std::vector<reco::TransientTrack> outInTrk,
32  const edm::Handle<reco::TrackCollection>& outInTrkHandle,
33  const edm::Handle<reco::TrackCaloClusterPtrAssociation>& outInTrackSCAssH,
34  std::vector<reco::TransientTrack> inOutTrk,
35  const edm::Handle<reco::TrackCollection>& inOutTrkHandle,
36  const edm::Handle<reco::TrackCaloClusterPtrAssociation>& inOutTrackSCAssH )
37 {
38 
39 
40  LogDebug("ConversionTrackPairFinder") << "ConversionTrackPairFinder::run " << "\n";
41 
42  std::vector<reco::TransientTrack> selectedOutInTk;
43  std::vector<reco::TransientTrack> selectedInOutTk;
44  std::vector<reco::TransientTrack> allSelectedTk;
45  std::map<reco::TransientTrack, reco::CaloClusterPtr,CompareTwoTracks> scTrkAssocMap;
46  std::multimap<int,reco::TransientTrack,std::greater<int> > auxMap;
47 
48  bool oneLeg=false;
49  bool noTrack=false;
50 
51 
52 
53  int iTrk=0;
54  for( std::vector<reco::TransientTrack>::iterator iTk = outInTrk.begin(); iTk != outInTrk.end(); iTk++) {
55  edm::Ref<reco::TrackCollection> trackRef(outInTrkHandle, iTrk );
56  iTrk++;
57 
58  if ( iTk->numberOfValidHits() <3 || iTk->normalizedChi2() > 5000 ) continue;
59  if ( fabs(iTk->impactPointState().globalPosition().x()) > 110 ||
60  fabs(iTk->impactPointState().globalPosition().y()) > 110 ||
61  fabs(iTk->impactPointState().globalPosition().z()) > 280 ) continue;
62 
63  // std::cout << " Out In Track charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner pt " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";
64  const reco::TrackTransientTrack* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
65  reco::TrackRef myTkRef= ttt->persistentTrackRef();
66  //std::cout << " ConversionTrackPairFinder persistent track ref hits " << myTkRef->recHitsSize() << " inner pt " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";
67  // std::cout << " ConversionTrackPairFinder track from handle hits " << trackRef->recHitsSize() << " inner pt " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";
68 
69  const reco::CaloClusterPtr aClus = (*outInTrackSCAssH)[trackRef];
70 
71  // std::cout << "ConversionTrackPairFinder Reading the OutIn Map " << *outInTrackSCAss[trackRef] << " " << &outInTrackSCAss[trackRef] << std::endl;
72  // std::cout << "ConversionTrackPairFinder Out In track belonging to SC with energy " << aClus->energy() << "\n";
73 
74  int nHits=iTk->recHitsSize();
75  scTrkAssocMap[*iTk]= aClus;
76  auxMap.insert(std::pair<int,reco::TransientTrack >(nHits,(*iTk)) );
77  selectedOutInTk.push_back(*iTk);
78  allSelectedTk.push_back(*iTk);
79 
80 
81 
82  }
83 
84 
85  iTrk=0;
86  for( std::vector<reco::TransientTrack>::iterator iTk = inOutTrk.begin(); iTk != inOutTrk.end(); iTk++) {
87  edm::Ref<reco::TrackCollection> trackRef(inOutTrkHandle, iTrk );
88  iTrk++;
89 
90  if ( iTk->numberOfValidHits() <3 || iTk->normalizedChi2() >5000 ) continue;
91  if ( fabs(iTk->impactPointState().globalPosition().x()) > 110 ||
92  fabs(iTk->impactPointState().globalPosition().y()) > 110 ||
93  fabs(iTk->impactPointState().globalPosition().z()) > 280 ) continue;
94 
95  // std::cout << " In Out Track charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner pt " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";
96  const reco::TrackTransientTrack* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
97  reco::TrackRef myTkRef= ttt->persistentTrackRef();
98  // std::cout << " ConversionTrackPairFinder persistent track ref hits " << myTkRef->recHitsSize() << " inner pt " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";
99  // std::cout << " ConversionTrackPairFinder track from handle hits " << trackRef->recHitsSize() << " inner pt " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";
100 
101  const reco::CaloClusterPtr aClus = (*inOutTrackSCAssH)[trackRef];
102 
103  // std::cout << "ConversionTrackPairFinder Filling the InOut Map " << &(*inOutTrackSCAss[trackRef]) << " " << &inOutTrackSCAss[trackRef] << std::endl;
104  // std::cout << "ConversionTrackPairFinder In Out track belonging to SC with energy " << aClus.energy() << "\n";
105 
106  scTrkAssocMap[*iTk]= aClus;
107  int nHits=iTk->recHitsSize();
108  auxMap.insert(std::pair<int,reco::TransientTrack >(nHits,(*iTk)) );
109  selectedInOutTk.push_back(*iTk);
110  allSelectedTk.push_back(*iTk);
111 
112 
113 
114  }
115 
116 
117  // std::cout << " ConversionTrackPairFinder allSelectedTk size " << allSelectedTk.size() << " scTrkAssocMap size " << scTrkAssocMap.size() << "\n";
118 
119  // Sort tracks in decreasing number of hits
120  if(selectedOutInTk.size() > 0)
121  std::stable_sort(selectedOutInTk.begin(), selectedOutInTk.end(), ByNumOfHits());
122  if(selectedInOutTk.size() > 0)
123  std::stable_sort(selectedInOutTk.begin(), selectedInOutTk.end(), ByNumOfHits());
124  if(allSelectedTk.size() > 0)
125  std::stable_sort(allSelectedTk.begin(), allSelectedTk.end(), ByNumOfHits());
126 
127 
128 
129  // for( std::vector<reco::TransientTrack>::const_iterator iTk = selectedOutInTk.begin(); iTk != selectedOutInTk.end(); iTk++) {
130  // std::cout << " Selected Out In Tracks charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << "\n";
131  //}
132  // for( std::vector<reco::TransientTrack>::const_iterator iTk = selectedInOutTk.begin(); iTk != selectedInOutTk.end(); iTk++) {
133  // std::cout << " Selected In Out Tracks charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << "\n";
134  // }
135  // for( std::vector<reco::TransientTrack>::const_iterator iTk = allSelectedTk.begin(); iTk != allSelectedTk.end(); iTk++) {
136  // std::cout << " All Selected Tracks charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " chi2 " << iTk->normalizedChi2() << " pt " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";
137  //}
138 
139 
140 
141 
142  std::vector<reco::TransientTrack > thePair(2);
143  std::vector<std::vector<reco::TransientTrack> > allPairs;
144 
145  std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors > allPairSCAss;
146  std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors> allPairOrdInPtSCAss;
147  std::map<reco::TransientTrack, reco::CaloClusterPtr>::const_iterator iMap1;
148  std::map<reco::TransientTrack, reco::CaloClusterPtr>::const_iterator iMap2;
149 
150  for( iMap1 = scTrkAssocMap.begin(); iMap1 != scTrkAssocMap.end(); ++iMap1) {
151 
152  // std::cout << " Ass map track charge " << (iMap1->first).charge() <<" pt " << sqrt(((iMap1->first)).track().innerMomentum().Perp2()) << " SC E " << (iMap1->second)->energy() << " SC eta " << (iMap1->second)->eta() << " SC phi " << (iMap1->second)->phi() << std::endl;
153  }
154 
155 
156  std::multimap<int, reco::TransientTrack>::const_iterator iAux;
157 
158 
159  // for( iAux = auxMap.begin(); iAux!= auxMap.end(); ++iAux) {
160  // // std::cout << " Aux Map " << (iAux->first) <<" pt " << sqrt(((iAux->second)).track().innerMomentum().Perp2()) << std::endl;
161  // for( iMap1 = scTrkAssocMap.begin(); iMap1 != scTrkAssocMap.end(); ++iMap1) {
162  // if ( (iMap1->first) == (iAux->second) ) std::cout << " ass SC " << (iMap1->second)->energy() << std::endl;
163  // }
164  // }
165 
166  if ( scTrkAssocMap.size() > 2 ){
167 
168 
169  for( iMap1 = scTrkAssocMap.begin(); iMap1 != scTrkAssocMap.end(); ++iMap1) {
170  for( iMap2 = iMap1; iMap2 != scTrkAssocMap.end(); ++iMap2) {
171  // consider only tracks associated to the same SC
172 
173  if ( (iMap1->second) != (iMap2->second) ) continue;
174 
175  if ( ((iMap1->first)).charge() * ((iMap2->first)).charge() < 0 ) {
176 
177  // std::cout << " ConversionTrackPairFinde All selected from the map First Track charge " << (iMap1->first).charge() << " Num of RecHits " << ((iMap1->first)).recHitsSize() << " inner pt " << sqrt(((iMap1->first)).track().innerMomentum().Perp2()) << " Ass SC " << (iMap1->second)->energy() << "\n";
178 
179  // std::cout << " ConversionTrackPairFinde All selected from the map Second Track charge " << ((iMap2->first)).charge() << " Num of RecHits " << ((iMap2->first)).recHitsSize() << " inner pt " << sqrt(((iMap2->first)).track().innerMomentum().Perp2()) << " Ass SC " << (iMap2->second)->energy() << "\n";
180 
181 
182 
183  thePair.clear();
184  thePair.push_back( iMap1->first );
185  thePair.push_back( iMap2->first );
186  allPairs.push_back ( thePair );
187  allPairSCAss[thePair]= iMap1->second;
188 
189  }
190  }
191  }
192 
193  // std::cout << " ConversionTrackPairFinder INTERMIDIATE allPairSCAss size " << allPairSCAss.size() << "\n";
194 
195  if ( allPairSCAss.size() == 0) {
196  // std::cout << " All Tracks had the same charge: Need to send out a single track " << "\n";
197 
198  for( iMap1 = scTrkAssocMap.begin(); iMap1 != scTrkAssocMap.end(); ++iMap1) {
199 
200  thePair.clear();
201  thePair.push_back(iMap1->first);
202  allPairs.push_back ( thePair );
203  allPairSCAss[thePair]= iMap1->second;
204 
205  }
206 
207  }
208 
209 
210 
211 
212 
213  } else if ( (scTrkAssocMap.size() ==2) ) {
214 
215  iMap1=scTrkAssocMap.begin();//get the first
216  iMap2=iMap1;
217  iMap2++;//get the second
218  if ( (iMap1->second) == (iMap2->second) ) {
219  if ( (iMap1->first).charge() * (iMap2->first).charge() < 0 ) {
220 
221  // std::cout << " ConversionTrackPairFinder Case when (scTrkAssocMap.size() ==2) " << (iMap1->first).charge() << std::endl;
222  //std::cout << " Num of RecHits " << ((iMap1->first)).recHitsSize() << std::endl;
223  // std::cout << " inner pt " << sqrt(((iMap1->first)).track().innerMomentum().Perp2()) << std::endl;
224  //std::cout << " Ass SC " << (iMap1->second)->energy() << "\n";
225 
226  // std::cout << " ConversionTrackPairFinder Case when (scTrkAssocMap.size() ==2) " << (iMap2->first).charge() << std::endl;
227  // std::cout << " Num of RecHits " << ((iMap2->first)).recHitsSize() << std::endl;
228  //std::cout << " inner pt " << sqrt(((iMap2->first)).track().innerMomentum().Perp2()) << std::endl;
229  //std::cout << " Ass SC " << (iMap2->second)->energy() << "\n";
230 
231  thePair.clear();
232  thePair.push_back( iMap1->first );
233  thePair.push_back( iMap2->first );
234  allPairs.push_back ( thePair );
235 
236  allPairSCAss[thePair]= iMap1->second;
237 
238 
239  } else {
240  //std::cout << " ConversionTrackPairFinder oneLeg case when 2 tracks with same sign Pick up the longest one" << std::endl;
241  if ( ((iMap1->first)).recHitsSize() > ((iMap2->first)).recHitsSize() ) {
242  thePair.clear();
243  thePair.push_back(iMap1->first);
244  allPairs.push_back ( thePair );
245  allPairSCAss[thePair]= iMap1->second;
246  } else {
247  thePair.clear();
248  thePair.push_back(iMap2->first);
249  allPairs.push_back ( thePair );
250  allPairSCAss[thePair]= iMap2->second;
251  }
252  }
253 
254  }
255 
256  } else if (scTrkAssocMap.size() ==1 ) {
257  // std::cout << " ConversionTrackPairFinder oneLeg case when 1 track only " << std::endl;
258  oneLeg=true;
259  } else {
260  noTrack=true;
261  }
262 
263 
264  if ( oneLeg ) {
265  thePair.clear();
266  // std::cout << " ConversionTrackPairFinder oneLeg case charge " << std::endl;
267 
268 
269  iMap1=scTrkAssocMap.begin();
270 
271  //std::cout << " ConversionTrackPairFinder oneLeg case charge " << (iMap1->first).charge() << " Num of RecHits " << ((iMap1->first)).recHitsSize() << " inner pt " << sqrt(((iMap1->first)).track().innerMomentum().Perp2()) << " Ass SC " << (iMap1->second)->energy() << "\n";
272 
273  thePair.push_back(iMap1->first);
274  allPairs.push_back ( thePair );
275  allPairSCAss[thePair]= iMap1->second;
276 
277  // std::cout << " WARNING ConversionTrackPairFinder::tracks The candidate has just one leg. Need to find another way to evaltuate the vertex !!! " << "\n";
278  }
279 
280  if ( noTrack) {
281  // std::cout << " WARNING ConversionTrackPairFinder::tracks case noTrack " << "\n";
282  thePair.clear();
283  allPairSCAss.clear();
284  }
285 
286 
287 
289  for( iMap1 = scTrkAssocMap.begin(); iMap1 != scTrkAssocMap.end(); ++iMap1) {
290 
291  int nFound=0;
292  for ( std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr>::const_iterator iPair= allPairSCAss.begin(); iPair!= allPairSCAss.end(); ++iPair ) {
293  if ( (iMap1->second) == (iPair->second) ) nFound++;
294  }
295 
296  if ( nFound == 0) {
297  // std::cout << " nFound zero case " << std::endl;
298  int iList=0;
299  for( iAux = auxMap.begin(); iAux!= auxMap.end(); ++iAux) {
300  if ( (iMap1->first) == (iAux->second) && iList==0 ) {
301  thePair.clear();
302  thePair.push_back(iAux->second);
303  allPairSCAss[thePair]= iMap1->second;
304 
305  }
306 
307  iList++;
308  }
309  }
310 
311  }
312 
313 
314 
315 
316  // order the tracks in the pair in order of decreasing pt
317  for ( std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr>::const_iterator iPair= allPairSCAss.begin(); iPair!= allPairSCAss.end(); ++iPair ) {
318  thePair.clear();
319  if ( (iPair->first).size() ==2 ) {
320  if ( sqrt((iPair->first)[0].track().innerMomentum().perp2()) > sqrt((iPair->first)[1].track().innerMomentum().perp2()) ) {
321  thePair.push_back((iPair->first)[0]);
322  thePair.push_back((iPair->first)[1]);
323  } else {
324  thePair.push_back((iPair->first)[1]);
325  thePair.push_back((iPair->first)[0]);
326  }
327  } else {
328  thePair.push_back((iPair->first)[0]);
329  }
330 
331  allPairOrdInPtSCAss[thePair]=iPair->second;
332  }
333 
334 
335  // std::cout << " ConversionTrackPairFinder FINAL allPairOrdInPtSCAss size " << allPairOrdInPtSCAss.size() << "\n";
336  // for ( std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr>::const_iterator iPair= allPairOrdInPtSCAss.begin(); iPair!= allPairOrdInPtSCAss.end(); ++iPair ) {
337  // std::cout << " ConversionTrackPairFindder FINAL allPairOrdInPtSCAss " << (iPair->first).size() << " SC Energy " << (iPair->second)->energy() << " eta " << (iPair->second)->eta() << " phi " << (iPair->second)->phi() << "\n";
338  // std::cout << " ConversionTrackPairFindder FINAL allPairOrdInPtSCAss (iPair->first).size() " << (iPair->first).size() << std::endl;
339 
340  // for ( std::vector<reco::TransientTrack>::const_iterator iTk=(iPair->first).begin(); iTk!= (iPair->first).end(); ++iTk) {
341  // std::cout << " ConversionTrackPair ordered track pt " << sqrt(iTk->track().innerMomentum().perp2()) << std::endl;
342  // }
343  //}
344 
345 
346 
347  return allPairOrdInPtSCAss;
348 
349 
350 }
351 
352 
353 
354 
355 
356 
357 
358 
LogDebug
#define LogDebug(id)
Definition: MessageLogger.h:501
ConversionTrackPairFinder::run
std::map< std::vector< reco::TransientTrack >, reco::CaloClusterPtr, CompareTwoTracksVectors > run(std::vector< reco::TransientTrack > outIn, const edm::Handle< reco::TrackCollection > &outInTrkHandle, const edm::Handle< reco::TrackCaloClusterPtrAssociation > &outInTrackSCAssH, std::vector< reco::TransientTrack > inOut, const edm::Handle< reco::TrackCollection > &inOutTrkHandle, const edm::Handle< reco::TrackCaloClusterPtrAssociation > &inOutTrackSCAssH)
Definition: ConversionTrackPairFinder.cc:31
reco::CaloClusterPtr
edm::Ptr< CaloCluster > CaloClusterPtr
Definition: CaloClusterFwd.h:17
DeDxDiscriminatorTools::charge
double charge(const std::vector< uint8_t > &Ampls)
Definition: DeDxDiscriminatorTools.cc:43
Association.map
dictionary map
Definition: Association.py:205
reco::TrackTransientTrack::persistentTrackRef
TrackRef persistentTrackRef() const
Definition: TrackTransientTrack.h:62
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
Generated for CMSSW Reference Manual by doxygen 1.8.5