ATLAS Offline Software
Loading...
Searching...
No Matches
VertexingAlgs.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
3*/
4
5// Header include
6#include "VrtSecInclusive/VrtSecInclusive.h"
7#include "VrtSecInclusive/NtupleVars.h"
8#include "VrtSecInclusive/Tools.h"
9#include "AthContainers/Accessor.h"
10
11#include "TrkTrackSummary/TrackSummary.h"
12
13#include "EventPrimitives/EventPrimitivesHelpers.h"
14
15#include "xAODEgamma/ElectronxAODHelpers.h"
16
17#include "TH1F.h"
18#include "TH2F.h"
19#include "TNtuple.h"
20#include "TTree.h"
21#include "TROOT.h"
22#include "TLorentzVector.h"
23
24#include <iostream>
25#include "TrkVKalVrtCore/PGraph.h"
26#include <algorithm>
27#include <array>
28
29//-------------------------------------------------
30
31using namespace std;
32
33
34
35namespace VKalVrtAthena {
36
37
38 //____________________________________________________________________________________________________
39 StatusCode VrtSecInclusive::extractIncompatibleTrackPairs( const EventContext& ctx,
40 std::vector<WrkVrt>* workVerticesContainer )
41 {
42
43 // Output SVs as xAOD::Vertex
44 // Needs a conversion function from WrkVrtSet to xAOD::Vertex here.
45 // The supposed form of the function will be as follows:
46
47 SG::WriteHandle<xAOD::VertexContainer> trackHandle;
48 xAOD::VertexContainer *twoTrksVertexContainer{};
49 if( m_FillIntermediateVertices ) {
50 trackHandle = SG::makeHandle( m_twoTrksVertexKey, ctx );
51 ATH_CHECK( trackHandle.record(std::make_unique<xAOD::VertexContainer>(),
52 std::make_unique<xAOD::VertexAuxContainer>()) );
53 twoTrksVertexContainer = trackHandle.ptr();
54 m_vertexCollectionsDefinitions[m_twoTrksVertexKey.key()] = true;
55 }
56
57 m_incomp.clear();
58
59 // Work variables
60 std::vector<const xAOD::TrackParticle*> baseTracks;
61 std::vector<const xAOD::NeutralParticle*> dummyNeutrals;
62
63 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Selected Tracks = "<< m_selectedTracks.size());
64 if( m_FillHist ) { m_hists["selTracksDist"]->Fill( m_selectedTracks.size() ); }
65
66 std::string msg;
67
68 enum recoStep { kStart, kInitVtxPosition, kImpactParamCheck, kVKalVrtFit, kChi2, kVposCut, kPatternMatch };
69
70 const double maxR { 563. }; // r = 563 mm is the TRT inner surface
71 double roughD0Cut = 100.;
72 double roughZ0Cut = 50.;
73 if(m_doDisappearingTrackVertexing){
74 roughD0Cut = 1000.;
75 roughZ0Cut = 1000.;
76 }
77
78 // Truth match map
79 std::map<const xAOD::TruthVertex*, bool> matchMap;
80 std::unique_ptr<Trk::IVKalState> state = m_fitSvc->makeState();
81 // first make all 2-track vertices
82 for( auto itrk = m_selectedTracks.begin(); itrk != m_selectedTracks.end(); ++itrk ) {
83 for( auto jtrk = std::next(itrk); jtrk != m_selectedTracks.end(); ++jtrk ) {
84
85 // avoid both tracks are too close to the beam line
86
87 const int itrk_id = itrk - m_selectedTracks.begin();
88 const int jtrk_id = jtrk - m_selectedTracks.begin();
89
90 WrkVrt wrkvrt;
91 wrkvrt.selectedTrackIndices.emplace_back( itrk_id );
92 wrkvrt.selectedTrackIndices.emplace_back( jtrk_id );
93
94 // Attempt to think the combination is incompatible by default
95 m_incomp.emplace_back( itrk_id, jtrk_id );
96
97 if(m_doDisappearingTrackVertexing) {
98
99 const auto* cont_i = dynamic_cast<const xAOD::TrackParticleContainer*>( (*itrk)->container() );
100 const auto* cont_j = dynamic_cast<const xAOD::TrackParticleContainer*>( (*jtrk)->container() );
101
102 if ( !cont_i || !cont_j ) {
103 ATH_MSG_DEBUG(" one of the track containers is null");
104 continue;
105 }
106
107 ElementLink<xAOD::TrackParticleContainer> link_i, link_j;
108 link_i.toIndexedElement( *cont_i, (*itrk)->index() );
109 link_j.toIndexedElement( *cont_j, (*jtrk)->index() );
110
111 if (!link_i.isValid() || !link_j.isValid()) {
112 ATH_MSG_DEBUG(" link itrk (" << (*itrk)->index() << ") or jtrk (" << (*jtrk)->index() << ") is not valid");
113 }
114 else {
115 if( link_i.dataID() == link_j.dataID() ) {
116 continue;
117 }
118 }
119 }
120
121
122 if( std::abs( (*itrk)->d0() ) < m_twoTrkVtxFormingD0Cut && std::abs( (*jtrk)->d0() ) < m_twoTrkVtxFormingD0Cut ) continue;
123
124 baseTracks.clear();
125 baseTracks.emplace_back( *itrk );
126 baseTracks.emplace_back( *jtrk );
127
128 if( m_FillHist ) m_hists["incompMonitor"]->Fill( kStart );
129
130 // new code to find initial approximate vertex
131 Amg::Vector3D initVertex;
132
133 StatusCode sc = m_fitSvc->VKalVrtFitFast( baseTracks, initVertex, *state );/* Fast crude estimation */
134 if( sc.isFailure() ) {
135 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": fast crude estimation fails ");
136 continue;
137 }
138
139 if( initVertex.perp() > maxR ) {
140 continue;
141 }
142 if( m_doDisappearingTrackVertexing && initVertex.perp() <m_twoTrVrtMinRadius){
143 continue;
144 }
145 if( m_FillHist ) m_hists["incompMonitor"]->Fill( kInitVtxPosition );
146
147 std::vector<double> impactParameters;
148 std::vector<double> impactParErrors;
149
150 if( !getSVImpactParameters( *itrk, initVertex, impactParameters, impactParErrors) ) continue;
151 const auto roughD0_itrk = impactParameters.at(TrkParameter::k_d0);
152 const auto roughZ0_itrk = impactParameters.at(TrkParameter::k_z0);
153 if( fabs( impactParameters.at(0)) > roughD0Cut || fabs( impactParameters.at(1) ) > roughZ0Cut ) {
154 continue;
155 }
156
157 if( !getSVImpactParameters( *jtrk, initVertex, impactParameters, impactParErrors) ) continue;
158 const auto roughD0_jtrk = impactParameters.at(TrkParameter::k_d0);
159 const auto roughZ0_jtrk = impactParameters.at(TrkParameter::k_z0);
160 if( fabs( impactParameters.at(0) ) > roughD0Cut || fabs( impactParameters.at(1) ) > roughZ0Cut ) {
161 continue;
162 }
163 if( m_FillHist ) m_hists["incompMonitor"]->Fill( kImpactParamCheck );
164
165 m_fitSvc->setApproximateVertex( initVertex.x(), initVertex.y(), initVertex.z(), *state );
166
167
168
169 // Vertex VKal Fitting
170 sc = m_fitSvc->VKalVrtFit( baseTracks,
171 dummyNeutrals,
172 wrkvrt.vertex, wrkvrt.vertexMom, wrkvrt.Charge,
173 wrkvrt.vertexCov, wrkvrt.Chi2PerTrk,
174 wrkvrt.TrkAtVrt, wrkvrt.Chi2, *state );
175
176 if( sc.isFailure() ) {
177 continue; /* No fit */
178 }
179 if( m_FillHist ) m_hists["incompMonitor"]->Fill( kVKalVrtFit );
180
181 // Compatibility to the primary vertex.
182 Amg::Vector3D vDist = wrkvrt.vertex - m_thePV->position();
183 const double vPos = ( vDist.x()*wrkvrt.vertexMom.Px()+vDist.y()*wrkvrt.vertexMom.Py()+vDist.z()*wrkvrt.vertexMom.Pz() )/wrkvrt.vertexMom.Rho();
184 const double vPosMomAngT = ( vDist.x()*wrkvrt.vertexMom.Px()+vDist.y()*wrkvrt.vertexMom.Py() ) / vDist.perp() / wrkvrt.vertexMom.Pt();
185 const double vPosMomAng3D = ( vDist.x()*wrkvrt.vertexMom.Px()+vDist.y()*wrkvrt.vertexMom.Py()+vDist.z()*wrkvrt.vertexMom.Pz() ) / (vDist.norm() * wrkvrt.vertexMom.Rho());
186
187 double dphi1 = TVector2::Phi_mpi_pi(vDist.phi() - (*itrk)->phi());
188 double dphi2 = TVector2::Phi_mpi_pi(vDist.phi() - (*jtrk)->phi());
189
190 const double dist_fromPV = vDist.norm();
191 if( m_FillHist ) m_hists["2trkVtxDistFromPV"]->Fill( dist_fromPV );
192
193 if( m_FillNtuple ) {
194 // Fill the 2-track vertex properties to AANT
195 m_ntupleVars->get<unsigned int>( "All2TrkVrtNum" )++;
196 m_ntupleVars->get< std::vector<double> >( "All2TrkVrtMass" ) .emplace_back(wrkvrt.vertexMom.M());
197 m_ntupleVars->get< std::vector<double> >( "All2TrkVrtPt" ) .emplace_back(wrkvrt.vertexMom.Perp());
198 m_ntupleVars->get< std::vector<int> > ( "All2TrkVrtCharge" ) .emplace_back(wrkvrt.Charge);
199 m_ntupleVars->get< std::vector<double> >( "All2TrkVrtX" ) .emplace_back(wrkvrt.vertex.x());
200 m_ntupleVars->get< std::vector<double> >( "All2TrkVrtY" ) .emplace_back(wrkvrt.vertex.y());
201 m_ntupleVars->get< std::vector<double> >( "All2TrkVrtZ" ) .emplace_back(wrkvrt.vertex.z());
202 m_ntupleVars->get< std::vector<double> >( "All2TrkVrtChiSq" ) .emplace_back(wrkvrt.Chi2);
203 }
204
205
206 // Create a xAOD::Vertex instance
207 xAOD::Vertex *vertex{};
208
209 if( m_FillIntermediateVertices ) {
210 vertex = new xAOD::Vertex;
211 twoTrksVertexContainer->emplace_back( vertex );
212
213 for( const auto *trk: baseTracks ) {
214
215 // Acquire link to the track
216 ElementLink<xAOD::TrackParticleContainer> trackElementLink( *( dynamic_cast<const xAOD::TrackParticleContainer*>( trk->container() ) ), trk->index() );
217
218 // Register link to the vertex
219 vertex->addTrackAtVertex( trackElementLink, 1. );
220 }
221
222 vertex->setVertexType( xAOD::VxType::SecVtx );
223 vertex->setPosition( wrkvrt.vertex );
224 vertex->setFitQuality( wrkvrt.Chi2, 1 ); // Ndof is always 1
225
226 static const SG::Accessor<float> massAcc("mass");
227 static const SG::Accessor<float> pTAcc("pT");
228 static const SG::Accessor<float> chargeAcc("charge");
229 static const SG::Accessor<float> vPosAcc("vPos");
230 static const SG::Accessor<bool> isFakeAcc("isFake");
231 massAcc(*vertex) = wrkvrt.vertexMom.M();
232 pTAcc(*vertex) = wrkvrt.vertexMom.Perp();
233 chargeAcc(*vertex) = wrkvrt.Charge;
234 vPosAcc(*vertex) = vPos;
235 isFakeAcc(*vertex) = true;
236 }
237
238
240
241 uint8_t trkiBLHit,trkjBLHit;
242 if( !((*itrk)->summaryValue( trkiBLHit,xAOD::numberOfInnermostPixelLayerHits))) trkiBLHit=0;
243 if( !((*jtrk)->summaryValue( trkjBLHit,xAOD::numberOfInnermostPixelLayerHits))) trkjBLHit=0;
244
245 if( m_FillNtuple ) m_ntupleVars->get< std::vector<int> >( "All2TrkSumBLHits" ).emplace_back( trkiBLHit + trkjBLHit );
246
247 // track chi2 cut
248 if( m_FillHist ) m_hists["2trkChi2Dist"]->Fill( log10( wrkvrt.Chi2 ) );
249
250 if( wrkvrt.fitQuality() > m_SelVrtChi2Cut) {
251 ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": failed to pass chi2 threshold." );
252 continue; /* Bad Chi2 */
253 }
254 if( m_FillHist ) m_hists["incompMonitor"]->Fill( kChi2 );
255
256
257 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": attempting form vertex from ( " << itrk_id << ", " << jtrk_id << " )." );
258 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": candidate vertex: "
259 << " isGood = " << (wrkvrt.isGood? "true" : "false")
260 << ", #ntrks = " << wrkvrt.nTracksTotal()
261 << ", #selectedTracks = " << wrkvrt.selectedTrackIndices.size()
262 << ", #associatedTracks = " << wrkvrt.associatedTrackIndices.size()
263 << ", chi2/ndof = " << wrkvrt.fitQuality()
264 << ", (r, z) = (" << wrkvrt.vertex.perp()
265 <<", " << wrkvrt.vertex.z() << ")" );
266
267 for( const auto* truthVertex : m_tracingTruthVertices ) {
268 Amg::Vector3D vTruth( truthVertex->x(), truthVertex->y(), truthVertex->z() );
269 Amg::Vector3D vReco ( wrkvrt.vertex.x(), wrkvrt.vertex.y(), wrkvrt.vertex.z() );
270
271 const auto distance = vReco - vTruth;
272
273 AmgSymMatrix(3) cov;
274 cov.fillSymmetric( 0, 0, wrkvrt.vertexCov.at(0) );
275 cov.fillSymmetric( 1, 0, wrkvrt.vertexCov.at(1) );
276 cov.fillSymmetric( 1, 1, wrkvrt.vertexCov.at(2) );
277 cov.fillSymmetric( 2, 0, wrkvrt.vertexCov.at(3) );
278 cov.fillSymmetric( 2, 1, wrkvrt.vertexCov.at(4) );
279 cov.fillSymmetric( 2, 2, wrkvrt.vertexCov.at(5) );
280
281 const double s2 = distance.transpose() * cov.inverse() * distance;
282
283 if( distance.norm() < 2.0 || s2 < 100. ) {
284 ATH_MSG_DEBUG ( " > " << __FUNCTION__ << ": truth-matched candidate! : signif^2 = " << s2 );
285 matchMap.emplace( truthVertex, true );
286 }
287 }
288
289 if( m_FillHist ) {
290 dynamic_cast<TH2F*>( m_hists["vPosDist"] )->Fill( wrkvrt.vertex.perp(), vPos );
291 dynamic_cast<TH2F*>( m_hists["vPosMomAngTDist"] )->Fill( wrkvrt.vertex.perp(), vPosMomAngT );
292 m_hists["vPosMomAngT"] ->Fill( vPosMomAngT );
293 m_hists["vPosMomAng3D"] ->Fill( vPosMomAng3D );
294 }
295
296 if( m_doTwoTrSoftBtag ){
297 if(dist_fromPV < m_twoTrVrtMinDistFromPV ){
298 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": failed to pass the 2tr vertex min distance from PV cut." );
299 continue;
300 }
301
302 if( vPosMomAng3D < m_twoTrVrtAngleCut ){
303 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": failed to pass the vertex angle cut." );
304 continue;
305 }
306 }
307
308 if( m_doPVcompatibilityCut ) {
309 if( cos( dphi1 ) < -0.8 && cos( dphi2 ) < -0.8 ) {
310 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": failed to pass the vPos cut. (both tracks are opposite against the vertex pos)" );
311 continue;
312 }
313 if (m_doTightPVcompatibilityCut && (cos( dphi1 ) < -0.8 || cos( dphi2 ) < -0.8)){
314 ATH_MSG_DEBUG(" > "<< __FUNCTION__ << ": failed to pass the tightened vPos cut. (at least one track is opposite against the vertex pos)" );
315 continue;
316 }
317 if( vPosMomAngT < -0.8 ) {
318 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": failed to pass the vPos cut. (pos-mom directions are opposite)" );
319 continue;
320 }
321 if( vPos < m_pvCompatibilityCut ) {
322 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": failed to pass the vPos cut." );
323 continue;
324 }
325 }
326 if( m_FillHist ) m_hists["incompMonitor"]->Fill( kVposCut );
327
328 // fake rejection cuts with track hit pattern consistencies
329 if( m_removeFakeVrt && !m_removeFakeVrtLate ) {
330 if( !this->passedFakeReject( wrkvrt.vertex, (*itrk), (*jtrk) ) ) {
331
332 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": failed to pass fake rejection algorithm." );
333 continue;
334 }
335 }
336 if( m_FillHist ) m_hists["incompMonitor"]->Fill( kPatternMatch );
337
338 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": passed fake rejection." );
339
340 if( m_FillNtuple ) {
341 // Fill AANT for vertices after fake rejection
342 m_ntupleVars->get< unsigned int >( "AfFakVrtNum" )++;
343 m_ntupleVars->get< std::vector<double> >( "AfFakVrtMass" ) .emplace_back(wrkvrt.vertexMom.M());
344 m_ntupleVars->get< std::vector<double> >( "AfFakVrtPt" ) .emplace_back(wrkvrt.vertexMom.Perp());
345 m_ntupleVars->get< std::vector<int> > ( "AfFakVrtCharge" ) .emplace_back(wrkvrt.Charge);
346 m_ntupleVars->get< std::vector<double> >( "AfFakVrtX" ) .emplace_back(wrkvrt.vertex.x());
347 m_ntupleVars->get< std::vector<double> >( "AfFakVrtY" ) .emplace_back(wrkvrt.vertex.y());
348 m_ntupleVars->get< std::vector<double> >( "AfFakVrtZ" ) .emplace_back(wrkvrt.vertex.z());
349 m_ntupleVars->get< std::vector<double> >( "AfFakVrtChiSq" ) .emplace_back(wrkvrt.Chi2);
350 }
351
352 // The vertex passed the quality cut: overwrite isFake to false
353 if( m_FillIntermediateVertices && vertex ) {
354 static const SG::Accessor<bool> isFakeAcc("isFake");
355 isFakeAcc(*vertex) = false;
356 }
357
358
359 // Now this vertex passed all criteria and considred to be a compatible vertices.
360 // Therefore the track pair is removed from the incompatibility list.
361 m_incomp.pop_back();
362
363 wrkvrt.isGood = true;
364
365 workVerticesContainer->emplace_back( wrkvrt );
366
367 msg += Form(" (%d, %d), ", itrk_id, jtrk_id );
368
369 if( m_FillHist ) {
370 m_hists["initVertexDispD0"]->Fill( roughD0_itrk, initVertex.perp() );
371 m_hists["initVertexDispD0"]->Fill( roughD0_jtrk, initVertex.perp() );
372 m_hists["initVertexDispZ0"]->Fill( roughZ0_itrk, initVertex.z() );
373 m_hists["initVertexDispZ0"]->Fill( roughZ0_jtrk, initVertex.z() );
374 }
375
376 }
377 }
378
379
380 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": compatible track pairs = " << msg );
381
382 if( m_FillNtuple ) m_ntupleVars->get<unsigned int>( "SizeIncomp" ) = m_incomp.size();
383
384 if( m_FillHist ) {
385 for( auto& pair: matchMap ) {
386 if( pair.second ) m_hists["nMatchedTruths"]->Fill( 1, pair.first->perp() );
387 }
388 }
389
390 return StatusCode::SUCCESS;
391 }
392
393
394 //____________________________________________________________________________________________________
395 StatusCode VrtSecInclusive::findNtrackVertices( const EventContext&,
396 std::vector<WrkVrt> *workVerticesContainer )
397 {
398 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": begin");
399 if(m_doDisappearingTrackVertexing){
400 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": skip");
401 return StatusCode::SUCCESS;
402 }
403
404
405 const auto compSize = m_selectedTracks.size()*(m_selectedTracks.size() - 1)/2 - m_incomp.size();
406 if( m_FillHist ) { m_hists["2trkVerticesDist"]->Fill( compSize ); }
407
408 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": compatible track pair size = " << compSize );
409 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": incompatible track pair size = " << m_incomp.size() );
410
411
412 if( not m_doFastMode ) {
413
414 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": incompatibility graph finder mode" );
415
416 // clear the container
417 workVerticesContainer->clear();
418
419 // Graph method: Trk::pgraphm_()
420 // used in order to find compatible sub-graphs from the incompatible graph
421
422 // List of edgeds between imcompatible nodes
423 // This weit is the data model of imcompatible graph used in Trk::pgraphm_().
424 std::vector<long int> weit;
425
426 for( auto& pair : m_incomp ) {
427 weit.emplace_back( pair.first + 1 ); /* +1 is needed for PGRAPH due to FORTRAN-style counting */
428 weit.emplace_back( pair.second + 1 ); /* +1 is needed for PGRAPH due to FORTRAN-style counting */
429 }
430
431 // Solution of the graph method routine (minimal covering of the graph)
432 // The size of the solution is returned by NPTR (see below)
433 std::vector<long int> solution( m_selectedTracks.size() );
434
435 // Number of edges in the list is the size of incompatibility track pairs.
436 long int nEdges = m_incomp.size();
437
438 // input number of nodes in the graph.
439 long int nTracks = static_cast<long int>( m_selectedTracks.size() );
440
441 // Input variable; the threshold. Solutions shorter than nth are not returned (ignored).
442 long int nth = 2; //VK some speed up
443
444 // NPTR: I/O variable (Destructive FORTRAN Style!!!)
445 // - on input: =0 for initialization, >0 to get next solution
446 // - on output: >0 : length of the solution stored in set; =0 : no more solutions can be found
447 long int solutionSize { 0 };
448
449 // This is just a unused strawman needed for m_fitSvc->VKalVrtFit()
450 std::vector<const xAOD::TrackParticle*> baseTracks;
451 std::vector<const xAOD::NeutralParticle*> dummyNeutrals;
452
453 std::unique_ptr<Trk::IVKalState> state = m_fitSvc->makeState();
454 auto pgraph = std::make_unique<Trk::PGraph>();
455 int iterationLimit(2000);
456 // Main iteration
457 while(true) {
458
459 // Find a solution from the given set of incompatible tracks (==weit)
460 pgraph->pgraphm_( weit.data(), nEdges, nTracks, solution.data(), &solutionSize, nth);
461
462 ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": Trk::pgraphm_() output: solutionSize = " << solutionSize );
463 if (0 == iterationLimit--){
464 ATH_MSG_WARNING("Iteration limit (2000) reached in VrtSecInclusive::findNtrackVertices, solution size = "<<solutionSize);
465 break;
466 }
467 if(solutionSize <= 0) break; // No more solutions ==> Exit
468 if(solutionSize == 1) continue; // i.e. single node ==> Not a good solution
469
470 baseTracks.clear();
471
472 std::string msg = "solution = [ ";
473 for( int i=0; i< solutionSize; i++) {
474 msg += Form( "%ld, ", solution[i]-1 );
475 }
476 msg += " ]";
477 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": " << msg );
478
479 // varaible of new vertex
480 WrkVrt wrkvrt;
481
482 // Try to compose a new vertex using the solution nodes
483 // Here the track ID is labelled with array
484 wrkvrt.isGood = true;
485 wrkvrt.selectedTrackIndices.clear();
486
487 for(long int i = 0; i<solutionSize; i++) {
488 wrkvrt.selectedTrackIndices.emplace_back(solution[i]-1);
489 baseTracks.emplace_back( m_selectedTracks.at(solution[i]-1) );
490 }
491
492 // Perform vertex fitting
493 Amg::Vector3D initVertex;
494
495 StatusCode sc = m_fitSvc->VKalVrtFitFast( baseTracks, initVertex, *state );/* Fast crude estimation */
496 if(sc.isFailure()) ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": fast crude estimation fails ");
497
498 m_fitSvc->setApproximateVertex( initVertex.x(), initVertex.y(), initVertex.z(), *state );
499
500 sc = m_fitSvc->VKalVrtFit(baseTracks, dummyNeutrals,
501 wrkvrt.vertex,
502 wrkvrt.vertexMom,
503 wrkvrt.Charge,
504 wrkvrt.vertexCov,
505 wrkvrt.Chi2PerTrk,
506 wrkvrt.TrkAtVrt,
507 wrkvrt.Chi2,
508 *state);
509
510 ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": FoundAppVrt=" << solutionSize << ", (r, z) = " << wrkvrt.vertex.perp() << ", " << wrkvrt.vertex.z() << ", chi2/ndof = " << wrkvrt.fitQuality() );
511
512 if( sc.isFailure() ) {
513
514 if( wrkvrt.selectedTrackIndices.size() <= 2 ) {
515 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": VKalVrtFit failed in 2-trk solution ==> give up.");
516 continue;
517 }
518
519 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": VKalVrtFit failed ==> retry...");
520
521 WrkVrt tmp;
522 tmp.isGood = false;
523
524 // Create 2-trk vertex combination and find any compatible vertex
525 for( auto& itrk: wrkvrt.selectedTrackIndices ) {
526 for( auto& jtrk: wrkvrt.selectedTrackIndices ) {
527 if( itrk == jtrk ) continue;
528 if( tmp.isGood ) continue;
529
530 tmp.selectedTrackIndices.clear();
531 tmp.selectedTrackIndices.emplace_back( itrk );
532 tmp.selectedTrackIndices.emplace_back( jtrk );
533
534 baseTracks.clear();
535 baseTracks.emplace_back( m_selectedTracks.at( itrk ) );
536 baseTracks.emplace_back( m_selectedTracks.at( jtrk ) );
537
538 // Perform vertex fitting
539 Amg::Vector3D initVertex;
540
541 sc = m_fitSvc->VKalVrtFitFast( baseTracks, initVertex, *state );
542 if( sc.isFailure() ) ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": fast crude estimation fails ");
543
544 m_fitSvc->setApproximateVertex( initVertex.x(), initVertex.y(), initVertex.z(), *state );
545
546 sc = m_fitSvc->VKalVrtFit(baseTracks, dummyNeutrals,
547 tmp.vertex,
548 tmp.vertexMom,
549 tmp.Charge,
550 tmp.vertexCov,
551 tmp.Chi2PerTrk,
552 tmp.TrkAtVrt,
553 tmp.Chi2,
554 *state);
555
556 if( sc.isFailure() ) continue;
557
558 tmp.isGood = true;
559
560 }
561 }
562
563 if( !tmp.isGood ) {
564 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Did not find any viable vertex in all 2-trk combinations. Give up.");
565 continue;
566 }
567
568 // Now, found at least one seed 2-track vertex. ==> attempt to attach other tracks
569 for( auto& itrk: wrkvrt.selectedTrackIndices ) {
570
571 if( std::find( tmp.selectedTrackIndices.begin(), tmp.selectedTrackIndices.end(), itrk ) != tmp.selectedTrackIndices.end() ) continue;
572
573 auto backup = tmp;
574
575 tmp.selectedTrackIndices.emplace_back( itrk );
576 baseTracks.clear();
577 for( auto& jtrk : tmp.selectedTrackIndices ) { baseTracks.emplace_back( m_selectedTracks.at(jtrk) ); }
578
579 // Perform vertex fitting
580 Amg::Vector3D initVertex;
581
582 sc = m_fitSvc->VKalVrtFitFast( baseTracks, initVertex, *state );/* Fast crude estimation */
583 if(sc.isFailure()) ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": fast crude estimation fails ");
584
585 m_fitSvc->setApproximateVertex( initVertex.x(), initVertex.y(), initVertex.z(), *state );
586
587 sc = m_fitSvc->VKalVrtFit(baseTracks, dummyNeutrals,
588 tmp.vertex,
589 tmp.vertexMom,
590 tmp.Charge,
591 tmp.vertexCov,
592 tmp.Chi2PerTrk,
593 tmp.TrkAtVrt,
594 tmp.Chi2,
595 *state);
596
597 if( sc.isFailure() ) {
598 tmp = backup;
599 continue;
600 }
601
602 }
603
604 wrkvrt = tmp;
605 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": VKalVrtFit succeeded; register the vertex to the list.");
606 wrkvrt.isGood = true;
607 wrkvrt.closestWrkVrtIndex = AlgConsts::invalidUnsigned;
608 wrkvrt.closestWrkVrtValue = AlgConsts::maxValue;
609 workVerticesContainer->emplace_back( wrkvrt );
610
611 } else {
612
613 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": VKalVrtFit succeeded; register the vertex to the list.");
614 wrkvrt.isGood = true;
615 wrkvrt.closestWrkVrtIndex = AlgConsts::invalidUnsigned;
616 wrkvrt.closestWrkVrtValue = AlgConsts::maxValue;
617 workVerticesContainer->emplace_back( wrkvrt );
618
619 }
620
621 }
622
623
624 } else {
625
626 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": rapid finder mode" );
627
628 struct Cluster {
629 Amg::Vector3D position;
630 std::set<long int> tracks;
631 };
632
633 std::vector<struct Cluster> clusters;
634
635 for( auto& wrkvrt : *workVerticesContainer ) {
636
637 bool foundCluster = false;
638
639 for( auto& cluster: clusters ) {
640 if( (wrkvrt.vertex - cluster.position).norm() < 1.0 ) {
641 for( auto& itrk : wrkvrt.selectedTrackIndices ) {
642 cluster.tracks.insert( itrk );
643 }
644 foundCluster = true;
645 break;
646 }
647 }
648
649 if( !foundCluster ) {
650 Cluster c;
651 c.position = wrkvrt.vertex;
652 for( auto& itrk : wrkvrt.selectedTrackIndices ) {
653 c.tracks.insert( itrk );
654 }
655 clusters.emplace_back( c );
656 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": added a new cluster" );
657 }
658
659 }
660
661 // This is just a unused strawman needed for m_fitSvc->VKalVrtFit()
662 std::vector<const xAOD::TrackParticle*> baseTracks;
663 std::vector<const xAOD::NeutralParticle*> dummyNeutrals;
664
665 workVerticesContainer->clear();
666
667 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": found cluster size =" << clusters.size() );
668
669 std::unique_ptr<Trk::IVKalState> state = m_fitSvc->makeState();
670 for( auto& cluster : clusters ) {
671
672 // varaible of new vertex
673 WrkVrt wrkvrt;
674
675 // Try to compose a new vertex using the solution nodes
676 // Here the track ID is labelled with array
677 wrkvrt.isGood = true;
678 wrkvrt.selectedTrackIndices.clear();
679
680 for(const auto& index: cluster.tracks) {
681 wrkvrt.selectedTrackIndices.emplace_back( index );
682 baseTracks.emplace_back( m_selectedTracks.at( index ) );
683 }
684
685 // Perform vertex fitting
686 Amg::Vector3D initVertex;
687
688 StatusCode sc = m_fitSvc->VKalVrtFitFast( baseTracks, initVertex, *state );/* Fast crude estimation */
689 if(sc.isFailure()) ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": fast crude estimation fails ");
690
691 m_fitSvc->setApproximateVertex( initVertex.x(), initVertex.y(), initVertex.z(), *state );
692
693 sc = m_fitSvc->VKalVrtFit(baseTracks, dummyNeutrals,
694 wrkvrt.vertex,
695 wrkvrt.vertexMom,
696 wrkvrt.Charge,
697 wrkvrt.vertexCov,
698 wrkvrt.Chi2PerTrk,
699 wrkvrt.TrkAtVrt,
700 wrkvrt.Chi2,
701 *state);
702
703 if( sc.isFailure() ) {
704 continue;
705 }
706
707 workVerticesContainer->emplace_back( wrkvrt );
708 }
709
710 }
711
712 if (m_truncateWrkVertices){
713 if (workVerticesContainer->size() > m_maxWrkVertices){
714 m_vertexingStatus = 3;
715 workVerticesContainer->resize(m_maxWrkVertices);
716 }
717 }
718
719 //-------------------------------------------------------
720 // Iterative cleanup algorithm
721
722 //-Remove vertices fully contained in other vertices
723 ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": Remove vertices fully contained in other vertices .");
724 while( workVerticesContainer->size() > 1 ) {
725 size_t tmpN = workVerticesContainer->size();
726
727 size_t iv = 0;
728 for(; iv<tmpN-1; iv++) {
729 size_t jv = iv+1;
730 for(; jv<tmpN; jv++) {
731 const auto nTCom = nTrkCommon( workVerticesContainer, {iv, jv} );
732
733 if( nTCom == workVerticesContainer->at(iv).selectedTrackIndices.size() ) { workVerticesContainer->erase(workVerticesContainer->begin()+iv); break; }
734 else if( nTCom == workVerticesContainer->at(jv).selectedTrackIndices.size() ) { workVerticesContainer->erase(workVerticesContainer->begin()+jv); break; }
735
736 }
737 if(jv!=tmpN) break; // One vertex is erased. Restart check
738 }
739 if(iv==tmpN-1) break; // No vertex deleted
740 }
741
742 //-Identify remaining 2-track vertices with very bad Chi2 and mass (b-tagging)
743 ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": Identify remaining 2-track vertices with very bad Chi2 and mass (b-tagging).");
744 for( auto& wrkvrt : *workVerticesContainer ) {
745
746 if( TMath::Prob( wrkvrt.Chi2, wrkvrt.ndof() ) < m_improveChi2ProbThreshold ) wrkvrt.isGood = false;
747 if( wrkvrt.selectedTrackIndices.size() != 2 ) continue;
748 if( m_FillHist ) m_hists["NtrkChi2Dist"]->Fill( log10( wrkvrt.fitQuality() ) );
749 }
750
751 if( m_FillNtuple) m_ntupleVars->get<unsigned int>( "NumInitSecVrt" ) = workVerticesContainer->size();
752
753 return StatusCode::SUCCESS;
754 }
755
756
757 //____________________________________________________________________________________________________
758 StatusCode VrtSecInclusive::rearrangeTracks( const EventContext&,
759 std::vector<WrkVrt> *workVerticesContainer )
760 {
761 if(m_doDisappearingTrackVertexing){
762 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": skip");
763 return StatusCode::SUCCESS;
764 }
765 //
766 // Rearrangement of solutions
767 //
768
769 std::vector<long int> processedTracks;
770
771 unsigned mergeCounter { 0 };
772 unsigned brokenCounter { 0 };
773 unsigned removeTrackCounter { 0 };
774
775 while( true ) {
776
777 // worstChi2: unit in [chi2 per track]
778 long int maxSharedTrack;
779 long int worstMatchingVertex;
780 std::pair<unsigned, unsigned> indexPair { AlgConsts::invalidUnsigned, AlgConsts::invalidUnsigned };
781
782
783 // trackToVertexMap has IDs of each track which can contain array of vertices.
784 // e.g. TrkInVrt->at( track_id ).size() gives the number of vertices which use the track [track_id].
785
786 std::map<long int, std::vector<long int> > trackToVertexMap;
787
788 // Fill trackToVertexMap with vertex IDs of each track
789 trackClassification( workVerticesContainer, trackToVertexMap );
790
791
792 auto worstChi2 = findWorstChi2ofMaximallySharedTrack( workVerticesContainer, trackToVertexMap, maxSharedTrack, worstMatchingVertex );
793
794 if( worstChi2 == AlgConsts::invalidFloat ) {
795 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": no shared tracks are found --> exit the while loop." );
796 break;
797 }
798
799 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": vertex [" << worstMatchingVertex << "]: maximally shared track index = " << maxSharedTrack
800 << ", multiplicity = " << trackToVertexMap.at( maxSharedTrack ).size()
801 << ", worst chi2_trk = " << worstChi2 );
802
803 //Choice of action
804 if( worstChi2 < m_TrackDetachCut ) {
805
806 // Here, the max-shared track is well-associated and cannot be detached.
807 // The closest vertex should be merged.
808
809 std::vector< std::pair<unsigned, unsigned> > badPairs;
810
811 while( true ) {
812
813 // find the closest vertices pair that share the track of interest
814 double minSignificance { AlgConsts::maxValue };
815 unsigned nShared { 0 };
816
817 {
818 auto& vrtList = trackToVertexMap.at( maxSharedTrack );
819
820 auto nGood = std::count_if( vrtList.begin(), vrtList.end(), [&]( auto& v ) { return workVerticesContainer->at(v).isGood; } );
821 ATH_MSG_VERBOSE( " > " << __FUNCTION__ << ": size of good vertices = " << nGood );
822
823 std::vector< std::tuple< std::pair<unsigned, unsigned>, double, unsigned> > significanceTuple;
824 enum { kIndexPair, kSignificance, kNshared };
825
826 for( auto ivrt = vrtList.begin(); ivrt != vrtList.end(); ++ivrt ) {
827 for( auto jvrt = std::next( ivrt ); jvrt != vrtList.end(); ++jvrt ) {
828 auto pair = std::pair<unsigned, unsigned>( *ivrt, *jvrt );
829
830 if( !( workVerticesContainer->at(*ivrt).isGood ) ) continue;
831 if( !( workVerticesContainer->at(*jvrt).isGood ) ) continue;
832
833 // skip known bad pairs
834 if( std::find( badPairs.begin(), badPairs.end(), pair ) != badPairs.end() ) continue;
835
836 auto signif = significanceBetweenVertices( workVerticesContainer->at( *ivrt ), workVerticesContainer->at( *jvrt ) );
837
838 auto& ivrtTrks = workVerticesContainer->at(*ivrt).selectedTrackIndices;
839 auto& jvrtTrks = workVerticesContainer->at(*jvrt).selectedTrackIndices;
840
841 auto nSharedTracks = std::count_if( ivrtTrks.begin(), ivrtTrks.end(),
842 [&]( auto& index ) {
843 return std::find( jvrtTrks.begin(), jvrtTrks.end(), index ) != jvrtTrks.end();
844 } );
845
846 significanceTuple.emplace_back( pair, signif, nSharedTracks );
847 }
848 }
849
850 if( significanceTuple.empty() ) {
851 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": no vertex pairs are found --> exit the while loop." );
852 break;
853 }
854
855 auto minSignificanceTuple = std::min_element( significanceTuple.begin(), significanceTuple.end(), [&]( auto& t1, auto&t2 ) { return std::get<kSignificance>(t1) < std::get<kSignificance>(t2); } );
856
857 indexPair = std::get<kIndexPair> ( *minSignificanceTuple );
858 minSignificance = std::get<kSignificance> ( *minSignificanceTuple );
859 nShared = std::get<kNshared> ( *minSignificanceTuple );
860 }
861
862 ATH_MSG_VERBOSE( " > " << __FUNCTION__ << ": minSignificance = " << minSignificance );
863
864 if( minSignificance < m_VertexMergeCut || nShared >= 2 ) {
865
866 ATH_MSG_VERBOSE( " > " << __FUNCTION__ << ": attempt to merge vertices " << indexPair.first << " and " << indexPair.second );
867
868 WrkVrt vertex_backup1 = workVerticesContainer->at( indexPair.first );
869 WrkVrt vertex_backup2 = workVerticesContainer->at( indexPair.second );
870
871 StatusCode sc = mergeVertices( workVerticesContainer->at( indexPair.first ), workVerticesContainer->at( indexPair.second ) );
872
873 if( m_FillHist ) { m_hists["mergeType"]->Fill( RECONSTRUCT_NTRK ); }
874
875 if( sc.isFailure() ) {
876 // revert to the original
877 workVerticesContainer->at( indexPair.first ) = vertex_backup1;
878 workVerticesContainer->at( indexPair.second ) = vertex_backup2;
879 badPairs.emplace_back( indexPair );
880 }
881
882 // The second vertex is merged to the first.
883 // Explicity flag the second vertex is invalid.
884 workVerticesContainer->at( indexPair.second ).isGood = false;
885
886 // Now the vertex is merged and the bad pair record is outdated.
887 badPairs.clear();
888
889 mergeCounter++;
890
891 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Merged vertices " << indexPair.first << " and " << indexPair.second << ". merged vertex multiplicity = " << workVerticesContainer->at( indexPair.first ).selectedTrackIndices.size() );
892
893 } else {
894
895 // Here, the significance between closest vertices sharing the track is sufficiently distant
896 // and cannot be merged, while the track-association chi2 is small as well.
897 // In order to resolve the ambiguity anyway, remove the track from the worst-associated vertex.
898
899 auto& wrkvrt = workVerticesContainer->at( worstMatchingVertex );
900
901 auto end = std::remove_if( wrkvrt.selectedTrackIndices.begin(), wrkvrt.selectedTrackIndices.end(), [&]( auto& index ) { return index == maxSharedTrack; } );
902 wrkvrt.selectedTrackIndices.erase( end, wrkvrt.selectedTrackIndices.end() );
903
904 removeTrackCounter++;
905
906 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": removed track " << maxSharedTrack << " from vertex " << worstMatchingVertex );
907
908 if( wrkvrt.selectedTrackIndices.size() < 2 ) {
909 wrkvrt.isGood = false;
910 brokenCounter++;
911 break;
912 }
913
914 StatusCode sc = refitVertex( wrkvrt );
915 if( sc.isFailure() ) {
916 ATH_MSG_WARNING(" > " << __FUNCTION__ << ": detected vertex fitting failure!" );
917 }
918
919 break;
920
921 }
922 }
923
924 } else {
925
926 // Here, a bad track association is detected
927 // The track is detached from the worst-associated vertex and refit.
928
929 auto& wrkvrt = workVerticesContainer->at( worstMatchingVertex );
930
931 auto end = std::remove_if( wrkvrt.selectedTrackIndices.begin(), wrkvrt.selectedTrackIndices.end(), [&]( auto& index ) { return index == maxSharedTrack; } );
932 wrkvrt.selectedTrackIndices.erase( end, wrkvrt.selectedTrackIndices.end() );
933
934 if( wrkvrt.nTracksTotal() >=2 ) {
935
936 auto wrkvrt_backup = wrkvrt;
937 StatusCode sc = refitVertex( wrkvrt );
938 if( sc.isFailure() ) {
939 ATH_MSG_WARNING(" > " << __FUNCTION__ << ": detected vertex fitting failure!" );
940 wrkvrt = wrkvrt_backup;
941 }
942
943 } else {
944 wrkvrt.isGood = false;
945 brokenCounter++;
946 }
947
948 removeTrackCounter++;
949
950 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": removed track " << maxSharedTrack << " from vertex " << worstMatchingVertex );
951
952 }
953
954 }
955
956 //
957 // Try to improve vertices with big Chi2
958 for( auto& wrkvrt : *workVerticesContainer ) {
959
960 if(!wrkvrt.isGood ) continue; //don't work on wrkvrt which is already bad
961 if( wrkvrt.selectedTrackIndices.size() < 3 ) continue;
962
963 WrkVrt backup = wrkvrt;
964 improveVertexChi2( wrkvrt );
965 if( wrkvrt.fitQuality() > backup.fitQuality() ) wrkvrt = backup;
966
967 if( wrkvrt.nTracksTotal() < 2 ) wrkvrt.isGood = false;
968
969 }
970
971 if( m_FillNtuple ) {
972 m_ntupleVars->get<unsigned int>( "NumRearrSecVrt" )=workVerticesContainer->size();
973 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Size of Solution Set: "<< m_ntupleVars->get<unsigned int>( "NumRearrSecVrt" ));
974 }
975
976 ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
977 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Number of merges = " << mergeCounter << ", Number of track removal = " << removeTrackCounter << ", broken vertices = " << brokenCounter );
978 ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
979
980 return StatusCode::SUCCESS;
981 }
982
983
984 //____________________________________________________________________________________________________
985 StatusCode VrtSecInclusive::reassembleVertices( const EventContext&,
986 std::vector<WrkVrt>* workVerticesContainer )
987 {
988 // Here, the supposed issue is that, the position of the reconstructed vertex may be significantly
989 // displaced from its truth position, even if the constituent tracks are all from that truth.
990 // The fundamental reason of this is speculated that the VKalVrt vertex fitting could fall in
991 // a local minimum. This function attempts to improve the situation, given that N-track vertices
992 // are already reconstructed, by attempting to asociate a track of a small multiplicity vertex
993 // to another large multiplicity vertex.
994
995 unsigned reassembleCounter { 0 };
996
997 // First, sort WrkVrt by the track multiplicity
998 std::sort( workVerticesContainer->begin(), workVerticesContainer->end(), [](WrkVrt& v1, WrkVrt& v2) { return v1.selectedTrackIndices.size() < v2.selectedTrackIndices.size(); } );
999
1000 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": #vertices = " << workVerticesContainer->size() );
1001 // Loop over vertices (small -> large Ntrk order)
1002 for( auto& wrkvrt : *workVerticesContainer ) {
1003 if( !wrkvrt.isGood ) continue;
1004 if( wrkvrt.selectedTrackIndices.size() <= 1 ) continue;
1005
1006 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": vertex " << &wrkvrt << " #tracks = " << wrkvrt.selectedTrackIndices.size() );
1007 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": candidate vertex: "
1008 << " isGood = " << (wrkvrt.isGood? "true" : "false")
1009 << ", #ntrks = " << wrkvrt.nTracksTotal()
1010 << ", #selectedTracks = " << wrkvrt.selectedTrackIndices.size()
1011 << ", #associatedTracks = " << wrkvrt.associatedTrackIndices.size()
1012 << ", chi2/ndof = " << wrkvrt.fitQuality()
1013 << ", (r, z) = (" << wrkvrt.vertex.perp()
1014 <<", " << wrkvrt.vertex.z() << ")" );
1015
1016 std::map<unsigned, std::vector<WrkVrt>::reverse_iterator> mergiableVertex;
1017 std::set<std::vector<WrkVrt>::reverse_iterator> mergiableVerticesSet;
1018
1019 for( auto& index : wrkvrt.selectedTrackIndices ) {
1020
1021 const xAOD::TrackParticle* trk = m_selectedTracks.at( index );
1022
1023 mergiableVertex[index] = workVerticesContainer->rend();
1024
1025 std::vector<double> distances;
1026
1027 // Reverse iteration: large Ntrk -> small Ntrk order
1028 for( auto ritr = workVerticesContainer->rbegin(); ritr != workVerticesContainer->rend(); ++ritr ) {
1029 auto& targetVertex = *ritr;
1030
1031 if( &wrkvrt == &targetVertex ) continue;
1032 if( wrkvrt.selectedTrackIndices.size() >= targetVertex.selectedTrackIndices.size() ) continue;
1033
1034 // Get the closest approach
1035 std::vector<double> impactParameters;
1036 std::vector<double> impactParErrors;
1037
1038 if( !getSVImpactParameters(trk,targetVertex.vertex,impactParameters,impactParErrors) ) continue;
1039
1040 const auto& distance = hypot( impactParameters.at(0), impactParameters.at(1) );
1041 distances.emplace_back( distance );
1042
1043 if( std::abs( impactParameters.at(0) ) > m_reassembleMaxImpactParameterD0 ) continue;
1044 if( std::abs( impactParameters.at(1) ) > m_reassembleMaxImpactParameterZ0 ) continue;
1045
1046 mergiableVertex[index] = ritr;
1047 mergiableVerticesSet.emplace( ritr );
1048
1049 }
1050
1051 auto min_distance = !distances.empty() ? *(std::min_element( distances.begin(), distances.end() )) : AlgConsts::invalidFloat;
1052
1053 if( mergiableVertex[index] == workVerticesContainer->rend() ) {
1054 ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": track " << trk << " --> none : min distance = " << min_distance );
1055 } else {
1056 ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": track " << trk << " --> " << &( *(mergiableVertex[index]) ) << " --> size = " << mergiableVertex[index]->selectedTrackIndices.size() << ": min distance = " << min_distance );
1057 }
1058
1059 }
1060
1061 size_t count_mergiable = std::count_if( mergiableVertex.begin(), mergiableVertex.end(),
1062 [&](const std::pair<unsigned, std::vector<WrkVrt>::reverse_iterator>& p ) {
1063 return p.second != workVerticesContainer->rend(); } );
1064
1065 if( mergiableVerticesSet.size() == 1 && count_mergiable == wrkvrt.selectedTrackIndices.size() ) {
1066
1067 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": identified a unique association destination vertex" );
1068
1069 WrkVrt& destination = *( mergiableVertex.begin()->second );
1070 ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": destination #tracks before merging = " << destination.selectedTrackIndices.size() );
1071
1072 StatusCode sc = mergeVertices( destination, wrkvrt );
1073 if( sc.isFailure() ) {
1074 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": failure in vertex merging" );
1075 }
1076
1077 improveVertexChi2( destination );
1078
1079 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": merged destination vertex: "
1080 << " isGood = " << (destination.isGood? "true" : "false")
1081 << ", #ntrks = " << destination.nTracksTotal()
1082 << ", #selectedTracks = " << destination.selectedTrackIndices.size()
1083 << ", #associatedTracks = " << destination.associatedTrackIndices.size()
1084 << ", chi2/ndof = " << destination.fitQuality()
1085 << ", (r, z) = (" << destination.vertex.perp()
1086 <<", " << destination.vertex.z() << ")" );
1087
1088 if( m_FillHist ) { m_hists["mergeType"]->Fill( REASSEMBLE ); }
1089
1090 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": destination #tracks after merging = " << destination.selectedTrackIndices.size() );
1091
1092 reassembleCounter++;
1093
1094 }
1095
1096 }
1097
1098 ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
1099 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": reassembled vertices = " << reassembleCounter );
1100 ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
1101
1102 return StatusCode::SUCCESS;
1103 }
1104
1105
1106 //____________________________________________________________________________________________________
1107 StatusCode VrtSecInclusive::associateNonSelectedTracks( const EventContext& ctx,
1108 std::vector<WrkVrt>* workVerticesContainer )
1109 {
1110 SG::ReadHandle<xAOD::TrackParticleContainer> trackHandle = SG::makeHandle( m_TrackLocation, ctx );
1111 ATH_CHECK( trackHandle.isValid() );
1112 const xAOD::TrackParticleContainer *allTracks = trackHandle.cptr();
1113
1114 SG::ReadHandle<xAOD::VertexContainer> primVtxHandle = SG::makeHandle( m_PrimVrtLocation, ctx );
1115 ATH_CHECK( primVtxHandle.isValid() );
1116 const xAOD::VertexContainer *pvs = primVtxHandle.cptr();
1117
1118 if( !m_decor_isAssociated ) {
1119 m_decor_isAssociated.emplace ( "is_associated" + m_augVerString );
1120 }
1121
1122 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": #verticess = " << workVerticesContainer->size() );
1123
1124 unsigned associateCounter { 0 };
1125
1126 // Loop over vertices
1127 for( auto& wrkvrt : *workVerticesContainer ) {
1128
1129 if( !wrkvrt.isGood ) continue;
1130 if( wrkvrt.selectedTrackIndices.size() <= 1 ) continue;
1131
1132 improveVertexChi2( wrkvrt );
1133
1134 wrkvrt.Chi2_core = wrkvrt.Chi2;
1135
1136 auto& vertexPos = wrkvrt.vertex;
1137
1138 std::vector<double> distanceToPVs;
1139
1140 for( const auto* pv : *pvs ) {
1141 distanceToPVs.emplace_back( VKalVrtAthena::vtxVtxDistance( vertexPos, pv->position() ) );
1142 }
1143 const auto& minDistance = *( std::min_element( distanceToPVs.begin(), distanceToPVs.end() ) );
1144
1145 if( minDistance < m_associateMinDistanceToPV ) continue;
1146
1147
1148 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": vertex pos = (" << vertexPos.x() << ", " << vertexPos.y() << ", " << vertexPos.z() << "), "
1149 "#selected = " << wrkvrt.selectedTrackIndices.size() << ", #assoc = " << wrkvrt.associatedTrackIndices.size() );
1150
1151 std::vector<const xAOD::TrackParticle*> candidates;
1152
1153 // Search for candidate tracks
1154 for( auto itr = allTracks->begin(); itr != allTracks->end(); ++itr ) {
1155 const auto* trk = *itr;
1156
1157 // If the track is already used for any DV candidate, reject.
1158 {
1159 auto result = std::find_if( workVerticesContainer->begin(), workVerticesContainer->end(),
1160 [&] ( WrkVrt& wrkvrt ) {
1161 auto found = std::find_if( wrkvrt.selectedTrackIndices.begin(), wrkvrt.selectedTrackIndices.end(),
1162 [&]( long int index ) {
1163 // when using selected tracks from electrons, also check the orginal track particle from GSF to see if InDetTrackParticle (trk) is an electron that is already in the vertex
1164 if (m_doSelectTracksFromElectrons || m_doSelectIDAndGSFTracks) {
1165 const xAOD::TrackParticle *id_tr;
1166 id_tr = xAOD::EgammaHelpers::getOriginalTrackParticleFromGSF(m_selectedTracks.at(index));
1167 return trk == m_selectedTracks.at(index) or trk == id_tr;
1168 }
1169 else{
1170 return trk == m_selectedTracks.at(index);
1171 }
1172 } );
1173 return found != wrkvrt.selectedTrackIndices.end();
1174 } );
1175 if( result != workVerticesContainer->end() ) continue;
1176 }
1177
1178 // If the track is already registered to the associated track list, reject.
1179 {
1180 auto result = std::find_if( m_associatedTracks.begin(), m_associatedTracks.end(),
1181 [&] (const auto* atrk) { return trk == atrk; } );
1182 if( result != m_associatedTracks.end() ) continue;
1183 }
1184
1185 // Reject PV-associated tracks
1186 // if( !selectTrack_notPVassociated( trk ) ) continue;
1187
1188 // pT selection
1189 if( trk->pt() < m_associatePtCut ) continue;
1190
1191 // chi2 selection
1192 if( trk->chiSquared() / trk->numberDoF() > m_associateChi2Cut ) continue;
1193
1194 // Hit pattern consistentcy requirement
1195 if( !checkTrackHitPatternToVertexOuterOnly( trk, vertexPos ) ) continue;
1196
1197 // Get the closest approach
1198 std::vector<double> impactParameters;
1199 std::vector<double> impactParErrors;
1200
1201 if( !getSVImpactParameters( trk, vertexPos, impactParameters, impactParErrors) ) continue;
1202
1203 if( std::abs( impactParameters.at(0) ) / sqrt( impactParErrors.at(0) ) > m_associateMaxD0Signif ) continue;
1204 if( std::abs( impactParameters.at(1) ) / sqrt( impactParErrors.at(1) ) > m_associateMaxZ0Signif ) continue;
1205
1206 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": trk " << trk
1207 << ": d0 to vtx = " << impactParameters.at(k_d0)
1208 << ", z0 to vtx = " << impactParameters.at(k_z0)
1209 << ", distance to vtx = " << hypot( impactParameters.at(k_d0), impactParameters.at(k_z0) ) );
1210
1211 candidates.emplace_back( trk );
1212
1213 }
1214
1215 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": number of candidate tracks = " << candidates.size() );
1216
1217 std::unique_ptr<Trk::IVKalState> state = m_fitSvc->makeState();
1218 // Attempt to add the track to the vertex and try fitting
1219 for( const auto* trk : candidates ) {
1220
1221 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": attempting to associate track = " << trk );
1222
1223 // Backup the current vertes status
1224 WrkVrt wrkvrt_backup = wrkvrt;
1225
1226 m_fitSvc->setApproximateVertex( vertexPos.x(), vertexPos.y(), vertexPos.z(), *state );
1227
1228 std::vector<const xAOD::TrackParticle*> baseTracks;
1229 std::vector<const xAOD::NeutralParticle*> dummyNeutrals;
1230
1231 wrkvrt.Chi2PerTrk.clear();
1232
1233 for( const auto& index : wrkvrt.selectedTrackIndices ) {
1234 baseTracks.emplace_back( m_selectedTracks.at( index ) );
1235 wrkvrt.Chi2PerTrk.emplace_back( AlgConsts::chi2PerTrackInitValue );
1236 }
1237 for( const auto& index : wrkvrt.associatedTrackIndices ) {
1238 baseTracks.emplace_back( m_associatedTracks.at( index ) );
1239 wrkvrt.Chi2PerTrk.emplace_back( AlgConsts::chi2PerTrackInitValue );
1240 }
1241
1242 baseTracks.emplace_back( trk );
1243 wrkvrt.Chi2PerTrk.emplace_back( AlgConsts::chi2PerTrackInitValue );
1244
1245 Amg::Vector3D initPos;
1246
1247 {
1248 StatusCode sc = m_fitSvc->VKalVrtFitFast( baseTracks, initPos, *state );/* Fast crude estimation */
1249
1250 if( sc.isFailure() ) ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": fast crude estimation failed.");
1251
1252 const auto& diffPos = initPos - vertexPos;
1253
1254 if( diffPos.norm() > 10. ) {
1255
1256 ATH_MSG_VERBOSE( " > " << __FUNCTION__ << ": approx vertex as original" );
1257 m_fitSvc->setApproximateVertex( vertexPos.x(), vertexPos.y(), vertexPos.z(), *state );
1258
1259 } else {
1260
1261 ATH_MSG_VERBOSE( " > " << __FUNCTION__ << ": approx vertex set to (" << initPos.x() << ", " << initPos.y() << ", " << initPos.z() << ")" );
1262 m_fitSvc->setApproximateVertex( initPos.x(), initPos.y(), initPos.z(), *state );
1263
1264 }
1265 }
1266
1267
1268 ATH_MSG_VERBOSE( " > " << __FUNCTION__ << ": now vertex fitting..." );
1269
1270 StatusCode sc = m_fitSvc->VKalVrtFit(baseTracks, dummyNeutrals,
1271 wrkvrt.vertex,
1272 wrkvrt.vertexMom,
1273 wrkvrt.Charge,
1274 wrkvrt.vertexCov,
1275 wrkvrt.Chi2PerTrk,
1276 wrkvrt.TrkAtVrt,
1277 wrkvrt.Chi2,
1278 *state);
1279
1280 if( sc.isFailure() ) {
1281 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": VKalVrtFit failure. Revert to backup");
1282 wrkvrt = wrkvrt_backup;
1283
1284 if( m_FillHist ) m_hists["associateMonitor"]->Fill( 1 );
1285
1286 continue;
1287 }
1288
1289
1290 if( m_FillHist ) m_hists["associateMonitor"]->Fill( 0 );
1291
1292 auto& cov = wrkvrt.vertexCov;
1293
1294 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": succeeded in associating. New vertex pos = (" << vertexPos.perp() << ", " << vertexPos.z() << ", " << vertexPos.perp()*vertexPos.phi() << ")" );
1295 ATH_MSG_VERBOSE( " > " << __FUNCTION__ << ": New vertex cov = (" << cov.at(0) << ", " << cov.at(1) << ", " << cov.at(2) << ", " << cov.at(3) << ", " << cov.at(4) << ", " << cov.at(5) << ")" );
1296
1297 associateCounter++;
1298
1299 wrkvrt.associatedTrackIndices.emplace_back( m_associatedTracks.size() );
1300
1301 m_associatedTracks.emplace_back( trk );
1302 (*m_decor_isAssociated)( *trk ) = true;
1303
1304 }
1305
1306 }
1307
1308 ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
1309 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": total associated number of tracks = " << associateCounter );
1310 ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
1311
1312 return StatusCode::SUCCESS;
1313 }
1314
1315
1316 //____________________________________________________________________________________________________
1317 StatusCode VrtSecInclusive::mergeByShuffling( const EventContext&,
1318 std::vector<WrkVrt> *workVerticesContainer )
1319 {
1320
1321 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": #verticess = " << workVerticesContainer->size() );
1322
1323 unsigned mergeCounter { 0 };
1324
1325 // First, sort WrkVrt by the track multiplicity
1326 std::sort( workVerticesContainer->begin(), workVerticesContainer->end(), [](WrkVrt& v1, WrkVrt& v2) { return v1.selectedTrackIndices.size() < v2.selectedTrackIndices.size(); } );
1327
1328 // Loop over vertices (small -> large Ntrk order)
1329 for( auto& wrkvrt : *workVerticesContainer ) {
1330 if( !wrkvrt.isGood ) continue;
1331 if( wrkvrt.selectedTrackIndices.size() <= 1 ) continue;
1332
1333 // Reverse iteration: large Ntrk -> small Ntrk order
1334 for( auto ritr = workVerticesContainer->rbegin(); ritr != workVerticesContainer->rend(); ++ritr ) {
1335 auto& vertexToMerge = *ritr;
1336
1337 if( !vertexToMerge.isGood ) continue;
1338 if( vertexToMerge.selectedTrackIndices.size() <= 1 ) continue;
1339 if( &wrkvrt == &vertexToMerge ) continue;
1340 if( vertexToMerge.selectedTrackIndices.size() < wrkvrt.selectedTrackIndices.size() ) continue;
1341
1342 const double& significance = significanceBetweenVertices( wrkvrt, vertexToMerge );
1343
1344 if( significance > m_mergeByShufflingMaxSignificance ) continue;
1345
1346 bool mergeFlag { false };
1347
1348 ATH_MSG_DEBUG(" > " << __FUNCTION__
1349 << ": vertex " << &wrkvrt << " #tracks = " << wrkvrt.selectedTrackIndices.size()
1350 << " --> to Merge : " << &vertexToMerge << ", #tracks = " << vertexToMerge.selectedTrackIndices.size()
1351 << " significance = " << significance );
1352
1353 double min_signif = AlgConsts::maxValue;
1354
1355 // Method 1. Assume that the solution is somewhat wrong, and the solution gets correct if it starts from the other vertex position
1356 if( m_doSuggestedRefitOnMerging && !mergeFlag ) {
1357 WrkVrt testVertex = wrkvrt;
1358 StatusCode sc = refitVertexWithSuggestion( testVertex, vertexToMerge.vertex );
1359 if( sc.isFailure() ) {
1360 //ATH_MSG_WARNING(" > " << __FUNCTION__ << ": detected vertex fitting failure!" );
1361 } else {
1362
1363 const auto signif = significanceBetweenVertices( testVertex, vertexToMerge );
1364 if( signif < min_signif ) min_signif = signif;
1365
1366 if( signif < m_mergeByShufflingAllowance ) {
1367 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": method1: vertexToMerge " << &vertexToMerge << ": test signif = " << signif );
1368 mergeFlag = true;
1369
1370 }
1371
1372 if( m_FillHist && min_signif > 0. ) m_hists["shuffleMinSignif1"]->Fill( log10( min_signif ) );
1373 if( m_FillHist && mergeFlag ) { m_hists["mergeType"]->Fill( SHUFFLE1 ); }
1374 }
1375 }
1376
1377 // Method 2. magnet merging: borrowing another track from the target vertex to merge
1378 if( m_doMagnetMerging && !mergeFlag ) {
1379
1380 // Loop over tracks in vertexToMerge
1381 for( auto& index : vertexToMerge.selectedTrackIndices ) {
1382
1383 WrkVrt testVertex = wrkvrt;
1384 testVertex.selectedTrackIndices.emplace_back( index );
1385
1386 StatusCode sc = refitVertexWithSuggestion( testVertex, vertexToMerge.vertex );
1387 if( sc.isFailure() ) {
1388 //ATH_MSG_WARNING(" > " << __FUNCTION__ << ": detected vertex fitting failure!" );
1389 } else {
1390
1391 const auto signif = significanceBetweenVertices( testVertex, vertexToMerge );
1392 if( signif < min_signif ) min_signif = signif;
1393
1394 if( signif < m_mergeByShufflingAllowance ) {
1395 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": method2: vertexToMerge " << &vertexToMerge << " track index " << index << ": test signif = " << signif );
1396 mergeFlag = true;
1397 }
1398
1399 }
1400 }
1401
1402 if( m_FillHist && min_signif > 0. ) m_hists["shuffleMinSignif2"]->Fill( log10( min_signif ) );
1403
1404 if( m_FillHist && mergeFlag ) { m_hists["mergeType"]->Fill( SHUFFLE2 ); }
1405 }
1406
1407 // Method 3. Attempt to force merge
1408 if( m_doWildMerging && !mergeFlag ) {
1409
1410 WrkVrt testVertex = wrkvrt;
1411
1412 for( auto& index : vertexToMerge.selectedTrackIndices ) {
1413 testVertex.selectedTrackIndices.emplace_back( index );
1414 }
1415
1416 StatusCode sc = refitVertexWithSuggestion( testVertex, vertexToMerge.vertex );
1417 if( sc.isFailure() ) {
1418 //ATH_MSG_WARNING(" > " << __FUNCTION__ << ": detected vertex fitting failure!" );
1419 } else {
1420
1421 const auto signif = significanceBetweenVertices( testVertex, vertexToMerge );
1422 if( signif < min_signif ) min_signif = signif;
1423
1424 if( signif < m_mergeByShufflingAllowance ) {
1425 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": method3: vertexToMerge " << &vertexToMerge << ": test signif = " << signif );
1426 mergeFlag = true;
1427 }
1428
1429 if( m_FillHist && min_signif > 0. ) m_hists["shuffleMinSignif3"]->Fill( log10( min_signif ) );
1430 if( m_FillHist && mergeFlag ) { m_hists["mergeType"]->Fill( SHUFFLE3 ); }
1431
1432 }
1433 }
1434
1435
1436 if( mergeFlag ) {
1437 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": vertexToMerge " << &vertexToMerge << " ==> min signif = " << min_signif << " judged to merge" );
1438
1439 auto vertexToMerge_backup = vertexToMerge;
1440 auto wrkvrt_backup = wrkvrt;
1441
1442 StatusCode sc = mergeVertices( vertexToMerge, wrkvrt );
1443 if( sc.isFailure() ) {
1444 vertexToMerge = vertexToMerge_backup;
1445 wrkvrt = wrkvrt_backup;
1446 continue;
1447 }
1448
1449 improveVertexChi2( wrkvrt );
1450
1451 mergeCounter++;
1452 }
1453
1454 }
1455
1456 }
1457
1458 ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
1459 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Number of merges = " << mergeCounter );
1460 ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
1461
1462 return StatusCode::SUCCESS;
1463 }
1464
1465
1466 //____________________________________________________________________________________________________
1467 StatusCode VrtSecInclusive::mergeFinalVertices( const EventContext&,
1468 std::vector<WrkVrt> *workVerticesContainer )
1469 {
1470
1471 unsigned mergeCounter { 0 };
1472
1473 while (true) {
1474 //
1475 // Minimal vertex-vertex distance
1476 //
1477 for( auto& wrkvrt : *workVerticesContainer) {
1478 wrkvrt.closestWrkVrtIndex = AlgConsts::invalidUnsigned;
1479 wrkvrt.closestWrkVrtValue = AlgConsts::maxValue;
1480 }
1481
1482 std::pair<unsigned, unsigned> indexPair { AlgConsts::invalidUnsigned, AlgConsts::invalidUnsigned };
1483 auto minDistance = findMinVerticesPair( workVerticesContainer, indexPair, &VrtSecInclusive::distanceBetweenVertices );
1484
1485 if( minDistance == AlgConsts::maxValue ) break;
1486 if( indexPair.first == AlgConsts::invalidUnsigned ) break;
1487 if( indexPair.second == AlgConsts::invalidUnsigned ) break;
1488
1489 auto& v1 = workVerticesContainer->at(indexPair.first);
1490 auto& v2 = workVerticesContainer->at(indexPair.second);
1491
1492 const double averageRadius = ( v1.vertex.perp() + v2.vertex.perp() ) / 2.0;
1493
1494 if( minDistance > m_VertexMergeFinalDistCut + m_VertexMergeFinalDistScaling * averageRadius ) {
1495 ATH_MSG_DEBUG( "Vertices " << indexPair.first << " and " << indexPair.second
1496 <<" are separated by distance " << minDistance );
1497 break;
1498 }
1499
1500 ATH_MSG_DEBUG( "Merging FINAL vertices " << indexPair.first << " and " << indexPair.second
1501 <<" which are separated by distance "<< minDistance );
1502
1503 StatusCode sc = mergeVertices( v1, v2 );
1504 if( sc.isFailure() ) {}
1505 if( m_FillHist ) { m_hists["mergeType"]->Fill( FINAL ); }
1506
1507 improveVertexChi2( v1 );
1508
1509 mergeCounter++;
1510
1511 }
1512
1513 ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
1514 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Number of merges = " << mergeCounter );
1515 ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
1516
1517 return StatusCode::SUCCESS;
1518
1519 } // end of mergeFinalVertices
1520
1521
1522
1523 //____________________________________________________________________________________________________
1524 StatusCode VrtSecInclusive::refitAndSelectGoodQualityVertices( const EventContext& ctx,
1525 std::vector<WrkVrt> *workVerticesContainer )
1526 {
1527 // Output SVs as xAOD::Vertex
1528 // Needs a conversion function from workVerticesContainer to xAOD::Vertex here.
1529 // The supposed form of the function will be as follows:
1530
1531 try {
1532 SG::WriteHandle<xAOD::VertexContainer> secVtxHandle = SG::makeHandle( m_vertexKey, ctx );
1533 ATH_CHECK( secVtxHandle.record( std::make_unique<xAOD::VertexContainer>(),
1534 std::make_unique<xAOD::VertexAuxContainer>() ) );
1535 xAOD::VertexContainer *secondaryVertexContainer = secVtxHandle.ptr();
1536 m_vertexCollectionsDefinitions[m_vertexKey.key()] = true;
1537
1538 enum { kPt, kEta, kPhi, kD0, kZ0, kErrP, kErrD0, kErrZ0, kChi2SV };
1539 if( m_trkDecors.empty() ) {
1540 m_trkDecors.emplace( kPt, SG::AuxElement::Decorator<float>("pt_wrtSV" + m_augVerString) );
1541 m_trkDecors.emplace( kEta, SG::AuxElement::Decorator<float>("eta_wrtSV" + m_augVerString) );
1542 m_trkDecors.emplace( kPhi, SG::AuxElement::Decorator<float>("phi_wrtSV" + m_augVerString) );
1543 m_trkDecors.emplace( kD0, SG::AuxElement::Decorator<float>("d0_wrtSV" + m_augVerString) );
1544 m_trkDecors.emplace( kZ0, SG::AuxElement::Decorator<float>("z0_wrtSV" + m_augVerString) );
1545 m_trkDecors.emplace( kErrP, SG::AuxElement::Decorator<float>("errP_wrtSV" + m_augVerString) );
1546 m_trkDecors.emplace( kErrD0, SG::AuxElement::Decorator<float>("errd0_wrtSV" + m_augVerString) );
1547 m_trkDecors.emplace( kErrZ0, SG::AuxElement::Decorator<float>("errz0_wrtSV" + m_augVerString) );
1548 m_trkDecors.emplace( kChi2SV, SG::AuxElement::Decorator<float>("chi2_toSV" + m_augVerString) );
1549 }
1550 if( !m_decor_is_svtrk_final ) {
1551 m_decor_is_svtrk_final.emplace ( "is_svtrk_final" + m_augVerString );
1552 }
1553
1554 std::map<const WrkVrt*, const xAOD::Vertex*> wrkvrtLinkMap;
1555
1556 //----------------------------------------------------------
1557
1558 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": input #vertices = " << workVerticesContainer->size() );
1559
1560 // Loop over vertices
1561 for( auto& wrkvrt : *workVerticesContainer ) {
1562
1563 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": candidate vertex: "
1564 << " isGood = " << (wrkvrt.isGood? "true" : "false")
1565 << ", #ntrks = " << wrkvrt.nTracksTotal()
1566 << ", #selectedTracks = " << wrkvrt.selectedTrackIndices.size()
1567 << ", #associatedTracks = " << wrkvrt.associatedTrackIndices.size()
1568 << ", chi2/ndof = " << wrkvrt.Chi2 / ( wrkvrt.ndof() + AlgConsts::infinitesimal )
1569 << ", (r, z) = (" << wrkvrt.vertex.perp()
1570 <<", " << wrkvrt.vertex.z() << ")" );
1571
1572 if( m_FillHist ) m_hists["finalCutMonitor"]->Fill( 0 );
1573
1574 if( m_removeFakeVrt && m_removeFakeVrtLate ) {
1575 removeInconsistentTracks( wrkvrt );
1576 }
1577
1578 if( wrkvrt.nTracksTotal() < 2 ) {
1579 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": ntrk < 2 --> rejected." );
1580 continue; /* Bad vertices */
1581 }
1582
1583 if( m_FillHist ) m_hists["finalCutMonitor"]->Fill( 1 );
1584
1585
1586 // Remove track if the vertex is inner than IBL and the track does not have pixel hits!
1587 if( wrkvrt.vertex.perp() < 31.0 ) {
1588
1589 // for selected tracks
1590 wrkvrt.selectedTrackIndices.erase( std::remove_if( wrkvrt.selectedTrackIndices.begin(), wrkvrt.selectedTrackIndices.end(),
1591 [&]( auto& index ) {
1592 auto* trk = m_selectedTracks.at( index );
1593 uint8_t nPixelHits { 0 }; trk->summaryValue( nPixelHits, xAOD::numberOfPixelHits );
1594 return ( nPixelHits < 3 );
1595 } ),
1596 wrkvrt.selectedTrackIndices.end() );
1597
1598 // for associated tracks
1599 wrkvrt.associatedTrackIndices.erase( std::remove_if( wrkvrt.associatedTrackIndices.begin(), wrkvrt.associatedTrackIndices.end(),
1600 [&]( auto& index ) {
1601 auto* trk = m_associatedTracks.at( index );
1602 uint8_t nPixelHits { 0 }; trk->summaryValue( nPixelHits, xAOD::numberOfPixelHits );
1603 return ( nPixelHits < 3 );
1604 } ),
1605 wrkvrt.associatedTrackIndices.end() );
1606
1607 auto statusCode = refitVertex( wrkvrt );
1608 if( statusCode.isFailure() ) {}
1609
1610 }
1611
1612
1613 if( m_doFinalImproveChi2 ) {
1614
1615 WrkVrt backup = wrkvrt;
1616
1617 improveVertexChi2( wrkvrt );
1618
1619 if( wrkvrt.fitQuality() > backup.fitQuality() ) wrkvrt = backup;
1620
1621 }
1622
1623 // If the number of remaining tracks is less than 2, drop.
1624 if( wrkvrt.nTracksTotal() < 2 ) continue;
1625
1626 // Select only vertices with keeping more than 2 selectedTracks
1627 if( wrkvrt.selectedTrackIndices.size() < 2 ) continue;
1628
1629
1630 if( m_FillHist ) m_hists["finalCutMonitor"]->Fill( 2 );
1631
1632
1633 {
1634 WrkVrt backup = wrkvrt;
1635
1636 StatusCode sc = refitVertex( wrkvrt );
1637 if( sc.isFailure() ) {
1638
1639 auto indices = wrkvrt.associatedTrackIndices;
1640
1641 wrkvrt.associatedTrackIndices.clear();
1642 sc = refitVertex( wrkvrt );
1643 if( sc.isFailure() ) {
1644 ATH_MSG_WARNING(" > " << __FUNCTION__ << ": detected vertex fitting failure!" );
1645 wrkvrt = backup;
1646 }
1647 if( wrkvrt.fitQuality() > backup.fitQuality() ) wrkvrt = backup;
1648
1649 for( auto& index : indices ) {
1650 backup = wrkvrt;
1651 wrkvrt.associatedTrackIndices.emplace_back( index );
1652 sc = refitVertex( wrkvrt );
1653 if( sc.isFailure() || TMath::Prob( wrkvrt.Chi2, wrkvrt.ndof() ) < m_improveChi2ProbThreshold ) {
1654 ATH_MSG_WARNING(" > " << __FUNCTION__ << ": detected vertex fitting failure!" );
1655 wrkvrt = backup;
1656 continue;
1657 }
1658 }
1659
1660 } else {
1661 if( wrkvrt.fitQuality() > backup.fitQuality() ) wrkvrt = backup;
1662 }
1663 }
1664
1665 if( m_FillHist ) m_hists["finalCutMonitor"]->Fill( 3 );
1666
1667 //
1668 // Store good vertices into StoreGate
1669 //
1670 if( m_FillNtuple ) m_ntupleVars->get<unsigned int>( "NumSecVrt" )++;
1671
1672 TLorentzVector sumP4_pion;
1673 TLorentzVector sumP4_electron;
1674 TLorentzVector sumP4_proton;
1675
1676 // Pre-check before storing vertex if the SV perigee is available
1677 bool good_flag = true;
1678
1679 std::map<const std::deque<long int>*, const std::vector<const xAOD::TrackParticle*>&> indicesSet
1680 = {
1681 { &(wrkvrt.selectedTrackIndices), m_selectedTracks },
1682 { &(wrkvrt.associatedTrackIndices), m_associatedTracks }
1683 };
1684
1685 for( auto& pair : indicesSet ) {
1686
1687 const auto* indices = pair.first;
1688 const auto& tracks = pair.second;
1689
1690 for( const auto& itrk : *indices ) {
1691 const auto* trk = tracks.at( itrk );
1692 auto sv_perigee = m_trackToVertexTool->perigeeAtVertex(ctx, *trk, wrkvrt.vertex );
1693 if( !sv_perigee ) {
1694 ATH_MSG_INFO(" > " << __FUNCTION__ << ": > Track index " << trk->index() << ": Failed in obtaining the SV perigee!" );
1695 good_flag = false;
1696 }
1697 }
1698
1699 }
1700
1701 if( !good_flag ) {
1702 ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": sv perigee could not be obtained --> rejected" );
1703 continue;
1704 }
1705
1706 if( m_FillHist ) m_hists["finalCutMonitor"]->Fill( 4 );
1707
1708
1709 std::vector<const xAOD::TrackParticle*> tracks;
1710 std::vector< std::pair<const xAOD::TrackParticle*, double> > trackChi2Pairs;
1711
1712 {
1713
1714 for( auto& pair : indicesSet ) {
1715 for( const auto& index : *pair.first ) tracks.emplace_back( pair.second.at( index ) );
1716 }
1717
1718 auto trkitr = tracks.begin();
1719 auto chi2itr = wrkvrt.Chi2PerTrk.begin();
1720
1721 for( ; ( trkitr!=tracks.end() && chi2itr!=wrkvrt.Chi2PerTrk.end() ); ++trkitr, ++chi2itr ) {
1722 trackChi2Pairs.emplace_back( *trkitr, *chi2itr );
1723 }
1724
1725 }
1726
1727
1728 TLorentzVector sumP4_selected;
1729
1730 bool badIPflag { false };
1731
1732 // loop over vertex tracks
1733 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Track loop: size = " << tracks.size() );
1734 for( auto& pair : trackChi2Pairs ) {
1735
1736 const auto* trk = pair.first;
1737 const auto& chi2AtSV = pair.second;
1738
1739 ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": > Track index " << trk->index() << ": start." );
1740
1741 track_summary trk_summary;
1742 fillTrackSummary( trk_summary, trk );
1743
1744 //
1745 // calculate mass/pT of tracks and track parameters
1746 //
1747
1748 double trk_pt = trk->pt();
1749 double trk_eta = trk->eta();
1750 double trk_phi = trk->phi();
1751
1752 ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": > Track index " << trk->index() << ": in vrt chg/pt/phi/eta = "
1753 << trk->charge() <<","
1754 <<trk_pt<<","
1755 <<trk_phi<<","
1756 <<trk_eta);
1757
1759 // Get the perigee of the track at the vertex
1760 ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": > Track index " << trk->index() << ": Get the prigee of the track at the vertex." );
1761
1762 auto sv_perigee = m_trackToVertexTool->perigeeAtVertex(ctx, *trk, wrkvrt.vertex );
1763 if( !sv_perigee ) {
1764 ATH_MSG_WARNING(" > " << __FUNCTION__ << ": > Track index " << trk->index() << ": Failed in obtaining the SV perigee!" );
1765
1766 for( auto& pair : m_trkDecors ) {
1767 pair.second( *trk ) = AlgConsts::invalidFloat;
1768 }
1769 (*m_decor_is_svtrk_final)( *trk ) = true;
1770 continue;
1771 }
1772
1773 double qOverP_wrtSV = sv_perigee->parameters() [Trk::qOverP];
1774 double theta_wrtSV = sv_perigee->parameters() [Trk::theta];
1775 double p_wrtSV = 1.0 / std::abs( qOverP_wrtSV );
1776 double pt_wrtSV = p_wrtSV * sin( theta_wrtSV );
1777 double eta_wrtSV = -log( tan( theta_wrtSV/2. ) );
1778 double phi_wrtSV = sv_perigee->parameters() [Trk::phi];
1779 double d0_wrtSV = sv_perigee->parameters() [Trk::d0];
1780 double z0_wrtSV = sv_perigee->parameters() [Trk::z0];
1781 double errd0_wrtSV = (*sv_perigee->covariance())( Trk::d0, Trk::d0 );
1782 double errz0_wrtSV = (*sv_perigee->covariance())( Trk::z0, Trk::z0 );
1783 double errP_wrtSV = (*sv_perigee->covariance())( Trk::qOverP, Trk::qOverP );
1784
1785 // xAOD::Track augmentation
1786 ( m_trkDecors.at(kPt) )( *trk ) = pt_wrtSV;
1787 ( m_trkDecors.at(kEta) )( *trk ) = eta_wrtSV;
1788 ( m_trkDecors.at(kPhi) )( *trk ) = phi_wrtSV;
1789 ( m_trkDecors.at(kD0) )( *trk ) = d0_wrtSV;
1790 ( m_trkDecors.at(kZ0) )( *trk ) = z0_wrtSV;
1791 ( m_trkDecors.at(kErrP) )( *trk ) = errP_wrtSV;
1792 ( m_trkDecors.at(kErrD0) )( *trk ) = errd0_wrtSV;
1793 ( m_trkDecors.at(kErrZ0) )( *trk ) = errz0_wrtSV;
1794 ( m_trkDecors.at(kChi2SV))( *trk ) = chi2AtSV;
1795
1796 (*m_decor_is_svtrk_final)( *trk ) = true;
1797
1798 TLorentzVector p4wrtSV_pion;
1799 TLorentzVector p4wrtSV_electron;
1800 TLorentzVector p4wrtSV_proton;
1801
1802 p4wrtSV_pion .SetPtEtaPhiM( pt_wrtSV, eta_wrtSV, phi_wrtSV, PhysConsts::mass_chargedPion );
1803 p4wrtSV_electron.SetPtEtaPhiM( pt_wrtSV, eta_wrtSV, phi_wrtSV, PhysConsts::mass_electron );
1804
1805 // for selected tracks only
1806 static const SG::ConstAccessor<char> is_associatedAcc("is_associated" + m_augVerString);
1807 if( is_associatedAcc.isAvailable(*trk) ) {
1808 if( !is_associatedAcc(*trk) ) {
1809 sumP4_selected += p4wrtSV_pion;
1810 }
1811 } else {
1812 sumP4_selected += p4wrtSV_pion;
1813 }
1814
1815 sumP4_pion += p4wrtSV_pion;
1816 sumP4_electron += p4wrtSV_electron;
1817 sumP4_proton += p4wrtSV_proton;
1818
1819 ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": > Track index " << trk->index() << ": end." );
1820 } // loop over tracks in vertex
1821
1822 ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": Track loop end. ");
1823
1824 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Final Sec.Vertex=" << wrkvrt.nTracksTotal() <<", "
1825 <<wrkvrt.vertex.perp() <<", "<<wrkvrt.vertex.z() <<", "
1826 <<wrkvrt.vertex.phi() <<", mass = "<< sumP4_pion.M() << "," << sumP4_electron.M() );
1827
1828 // Save the perigee parameters for the first two tracks
1829 float perigee_x_trk1 = 0.0;
1830 float perigee_y_trk1 = 0.0;
1831 float perigee_z_trk1 = 0.0;
1832 float perigee_x_trk2 = 0.0;
1833 float perigee_y_trk2 = 0.0;
1834 float perigee_z_trk2 = 0.0;
1835 float perigee_px_trk1 = 0.0;
1836 float perigee_py_trk1 = 0.0;
1837 float perigee_pz_trk1 = 0.0;
1838 float perigee_px_trk2 = 0.0;
1839 float perigee_py_trk2 = 0.0;
1840 float perigee_pz_trk2 = 0.0;
1841 float perigee_cov_xx_trk1 = 0.0;
1842 float perigee_cov_xy_trk1 = 0.0;
1843 float perigee_cov_xz_trk1 = 0.0;
1844 float perigee_cov_yy_trk1 = 0.0;
1845 float perigee_cov_yz_trk1 = 0.0;
1846 float perigee_cov_zz_trk1 = 0.0;
1847 float perigee_cov_xx_trk2 = 0.0;
1848 float perigee_cov_xy_trk2 = 0.0;
1849 float perigee_cov_xz_trk2 = 0.0;
1850 float perigee_cov_yy_trk2 = 0.0;
1851 float perigee_cov_yz_trk2 = 0.0;
1852 float perigee_cov_zz_trk2 = 0.0;
1853 float perigee_d0_trk1 = 0.0;
1854 float perigee_d0_trk2 = 0.0;
1855 float perigee_z0_trk1 = 0.0;
1856 float perigee_z0_trk2 = 0.0;
1857 float perigee_qOverP_trk1 = 0.0;
1858 float perigee_qOverP_trk2 = 0.0;
1859 float perigee_theta_trk1 = 0.0;
1860 float perigee_theta_trk2 = 0.0;
1861 float perigee_phi_trk1 = 0.0;
1862 float perigee_phi_trk2 = 0.0;
1863 int perigee_charge_trk1 = 0;
1864 int perigee_charge_trk2 = 0;
1865 float perigee_distance = 9999.0;
1866
1867 Amg::Vector3D vDist = wrkvrt.vertex - m_thePV->position();
1868 float vPos = (vDist.x() * wrkvrt.vertexMom.Px() + vDist.y() * wrkvrt.vertexMom.Py() + vDist.z() * wrkvrt.vertexMom.Pz()) / wrkvrt.vertexMom.Rho();
1869 float vPosMomAngT = (vDist.x() * wrkvrt.vertexMom.Px() + vDist.y() * wrkvrt.vertexMom.Py()) / vDist.perp() / wrkvrt.vertexMom.Pt();
1870 float vPosMomAng3D = (vDist.x() * wrkvrt.vertexMom.Px() + vDist.y() * wrkvrt.vertexMom.Py() + vDist.z() * wrkvrt.vertexMom.Pz()) / (vDist.norm() * wrkvrt.vertexMom.Rho());
1871 float dphi_trk1 = 0.0;
1872 float dphi_trk2 = 0.0;
1873
1874 if (m_doDisappearingTrackVertexing){
1875 // Process track1
1876 const auto* track1 = trackChi2Pairs[0].first;
1877 dphi_trk1 = TVector2::Phi_mpi_pi(vDist.phi() - track1->phi());
1878 auto sv_perigee1 = m_trackToVertexTool->perigeeAtVertex(ctx, *track1, wrkvrt.vertex);
1879 if (sv_perigee1) {
1880 perigee_x_trk1 = sv_perigee1->position().x();
1881 perigee_y_trk1 = sv_perigee1->position().y();
1882 perigee_z_trk1 = sv_perigee1->position().z();
1883 perigee_px_trk1 = sv_perigee1->momentum().x();
1884 perigee_py_trk1 = sv_perigee1->momentum().y();
1885 perigee_pz_trk1 = sv_perigee1->momentum().z();
1886 perigee_cov_xx_trk1 = (*sv_perigee1->covariance())(0, 0);
1887 perigee_cov_xy_trk1 = (*sv_perigee1->covariance())(0, 1);
1888 perigee_cov_xz_trk1 = (*sv_perigee1->covariance())(0, 2);
1889 perigee_cov_yy_trk1 = (*sv_perigee1->covariance())(1, 1);
1890 perigee_cov_yz_trk1 = (*sv_perigee1->covariance())(1, 2);
1891 perigee_cov_zz_trk1 = (*sv_perigee1->covariance())(2, 2);
1892 perigee_d0_trk1 = sv_perigee1->parameters()[Trk::d0];
1893 perigee_z0_trk1 = sv_perigee1->parameters()[Trk::z0];
1894 perigee_qOverP_trk1 = sv_perigee1->parameters()[Trk::qOverP];
1895 perigee_theta_trk1 = sv_perigee1->parameters()[Trk::theta];
1896 perigee_phi_trk1 = sv_perigee1->parameters()[Trk::phi];
1897 perigee_charge_trk1 = sv_perigee1->parameters()[Trk::qOverP] > 0 ? 1 : -1;
1898 }else{
1899 ATH_MSG_DEBUG("Failed to obtain perigee for track1 at vertex.");
1900 }
1901
1902 //Process track2
1903 const auto* track2 = trackChi2Pairs[1].first;
1904 dphi_trk2 = TVector2::Phi_mpi_pi(vDist.phi() - track2->phi());
1905 auto sv_perigee2 = m_trackToVertexTool->perigeeAtVertex(ctx, *track2, wrkvrt.vertex);
1906 if (sv_perigee2) {
1907 perigee_x_trk2 = sv_perigee2->position().x();
1908 perigee_y_trk2 = sv_perigee2->position().y();
1909 perigee_z_trk2 = sv_perigee2->position().z();
1910 perigee_px_trk2 = sv_perigee2->momentum().x();
1911 perigee_py_trk2 = sv_perigee2->momentum().y();
1912 perigee_pz_trk2 = sv_perigee2->momentum().z();
1913 perigee_cov_xx_trk2 = (*sv_perigee2->covariance())(0, 0);
1914 perigee_cov_xy_trk2 = (*sv_perigee2->covariance())(0, 1);
1915 perigee_cov_xz_trk2 = (*sv_perigee2->covariance())(0, 2);
1916 perigee_cov_yy_trk2 = (*sv_perigee2->covariance())(1, 1);
1917 perigee_cov_yz_trk2 = (*sv_perigee2->covariance())(1, 2);
1918 perigee_cov_zz_trk2 = (*sv_perigee2->covariance())(2, 2);
1919 perigee_d0_trk2 = sv_perigee2->parameters()[Trk::d0];
1920 perigee_z0_trk2 = sv_perigee2->parameters()[Trk::z0];
1921 perigee_qOverP_trk2 = sv_perigee2->parameters()[Trk::qOverP];
1922 perigee_theta_trk2 = sv_perigee2->parameters()[Trk::theta];
1923 perigee_phi_trk2 = sv_perigee2->parameters()[Trk::phi];
1924 perigee_charge_trk2 = sv_perigee2->parameters()[Trk::qOverP] > 0 ? 1 : -1;
1925 }else{
1926 ATH_MSG_DEBUG("Failed to obtain perigee for track2 at vertex.");
1927 }
1928
1929 if(sv_perigee1 && sv_perigee2){
1930 perigee_distance = sqrt(
1931 (perigee_x_trk1 - perigee_x_trk2) * (perigee_x_trk1 - perigee_x_trk2) +
1932 (perigee_y_trk1 - perigee_y_trk2) * (perigee_y_trk1 - perigee_y_trk2) +
1933 (perigee_z_trk1 - perigee_z_trk2) * (perigee_z_trk1 - perigee_z_trk2)
1934 );
1935 }
1936 if(perigee_distance > m_twoTrVrtMaxPerigeeDist) continue;
1937 }
1938
1939 //
1940 // calculate opening angle between all 2-track pairs, and store the minimum
1941 //
1942 double minOpAng = AlgConsts::invalidFloat;
1943 std::vector<double> opAngles;
1944
1945 for( auto itr1 = tracks.begin(); itr1 != tracks.end(); ++itr1 ) {
1946 for( auto itr2 = std::next( itr1 ); itr2 != tracks.end(); ++itr2 ) {
1947 const auto& p1 = (*itr1)->p4().Vect();
1948 const auto& p2 = (*itr2)->p4().Vect();
1949 auto cos = p1 * p2 / p1.Mag() / p2.Mag();
1950 opAngles.emplace_back( cos );
1951 }
1952 }
1953 minOpAng = *( std::max_element( opAngles.begin(), opAngles.end() ) );
1954 if( m_FillNtuple ) m_ntupleVars->get< vector<double> >( "SecVtx_MinOpAng" ).emplace_back(minOpAng);
1955
1956
1957 if( m_FillHist ) m_hists["finalCutMonitor"]->Fill( 5 );
1958
1959 if( badIPflag ) {
1960 ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Bad impact parameter signif wrt SV was flagged." );
1961 }
1962
1963 if (m_doRemoveNonLeptonVertices) {
1964
1965 bool oneLepMatchTrack = false;
1966 for (const auto *trk: tracks) {
1967 if ( std::find(m_leptonicTracks.begin(), m_leptonicTracks.end(), trk) != m_leptonicTracks.end() ) {
1968 oneLepMatchTrack = true;
1969 break;
1970 }
1971 }
1972
1973 // If there are no tracks matched to leptons, do not save the container to the output.
1974 if (!oneLepMatchTrack) continue;
1975 }
1976
1978 // Data filling to xAOD container
1979
1980 wrkvrt.isGood = true;
1981
1982 // Firstly store the new vertex to the container before filling properties.
1983 // (This is the feature of xAOD.)
1984 xAOD::Vertex* vertex = new xAOD::Vertex;
1985 secondaryVertexContainer->emplace_back( vertex );
1986
1987 // Registering the vertex position to xAOD::Vertex
1988 vertex->setPosition( wrkvrt.vertex );
1989
1990 // Registering the vertex type: SV
1991 vertex->setVertexType( xAOD::VxType::SecVtx );
1992
1993 // Registering the vertex chi2 and Ndof
1994 // Here, we register the core chi2 of the core (before track association)
1995 vertex->setFitQuality( wrkvrt.Chi2_core, wrkvrt.ndof_core() );
1996
1997 // Registering the vertex covariance matrix
1998 std::vector<float> fCov(wrkvrt.vertexCov.cbegin(), wrkvrt.vertexCov.cend());
1999 vertex->setCovariance(fCov);
2000
2001 // Registering the vertex momentum and charge
2002 static const SG::Accessor<float> vtx_pxAcc("vtx_px");
2003 static const SG::Accessor<float> vtx_pyAcc("vtx_py");
2004 static const SG::Accessor<float> vtx_pzAcc("vtx_pz");
2005 static const SG::Accessor<float> vtx_massAcc("vtx_mass");
2006 static const SG::Accessor<float> vtx_chargeAcc("vtx_charge");
2007 static const SG::Accessor<float> chi2_coreAcc("chi2_core");
2008 static const SG::Accessor<float> ndof_coreAcc("ndof_core");
2009 static const SG::Accessor<float> chi2_assocAcc("chi2_assoc");
2010 static const SG::Accessor<float> ndof_assocAcc("ndof_assoc");
2011 static const SG::Accessor<float> massAcc("mass");
2012 static const SG::Accessor<float> mass_eAcc("mass_e");
2013 static const SG::Accessor<float> mass_selectedTracksAcc("mass_selectedTracks");
2014 static const SG::Accessor<float> minOpAngAcc("minOpAng");
2015 static const SG::Accessor<int> num_trksAcc("num_trks");
2016 static const SG::Accessor<int> num_selectedTracksAcc("num_selectedTracks");
2017 static const SG::Accessor<int> num_associatedTracksAcc("num_associatedTracks");
2018 static const SG::Accessor<float> dCloseVrtAcc("dCloseVrt");
2019
2020 vtx_pxAcc(*vertex) = wrkvrt.vertexMom.Px();
2021 vtx_pyAcc(*vertex) = wrkvrt.vertexMom.Py();
2022 vtx_pzAcc(*vertex) = wrkvrt.vertexMom.Pz();
2023
2024 vtx_massAcc(*vertex) = wrkvrt.vertexMom.M();
2025 vtx_chargeAcc(*vertex) = wrkvrt.Charge;
2026
2027 chi2_coreAcc(*vertex) = wrkvrt.Chi2_core;
2028 ndof_coreAcc(*vertex) = wrkvrt.ndof_core();
2029 chi2_assocAcc(*vertex) = wrkvrt.Chi2;
2030 ndof_assocAcc(*vertex) = wrkvrt.ndof();
2031 // Other SV properties
2032 massAcc(*vertex) = sumP4_pion.M();
2033 mass_eAcc(*vertex) = sumP4_electron.M();
2034 mass_selectedTracksAcc(*vertex) = sumP4_selected.M();
2035 minOpAngAcc(*vertex) = minOpAng;
2036 num_trksAcc(*vertex) = wrkvrt.nTracksTotal();
2037 num_selectedTracksAcc(*vertex) = wrkvrt.selectedTrackIndices.size();
2038 num_associatedTracksAcc(*vertex) = wrkvrt.associatedTrackIndices.size();
2039 dCloseVrtAcc(*vertex) = wrkvrt.closestWrkVrtValue;
2040
2041 // Registering the vertex momentum and charge
2042 if (m_doDisappearingTrackVertexing){
2043 static const SG::Accessor<float> perigee_x_trk1Acc("perigee_x_trk1");
2044 static const SG::Accessor<float> perigee_y_trk1Acc("perigee_y_trk1");
2045 static const SG::Accessor<float> perigee_z_trk1Acc("perigee_z_trk1");
2046 static const SG::Accessor<float> perigee_x_trk2Acc("perigee_x_trk2");
2047 static const SG::Accessor<float> perigee_y_trk2Acc("perigee_y_trk2");
2048 static const SG::Accessor<float> perigee_z_trk2Acc("perigee_z_trk2");
2049 static const SG::Accessor<float> perigee_px_trk1Acc("perigee_px_trk1");
2050 static const SG::Accessor<float> perigee_py_trk1Acc("perigee_py_trk1");
2051 static const SG::Accessor<float> perigee_pz_trk1Acc("perigee_pz_trk1");
2052 static const SG::Accessor<float> perigee_px_trk2Acc("perigee_px_trk2");
2053 static const SG::Accessor<float> perigee_py_trk2Acc("perigee_py_trk2");
2054 static const SG::Accessor<float> perigee_pz_trk2Acc("perigee_pz_trk2");
2055 static const SG::Accessor<float> perigee_cov_xx_trk1Acc("perigee_cov_xx_trk1");
2056 static const SG::Accessor<float> perigee_cov_xy_trk1Acc("perigee_cov_xy_trk1");
2057 static const SG::Accessor<float> perigee_cov_xz_trk1Acc("perigee_cov_xz_trk1");
2058 static const SG::Accessor<float> perigee_cov_yy_trk1Acc("perigee_cov_yy_trk1");
2059 static const SG::Accessor<float> perigee_cov_yz_trk1Acc("perigee_cov_yz_trk1");
2060 static const SG::Accessor<float> perigee_cov_zz_trk1Acc("perigee_cov_zz_trk1");
2061 static const SG::Accessor<float> perigee_cov_xx_trk2Acc("perigee_cov_xx_trk2");
2062 static const SG::Accessor<float> perigee_cov_xy_trk2Acc("perigee_cov_xy_trk2");
2063 static const SG::Accessor<float> perigee_cov_xz_trk2Acc("perigee_cov_xz_trk2");
2064 static const SG::Accessor<float> perigee_cov_yy_trk2Acc("perigee_cov_yy_trk2");
2065 static const SG::Accessor<float> perigee_cov_yz_trk2Acc("perigee_cov_yz_trk2");
2066 static const SG::Accessor<float> perigee_cov_zz_trk2Acc("perigee_cov_zz_trk2");
2067 static const SG::Accessor<float> perigee_d0_trk1Acc("perigee_d0_trk1");
2068 static const SG::Accessor<float> perigee_d0_trk2Acc("perigee_d0_trk2");
2069 static const SG::Accessor<float> perigee_z0_trk1Acc("perigee_z0_trk1");
2070 static const SG::Accessor<float> perigee_z0_trk2Acc("perigee_z0_trk2");
2071 static const SG::Accessor<float> perigee_qOverP_trk1Acc("perigee_qOverP_trk1");
2072 static const SG::Accessor<float> perigee_qOverP_trk2Acc("perigee_qOverP_trk2");
2073 static const SG::Accessor<float> perigee_theta_trk1Acc("perigee_theta_trk1");
2074 static const SG::Accessor<float> perigee_theta_trk2Acc("perigee_theta_trk2");
2075 static const SG::Accessor<float> perigee_phi_trk1Acc("perigee_phi_trk1");
2076 static const SG::Accessor<float> perigee_phi_trk2Acc("perigee_phi_trk2");
2077 static const SG::Accessor<int> perigee_charge_trk1Acc("perigee_charge_trk1");
2078 static const SG::Accessor<int> perigee_charge_trk2Acc("perigee_charge_trk2");
2079 static const SG::Accessor<float> vPosAcc("vPos");
2080 static const SG::Accessor<float> vPosMomAngTAcc("vPosMomAngT");
2081 static const SG::Accessor<float> vPosMomAng3DAcc("vPosMomAng3D");
2082 static const SG::Accessor<float> dphi_trk1Acc("dphi_trk1");
2083 static const SG::Accessor<float> dphi_trk2Acc("dphi_trk2");
2084 perigee_x_trk1Acc(*vertex) = perigee_x_trk1;
2085 perigee_y_trk1Acc(*vertex) = perigee_y_trk1;
2086 perigee_z_trk1Acc(*vertex) = perigee_z_trk1;
2087 perigee_x_trk2Acc(*vertex) = perigee_x_trk2;
2088 perigee_y_trk2Acc(*vertex) = perigee_y_trk2;
2089 perigee_z_trk2Acc(*vertex) = perigee_z_trk2;
2090 perigee_px_trk1Acc(*vertex) = perigee_px_trk1;
2091 perigee_py_trk1Acc(*vertex) = perigee_py_trk1;
2092 perigee_pz_trk1Acc(*vertex) = perigee_pz_trk1;
2093 perigee_px_trk2Acc(*vertex) = perigee_px_trk2;
2094 perigee_py_trk2Acc(*vertex) = perigee_py_trk2;
2095 perigee_pz_trk2Acc(*vertex) = perigee_pz_trk2;
2096 perigee_cov_xx_trk1Acc(*vertex) = perigee_cov_xx_trk1;
2097 perigee_cov_xy_trk1Acc(*vertex) = perigee_cov_xy_trk1;
2098 perigee_cov_xz_trk1Acc(*vertex) = perigee_cov_xz_trk1;
2099 perigee_cov_yy_trk1Acc(*vertex) = perigee_cov_yy_trk1;
2100 perigee_cov_yz_trk1Acc(*vertex) = perigee_cov_yz_trk1;
2101 perigee_cov_zz_trk1Acc(*vertex) = perigee_cov_zz_trk1;
2102 perigee_cov_xx_trk2Acc(*vertex) = perigee_cov_xx_trk2;
2103 perigee_cov_xy_trk2Acc(*vertex) = perigee_cov_xy_trk2;
2104 perigee_cov_xz_trk2Acc(*vertex) = perigee_cov_xz_trk2;
2105 perigee_cov_yy_trk2Acc(*vertex) = perigee_cov_yy_trk2;
2106 perigee_cov_yz_trk2Acc(*vertex) = perigee_cov_yz_trk2;
2107 perigee_cov_zz_trk2Acc(*vertex) = perigee_cov_zz_trk2;
2108 perigee_d0_trk1Acc(*vertex) = perigee_d0_trk1;
2109 perigee_d0_trk2Acc(*vertex) = perigee_d0_trk2;
2110 perigee_z0_trk1Acc(*vertex) = perigee_z0_trk1;
2111 perigee_z0_trk2Acc(*vertex) = perigee_z0_trk2;
2112 perigee_qOverP_trk1Acc(*vertex) = perigee_qOverP_trk1;
2113 perigee_qOverP_trk2Acc(*vertex) = perigee_qOverP_trk2;
2114 perigee_theta_trk1Acc(*vertex) = perigee_theta_trk1;
2115 perigee_theta_trk2Acc(*vertex) = perigee_theta_trk2;
2116 perigee_phi_trk1Acc(*vertex) = perigee_phi_trk1;
2117 perigee_phi_trk2Acc(*vertex) = perigee_phi_trk2;
2118 perigee_charge_trk1Acc(*vertex) = perigee_charge_trk1;
2119 perigee_charge_trk2Acc(*vertex) = perigee_charge_trk2;
2120 vPosAcc(*vertex) = vPos;
2121 vPosMomAngTAcc(*vertex) = vPosMomAngT;
2122 vPosMomAng3DAcc(*vertex) = vPosMomAng3D;
2123 dphi_trk1Acc(*vertex) = dphi_trk1;
2124 dphi_trk2Acc(*vertex) = dphi_trk2;
2125 }
2126
2127 // Registering tracks comprising the vertex to xAOD::Vertex
2128 // loop over the tracks comprising the vertex
2129 for( auto trk_id : wrkvrt.selectedTrackIndices ) {
2130
2131 const xAOD::TrackParticle *trk = m_selectedTracks.at( trk_id );
2132
2133 // Acquire link the track to the vertex
2134 ElementLink<xAOD::TrackParticleContainer> link_trk( *( dynamic_cast<const xAOD::TrackParticleContainer*>( trk->container() ) ), static_cast<long unsigned int>(trk->index()) );
2135
2136 // Register the link to the vertex
2137 vertex->addTrackAtVertex( link_trk, 1. );
2138
2139 }
2140
2141 for( auto trk_id : wrkvrt.associatedTrackIndices ) {
2142
2143 const xAOD::TrackParticle *trk = m_associatedTracks.at( trk_id );
2144
2145 // Acquire link the track to the vertex
2146 ElementLink<xAOD::TrackParticleContainer> link_trk( *( dynamic_cast<const xAOD::TrackParticleContainer*>( trk->container() ) ), static_cast<long unsigned int>(trk->index()) );
2147
2148 // Register the link to the vertex
2149 vertex->addTrackAtVertex( link_trk, 1. );
2150
2151 }
2152
2153
2154 if( m_doMapToLocal ) {
2155 // Obtain the local mapping of the reconstructed vertex
2156 Trk::MappedVertex mappedVtx = m_vertexMapper->mapToLocal( wrkvrt.vertex );
2157 static const SG::Accessor<int> local_identifierHashAcc("local_identifierHash");
2158 static const SG::Accessor<int> local_layerIndexAcc("local_layerIndex");
2159 static const SG::Accessor<float> local_posXAcc("local_posX");
2160 static const SG::Accessor<float> local_posYAcc("local_posY");
2161 static const SG::Accessor<float> local_posZAcc("local_posZ");
2162 if( mappedVtx.valid ) {
2163 local_identifierHashAcc(*vertex) = mappedVtx.identifierHash;
2164 local_layerIndexAcc(*vertex) = mappedVtx.layerIndex;
2165 local_posXAcc(*vertex) = mappedVtx.localPosition.x();
2166 local_posYAcc(*vertex) = mappedVtx.localPosition.y();
2167 local_posZAcc(*vertex) = mappedVtx.localPosition.z();
2168 } else {
2169 local_identifierHashAcc(*vertex) = AlgConsts::invalidInt;
2170 local_layerIndexAcc(*vertex) = AlgConsts::invalidInt;
2171 local_posXAcc(*vertex) = AlgConsts::invalidFloat;
2172 local_posYAcc(*vertex) = AlgConsts::invalidFloat;
2173 local_posZAcc(*vertex) = AlgConsts::invalidFloat;
2174 }
2175 }
2176
2177
2178 // For MC, try to trace down to the truth particles,
2179 // and depending on the topology, categorize the label of the reconstructed vertex.
2180 if( m_doTruth ) {
2181 ATH_CHECK( categorizeVertexTruthTopology( vertex ) );
2182 }
2183
2184 // Keep the link between wrkvrt and vertex for later use
2185 wrkvrtLinkMap[&wrkvrt] = vertex;
2186
2187
2188 } // loop over vertices
2189
2190 if( m_FillNtuple ) {
2191 ATH_CHECK( fillAANT_SecondaryVertices( secondaryVertexContainer ) );
2192 }
2193
2194
2195 // Post process -- Additional augmentations
2196 if( m_doAugmentDVimpactParametersToMuons ) { ATH_CHECK( augmentDVimpactParametersToLeptons<xAOD::Muon> ( "Muons" ) ); }
2197 if( m_doAugmentDVimpactParametersToElectrons ) { ATH_CHECK( augmentDVimpactParametersToLeptons<xAOD::Electron>( "Electrons" ) ); }
2198
2199 } catch (const std::out_of_range& e) {
2200
2201 ATH_MSG_WARNING( " > " << __FUNCTION__ << ": out of range error is detected: " << e.what() );
2202
2203 return StatusCode::SUCCESS;
2204
2205 } catch( ... ) {
2206
2207 ATH_MSG_WARNING( " > " << __FUNCTION__ << ": some other error is detected." );
2208
2209 return StatusCode::SUCCESS;
2210
2211 }
2212
2213 return StatusCode::SUCCESS;
2214 }
2215
2216
2217 //____________________________________________________________________________________________________
2218 StatusCode VrtSecInclusive::monitorVertexingAlgorithmStep( const EventContext& ctx,
2219 std::vector<WrkVrt>* workVerticesContainer, const std::string& name, bool final ) {
2220 if( m_FillIntermediateVertices ) {
2221
2222 SG::WriteHandle<xAOD::VertexContainer> vertexHandle = SG::makeHandle( m_intermediateVertexKey[name], ctx );
2223 ATH_CHECK( vertexHandle.record(std::make_unique<xAOD::VertexContainer>(),
2224 std::make_unique<xAOD::VertexAuxContainer>()) );
2225 xAOD::VertexContainer* intermediateVertexContainer = vertexHandle.ptr();
2226 m_vertexCollectionsDefinitions[ m_intermediateVertexKey.at(name).key() ] = true;
2227
2228 for( auto& wrkvrt : *workVerticesContainer ) {
2229
2230 xAOD::Vertex* vertex = new xAOD::Vertex;
2231 intermediateVertexContainer->emplace_back( vertex );
2232
2233 // Registering the vertex position to xAOD::Vertex
2234 vertex->setPosition( wrkvrt.vertex );
2235
2236 // Registering the vertex type: SV
2237 vertex->setVertexType( xAOD::VxType::SecVtx );
2238
2239 // Registering the vertex chi2 and Ndof
2240 int ndof = wrkvrt.ndof();
2241 vertex->setFitQuality( wrkvrt.Chi2, ndof );
2242
2243 // Registering the vertex covariance matrix
2244 std::vector<float> fCov(wrkvrt.vertexCov.cbegin(), wrkvrt.vertexCov.cend());
2245 vertex->setCovariance(fCov);
2246
2247 // Registering tracks comprising the vertex to xAOD::Vertex
2248 // loop over the tracks comprising the vertex
2249 for( auto trk_id : wrkvrt.selectedTrackIndices ) {
2250
2251 const xAOD::TrackParticle *trk = m_selectedTracks.at( trk_id );
2252
2253 // Acquire link the track to the vertex
2254 ElementLink<xAOD::TrackParticleContainer> link_trk( *( dynamic_cast<const xAOD::TrackParticleContainer*>( trk->container() ) ), static_cast<long unsigned int>(trk->index()) );
2255
2256 // Register the link to the vertex
2257 vertex->addTrackAtVertex( link_trk, 1. );
2258
2259 }
2260
2261 for( auto trk_id : wrkvrt.associatedTrackIndices ) {
2262
2263 const xAOD::TrackParticle *trk = m_associatedTracks.at( trk_id );
2264
2265 // Acquire link the track to the vertex
2266 ElementLink<xAOD::TrackParticleContainer> link_trk( *( dynamic_cast<const xAOD::TrackParticleContainer*>( trk->container() ) ), static_cast<long unsigned int>(trk->index()) );
2267
2268 // Register the link to the vertex
2269 vertex->addTrackAtVertex( link_trk, 1. );
2270
2271 }
2272 }
2273
2274 }
2275
2276
2277
2278 if( !m_FillHist ) return StatusCode::SUCCESS;
2279
2280 printWrkSet( workVerticesContainer, Form("%s (step %u)", name.c_str(), m_vertexingAlgorithmStep) );
2281
2282 unsigned count = std::count_if( workVerticesContainer->begin(), workVerticesContainer->end(),
2283 []( WrkVrt& v ) { return ( v.selectedTrackIndices.size() + v.associatedTrackIndices.size() ) >= 2; } );
2284
2285 if( m_vertexingAlgorithmStep == 0 ) {
2286
2287 const auto compSize = m_selectedTracks.size()*(m_selectedTracks.size() - 1)/2 - m_incomp.size();
2288 m_hists["vertexYield"]->Fill( m_vertexingAlgorithmStep, compSize );
2289
2290 } else {
2291
2292 m_hists["vertexYield"]->Fill( m_vertexingAlgorithmStep, count );
2293
2294 }
2295
2296 m_hists["vertexYield"]->GetXaxis()->SetBinLabel( m_vertexingAlgorithmStep+1, name.c_str() );
2297
2298 for( auto& vertex : *workVerticesContainer ) {
2299 auto ntrk = vertex.selectedTrackIndices.size() + vertex.associatedTrackIndices.size();
2300 if( vertex.isGood && ntrk >= 2 ) {
2301 dynamic_cast<TH2F*>( m_hists["vertexYieldNtrk"] )->Fill( ntrk, m_vertexingAlgorithmStep );
2302 dynamic_cast<TH2F*>( m_hists["vertexYieldChi2"] )->Fill( vertex.Chi2/(vertex.ndof() + AlgConsts::infinitesimal), m_vertexingAlgorithmStep );
2303 }
2304 }
2305 m_hists["vertexYieldNtrk"]->GetYaxis()->SetBinLabel( m_vertexingAlgorithmStep+1, name.c_str() );
2306 m_hists["vertexYieldChi2"]->GetYaxis()->SetBinLabel( m_vertexingAlgorithmStep+1, name.c_str() );
2307
2308
2309 if( !final ) return StatusCode::SUCCESS;
2310
2311 for( auto& vertex : *workVerticesContainer ) {
2312 auto ntrk = vertex.selectedTrackIndices.size() + vertex.associatedTrackIndices.size();
2313 if( vertex.isGood && ntrk >= 2 ) {
2314 m_hists["finalVtxNtrk"] ->Fill( ntrk );
2315 m_hists["finalVtxR"] ->Fill( vertex.vertex.perp() );
2316 dynamic_cast<TH2F*>( m_hists["finalVtxNtrkR"] )->Fill( ntrk, vertex.vertex.perp() );
2317 }
2318 }
2319
2320 return StatusCode::SUCCESS;
2321 }
2322
2323 //____________________________________________________________________________________________________
2324 bool VrtSecInclusive::getSVImpactParameters(const xAOD::TrackParticle* trk, const Amg::Vector3D& vertex,
2325 std::vector<double>& impactParameters,
2326 std::vector<double>& impactParErrors){
2327
2328 impactParameters.clear();
2329 impactParErrors.clear();
2330
2331 if( m_trkExtrapolator==1 ){
2332 m_fitSvc->VKalGetImpact(trk, vertex, static_cast<int>( trk->charge() ), impactParameters, impactParErrors);
2333 }
2334 else if( m_trkExtrapolator==2 ){
2335 auto sv_perigee = m_trackToVertexTool->perigeeAtVertex(Gaudi::Hive::currentContext(), *trk, vertex );
2336 if( !sv_perigee ) return false;
2337 impactParameters.push_back(sv_perigee->parameters() [Trk::d0]);
2338 impactParameters.push_back(sv_perigee->parameters() [Trk::z0]);
2339 impactParErrors.push_back((*sv_perigee->covariance())( Trk::d0, Trk::d0 ));
2340 impactParErrors.push_back((*sv_perigee->covariance())( Trk::z0, Trk::z0 ));
2341 }
2342 else{
2343 ATH_MSG_WARNING( " > " << __FUNCTION__ << ": Unknown track extrapolator " << m_trkExtrapolator );
2344 return false;
2345 }
2346
2347 return true;
2348
2349 } // getSVImpactParameters
2350
2351
2352
2353} // end of namespace VKalVrtAthena
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition AthMsgStreamMacros.h:31
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition AthMsgStreamMacros.h:28
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
Accessor.h
Helper class to provide type-safe access to aux data.
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition EventPrimitives.h:50
sc
static Double_t sc
Definition LArPhysWaveHECTool.cxx:37
AthCommonMsg< Algorithm >::msg
MsgStream & msg() const
Definition AthCommonMsg.h:24
DataVector::emplace_back
value_type emplace_back(value_type pElem)
Add an element to the end of the collection.
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
SG::Accessor
Helper class to provide type-safe access to aux data.
Definition Control/AthContainers/AthContainers/Accessor.h:68
SG::AuxElement::Decorator
SG::Decorator< T, ALLOC > Decorator
Definition AuxElement.h:576
SG::AuxElement::container
const SG::AuxVectorData * container() const
Return the container holding this element.
SG::AuxElement::index
size_t index() const
Return the index of this element within its container.
SG::ConstAccessor
Helper class to provide constant type-safe access to aux data.
Definition ConstAccessor.h:55
SG::ConstAccessor::isAvailable
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
SG::ReadHandle::cptr
const_pointer_type cptr()
Dereference the pointer.
SG::WriteHandle::record
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
SG::WriteHandle::ptr
pointer_type ptr()
Dereference the pointer.
VKalVrtAthena::VrtSecInclusive::m_doRemoveNonLeptonVertices
Gaudi::Property< bool > m_doRemoveNonLeptonVertices
Definition VrtSecInclusive.h:228
VKalVrtAthena::VrtSecInclusive::m_pvCompatibilityCut
Gaudi::Property< double > m_pvCompatibilityCut
Definition VrtSecInclusive.h:192
VKalVrtAthena::VrtSecInclusive::m_associateChi2Cut
Gaudi::Property< double > m_associateChi2Cut
Definition VrtSecInclusive.h:214
VKalVrtAthena::VrtSecInclusive::m_doWildMerging
Gaudi::Property< bool > m_doWildMerging
Definition VrtSecInclusive.h:188
VKalVrtAthena::VrtSecInclusive::m_SelVrtChi2Cut
Gaudi::Property< double > m_SelVrtChi2Cut
Definition VrtSecInclusive.h:193
VKalVrtAthena::VrtSecInclusive::refitVertexWithSuggestion
StatusCode refitVertexWithSuggestion(WrkVrt &, const Amg::Vector3D &)
refit the vertex with suggestion
Definition Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:585
VKalVrtAthena::VrtSecInclusive::passedFakeReject
bool passedFakeReject(const Amg::Vector3D &FitVertex, const xAOD::TrackParticle *itrk, const xAOD::TrackParticle *jtrk)
Flag false if the consistituent tracks are not consistent with the vertex position.
Definition Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:2184
VKalVrtAthena::VrtSecInclusive::m_trkExtrapolator
Gaudi::Property< int > m_trkExtrapolator
Definition VrtSecInclusive.h:177
VKalVrtAthena::VrtSecInclusive::printWrkSet
void printWrkSet(const std::vector< WrkVrt > *WrkVrtSet, const std::string &name)
print the contents of reconstructed vertices
Definition Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:862
VKalVrtAthena::VrtSecInclusive::m_twoTrVrtMaxPerigeeDist
Gaudi::Property< double > m_twoTrVrtMaxPerigeeDist
Definition VrtSecInclusive.h:232
VKalVrtAthena::VrtSecInclusive::m_doTightPVcompatibilityCut
Gaudi::Property< bool > m_doTightPVcompatibilityCut
Definition VrtSecInclusive.h:181
VKalVrtAthena::VrtSecInclusive::m_vertexKey
SG::WriteHandleKey< xAOD::VertexContainer > m_vertexKey
Definition VrtSecInclusive.h:299
VKalVrtAthena::VrtSecInclusive::m_reassembleMaxImpactParameterZ0
Gaudi::Property< double > m_reassembleMaxImpactParameterZ0
Definition VrtSecInclusive.h:217
VKalVrtAthena::VrtSecInclusive::m_mergeByShufflingAllowance
Gaudi::Property< double > m_mergeByShufflingAllowance
Definition VrtSecInclusive.h:219
VKalVrtAthena::VrtSecInclusive::trackClassification
void trackClassification(std::vector< WrkVrt > *, std::map< long int, std::vector< long int > > &)
Definition Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:775
VKalVrtAthena::VrtSecInclusive::m_doPVcompatibilityCut
Gaudi::Property< bool > m_doPVcompatibilityCut
Definition VrtSecInclusive.h:180
VKalVrtAthena::VrtSecInclusive::m_trkDecors
std::map< unsigned, SG::Decorator< float > > m_trkDecors
Definition VrtSecInclusive.h:295
VKalVrtAthena::VrtSecInclusive::m_vertexMapper
PublicToolHandle< Trk::IVertexMapper > m_vertexMapper
Definition VrtSecInclusive.h:277
VKalVrtAthena::VrtSecInclusive::findWorstChi2ofMaximallySharedTrack
double findWorstChi2ofMaximallySharedTrack(std::vector< WrkVrt > *, std::map< long int, std::vector< long int > > &, long int &, long int &)
Definition Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:807
VKalVrtAthena::VrtSecInclusive::m_associateMaxD0Signif
Gaudi::Property< double > m_associateMaxD0Signif
Definition VrtSecInclusive.h:211
VKalVrtAthena::VrtSecInclusive::m_ntupleVars
std::unique_ptr< NtupleVars > m_ntupleVars
Definition VrtSecInclusive.h:318
VKalVrtAthena::VrtSecInclusive::m_doFastMode
Gaudi::Property< bool > m_doFastMode
Definition VrtSecInclusive.h:135
VKalVrtAthena::VrtSecInclusive::m_truncateWrkVertices
Gaudi::Property< bool > m_truncateWrkVertices
Definition VrtSecInclusive.h:207
VKalVrtAthena::VrtSecInclusive::m_associateMinDistanceToPV
Gaudi::Property< double > m_associateMinDistanceToPV
Definition VrtSecInclusive.h:210
VKalVrtAthena::VrtSecInclusive::getSVImpactParameters
bool getSVImpactParameters(const xAOD::TrackParticle *trk, const Amg::Vector3D &vertex, std::vector< double > &impactParameters, std::vector< double > &impactParErrors)
get secondary vertex impact parameters
Definition VertexingAlgs.cxx:2324
VKalVrtAthena::VrtSecInclusive::categorizeVertexTruthTopology
StatusCode categorizeVertexTruthTopology(xAOD::Vertex *vertex)
Definition TruthAlgs.cxx:44
VKalVrtAthena::VrtSecInclusive::m_vertexCollectionsDefinitions
std::unordered_map< std::string, bool > m_vertexCollectionsDefinitions
Definition VrtSecInclusive.h:583
VKalVrtAthena::VrtSecInclusive::removeInconsistentTracks
void removeInconsistentTracks(WrkVrt &)
Remove inconsistent tracks from vertices.
Definition Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:2198
VKalVrtAthena::VrtSecInclusive::m_doSuggestedRefitOnMerging
Gaudi::Property< bool > m_doSuggestedRefitOnMerging
Definition VrtSecInclusive.h:186
VKalVrtAthena::VrtSecInclusive::m_twoTrVrtMinDistFromPV
Gaudi::Property< double > m_twoTrVrtMinDistFromPV
Definition VrtSecInclusive.h:201
VKalVrtAthena::VrtSecInclusive::associateNonSelectedTracks
StatusCode associateNonSelectedTracks(const EventContext &ctx, std::vector< WrkVrt > *)
in addition to selected tracks, associate as much tracks as possible
Definition VertexingAlgs.cxx:1107
VKalVrtAthena::VrtSecInclusive::m_twoTrVrtMinRadius
Gaudi::Property< double > m_twoTrVrtMinRadius
Definition VrtSecInclusive.h:233
VKalVrtAthena::VrtSecInclusive::m_VertexMergeFinalDistCut
Gaudi::Property< double > m_VertexMergeFinalDistCut
Definition VrtSecInclusive.h:194
VKalVrtAthena::VrtSecInclusive::m_associatePtCut
Gaudi::Property< double > m_associatePtCut
Definition VrtSecInclusive.h:213
VKalVrtAthena::VrtSecInclusive::m_FillHist
Gaudi::Property< bool > m_FillHist
Definition VrtSecInclusive.h:126
VKalVrtAthena::VrtSecInclusive::checkTrackHitPatternToVertexOuterOnly
bool checkTrackHitPatternToVertexOuterOnly(const xAOD::TrackParticle *trk, const Amg::Vector3D &vertex)
A classical method with hard-coded geometry.
Definition Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:2057
VKalVrtAthena::VrtSecInclusive::m_associateMaxZ0Signif
Gaudi::Property< double > m_associateMaxZ0Signif
Definition VrtSecInclusive.h:212
VKalVrtAthena::VrtSecInclusive::m_TrackDetachCut
Gaudi::Property< double > m_TrackDetachCut
Definition VrtSecInclusive.h:197
VKalVrtAthena::VrtSecInclusive::m_TrackLocation
SG::ReadHandleKey< xAOD::TrackParticleContainer > m_TrackLocation
Definition VrtSecInclusive.h:112
VKalVrtAthena::VrtSecInclusive::m_hists
std::map< std::string, TH1 * > m_hists
Definition VrtSecInclusive.h:321
VKalVrtAthena::VrtSecInclusive::fillTrackSummary
static void fillTrackSummary(track_summary &summary, const xAOD::TrackParticle *trk)
retrieve the track hit information
Definition Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:954
VKalVrtAthena::VrtSecInclusive::improveVertexChi2
double improveVertexChi2(WrkVrt &)
attempt to improve the vertex chi2 by removing the most-outlier track one by one until the vertex chi...
Definition Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:265
VKalVrtAthena::VrtSecInclusive::distanceBetweenVertices
double distanceBetweenVertices(const WrkVrt &, const WrkVrt &) const
calculate the physical distance
Definition Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:112
VKalVrtAthena::VrtSecInclusive::reassembleVertices
StatusCode reassembleVertices(const EventContext &ctx, std::vector< WrkVrt > *)
attempt to merge vertices when all tracks of a vertex A is close to vertex B in terms of impact param...
Definition VertexingAlgs.cxx:985
VKalVrtAthena::VrtSecInclusive::mergeFinalVertices
StatusCode mergeFinalVertices(const EventContext &ctx, std::vector< WrkVrt > *)
attempt to merge vertices by lookng at the distance between two vertices
Definition VertexingAlgs.cxx:1467
VKalVrtAthena::VrtSecInclusive::nTrkCommon
static size_t nTrkCommon(std::vector< WrkVrt > *WrkVrtSet, const std::pair< unsigned, unsigned > &pairIndex)
returns the number of tracks commonly present in both vertices
Definition Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:662
VKalVrtAthena::VrtSecInclusive::m_FillNtuple
Gaudi::Property< bool > m_FillNtuple
Definition VrtSecInclusive.h:127
VKalVrtAthena::VrtSecInclusive::m_decor_is_svtrk_final
std::optional< SG::Decorator< char > > m_decor_is_svtrk_final
Definition VrtSecInclusive.h:294
VKalVrtAthena::VrtSecInclusive::findMinVerticesPair
double findMinVerticesPair(std::vector< WrkVrt > *, std::pair< unsigned, unsigned > &, const AlgForVerticesPair &)
returns the pair of vertices that give minimum in terms of some observable (e.g.
Definition Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:374
VKalVrtAthena::VrtSecInclusive::m_doFinalImproveChi2
Gaudi::Property< bool > m_doFinalImproveChi2
Definition VrtSecInclusive.h:191
VKalVrtAthena::VrtSecInclusive::m_doMagnetMerging
Gaudi::Property< bool > m_doMagnetMerging
Definition VrtSecInclusive.h:187
VKalVrtAthena::VrtSecInclusive::m_VertexMergeFinalDistScaling
Gaudi::Property< double > m_VertexMergeFinalDistScaling
Definition VrtSecInclusive.h:195
VKalVrtAthena::VrtSecInclusive::m_twoTrkVtxFormingD0Cut
Gaudi::Property< double > m_twoTrkVtxFormingD0Cut
Definition VrtSecInclusive.h:159
VKalVrtAthena::VrtSecInclusive::m_improveChi2ProbThreshold
Gaudi::Property< double > m_improveChi2ProbThreshold
Definition VrtSecInclusive.h:221
VKalVrtAthena::VrtSecInclusive::refitAndSelectGoodQualityVertices
StatusCode refitAndSelectGoodQualityVertices(const EventContext &ctx, std::vector< WrkVrt > *)
finalization of the vertex and store to xAOD::VertexContainer
Definition VertexingAlgs.cxx:1524
VKalVrtAthena::VrtSecInclusive::m_doTwoTrSoftBtag
Gaudi::Property< bool > m_doTwoTrSoftBtag
Definition VrtSecInclusive.h:199
VKalVrtAthena::VrtSecInclusive::mergeVertices
StatusCode mergeVertices(WrkVrt &destination, WrkVrt &source)
the 2nd vertex is merged into the 1st vertex.
Definition Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:452
VKalVrtAthena::VrtSecInclusive::m_mergeByShufflingMaxSignificance
Gaudi::Property< double > m_mergeByShufflingMaxSignificance
Definition VrtSecInclusive.h:218
VKalVrtAthena::VrtSecInclusive::m_associatedTracks
std::vector< const xAOD::TrackParticle * > m_associatedTracks
Definition VrtSecInclusive.h:261
VKalVrtAthena::VrtSecInclusive::m_trackToVertexTool
PublicToolHandle< Reco::ITrackToVertex > m_trackToVertexTool
get a handle on the Track to Vertex tool
Definition VrtSecInclusive.h:274
VKalVrtAthena::VrtSecInclusive::m_fitSvc
ToolHandle< Trk::ITrkVKalVrtFitter > m_fitSvc
Definition VrtSecInclusive.h:270
VKalVrtAthena::VrtSecInclusive::m_PrimVrtLocation
SG::ReadHandleKey< xAOD::VertexContainer > m_PrimVrtLocation
Definition VrtSecInclusive.h:115
VKalVrtAthena::VrtSecInclusive::m_selectedTracks
std::vector< const xAOD::TrackParticle * > m_selectedTracks
Definition VrtSecInclusive.h:260
VKalVrtAthena::VrtSecInclusive::m_maxWrkVertices
Gaudi::Property< size_t > m_maxWrkVertices
Definition VrtSecInclusive.h:208
VKalVrtAthena::VrtSecInclusive::rearrangeTracks
StatusCode rearrangeTracks(const EventContext &ctx, std::vector< WrkVrt > *)
Definition VertexingAlgs.cxx:758
VKalVrtAthena::VrtSecInclusive::m_doMapToLocal
Gaudi::Property< bool > m_doMapToLocal
Definition VrtSecInclusive.h:130
VKalVrtAthena::VrtSecInclusive::m_removeFakeVrtLate
Gaudi::Property< bool > m_removeFakeVrtLate
Definition VrtSecInclusive.h:183
VKalVrtAthena::VrtSecInclusive::m_reassembleMaxImpactParameterD0
Gaudi::Property< double > m_reassembleMaxImpactParameterD0
Definition VrtSecInclusive.h:216
VKalVrtAthena::VrtSecInclusive::m_removeFakeVrt
Gaudi::Property< bool > m_removeFakeVrt
Definition VrtSecInclusive.h:182
VKalVrtAthena::VrtSecInclusive::m_augVerString
Gaudi::Property< std::string > m_augVerString
Definition VrtSecInclusive.h:118
VKalVrtAthena::VrtSecInclusive::m_doDisappearingTrackVertexing
Gaudi::Property< bool > m_doDisappearingTrackVertexing
Definition VrtSecInclusive.h:231
VKalVrtAthena::VrtSecInclusive::track_summary
struct VKalVrtAthena::VrtSecInclusive::track_summary_properties track_summary
VKalVrtAthena::VrtSecInclusive::m_leptonicTracks
std::vector< const xAOD::TrackParticle * > m_leptonicTracks
Definition VrtSecInclusive.h:262
VKalVrtAthena::VrtSecInclusive::monitorVertexingAlgorithmStep
StatusCode monitorVertexingAlgorithmStep(const EventContext &ctx, std::vector< WrkVrt > *, const std::string &name, bool final=false)
monitor the intermediate status of vertexing
Definition VertexingAlgs.cxx:2218
VKalVrtAthena::VrtSecInclusive::m_incomp
std::vector< std::pair< int, int > > m_incomp
Definition VrtSecInclusive.h:374
VKalVrtAthena::VrtSecInclusive::m_FillIntermediateVertices
Gaudi::Property< bool > m_FillIntermediateVertices
Definition VrtSecInclusive.h:128
VKalVrtAthena::VrtSecInclusive::mergeByShuffling
StatusCode mergeByShuffling(const EventContext &ctx, std::vector< WrkVrt > *)
attempt to merge splitted vertices when they are significantly distant due to the long-tail behavior ...
Definition VertexingAlgs.cxx:1317
VKalVrtAthena::VrtSecInclusive::m_intermediateVertexKey
std::map< std::string, SG::WriteHandleKey< xAOD::VertexContainer > > m_intermediateVertexKey
Definition VrtSecInclusive.h:301
VKalVrtAthena::VrtSecInclusive::m_decor_isAssociated
std::optional< SG::Decorator< char > > m_decor_isAssociated
Definition VrtSecInclusive.h:293
VKalVrtAthena::VrtSecInclusive::m_tracingTruthVertices
std::vector< const xAOD::TruthVertex * > m_tracingTruthVertices
Definition VrtSecInclusive.h:567
VKalVrtAthena::VrtSecInclusive::extractIncompatibleTrackPairs
StatusCode extractIncompatibleTrackPairs(const EventContext &ctx, std::vector< WrkVrt > *)
related to the graph method and verte finding
Definition VertexingAlgs.cxx:39
VKalVrtAthena::VrtSecInclusive::m_doTruth
Gaudi::Property< bool > m_doTruth
Definition VrtSecInclusive.h:125
VKalVrtAthena::VrtSecInclusive::m_twoTrksVertexKey
SG::WriteHandleKey< xAOD::VertexContainer > m_twoTrksVertexKey
Definition VrtSecInclusive.h:300
VKalVrtAthena::VrtSecInclusive::findNtrackVertices
StatusCode findNtrackVertices(const EventContext &ctx, std::vector< WrkVrt > *)
Definition VertexingAlgs.cxx:395
VKalVrtAthena::VrtSecInclusive::m_twoTrVrtAngleCut
Gaudi::Property< double > m_twoTrVrtAngleCut
Definition VrtSecInclusive.h:200
VKalVrtAthena::VrtSecInclusive::significanceBetweenVertices
double significanceBetweenVertices(const WrkVrt &, const WrkVrt &) const
calculate the significance (Mahalanobis distance) between two reconstructed vertices
Definition Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:61
xAOD::TrackParticle_v1::charge
float charge() const
Returns the charge.
Definition TrackParticle_v1.cxx:143
count
int count(std::string s, const std::string &regx)
count how many occurances of a regx are in a string
Definition hcg.cxx:146
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
SG::makeHandle
SG::ReadCondHandle< T > makeHandle(const SG::ReadCondHandleKey< T > &key, const EventContext &ctx=Gaudi::Hive::currentContext())
Definition ReadCondHandle.h:269
Trk::theta
@ theta
Definition ParamDefs.h:66
Trk::qOverP
@ qOverP
perigee
Definition ParamDefs.h:67
Trk::phi
@ phi
Definition ParamDefs.h:75
Trk::d0
@ d0
Definition ParamDefs.h:63
Trk::z0
@ z0
Definition ParamDefs.h:64
VKalVrtAthena::vtxVtxDistance
double vtxVtxDistance(const Amg::Vector3D &v1, const Amg::Vector3D &v2)
Definition Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:55
index
Definition index.py:1
std
STL namespace.
std::sort
void sort(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end)
Specialization of sort for DataVector/List.
Definition DVL_algorithms.h:554
std::remove_if
DataModel_detail::iterator< DVL > remove_if(typename DataModel_detail::iterator< DVL > beg, typename DataModel_detail::iterator< DVL > end, Predicate pred)
Specialization of remove_if for DataVector/List.
Definition DVL_algorithms.h:70
xAOD::VxType::SecVtx
@ SecVtx
Secondary vertex.
Definition TrackingPrimitives.h:573
xAOD::TrackParticle
TrackParticle_v1 TrackParticle
Reference the current persistent version:
Definition Event/xAOD/xAODTracking/xAODTracking/TrackParticle.h:13
xAOD::VertexContainer
VertexContainer_v1 VertexContainer
Definition of the current "Vertex container version".
Definition VertexContainer.h:14
xAOD::Vertex
Vertex_v1 Vertex
Define the latest version of the vertex class.
Definition Event/xAOD/xAODTracking/xAODTracking/Vertex.h:16
xAOD::TrackParticleContainer
TrackParticleContainer_v1 TrackParticleContainer
Definition of the current "TrackParticle container version".
Definition Event/xAOD/xAODTracking/xAODTracking/TrackParticleContainer.h:14
xAOD::numberOfInnermostPixelLayerHits
@ numberOfInnermostPixelLayerHits
these are the hits in the 0th pixel barrel layer
Definition TrackingPrimitives.h:238
xAOD::numberOfPixelHits
@ numberOfPixelHits
these are the pixel hits, including the b-layer [unit8_t].
Definition TrackingPrimitives.h:260
Trk::MappedVertex::localPosition
Amg::Vector3D localPosition
Definition IVertexMapper.h:25
VKalVrtAthena::VrtSecInclusive::WrkVrt::Chi2
double Chi2
VKalVrt fit covariance.
Definition VrtSecInclusive.h:349
VKalVrtAthena::VrtSecInclusive::WrkVrt::vertexCov
std::vector< double > vertexCov
VKalVrt fit vertex 4-momentum.
Definition VrtSecInclusive.h:348
VKalVrtAthena::VrtSecInclusive::WrkVrt::Chi2PerTrk
std::vector< double > Chi2PerTrk
VKalVrt fit chi2 result.
Definition VrtSecInclusive.h:351
VKalVrtAthena::VrtSecInclusive::WrkVrt::selectedTrackIndices
std::deque< long int > selectedTrackIndices
flagged true for good vertex candidates
Definition VrtSecInclusive.h:344
VKalVrtAthena::VrtSecInclusive::WrkVrt::associatedTrackIndices
std::deque< long int > associatedTrackIndices
list if indices in TrackParticleContainer for selectedBaseTracks
Definition VrtSecInclusive.h:345
VKalVrtAthena::VrtSecInclusive::WrkVrt::Charge
long int Charge
list of VKalVrt fit chi2 for each track
Definition VrtSecInclusive.h:352
VKalVrtAthena::VrtSecInclusive::WrkVrt::vertex
Amg::Vector3D vertex
list if indices in TrackParticleContainer for associatedTracks
Definition VrtSecInclusive.h:346
VKalVrtAthena::VrtSecInclusive::WrkVrt::vertexMom
TLorentzVector vertexMom
VKalVrt fit vertex position.
Definition VrtSecInclusive.h:347
VKalVrtAthena::VrtSecInclusive::WrkVrt::closestWrkVrtValue
double closestWrkVrtValue
stores the index of the closest WrkVrt in std::vector<WrkVrt>
Definition VrtSecInclusive.h:355
VKalVrtAthena::VrtSecInclusive::WrkVrt::TrkAtVrt
std::vector< std::vector< double > > TrkAtVrt
total charge of the vertex
Definition VrtSecInclusive.h:353
VKalVrtAthena::VrtSecInclusive::WrkVrt::closestWrkVrtIndex
unsigned long closestWrkVrtIndex
list of track parameters wrt the reconstructed vertex
Definition VrtSecInclusive.h:354
result
Generated on for ATLAS Offline Software by doxygen 1.14.0