108 {
109
110 using TLC = std::vector<ElementLink<DataVector<xAOD::TrackParticle_v1>>>;
111 using TL = ElementLink<DataVector<xAOD::TrackParticle_v1>>;
112
113
114 SG::ReadDecorHandle<xAOD::JetContainer, TLC> trackLinksHandle(
m_trackLinksKey, ctx);
115 SG::ReadDecorHandle<xAOD::JetContainer, std::vector<char>> trackOriginsHandle(
m_trackOriginsKey, ctx);
116 SG::ReadDecorHandle<xAOD::JetContainer, std::vector<char>> vertexLinksHandle(
m_vertexLinksKey, ctx);
117 SG::WriteDecorHandle<xAOD::JetContainer, std::vector<ElementLink<xAOD::VertexContainer>>>
119
120
121 for (
const xAOD::Jet* jet : *inJetContainer) {
122
123
124 auto vertexCollection = vertexLinksHandle(*jet);
125 auto trackCollection = trackLinksHandle(*jet);
126 auto trackOriginCollection = trackOriginsHandle(*jet);
127
128 using indexList = std::vector<int>;
129 using trackCountMap = std::map<char, std::set<TL>>;
130
131 indexList iList(vertexCollection.size());
132
133 trackCountMap allTracksMap;
134 trackCountMap heavyFlavourVertexMap;
135 trackCountMap primaryVertexMap;
136 trackCountMap fittingMap;
137
138 fittingMap.clear();
139 heavyFlavourVertexMap.clear();
140 primaryVertexMap.clear();
141 allTracksMap.clear();
142
143
144 for (
int index=0;
index <
int(vertexCollection.size());
index++){
145
147 auto trackOrigin = trackOriginCollection[
index];
148 auto trackLink = trackCollection[
index];
149
150 allTracksMap[
vertex].insert(trackLink);
151
152
155 heavyFlavourVertexMap[
vertex].insert(trackLink);
156 }
157
159 primaryVertexMap[
vertex].insert(trackLink);
160 }
161 };
162
163
164 auto pvCandidate = std::max_element(
165 primaryVertexMap.begin(), primaryVertexMap.end(),
166 [](
const auto&
a,
const auto& b) { return a.second.size() < b.second.size(); }
167 )->first;
168
169
170 for (const auto &vertexTrackPair : allTracksMap) {
171
172 const auto &[
vertex, trackLinks] = vertexTrackPair;
173
175 if(vertex == pvCandidate) {
176 continue;
177 }
178 }
179
181 if (heavyFlavourVertexMap.find(vertex) == heavyFlavourVertexMap.end()) {
182 continue;
183 }
184 }
185
187 if (heavyFlavourVertexMap.find(vertex) != heavyFlavourVertexMap.end() &&
188 (
static_cast<float>(heavyFlavourVertexMap[vertex].size()) / trackLinks.size()) <
m_HFRatioThres) {
189 continue;
190 }
191 }
192
193 fittingMap.insert(std::pair<char, std::set<TL>>(vertex, trackLinks));
194 }
195
197
198 const char mergedVertexKey = 0;
199 std::set<TL> mergedSet;
200 std::set<TL> combinedHFTracks;
201
202
203 for (const auto &vertexTrackPair : fittingMap) {
204
205 mergedSet.insert(vertexTrackPair.second.begin(), vertexTrackPair.second.end());
206 combinedHFTracks.insert(heavyFlavourVertexMap[vertexTrackPair.first].begin(), heavyFlavourVertexMap[vertexTrackPair.first].end());
207 }
208
209
210 fittingMap.clear();
211 fittingMap[mergedVertexKey] = std::move(mergedSet);
212 heavyFlavourVertexMap.clear();
213 heavyFlavourVertexMap[mergedVertexKey] = std::move(combinedHFTracks);
214 }
215
216
217 std::unique_ptr<workVectorArrxAOD> xAODwrk (new workVectorArrxAOD());
218 SG::ReadCondHandle<InDet::BeamSpotData> beamSpotHandle{
m_beamSpotKey, ctx};
219
220
221 xAODwrk->beamX = beamSpotHandle->beamPos().x();
222 xAODwrk->beamY = beamSpotHandle->beamPos().y();
223 xAODwrk->beamZ = beamSpotHandle->beamPos().z();
224 xAODwrk->tanBeamTiltX =
tan(beamSpotHandle->beamTilt(0));
225 xAODwrk->tanBeamTiltY =
tan(beamSpotHandle->beamTilt(1));
226
227 std::unique_ptr<std::vector<WrkVrt>> wrkVrtSet = std::make_unique<std::vector<WrkVrt>>();
231 std::vector<const xAOD::NeutralParticle *> neutralPartDummy(0);
233
234 for (const auto &pair : fittingMap) {
235
236 if (pair.second.size() >= 2) {
237 int NTRKS = pair.second.size();
238 std::vector<double> InpMass(NTRKS,
m_massPi);
240
241 xAODwrk->listSelTracks.clear();
242
244 xAODwrk->listSelTracks.push_back(*
TrackLink);
245 }
246
248 TLorentzVector jetDir(jet->p4());
249
250
252
253 if (
sc.isFailure() || FitVertex.perp() >
m_maxLxy) {
254 IniVrt = primVrt.position();
256 IniVrt.setZero();
257 } else {
258 vDist = FitVertex - primVrt.position();
259 double JetVrtDir = jetDir.Px() * vDist.x() + jetDir.Py() * vDist.y() + jetDir.Pz() * vDist.z();
261 JetVrtDir = fabs(JetVrtDir);
262 if (JetVrtDir > 0.)
263 IniVrt = FitVertex;
264 else
265 IniVrt = primVrt.position();
266 }
267
269
270
271 sc = (
m_vertexFitterTool->VKalVrtFit(xAODwrk->listSelTracks, neutralPartDummy, newvrt.vertex, newvrt.vertexMom,
272 newvrt.vertexCharge, newvrt.vertexCov, newvrt.chi2PerTrk, newvrt.trkAtVrt,
273 newvrt.chi2, *state, false));
275 continue;
276
277
278 auto NDOF = 2 * (newvrt.trkAtVrt.size()) - 3.0;
279
281 continue;
282 ATH_MSG_DEBUG(
"Found IniVertex=" << newvrt.vertex[0] <<
", " << newvrt.vertex[1] <<
", " << newvrt.vertex[2]
283 << " trks " << newvrt.trkAtVrt.size());
284
286
287 Amg::Vector3D jetVrtDir(jet->p4().Px(), jet->p4().Py(), jet->p4().Pz());
288
289 double vPos =
290 (vDir.x() * newvrt.vertexMom.Px() + vDir.y() * newvrt.vertexMom.Py() + vDir.z() * newvrt.vertexMom.Pz()) /
291 newvrt.vertexMom.Rho();
292
293 double Lxy = sqrt(vDir[0] * vDir[0] + vDir[1] * vDir[1]);
294 double Lxyz = sqrt(vDir[0] * vDir[0] + vDir[1] * vDir[1] + vDir[2] * vDir[2]);
296
298
299 int ntrk = newvrt.trkAtVrt.size();
300
301 TLorentzVector MomentumVtx =
TotalMom(xAODwrk->listSelTracks);
302
303 double eRatio = MomentumVtx.E() / jet->p4().E();
304 double signif3D;
305 [[maybe_unused]]
double distToPV =
vrtVrtDist(primVrt, newvrt.vertex, newvrt.vertexCov, signif3D);
306
307
309
311 continue;
312
314 continue;
315
317 continue;
318 }
319
320
322
323 for (const auto *trk : xAODwrk->listSelTracks) {
324 ElementLink<xAOD::TrackParticleContainer> link_trk(
326 trk->index());
327 GNNvertex->addTrackAtVertex(link_trk, 1.);
328 }
329
330
332 GNNvertex->setPosition(newvrt.vertex);
333 GNNvertex->setFitQuality(newvrt.chi2, NDOF);
334
336 m_deco_pt(*GNNvertex) = newvrt.vertexMom.Perp();
346
347 ElementLink<xAOD::VertexContainer> linkVertex;
350 jetWriteDecorHandleVertexLink(*jet).push_back(linkVertex);
351
352 }
353 }
354 }
355 return StatusCode::SUCCESS;
356}
ElementLink< xAOD::TrackParticleContainer > TrackLink
value_type emplace_back(value_type pElem)
Add an element to the end of the collection.
bool setElement(ElementType element)
Set to point to an element.
bool setStorableObject(BaseConstReference data, bool replace=false, IProxyDict *sg=0)
Set link to point to a new container (storable).
void setVertexType(VxType::VertexType vType)
Set the type of the vertex.
Eigen::Matrix< double, 3, 1 > Vector3D
double deltaR(const Amg::Vector3D &v1, const Amg::Vector3D &v2)
::StatusCode StatusCode
StatusCode definition for legacy code.
@ SecVtx
Secondary vertex.
Jet_v1 Jet
Definition of the current "jet version".
Vertex_v1 Vertex
Define the latest version of the vertex class.
TrackParticleContainer_v1 TrackParticleContainer
Definition of the current "TrackParticle container version".
1.14.0