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