ATLAS Offline Software
VertexingAlgs.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 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 ) m_decor_isAssociated = std::make_unique< SG::AuxElement::Decorator<char> >( "is_associated" + m_jp.augVerString );
1071 
1072  ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": #verticess = " << workVerticesContainer->size() );
1073 
1074  unsigned associateCounter { 0 };
1075 
1076  // Loop over vertices
1077  for( auto& wrkvrt : *workVerticesContainer ) {
1078 
1079  if( !wrkvrt.isGood ) continue;
1080  if( wrkvrt.selectedTrackIndices.size() <= 1 ) continue;
1081 
1082  improveVertexChi2( wrkvrt );
1083 
1084  wrkvrt.Chi2_core = wrkvrt.Chi2;
1085 
1086  auto& vertexPos = wrkvrt.vertex;
1087 
1088  std::vector<double> distanceToPVs;
1089 
1090  for( const auto* pv : *pvs ) {
1091  distanceToPVs.emplace_back( VKalVrtAthena::vtxVtxDistance( vertexPos, pv->position() ) );
1092  }
1093  const auto& minDistance = *( std::min_element( distanceToPVs.begin(), distanceToPVs.end() ) );
1094 
1095  if( minDistance < m_jp.associateMinDistanceToPV ) continue;
1096 
1097 
1098  ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": vertex pos = (" << vertexPos.x() << ", " << vertexPos.y() << ", " << vertexPos.z() << "), "
1099  "#selected = " << wrkvrt.selectedTrackIndices.size() << ", #assoc = " << wrkvrt.associatedTrackIndices.size() );
1100 
1101  std::vector<const xAOD::TrackParticle*> candidates;
1102 
1103  // Search for candidate tracks
1104  for( auto itr = allTracks->begin(); itr != allTracks->end(); ++itr ) {
1105  const auto* trk = *itr;
1106 
1107  // If the track is already used for any DV candidate, reject.
1108  {
1109  auto result = std::find_if( workVerticesContainer->begin(), workVerticesContainer->end(),
1110  [&] ( WrkVrt& wrkvrt ) {
1111  auto found = std::find_if( wrkvrt.selectedTrackIndices.begin(), wrkvrt.selectedTrackIndices.end(),
1112  [&]( long int index ) {
1113  // 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
1114  if (m_jp.doSelectTracksFromElectrons || m_jp.doSelectIDAndGSFTracks) {
1115  const xAOD::TrackParticle *id_tr;
1116  id_tr = xAOD::EgammaHelpers::getOriginalTrackParticleFromGSF(m_selectedTracks->at(index));
1117  return trk == m_selectedTracks->at(index) or trk == id_tr;
1118  }
1119  else{
1120  return trk == m_selectedTracks->at(index);
1121  }
1122  } );
1123  return found != wrkvrt.selectedTrackIndices.end();
1124  } );
1125  if( result != workVerticesContainer->end() ) continue;
1126  }
1127 
1128  // If the track is already registered to the associated track list, reject.
1129  {
1130  auto result = std::find_if( m_associatedTracks->begin(), m_associatedTracks->end(),
1131  [&] (const auto* atrk) { return trk == atrk; } );
1132  if( result != m_associatedTracks->end() ) continue;
1133  }
1134 
1135  // Reject PV-associated tracks
1136  // if( !selectTrack_notPVassociated( trk ) ) continue;
1137 
1138  // pT selection
1139  if( trk->pt() < m_jp.associatePtCut ) continue;
1140 
1141  // chi2 selection
1142  if( trk->chiSquared() / trk->numberDoF() > m_jp.associateChi2Cut ) continue;
1143 
1144  // Hit pattern consistentcy requirement
1145  if( !checkTrackHitPatternToVertexOuterOnly( trk, vertexPos ) ) continue;
1146 
1147  // Get the closest approach
1148  std::vector<double> impactParameters;
1149  std::vector<double> impactParErrors;
1150 
1151  if( !getSVImpactParameters( trk, vertexPos, impactParameters, impactParErrors) ) continue;
1152 
1153  if( std::abs( impactParameters.at(0) ) / sqrt( impactParErrors.at(0) ) > m_jp.associateMaxD0Signif ) continue;
1154  if( std::abs( impactParameters.at(1) ) / sqrt( impactParErrors.at(1) ) > m_jp.associateMaxZ0Signif ) continue;
1155 
1156  ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": trk " << trk
1157  << ": d0 to vtx = " << impactParameters.at(k_d0)
1158  << ", z0 to vtx = " << impactParameters.at(k_z0)
1159  << ", distance to vtx = " << hypot( impactParameters.at(k_d0), impactParameters.at(k_z0) ) );
1160 
1161  candidates.emplace_back( trk );
1162 
1163  }
1164 
1165  ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": number of candidate tracks = " << candidates.size() );
1166 
1167  std::unique_ptr<Trk::IVKalState> state = m_fitSvc->makeState();
1168  // Attempt to add the track to the vertex and try fitting
1169  for( const auto* trk : candidates ) {
1170 
1171  ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": attempting to associate track = " << trk );
1172 
1173  // Backup the current vertes status
1174  WrkVrt wrkvrt_backup = wrkvrt;
1175 
1176  m_fitSvc->setApproximateVertex( vertexPos.x(), vertexPos.y(), vertexPos.z(), *state );
1177 
1178  std::vector<const xAOD::TrackParticle*> baseTracks;
1179  std::vector<const xAOD::NeutralParticle*> dummyNeutrals;
1180 
1181  wrkvrt.Chi2PerTrk.clear();
1182 
1183  for( const auto& index : wrkvrt.selectedTrackIndices ) {
1184  baseTracks.emplace_back( m_selectedTracks->at( index ) );
1185  wrkvrt.Chi2PerTrk.emplace_back( AlgConsts::chi2PerTrackInitValue );
1186  }
1187  for( const auto& index : wrkvrt.associatedTrackIndices ) {
1188  baseTracks.emplace_back( m_associatedTracks->at( index ) );
1189  wrkvrt.Chi2PerTrk.emplace_back( AlgConsts::chi2PerTrackInitValue );
1190  }
1191 
1192  baseTracks.emplace_back( trk );
1193  wrkvrt.Chi2PerTrk.emplace_back( AlgConsts::chi2PerTrackInitValue );
1194 
1195  Amg::Vector3D initPos;
1196 
1197  {
1198  StatusCode sc = m_fitSvc->VKalVrtFitFast( baseTracks, initPos, *state );/* Fast crude estimation */
1199 
1200  if( sc.isFailure() ) ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": fast crude estimation failed.");
1201 
1202  const auto& diffPos = initPos - vertexPos;
1203 
1204  if( diffPos.norm() > 10. ) {
1205 
1206  ATH_MSG_VERBOSE( " > " << __FUNCTION__ << ": approx vertex as original" );
1207  m_fitSvc->setApproximateVertex( vertexPos.x(), vertexPos.y(), vertexPos.z(), *state );
1208 
1209  } else {
1210 
1211  ATH_MSG_VERBOSE( " > " << __FUNCTION__ << ": approx vertex set to (" << initPos.x() << ", " << initPos.y() << ", " << initPos.z() << ")" );
1212  m_fitSvc->setApproximateVertex( initPos.x(), initPos.y(), initPos.z(), *state );
1213 
1214  }
1215  }
1216 
1217 
1218  ATH_MSG_VERBOSE( " > " << __FUNCTION__ << ": now vertex fitting..." );
1219 
1220  StatusCode sc = m_fitSvc->VKalVrtFit(baseTracks, dummyNeutrals,
1221  wrkvrt.vertex,
1222  wrkvrt.vertexMom,
1223  wrkvrt.Charge,
1224  wrkvrt.vertexCov,
1225  wrkvrt.Chi2PerTrk,
1226  wrkvrt.TrkAtVrt,
1227  wrkvrt.Chi2,
1228  *state);
1229 
1230  if( sc.isFailure() ) {
1231  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": VKalVrtFit failure. Revert to backup");
1232  wrkvrt = wrkvrt_backup;
1233 
1234  if( m_jp.FillHist ) m_hists["associateMonitor"]->Fill( 1 );
1235 
1236  continue;
1237  }
1238 
1239 
1240  if( m_jp.FillHist ) m_hists["associateMonitor"]->Fill( 0 );
1241 
1242  auto& cov = wrkvrt.vertexCov;
1243 
1244  ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": succeeded in associating. New vertex pos = (" << vertexPos.perp() << ", " << vertexPos.z() << ", " << vertexPos.perp()*vertexPos.phi() << ")" );
1245  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) << ")" );
1246 
1247  associateCounter++;
1248 
1249  wrkvrt.associatedTrackIndices.emplace_back( m_associatedTracks->size() );
1250 
1251  m_associatedTracks->emplace_back( trk );
1252  (*m_decor_isAssociated)( *trk ) = true;
1253 
1254  }
1255 
1256  }
1257 
1258  ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
1259  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": total associated number of tracks = " << associateCounter );
1260  ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
1261 
1262  return StatusCode::SUCCESS;
1263  }
1264 
1265 
1266  //____________________________________________________________________________________________________
1267  StatusCode VrtSecInclusive::mergeByShuffling( std::vector<WrkVrt> *workVerticesContainer )
1268  {
1269 
1270  ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": #verticess = " << workVerticesContainer->size() );
1271 
1272  unsigned mergeCounter { 0 };
1273 
1274  // First, sort WrkVrt by the track multiplicity
1275  std::sort( workVerticesContainer->begin(), workVerticesContainer->end(), [](WrkVrt& v1, WrkVrt& v2) { return v1.selectedTrackIndices.size() < v2.selectedTrackIndices.size(); } );
1276 
1277  // Loop over vertices (small -> large Ntrk order)
1278  for( auto& wrkvrt : *workVerticesContainer ) {
1279  if( !wrkvrt.isGood ) continue;
1280  if( wrkvrt.selectedTrackIndices.size() <= 1 ) continue;
1281 
1282  // Reverse iteration: large Ntrk -> small Ntrk order
1283  for( auto ritr = workVerticesContainer->rbegin(); ritr != workVerticesContainer->rend(); ++ritr ) {
1284  auto& vertexToMerge = *ritr;
1285 
1286  if( !vertexToMerge.isGood ) continue;
1287  if( vertexToMerge.selectedTrackIndices.size() <= 1 ) continue;
1288  if( &wrkvrt == &vertexToMerge ) continue;
1289  if( vertexToMerge.selectedTrackIndices.size() < wrkvrt.selectedTrackIndices.size() ) continue;
1290 
1291  const double& significance = significanceBetweenVertices( wrkvrt, vertexToMerge );
1292 
1293  if( significance > m_jp.mergeByShufflingMaxSignificance ) continue;
1294 
1295  bool mergeFlag { false };
1296 
1297  ATH_MSG_DEBUG(" > " << __FUNCTION__
1298  << ": vertex " << &wrkvrt << " #tracks = " << wrkvrt.selectedTrackIndices.size()
1299  << " --> to Merge : " << &vertexToMerge << ", #tracks = " << vertexToMerge.selectedTrackIndices.size()
1300  << " significance = " << significance );
1301 
1302  double min_signif = AlgConsts::maxValue;
1303 
1304  // Method 1. Assume that the solution is somewhat wrong, and the solution gets correct if it starts from the other vertex position
1305  if( m_jp.doSuggestedRefitOnMerging && !mergeFlag ) {
1306  WrkVrt testVertex = wrkvrt;
1307  StatusCode sc = refitVertexWithSuggestion( testVertex, vertexToMerge.vertex );
1308  if( sc.isFailure() ) {
1309  //ATH_MSG_WARNING(" > " << __FUNCTION__ << ": detected vertex fitting failure!" );
1310  } else {
1311 
1312  const auto signif = significanceBetweenVertices( testVertex, vertexToMerge );
1313  if( signif < min_signif ) min_signif = signif;
1314 
1315  if( signif < m_jp.mergeByShufflingAllowance ) {
1316  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": method1: vertexToMerge " << &vertexToMerge << ": test signif = " << signif );
1317  mergeFlag = true;
1318 
1319  }
1320 
1321  if( m_jp.FillHist && min_signif > 0. ) m_hists["shuffleMinSignif1"]->Fill( log10( min_signif ) );
1322  if( m_jp.FillHist && mergeFlag ) { m_hists["mergeType"]->Fill( SHUFFLE1 ); }
1323  }
1324  }
1325 
1326  // Method 2. magnet merging: borrowing another track from the target vertex to merge
1327  if( m_jp.doMagnetMerging && !mergeFlag ) {
1328 
1329  // Loop over tracks in vertexToMerge
1330  for( auto& index : vertexToMerge.selectedTrackIndices ) {
1331 
1332  WrkVrt testVertex = wrkvrt;
1333  testVertex.selectedTrackIndices.emplace_back( index );
1334 
1335  StatusCode sc = refitVertexWithSuggestion( testVertex, vertexToMerge.vertex );
1336  if( sc.isFailure() ) {
1337  //ATH_MSG_WARNING(" > " << __FUNCTION__ << ": detected vertex fitting failure!" );
1338  } else {
1339 
1340  const auto signif = significanceBetweenVertices( testVertex, vertexToMerge );
1341  if( signif < min_signif ) min_signif = signif;
1342 
1343  if( signif < m_jp.mergeByShufflingAllowance ) {
1344  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": method2: vertexToMerge " << &vertexToMerge << " track index " << index << ": test signif = " << signif );
1345  mergeFlag = true;
1346  }
1347 
1348  }
1349  }
1350 
1351  if( m_jp.FillHist && min_signif > 0. ) m_hists["shuffleMinSignif2"]->Fill( log10( min_signif ) );
1352 
1353  if( m_jp.FillHist && mergeFlag ) { m_hists["mergeType"]->Fill( SHUFFLE2 ); }
1354  }
1355 
1356  // Method 3. Attempt to force merge
1357  if( m_jp.doWildMerging && !mergeFlag ) {
1358 
1359  WrkVrt testVertex = wrkvrt;
1360 
1361  for( auto& index : vertexToMerge.selectedTrackIndices ) {
1362  testVertex.selectedTrackIndices.emplace_back( index );
1363  }
1364 
1365  StatusCode sc = refitVertexWithSuggestion( testVertex, vertexToMerge.vertex );
1366  if( sc.isFailure() ) {
1367  //ATH_MSG_WARNING(" > " << __FUNCTION__ << ": detected vertex fitting failure!" );
1368  } else {
1369 
1370  const auto signif = significanceBetweenVertices( testVertex, vertexToMerge );
1371  if( signif < min_signif ) min_signif = signif;
1372 
1373  if( signif < m_jp.mergeByShufflingAllowance ) {
1374  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": method3: vertexToMerge " << &vertexToMerge << ": test signif = " << signif );
1375  mergeFlag = true;
1376  }
1377 
1378  if( m_jp.FillHist && min_signif > 0. ) m_hists["shuffleMinSignif3"]->Fill( log10( min_signif ) );
1379  if( m_jp.FillHist && mergeFlag ) { m_hists["mergeType"]->Fill( SHUFFLE3 ); }
1380 
1381  }
1382  }
1383 
1384 
1385  if( mergeFlag ) {
1386  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": vertexToMerge " << &vertexToMerge << " ==> min signif = " << min_signif << " judged to merge" );
1387 
1388  auto vertexToMerge_backup = vertexToMerge;
1389  auto wrkvrt_backup = wrkvrt;
1390 
1391  StatusCode sc = mergeVertices( vertexToMerge, wrkvrt );
1392  if( sc.isFailure() ) {
1393  vertexToMerge = vertexToMerge_backup;
1394  wrkvrt = wrkvrt_backup;
1395  continue;
1396  }
1397 
1398  improveVertexChi2( wrkvrt );
1399 
1400  mergeCounter++;
1401  }
1402 
1403  }
1404 
1405  }
1406 
1407  ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
1408  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Number of merges = " << mergeCounter );
1409  ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
1410 
1411  return StatusCode::SUCCESS;
1412  }
1413 
1414 
1415  //____________________________________________________________________________________________________
1416  StatusCode VrtSecInclusive::mergeFinalVertices( std::vector<WrkVrt> *workVerticesContainer )
1417  {
1418 
1419  unsigned mergeCounter { 0 };
1420 
1421  while (true) {
1422  //
1423  // Minimal vertex-vertex distance
1424  //
1425  for( auto& wrkvrt : *workVerticesContainer) {
1426  wrkvrt.closestWrkVrtIndex = AlgConsts::invalidUnsigned;
1427  wrkvrt.closestWrkVrtValue = AlgConsts::maxValue;
1428  }
1429 
1430  std::pair<unsigned, unsigned> indexPair { AlgConsts::invalidUnsigned, AlgConsts::invalidUnsigned };
1431  auto minDistance = findMinVerticesPair( workVerticesContainer, indexPair, &VrtSecInclusive::distanceBetweenVertices );
1432 
1433  if( minDistance == AlgConsts::maxValue ) break;
1434  if( indexPair.first == AlgConsts::invalidUnsigned ) break;
1435  if( indexPair.second == AlgConsts::invalidUnsigned ) break;
1436 
1437  auto& v1 = workVerticesContainer->at(indexPair.first);
1438  auto& v2 = workVerticesContainer->at(indexPair.second);
1439 
1440  const double averageRadius = ( v1.vertex.perp() + v2.vertex.perp() ) / 2.0;
1441 
1442  if( minDistance > m_jp.VertexMergeFinalDistCut + m_jp.VertexMergeFinalDistScaling * averageRadius ) {
1443  ATH_MSG_DEBUG( "Vertices " << indexPair.first << " and " << indexPair.second
1444  <<" are separated by distance " << minDistance );
1445  break;
1446  }
1447 
1448  ATH_MSG_DEBUG( "Merging FINAL vertices " << indexPair.first << " and " << indexPair.second
1449  <<" which are separated by distance "<< minDistance );
1450 
1451  StatusCode sc = mergeVertices( v1, v2 );
1452  if( sc.isFailure() ) {}
1453  if( m_jp.FillHist ) { m_hists["mergeType"]->Fill( FINAL ); }
1454 
1455  improveVertexChi2( v1 );
1456 
1457  mergeCounter++;
1458 
1459  }
1460 
1461  ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
1462  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Number of merges = " << mergeCounter );
1463  ATH_MSG_DEBUG(" > " << __FUNCTION__ << "----------------------------------------------" );
1464 
1465  return StatusCode::SUCCESS;
1466 
1467  } // end of mergeFinalVertices
1468 
1469 
1470 
1471  //____________________________________________________________________________________________________
1472  StatusCode VrtSecInclusive::refitAndSelectGoodQualityVertices( std::vector<WrkVrt> *workVerticesContainer )
1473  {
1474 
1475  // Output SVs as xAOD::Vertex
1476  // Needs a conversion function from workVerticesContainer to xAOD::Vertex here.
1477  // The supposed form of the function will be as follows:
1478 
1479  try {
1480 
1481  xAOD::VertexContainer *secondaryVertexContainer( nullptr );
1482  ATH_CHECK( evtStore()->retrieve( secondaryVertexContainer, "VrtSecInclusive_" + m_jp.secondaryVerticesContainerName + m_jp.augVerString ) );
1483 
1484  const xAOD::TrackParticleContainer* trackParticleContainer ( nullptr );
1485  ATH_CHECK( evtStore()->retrieve( trackParticleContainer, m_jp.TrackLocation) );
1486 
1487  enum { kPt, kEta, kPhi, kD0, kZ0, kErrP, kErrD0, kErrZ0, kChi2SV };
1488  if( m_trkDecors.empty() ) {
1489  m_trkDecors.emplace( kPt, SG::AuxElement::Decorator<float>("pt_wrtSV" + m_jp.augVerString) );
1490  m_trkDecors.emplace( kEta, SG::AuxElement::Decorator<float>("eta_wrtSV" + m_jp.augVerString) );
1491  m_trkDecors.emplace( kPhi, SG::AuxElement::Decorator<float>("phi_wrtSV" + m_jp.augVerString) );
1492  m_trkDecors.emplace( kD0, SG::AuxElement::Decorator<float>("d0_wrtSV" + m_jp.augVerString) );
1493  m_trkDecors.emplace( kZ0, SG::AuxElement::Decorator<float>("z0_wrtSV" + m_jp.augVerString) );
1494  m_trkDecors.emplace( kErrP, SG::AuxElement::Decorator<float>("errP_wrtSV" + m_jp.augVerString) );
1495  m_trkDecors.emplace( kErrD0, SG::AuxElement::Decorator<float>("errd0_wrtSV" + m_jp.augVerString) );
1496  m_trkDecors.emplace( kErrZ0, SG::AuxElement::Decorator<float>("errz0_wrtSV" + m_jp.augVerString) );
1497  m_trkDecors.emplace( kChi2SV, SG::AuxElement::Decorator<float>("chi2_toSV" + m_jp.augVerString) );
1498  }
1499  if( !m_decor_is_svtrk_final ) m_decor_is_svtrk_final = std::make_unique< SG::AuxElement::Decorator<char> >( "is_svtrk_final" + m_jp.augVerString );
1500 
1501  std::map<const WrkVrt*, const xAOD::Vertex*> wrkvrtLinkMap;
1502 
1503  //----------------------------------------------------------
1504  const auto& ctx = Gaudi::Hive::currentContext();
1505 
1506  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": input #vertices = " << workVerticesContainer->size() );
1507 
1508  // Loop over vertices
1509  for( auto& wrkvrt : *workVerticesContainer ) {
1510 
1511  ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": candidate vertex: "
1512  << " isGood = " << (wrkvrt.isGood? "true" : "false")
1513  << ", #ntrks = " << wrkvrt.nTracksTotal()
1514  << ", #selectedTracks = " << wrkvrt.selectedTrackIndices.size()
1515  << ", #associatedTracks = " << wrkvrt.associatedTrackIndices.size()
1516  << ", chi2/ndof = " << wrkvrt.Chi2 / ( wrkvrt.ndof() + AlgConsts::infinitesimal )
1517  << ", (r, z) = (" << wrkvrt.vertex.perp()
1518  <<", " << wrkvrt.vertex.z() << ")" );
1519 
1520  if( m_jp.FillHist ) m_hists["finalCutMonitor"]->Fill( 0 );
1521 
1522  if( m_jp.removeFakeVrt && m_jp.removeFakeVrtLate ) {
1523  removeInconsistentTracks( wrkvrt );
1524  }
1525 
1526  if( wrkvrt.nTracksTotal() < 2 ) {
1527  ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": ntrk < 2 --> rejected." );
1528  continue; /* Bad vertices */
1529  }
1530 
1531  if( m_jp.FillHist ) m_hists["finalCutMonitor"]->Fill( 1 );
1532 
1533 
1534  // Remove track if the vertex is inner than IBL and the track does not have pixel hits!
1535  if( wrkvrt.vertex.perp() < 31.0 ) {
1536 
1537  // for selected tracks
1538  wrkvrt.selectedTrackIndices.erase( std::remove_if( wrkvrt.selectedTrackIndices.begin(), wrkvrt.selectedTrackIndices.end(),
1539  [&]( auto& index ) {
1540  auto* trk = m_selectedTracks->at( index );
1541  uint8_t nPixelHits { 0 }; trk->summaryValue( nPixelHits, xAOD::numberOfPixelHits );
1542  return ( nPixelHits < 3 );
1543  } ),
1544  wrkvrt.selectedTrackIndices.end() );
1545 
1546  // for associated tracks
1547  wrkvrt.associatedTrackIndices.erase( std::remove_if( wrkvrt.associatedTrackIndices.begin(), wrkvrt.associatedTrackIndices.end(),
1548  [&]( auto& index ) {
1549  auto* trk = m_associatedTracks->at( index );
1550  uint8_t nPixelHits { 0 }; trk->summaryValue( nPixelHits, xAOD::numberOfPixelHits );
1551  return ( nPixelHits < 3 );
1552  } ),
1553  wrkvrt.associatedTrackIndices.end() );
1554 
1555  auto statusCode = refitVertex( wrkvrt );
1556  if( statusCode.isFailure() ) {}
1557 
1558  }
1559 
1560 
1561  if( m_jp.doFinalImproveChi2 ) {
1562 
1563  WrkVrt backup = wrkvrt;
1564 
1565  improveVertexChi2( wrkvrt );
1566 
1567  if( wrkvrt.fitQuality() > backup.fitQuality() ) wrkvrt = backup;
1568 
1569  }
1570 
1571  // If the number of remaining tracks is less than 2, drop.
1572  if( wrkvrt.nTracksTotal() < 2 ) continue;
1573 
1574  // Select only vertices with keeping more than 2 selectedTracks
1575  if( wrkvrt.selectedTrackIndices.size() < 2 ) continue;
1576 
1577 
1578  if( m_jp.FillHist ) m_hists["finalCutMonitor"]->Fill( 2 );
1579 
1580 
1581  {
1582  WrkVrt backup = wrkvrt;
1583 
1584  StatusCode sc = refitVertex( wrkvrt );
1585  if( sc.isFailure() ) {
1586 
1587  auto indices = wrkvrt.associatedTrackIndices;
1588 
1589  wrkvrt.associatedTrackIndices.clear();
1590  sc = refitVertex( wrkvrt );
1591  if( sc.isFailure() ) {
1592  ATH_MSG_WARNING(" > " << __FUNCTION__ << ": detected vertex fitting failure!" );
1593  wrkvrt = backup;
1594  }
1595  if( wrkvrt.fitQuality() > backup.fitQuality() ) wrkvrt = backup;
1596 
1597  for( auto& index : indices ) {
1598  backup = wrkvrt;
1599  wrkvrt.associatedTrackIndices.emplace_back( index );
1600  sc = refitVertex( wrkvrt );
1601  if( sc.isFailure() || TMath::Prob( wrkvrt.Chi2, wrkvrt.ndof() ) < m_jp.improveChi2ProbThreshold ) {
1602  ATH_MSG_WARNING(" > " << __FUNCTION__ << ": detected vertex fitting failure!" );
1603  wrkvrt = backup;
1604  continue;
1605  }
1606  }
1607 
1608  } else {
1609  if( wrkvrt.fitQuality() > backup.fitQuality() ) wrkvrt = backup;
1610  }
1611  }
1612 
1613  if( m_jp.FillHist ) m_hists["finalCutMonitor"]->Fill( 3 );
1614 
1615  //
1616  // Store good vertices into StoreGate
1617  //
1618  if( m_jp.FillNtuple ) m_ntupleVars->get<unsigned int>( "NumSecVrt" )++;
1619 
1620  TLorentzVector sumP4_pion;
1621  TLorentzVector sumP4_electron;
1622  TLorentzVector sumP4_proton;
1623 
1624  // Pre-check before storing vertex if the SV perigee is available
1625  bool good_flag = true;
1626 
1627  std::map<const std::deque<long int>*, const std::vector<const xAOD::TrackParticle*>&> indicesSet
1628  = {
1629  { &(wrkvrt.selectedTrackIndices), *m_selectedTracks },
1630  { &(wrkvrt.associatedTrackIndices), *m_associatedTracks }
1631  };
1632 
1633  for( auto& pair : indicesSet ) {
1634 
1635  const auto* indices = pair.first;
1636  const auto& tracks = pair.second;
1637 
1638  for( const auto& itrk : *indices ) {
1639  const auto* trk = tracks.at( itrk );
1640  auto sv_perigee = m_trackToVertexTool->perigeeAtVertex(ctx, *trk, wrkvrt.vertex );
1641  if( !sv_perigee ) {
1642  ATH_MSG_INFO(" > " << __FUNCTION__ << ": > Track index " << trk->index() << ": Failed in obtaining the SV perigee!" );
1643  good_flag = false;
1644  }
1645  }
1646 
1647  }
1648 
1649  if( !good_flag ) {
1650  ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": sv perigee could not be obtained --> rejected" );
1651  continue;
1652  }
1653 
1654  if( m_jp.FillHist ) m_hists["finalCutMonitor"]->Fill( 4 );
1655 
1656 
1657  std::vector<const xAOD::TrackParticle*> tracks;
1658  std::vector< std::pair<const xAOD::TrackParticle*, double> > trackChi2Pairs;
1659 
1660  {
1661 
1662  for( auto& pair : indicesSet ) {
1663  for( const auto& index : *pair.first ) tracks.emplace_back( pair.second.at( index ) );
1664  }
1665 
1666  auto trkitr = tracks.begin();
1667  auto chi2itr = wrkvrt.Chi2PerTrk.begin();
1668 
1669  for( ; ( trkitr!=tracks.end() && chi2itr!=wrkvrt.Chi2PerTrk.end() ); ++trkitr, ++chi2itr ) {
1670  trackChi2Pairs.emplace_back( *trkitr, *chi2itr );
1671  }
1672 
1673  }
1674 
1675 
1676  TLorentzVector sumP4_selected;
1677 
1678  bool badIPflag { false };
1679 
1680  // loop over vertex tracks
1681  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Track loop: size = " << tracks.size() );
1682  for( auto& pair : trackChi2Pairs ) {
1683 
1684  const auto* trk = pair.first;
1685  const auto& chi2AtSV = pair.second;
1686 
1687  ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": > Track index " << trk->index() << ": start." );
1688 
1689  track_summary trk_summary;
1690  fillTrackSummary( trk_summary, trk );
1691 
1692  //
1693  // calculate mass/pT of tracks and track parameters
1694  //
1695 
1696  double trk_pt = trk->pt();
1697  double trk_eta = trk->eta();
1698  double trk_phi = trk->phi();
1699 
1700  ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": > Track index " << trk->index() << ": in vrt chg/pt/phi/eta = "
1701  << trk->charge() <<","
1702  <<trk_pt<<","
1703  <<trk_phi<<","
1704  <<trk_eta);
1705 
1707  // Get the perigee of the track at the vertex
1708  ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": > Track index " << trk->index() << ": Get the prigee of the track at the vertex." );
1709 
1710  auto sv_perigee = m_trackToVertexTool->perigeeAtVertex(ctx, *trk, wrkvrt.vertex );
1711  if( !sv_perigee ) {
1712  ATH_MSG_WARNING(" > " << __FUNCTION__ << ": > Track index " << trk->index() << ": Failed in obtaining the SV perigee!" );
1713 
1714  for( auto& pair : m_trkDecors ) {
1715  pair.second( *trk ) = AlgConsts::invalidFloat;
1716  }
1717  (*m_decor_is_svtrk_final)( *trk ) = true;
1718  continue;
1719  }
1720 
1721  double qOverP_wrtSV = sv_perigee->parameters() [Trk::qOverP];
1722  double theta_wrtSV = sv_perigee->parameters() [Trk::theta];
1723  double p_wrtSV = 1.0 / std::abs( qOverP_wrtSV );
1724  double pt_wrtSV = p_wrtSV * sin( theta_wrtSV );
1725  double eta_wrtSV = -log( tan( theta_wrtSV/2. ) );
1726  double phi_wrtSV = sv_perigee->parameters() [Trk::phi];
1727  double d0_wrtSV = sv_perigee->parameters() [Trk::d0];
1728  double z0_wrtSV = sv_perigee->parameters() [Trk::z0];
1729  double errd0_wrtSV = (*sv_perigee->covariance())( Trk::d0, Trk::d0 );
1730  double errz0_wrtSV = (*sv_perigee->covariance())( Trk::z0, Trk::z0 );
1731  double errP_wrtSV = (*sv_perigee->covariance())( Trk::qOverP, Trk::qOverP );
1732 
1733  // xAOD::Track augmentation
1734  ( m_trkDecors.at(kPt) )( *trk ) = pt_wrtSV;
1735  ( m_trkDecors.at(kEta) )( *trk ) = eta_wrtSV;
1736  ( m_trkDecors.at(kPhi) )( *trk ) = phi_wrtSV;
1737  ( m_trkDecors.at(kD0) )( *trk ) = d0_wrtSV;
1738  ( m_trkDecors.at(kZ0) )( *trk ) = z0_wrtSV;
1739  ( m_trkDecors.at(kErrP) )( *trk ) = errP_wrtSV;
1740  ( m_trkDecors.at(kErrD0) )( *trk ) = errd0_wrtSV;
1741  ( m_trkDecors.at(kErrZ0) )( *trk ) = errz0_wrtSV;
1742  ( m_trkDecors.at(kChi2SV))( *trk ) = chi2AtSV;
1743 
1744  (*m_decor_is_svtrk_final)( *trk ) = true;
1745 
1746  TLorentzVector p4wrtSV_pion;
1747  TLorentzVector p4wrtSV_electron;
1748  TLorentzVector p4wrtSV_proton;
1749 
1750  p4wrtSV_pion .SetPtEtaPhiM( pt_wrtSV, eta_wrtSV, phi_wrtSV, PhysConsts::mass_chargedPion );
1751  p4wrtSV_electron.SetPtEtaPhiM( pt_wrtSV, eta_wrtSV, phi_wrtSV, PhysConsts::mass_electron );
1752 
1753  // for selected tracks only
1754  static const SG::ConstAccessor<char> is_associatedAcc("is_associated" + m_jp.augVerString);
1755  if( is_associatedAcc.isAvailable(*trk) ) {
1756  if( !is_associatedAcc(*trk) ) {
1757  sumP4_selected += p4wrtSV_pion;
1758  }
1759  } else {
1760  sumP4_selected += p4wrtSV_pion;
1761  }
1762 
1763  sumP4_pion += p4wrtSV_pion;
1764  sumP4_electron += p4wrtSV_electron;
1765  sumP4_proton += p4wrtSV_proton;
1766 
1767  ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": > Track index " << trk->index() << ": end." );
1768  } // loop over tracks in vertex
1769 
1770  ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": Track loop end. ");
1771 
1772  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Final Sec.Vertex=" << wrkvrt.nTracksTotal() <<", "
1773  <<wrkvrt.vertex.perp() <<", "<<wrkvrt.vertex.z() <<", "
1774  <<wrkvrt.vertex.phi() <<", mass = "<< sumP4_pion.M() << "," << sumP4_electron.M() );
1775 
1776 
1777  //
1778  // calculate opening angle between all 2-track pairs, and store the minimum
1779  //
1780  double minOpAng = AlgConsts::invalidFloat;
1781  std::vector<double> opAngles;
1782 
1783  for( auto itr1 = tracks.begin(); itr1 != tracks.end(); ++itr1 ) {
1784  for( auto itr2 = std::next( itr1 ); itr2 != tracks.end(); ++itr2 ) {
1785  const auto& p1 = (*itr1)->p4().Vect();
1786  const auto& p2 = (*itr2)->p4().Vect();
1787  auto cos = p1 * p2 / p1.Mag() / p2.Mag();
1788  opAngles.emplace_back( cos );
1789  }
1790  }
1791  minOpAng = *( std::max_element( opAngles.begin(), opAngles.end() ) );
1792  if( m_jp.FillNtuple ) m_ntupleVars->get< vector<double> >( "SecVtx_MinOpAng" ).emplace_back(minOpAng);
1793 
1794 
1795  if( m_jp.FillHist ) m_hists["finalCutMonitor"]->Fill( 5 );
1796 
1797  if( badIPflag ) {
1798  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Bad impact parameter signif wrt SV was flagged." );
1799  }
1800 
1801  if (m_jp.doRemoveNonLeptonVertices) {
1802 
1803  bool oneLepMatchTrack = false;
1804  for (const auto *trk: tracks) {
1805  if ( std::find(m_leptonicTracks->begin(), m_leptonicTracks->end(), trk) != m_leptonicTracks->end() ) {
1806  oneLepMatchTrack = true;
1807  break;
1808  }
1809  }
1810 
1811  // If there are no tracks matched to leptons, do not save the container to the output.
1812  if (!oneLepMatchTrack) continue;
1813  }
1814 
1816  // Data filling to xAOD container
1817 
1818  wrkvrt.isGood = true;
1819 
1820  // Firstly store the new vertex to the container before filling properties.
1821  // (This is the feature of xAOD.)
1822  xAOD::Vertex* vertex = new xAOD::Vertex;
1823  secondaryVertexContainer->emplace_back( vertex );
1824 
1825  // Registering the vertex position to xAOD::Vertex
1826  vertex->setPosition( wrkvrt.vertex );
1827 
1828  // Registering the vertex type: SV
1829  vertex->setVertexType( xAOD::VxType::SecVtx );
1830 
1831  // Registering the vertex chi2 and Ndof
1832  // Here, we register the core chi2 of the core (before track association)
1833  vertex->setFitQuality( wrkvrt.Chi2_core, wrkvrt.ndof_core() );
1834 
1835  // Registering the vertex covariance matrix
1836  std::vector<float> fCov(wrkvrt.vertexCov.cbegin(), wrkvrt.vertexCov.cend());
1837  vertex->setCovariance(fCov);
1838 
1839  // Registering the vertex momentum and charge
1840  static const SG::Accessor<float> vtx_pxAcc("vtx_px");
1841  static const SG::Accessor<float> vtx_pyAcc("vtx_py");
1842  static const SG::Accessor<float> vtx_pzAcc("vtx_pz");
1843  static const SG::Accessor<float> vtx_massAcc("vtx_mass");
1844  static const SG::Accessor<float> vtx_chargeAcc("vtx_charge");
1845  static const SG::Accessor<float> chi2_coreAcc("chi2_core");
1846  static const SG::Accessor<float> ndof_coreAcc("ndof_core");
1847  static const SG::Accessor<float> chi2_assocAcc("chi2_assoc");
1848  static const SG::Accessor<float> ndof_assocAcc("ndof_assoc");
1849  static const SG::Accessor<float> massAcc("mass");
1850  static const SG::Accessor<float> mass_eAcc("mass_e");
1851  static const SG::Accessor<float> mass_selectedTracksAcc("mass_selectedTracks");
1852  static const SG::Accessor<float> minOpAngAcc("minOpAng");
1853  static const SG::Accessor<int> num_trksAcc("num_trks");
1854  static const SG::Accessor<int> num_selectedTracksAcc("num_selectedTracks");
1855  static const SG::Accessor<int> num_associatedTracksAcc("num_associatedTracks");
1856  static const SG::Accessor<float> dCloseVrtAcc("dCloseVrt");
1857 
1858  vtx_pxAcc(*vertex) = wrkvrt.vertexMom.Px();
1859  vtx_pyAcc(*vertex) = wrkvrt.vertexMom.Py();
1860  vtx_pzAcc(*vertex) = wrkvrt.vertexMom.Pz();
1861 
1862  vtx_massAcc(*vertex) = wrkvrt.vertexMom.M();
1863  vtx_chargeAcc(*vertex) = wrkvrt.Charge;
1864 
1865  chi2_coreAcc(*vertex) = wrkvrt.Chi2_core;
1866  ndof_coreAcc(*vertex) = wrkvrt.ndof_core();
1867  chi2_assocAcc(*vertex) = wrkvrt.Chi2;
1868  ndof_assocAcc(*vertex) = wrkvrt.ndof();
1869  // Other SV properties
1870  massAcc(*vertex) = sumP4_pion.M();
1871  mass_eAcc(*vertex) = sumP4_electron.M();
1872  mass_selectedTracksAcc(*vertex) = sumP4_selected.M();
1873  minOpAngAcc(*vertex) = minOpAng;
1874  num_trksAcc(*vertex) = wrkvrt.nTracksTotal();
1875  num_selectedTracksAcc(*vertex) = wrkvrt.selectedTrackIndices.size();
1876  num_associatedTracksAcc(*vertex) = wrkvrt.associatedTrackIndices.size();
1877  dCloseVrtAcc(*vertex) = wrkvrt.closestWrkVrtValue;
1878 
1879  // Registering tracks comprising the vertex to xAOD::Vertex
1880  // loop over the tracks comprising the vertex
1881  for( auto trk_id : wrkvrt.selectedTrackIndices ) {
1882 
1883  const xAOD::TrackParticle *trk = m_selectedTracks->at( trk_id );
1884 
1885  // Acquire link the track to the vertex
1886  ElementLink<xAOD::TrackParticleContainer> link_trk( *( dynamic_cast<const xAOD::TrackParticleContainer*>( trk->container() ) ), static_cast<long unsigned int>(trk->index()) );
1887 
1888  // Register the link to the vertex
1889  vertex->addTrackAtVertex( link_trk, 1. );
1890 
1891  }
1892 
1893  for( auto trk_id : wrkvrt.associatedTrackIndices ) {
1894 
1895  const xAOD::TrackParticle *trk = m_associatedTracks->at( trk_id );
1896 
1897  // Acquire link the track to the vertex
1898  ElementLink<xAOD::TrackParticleContainer> link_trk( *( dynamic_cast<const xAOD::TrackParticleContainer*>( trk->container() ) ), static_cast<long unsigned int>(trk->index()) );
1899 
1900  // Register the link to the vertex
1901  vertex->addTrackAtVertex( link_trk, 1. );
1902 
1903  }
1904 
1905 
1906  if( m_jp.doMapToLocal ) {
1907  // Obtain the local mapping of the reconstructed vertex
1908  Trk::MappedVertex mappedVtx = m_vertexMapper->mapToLocal( wrkvrt.vertex );
1909  static const SG::Accessor<int> local_identifierHashAcc("local_identifierHash");
1910  static const SG::Accessor<int> local_layerIndexAcc("local_layerIndex");
1911  static const SG::Accessor<float> local_posXAcc("local_posX");
1912  static const SG::Accessor<float> local_posYAcc("local_posY");
1913  static const SG::Accessor<float> local_posZAcc("local_posZ");
1914  if( mappedVtx.valid ) {
1915  local_identifierHashAcc(*vertex) = mappedVtx.identifierHash;
1916  local_layerIndexAcc(*vertex) = mappedVtx.layerIndex;
1917  local_posXAcc(*vertex) = mappedVtx.localPosition.x();
1918  local_posYAcc(*vertex) = mappedVtx.localPosition.y();
1919  local_posZAcc(*vertex) = mappedVtx.localPosition.z();
1920  } else {
1921  local_identifierHashAcc(*vertex) = AlgConsts::invalidInt;
1922  local_layerIndexAcc(*vertex) = AlgConsts::invalidInt;
1923  local_posXAcc(*vertex) = AlgConsts::invalidFloat;
1924  local_posYAcc(*vertex) = AlgConsts::invalidFloat;
1925  local_posZAcc(*vertex) = AlgConsts::invalidFloat;
1926  }
1927  }
1928 
1929 
1930  // For MC, try to trace down to the truth particles,
1931  // and depending on the topology, categorize the label of the reconstructed vertex.
1932  if( m_jp.doTruth ) {
1933  ATH_CHECK( categorizeVertexTruthTopology( vertex ) );
1934  }
1935 
1936  // Keep the link between wrkvrt and vertex for later use
1937  wrkvrtLinkMap[&wrkvrt] = vertex;
1938 
1939 
1940  } // loop over vertices
1941 
1942  if( m_jp.FillNtuple ) {
1943  ATH_CHECK( fillAANT_SecondaryVertices( secondaryVertexContainer ) );
1944  }
1945 
1946 
1947  // Post process -- Additional augmentations
1948  if( m_jp.doAugmentDVimpactParametersToMuons ) { ATH_CHECK( augmentDVimpactParametersToLeptons<xAOD::Muon> ( "Muons" ) ); }
1949  if( m_jp.doAugmentDVimpactParametersToElectrons ) { ATH_CHECK( augmentDVimpactParametersToLeptons<xAOD::Electron>( "Electrons" ) ); }
1950 
1951  } catch (const std::out_of_range& e) {
1952 
1953  ATH_MSG_WARNING( " > " << __FUNCTION__ << ": out of range error is detected: " << e.what() );
1954 
1955  return StatusCode::SUCCESS;
1956 
1957  } catch( ... ) {
1958 
1959  ATH_MSG_WARNING( " > " << __FUNCTION__ << ": some other error is detected." );
1960 
1961  return StatusCode::SUCCESS;
1962 
1963  }
1964 
1965  return StatusCode::SUCCESS;
1966  }
1967 
1968 
1969  //____________________________________________________________________________________________________
1970  StatusCode VrtSecInclusive::monitorVertexingAlgorithmStep( std::vector<WrkVrt>* workVerticesContainer, const std::string& name, bool final ) {
1971 
1972  if( m_jp.FillIntermediateVertices ) {
1973 
1974  const xAOD::TrackParticleContainer* trackParticleContainer ( nullptr );
1975  ATH_CHECK( evtStore()->retrieve( trackParticleContainer, m_jp.TrackLocation) );
1976 
1977  xAOD::VertexContainer* intermediateVertexContainer { nullptr };
1978 
1979  ATH_CHECK( evtStore()->retrieve( intermediateVertexContainer, "VrtSecInclusive_IntermediateVertices_" + name + m_jp.augVerString ) );
1980 
1981  for( auto& wrkvrt : *workVerticesContainer ) {
1982 
1983  xAOD::Vertex* vertex = new xAOD::Vertex;
1984  intermediateVertexContainer->emplace_back( vertex );
1985 
1986  // Registering the vertex position to xAOD::Vertex
1987  vertex->setPosition( wrkvrt.vertex );
1988 
1989  // Registering the vertex type: SV
1990  vertex->setVertexType( xAOD::VxType::SecVtx );
1991 
1992  // Registering the vertex chi2 and Ndof
1993  int ndof = wrkvrt.ndof();
1994  vertex->setFitQuality( wrkvrt.Chi2, ndof );
1995 
1996  // Registering the vertex covariance matrix
1997  std::vector<float> fCov(wrkvrt.vertexCov.cbegin(), wrkvrt.vertexCov.cend());
1998  vertex->setCovariance(fCov);
1999 
2000  // Registering tracks comprising the vertex to xAOD::Vertex
2001  // loop over the tracks comprising the vertex
2002  for( auto trk_id : wrkvrt.selectedTrackIndices ) {
2003 
2004  const xAOD::TrackParticle *trk = m_selectedTracks->at( trk_id );
2005 
2006  // Acquire link the track to the vertex
2007  ElementLink<xAOD::TrackParticleContainer> link_trk( *( dynamic_cast<const xAOD::TrackParticleContainer*>( trk->container() ) ), static_cast<long unsigned int>(trk->index()) );
2008 
2009  // Register the link to the vertex
2010  vertex->addTrackAtVertex( link_trk, 1. );
2011 
2012  }
2013 
2014  for( auto trk_id : wrkvrt.associatedTrackIndices ) {
2015 
2016  const xAOD::TrackParticle *trk = m_associatedTracks->at( trk_id );
2017 
2018  // Acquire link the track to the vertex
2019  ElementLink<xAOD::TrackParticleContainer> link_trk( *( dynamic_cast<const xAOD::TrackParticleContainer*>( trk->container() ) ), static_cast<long unsigned int>(trk->index()) );
2020 
2021  // Register the link to the vertex
2022  vertex->addTrackAtVertex( link_trk, 1. );
2023 
2024  }
2025  }
2026 
2027  }
2028 
2029 
2030 
2031  if( !m_jp.FillHist ) return StatusCode::SUCCESS;
2032 
2033  printWrkSet( workVerticesContainer, Form("%s (step %u)", name.c_str(), m_vertexingAlgorithmStep) );
2034 
2035  unsigned count = std::count_if( workVerticesContainer->begin(), workVerticesContainer->end(),
2036  []( WrkVrt& v ) { return ( v.selectedTrackIndices.size() + v.associatedTrackIndices.size() ) >= 2; } );
2037 
2038  if( m_vertexingAlgorithmStep == 0 ) {
2039 
2040  const auto compSize = m_selectedTracks->size()*(m_selectedTracks->size() - 1)/2 - m_incomp.size();
2041  m_hists["vertexYield"]->Fill( m_vertexingAlgorithmStep, compSize );
2042 
2043  } else {
2044 
2045  m_hists["vertexYield"]->Fill( m_vertexingAlgorithmStep, count );
2046 
2047  }
2048 
2049  m_hists["vertexYield"]->GetXaxis()->SetBinLabel( m_vertexingAlgorithmStep+1, name.c_str() );
2050 
2051  for( auto& vertex : *workVerticesContainer ) {
2052  auto ntrk = vertex.selectedTrackIndices.size() + vertex.associatedTrackIndices.size();
2053  if( vertex.isGood && ntrk >= 2 ) {
2054  dynamic_cast<TH2F*>( m_hists["vertexYieldNtrk"] )->Fill( ntrk, m_vertexingAlgorithmStep );
2055  dynamic_cast<TH2F*>( m_hists["vertexYieldChi2"] )->Fill( vertex.Chi2/(vertex.ndof() + AlgConsts::infinitesimal), m_vertexingAlgorithmStep );
2056  }
2057  }
2058  m_hists["vertexYieldNtrk"]->GetYaxis()->SetBinLabel( m_vertexingAlgorithmStep+1, name.c_str() );
2059  m_hists["vertexYieldChi2"]->GetYaxis()->SetBinLabel( m_vertexingAlgorithmStep+1, name.c_str() );
2060 
2061 
2062  if( !final ) return StatusCode::SUCCESS;
2063 
2064  for( auto& vertex : *workVerticesContainer ) {
2065  auto ntrk = vertex.selectedTrackIndices.size() + vertex.associatedTrackIndices.size();
2066  if( vertex.isGood && ntrk >= 2 ) {
2067  m_hists["finalVtxNtrk"] ->Fill( ntrk );
2068  m_hists["finalVtxR"] ->Fill( vertex.vertex.perp() );
2069  dynamic_cast<TH2F*>( m_hists["finalVtxNtrkR"] )->Fill( ntrk, vertex.vertex.perp() );
2070  }
2071  }
2072 
2073  return StatusCode::SUCCESS;
2074  }
2075 
2076  //____________________________________________________________________________________________________
2077  bool VrtSecInclusive::getSVImpactParameters(const xAOD::TrackParticle* trk, const Amg::Vector3D& vertex,
2078  std::vector<double>& impactParameters,
2079  std::vector<double>& impactParErrors){
2080 
2081  impactParameters.clear();
2082  impactParErrors.clear();
2083 
2084  if( m_jp.trkExtrapolator==1 ){
2085  m_fitSvc->VKalGetImpact(trk, vertex, static_cast<int>( trk->charge() ), impactParameters, impactParErrors);
2086  }
2087  else if( m_jp.trkExtrapolator==2 ){
2088  auto sv_perigee = m_trackToVertexTool->perigeeAtVertex(Gaudi::Hive::currentContext(), *trk, vertex );
2089  if( !sv_perigee ) return false;
2090  impactParameters.push_back(sv_perigee->parameters() [Trk::d0]);
2091  impactParameters.push_back(sv_perigee->parameters() [Trk::z0]);
2092  impactParErrors.push_back((*sv_perigee->covariance())( Trk::d0, Trk::d0 ));
2093  impactParErrors.push_back((*sv_perigee->covariance())( Trk::z0, Trk::z0 ));
2094  }
2095  else{
2096  ATH_MSG_WARNING( " > " << __FUNCTION__ << ": Unknown track extrapolator " << m_jp.trkExtrapolator );
2097  return false;
2098  }
2099 
2100  return true;
2101 
2102  } // getSVImpactParameters
2103 
2104 
2105 
2106 } // 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:343
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:349
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:342
maxValue
#define maxValue(current, test)
Definition: CompoundLayerMaterialCreator.h:22
VKalVrtAthena::VrtSecInclusive::WrkVrt::isGood
bool isGood
Definition: VrtSecInclusive.h:340
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:344
Trk::MappedVertex
Definition: IVertexMapper.h:23
Trk::z0
@ z0
Definition: ParamDefs.h:64
VKalVrtAthena::VrtSecInclusive::track_summary_properties
Definition: VrtSecInclusive.h:492
xAOD::numberOfPixelHits
@ numberOfPixelHits
these are the pixel hits, including the b-layer [unit8_t].
Definition: TrackingPrimitives.h:259
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:348
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:572
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:345
VKalVrtAthena::VrtSecInclusive::WrkVrt::Chi2
double Chi2
VKalVrt fit covariance.
Definition: VrtSecInclusive.h:346
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:351
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:228
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.
VKalVrtAthena::VrtSecInclusive::WrkVrt::selectedTrackIndices
std::deque< long int > selectedTrackIndices
flagged true for good vertex candidates
Definition: VrtSecInclusive.h:341
VKalVrtAthena::VrtSecInclusive::WrkVrt
Definition: VrtSecInclusive.h:339
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
VKalVrtAthena::VrtSecInclusive::WrkVrt::nTracksTotal
unsigned nTracksTotal() const
Definition: VrtSecInclusive.h:356
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:352
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
CondAlgsOpts.found
int found
Definition: CondAlgsOpts.py:101
VKalVrtAthena::VrtSecInclusive::WrkVrt::fitQuality
double fitQuality() const
Definition: VrtSecInclusive.h:357
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:350
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:237
Generated on Sun Dec 22 2024 21:22:21 for ATLAS Offline Software by   doxygen 1.8.18