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