|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
|
#include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
#include "Calibration/IsolatedParticles/interface/MatchingSimTrack.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <sstream>
namespace spr {
edm::SimTrackContainer::const_iterator matchedSimTrack(const edm::Event& iEvent,
edm::Handle<edm::SimTrackContainer>& SimTk,
edm::Handle<edm::SimVertexContainer>& SimVtx,
const reco::Track* pTrack,
TrackerHitAssociator& associate,
bool debug) {
edm::SimTrackContainer::const_iterator itr = SimTk->end();
;
//Get the vector of PsimHits associated to TrackerRecHits and select the
//matching SimTrack on the basis of maximum occurance of trackIds
std::vector<unsigned int> trkId, trkOcc;
for (auto const& trkHit : pTrack->recHits()) {
std::vector<PSimHit> matchedSimIds = associate.associateHit(*trkHit);
for (unsigned int isim = 0; isim < matchedSimIds.size(); isim++) {
unsigned tkId = matchedSimIds[isim].trackId();
bool found = false;
for (unsigned int j = 0; j < trkId.size(); j++) {
if (tkId == trkId[j]) {
trkOcc[j]++;
found = true;
break;
}
}
if (!found) {
trkId.push_back(tkId);
trkOcc.push_back(1);
}
}
}
if (debug) {
std::ostringstream st1;
for (unsigned int isim = 0; isim < trkId.size(); isim++) {
st1 << "\n trkId " << trkId[isim] << " Occurance " << trkOcc[isim] << ", ";
}
edm::LogVerbatim("IsoTrack") << st1.str();
}
int matchedId = 0;
unsigned int matchSimTrk = 0;
if (!trkOcc.empty()) {
unsigned int maxTrkOcc = 0, idxMax = 0;
for (unsigned int j = 0; j < trkOcc.size(); j++) {
if (trkOcc[j] > maxTrkOcc) {
maxTrkOcc = trkOcc[j];
idxMax = j;
}
}
matchSimTrk = trkId[idxMax];
for (auto simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
if (simTrkItr->trackId() == matchSimTrk) {
matchedId = simTrkItr->type();
if (debug)
edm::LogVerbatim("IsoTrack") << "matched trackId (maximum occurance) " << matchSimTrk << " type "
<< matchedId;
itr = simTrkItr;
break;
}
}
}
if (matchedId == 0 && debug) {
edm::LogVerbatim("IsoTrack") << "Could not find matched SimTrk and track history now ";
}
return itr;
}
//Returns a vector of TrackId originating from the matching SimTrack
std::vector<int> matchedSimTrackId(const edm::Event& iEvent,
edm::Handle<edm::SimTrackContainer>& SimTk,
edm::Handle<edm::SimVertexContainer>& SimVtx,
const reco::Track* pTrack,
TrackerHitAssociator& associate,
bool debug) {
// get the matching SimTrack
edm::SimTrackContainer::const_iterator trkInfo =
spr::matchedSimTrack(iEvent, SimTk, SimVtx, pTrack, associate, debug);
unsigned int matchSimTrk = trkInfo->trackId();
if (debug)
edm::LogVerbatim("IsoTrack") << "matchedSimTrackId finds the SimTrk ID of the current track to be "
<< matchSimTrk;
std::vector<int> matchTkid;
if (trkInfo->type() != 0) {
for (auto simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
if (validSimTrack(matchSimTrk, simTrkItr, SimTk, SimVtx, false))
matchTkid.push_back(static_cast<int>(simTrkItr->trackId()));
}
}
return matchTkid;
}
spr::simTkInfo matchedSimTrackInfo(unsigned int simTkId,
edm::Handle<edm::SimTrackContainer>& SimTk,
edm::Handle<edm::SimVertexContainer>& SimVtx,
bool debug) {
spr::simTkInfo info;
for (edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
if (simTkId == simTrkItr->trackId()) {
if (spr::validSimTrack(simTkId, simTrkItr, SimTk, SimVtx, debug)) {
info.found = true;
info.pdgId = simTrkItr->type();
info.charge = simTrkItr->charge();
} else {
edm::SimTrackContainer::const_iterator parentItr = spr::parentSimTrack(simTrkItr, SimTk, SimVtx, debug);
if (debug) {
if (parentItr != SimTk->end())
edm::LogVerbatim("IsoTrack") << "original parent of " << simTrkItr->trackId() << " "
<< parentItr->trackId() << ", " << parentItr->type();
else
edm::LogVerbatim("IsoTrack") << "original parent of " << simTrkItr->trackId() << " not found";
}
if (parentItr != SimTk->end()) {
info.found = true;
info.pdgId = parentItr->type();
info.charge = parentItr->charge();
}
}
break;
}
}
return info;
}
// Checks if this SimTrack=thisTrkItr originates from the one with trackId=simTkId
bool validSimTrack(unsigned int simTkId,
edm::SimTrackContainer::const_iterator thisTrkItr,
edm::Handle<edm::SimTrackContainer>& SimTk,
edm::Handle<edm::SimVertexContainer>& SimVtx,
bool debug) {
if (debug)
edm::LogVerbatim("IsoTrack") << "Inside validSimTrack: trackId " << thisTrkItr->trackId() << " vtxIndex "
<< thisTrkItr->vertIndex() << " to be matched to " << simTkId;
//This track originates from simTkId
if (thisTrkItr->trackId() == simTkId)
return true;
//Otherwise trace back the history using SimTracks and SimVertices
int vertIndex = thisTrkItr->vertIndex();
if (vertIndex == -1 || vertIndex >= static_cast<int>(SimVtx->size()))
return false;
edm::SimVertexContainer::const_iterator simVtxItr = SimVtx->begin();
for (int iv = 0; iv < vertIndex; iv++)
simVtxItr++;
int parent = simVtxItr->parentIndex();
if (debug)
edm::LogVerbatim("IsoTrack") << "validSimTrack:: parent index " << parent << " ";
if (parent < 0 && simVtxItr != SimVtx->begin()) {
const math::XYZTLorentzVectorD pos1 = simVtxItr->position();
for (simVtxItr = SimVtx->begin(); simVtxItr != SimVtx->end(); ++simVtxItr) {
if (simVtxItr->parentIndex() > 0) {
const math::XYZTLorentzVectorD pos2 = pos1 - simVtxItr->position();
double dist = pos2.P();
if (dist < 0.001) {
parent = simVtxItr->parentIndex();
break;
}
}
}
}
if (debug)
edm::LogVerbatim("IsoTrack") << "final index " << parent;
for (edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
if (static_cast<int>(simTrkItr->trackId()) == parent && simTrkItr != thisTrkItr)
return validSimTrack(simTkId, simTrkItr, SimTk, SimVtx, debug);
}
return false;
}
//Returns the parent of a SimTrack
edm::SimTrackContainer::const_iterator parentSimTrack(edm::SimTrackContainer::const_iterator thisTrkItr,
edm::Handle<edm::SimTrackContainer>& SimTk,
edm::Handle<edm::SimVertexContainer>& SimVtx,
bool debug) {
edm::SimTrackContainer::const_iterator itr = SimTk->end();
int vertIndex = thisTrkItr->vertIndex();
if (debug)
edm::LogVerbatim("IsoTrack") << "SimTrackParent " << thisTrkItr->trackId() << " Vertex " << vertIndex << " Type "
<< thisTrkItr->type() << " Charge " << static_cast<int>(thisTrkItr->charge());
if (vertIndex == -1)
return thisTrkItr;
else if (vertIndex >= static_cast<int>(SimVtx->size()))
return itr;
edm::SimVertexContainer::const_iterator simVtxItr = SimVtx->begin();
for (int iv = 0; iv < vertIndex; iv++)
simVtxItr++;
int parent = simVtxItr->parentIndex();
if (parent < 0 && simVtxItr != SimVtx->begin()) {
const math::XYZTLorentzVectorD pos1 = simVtxItr->position();
for (simVtxItr = SimVtx->begin(); simVtxItr != SimVtx->end(); ++simVtxItr) {
if (simVtxItr->parentIndex() > 0) {
const math::XYZTLorentzVectorD pos2 = pos1 - simVtxItr->position();
double dist = pos2.P();
if (dist < 0.001) {
parent = simVtxItr->parentIndex();
break;
}
}
}
}
for (edm::SimTrackContainer::const_iterator simTrkItr = SimTk->begin(); simTrkItr != SimTk->end(); simTrkItr++) {
if (static_cast<int>(simTrkItr->trackId()) == parent && simTrkItr != thisTrkItr)
return parentSimTrack(simTrkItr, SimTk, SimVtx, debug);
}
return thisTrkItr;
}
} // namespace spr
|