ATLAS Offline Software
TrigVrtSecInclusive.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 
4  * Trigger specific modification of VSI, that is aimed at reconstructing displaced vertex
5  * author Kunihiro Nagano <kunihiro.nagano@cern.ch> - KEK
6  * Takane Sano <takane.sano@cern.ch> - Kyoto University
7 */
8 #include "TrigVrtSecInclusive/KDPoint.h"
9 #include "TrigVrtSecInclusive/VtxMap.h"
10 #include "TrigVrtSecInclusive/Coordinate.h"
11 #include "TrigVrtSecInclusive/Constants.h"
12 #include "TrigVrtSecInclusive/TrigVrtSecInclusive.h"
13 
14 #include <TFile.h>
15 #include "PathResolver/PathResolver.h"
16 #include "AthenaMonitoringKernel/Monitored.h"
17 #include "TrkParticleBase/TrackParticleBase.h"
18 #include "TrkTrackLink/ITrackLink.h"
19 #include "TrkTrack/LinkToTrack.h"
20 #include "TrkParticleBase/LinkToTrackParticleBase.h"
21 #include "TrkLinks/LinkToXAODTrackParticle.h"
22 #include "TrkVKalVrtCore/PGraph.h"
23 #include "AthContainers/Decorator.h"
24 
25 #include <limits>
26 
27 namespace TrigVSI {
28 
29 TrigVrtSecInclusive::TrigVrtSecInclusive(const std::string& name,ISvcLocator* pSvcLocator)
30  : AthReentrantAlgorithm(name, pSvcLocator),
31  m_materialMapInnerFileName("MaterialMap_v3.2_Inner.root"),
32  m_materialMapInnerHistName("FinalMap_inner"),
33  m_materialMapInnerMatrixName("FoldingInfo"),
34  m_materialMapOuterFileName("MaterialMap_v3_Outer.root"),
35  m_materialMapOuterHistName("matmap_outer")
36 {}
37 
38 StatusCode TrigVrtSecInclusive::initialize()
39 {
40  // read material maps
41  if (m_doMaterialMapVeto) {
42 
43  std::string inFileName = PathResolver::find_file(m_materialMapInnerFileName, "DATAPATH");
44  if ( inFileName.empty()) {
45  ATH_MSG_ERROR("Cannot find material map inner file: " << m_materialMapInnerFileName);
46  return StatusCode::FAILURE;
47  }
48  std::unique_ptr<TFile> materialMapInnerFile = std::make_unique<TFile>(inFileName.c_str(), "READ");
49  if(materialMapInnerFile->IsOpen()){
50  materialMapInnerFile->GetObject(m_materialMapInnerHistName.c_str(), m_materialMapInner);
51  materialMapInnerFile->GetObject(m_materialMapInnerMatrixName.c_str(), m_materialMapMatrix);
52  if(m_materialMapInner) m_materialMapInner->SetDirectory(0);
53  }
54  materialMapInnerFile->Close();
55 
56  std::string outFileName = PathResolver::find_file(m_materialMapOuterFileName, "DATAPATH");
57  if ( outFileName.empty()) {
58  ATH_MSG_ERROR("Cannot find material map outer file: " << m_materialMapOuterFileName);
59  return StatusCode::FAILURE;
60  }
61  std::unique_ptr<TFile> materialMapOuterFile = std::make_unique<TFile>(outFileName.c_str(), "READ");
62  if(materialMapOuterFile->IsOpen()){
63  materialMapOuterFile->GetObject(m_materialMapOuterHistName.c_str(), m_materialMapOuter);
64  if(m_materialMapOuter) m_materialMapOuter->SetDirectory(0);
65  }
66  materialMapOuterFile->Close();
67  ATH_MSG_DEBUG("material maps are read correctly" );
68 
69  }
70 
71  // monitoring tool
72  if ( !m_monTool.empty() ) ATH_CHECK(m_monTool.retrieve());
73 
74  // vertex fitters
75  ATH_CHECK( m_fitSvc.retrieve() );
76  ATH_CHECK( m_vertexPointEstimator.retrieve() );
77 
78  // collection names
79  ATH_CHECK(m_firstPassTracksName.initialize());
80  ATH_CHECK(m_secondPassTracksName.initialize());
81  ATH_CHECK(m_vxCandidatesOutputName.initialize());
82  ATH_CHECK(m_trkPairOutputName.initialize());
83  ATH_CHECK(m_PrimaryVxInputName.initialize());
84 
85  //
86  ATH_MSG_INFO( "Initialization completed successfully" );
87  return StatusCode::SUCCESS;
88 }
89 
90 StatusCode TrigVrtSecInclusive::execute(const EventContext& ctx) const
91 {
92  // Setup output container for vertex
93  SG::WriteHandle<xAOD::VertexContainer> outputVertices (m_vxCandidatesOutputName, ctx);
94 
95  std::unique_ptr<xAOD::VertexContainer> theXAODContainer = std::make_unique<xAOD::VertexContainer>();
96  std::unique_ptr<xAOD::VertexAuxContainer> theXAODAuxContainer = std::make_unique<xAOD::VertexAuxContainer>();
97  theXAODContainer->setStore( theXAODAuxContainer.get() );
98  xAODContainers theXAODContainers = std::make_pair( std::move(theXAODContainer), std::move(theXAODAuxContainer) );
99 
100  // Setup output container for track pair
101  SG::WriteHandle<xAOD::VertexContainer> outputTrkPairs (m_trkPairOutputName, ctx);
102 
103  std::unique_ptr<xAOD::VertexContainer> theXAODTrkPairContainer = std::make_unique<xAOD::VertexContainer>();
104  std::unique_ptr<xAOD::VertexAuxContainer> theXAODTrkPairAuxContainer = std::make_unique<xAOD::VertexAuxContainer>();
105  theXAODTrkPairContainer->setStore( theXAODTrkPairAuxContainer.get() );
106  xAODContainers theXAODTrkPairContainers = std::make_pair( std::move(theXAODTrkPairContainer), std::move(theXAODTrkPairAuxContainer) );
107 
108  // retrieve primary vertices
109  //_________________________________________________________________________
110  SG::ReadHandle<xAOD::VertexContainer> vtxCont(m_PrimaryVxInputName, ctx);
111  const xAOD::Vertex* privtx = nullptr;
112  if( vtxCont.isValid() ){
113  const xAOD::VertexContainer* vertex_container = vtxCont.get();
114  privtx = vertex_container->at(0);
115  if ( !( privtx->vertexType() == xAOD::VxType::PriVtx
116  && privtx->nTrackParticles() >= 2 ) ){
117  ATH_MSG_WARNING( "Illed primary vertex, keeping null privtx" );
118  privtx = nullptr;
119  } else {
120  ATH_MSG_DEBUG( "primary vertex successfully retrieved" );
121  }
122  }
123  else {
124  ATH_MSG_WARNING( "couldn't retrieve primary vertex, keeping null privtx" );
125  }
126 
127  // retrieve the tracks
128  //_________________________________________________________________________
129  SG::ReadHandle<xAOD::TrackParticleContainer> firstPassTrackParticleCollection(m_firstPassTracksName, ctx);
130  SG::ReadHandle<xAOD::TrackParticleContainer> secondPassTrackParticleCollection(m_secondPassTracksName, ctx);
131  bool isFirstTrackValid = firstPassTrackParticleCollection.isValid();
132  bool isSecondTrackValid = secondPassTrackParticleCollection.isValid();
133  if( ! isFirstTrackValid ) ATH_MSG_WARNING( "couldn't retrieve first pass tracks: " << m_firstPassTracksName);
134  if( ! isSecondTrackValid ) ATH_MSG_WARNING( "couldn't retrieve second pass tracks: " << m_secondPassTracksName);
135 
136  // monitoring
137  //_________________________________________________________________________
138  auto timerTrackSel = Monitored::Timer<std::chrono::nanoseconds> ("TIME_TrackSel");
139  auto timerTwoTrackVertex = Monitored::Timer<std::chrono::milliseconds> ("TIME_TwoTrackVertex");
140  auto timerMapClustering = Monitored::Timer<std::chrono::microseconds> ("TIME_MapClustering");
141  auto timerNTrackVertex = Monitored::Timer<std::chrono::milliseconds> ("TIME_NTrackVertex");
142  auto timerNTrackVtxOffVSI = Monitored::Timer<std::chrono::milliseconds> ("TIME_NTrackVtxOffVSI");
143  auto timerCleanUp = Monitored::Timer<std::chrono::nanoseconds> ("TIME_CleanUp");
144  auto timerWriteVertices = Monitored::Timer<std::chrono::milliseconds> ("TIME_WriteVertices");
145  auto timerOverall = Monitored::Timer<std::chrono::milliseconds> ("TIME_Overall");
146 
147  auto monTime = Monitored::Group(m_monTool, timerTrackSel, timerTwoTrackVertex, timerMapClustering, timerNTrackVertex, timerNTrackVtxOffVSI, timerCleanUp, timerWriteVertices, timerOverall);
148 
149  timerOverall.start();
150 
151  // Reset incompatible pair list and selectedTracks
152  //_________________________________________________________________________
153  std::vector<std::pair<size_t,size_t>> v_incomp;
154  std::vector<const xAOD::TrackParticle*> v_selectedTracks;
155 
156  // tracks selection
157  //_________________________________________________________________________
158  timerTrackSel.start();
159 
160  if( isFirstTrackValid ) {
161  if( isSecondTrackValid ) {
162  ATH_CHECK( trackSelection( firstPassTrackParticleCollection.cptr(), secondPassTrackParticleCollection.cptr(), v_selectedTracks ) );
163  }
164  else {
165  ATH_CHECK( trackSelection( firstPassTrackParticleCollection.cptr(), nullptr, v_selectedTracks ) );
166  }
167  }
168  timerTrackSel.stop();
169 
170 
171  // Vertex reconstruction
172  //_________________________________________________________________________
173 
174  // Find compatible track pairs and record as di-trakc vertex
175  // ===============================
176  WrkVrtContainer v_diwrkvrt;
177  v_diwrkvrt.reserve(5000);
178 
179  if ( m_vtxAlgorithm == 0 ) {
180  // TrigVSI
181  timerTwoTrackVertex.start();
182  ATH_CHECK( findDiTrackVertex( v_diwrkvrt, v_incomp, v_selectedTracks, ctx, privtx ) );
183  timerTwoTrackVertex.stop();
184  } else {
185  // Offline VSI like
186  timerTwoTrackVertex.start();
187  ATH_CHECK( findDiTrackVertexVSI( v_diwrkvrt, v_incomp, v_selectedTracks, ctx, privtx ) );
188  timerTwoTrackVertex.stop();
189  }
190 
191  // Reconstruct N track vertices
192  // ===============================
193  // Ensure that the adress of each di-track wrkvrt can't be changed.
194  const WrkVrtContainer v_diTrkVertices = std::move(v_diwrkvrt);
195 
196  // define vertex map
197  namespace vsic = TrigVSI::AlgConsts;
198  TrigVSI::VtxMap<WrkVrt, TrigVSI::Coordinate::Pseudo> vtxmap( std::make_unique<TH3D>("TrigVrtSecInclusive_VtxMap", "TrigVrtSecInclusive_VtxMap", vsic::binNumR, vsic::IDrInner, vsic::IDrInner+vsic::mapBinWid*vsic::binNumR, vsic::nEta, -vsic::maxEta, vsic::maxEta, vsic::nPhi, -TMath::Pi(), TMath::Pi()) );
199 
200  // Container for reconstructed N-track vertices
201  WrkVrtContainer v_wrkvrt;
202  v_wrkvrt.reserve(1000);
203 
204  if ( m_vtxAlgorithm == 0 ) {
205  // TrigVSI
206  int trkpair_cnt = 0;
207 
208  // Fill vertex map
209  timerMapClustering.start();
210  for (auto wrkvrt_iter = v_diTrkVertices.begin(); wrkvrt_iter != v_diTrkVertices.end(); ++wrkvrt_iter ) {
211  if ( wrkvrt_iter->isGood && wrkvrt_iter->cuts.isFullPass() ) {
212  vtxmap.Fill( &(*wrkvrt_iter) );
213  trkpair_cnt++;
214  }
215  }
216  vtxmap.lock();
217  ATH_MSG_DEBUG( "filled vertex map with " << trkpair_cnt << " pairs" );
218  vtxmap.ClusterizeCells(1.,1);
219  timerMapClustering.stop();
220  ATH_MSG_DEBUG( "vertex map clustering successfully completed : " << vtxmap.nClusters() << " clusters" );
221 
222  // N track vertexing
223  timerNTrackVertex.start();
224  ATH_CHECK( findNTrackVertex( v_wrkvrt, vtxmap, v_selectedTracks, ctx ) );
225  timerNTrackVertex.stop();
226 
227  // Cleanup vertex fully included in other vertex.
228  timerCleanUp.start();
229  size_t n_pre_clean = v_wrkvrt.size();
230  ATH_CHECK( cleanUp(v_wrkvrt) );
231  size_t n_post_clean = v_wrkvrt.size();
232  timerCleanUp.stop();
233  ATH_MSG_DEBUG( "cleaned up " << n_pre_clean - n_post_clean << " out of " << n_pre_clean );
234 
235  } else {
236  // Offline VSI like
237  timerNTrackVtxOffVSI.start();
238  ATH_CHECK( findNtrackVerticesVSI( v_wrkvrt, v_incomp, v_selectedTracks, ctx ) );
239  timerNTrackVtxOffVSI.stop();
240 
241  }
242 
243  // Fill Vertex in the VertexContainer
244  timerWriteVertices.start();
245 
246  ATH_CHECK( fillVtxContainer( theXAODContainers, v_wrkvrt, v_selectedTracks ) );
247  if (m_recordTrkPair) {
248  ATH_CHECK( fillVtxContainer( theXAODTrkPairContainers, v_diTrkVertices, v_selectedTracks ) );
249  }
250 
251  // record AOD object
252  ATH_MSG_DEBUG( "Record vertices into container." );
253  ATH_CHECK(outputVertices.record(std::move(theXAODContainers.first), std::move(theXAODContainers.second)));
254  ATH_CHECK(outputTrkPairs.record(std::move(theXAODTrkPairContainers.first), std::move(theXAODTrkPairContainers.second)));
255  ATH_MSG_DEBUG( "Recorded Vertices with key: " << m_vxCandidatesOutputName.key() );
256 
257  timerWriteVertices.stop();
258 
259  timerOverall.stop();
260  ATH_MSG_DEBUG( "::execute() finished successfully." );
261 
262  return StatusCode::SUCCESS;
263 }
264 
265 
266 StatusCode TrigVrtSecInclusive::trackSelection( const xAOD::TrackParticleContainer* firstPassTrackParticles, const xAOD::TrackParticleContainer* secondPassTrackParticles, std::vector<const xAOD::TrackParticle*>& selectedTracks ) const
267 {
268  // Pack TrackParticleContainers to one vector
269  std::vector<const xAOD::TrackParticleContainer*> v_containers;
270  if( firstPassTrackParticles != nullptr ) v_containers.push_back(firstPassTrackParticles);
271  if( secondPassTrackParticles != nullptr ) v_containers.push_back(secondPassTrackParticles);
272 
273  typedef DataVector<xAOD::TrackParticle>::const_iterator TrackParticleDataVecIter;
274 
275  for(unsigned int i=0; i<v_containers.size(); i++) {
276  const xAOD::TrackParticleContainer* trackParticles = v_containers[i];
277  size_t n_trk = 0;
278  for (TrackParticleDataVecIter itr = trackParticles->begin(); itr != trackParticles->end(); ++itr) {
279  if( selectTrack(*itr) ) { selectedTracks.push_back(*itr); n_trk++; }
280  }
281  ATH_MSG_DEBUG("Of " << trackParticles->size() << " tracks " << n_trk << " survived the preselection.");
282  }
283  ATH_MSG_DEBUG(selectedTracks.size() << " tracks in total passed the selection.");
284  return StatusCode::SUCCESS;
285 }
286 
287 
288 StatusCode TrigVrtSecInclusive::fillVtxContainer( xAODContainers& theXAODContainers, const WrkVrtContainer& workVerticesContainer, std::vector<const xAOD::TrackParticle*>& selectedTracks ) const
289 {
290  xAOD::VertexContainer* theVertexContainer = theXAODContainers.first.get();
291 
292  for (const auto& wrkvrt : workVerticesContainer) {
293 
294  if ( !wrkvrt.isGood && !wrkvrt.isPair ) continue;
295 
296  std::vector<const xAOD::TrackParticle*> baseTracks;
297  for ( auto i : wrkvrt.selectedTrackIndices() ) {
298  baseTracks.emplace_back( selectedTracks.at(i) );
299  }
300 
301  // Create a xAOD::Vertex instance
302  xAOD::Vertex* vertex = new xAOD::Vertex();
303  vertex->makePrivateStore();
304  theVertexContainer->emplace_back( vertex );
305 
306  for( auto *trk: baseTracks ) {
307  // Acquire link to the track
308  ElementLink<xAOD::TrackParticleContainer> trackElementLink( *( dynamic_cast<const xAOD::TrackParticleContainer*>( trk->container() ) ), trk->index() );
309  // Register link to the vertex
310  vertex->addTrackAtVertex( trackElementLink, 1. );
311  }
312 
313  vertex->setVertexType( xAOD::VxType::SecVtx );
314  vertex->setPosition( wrkvrt.vertex );
315  vertex->setFitQuality( wrkvrt.chi2, 1 ); // Ndof is always 1
316 
317  static const SG::Accessor<float> vsi_massAcc("vsi_mass");
318  static const SG::Accessor<float> vsi_pTAcc("vsi_pT");
319  static const SG::Accessor<float> vsi_chargeAcc("vsi_charge");
320  static const SG::Accessor<char> vsi_isFakeAcc("vsi_isFake");
321  vsi_massAcc(*vertex) = wrkvrt.vertexMom.M();
322  vsi_pTAcc(*vertex) = wrkvrt.vertexMom.Perp();
323  vsi_chargeAcc(*vertex) = wrkvrt.charge;
324  vsi_isFakeAcc(*vertex) = 0;
325 
326  static const SG::Accessor<float> vsi_twoCirc_drAcc("vsi_twoCirc_dr");
327  static const SG::Accessor<float> vsi_twoCirc_dphiAcc("vsi_twoCirc_dphi");
328  static const SG::Accessor<float> vsi_twoCirc_int_rAcc("vsi_twoCirc_int_r");
329  static const SG::Accessor<float> vsi_vrtFast_rAcc("vsi_vrtFast_r");
330  static const SG::Accessor<float> vsi_vrtFast_etaAcc("vsi_vrtFast_eta");
331  static const SG::Accessor<float> vsi_vrtFast_phiAcc("vsi_vrtFast_phi");
332  static const SG::Accessor< std::vector<float> > vsi_vrtFast_trkd0Acc("vsi_vrtFast_trkd0");
333  static const SG::Accessor< std::vector<float> > vsi_vrtFast_trkz0Acc("vsi_vrtFast_trkz0");
334  static const SG::Accessor<float> vsi_vrtFit_rAcc("vsi_vrtFit_r");
335  static const SG::Accessor<float> vsi_vrtFit_chi2Acc("vsi_vrtFit_chi2");
336  static const SG::Accessor<float> vsi_vPosAcc("vsi_vPos");
337  static const SG::Accessor<float> vsi_vPosMomAngTAcc("vsi_vPosMomAngT");
338  static const SG::Accessor<float> vsi_dphi1Acc("vsi_dphi1");
339  static const SG::Accessor<float> vsi_dphi2Acc("vsi_dphi2");
340  static const SG::Accessor<char> vsi_isPassMMVAcc("vsi_isPassMMV");
341  vsi_twoCirc_drAcc(*vertex) = wrkvrt.param.twoCirc_dr;
342  vsi_twoCirc_dphiAcc(*vertex) = wrkvrt.param.twoCirc_dphi;
343  vsi_twoCirc_int_rAcc(*vertex) = wrkvrt.param.twoCirc_int_r;
344  vsi_vrtFast_rAcc(*vertex) = wrkvrt.param.vrtFast_r;
345  vsi_vrtFast_etaAcc(*vertex) = wrkvrt.param.vrtFast_eta;
346  vsi_vrtFast_phiAcc(*vertex) = wrkvrt.param.vrtFast_phi;
347  vsi_vrtFast_trkd0Acc(*vertex) = wrkvrt.param.vrtFast_trkd0;
348  vsi_vrtFast_trkz0Acc(*vertex) = wrkvrt.param.vrtFast_trkz0;
349  vsi_vrtFit_rAcc(*vertex) = wrkvrt.vertex.perp();
350  vsi_vrtFit_chi2Acc(*vertex) = wrkvrt.chi2;
351  vsi_vPosAcc(*vertex) = wrkvrt.param.vPos;
352  vsi_vPosMomAngTAcc(*vertex) = wrkvrt.param.vPosMomAngT;
353  vsi_dphi1Acc(*vertex) = wrkvrt.param.dphi1;
354  vsi_dphi2Acc(*vertex) = wrkvrt.param.dphi2;
355  vsi_isPassMMVAcc(*vertex) = wrkvrt.param.isPassMMV ? 1 : 0;
356 
357  static const SG::Accessor<char> vsi_trkd0cutAcc("vsi_trkd0cut");
358  static const SG::Accessor<char> vsi_twoCircErrcutAcc("vsi_twoCircErrcut");
359  static const SG::Accessor<char> vsi_twoCircRcutAcc("vsi_twoCircRcut");
360  static const SG::Accessor<char> vsi_fastErrcutAcc("vsi_fastErrcut");
361  static const SG::Accessor<char> vsi_fastRcutAcc("vsi_fastRcut");
362  static const SG::Accessor<char> vsi_fitErrcutAcc("vsi_fitErrcut");
363  static const SG::Accessor<char> vsi_chi2cutAcc("vsi_chi2cut");
364  vsi_trkd0cutAcc(*vertex) = wrkvrt.cuts.trkd0cut ? 1 : 0;
365  vsi_twoCircErrcutAcc(*vertex) = wrkvrt.cuts.twoCircErrcut ? 1 : 0;
366  vsi_twoCircRcutAcc(*vertex) = wrkvrt.cuts.twoCircRcut ? 1 : 0;
367  vsi_fastErrcutAcc(*vertex) = wrkvrt.cuts.fastErrcut ? 1 : 0;
368  vsi_fastRcutAcc(*vertex) = wrkvrt.cuts.fastRcut ? 1 : 0;
369  vsi_fitErrcutAcc(*vertex) = wrkvrt.cuts.fitErrcut ? 1 : 0;
370  vsi_chi2cutAcc(*vertex) = wrkvrt.cuts.chi2cut ? 1 : 0;
371 
372  }
373  //
374  return StatusCode::SUCCESS;
375 }
376 
377 
378 StatusCode TrigVrtSecInclusive::findDiTrackVertex
379 ( WrkVrtContainer& workVerticesContainer, std::vector<std::pair<size_t,size_t>>& incomp, std::vector<const xAOD::TrackParticle*>& selectedTracks, const EventContext& ctx, const xAOD::Vertex *primVertex ) const
380 {
381  // monitoring
382  auto timerTwoCircIntsect = Monitored::Timer<std::chrono::nanoseconds> ("TIME_TwoCircIntsect");
383  auto timerVrtFitFast = Monitored::Timer<std::chrono::nanoseconds> ("TIME_VrtFitFast");
384  auto timerVrtFit = Monitored::Timer<std::chrono::microseconds>("TIME_VrtFit");
385 
386  std::vector<float> mnt_pair_dphi;
387  std::vector<float> mnt_pair_dr;
388  std::vector<float> mnt_intersect_r;
389  std::vector<float> mnt_init_r;
390  std::vector<float> mnt_init_trkd0;
391  std::vector<float> mnt_init_trkz0;
392  std::vector<float> mnt_vtxfit_chi2;
393  std::vector<float> mnt_vtxfit_r;
394  std::vector<float> mnt_vtx_chi2;
395  std::vector<float> mnt_vtx_mass;
396  std::vector<float> mnt_vtx_mass_high;
397  std::vector<float> mnt_vtx_pt;
398  std::vector<float> mnt_vtx_charge;
399  std::vector<float> mnt_vtx_r;
400  auto mon_pair_dphi = Monitored::Collection("pair_dphi", mnt_pair_dphi);
401  auto mon_pair_dr = Monitored::Collection("pair_dr", mnt_pair_dr);
402  auto mon_intersect_r = Monitored::Collection("intersect_r", mnt_intersect_r);
403  auto mon_init_r = Monitored::Collection("init_r", mnt_init_r);
404  auto mon_init_trkd0 = Monitored::Collection("init_trkd0", mnt_init_trkd0);
405  auto mon_init_trkz0 = Monitored::Collection("init_trkz0", mnt_init_trkz0);
406  auto mon_vtxfit_chi2 = Monitored::Collection("vtxfit_chi2", mnt_vtxfit_chi2);
407  auto mon_vtxfit_r = Monitored::Collection("vtxfit_r", mnt_vtxfit_r);
408  auto mon_vtx_chi2 = Monitored::Collection("vtx_chi2", mnt_vtx_chi2);
409  auto mon_vtx_mass = Monitored::Collection("vtx_mass", mnt_vtx_mass);
410  auto mon_vtx_mass_high= Monitored::Collection("vtx_mass_high",mnt_vtx_mass);
411  auto mon_vtx_pt = Monitored::Collection("vtx_pt", mnt_vtx_pt);
412  auto mon_vtx_charge = Monitored::Collection("vtx_charge", mnt_vtx_charge);
413  auto mon_vtx_r = Monitored::Collection("vtx_r", mnt_vtx_r);
414  auto monVertex = Monitored::Group(m_monTool, mon_pair_dphi, mon_pair_dr, mon_intersect_r, mon_init_r, mon_init_trkd0, mon_init_trkz0, mon_vtxfit_chi2, mon_vtxfit_r,
415  mon_vtx_chi2,mon_vtx_mass,mon_vtx_mass_high,mon_vtx_pt,mon_vtx_charge,mon_vtx_r);
416 
417  // vertexing
418  unsigned int nPairsAll = 0;
419  unsigned int nPairsTrkd0 = 0;
420  unsigned int nPairsIntersect = 0;
421  unsigned int nPairsIntersectPos = 0;
422  unsigned int nPairsInitPos = 0;
423  unsigned int nPairsInitTrkd0z0 = 0;
424  unsigned int nPairsVtxFit = 0;
425  unsigned int nPairsVtxChi2 = 0;
426  unsigned int nPairsVtxComp = 0;
427  unsigned int nPairsVtxMatveto = 0;
428 
429  for( auto itrkIter = selectedTracks.begin(); itrkIter != selectedTracks.end(); ++itrkIter ) {
430  const xAOD::TrackParticle* itrk = *itrkIter;
431  for( auto jtrkIter = std::next(itrkIter); jtrkIter != selectedTracks.end(); ++jtrkIter ) { // for starts from here
432  const xAOD::TrackParticle* jtrk = *jtrkIter;
433  nPairsAll++;
434 
435  auto monTimeStep = Monitored::Group(m_monTool, timerTwoCircIntsect, timerVrtFitFast, timerVrtFit);
436 
437  //
438  const int itrk_id = std::distance(selectedTracks.begin(), itrkIter);
439  const int jtrk_id = std::distance(selectedTracks.begin(), jtrkIter);
440 
441  WrkVrt::AlgCuts wrkcuts;
442  WrkVrt::AlgParam wrkprm;
443  WrkVrt wrkvrt;
444 
445  // Set track indices
446  wrkvrt.selectedTrackIndices().clear();
447  wrkvrt.selectedTrackIndices().emplace_back(itrk_id);
448  wrkvrt.selectedTrackIndices().emplace_back(jtrk_id);
449 
450  // trk d0
451  wrkcuts.trkd0cut = std::abs( itrk->d0() ) < m_twoTrkVtxFormingD0Cut && std::abs( jtrk->d0() ) < m_twoTrkVtxFormingD0Cut;
452  nPairsTrkd0++;
453 
454  // Attempt to think the combination is incompatible by default
455  incomp.emplace_back( std::pair<size_t, size_t>(itrk_id, jtrk_id) );
456 
457  // two circles intersection
458  const Trk::Perigee& perigee1 = itrk->perigeeParameters();
459  const Trk::Perigee& perigee2 = jtrk->perigeeParameters();
460  int flag = 0;
461  int errorcode = 0;
462  InDet::VertexPointEstimator::Values_t decors;
463  timerTwoCircIntsect.start();
464  Amg::Vector3D circIntersect = m_vertexPointEstimator->getCirclesIntersectionPoint(&perigee1, &perigee2, flag, errorcode, decors);
465  timerTwoCircIntsect.stop();
466  if (errorcode != 0) {
467  ATH_MSG_VERBOSE( ": two circles intersect failed");
468  wrkcuts.twoCircErrcut = true;
469  wrkvrt.param = wrkprm;
470  wrkvrt.cuts = wrkcuts;
471  wrkvrt.isPair = true;
472 
473  workVerticesContainer.push_back( std::move(wrkvrt) );
474  continue;
475  }
476  float dphi = std::abs(decors["deltaPhiTracks"]);
477  float dr = std::abs(decors["DR1R2"]);
478  float intersect_r = circIntersect.perp();
479  mnt_pair_dphi.push_back(dphi);
480  mnt_pair_dr.push_back(dr);
481  mnt_intersect_r.push_back(intersect_r);
482 
483  wrkprm.twoCirc_dphi = dphi;
484  wrkprm.twoCirc_dr = dr;
485  wrkprm.twoCirc_int_r = intersect_r;
486 
487  nPairsIntersect++;
488  wrkcuts.twoCircRcut = intersect_r > m_maxR;
489  if (m_doTwoCircRCut && intersect_r > m_maxR) continue;
490  nPairsIntersectPos++;
491 
492  // initial fast fitting
493  std::vector<const xAOD::TrackParticle*> baseTracks;
494  baseTracks.emplace_back( itrk );
495  baseTracks.emplace_back( jtrk );
496  Amg::Vector3D initVertex;
497 
498  timerVrtFitFast.start();
499  auto fitterState = m_fitSvc->makeState(ctx);
500  m_fitSvc->setApproximateVertex( circIntersect.x(), circIntersect.y(), circIntersect.z(), *fitterState );
501 
502  StatusCode sc = m_fitSvc->VKalVrtFitFast( baseTracks, initVertex, *fitterState );
503  timerVrtFitFast.stop();
504  if( sc.isFailure() ) {
505  ATH_MSG_VERBOSE( ": fast crude estimation fails ");
506  wrkcuts.fastErrcut = true;
507  wrkvrt.param = wrkprm;
508  wrkvrt.cuts = wrkcuts;
509  wrkvrt.isPair = true;
510 
511  workVerticesContainer.push_back( std::move(wrkvrt) );
512  continue;
513  }
514 
515  // position of initial fast fitting
516  mnt_init_r.push_back(initVertex.perp());
517  wrkprm.vrtFast_r = initVertex.perp();
518  wrkprm.vrtFast_eta = initVertex.eta();
519  wrkprm.vrtFast_phi = initVertex.phi();
520 
521  wrkcuts.fastRcut = initVertex.perp() > m_maxR;
522  if( m_doFastRCut && initVertex.perp() > m_maxR ) continue;
523  nPairsInitPos++;
524 
525  // impact parameter on initial fast fitting
526  std::vector<double> impactParameters;
527  std::vector<double> impactParErrors;
528  if( ! m_fitSvc->VKalGetImpact(itrk, initVertex, static_cast<int>( itrk->charge() ), impactParameters, impactParErrors) ) continue;
529  const auto roughD0_itrk = impactParameters.at(TrkParameter::k_d0);
530  const auto roughZ0_itrk = impactParameters.at(TrkParameter::k_z0);
531  mnt_init_trkd0.push_back(std::abs(roughD0_itrk));
532  mnt_init_trkz0.push_back(std::abs(roughZ0_itrk));
533  wrkprm.vrtFast_trkd0.push_back(std::abs(roughD0_itrk));
534  wrkprm.vrtFast_trkz0.push_back(std::abs(roughZ0_itrk));
535 
536  if( ! m_fitSvc->VKalGetImpact(jtrk, initVertex, static_cast<int>( jtrk->charge() ), impactParameters, impactParErrors) ) continue;
537  const auto roughD0_jtrk = impactParameters.at(TrkParameter::k_d0);
538  const auto roughZ0_jtrk = impactParameters.at(TrkParameter::k_z0);
539  mnt_init_trkd0.push_back(std::abs(roughD0_jtrk));
540  mnt_init_trkz0.push_back(std::abs(roughZ0_jtrk));
541  wrkprm.vrtFast_trkd0.push_back(std::abs(roughD0_jtrk));
542  wrkprm.vrtFast_trkz0.push_back(std::abs(roughZ0_jtrk));
543 
544  if ( std::min(roughD0_itrk, roughD0_jtrk) > m_fastD0minCut || std::abs(roughD0_itrk - roughD0_jtrk) > m_fastD0deltaCut ) continue;
545  if ( std::min(roughZ0_itrk, roughZ0_jtrk) > m_fastZ0minCut || std::abs(roughZ0_itrk - roughZ0_jtrk) > m_fastZ0deltaCut ) continue;
546 
547  nPairsInitTrkd0z0++;
548 
549  // Vertex VKal Fitting
550  std::vector<const xAOD::NeutralParticle*> dummyNeutrals;
551  timerVrtFit.start();
552  m_fitSvc->setApproximateVertex( initVertex.x(), initVertex.y(), initVertex.z(), *fitterState );
553  sc = m_fitSvc->VKalVrtFit( baseTracks,
554  dummyNeutrals,
555  wrkvrt.vertex, wrkvrt.vertexMom, wrkvrt.charge,
556  wrkvrt.vertexCov, wrkvrt.chi2PerTrk,
557  wrkvrt.trkAtVrt, wrkvrt.chi2, *fitterState, false );
558  timerVrtFit.stop();
559  if( sc.isFailure() ) {
560  wrkcuts.fitErrcut = true;
561  wrkvrt.param = wrkprm;
562  wrkvrt.cuts = wrkcuts;
563  wrkvrt.isPair = true;
564 
565  workVerticesContainer.push_back( std::move(wrkvrt) );
566  continue;
567  }
568  nPairsVtxFit++;
569 
570  // Chi2 cut
571  mnt_vtxfit_chi2.push_back(wrkvrt.chi2);
572  if( wrkvrt.chi2 > m_selVrtChi2Cut) {
573  ATH_MSG_VERBOSE( ": failed to pass chi2 threshold." );
574  wrkcuts.chi2cut = true;
575  wrkvrt.param = wrkprm;
576  wrkvrt.cuts = wrkcuts;
577  wrkvrt.isPair = true;
578 
579  workVerticesContainer.push_back( std::move(wrkvrt) );
580  continue;
581  }
582  nPairsVtxChi2++;
583  mnt_vtxfit_r.push_back(wrkvrt.vertex.perp());
584 
585  // Compatibility to the primary vertex
586  Amg::Vector3D vDist = (primVertex!=nullptr)? wrkvrt.vertex - primVertex->position() : wrkvrt.vertex;
587  const double vPos = ( vDist.x()*wrkvrt.vertexMom.Px()+vDist.y()*wrkvrt.vertexMom.Py()+vDist.z()*wrkvrt.vertexMom.Pz() )/wrkvrt.vertexMom.Rho();
588  const double vPosMomAngT = ( vDist.x()*wrkvrt.vertexMom.Px()+vDist.y()*wrkvrt.vertexMom.Py() ) / vDist.perp() / wrkvrt.vertexMom.Pt();
589 
590  double dphi1 = vDist.phi() - itrk->phi(); while( dphi1 > TMath::Pi() ) { dphi1 -= TMath::TwoPi(); } while( dphi1 < -TMath::Pi() ) { dphi1 += TMath::TwoPi(); }
591  double dphi2 = vDist.phi() - jtrk->phi(); while( dphi2 > TMath::Pi() ) { dphi2 -= TMath::TwoPi(); } while( dphi2 < -TMath::Pi() ) { dphi2 += TMath::TwoPi(); }
592 
593  wrkprm.vPos = vPos;
594  wrkprm.vPosMomAngT = vPosMomAngT;
595  wrkprm.dphi1 = dphi1;
596  wrkprm.dphi2 = dphi2;
597 
598  if (m_doPVCompatibilityCut) {
599 
600  if( std::cos( dphi1 ) < m_dphiPVCut && std::cos( dphi2 ) < m_dphiPVCut ) {
601  ATH_MSG_DEBUG( ": failed to pass the vPos cut. (both tracks are opposite against the vertex pos)" );
602  continue;
603  }
604  if( vPosMomAngT < m_dphiPVCut ) {
605  ATH_MSG_DEBUG( ": failed to pass the vPos cut. (pos-mom directions are opposite)" );
606  continue;
607  }
608 
609  if( vPos < m_pvCompatibilityCut ) {
610  ATH_MSG_DEBUG( ": failed to pass the vPos cut." );
611  continue;
612  }
613 
614  }
615  nPairsVtxComp++;
616 
617  // material map veto
618  wrkprm.isPassMMV = true;
619  if (m_doMaterialMapVeto) {
620  if (wrkvrt.vertex.perp() > 150) {
621  wrkprm.isPassMMV = ( m_materialMapOuter->GetBinContent(m_materialMapOuter->FindFixBin( wrkvrt.vertex.perp(), wrkvrt.vertex.phi(), wrkvrt.vertex.z() ) ) == 0);
622  }
623  else {
624  const TMatrixT<double>& mmm = *m_materialMapMatrix;
625  for (int i=0;i<5;i++){
626  if (wrkvrt.vertex.perp() < mmm[i][0]) {
627  float test_x = wrkvrt.vertex.x() + mmm[i][1];
628  float test_y = wrkvrt.vertex.y() + mmm[i][2];
629  double calc_phi = std::fmod( TMath::ATan2(test_y,test_x),TMath::Pi()/mmm[i][3] );
630  if (calc_phi <0) calc_phi = calc_phi + TMath::Pi()/mmm[i][3];
631  wrkprm.isPassMMV = ( m_materialMapInner->GetBinContent(m_materialMapInner->FindFixBin(sqrt(test_x*test_x + test_y*test_y), calc_phi, wrkvrt.vertex.z() ) ) == 0);
632  }
633  }
634  }
635  if( ! wrkprm.isPassMMV ) continue;
636  }
637  nPairsVtxMatveto++;
638 
639  // Now this vertex passed all criteria and considred to be a compatible vertices.
640  // Therefore the track pair is removed from the incompatibility list.
641  if (wrkcuts.isFullPass()) incomp.pop_back();
642 
643  wrkvrt.param = wrkprm;
644  wrkvrt.cuts = wrkcuts;
645  wrkvrt.isGood = true;
646  wrkvrt.isPair = true;
647 
648  workVerticesContainer.push_back( std::move(wrkvrt) );
649  }
650  }
651  ATH_MSG_DEBUG("Of " << nPairsAll << " pairs : trkd0" << " / " << "intersect" << " / " << "intersectPos" << " / " << "initPos" << " / " << "initD0Z0" << " / " << "vtxFit" << " / " << "vtxChi2" << " / " << "vtxComp" << " / " << "matVeto = "
652  << nPairsTrkd0 << " / " << nPairsIntersect << " / " << nPairsIntersectPos << " / " << nPairsInitPos << " / " << nPairsInitTrkd0z0 << " / " << nPairsVtxFit << " / " << nPairsVtxChi2 << " / " << nPairsVtxComp << " / " << nPairsVtxMatveto);
653  //
654  return StatusCode::SUCCESS;
655 }
656 
657 
658 StatusCode TrigVrtSecInclusive::findDiTrackVertexVSI
659 ( WrkVrtContainer& workVerticesContainer, std::vector<std::pair<size_t,size_t>>& incomp, std::vector<const xAOD::TrackParticle*>& selectedTracks, const EventContext& ctx, const xAOD::Vertex *primVertex ) const
660 {
661  // monitoring
662  auto timerTwoCircIntsect = Monitored::Timer<std::chrono::nanoseconds> ("TIME_TwoCircIntsect");
663  auto timerVrtFitFast = Monitored::Timer<std::chrono::nanoseconds> ("TIME_VrtFitFast");
664  auto timerVrtFit = Monitored::Timer<std::chrono::microseconds>("TIME_VrtFit");
665 
666  std::vector<float> mnt_pair_dphi;
667  std::vector<float> mnt_pair_dr;
668  std::vector<float> mnt_intersect_r;
669  std::vector<float> mnt_init_r;
670  std::vector<float> mnt_init_trkd0;
671  std::vector<float> mnt_init_trkz0;
672  std::vector<float> mnt_vtxfit_chi2;
673  std::vector<float> mnt_vtxfit_r;
674  std::vector<float> mnt_vtx_chi2;
675  std::vector<float> mnt_vtx_mass;
676  std::vector<float> mnt_vtx_mass_high;
677  std::vector<float> mnt_vtx_pt;
678  std::vector<float> mnt_vtx_charge;
679  std::vector<float> mnt_vtx_r;
680  auto mon_pair_dphi = Monitored::Collection("pair_dphi", mnt_pair_dphi);
681  auto mon_pair_dr = Monitored::Collection("pair_dr", mnt_pair_dr);
682  auto mon_intersect_r = Monitored::Collection("intersect_r", mnt_intersect_r);
683  auto mon_init_r = Monitored::Collection("init_r", mnt_init_r);
684  auto mon_init_trkd0 = Monitored::Collection("init_trkd0", mnt_init_trkd0);
685  auto mon_init_trkz0 = Monitored::Collection("init_trkz0", mnt_init_trkz0);
686  auto mon_vtxfit_chi2 = Monitored::Collection("vtxfit_chi2", mnt_vtxfit_chi2);
687  auto mon_vtxfit_r = Monitored::Collection("vtxfit_r", mnt_vtxfit_r);
688  auto mon_vtx_chi2 = Monitored::Collection("vtx_chi2", mnt_vtx_chi2);
689  auto mon_vtx_mass = Monitored::Collection("vtx_mass", mnt_vtx_mass);
690  auto mon_vtx_mass_high= Monitored::Collection("vtx_mass_high",mnt_vtx_mass);
691  auto mon_vtx_pt = Monitored::Collection("vtx_pt", mnt_vtx_pt);
692  auto mon_vtx_charge = Monitored::Collection("vtx_charge", mnt_vtx_charge);
693  auto mon_vtx_r = Monitored::Collection("vtx_r", mnt_vtx_r);
694  auto monVertex = Monitored::Group(m_monTool, mon_pair_dphi, mon_pair_dr, mon_intersect_r, mon_init_r, mon_init_trkd0, mon_init_trkz0, mon_vtxfit_chi2, mon_vtxfit_r,
695  mon_vtx_chi2,mon_vtx_mass,mon_vtx_mass_high,mon_vtx_pt,mon_vtx_charge,mon_vtx_r);
696 
697  // vertexing
698  const double roughD0Cut{ 100. };
699  const double roughZ0Cut{ 50. };
700 
701  unsigned int nPairsAll = 0;
702  unsigned int nPairsTrkd0 = 0;
703  unsigned int nPairsIntersect = 0;
704  unsigned int nPairsIntersectPos = 0;
705  unsigned int nPairsInitPos = 0;
706  unsigned int nPairsInitTrkd0z0 = 0;
707  unsigned int nPairsVtxFit = 0;
708  unsigned int nPairsVtxChi2 = 0;
709  unsigned int nPairsVtxComp = 0;
710 
711  for( auto itrkIter = selectedTracks.begin(); itrkIter != selectedTracks.end(); ++itrkIter ) {
712  const xAOD::TrackParticle* itrk = *itrkIter;
713  for( auto jtrkIter = std::next(itrkIter); jtrkIter != selectedTracks.end(); ++jtrkIter ) { // for starts from here
714  const xAOD::TrackParticle* jtrk = *jtrkIter;
715  nPairsAll++;
716 
717  auto monTimeStep = Monitored::Group(m_monTool, timerTwoCircIntsect, timerVrtFitFast, timerVrtFit);
718 
719  //
720  const int itrk_id = std::distance(selectedTracks.begin(), itrkIter);
721  const int jtrk_id = std::distance(selectedTracks.begin(), jtrkIter);
722 
723  WrkVrt::AlgCuts wrkcuts;
724  WrkVrt::AlgParam wrkprm;
725  WrkVrt wrkvrt;
726 
727  // Set track indices
728  wrkvrt.selectedTrackIndices().clear();
729  wrkvrt.selectedTrackIndices().emplace_back(itrk_id);
730  wrkvrt.selectedTrackIndices().emplace_back(jtrk_id);
731 
732  // trk d0
733  if( std::abs( itrk->d0() ) < m_twoTrkVtxFormingD0Cut && std::abs( jtrk->d0() ) < m_twoTrkVtxFormingD0Cut ) continue;
734  nPairsTrkd0++;
735 
736  // Attempt to think the combination is incompatible by default
737  incomp.emplace_back( std::pair<size_t, size_t>(itrk_id, jtrk_id) );
738 
739  // two circles intersection
740  const Trk::Perigee& perigee1 = itrk->perigeeParameters();
741  const Trk::Perigee& perigee2 = jtrk->perigeeParameters();
742  int flag = 0;
743  int errorcode = 0;
744  InDet::VertexPointEstimator::Values_t decors;
745  timerTwoCircIntsect.start();
746  Amg::Vector3D circIntersect = m_vertexPointEstimator->getCirclesIntersectionPoint(&perigee1, &perigee2, flag, errorcode, decors);
747  timerTwoCircIntsect.stop();
748  if (errorcode != 0) {
749  ATH_MSG_VERBOSE( ": two circles intersect failed");
750  wrkcuts.twoCircErrcut = true;
751  wrkvrt.param = wrkprm;
752  wrkvrt.cuts = wrkcuts;
753  wrkvrt.isPair = true;
754 
755  workVerticesContainer.push_back( std::move(wrkvrt) );
756  continue;
757  }
758  float dphi = std::abs(decors["deltaPhiTracks"]);
759  float dr = std::abs(decors["DR1R2"]);
760  float intersect_r = circIntersect.perp();
761  mnt_pair_dphi.push_back(dphi);
762  mnt_pair_dr.push_back(dr);
763  mnt_intersect_r.push_back(intersect_r);
764 
765  wrkprm.twoCirc_dphi = dphi;
766  wrkprm.twoCirc_dr = dr;
767  wrkprm.twoCirc_int_r = intersect_r;
768 
769  nPairsIntersect++;
770  wrkcuts.twoCircRcut = intersect_r > m_maxR;
771  nPairsIntersectPos++;
772 
773  // initial fast fitting
774  std::vector<const xAOD::TrackParticle*> baseTracks;
775  baseTracks.emplace_back( itrk );
776  baseTracks.emplace_back( jtrk );
777  Amg::Vector3D initVertex;
778  timerVrtFitFast.start();
779  auto fitterState = m_fitSvc->makeState(ctx);
780  m_fitSvc->setApproximateVertex( 0., 0., 0., *fitterState );
781  StatusCode sc = m_fitSvc->VKalVrtFitFast( baseTracks, initVertex, *fitterState );
782  timerVrtFitFast.stop();
783  if( sc.isFailure() ) {
784  ATH_MSG_VERBOSE( ": fast crude estimation fails ");
785  wrkcuts.fastErrcut = true;
786  wrkvrt.param = wrkprm;
787  wrkvrt.cuts = wrkcuts;
788  wrkvrt.isPair = true;
789 
790  workVerticesContainer.push_back( std::move(wrkvrt) );
791  continue;
792  }
793 
794  // position of initial fast fitting
795  mnt_init_r.push_back(initVertex.perp());
796  wrkprm.vrtFast_r = initVertex.perp();
797  wrkprm.vrtFast_eta = initVertex.eta();
798  wrkprm.vrtFast_phi = initVertex.phi();
799  wrkcuts.fastRcut = initVertex.perp() > m_maxR;
800 
801  nPairsInitPos++;
802 
803  // impact parameter on initial fast fitting
804  std::vector<double> impactParameters;
805  std::vector<double> impactParErrors;
806  if( ! m_fitSvc->VKalGetImpact(itrk, initVertex, static_cast<int>( itrk->charge() ), impactParameters, impactParErrors) ) continue;
807  const auto roughD0_itrk = impactParameters.at(TrkParameter::k_d0);
808  const auto roughZ0_itrk = impactParameters.at(TrkParameter::k_z0);
809  mnt_init_trkd0.push_back(std::abs(roughD0_itrk));
810  mnt_init_trkz0.push_back(std::abs(roughZ0_itrk));
811  wrkprm.vrtFast_trkd0.push_back(std::abs(roughD0_itrk));
812  wrkprm.vrtFast_trkz0.push_back(std::abs(roughZ0_itrk));
813  if( std::abs(roughD0_itrk) > roughD0Cut || std::abs(roughZ0_itrk) > roughZ0Cut ) continue;
814 
815  if( ! m_fitSvc->VKalGetImpact(jtrk, initVertex, static_cast<int>( jtrk->charge() ), impactParameters, impactParErrors) ) continue;
816  const auto roughD0_jtrk = impactParameters.at(TrkParameter::k_d0);
817  const auto roughZ0_jtrk = impactParameters.at(TrkParameter::k_z0);
818  mnt_init_trkd0.push_back(std::abs(roughD0_jtrk));
819  mnt_init_trkz0.push_back(std::abs(roughZ0_jtrk));
820  wrkprm.vrtFast_trkd0.push_back(std::abs(roughD0_jtrk));
821  wrkprm.vrtFast_trkz0.push_back(std::abs(roughZ0_jtrk));
822  if( std::abs(roughD0_jtrk) > roughD0Cut || std::abs(roughZ0_jtrk) > roughZ0Cut ) continue;
823 
824  nPairsInitTrkd0z0++;
825 
826  // Vertex VKal Fitting
827  std::vector<const xAOD::NeutralParticle*> dummyNeutrals;
828  timerVrtFit.start();
829  m_fitSvc->setApproximateVertex( initVertex.x(), initVertex.y(), initVertex.z(), *fitterState );
830  sc = m_fitSvc->VKalVrtFit( baseTracks,
831  dummyNeutrals,
832  wrkvrt.vertex, wrkvrt.vertexMom, wrkvrt.charge,
833  wrkvrt.vertexCov, wrkvrt.chi2PerTrk,
834  wrkvrt.trkAtVrt, wrkvrt.chi2, *fitterState, false );
835  timerVrtFit.stop();
836  if( sc.isFailure() ) {
837  wrkcuts.fitErrcut = true;
838  wrkvrt.param = wrkprm;
839  wrkvrt.cuts = wrkcuts;
840  wrkvrt.isPair = true;
841 
842  workVerticesContainer.push_back( std::move(wrkvrt) );
843  continue;
844  }
845 
846  nPairsVtxFit++;
847 
848  // Chi2 cut
849  mnt_vtxfit_chi2.push_back(wrkvrt.chi2);
850  if( wrkvrt.chi2 > m_selVrtChi2Cut) {
851  ATH_MSG_VERBOSE( ": failed to pass chi2 threshold." );
852  wrkcuts.chi2cut = true;
853  wrkvrt.param = wrkprm;
854  wrkvrt.cuts = wrkcuts;
855  wrkvrt.isPair = true;
856 
857  workVerticesContainer.push_back( std::move(wrkvrt) );
858  continue;
859  }
860  nPairsVtxChi2++;
861  mnt_vtxfit_r.push_back(wrkvrt.vertex.perp());
862 
863  // Compatibility to the primary vertex
864  Amg::Vector3D vDist = (primVertex!=nullptr)? wrkvrt.vertex - primVertex->position() : wrkvrt.vertex;
865  const double vPos = ( vDist.x()*wrkvrt.vertexMom.Px()+vDist.y()*wrkvrt.vertexMom.Py()+vDist.z()*wrkvrt.vertexMom.Pz() )/wrkvrt.vertexMom.Rho();
866  const double vPosMomAngT = ( vDist.x()*wrkvrt.vertexMom.Px()+vDist.y()*wrkvrt.vertexMom.Py() ) / vDist.perp() / wrkvrt.vertexMom.Pt();
867 
868  double dphi1 = vDist.phi() - itrk->phi(); while( dphi1 > TMath::Pi() ) { dphi1 -= TMath::TwoPi(); } while( dphi1 < -TMath::Pi() ) { dphi1 += TMath::TwoPi(); }
869  double dphi2 = vDist.phi() - jtrk->phi(); while( dphi2 > TMath::Pi() ) { dphi2 -= TMath::TwoPi(); } while( dphi2 < -TMath::Pi() ) { dphi2 += TMath::TwoPi(); }
870 
871  wrkprm.vPos = vPos;
872  wrkprm.vPosMomAngT = vPosMomAngT;
873  wrkprm.dphi1 = dphi1;
874  wrkprm.dphi2 = dphi2;
875 
876  if( std::cos( dphi1 ) < m_dphiPVCut && std::cos( dphi2 ) < m_dphiPVCut ) {
877  ATH_MSG_DEBUG( ": failed to pass the vPos cut. (both tracks are opposite against the vertex pos)" );
878  continue;
879  }
880  if( vPosMomAngT < m_dphiPVCut ) {
881  ATH_MSG_DEBUG( ": failed to pass the vPos cut. (pos-mom directions are opposite)" );
882  continue;
883  }
884  if( vPos < m_pvCompatibilityCut ) {
885  ATH_MSG_DEBUG( ": failed to pass the vPos cut." );
886  continue;
887  }
888  nPairsVtxComp++;
889 
890  // Now this vertex passed all criteria and considred to be a compatible vertices.
891  // Therefore the track pair is removed from the incompatibility list.
892  if (wrkcuts.isFullPass()) incomp.pop_back();
893 
894  wrkvrt.param = wrkprm;
895  wrkvrt.cuts = wrkcuts;
896  wrkvrt.isGood = true;
897  wrkvrt.isPair = true;
898 
899  workVerticesContainer.push_back( std::move(wrkvrt) );
900  }
901  }
902  ATH_MSG_DEBUG("Of " << nPairsAll << " pairs : trkd0" << " / " << "intersect" << " / " << "intersectPos" << " / " << "initPos" << " / " << "initD0Z0" << " / " << "vtxFit" << " / " << "vtxChi2" << " / " << "vtxComp = "
903  << nPairsTrkd0 << " / " << nPairsIntersect << " / " << nPairsIntersectPos << " / " << nPairsInitPos << " / " << nPairsInitTrkd0z0 << " / " << nPairsVtxFit << " / " << nPairsVtxChi2 << " / " << nPairsVtxComp );
904  //
905  return StatusCode::SUCCESS;
906 }
907 
908 
909 StatusCode TrigVrtSecInclusive::findNtrackVerticesVSI
910 ( WrkVrtContainer& workVerticesContainer, std::vector<std::pair<size_t,size_t>>& incomp, std::vector<const xAOD::TrackParticle*>& selectedTracks, const EventContext& ctx ) const
911 {
912  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": begin");
913 
914  const auto compSize = selectedTracks.size()*(selectedTracks.size() - 1)/2 - incomp.size();
915 
916  auto pgraph = std::make_unique<Trk::PGraph>();
917 
918  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": compatible track pair size = " << compSize );
919  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": incompatible track pair size = " << incomp.size() );
920 
921  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": incompatibility graph finder mode" );
922 
923  // clear the container
924  workVerticesContainer.clear();
925 
926  // Graph method: Trk::pgraphm_()
927  // used in order to find compatible sub-graphs from the incompatible graph
928 
929  // List of edges between incompatible nodes
930  // This weight is the data model of incompatible graph used in Trk::pgraphm_().
931  std::vector<long int> weight;
932 
933  for( auto& pair : incomp ) {
934  weight.emplace_back( pair.first + 1 ); /* +1 is needed for PGRAPH due to FORTRAN-style counting */
935  weight.emplace_back( pair.second + 1 ); /* +1 is needed for PGRAPH due to FORTRAN-style counting */
936  }
937 
938  // Solution of the graph method routine (minimal covering of the graph)
939  // The size of the solution is returned by NPTR (see below)
940  std::vector<long int> solution( selectedTracks.size() );
941 
942  // Number of edges in the list is the size of incompatibility track pairs.
943  long int nEdges = incomp.size();
944 
945  // input number of nodes in the graph.
946  long int nTracks = static_cast<long int>( selectedTracks.size() );
947 
948  // Input variable: the threshold. Solutions shorter than nth are not returned (ignored).
949  long int nth = 2; //VK some speed up
950 
951  // NPTR: I/O variable (Destructive FORTRAN Style!!!)
952  // - on input: =0 for initialization, >0 to get next solution
953  // - on output: >0 : length of the solution stored in set; =0 : no more solutions can be found
954  long int solutionSize { 0 };
955 
956  // This is just a unused strawman needed for m_fitSvc->VKalVrtFit()
957  std::vector<const xAOD::TrackParticle*> baseTracks;
958  std::vector<const xAOD::NeutralParticle*> dummyNeutrals;
959 
960  // Main iteration
961  while(true) {
962  // Find a solution from the given set of incompatible tracks (==weight)
963  pgraph->pgraphm_( weight.data(), nEdges, nTracks, solution.data(), &solutionSize, nth);
964 
965  ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": Trk::pgraphm_() output: solutionSize = " << solutionSize );
966 
967  if(solutionSize <= 0) break; // No more solutions ==> Exit
968  if(solutionSize == 1) continue; // i.e. single node ==> Not a good solution
969 
970  baseTracks.clear();
971 
972  std::string msg = "solution = [ ";
973  for( int i=0; i< solutionSize; i++) {
974  msg += Form( "%ld, ", solution[i]-1 );
975  }
976  msg += " ]";
977  ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": " << msg );
978 
979  // varaible of new vertex
980  WrkVrt wrkvrt;
981 
982  // Try to compose a new vertex using the solution nodes
983  // Here the track ID is labelled with array
984  wrkvrt.isGood = true;
985  wrkvrt.selectedTrackIndices().clear();
986 
987  for(long int i = 0; i<solutionSize; i++) {
988  wrkvrt.selectedTrackIndices().emplace_back(solution[i]-1);
989  baseTracks.emplace_back( selectedTracks.at(solution[i]-1) );
990  }
991 
992  // Perform vertex fitting
993  Amg::Vector3D initVertex;
994  auto fitterState = m_fitSvc->makeState(ctx);
995 
996  StatusCode sc = m_fitSvc->VKalVrtFitFast( baseTracks, initVertex, *fitterState );/* Fast crude estimation */
997  if(sc.isFailure()) ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": fast crude estimation fails ");
998 
999  m_fitSvc->setApproximateVertex( initVertex.x(), initVertex.y(), initVertex.z(), *fitterState );
1000 
1001  sc = m_fitSvc->VKalVrtFit(baseTracks, dummyNeutrals,
1002  wrkvrt.vertex,
1003  wrkvrt.vertexMom,
1004  wrkvrt.charge,
1005  wrkvrt.vertexCov,
1006  wrkvrt.chi2PerTrk,
1007  wrkvrt.trkAtVrt,
1008  wrkvrt.chi2,
1009  *fitterState, false);
1010 
1011  ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": FoundAppVrt=" << solutionSize << ", (r, z) = " << wrkvrt.vertex.perp() << ", " << wrkvrt.vertex.z() << ", chi2/ndof = " << wrkvrt.fitQuality() );
1012 
1013  if( sc.isFailure() ) {
1014 
1015  if( wrkvrt.selectedTrackIndices().size() <= 2 ) {
1016  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": VKalVrtFit failed in 2-trk solution ==> give up.");
1017  continue;
1018  }
1019 
1020  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": VKalVrtFit failed ==> retry...");
1021 
1022  WrkVrt tmp;
1023  tmp.isGood = false;
1024 
1025  // Create 2-trk vertex combination and find any compatible vertex
1026  for( auto& itrk: wrkvrt.selectedTrackIndices() ) {
1027  for( auto& jtrk: wrkvrt.selectedTrackIndices() ) {
1028  if( itrk == jtrk ) continue;
1029  if( tmp.isGood ) break;
1030 
1031  tmp.selectedTrackIndices().clear();
1032  tmp.selectedTrackIndices().emplace_back( itrk );
1033  tmp.selectedTrackIndices().emplace_back( jtrk );
1034 
1035  baseTracks.clear();
1036  baseTracks.emplace_back( selectedTracks.at( itrk ) );
1037  baseTracks.emplace_back( selectedTracks.at( jtrk ) );
1038 
1039  // Perform vertex fitting
1040  Amg::Vector3D initVertex;
1041 
1042  sc = m_fitSvc->VKalVrtFitFast( baseTracks, initVertex, *fitterState );
1043  if( sc.isFailure() ) ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": fast crude estimation fails ");
1044 
1045  m_fitSvc->setApproximateVertex( initVertex.x(), initVertex.y(), initVertex.z(), *fitterState );
1046 
1047  sc = m_fitSvc->VKalVrtFit(baseTracks, dummyNeutrals,
1048  tmp.vertex,
1049  tmp.vertexMom,
1050  tmp.charge,
1051  tmp.vertexCov,
1052  tmp.chi2PerTrk,
1053  tmp.trkAtVrt,
1054  tmp.chi2,
1055  *fitterState, false);
1056 
1057  if( sc.isFailure() ) continue;
1058 
1059  tmp.isGood = true;
1060 
1061  }
1062  }
1063 
1064  if( !tmp.isGood ) {
1065  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Did not find any viable vertex in all 2-trk combinations. Give up.");
1066  continue;
1067  }
1068 
1069  // Now, found at least one seed 2-track vertex. ==> attempt to attach other tracks
1070  for( auto& itrk: wrkvrt.selectedTrackIndices() ) {
1071 
1072  if( std::find( tmp.selectedTrackIndices().begin(), tmp.selectedTrackIndices().end(), itrk ) != tmp.selectedTrackIndices().end() ) continue;
1073 
1074  auto backup = tmp;
1075 
1076  tmp.selectedTrackIndices().emplace_back( itrk );
1077  baseTracks.clear();
1078  for( auto& jtrk : tmp.selectedTrackIndices() ) { baseTracks.emplace_back( selectedTracks.at(jtrk) ); }
1079 
1080  // Perform vertex fitting
1081  Amg::Vector3D initVertex;
1082 
1083  sc = m_fitSvc->VKalVrtFitFast( baseTracks, initVertex, *fitterState );/* Fast crude estimation */
1084  if(sc.isFailure()) ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": fast crude estimation fails ");
1085 
1086  m_fitSvc->setApproximateVertex( initVertex.x(), initVertex.y(), initVertex.z(), *fitterState );
1087  sc = m_fitSvc->VKalVrtFit(baseTracks, dummyNeutrals,
1088  tmp.vertex,
1089  tmp.vertexMom,
1090  tmp.charge,
1091  tmp.vertexCov,
1092  tmp.chi2PerTrk,
1093  tmp.trkAtVrt,
1094  tmp.chi2,
1095  *fitterState, false);
1096 
1097  if( sc.isFailure() ) {
1098  tmp = backup;
1099  continue;
1100  }
1101 
1102  }
1103 
1104  wrkvrt = tmp;
1105  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": VKalVrtFit succeeded; register the vertex to the list.");
1106  wrkvrt.isGood = true;
1107  wrkvrt.closestWrkVrtIndex = std::numeric_limits<unsigned>::max();
1108  wrkvrt.closestWrkVrtValue = std::numeric_limits<double>::max();
1109  workVerticesContainer.emplace_back( wrkvrt );
1110 
1111  } else {
1112 
1113  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": VKalVrtFit succeeded; register the vertex to the list.");
1114  wrkvrt.isGood = true;
1115  wrkvrt.closestWrkVrtIndex = std::numeric_limits<unsigned>::max();
1116  wrkvrt.closestWrkVrtValue = std::numeric_limits<double>::max();
1117  workVerticesContainer.emplace_back( wrkvrt );
1118 
1119  }
1120 
1121  }
1122 
1123 
1124 
1125 
1126  if (m_truncateWrkVertices){
1127  if (workVerticesContainer.size() > m_maxWrkVertices){
1128  ATH_MSG_INFO("WrkVertex truncated. Too many vertices");
1129  workVerticesContainer.resize(m_maxWrkVertices);
1130  }
1131  }
1132 
1133  ATH_CHECK( cleanUp(workVerticesContainer) );
1134 
1135  //-Identify remaining 2-track vertices with very bad Chi2 and mass (b-tagging)
1136  ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": Identify remaining 2-track vertices with very bad Chi2 and mass (b-tagging).");
1137  for( auto& wrkvrt : workVerticesContainer ) {
1138 
1139  if( TMath::Prob( wrkvrt.chi2, wrkvrt.ndof() ) < m_improveChi2ProbThreshold ) wrkvrt.isGood = false;
1140  if( wrkvrt.selectedTrackIndices().size() != 2 ) continue;
1141  }
1142 
1143  return StatusCode::SUCCESS;
1144 
1145 }
1146 
1147 
1148 
1149 template<typename VrtType, typename Coord>
1150 StatusCode TrigVrtSecInclusive::findNTrackVertex
1151 ( TrigVrtSecInclusive::WrkVrtContainer& workVerticesContainer, TrigVSI::VtxMap<VrtType,Coord>& vtxmap, const std::vector<const xAOD::TrackParticle*>& selectedTracks, const EventContext& ctx ) const
1152 {
1153  ATH_MSG_DEBUG( "findNTrackVertex start" );
1154 
1155  // monitoring
1156  auto timerRetrvFromMap = Monitored::Timer<std::chrono::nanoseconds> ("TIME_RetrvFromMap");
1157  auto timerMergeParGraph = Monitored::Timer<std::chrono::microseconds> ("TIME_MergeParGraph");
1158  auto timerMergeSimple = Monitored::Timer<std::chrono::microseconds> ("TIME_MergeSimple");
1159 
1160  size_t n_clst = vtxmap.nClusters();
1161  size_t n_vtx_pargraph = 0;
1162  size_t n_vtx_simple = 0;
1163 
1164  for (size_t i_clst = 0; i_clst < n_clst; i_clst++) {
1165  auto monTimeClust = Monitored::Group(m_monTool, timerRetrvFromMap, timerMergeParGraph, timerMergeSimple);
1166  timerRetrvFromMap.start();
1167 
1168  ATH_MSG_DEBUG( "Retrieve cluster, track indices and incompatible pairs" );
1169 
1170  auto clst = vtxmap.getCluster(i_clst);
1171  auto selected_track_indices = clst.getSelectedTrackIndices();
1172  auto incomp = clst.getIncompIndices();
1173 
1174  ATH_MSG_DEBUG( "Convert set to vector" );
1175 
1176  std::vector<size_t> v_track_indices( selected_track_indices.begin(), selected_track_indices.end() );
1177 
1178  timerRetrvFromMap.stop();
1179 
1180  ATH_MSG_DEBUG( "Cluster number : " << i_clst << " | " << incomp.size() << " incompatible pairs | " << v_track_indices.size() << " tracks" );
1181  ATH_MSG_DEBUG( "Pos avr : ( R: " << clst.PosCoord().at(0) << ", eta: " << clst.PosCoord().at(1) << ", phi: " << clst.PosCoord().at(2) << " )" );
1182  ATH_MSG_DEBUG( "Number of cells : " << clst.nPoints() );
1183 
1184  if ( v_track_indices.size() > m_maxTrks ) {
1185  ATH_MSG_DEBUG( "Skipped the cluster. Too many tracks" );
1186  continue;
1187  }
1188 
1189  if ( clst.nVtx() >= m_minTrkPairsMerge || selected_track_indices.size() >= m_minTrksMerge ) {
1190 
1191  timerMergeParGraph.start();
1192  ATH_CHECK( mergeVertexFromDiTrkVrt( workVerticesContainer, incomp, v_track_indices, selectedTracks, ctx ) );
1193  timerMergeParGraph.stop();
1194 
1195  n_vtx_pargraph++;
1196 
1197  } else {
1198 
1199  timerMergeSimple.start();
1200 
1201  WrkVrt wrkvrt;
1202  wrkvrt.isGood = true;
1203  wrkvrt.selectedTrackIndices().clear();
1204 
1205  for (size_t i_trk = 0; i_trk < v_track_indices.size(); i_trk++) {
1206  wrkvrt.selectedTrackIndices().emplace_back( v_track_indices.at(i_trk) );
1207  }
1208 
1209  if ( fitVertexFromTracks(wrkvrt, selectedTracks, ctx) == StatusCode::SUCCESS ) {
1210  workVerticesContainer.emplace_back( std::move(wrkvrt) );
1211  }
1212 
1213  timerMergeSimple.stop();
1214 
1215  n_vtx_simple++;
1216 
1217  }
1218 
1219  }
1220 
1221  ATH_MSG_DEBUG( "Partial graph method / Simple method / All vertex = " << n_vtx_pargraph << " / " << n_vtx_simple << " / " << n_vtx_pargraph+n_vtx_simple );
1222  return StatusCode::SUCCESS;
1223 
1224 }
1225 
1226 template StatusCode TrigVrtSecInclusive::findNTrackVertex ( WrkVrtContainer&, TrigVSI::VtxMap<TrigVrtSecInclusive::WrkVrt,TrigVSI::Coordinate::Pseudo>&, const std::vector<const xAOD::TrackParticle*>&, const EventContext& ) const;
1227 
1228 
1229 
1243 StatusCode TrigVrtSecInclusive::mergeVertexFromDiTrkVrt
1244 ( TrigVrtSecInclusive::WrkVrtContainer& workVerticesContainer, const std::vector<std::pair<size_t,size_t>>& incomp, const std::vector<size_t>& trkIdx, const std::vector<const xAOD::TrackParticle*>& selectedTracks, const EventContext& ctx ) const
1245 {
1246 
1247  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": begin");
1248  ATH_MSG_DEBUG("TrigVrtSecInclusive::mergeVertexFromDiTrkVrt : start");
1249 
1250  // Graph method: Trk::pgraphm_()
1251  // used in order to find compatible sub-graphs from the incompatible graph
1252  auto pgraph = std::make_unique<Trk::PGraph>();
1253 
1254  std::unordered_map<size_t,size_t> dict_trk_idx;
1255  size_t n_trk = trkIdx.size(); // Dictionary: global trkId -> local trkId
1256 
1257  for (size_t i = 0; i < n_trk; i++) {
1258  dict_trk_idx.emplace( trkIdx.at(i), i );
1259  }
1260 
1261  const auto compSize = (n_trk*(n_trk - 1))/2 - incomp.size();
1262 
1263  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": compatible track pair size = " << compSize );
1264  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": incompatible track pair size = " << incomp.size() );
1265 
1266  // List of edges between incompatible nodes
1267  // This weight is the data model of incompatible graph used in Trk::pgraphm_().
1268  std::vector<long int> weight;
1269 
1270  for( auto& pair : incomp ) {
1271  weight.emplace_back( dict_trk_idx[pair.first] + 1 ); /* +1 is needed for PGRAPH due to FORTRAN-style counting */
1272  weight.emplace_back( dict_trk_idx[pair.second] + 1 ); /* +1 is needed for PGRAPH due to FORTRAN-style counting */
1273  }
1274 
1275  // Solution of the graph method routine (minimal covering of the graph)
1276  // The size of the solution is returned by NPTR (see below)
1277  std::vector<long int> solution( n_trk );
1278 
1279  // Number of edges in the list is the size of incompatibility track pairs.
1280  long int nEdges = incomp.size();
1281 
1282  // input number of nodes in the graph.
1283  long int nTracks = static_cast<long int>( n_trk );
1284 
1285  // Input variable; the threshold. Solutions shorter than nth are not returned (ignored).
1286  long int nth = 2; //VK some speed up
1287 
1288  // NPTR: I/O variable (Destructive FORTRAN Style!!!)
1289  // - on input: =0 for initialization, >0 to get next solution
1290  // - on output: >0 : length of the solution stored in set; =0 : no more solutions can be found
1291  long int solutionSize { 0 };
1292 
1293  ATH_MSG_DEBUG( "TrigVrtSecInclusive::mergeVertexFromDiTrkVrt : while loop start" );
1294 
1295  // Main iteration
1296  while(true) {
1297 
1298  // Find a solution from the given set of incompatible tracks (==weight)
1299  pgraph->pgraphm_( weight.data(), nEdges, nTracks, solution.data(), &solutionSize, nth);
1300 
1301  ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": Trk::pgraphm_() output: solutionSize = " << solutionSize );
1302  ATH_MSG_DEBUG( "TrigVrtSecInclusive::mergeVertexFromDiTrkVrt : Trk::pgraphm_() output: solutionSize = " << solutionSize );
1303 
1304  if(solutionSize <= 0) break; // No more solutions ==> Exit
1305  if(solutionSize == 1) continue; // i.e. single node ==> Not a good solution
1306 
1307  // Print solution
1308  std::string msg = "solution = [ ";
1309  for( int i=0; i< solutionSize; i++) {
1310  msg += Form( "%ld, ", trkIdx[ solution[i]-1 ] );
1311  }
1312  msg += " ]";
1313  ATH_MSG_DEBUG( " > " << __FUNCTION__ << ": " << msg );
1314 
1315  // varaible of new vertex
1316  WrkVrt wrkvrt;
1317 
1318  // Try to compose a new vertex using the solution nodes
1319  // Here the track ID is labelled with array
1320  wrkvrt.isGood = true;
1321  wrkvrt.selectedTrackIndices().clear();
1322 
1323  for(long int i = 0; i<solutionSize; i++) {
1324  wrkvrt.selectedTrackIndices().emplace_back( trkIdx[ solution[i]-1 ] );
1325  }
1326 
1327  if ( fitVertexFromTracks(wrkvrt, selectedTracks, ctx) == StatusCode::SUCCESS ) {
1328  workVerticesContainer.emplace_back( std::move(wrkvrt) );
1329  }
1330 
1331  }
1332 
1333  return StatusCode::SUCCESS;
1334 }
1335 
1336 
1350 StatusCode TrigVrtSecInclusive::fitVertexFromTracks
1351 ( WrkVrt& wrkvrt, const std::vector<const xAOD::TrackParticle*>& selectedTracks, const EventContext& ctx ) const
1352 {
1353 
1354  std::vector<const xAOD::TrackParticle*> baseTracks;
1355  std::vector<const xAOD::NeutralParticle*> dummyNeutrals; // This is just a unused strawman needed for m_fitSvc->VKalVrtFit()
1356 
1357  baseTracks.clear();
1358 
1359  // Try to compose a new vertex using the solution nodes
1360  // Here the track ID is labelled with array
1361  wrkvrt.isGood = true;
1362 
1363  size_t n_trk = wrkvrt.selectedTrackIndices().size();
1364 
1365  for(size_t i = 0; i<n_trk; i++) {
1366  baseTracks.emplace_back( selectedTracks.at( wrkvrt.selectedTrackIndices().at(i) ) );
1367  }
1368 
1369  // Perform vertex fitting
1370  Amg::Vector3D initVertex;
1371  auto fitterState = m_fitSvc->makeState(ctx);
1372 
1373  StatusCode sc = m_fitSvc->VKalVrtFitFast( baseTracks, initVertex, *fitterState );/* Fast crude estimation */
1374  if(sc.isFailure()) ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": fast crude estimation fails ");
1375 
1376  m_fitSvc->setApproximateVertex( initVertex.x(), initVertex.y(), initVertex.z(), *fitterState );
1377 
1378  sc = m_fitSvc->VKalVrtFit(baseTracks, dummyNeutrals,
1379  wrkvrt.vertex,
1380  wrkvrt.vertexMom,
1381  wrkvrt.charge,
1382  wrkvrt.vertexCov,
1383  wrkvrt.chi2PerTrk,
1384  wrkvrt.trkAtVrt,
1385  wrkvrt.chi2,
1386  *fitterState, false);
1387 
1388  ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": FoundAppVrt=" << n_trk << ", (r, z) = " << wrkvrt.vertex.perp() << ", " << wrkvrt.vertex.z() << ", chi2/ndof = " << wrkvrt.fitQuality() );
1389 
1390  if( sc.isFailure() ) {
1391 
1392  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": VKalVrtFit failed ==> retry...");
1393 
1394  WrkVrt tmp;
1395  tmp.isGood = false;
1396 
1397  // Create 2-trk vertex combination and find any compatible vertex
1398  for( auto& itrk: wrkvrt.selectedTrackIndices() ) {
1399  for( auto& jtrk: wrkvrt.selectedTrackIndices() ) {
1400  if( itrk == jtrk ) continue;
1401  if( tmp.isGood ) break;
1402 
1403  tmp.selectedTrackIndices().clear();
1404  tmp.selectedTrackIndices().emplace_back( itrk );
1405  tmp.selectedTrackIndices().emplace_back( jtrk );
1406 
1407  baseTracks.clear();
1408  baseTracks.emplace_back( selectedTracks.at( itrk ) );
1409  baseTracks.emplace_back( selectedTracks.at( jtrk ) );
1410 
1411  // Perform vertex fitting
1412  Amg::Vector3D initVertex;
1413 
1414  sc = m_fitSvc->VKalVrtFitFast( baseTracks, initVertex, *fitterState );
1415  if( sc.isFailure() ) ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": fast crude estimation fails ");
1416 
1417  m_fitSvc->setApproximateVertex( initVertex.x(), initVertex.y(), initVertex.z(), *fitterState );
1418 
1419  sc = m_fitSvc->VKalVrtFit(baseTracks, dummyNeutrals,
1420  tmp.vertex,
1421  tmp.vertexMom,
1422  tmp.charge,
1423  tmp.vertexCov,
1424  tmp.chi2PerTrk,
1425  tmp.trkAtVrt,
1426  tmp.chi2,
1427  *fitterState, false);
1428 
1429  if( sc.isFailure() ) continue;
1430  if( tmp.chi2 > m_selVrtChi2Cut ) continue;
1431 
1432  tmp.isGood = true;
1433 
1434  }
1435  }
1436 
1437  if( !tmp.isGood ) {
1438  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": Did not find any viable vertex in all 2-trk combinations. Give up.");
1439  wrkvrt = std::move(tmp);
1440  return StatusCode::FAILURE;
1441  }
1442 
1443  // Now, found at least one seed 2-track vertex. ==> attempt to attach other tracks
1444  for( auto& itrk: wrkvrt.selectedTrackIndices() ) {
1445 
1446  if( std::find( tmp.selectedTrackIndices().begin(), tmp.selectedTrackIndices().end(), itrk ) != tmp.selectedTrackIndices().end() ) continue;
1447 
1448  auto backup = tmp;
1449 
1450  tmp.selectedTrackIndices().emplace_back( itrk );
1451  baseTracks.clear();
1452  for( auto& jtrk : tmp.selectedTrackIndices() ) { baseTracks.emplace_back( selectedTracks.at(jtrk) ); }
1453 
1454  // Perform vertex fitting
1455  Amg::Vector3D initVertex;
1456 
1457  sc = m_fitSvc->VKalVrtFitFast( baseTracks, initVertex, *fitterState );/* Fast crude estimation */
1458  if(sc.isFailure()) ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": fast crude estimation fails ");
1459 
1460  m_fitSvc->setApproximateVertex( initVertex.x(), initVertex.y(), initVertex.z(), *fitterState );
1461  sc = m_fitSvc->VKalVrtFit(baseTracks, dummyNeutrals,
1462  tmp.vertex,
1463  tmp.vertexMom,
1464  tmp.charge,
1465  tmp.vertexCov,
1466  tmp.chi2PerTrk,
1467  tmp.trkAtVrt,
1468  tmp.chi2,
1469  *fitterState, false);
1470 
1471  if( sc.isFailure() ) {
1472  tmp = backup;
1473  continue;
1474  }
1475 
1476  }
1477 
1478  wrkvrt = std::move(tmp);
1479 
1480  }
1481 
1482  ATH_MSG_DEBUG(" > " << __FUNCTION__ << ": VKalVrtFit succeeded; register the vertex to the list.");
1483  wrkvrt.isGood = true;
1484  wrkvrt.closestWrkVrtIndex = std::numeric_limits<unsigned>::max();
1485  wrkvrt.closestWrkVrtValue = std::numeric_limits<double>::max();
1486 
1487  return StatusCode::SUCCESS;
1488 }
1489 
1490 
1491 bool TrigVrtSecInclusive::selectTrack( const xAOD::TrackParticle* trk ) const
1492 {
1493  if( ! selectTrack_pTCut(trk) ) return false;
1494  if( ! selectTrack_d0Cut(trk) ) return false;
1495  if( ! selectTrack_z0Cut(trk) ) return false;
1496  if( ! selectTrack_chi2Cut(trk) ) return false;
1497  if( ! selectTrack_hitPattern(trk) ) return false;
1498  return true;
1499 }
1500 
1501 
1502 bool TrigVrtSecInclusive::selectTrack_d0Cut ( const xAOD::TrackParticle* trk ) const { return ( std::abs( trk->d0() ) > m_d0TrkPVDstMinCut && std::abs( trk->d0() ) < m_d0TrkPVDstMaxCut ); }
1503 bool TrigVrtSecInclusive::selectTrack_z0Cut ( const xAOD::TrackParticle* trk ) const { return ( std::abs( trk->z0() ) > m_z0TrkPVDstMinCut && std::abs( trk->z0() ) < m_z0TrkPVDstMaxCut ); }
1504 bool TrigVrtSecInclusive::selectTrack_pTCut ( const xAOD::TrackParticle* trk ) const { return trk->pt() > m_trkPtCut; }
1505 bool TrigVrtSecInclusive::selectTrack_chi2Cut ( const xAOD::TrackParticle* trk ) const { return trk->chiSquared() / (trk->numberDoF()+1e-9) < m_trkChi2Cut; }
1506 
1507 
1508 bool TrigVrtSecInclusive::selectTrack_hitPattern( const xAOD::TrackParticle* trk ) const
1509 {
1510  uint8_t PixelHits = 0;
1511  uint8_t SCTHits = 0;
1512  uint8_t BLayHits = 0;
1513  uint8_t PixShare = 0;
1514  uint8_t SCTShare = 0;
1515  uint8_t TRTHits = 0;
1516 
1517  if( !(trk->summaryValue( PixelHits, xAOD::numberOfPixelHits ) ) ) PixelHits =0;
1518  if( !(trk->summaryValue( SCTHits, xAOD::numberOfSCTHits ) ) ) SCTHits =0;
1519  if( !(trk->summaryValue( BLayHits, xAOD::numberOfInnermostPixelLayerHits ) ) ) BLayHits =0;
1520  if( !(trk->summaryValue( PixShare, xAOD::numberOfPixelSharedHits ) ) ) PixShare =0;
1521  if( !(trk->summaryValue( SCTShare, xAOD::numberOfSCTSharedHits ) ) ) SCTShare =0;
1522  if( !(trk->summaryValue( TRTHits, xAOD::numberOfTRTHits ) ) ) TRTHits =0;
1523 
1524  uint8_t SharedHits = PixShare + SCTShare;
1525 
1526  // do Pixel/SCT/SiHits only if we exclude StandAlone TRT hits
1527  if(PixelHits < m_cutPixelHits) return false;
1528  if(SCTHits < m_cutSctHits) return false;
1529  if((PixelHits+SCTHits) < m_cutSiHits) return false;
1530  if(BLayHits < m_cutBLayHits) return false;
1531  if(SharedHits > m_cutSharedHits) return false;
1532 
1533  // passed
1534  return true;
1535 }
1536 
1537 
1538 size_t TrigVrtSecInclusive::nTrkCommon( WrkVrtContainer& workVerticesContainer, const std::pair<unsigned, unsigned>& pairIndex) const
1539 {
1540  //
1541  // Number of common tracks for 2 vertices
1542  //
1543 
1544  auto& trackIndices1 = workVerticesContainer.at( pairIndex.first ).selectedTrackIndices();
1545  auto& trackIndices2 = workVerticesContainer.at( pairIndex.second ).selectedTrackIndices();
1546 
1547  if( trackIndices1.size() < 2 ) return 0;
1548  if( trackIndices2.size() < 2 ) return 0;
1549 
1550  size_t nTrkCom = 0;
1551 
1552  for( auto& index : trackIndices1 ) {
1553  if( std::find(trackIndices2.begin(),trackIndices2.end(), index) != trackIndices2.end()) nTrkCom++;
1554  }
1555 
1556  return nTrkCom;
1557 }
1558 
1559 
1560 StatusCode TrigVrtSecInclusive::cleanUp
1561 ( WrkVrtContainer& workVerticesContainer ) const
1562 {
1563  //-Remove vertices fully contained in other vertices
1564  ATH_MSG_VERBOSE(" > " << __FUNCTION__ << ": Remove vertices fully contained in other vertices .");
1565  size_t n_vtx = workVerticesContainer.size();
1566 
1567  for (size_t iv = 0; iv < n_vtx; iv++) {
1568  for (size_t jv = iv+1; jv < n_vtx; jv++) {
1569 
1570  if ( !workVerticesContainer.at(iv).isGood ) continue;
1571  if ( !workVerticesContainer.at(jv).isGood ) continue;
1572 
1573  const auto nTCom = nTrkCommon( workVerticesContainer, {iv, jv} );
1574 
1575  if( nTCom == workVerticesContainer.at(iv).selectedTrackIndices().size() ) { workVerticesContainer.at(iv).isGood = false; continue; }
1576  else if( nTCom == workVerticesContainer.at(jv).selectedTrackIndices().size() ) { workVerticesContainer.at(jv).isGood = false; continue; }
1577 
1578  }
1579  }
1580 
1581  return StatusCode::SUCCESS;
1582 }
1583 
1584 
1585 StatusCode TrigVrtSecInclusive::finalize()
1586 {
1587  return StatusCode::SUCCESS;
1588 }
1589 
1590 TrigVrtSecInclusive::~TrigVrtSecInclusive() {}
1591 
1592 } // end of namespace TrigVSI
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgParam::isPassMMV
bool isPassMMV
Definition: TrigVrtSecInclusive.h:135
xAOD::TrackParticle_v1::pt
virtual double pt() const override final
The transverse momentum ( ) of the particle.
Definition: TrackParticle_v1.cxx:73
LinkToTrack.h
TrigVSI::TrigVrtSecInclusive::m_doMaterialMapVeto
Gaudi::Property< bool > m_doMaterialMapVeto
Definition: TrigVrtSecInclusive.h:87
TrigVSI::TrigVrtSecInclusive::findDiTrackVertex
StatusCode findDiTrackVertex(WrkVrtContainer &, std::vector< std::pair< size_t, size_t >> &, std::vector< const xAOD::TrackParticle * > &, const EventContext &, const xAOD::Vertex *) const
Definition: TrigVrtSecInclusive.cxx:379
DataModel_detail::const_iterator
Const iterator class for DataVector/DataList.
Definition: DVLIterator.h:82
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgParam::twoCirc_dr
double twoCirc_dr
Definition: TrigVrtSecInclusive.h:121
TrigVSI::TrigVrtSecInclusive::fillVtxContainer
StatusCode fillVtxContainer(xAODContainers &, const WrkVrtContainer &, std::vector< const xAOD::TrackParticle * > &) const
Definition: TrigVrtSecInclusive.cxx:288
TrigVSI::TrigVrtSecInclusive::TrigVrtSecInclusive
TrigVrtSecInclusive(const std::string &name, ISvcLocator *pSvcLocator)
Definition: TrigVrtSecInclusive.cxx:29
TrigDefs::Group
Group
Properties of a chain group.
Definition: GroupProperties.h:13
TrigVSI::TrigVrtSecInclusive::findDiTrackVertexVSI
StatusCode findDiTrackVertexVSI(WrkVrtContainer &, std::vector< std::pair< size_t, size_t >> &, std::vector< const xAOD::TrackParticle * > &, const EventContext &, const xAOD::Vertex *) const
Definition: TrigVrtSecInclusive.cxx:659
xAOD::Vertex_v1::nTrackParticles
size_t nTrackParticles() const
Get the number of tracks associated with this vertex.
Definition: Vertex_v1.cxx:270
TrigVSI::TrigVrtSecInclusive::m_materialMapInner
TH3S * m_materialMapInner
Definition: TrigVrtSecInclusive.h:177
TrigVSI::TrigVrtSecInclusive::WrkVrt::selectedTrackIndices
virtual std::deque< size_t > & selectedTrackIndices() override
Return indices of tracks associated with the vertex.
Definition: TrigVrtSecInclusive.h:153
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgParam::twoCirc_dphi
double twoCirc_dphi
Definition: TrigVrtSecInclusive.h:120
max
#define max(a, b)
Definition: cfImp.cxx:41
TrigVSI::TrigVrtSecInclusive::m_materialMapOuterFileName
std::string m_materialMapOuterFileName
Definition: TrigVrtSecInclusive.h:183
xAOD::numberOfSCTSharedHits
@ numberOfSCTSharedHits
number of SCT hits shared by several tracks [unit8_t].
Definition: TrackingPrimitives.h:272
xAOD::Vertex
Vertex_v1 Vertex
Define the latest version of the vertex class.
Definition: Event/xAOD/xAODTracking/xAODTracking/Vertex.h:16
TrigVSI::TrigVrtSecInclusive::WrkVrt::vertex
Amg::Vector3D vertex
flaged true for track pair vertex
Definition: TrigVrtSecInclusive.h:103
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:575
find
std::string find(const std::string &s)
return a remapped string
Definition: hcg.cxx:135
PathResolver::find_file
static std::string find_file(const std::string &logical_file_name, const std::string &search_path, SearchType search_type=LocalSearch)
Definition: PathResolver.cxx:251
SG::ReadHandle::cptr
const_pointer_type cptr()
Dereference the pointer.
TrigVSI::TrigVrtSecInclusive::fitVertexFromTracks
StatusCode fitVertexFromTracks(WrkVrt &, const std::vector< const xAOD::TrackParticle * > &, const EventContext &) const
Reconstruct vertex from given tracks.
Definition: TrigVrtSecInclusive.cxx:1351
SG::Accessor< float >
TrigVSI::TrigVrtSecInclusive::m_maxWrkVertices
Gaudi::Property< size_t > m_maxWrkVertices
Definition: TrigVrtSecInclusive.h:96
TrigVSI::TrigVrtSecInclusive::m_minTrksMerge
Gaudi::Property< size_t > m_minTrksMerge
Definition: TrigVrtSecInclusive.h:93
Constants.h
Constants for algorithm.
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgCuts::fastRcut
bool fastRcut
Definition: TrigVrtSecInclusive.h:143
xAOD::TrackParticle_v1::charge
float charge() const
Returns the charge.
Definition: TrackParticle_v1.cxx:150
TrigVSI::TrigVrtSecInclusive::m_materialMapInnerFileName
std::string m_materialMapInnerFileName
Definition: TrigVrtSecInclusive.h:180
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
index
Definition: index.py:1
Trk::ParametersT
Dummy class used to allow special convertors to be called for surfaces owned by a detector element.
Definition: EMErrorDetail.h:25
TrigVSI::TrigVrtSecInclusive::findNTrackVertex
StatusCode findNTrackVertex(WrkVrtContainer &, TrigVSI::VtxMap< VrtType, Coord > &, const std::vector< const xAOD::TrackParticle * > &, const EventContext &) const
Definition: TrigVrtSecInclusive.cxx:1151
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgCuts::twoCircErrcut
bool twoCircErrcut
Definition: TrigVrtSecInclusive.h:140
TrigVSI::TrigVrtSecInclusive::WrkVrt::cuts
AlgCuts cuts
Definition: TrigVrtSecInclusive.h:151
TrackParticleBase.h
PGraph.h
xAOD::TrackParticle_v1::summaryValue
bool summaryValue(uint8_t &value, const SummaryType &information) const
Accessor for TrackSummary values.
Definition: TrackParticle_v1.cxx:736
TrigVSI::TrigVrtSecInclusive::m_trkPtCut
Gaudi::Property< double > m_trkPtCut
Definition: TrigVrtSecInclusive.h:62
TrigVSI::TrigVrtSecInclusive::mergeVertexFromDiTrkVrt
StatusCode mergeVertexFromDiTrkVrt(WrkVrtContainer &, const std::vector< std::pair< size_t, size_t >> &, const std::vector< size_t > &, const std::vector< const xAOD::TrackParticle * > &, const EventContext &) const
Reconstruct multi-track vertices from incompatible track pair lists.
Definition: TrigVrtSecInclusive.cxx:1244
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgCuts::isFullPass
bool isFullPass() const
Definition: TrigVrtSecInclusive.h:147
TrigVSI::TrigVrtSecInclusive::m_doTwoCircRCut
Gaudi::Property< bool > m_doTwoCircRCut
Definition: TrigVrtSecInclusive.h:73
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgCuts::twoCircRcut
bool twoCircRcut
Definition: TrigVrtSecInclusive.h:141
TrigVSI::TrigVrtSecInclusive::m_minTrkPairsMerge
Gaudi::Property< size_t > m_minTrkPairsMerge
Definition: TrigVrtSecInclusive.h:92
TrigVSI::TrigVrtSecInclusive::m_z0TrkPVDstMaxCut
Gaudi::Property< double > m_z0TrkPVDstMaxCut
Definition: TrigVrtSecInclusive.h:67
xAOD::TrackParticle_v1::z0
float z0() const
Returns the parameter.
TrigVSI::VtxMap::lock
void lock()
Lock the map to prevent adding vertex after clustering.
Definition: VtxMap.h:169
xAOD::TrackParticle_v1::chiSquared
float chiSquared() const
Returns the of the overall track fit.
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgParam::dphi2
double dphi2
Definition: TrigVrtSecInclusive.h:127
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgParam::vPos
double vPos
Definition: TrigVrtSecInclusive.h:124
DataVector::get
const T * get(size_type n) const
Access an element, as an rvalue.
TrigVSI::VtxMap::Fill
void Fill(const WrkVrt *)
Fill vertex map with vertex from its pointer.
Definition: VtxMap.h:308
TrigVSI::TrigVrtSecInclusive::WrkVrt::closestWrkVrtValue
double closestWrkVrtValue
stores the index of the closest WrkVrt in std::vector<WrkVrt>
Definition: TrigVrtSecInclusive.h:112
TrigVSI::TrigVrtSecInclusive::m_doPVCompatibilityCut
Gaudi::Property< bool > m_doPVCompatibilityCut
Definition: TrigVrtSecInclusive.h:83
xAOD::numberOfPixelHits
@ numberOfPixelHits
these are the pixel hits, including the b-layer [unit8_t].
Definition: TrackingPrimitives.h:259
xAOD::numberOfTRTHits
@ numberOfTRTHits
number of TRT hits [unit8_t].
Definition: TrackingPrimitives.h:275
xAOD::Vertex_v1::position
const Amg::Vector3D & position() const
Returns the 3-pos.
TrigVSI::TrigVrtSecInclusive::WrkVrt::param
AlgParam param
Definition: TrigVrtSecInclusive.h:150
python.TurnDataReader.dr
dr
Definition: TurnDataReader.py:112
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
TrigVSI::TrigVrtSecInclusive::selectTrack_pTCut
bool selectTrack_pTCut(const xAOD::TrackParticle *trk) const
Definition: TrigVrtSecInclusive.cxx:1504
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
TrigVSI::TrigVrtSecInclusive::m_doFastRCut
Gaudi::Property< bool > m_doFastRCut
Definition: TrigVrtSecInclusive.h:74
TrigVSI::TrigVrtSecInclusive::xAODContainers
std::pair< std::unique_ptr< xAOD::VertexContainer >, std::unique_ptr< xAOD::VertexAuxContainer > > xAODContainers
Definition: TrigVrtSecInclusive.h:48
TrigVSI::TrigVrtSecInclusive::m_fastZ0deltaCut
Gaudi::Property< double > m_fastZ0deltaCut
Definition: TrigVrtSecInclusive.h:79
xAOD::Vertex_v1::vertexType
VxType::VertexType vertexType() const
The type of the vertex.
TrigVSI::TrigVrtSecInclusive::WrkVrt::charge
long int charge
list of VKalVrt fit chi2 for each track
Definition: TrigVrtSecInclusive.h:109
TrigVSI::AlgConsts::nPhi
constexpr int nPhi
Default bin number of phi for vertex map.
Definition: Trigger/TrigTools/TrigVrtSecInclusive/TrigVrtSecInclusive/Constants.h:27
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgParam::vrtFast_trkz0
std::vector< float > vrtFast_trkz0
Definition: TrigVrtSecInclusive.h:133
TrigVSI::TrigVrtSecInclusive::trackSelection
StatusCode trackSelection(const xAOD::TrackParticleContainer *, const xAOD::TrackParticleContainer *, std::vector< const xAOD::TrackParticle * > &) const
Definition: TrigVrtSecInclusive.cxx:266
xAOD::TrackParticle_v1::d0
float d0() const
Returns the parameter.
AthenaPoolTestRead.sc
sc
Definition: AthenaPoolTestRead.py:27
TrigVSI::TrigVrtSecInclusive::m_fitSvc
ToolHandle< Trk::TrkVKalVrtFitter > m_fitSvc
Definition: TrigVrtSecInclusive.h:201
PUfitVar::maxEta
constexpr float maxEta
Definition: GepMETPufitAlg.cxx:13
Monitored::Collection
ValuesCollection< T > Collection(std::string name, const T &collection)
Declare a monitored (double-convertible) collection.
Definition: MonitoredCollection.h:38
TrigVSI::TrigVrtSecInclusive::m_vertexPointEstimator
ToolHandle< InDet::VertexPointEstimator > m_vertexPointEstimator
Definition: TrigVrtSecInclusive.h:202
AthReentrantAlgorithm
An algorithm that can be simultaneously executed in multiple threads.
Definition: AthReentrantAlgorithm.h:83
TrigVSI::TrigVrtSecInclusive::m_firstPassTracksName
SG::ReadHandleKey< xAOD::TrackParticleContainer > m_firstPassTracksName
Definition: TrigVrtSecInclusive.h:195
TrigVSI::TrigVrtSecInclusive::m_trkPairOutputName
SG::WriteHandleKey< xAOD::VertexContainer > m_trkPairOutputName
Definition: TrigVrtSecInclusive.h:198
TrigVSI::TrigVrtSecInclusive::selectTrack_z0Cut
bool selectTrack_z0Cut(const xAOD::TrackParticle *trk) const
Definition: TrigVrtSecInclusive.cxx:1503
TrigVSI::TrigVrtSecInclusive::m_vtxAlgorithm
Gaudi::Property< int > m_vtxAlgorithm
Definition: TrigVrtSecInclusive.h:51
dqt_zlumi_pandas.weight
int weight
Definition: dqt_zlumi_pandas.py:200
InDet::VertexPointEstimator::Values_t
std::map< std::string, float > Values_t
Definition: VertexPointEstimator.h:33
DumpGeoConfig.outFileName
string outFileName
Definition: DumpGeoConfig.py:238
xAOD::TrackParticle_v1::perigeeParameters
const Trk::Perigee & perigeeParameters() const
Returns the Trk::MeasuredPerigee track parameters.
Definition: TrackParticle_v1.cxx:485
TrigVSI::TrigVrtSecInclusive::m_trkChi2Cut
Gaudi::Property< double > m_trkChi2Cut
Definition: TrigVrtSecInclusive.h:61
xAOD::numberOfPixelSharedHits
@ numberOfPixelSharedHits
number of Pixel all-layer hits shared by several tracks [unit8_t].
Definition: TrackingPrimitives.h:262
TrigVSI::TrigVrtSecInclusive::findNtrackVerticesVSI
StatusCode findNtrackVerticesVSI(WrkVrtContainer &, std::vector< std::pair< size_t, size_t >> &, std::vector< const xAOD::TrackParticle * > &, const EventContext &) const
Definition: TrigVrtSecInclusive.cxx:910
TrigVSI::TrigVrtSecInclusive::selectTrack_d0Cut
bool selectTrack_d0Cut(const xAOD::TrackParticle *trk) const
Definition: TrigVrtSecInclusive.cxx:1502
TrigVSI::TrigVrtSecInclusive::WrkVrt::chi2PerTrk
std::vector< double > chi2PerTrk
VKalVrt fit chi2 result.
Definition: TrigVrtSecInclusive.h:108
TrigVSI::TrigVrtSecInclusive::m_z0TrkPVDstMinCut
Gaudi::Property< double > m_z0TrkPVDstMinCut
Definition: TrigVrtSecInclusive.h:66
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
ConvertOldHistosToNewHistos.inFileName
inFileName
Definition: ConvertOldHistosToNewHistos.py:21
TrigVSI::AlgConsts::mapBinWid
constexpr float mapBinWid
Default bin width for vertex map.
Definition: Trigger/TrigTools/TrigVrtSecInclusive/TrigVrtSecInclusive/Constants.h:23
LinkToXAODTrackParticle.h
TrigVSI::TrigVrtSecInclusive::m_materialMapMatrix
TMatrixT< double > * m_materialMapMatrix
Definition: TrigVrtSecInclusive.h:179
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgParam::vrtFast_r
double vrtFast_r
Definition: TrigVrtSecInclusive.h:129
fillPileUpNoiseLumi.next
next
Definition: fillPileUpNoiseLumi.py:52
lumiFormat.i
int i
Definition: lumiFormat.py:92
TrigVSI::TrigVrtSecInclusive::m_truncateWrkVertices
Gaudi::Property< bool > m_truncateWrkVertices
Definition: TrigVrtSecInclusive.h:95
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgParam::dphi1
double dphi1
Definition: TrigVrtSecInclusive.h:126
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
SG::ReadHandle::get
const_pointer_type get() const
Dereference the pointer, but don't cache anything.
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
TrigVSI::TrigVrtSecInclusive::WrkVrt::trkAtVrt
std::vector< std::vector< double > > trkAtVrt
total charge of the vertex
Definition: TrigVrtSecInclusive.h:110
TrigVSI::TrigVrtSecInclusive::cleanUp
StatusCode cleanUp(WrkVrtContainer &) const
Definition: TrigVrtSecInclusive.cxx:1561
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgParam
Definition: TrigVrtSecInclusive.h:119
TrigVSI::TrigVrtSecInclusive::finalize
virtual StatusCode finalize() override
Definition: TrigVrtSecInclusive.cxx:1585
master.flag
bool flag
Definition: master.py:29
TrigVSI::TrigVrtSecInclusive::m_cutSharedHits
Gaudi::Property< int > m_cutSharedHits
Definition: TrigVrtSecInclusive.h:57
TrigVSI::TrigVrtSecInclusive::~TrigVrtSecInclusive
virtual ~TrigVrtSecInclusive()
Definition: TrigVrtSecInclusive.cxx:1590
TrigVSI::TrigVrtSecInclusive::m_improveChi2ProbThreshold
Gaudi::Property< double > m_improveChi2ProbThreshold
Definition: TrigVrtSecInclusive.h:97
xAOD::VxType::PriVtx
@ PriVtx
Primary vertex.
Definition: TrackingPrimitives.h:571
TrigVSI::VtxMap::nClusters
size_t nClusters() const
Return the number of the clusters.
Definition: VtxMap.h:188
TrigVSI::VtxMap::getCluster
CellCluster getCluster(size_t)
Retrieve clustering result as CellCluster object.
Definition: VtxMap.h:420
TrigVSI::TrigVrtSecInclusive::m_recordTrkPair
Gaudi::Property< bool > m_recordTrkPair
Definition: TrigVrtSecInclusive.h:53
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgParam::vrtFast_eta
double vrtFast_eta
Definition: TrigVrtSecInclusive.h:130
Coordinate.h
Coordinate policies.
TrigVSI::TrigVrtSecInclusive::initialize
virtual StatusCode initialize() override
Definition: TrigVrtSecInclusive.cxx:38
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgParam::twoCirc_int_r
double twoCirc_int_r
Definition: TrigVrtSecInclusive.h:122
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
TrigVSI::TrigVrtSecInclusive::WrkVrt::vertexCov
std::vector< double > vertexCov
VKalVrt fit vertex 4-momentum.
Definition: TrigVrtSecInclusive.h:105
TrigVSI::TrigVrtSecInclusive::m_materialMapOuterHistName
std::string m_materialMapOuterHistName
Definition: TrigVrtSecInclusive.h:184
TrigVSI::TrigVrtSecInclusive::execute
virtual StatusCode execute(const EventContext &ctx) const override
Definition: TrigVrtSecInclusive.cxx:90
SG::VarHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:103
DeMoUpdate.tmp
string tmp
Definition: DeMoUpdate.py:1167
TrigVSI::VtxMap::ClusterizeCells
void ClusterizeCells(double, size_t)
Generate clusters.
Definition: VtxMap.h:342
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgCuts::chi2cut
bool chi2cut
Definition: TrigVrtSecInclusive.h:145
TrigVSI::TrigVrtSecInclusive::selectTrack
bool selectTrack(const xAOD::TrackParticle *trk) const
Definition: TrigVrtSecInclusive.cxx:1491
DataVector
Derived DataVector<T>.
Definition: DataVector.h:581
Monitored.h
Header file to be included by clients of the Monitored infrastructure.
TrigVSI::TrigVrtSecInclusive::WrkVrt
Definition: TrigVrtSecInclusive.h:100
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgParam::vrtFast_trkd0
std::vector< float > vrtFast_trkd0
Definition: TrigVrtSecInclusive.h:132
min
#define min(a, b)
Definition: cfImp.cxx:40
TrigVSI::TrigVrtSecInclusive::selectTrack_hitPattern
bool selectTrack_hitPattern(const xAOD::TrackParticle *trk) const
Definition: TrigVrtSecInclusive.cxx:1508
TrigVSI::TrigVrtSecInclusive::m_vxCandidatesOutputName
SG::WriteHandleKey< xAOD::VertexContainer > m_vxCandidatesOutputName
Definition: TrigVrtSecInclusive.h:197
LinkToTrackParticleBase.h
VtxMap.h
point classes for clustering
TrigVSI::TrigVrtSecInclusive::m_d0TrkPVDstMaxCut
Gaudi::Property< double > m_d0TrkPVDstMaxCut
Definition: TrigVrtSecInclusive.h:65
TrigVSI::TrigVrtSecInclusive::WrkVrt::chi2
double chi2
VKalVrt fit covariance.
Definition: TrigVrtSecInclusive.h:106
PathResolver.h
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:192
TH3S::GetBinContent
double GetBinContent(int) const
Definition: rootspy.cxx:480
TrigVSI::TrigVrtSecInclusive::m_materialMapInnerMatrixName
std::string m_materialMapInnerMatrixName
Definition: TrigVrtSecInclusive.h:182
TrigVSI::VtxMap
The vertex map class to be used to find multi-track vertices.
Definition: VtxMap.h:40
TrigVSI::TrigVrtSecInclusive::m_cutSctHits
Gaudi::Property< int > m_cutSctHits
Definition: TrigVrtSecInclusive.h:56
TrigVSI::AlgConsts
Namespace for constants.
Definition: Trigger/TrigTools/TrigVrtSecInclusive/TrigVrtSecInclusive/Constants.h:19
TrigVSI::TrigVrtSecInclusive::m_monTool
ToolHandle< GenericMonitoringTool > m_monTool
Definition: TrigVrtSecInclusive.h:204
TrigVSI::TrigVrtSecInclusive::m_materialMapInnerHistName
std::string m_materialMapInnerHistName
Definition: TrigVrtSecInclusive.h:181
TrigVrtSecInclusive.h
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
TrigVSI::TrigVrtSecInclusive::m_cutPixelHits
Gaudi::Property< int > m_cutPixelHits
Definition: TrigVrtSecInclusive.h:55
TrigVSI::TrigVrtSecInclusive::m_dphiPVCut
Gaudi::Property< double > m_dphiPVCut
Definition: TrigVrtSecInclusive.h:84
TrigVSI::TrigVrtSecInclusive::m_fastD0deltaCut
Gaudi::Property< double > m_fastD0deltaCut
Definition: TrigVrtSecInclusive.h:77
TrigVSI::TrigVrtSecInclusive::m_fastD0minCut
Gaudi::Property< double > m_fastD0minCut
Definition: TrigVrtSecInclusive.h:76
KDPoint.h
point classes for clustering
DataVector::end
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
SG::WriteHandle
Definition: StoreGate/StoreGate/WriteHandle.h:76
Trk::vertex
@ vertex
Definition: MeasurementType.h:21
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
TrigVSI::AlgConsts::binNumR
constexpr int binNumR
Default bin number of R for vertex map.
Definition: Trigger/TrigTools/TrigVrtSecInclusive/TrigVrtSecInclusive/Constants.h:29
TrigVSI::TrigVrtSecInclusive::m_maxTrks
Gaudi::Property< size_t > m_maxTrks
Definition: TrigVrtSecInclusive.h:90
xAOD::Vertex_v1
Class describing a Vertex.
Definition: Vertex_v1.h:42
SG::WriteHandle::record
StatusCode record(std::unique_ptr< T > data)
Record a const object to the store.
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
TrigVSI::TrigVrtSecInclusive::WrkVrtContainer
std::vector< WrkVrt > WrkVrtContainer
Definition: TrigVrtSecInclusive.h:174
AthCommonMsg< Gaudi::Algorithm >::msg
MsgStream & msg() const
Definition: AthCommonMsg.h:24
TrigVSI::TrigVrtSecInclusive::selectTrack_chi2Cut
bool selectTrack_chi2Cut(const xAOD::TrackParticle *trk) const
Definition: TrigVrtSecInclusive.cxx:1505
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgParam::vrtFast_phi
double vrtFast_phi
Definition: TrigVrtSecInclusive.h:131
TrigVSI::TrigVrtSecInclusive::WrkVrt::fitQuality
double fitQuality() const
Definition: TrigVrtSecInclusive.h:117
TrigVSI::TrigVrtSecInclusive::m_cutBLayHits
Gaudi::Property< int > m_cutBLayHits
Definition: TrigVrtSecInclusive.h:59
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgCuts::fastErrcut
bool fastErrcut
Definition: TrigVrtSecInclusive.h:142
xAOD::numberOfSCTHits
@ numberOfSCTHits
number of hits in SCT [unit8_t].
Definition: TrackingPrimitives.h:268
TrigVSI::TrigVrtSecInclusive::m_selVrtChi2Cut
Gaudi::Property< double > m_selVrtChi2Cut
Definition: TrigVrtSecInclusive.h:81
TrigVSI::AlgConsts::IDrInner
constexpr float IDrInner
Minimum R of ID volume.
Definition: Trigger/TrigTools/TrigVrtSecInclusive/TrigVrtSecInclusive/Constants.h:21
xAOD::TrackParticle_v1::numberDoF
float numberDoF() const
Returns the number of degrees of freedom of the overall track or vertex fit as float.
xAOD::TrackParticle_v1
Class describing a TrackParticle.
Definition: TrackParticle_v1.h:43
TrigVSI::TrigVrtSecInclusive::nTrkCommon
size_t nTrkCommon(WrkVrtContainer &, const std::pair< unsigned, unsigned > &) const
Definition: TrigVrtSecInclusive.cxx:1538
DataVector::at
const T * at(size_type n) const
Access an element, as an rvalue.
TrigVSI::AlgConsts::nEta
constexpr int nEta
Default bin number of eta for vertex map.
Definition: Trigger/TrigTools/TrigVrtSecInclusive/TrigVrtSecInclusive/Constants.h:26
TrigVSI::TrigVrtSecInclusive::m_PrimaryVxInputName
SG::ReadHandleKey< xAOD::VertexContainer > m_PrimaryVxInputName
Definition: TrigVrtSecInclusive.h:199
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgParam::vPosMomAngT
double vPosMomAngT
Definition: TrigVrtSecInclusive.h:125
TrigVSI::TrigVrtSecInclusive::WrkVrt::isPair
bool isPair
flaged true for good vertex candidates
Definition: TrigVrtSecInclusive.h:102
TrigVSI::TrigVrtSecInclusive::m_twoTrkVtxFormingD0Cut
Gaudi::Property< double > m_twoTrkVtxFormingD0Cut
Definition: TrigVrtSecInclusive.h:68
TrigVSI::TrigVrtSecInclusive::m_materialMapOuter
TH3S * m_materialMapOuter
Definition: TrigVrtSecInclusive.h:178
TrigVSI::TrigVrtSecInclusive::WrkVrt::vertexMom
TLorentzVector vertexMom
VKalVrt fit vertex position.
Definition: TrigVrtSecInclusive.h:104
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgCuts::fitErrcut
bool fitErrcut
Definition: TrigVrtSecInclusive.h:144
TrigVSI::TrigVrtSecInclusive::m_fastZ0minCut
Gaudi::Property< double > m_fastZ0minCut
Definition: TrigVrtSecInclusive.h:78
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
TrigVSI::TrigVrtSecInclusive::m_pvCompatibilityCut
Gaudi::Property< double > m_pvCompatibilityCut
Definition: TrigVrtSecInclusive.h:85
TrigVSI::TrigVrtSecInclusive::WrkVrt::isGood
bool isGood
Definition: TrigVrtSecInclusive.h:101
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
Decorator.h
Helper class to provide type-safe access to aux data.
Monitored::Timer
A monitored timer.
Definition: MonitoredTimer.h:32
TrigVSI::TrigVrtSecInclusive::m_cutSiHits
Gaudi::Property< int > m_cutSiHits
Definition: TrigVrtSecInclusive.h:58
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgCuts::trkd0cut
bool trkd0cut
Definition: TrigVrtSecInclusive.h:139
TrigVSI
Definition: TrigVrtSecInclusive.cxx:27
TrigVSI::TrigVrtSecInclusive::m_d0TrkPVDstMinCut
Gaudi::Property< double > m_d0TrkPVDstMinCut
Definition: TrigVrtSecInclusive.h:64
TrigVSI::TrigVrtSecInclusive::m_secondPassTracksName
SG::ReadHandleKey< xAOD::TrackParticleContainer > m_secondPassTracksName
Definition: TrigVrtSecInclusive.h:196
DataVector::begin
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
TrigVSI::TrigVrtSecInclusive::WrkVrt::AlgCuts
Definition: TrigVrtSecInclusive.h:138
xAOD::TrackParticle_v1::phi
virtual double phi() const override final
The azimuthal angle ( ) of the particle (has range to .)
xAOD::numberOfInnermostPixelLayerHits
@ numberOfInnermostPixelLayerHits
these are the hits in the 0th pixel barrel layer
Definition: TrackingPrimitives.h:237
TrigVSI::TrigVrtSecInclusive::WrkVrt::closestWrkVrtIndex
unsigned long closestWrkVrtIndex
list of track parameters wrt the reconstructed vertex
Definition: TrigVrtSecInclusive.h:111
LArGeo::ATan2
GeoGenfun::FunctionNoop ATan2(GeoGenfun::GENFUNCTION y, GeoGenfun::GENFUNCTION x)
Definition: BarrelAuxFunctions.cxx:50
TrigVSI::TrigVrtSecInclusive::m_maxR
Gaudi::Property< double > m_maxR
Definition: TrigVrtSecInclusive.h:71
Generated on Sun Jun 30 2024 21:31:08 for ATLAS Offline Software by   doxygen 1.8.18