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