ATLAS Offline Software
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 
31 using namespace std;
32 
33 
34 
35 namespace 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
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
VKalVrtAthena::VrtSecInclusive::WrkVrt::vertex
Amg::Vector3D vertex
list if indices in TrackParticleContainer for associatedTracks
Definition: VrtSecInclusive.h:344
VKalVrtAthena::AlgConsts::infinitesimal
constexpr double infinitesimal
Definition: Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/Constants.h:20
IDTPM::ndof
float ndof(const U &p)
Definition: TrackParametersHelper.h:134
VKalVrtAthena::VrtSecInclusive::WrkVrt::Charge
long int Charge
list of VKalVrt fit chi2 for each track
Definition: VrtSecInclusive.h:350
get_generator_info.result
result
Definition: get_generator_info.py:21
VKalVrtAthena::VrtSecInclusive::WrkVrt::associatedTrackIndices
std::deque< long int > associatedTrackIndices
list if indices in TrackParticleContainer for selectedBaseTracks
Definition: VrtSecInclusive.h:343
maxValue
#define maxValue(current, test)
Definition: CompoundLayerMaterialCreator.h:22
VKalVrtAthena::VrtSecInclusive::WrkVrt::isGood
bool isGood
Definition: VrtSecInclusive.h:341
xAOD::Vertex
Vertex_v1 Vertex
Define the latest version of the vertex class.
Definition: Event/xAOD/xAODTracking/xAODTracking/Vertex.h:16
DataVector::emplace_back
value_type emplace_back(value_type pElem)
Add an element to the end of the collection.
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
xAOD::uint8_t
uint8_t
Definition: Muon_v1.cxx:558
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
SG::Accessor< float >
xAOD::TrackParticle_v1::charge
float charge() const
Returns the charge.
Definition: TrackParticle_v1.cxx:143
index
Definition: index.py:1
drawFromPickle.candidates
candidates
Definition: drawFromPickle.py:271
Trk::indices
std::pair< long int, long int > indices
Definition: AlSymMatBase.h:24
PGraph.h
VKalVrtAthena::AlgConsts::chi2PerTrackInitValue
constexpr double chi2PerTrackInitValue
Definition: Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/Constants.h:22
Trk::MappedVertex::identifierHash
int identifierHash
Definition: IVertexMapper.h:27
EventPrimitivesHelpers.h
VKalVrtAthena::AlgConsts::invalidInt
constexpr int invalidInt
Definition: Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/Constants.h:25
TRTCalib_cfilter.p1
p1
Definition: TRTCalib_cfilter.py:130
ALFA_EventTPCnv_Dict::t1
std::vector< ALFA_RawDataCollection_p1 > t1
Definition: ALFA_EventTPCnvDict.h:43
plotBeamSpotVxVal.cov
cov
Definition: plotBeamSpotVxVal.py:200
ElectronxAODHelpers.h
VKalVrtAthena::VrtSecInclusive::WrkVrt::vertexMom
TLorentzVector vertexMom
VKalVrt fit vertex position.
Definition: VrtSecInclusive.h:345
Trk::MappedVertex
Definition: IVertexMapper.h:23
Trk::z0
@ z0
Definition: ParamDefs.h:64
VKalVrtAthena::VrtSecInclusive::track_summary_properties
Definition: VrtSecInclusive.h:493
xAOD::numberOfPixelHits
@ numberOfPixelHits
these are the pixel hits, including the b-layer [unit8_t].
Definition: TrackingPrimitives.h:260
SG::ConstAccessor< char >
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
dumpTruth.statusCode
statusCode
Definition: dumpTruth.py:85
NtupleVars.h
Phi_mpi_pi
__HOSTDEV__ double Phi_mpi_pi(double)
Definition: GeoRegion.cxx:10
python.TrigEgammaMonitorHelper.TH2F
def TH2F(name, title, nxbins, bins_par2, bins_par3, bins_par4, bins_par5=None, bins_par6=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:45
IDTPM::nPixelHits
float nPixelHits(const U &p)
Definition: TrackParametersHelper.h:354
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition: EventPrimitives.h:50
XMLtoHeader.count
count
Definition: XMLtoHeader.py:84
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
mergePhysValFiles.end
end
Definition: DataQuality/DataQualityUtils/scripts/mergePhysValFiles.py:92
VrtSecInclusive.h
VKalVrtAthena
Definition: AANT_Tools.cxx:25
Trk::MappedVertex::layerIndex
int layerIndex
Definition: IVertexMapper.h:30
TRTCalib_cfilter.p2
p2
Definition: TRTCalib_cfilter.py:131
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:209
SG::Decorator
Helper class to provide type-safe access to aux data.
Definition: Decorator.h:59
fillPileUpNoiseLumi.next
next
Definition: fillPileUpNoiseLumi.py:52
VKalVrtAthena::VrtSecInclusive::WrkVrt::Chi2PerTrk
std::vector< double > Chi2PerTrk
VKalVrt fit chi2 result.
Definition: VrtSecInclusive.h:349
lumiFormat.i
int i
Definition: lumiFormat.py:85
Trk::theta
@ theta
Definition: ParamDefs.h:66
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
xAOD::VxType::SecVtx
@ SecVtx
Secondary vertex.
Definition: TrackingPrimitives.h:573
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
VKalVrtAthena::VrtSecInclusive::WrkVrt::vertexCov
std::vector< double > vertexCov
VKalVrt fit vertex 4-momentum.
Definition: VrtSecInclusive.h:346
VKalVrtAthena::VrtSecInclusive::WrkVrt::Chi2
double Chi2
VKalVrt fit covariance.
Definition: VrtSecInclusive.h:347
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
std::remove_if
std::reverse_iterator< DataModel_detail::iterator< DVL > > remove_if(typename std::reverse_iterator< DataModel_detail::iterator< DVL > > beg, typename std::reverse_iterator< DataModel_detail::iterator< DVL > > end, Predicate pred)
Specialization of remove_if for DataVector/List.
Definition: DVL_algorithms.h:114
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
TrackSummary.h
VKalVrtAthena::vtxVtxDistance
double vtxVtxDistance(const Amg::Vector3D &v1, const Amg::Vector3D &v2)
Definition: Reconstruction/VKalVrt/VrtSecInclusive/src/Utilities.cxx:55
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
VKalVrtAthena::VrtSecInclusive::WrkVrt::closestWrkVrtIndex
unsigned long closestWrkVrtIndex
list of track parameters wrt the reconstructed vertex
Definition: VrtSecInclusive.h:352
python.compareNtuple.nGood
nGood
Definition: compareNtuple.py:55
DeMoUpdate.tmp
string tmp
Definition: DeMoUpdate.py:1167
SG::AuxElement::index
size_t index() const
Return the index of this element within its container.
DataVector
Derived DataVector<T>.
Definition: DataVector.h:795
Reflex::FINAL
@ FINAL
Definition: RootType.h:25
VKalVrtAthena::PhysConsts::mass_electron
constexpr double mass_electron
Definition: Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/Constants.h:15
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
Trk::d0
@ d0
Definition: ParamDefs.h:63
VKalVrtAthena::PhysConsts::mass_chargedPion
constexpr double mass_chargedPion
Definition: Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/Constants.h:14
Accessor.h
Helper class to provide type-safe access to aux data.
checkTriggerxAOD.found
found
Definition: checkTriggerxAOD.py:328
VKalVrtAthena::VrtSecInclusive::WrkVrt::selectedTrackIndices
std::deque< long int > selectedTrackIndices
flagged true for good vertex candidates
Definition: VrtSecInclusive.h:342
VKalVrtAthena::VrtSecInclusive::WrkVrt
Definition: VrtSecInclusive.h:340
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
VKalVrtAthena::VrtSecInclusive::WrkVrt::nTracksTotal
unsigned nTracksTotal() const
Definition: VrtSecInclusive.h:357
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
ReadCellNoiseFromCoolCompare.v2
v2
Definition: ReadCellNoiseFromCoolCompare.py:364
VKalVrtAthena::VrtSecInclusive::WrkVrt::closestWrkVrtValue
double closestWrkVrtValue
stores the index of the closest WrkVrt in std::vector<WrkVrt>
Definition: VrtSecInclusive.h:353
python.PyAthena.v
v
Definition: PyAthena.py:154
ALFA_EventTPCnv_Dict::t2
std::vector< ALFA_RawDataContainer_p1 > t2
Definition: ALFA_EventTPCnvDict.h:44
InDetDD::detail::kEta
@ kEta
Definition: PixelGeoUtils.h:22
Trk::vertex
@ vertex
Definition: MeasurementType.h:21
DeMoScan.index
string index
Definition: DeMoScan.py:362
VKalVrtAthena::AlgConsts::invalidFloat
constexpr double invalidFloat
Definition: Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/Constants.h:24
xAOD::Vertex_v1
Class describing a Vertex.
Definition: Vertex_v1.h:42
VKalVrtAthena::VrtSecInclusive::WrkVrt::fitQuality
double fitQuality() const
Definition: VrtSecInclusive.h:358
std::sort
void sort(typename std::reverse_iterator< DataModel_detail::iterator< DVL > > beg, typename std::reverse_iterator< DataModel_detail::iterator< DVL > > end, const Compare &comp)
Specialization of sort for DataVector/List.
Definition: DVL_algorithms.h:623
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
VKalVrtAthena::VrtSecInclusive::WrkVrt::TrkAtVrt
std::vector< std::vector< double > > TrkAtVrt
total charge of the vertex
Definition: VrtSecInclusive.h:351
ReadCellNoiseFromCoolCompare.s2
s2
Definition: ReadCellNoiseFromCoolCompare.py:379
Trk::qOverP
@ qOverP
perigee
Definition: ParamDefs.h:67
python.CaloCondTools.log
log
Definition: CaloCondTools.py:20
RunTileMonitoring.clusters
clusters
Definition: RunTileMonitoring.py:133
python.changerun.pv
pv
Definition: changerun.py:79
InDetDD::detail::kPhi
@ kPhi
Definition: PixelGeoUtils.h:22
Trk::MappedVertex::valid
bool valid
Definition: IVertexMapper.h:34
SG::ConstAccessor::isAvailable
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
Trk::phi
@ phi
Definition: ParamDefs.h:75
Tools.h
xAOD::TrackParticle_v1
Class describing a TrackParticle.
Definition: TrackParticle_v1.h:44
SG::AuxElement::container
const SG::AuxVectorData * container() const
Return the container holding this element.
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
Amg::distance
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
Definition: GeoPrimitivesHelpers.h:54
python.compressB64.c
def c
Definition: compressB64.py:93
python.AutoConfigFlags.msg
msg
Definition: AutoConfigFlags.py:7
VKalVrtAthena::AlgConsts::invalidUnsigned
constexpr unsigned invalidUnsigned
Definition: Reconstruction/VKalVrt/VrtSecInclusive/VrtSecInclusive/Constants.h:26
Trk::MappedVertex::localPosition
Amg::Vector3D localPosition
Definition: IVertexMapper.h:25
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
xAOD::numberOfInnermostPixelLayerHits
@ numberOfInnermostPixelLayerHits
these are the hits in the 0th pixel barrel layer
Definition: TrackingPrimitives.h:238
Generated on Thu Oct 30 2025 21:21:59 for ATLAS Offline Software by   doxygen 1.8.18