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