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