36 bool is_pv_associated =
false;
40 for(
const auto* vtx : *vertices ) {
41 for(
size_t iv = 0; iv < vtx->nTrackParticles(); iv++ ) {
42 const auto* pvtrk = vtx->trackParticle( iv );
43 if (pvtrk ==
nullptr)
continue;
45 if ( (trk_from_gsf == pvtrk) or (trk == pvtrk) ) {
46 is_pv_associated =
true;
51 return is_pv_associated;
56 return (v1-
v2).norm();
61 double VrtSecInclusive::significanceBetweenVertices(
const WrkVrt& v1,
const WrkVrt&
v2 )
const {
66 sumCov.fillSymmetric(0, 0, v1.
vertexCov.at(0) +
v2.vertexCov.at(0));
67 sumCov.fillSymmetric(1, 0, v1.
vertexCov.at(1) +
v2.vertexCov.at(1));
68 sumCov.fillSymmetric(1, 1, v1.
vertexCov.at(2) +
v2.vertexCov.at(2));
69 sumCov.fillSymmetric(2, 0, v1.
vertexCov.at(3) +
v2.vertexCov.at(3));
70 sumCov.fillSymmetric(2, 1, v1.
vertexCov.at(4) +
v2.vertexCov.at(4));
71 sumCov.fillSymmetric(2, 2, v1.
vertexCov.at(5) +
v2.vertexCov.at(5));
77 ATH_MSG_WARNING(
" >>> " << __FUNCTION__ <<
": detected covariance matrix broken exception" );
86 xAOD::Muon::InnerDetectorTrackParticle,
87 xAOD::Muon::MuonSpectrometerTrackParticle,
88 xAOD::Muon::CombinedTrackParticle,
89 xAOD::Muon::ExtrapolatedMuonSpectrometerTrackParticle
95 for(
size_t i=0;
i<
electron->nTrackParticles();
i++ ) trackTypes.emplace_back(
i );
102 return muon->trackParticle(
static_cast<xAOD::Muon::TrackParticleType
>( trackType ) );
107 return electron->trackParticle( trackType );
112 double VrtSecInclusive::distanceBetweenVertices(
const WrkVrt& v1,
const WrkVrt&
v2 )
const {
113 return (
v2.vertex - v1.
vertex).norm();
118 StatusCode VrtSecInclusive::disassembleVertex(std::vector<WrkVrt> *workVerticesContainer,
const unsigned& iv )
121 auto& wrkvrt = workVerticesContainer->at(iv);
123 ATH_MSG_VERBOSE(
" >> disassembleVertex(): begin: disassembling vertex[" << iv <<
"], workVerticesContainer.size() = " << workVerticesContainer->size() );
124 ATH_MSG_VERBOSE(
" >> disassembleVertex(): Vertex: r = " << wrkvrt.vertex.perp() <<
", z = " << wrkvrt.vertex.z() );
127 const auto& ntrk = wrkvrt.selectedTrackIndices.size();
131 if( wrkvrt.selectedTrackIndices.size() <= 2 )
return StatusCode::SUCCESS;
133 for(
auto&
index : wrkvrt.selectedTrackIndices ) {
136 ATH_MSG_VERBOSE(
" >> disassembleVertex(): > track at vertex[" << iv <<
"]: "
137 <<
"index = " << trk->
index()
138 <<
", pT = " << trk->
pt()
139 <<
", phi = " << trk->
phi()
140 <<
", d0 = " << trk->
d0()
141 <<
", z0 = " << trk->
z0());
145 const auto&
max = std::max_element( wrkvrt.Chi2PerTrk.begin(), wrkvrt.Chi2PerTrk.end() );
147 if(
max == wrkvrt.Chi2PerTrk.end() )
return StatusCode::SUCCESS;
149 maxChi2TrackIndex =
max - wrkvrt.Chi2PerTrk.begin();
156 vector<const xAOD::NeutralParticle*> dummyNeutrals;
158 vector<WrkVrt> new_vertices;
161 ATH_MSG_VERBOSE(
" >> disassembleVertex(): Loop over the tracks associated to the vertex other than the selected tracks.");
162 for(
size_t itrk=0; itrk<ntrk; itrk++) {
164 ATH_MSG_VERBOSE(
" >> disassembleVertex(): > Loop itrk = " << itrk <<
" / " << ntrk );
167 if( itrk == maxChi2TrackIndex ) {
172 const size_t this_trk_id = wrkvrt.selectedTrackIndices[itrk];
173 const size_t selected_trk_id = wrkvrt.selectedTrackIndices[maxChi2TrackIndex];
175 ATH_MSG_VERBOSE(
" >> disassembleVertex(): > this_trk_id = " << this_trk_id <<
", selected_trk_id = " << selected_trk_id <<
", alltrks_size = " << m_selectedTracks.size() );
176 if( this_trk_id >= m_selectedTracks.size() ) {
177 ATH_MSG_VERBOSE(
" >> disassembleVertex(): > this_trk_id is invalid. continue!" );
180 if( selected_trk_id >= m_selectedTracks.size() ) {
181 ATH_MSG_VERBOSE(
" >> disassembleVertex(): > selected_trk_id is invalid. continue!" );
185 ATH_MSG_VERBOSE(
" >> disassembleVertex(): > Storing tracks to ListBaseTracks" );
186 ATH_MSG_VERBOSE(
" >> disassembleVertex(): > m_selectedTracks.at( this_trk_id ) = " << m_selectedTracks.at( this_trk_id )->index() );
187 ATH_MSG_VERBOSE(
" >> disassembleVertex(): > m_selectedTracks.at( this_trk_id ) = " << m_selectedTracks.at( selected_trk_id )->index() );
189 vector<const xAOD::TrackParticle*> ListBaseTracks;
190 ListBaseTracks.emplace_back( m_selectedTracks.at( this_trk_id ) );
191 ListBaseTracks.emplace_back( m_selectedTracks.at( selected_trk_id ) );
193 ATH_MSG_VERBOSE(
" >> disassembleVertex(): > ListBaseTracks was stored." );
202 std::unique_ptr<Trk::IVKalState> state = m_fitSvc->makeState();
203 ATH_CHECK( m_fitSvc->VKalVrtFitFast( ListBaseTracks, newvrt.
vertex, *state ) );
209 m_fitSvc->setApproximateVertex( wrkvrt.vertex[0], wrkvrt.vertex[1], wrkvrt.vertex[2], *state );
213 m_fitSvc->setApproximateVertex( newvrt.
vertex[0], newvrt.
vertex[1], newvrt.
vertex[2], *state );
228 if(
sc.isFailure() )
continue;
234 ATH_MSG_VERBOSE(
" >> disassembleVertex(): > register the new vertex to the vertex list" );
235 new_vertices.emplace_back( newvrt );
239 wrkvrt.selectedTrackIndices.erase( wrkvrt.selectedTrackIndices.begin() + maxChi2TrackIndex );
240 ATH_MSG_VERBOSE(
" >> disassembleVertex(): removed the selected track from the original vertex. wrkvrt.selectedTrackIndices.size = " << wrkvrt.selectedTrackIndices.size() );
243 ATH_MSG_VERBOSE(
" >> disassembleVertex(): refit the original vertex" );
246 if(
sc.isFailure() ) {
249 return StatusCode::SUCCESS;
253 for(
const auto&
vertex : new_vertices ) {
254 ATH_MSG_VERBOSE(
" >> disassembleVertex(): > emplace_back new vertex" );
255 workVerticesContainer->emplace_back(
vertex );
258 ATH_MSG_VERBOSE(
" >> disassembleVertex(): end. workVerticesContainer.size() = " << workVerticesContainer->size() );
259 return StatusCode::SUCCESS;
271 auto fitQuality_begin =
vertex.fitQuality();
273 auto removeCounter = 0;
275 if(
vertex.nTracksTotal() <= 2 )
return 0.;
280 if(
sc.isFailure() ) {
286 double chi2Probability = TMath::Prob(
vertex.Chi2,
vertex.ndof() );
288 while (chi2Probability < m_jp.improveChi2ProbThreshold ) {
289 if(
vertex.nTracksTotal() == 2 )
return chi2Probability;
293 auto maxChi2 = std::max_element(
vertex.Chi2PerTrk.begin(),
vertex.Chi2PerTrk.end() );
294 size_t index = maxChi2 -
vertex.Chi2PerTrk.begin();
297 ATH_MSG_DEBUG(
" >>> " << __FUNCTION__ <<
": worst chi2 trk index = " <<
index <<
", #trks = " <<
vertex.Chi2PerTrk.size() );
300 vertex.selectedTrackIndices.erase(
vertex.selectedTrackIndices.begin() +
index );
304 if(
index >=
vertex.associatedTrackIndices.size() ) {
308 vertex.associatedTrackIndices.erase(
vertex.associatedTrackIndices.begin() +
index );
320 chi2Probability = TMath::Prob(
vertex.Chi2,
vertex.ndof() );
323 auto fitQuality_end =
vertex.fitQuality();
325 if( 0 == removeCounter ) {
326 ATH_MSG_DEBUG(
" >>> " << __FUNCTION__ <<
": no improvement was found." );
328 ATH_MSG_DEBUG(
" >>> " << __FUNCTION__ <<
": Removed " << removeCounter <<
" tracks; Fit quality improvement: " << fitQuality_begin <<
" ==> " << fitQuality_end );
331 return chi2Probability;
337 void VrtSecInclusive::removeTrackFromVertex(std::vector<WrkVrt> *workVerticesContainer,
338 std::vector< std::deque<long int> > *TrkInVrt,
339 const long int & trackIndexToRemove,
340 const long int & SelectedVertex)
343 auto& wrkvrt = workVerticesContainer->at(SelectedVertex);
344 auto& tracks = wrkvrt.selectedTrackIndices;
347 auto end =
std::remove_if( tracks.begin(), tracks.end(), [&](
long int index) { return index == trackIndexToRemove; } );
348 tracks.erase(
end, tracks.end() );
352 for(
auto& trks : *TrkInVrt ) {
353 auto end =
std::remove_if( trks.begin(), trks.end(), [&](
long int index) { return index == trackIndexToRemove; } );
354 trks.erase(
end, trks.end() );
360 if( wrkvrt.selectedTrackIndices.size() == 1 ) {
362 const auto& leftTrackIndex = *( tracks.begin() );
363 auto&
list = TrkInVrt->at(leftTrackIndex);
380 for(
auto& workVertex : *workVerticesContainer ) {
382 workVertex.closestWrkVrtIndex = 0;
387 for(
auto iv = workVerticesContainer->begin(); iv != workVerticesContainer->end(); ++iv ) {
388 if( (*iv).selectedTrackIndices.size()< 2)
continue;
390 auto i_index = iv - workVerticesContainer->begin();
392 for(
auto jv =
std::next(iv); jv != workVerticesContainer->end(); ++jv ) {
393 if( (*jv).selectedTrackIndices.size()< 2)
continue;
395 auto j_index = iv - workVerticesContainer->begin();
401 indexPair.first = i_index;
402 indexPair.second = j_index;
404 if(
value < (*iv).closestWrkVrtValue ) {(*iv).closestWrkVrtValue =
value; (*iv).closestWrkVrtIndex = j_index; }
405 if(
value < (*jv).closestWrkVrtValue ) {(*jv).closestWrkVrtValue =
value; (*jv).closestWrkVrtIndex = i_index; }
414 double VrtSecInclusive::findMinVerticesNextPair( std::vector<WrkVrt> *workVerticesContainer, std::pair<unsigned, unsigned>& indexPair )
421 indexPair.second = 0;
425 for(
unsigned iv=0; iv<workVerticesContainer->size()-1; iv++) {
426 auto&
vertex = workVerticesContainer->at(iv);
428 if(
vertex.selectedTrackIndices.size() < 2)
continue;
429 if(
vertex.closestWrkVrtIndex == 0 )
continue;
433 unsigned jv =
vertex.closestWrkVrtIndex;
436 if( workVerticesContainer->at(jv).closestWrkVrtIndex == 0 )
continue;
440 indexPair.first = iv;
441 indexPair.second = jv;
471 v2.selectedTrackIndices.clear();
473 v2.closestWrkVrtIndex = 0;
481 if(
sc.isFailure() ) {
485 ATH_MSG_DEBUG(
" >>> " << __FUNCTION__ <<
": failure in merging" );
487 return StatusCode::FAILURE;
490 return StatusCode::SUCCESS;
498 std::unique_ptr<Trk::IVKalState> state = m_fitSvc->makeState();
499 return refitVertex (workVertex, *state);
510 vector<const xAOD::NeutralParticle*> dummyNeutrals;
515 workVertex.
isGood =
false;
516 return StatusCode::SUCCESS;
519 vector<const xAOD::TrackParticle*> ListBaseTracks;
524 ListBaseTracks.emplace_back( m_selectedTracks.at(
index ) );
529 ListBaseTracks.emplace_back( m_associatedTracks.at(
index ) );
533 auto& vertexPos = workVertex.
vertex;
535 m_fitSvc->setApproximateVertex( vertexPos.x(), vertexPos.y(), vertexPos.z(), istate );
537 ATH_MSG_VERBOSE(
" >>> refitVertex: ListBaseTracks.size = " << ListBaseTracks.size()
540 for(
const auto *trk : ListBaseTracks ) {
548 StatusCode sc = m_fitSvc->VKalVrtFitFast( ListBaseTracks, initVertex, istate );
549 if(
sc.isFailure())
ATH_MSG_DEBUG(
" >>> refitVertex: fast crude estimation failed.");
550 ATH_MSG_VERBOSE(
" >>> refitVertex: Fast VKalVrtFit succeeded. vertex (r,z) = (" << initVertex.perp() <<
", " << initVertex.z() <<
", " <<
")" );
554 m_fitSvc->setApproximateVertex( vertexPos.x(), vertexPos.y(), vertexPos.z(), istate );
558 m_fitSvc->setApproximateVertex( initVertex.x(), initVertex.y(), initVertex.z(), istate );
562 ATH_MSG_VERBOSE(
" >>> refitVertex: approx vertex is set. Now going to perform fitting..." );
564 StatusCode SC=m_fitSvc->VKalVrtFit(ListBaseTracks,dummyNeutrals,
576 if(
SC.isFailure())
ATH_MSG_DEBUG(
" >>> refitVertex: SC in refitVertex returned failure ");
578 ATH_MSG_VERBOSE(
" >>> refitVertex: succeeded in fitting. New vertex pos (r,z) = (" << vertexPos.perp() <<
", " << vertexPos.z() <<
")" );
579 ATH_MSG_VERBOSE(
" >>> refitVertex: New vertex cov = (" <<
cov.at(0) <<
", " <<
cov.at(1) <<
", " <<
cov.at(2) <<
", " <<
cov.at(3) <<
", " <<
cov.at(4) <<
", " <<
cov.at(5) <<
")" );
587 std::unique_ptr<Trk::IVKalState> state = m_fitSvc->makeState();
588 return refitVertexWithSuggestion (workVertex, suggestedPosition, *state);
599 vector<const xAOD::NeutralParticle*> dummyNeutrals;
604 workVertex.
isGood =
false;
605 return StatusCode::SUCCESS;
608 vector<const xAOD::TrackParticle*> ListBaseTracks;
613 ListBaseTracks.emplace_back( m_selectedTracks.at(
index ) );
618 ListBaseTracks.emplace_back( m_associatedTracks.at(
index ) );
622 auto& vertexPos = workVertex.
vertex;
624 m_fitSvc->setApproximateVertex( suggestedPosition.x(), suggestedPosition.y(), suggestedPosition.z(), istate );
626 ATH_MSG_VERBOSE(
" >>> " << __FUNCTION__ <<
": ListBaseTracks.size = " << ListBaseTracks.size()
629 for(
const auto *trk : ListBaseTracks ) {
630 ATH_MSG_VERBOSE(
" >>> " << __FUNCTION__ <<
": track index = " << trk->index() );
635 ATH_MSG_VERBOSE(
" >>> " << __FUNCTION__ <<
": approx vertex is set. Now going to perform fitting..." );
637 StatusCode SC=m_fitSvc->VKalVrtFit(ListBaseTracks,dummyNeutrals,
649 if(
SC.isFailure() )
ATH_MSG_VERBOSE(
" >>> " << __FUNCTION__ <<
": SC in refitVertex returned failure ");
652 if(
SC.isSuccess() ) {
653 ATH_MSG_VERBOSE(
" >>> " << __FUNCTION__ <<
": succeeded in fitting. New vertex pos = (" << vertexPos.x() <<
", " << vertexPos.y() <<
", " << vertexPos.z() <<
")" );
654 ATH_MSG_VERBOSE(
" >>> " << __FUNCTION__ <<
": New vertex cov = (" <<
cov.at(0) <<
", " <<
cov.at(1) <<
", " <<
cov.at(2) <<
", " <<
cov.at(3) <<
", " <<
cov.at(4) <<
", " <<
cov.at(5) <<
")" );
662 size_t VrtSecInclusive::nTrkCommon( std::vector<WrkVrt> *workVerticesContainer,
const std::pair<unsigned, unsigned>& pairIndex)
668 auto& trackIndices1 = workVerticesContainer->at( pairIndex.first ).selectedTrackIndices;
669 auto& trackIndices2 = workVerticesContainer->at( pairIndex.second ).selectedTrackIndices;
671 if( trackIndices1.size() < 2 )
return 0;
672 if( trackIndices2.size() < 2 )
return 0;
676 for(
auto&
index : trackIndices1 ) {
677 if(
std::find(trackIndices2.begin(),trackIndices2.end(),
index) != trackIndices2.end()) nTrkCom++;
690 declareProperty(
"TrackLocation", m_jp.TrackLocation =
"InDetTrackParticles" );
691 declareProperty(
"MuonLocation", m_jp.MuonLocation =
"Muons" );
692 declareProperty(
"ElectronLocation", m_jp.ElectronLocation =
"Electrons" );
693 declareProperty(
"PrimVrtLocation", m_jp.PrimVrtLocation =
"PrimaryVertices" );
694 declareProperty(
"McParticleContainer", m_jp.truthParticleContainerName =
"TruthParticles" );
695 declareProperty(
"MCEventContainer", m_jp.mcEventContainerName =
"TruthEvents" );
696 declareProperty(
"AugmentingVersionString", m_jp.augVerString =
"_VSI" );
697 declareProperty(
"TruthParticleFilter", m_jp.truthParticleFilter =
"Rhadron" );
699 declareProperty(
"All2trkVerticesContainerName", m_jp.all2trksVerticesContainerName =
"All2TrksVertices" );
700 declareProperty(
"SecondaryVerticesContainerName", m_jp.secondaryVerticesContainerName =
"SecondaryVertices" );
702 declareProperty(
"FillHist", m_jp.FillHist =
false );
703 declareProperty(
"FillNtuple", m_jp.FillNtuple =
false );
704 declareProperty(
"FillIntermediateVertices", m_jp.FillIntermediateVertices =
false );
705 declareProperty(
"DoIntersectionPos", m_jp.doIntersectionPos =
false );
706 declareProperty(
"DoMapToLocal", m_jp.doMapToLocal =
false );
707 declareProperty(
"DoTruth", m_jp.doTruth =
false );
708 declareProperty(
"DoPVcompatibility", m_jp.doPVcompatibilityCut =
true );
709 declareProperty(
"DoTightPVcompatibility", m_jp.doTightPVcompatibilityCut =
false );
710 declareProperty(
"RemoveFake2TrkVrt", m_jp.removeFakeVrt =
true );
711 declareProperty(
"DoDelayedFakeReject", m_jp.removeFakeVrtLate =
false );
712 declareProperty(
"CheckHitPatternStrategy", m_checkPatternStrategy =
"Classical" );
713 declareProperty(
"MCTrackResolution", m_jp.mcTrkResolution = 0.06 );
714 declareProperty(
"TruthTrkLen", m_jp.TruthTrkLen = 1000 );
715 declareProperty(
"ExtrapPV", m_jp.extrapPV =
false );
716 declareProperty(
"PassThroughTrackSelection", m_jp.passThroughTrackSelection =
false );
717 declareProperty(
"DoFastMode", m_jp.doFastMode =
false );
720 declareProperty(
"DoTwoTrSoftBtag", m_jp.doTwoTrSoftBtag =
false );
721 declareProperty(
"TwoTrVrtAngleCut", m_jp.twoTrVrtAngleCut = -10 );
722 declareProperty(
"TwoTrVrtMinDistFromPVCut", m_jp.twoTrVrtMinDistFromPV = 0. );
724 declareProperty(
"TruncateListOfWorkingVertices", m_jp.truncateWrkVertices =
true );
725 declareProperty(
"MaxNumberOfWorkingVertices", m_jp.maxWrkVertices = 1500 );
728 declareProperty(
"do_PVvetoCut", m_jp.do_PVvetoCut =
true );
729 declareProperty(
"do_d0Cut", m_jp.do_d0Cut =
true );
730 declareProperty(
"do_z0Cut", m_jp.do_z0Cut =
true );
731 declareProperty(
"do_d0errCut", m_jp.do_d0errCut =
false );
732 declareProperty(
"do_z0errCut", m_jp.do_z0errCut =
false );
733 declareProperty(
"do_d0signifCut", m_jp.do_d0signifCut =
false );
734 declareProperty(
"do_z0signifCut", m_jp.do_z0signifCut =
false );
736 declareProperty(
"ImpactWrtBL", m_jp.ImpactWrtBL =
true );
737 declareProperty(
"a0TrkPVDstMinCut", m_jp.d0TrkPVDstMinCut = 0. );
738 declareProperty(
"a0TrkPVDstMaxCut", m_jp.d0TrkPVDstMaxCut = 1000. );
739 declareProperty(
"a0TrkPVSignifCut", m_jp.d0TrkPVSignifCut = 0. );
740 declareProperty(
"twoTrkVtxFormingD0Cut", m_jp.twoTrkVtxFormingD0Cut = 1. );
741 declareProperty(
"zTrkPVDstMinCut", m_jp.z0TrkPVDstMinCut = 0. );
742 declareProperty(
"zTrkPVDstMaxCut", m_jp.z0TrkPVDstMaxCut = 1000. );
743 declareProperty(
"zTrkPVSignifCut", m_jp.z0TrkPVSignifCut = 0. );
744 declareProperty(
"TrkA0ErrCut", m_jp.d0TrkErrorCut = 10000 );
745 declareProperty(
"TrkZErrCut", m_jp.z0TrkErrorCut = 20000 );
747 declareProperty(
"SelTrkMaxCutoff", m_jp.SelTrkMaxCutoff = 50 );
748 declareProperty(
"TrkPtCut", m_jp.TrkPtCut = 1000. );
749 declareProperty(
"TrkChi2Cut", m_jp.TrkChi2Cut = 3. );
750 declareProperty(
"PVcompatibilityCut", m_jp.pvCompatibilityCut = -20. );
751 declareProperty(
"SelVrtChi2Cut", m_jp.SelVrtChi2Cut = 4.5 );
753 declareProperty(
"CutSctHits", m_jp.CutSctHits = 0 );
754 declareProperty(
"CutPixelHits", m_jp.CutPixelHits = 0 );
755 declareProperty(
"CutSiHits", m_jp.CutSiHits = 0 );
756 declareProperty(
"DoSAloneTRT", m_jp.SAloneTRT =
false );
757 declareProperty(
"CutBLayHits", m_jp.CutBLayHits = 0 );
758 declareProperty(
"CutSharedHits", m_jp.CutSharedHits = 0 );
759 declareProperty(
"doTRTPixCut", m_jp.doTRTPixCut =
false );
760 declareProperty(
"CutTRTHits", m_jp.CutTRTHits = 0 );
761 declareProperty(
"CutTightSCTHits", m_jp.CutTightSCTHits = 7 );
762 declareProperty(
"CutTightTRTHits", m_jp.CutTightTRTHits = 20 );
764 declareProperty(
"TrkExtrapolator", m_jp.trkExtrapolator = 2 );
766 declareProperty(
"doReassembleVertices", m_jp.doReassembleVertices =
false );
767 declareProperty(
"doMergeByShuffling", m_jp.doMergeByShuffling =
false );
768 declareProperty(
"doSuggestedRefitOnMerging", m_jp.doSuggestedRefitOnMerging =
true );
769 declareProperty(
"doMagnetMerging", m_jp.doMagnetMerging =
true );
770 declareProperty(
"doWildMerging", m_jp.doWildMerging =
true );
771 declareProperty(
"doMergeFinalVerticesDistance", m_jp.doMergeFinalVerticesDistance =
false );
772 declareProperty(
"doAssociateNonSelectedTracks", m_jp.doAssociateNonSelectedTracks =
false );
773 declareProperty(
"doFinalImproveChi2", m_jp.doFinalImproveChi2 =
false );
775 declareProperty(
"VertexMergeCut", m_jp.VertexMergeCut = 3 );
776 declareProperty(
"TrackDetachCut", m_jp.TrackDetachCut = 6 );
777 declareProperty(
"associateMinDistanceToPV", m_jp.associateMinDistanceToPV = 0.5 );
778 declareProperty(
"associateMaxD0Signif", m_jp.associateMaxD0Signif = 5. );
779 declareProperty(
"associateMaxZ0Signif", m_jp.associateMaxZ0Signif = 5. );
780 declareProperty(
"associatePtCut", m_jp.associatePtCut = 0. );
781 declareProperty(
"associateChi2Cut", m_jp.associateChi2Cut = 20. );
782 declareProperty(
"reassembleMaxImpactParameterD0", m_jp.reassembleMaxImpactParameterD0 = 1. );
783 declareProperty(
"reassembleMaxImpactParameterZ0", m_jp.reassembleMaxImpactParameterZ0 = 5. );
784 declareProperty(
"mergeByShufflingMaxSignificance", m_jp.mergeByShufflingMaxSignificance = 100. );
785 declareProperty(
"mergeByShufflingAllowance", m_jp.mergeByShufflingAllowance = 4. );
786 declareProperty(
"VertexMergeFinalDistCut", m_jp.VertexMergeFinalDistCut = 1. );
787 declareProperty(
"VertexMergeFinalDistScaling", m_jp.VertexMergeFinalDistScaling = 0. );
788 declareProperty(
"improveChi2ProbThreshold", m_jp.improveChi2ProbThreshold = 1.e-4 );
791 declareProperty(
"doSelectTracksFromMuons", m_jp.doSelectTracksFromMuons =
false );
792 declareProperty(
"doRemoveCaloTaggedMuons", m_jp.doRemoveCaloTaggedMuons =
false );
793 declareProperty(
"doSelectTracksFromElectrons", m_jp.doSelectTracksFromElectrons =
false );
794 declareProperty(
"doSelectIDAndGSFTracks", m_jp.doSelectIDAndGSFTracks =
false );
795 declareProperty(
"doRemoveNonLeptonVertices", m_jp.doRemoveNonLeptonVertices =
false );
798 declareProperty(
"doSelectTracksWithLRTCuts", m_jp.doSelectTracksWithLRTCuts =
false );
801 declareProperty(
"doAugmentDVimpactParametersToMuons", m_jp.doAugmentDVimpactParametersToMuons =
false );
802 declareProperty(
"doAugmentDVimpactParametersToElectrons", m_jp.doAugmentDVimpactParametersToElectrons =
false );
805 declareProperty(
"VertexFitterTool", m_fitSvc,
" Private TrkVKalVrtFitter" );
806 declareProperty(
"Extrapolator", m_extrapolator );
807 declareProperty(
"TrackToVertexTool", m_trackToVertexTool );
808 declareProperty(
"TrackToVertexIPEstimatorTool", m_trackToVertexIPEstimatorTool );
809 declareProperty(
"VertexMapper", m_vertexMapper );
810 declareProperty(
"TruthToTrack", m_truthToTrack );
826 if( m_jp.FillNtuple ) m_ntupleVars->get<
unsigned int>(
"NumPV" ) = 0;
829 ATH_MSG_DEBUG(
"processPrimaryVertices(): pv_size = " << m_primaryVertices->size() );
833 for(
const auto *
vertex : *m_primaryVertices ) {
844 if( m_jp.FillNtuple ) {
846 if( 0 == m_ntupleVars->get<
unsigned int>(
"NumPV" ) ) {
848 m_ntupleVars->get<
double>(
"PVX" ) =
vertex->x();
849 m_ntupleVars->get<
double>(
"PVY" ) =
vertex->y();
850 m_ntupleVars->get<
double>(
"PVZ" ) =
vertex->z();
851 m_ntupleVars->get<
unsigned int>(
"PVType" ) =
vertex->vertexType();
854 m_ntupleVars->get<
unsigned int>(
"NTrksPV" ) =
vertex->nTrackParticles();
857 m_ntupleVars->get<
unsigned int>(
"NumPV" )++;
859 m_ntupleVars->get< vector<int> > (
"NdofTrksPV" ) .emplace_back(
vertex->numberDoF() );
860 m_ntupleVars->get< vector<double> >(
"PVZpile" ) .emplace_back(
vertex->position().z() );
873 ATH_MSG_DEBUG(
"No Reconstructed PV was found. Using the dummy PV instead.");
874 for(
const auto *
vertex : *m_primaryVertices ) {
877 if( m_jp.FillNtuple ) {
879 if( 0 == m_ntupleVars->get<
unsigned int>(
"NumPV" ) ) {
882 m_ntupleVars->get<
double>(
"PVX" ) =
vertex->x();
883 m_ntupleVars->get<
double>(
"PVY" ) =
vertex->y();
884 m_ntupleVars->get<
double>(
"PVZ" ) =
vertex->z();
885 m_ntupleVars->get<
unsigned int>(
"PVType" ) =
vertex->vertexType();
888 m_ntupleVars->get<
unsigned int>(
"NTrksPV" ) =
vertex->nTrackParticles();
897 return StatusCode::FAILURE;
900 ATH_MSG_DEBUG(
" Primary vertex successful. thePV = " << m_thePV );
902 return StatusCode::SUCCESS;
907 void VrtSecInclusive::trackClassification(std::vector<WrkVrt> *workVerticesContainer, std::map<
long int, std::vector<long int> >& trackToVertexMap )
911 trackToVertexMap.clear();
913 for(
size_t iv = 0; iv<workVerticesContainer->size(); iv++ ) {
917 auto& trackIndices =
vertex.selectedTrackIndices;
918 if( !
vertex.isGood )
continue;
919 if( trackIndices.size() < 2 )
continue;
921 for(
auto&
index : trackIndices ) {
922 trackToVertexMap[
index].emplace_back( iv );
926 for(
auto& pair: trackToVertexMap ) {
927 std::string
msg = Form(
"track index [%ld]: vertices = (", pair.first);
928 for(
auto&
v : pair.second ) {
929 msg += Form(
"%ld, ",
v);
939 double VrtSecInclusive::findWorstChi2ofMaximallySharedTrack(std::vector<WrkVrt> *workVerticesContainer,
940 std::map<
long int, std::vector<long int> >& trackToVertexMap,
941 long int & maxSharedTrack,
942 long int & worstMatchingVertex)
948 auto maxSharedTrackToVertices = std::max_element( trackToVertexMap.begin(), trackToVertexMap.end(), [](
auto&
p1,
auto&
p2 ) { return p1.second.size() < p2.second.size(); } );
950 if( maxSharedTrackToVertices == trackToVertexMap.end() )
return worstChi2;
952 ATH_MSG_VERBOSE(
" > " << __FUNCTION__ <<
": max-shared track index = " << maxSharedTrackToVertices->first <<
", number of shared vertices = " << maxSharedTrackToVertices->second.size() );
954 if( maxSharedTrackToVertices->second.size() < 2 )
return worstChi2;
957 std::map<long int, double> vrtChi2Map;
960 for(
auto& iv : maxSharedTrackToVertices->second ) {
961 ATH_MSG_VERBOSE(
" > " << __FUNCTION__ <<
": loop over vertices: vertex index " << iv );
963 auto& wrkvrt = workVerticesContainer->at( iv );
964 auto& trackIndices = wrkvrt.selectedTrackIndices;
967 auto index = std::find_if( trackIndices.begin(), trackIndices.end(), [&](
auto&
index ) { return index == maxSharedTrackToVertices->first; } );
968 if(
index == trackIndices.end() ) {
969 ATH_MSG_WARNING(
" >> " << __FUNCTION__ <<
": index not found (algorithm inconsistent)" );
973 auto&
chi2 = wrkvrt.Chi2PerTrk.at(
index - trackIndices.begin() );
975 vrtChi2Map.emplace( std::pair<long int, double>(iv,
chi2) );
978 auto worstVrtChi2Pair = std::max_element( vrtChi2Map.begin(), vrtChi2Map.end(), [](
auto&
p1,
auto&
p2 ) { return p1.second < p2.second; } );
980 if( worstVrtChi2Pair == vrtChi2Map.end() ) {
981 ATH_MSG_WARNING(
" >> " << __FUNCTION__ <<
": max_element of vrtChi2Map not found" );
985 maxSharedTrack = maxSharedTrackToVertices->first;
986 worstMatchingVertex = worstVrtChi2Pair->first;
987 worstChi2 = worstVrtChi2Pair->second;
994 void VrtSecInclusive::printWrkSet(
const std::vector<WrkVrt> *workVerticesContainer,
const std::string&
name)
996 ATH_MSG_DEBUG(
" >> " << __FUNCTION__ <<
": ===============================================================" );
997 ATH_MSG_DEBUG(
" >> " << __FUNCTION__ <<
": " <<
name <<
": #vertices = " << workVerticesContainer->size() );
999 std::set<const xAOD::TrackParticle*> usedTracks;
1001 auto concatenateIndicesToString = [](
auto indices,
auto& collection ) -> std::string {
1002 if( 0 ==
indices.size() )
return "";
1004 [&collection](
const std::string&
str,
auto&
index ) { return str +
", " + std::to_string( collection.at(index)->index() ); } );
1007 std::map<const xAOD::TruthVertex*, bool>
previous;
1009 for(
auto& pair : m_matchMap ) {
previous.emplace( pair.first, pair.second ); }
1012 for(
const auto* truthVertex : m_tracingTruthVertices ) { m_matchMap.emplace( truthVertex,
false ); }
1014 for(
size_t iv=0; iv<workVerticesContainer->size(); iv++) {
1015 const auto& wrkvrt = workVerticesContainer->at(iv);
1017 if( wrkvrt.nTracksTotal() < 2 )
continue;
1019 std::string sels = concatenateIndicesToString( wrkvrt.selectedTrackIndices, m_selectedTracks );
1020 std::string assocs = concatenateIndicesToString( wrkvrt.associatedTrackIndices, m_associatedTracks );
1022 for(
const auto&
index : wrkvrt.selectedTrackIndices ) { usedTracks.insert( m_selectedTracks.at(
index) ); }
1023 for(
const auto&
index : wrkvrt.associatedTrackIndices ) { usedTracks.insert( m_associatedTracks.at(
index) ); }
1025 ATH_MSG_DEBUG(
" >> " << __FUNCTION__ <<
": " <<
name <<
" vertex [" << iv <<
"]: " << &wrkvrt
1026 <<
", isGood = " << (wrkvrt.isGood?
"true" :
"false")
1027 <<
", #ntrks(tot, sel, assoc) = (" << wrkvrt.nTracksTotal() <<
", " << wrkvrt.selectedTrackIndices.size() <<
", " << wrkvrt.associatedTrackIndices.size() <<
"), "
1028 <<
", chi2/ndof = " << wrkvrt.fitQuality()
1029 <<
", (r, z) = (" << wrkvrt.vertex.perp()
1030 <<
", " << wrkvrt.vertex.z() <<
")"
1031 <<
", sels = { " << sels <<
" }"
1032 <<
", assocs = { " << assocs <<
" }" );
1035 for(
const auto* truthVertex : m_tracingTruthVertices ) {
1038 Amg::Vector3D vTruth( truthVertex->x(), truthVertex->y(), truthVertex->z() );
1039 Amg::Vector3D vReco ( wrkvrt.vertex.x(), wrkvrt.vertex.y(), wrkvrt.vertex.z() );
1041 const auto distance = vReco - vTruth;
1044 cov.fillSymmetric( 0, 0, wrkvrt.vertexCov.at(0) );
1045 cov.fillSymmetric( 1, 0, wrkvrt.vertexCov.at(1) );
1046 cov.fillSymmetric( 1, 1, wrkvrt.vertexCov.at(2) );
1047 cov.fillSymmetric( 2, 0, wrkvrt.vertexCov.at(3) );
1048 cov.fillSymmetric( 2, 1, wrkvrt.vertexCov.at(4) );
1049 cov.fillSymmetric( 2, 2, wrkvrt.vertexCov.at(5) );
1053 if(
distance.norm() < 2.0 ||
s2 < 100. ) m_matchMap.at( truthVertex ) =
true;
1059 ATH_MSG_DEBUG(
" >> " << __FUNCTION__ <<
": number of used tracks = " << usedTracks.size() );
1062 for(
auto& pair : m_matchMap ) {
1064 if( pair.second !=
previous.at( pair.first ) ) {
1065 ATH_MSG_DEBUG(
" >> " << __FUNCTION__ <<
": Match flag has changed: (r, z) = (" << pair.first->perp() <<
", " << pair.first->z() <<
")" );
1070 if( m_jp.FillHist ) {
1071 for(
auto& pair : m_matchMap ) {
1072 if( pair.second ) m_hists[
"nMatchedTruths"]->Fill( m_vertexingAlgorithmStep+2, pair.first->perp() );
1077 for(
const auto* trk : usedTracks ) {
msg += Form(
"%ld, ", trk->index() ); }
1080 ATH_MSG_DEBUG(
" >> " << __FUNCTION__ <<
": ===============================================================" );
1106 const EventContext& ctx = Gaudi::Hive::currentContext();
1107 std::vector<std::unique_ptr<Trk::TrackParameters>> paramsVector =
1114 for(
auto&
params : paramsVector ) {
1116 const TVector3 position(
params->position().x(),
params->position().y(),
params->position().z() );
1118 if( prevPos == position ) {
1124 const auto* detElement =
params->associatedSurface().associatedDetectorElement();
1130 const auto&
id = detElement->identify();
1133 if( m_atlasId->is_pixel(
id) ) {
1135 auto idHash = m_pixelId->wafer_hash(
id );
1136 good = m_pixelCondSummaryTool->isGood( idHash, ctx );
1138 pattern->emplace_back( std::make_tuple( position,
Pixel, m_pixelId->barrel_ec(
id), m_pixelId->layer_disk(
id),
good ) );
1140 }
else if( m_atlasId->is_sct(
id) ) {
1142 auto idHash = m_sctId->wafer_hash(
id );
1143 good = m_sctCondSummaryTool->isGood( idHash, ctx );
1145 pattern->emplace_back( std::make_tuple( position,
SCT, m_sctId->barrel_ec(
id), m_sctId->layer_disk(
id),
good ) );
1151 ATH_MSG_VERBOSE(
" >> " << __FUNCTION__ <<
", track " << trk <<
": position = (" << position.Perp() <<
", " << position.z() <<
", " << position.Phi() <<
"), detElement ID = " <<
id <<
", good = " <<
good
1152 <<
": (det, bec, layer) = (" << std::get<1>(
pattern->back() ) <<
", " << std::get<2>(
pattern->back() ) <<
", " << std::get<3>(
pattern->back() ) <<
")" );
1154 if( !
good ) nDisabled++;
1161 if( m_jp.FillHist ) {
1162 m_hists[
"disabledCount"]->Fill( nDisabled );
1174 if( m_extrapolatedPatternBank.find( trk ) == m_extrapolatedPatternBank.end() ) {
1176 std::unique_ptr<ExtrapolatedPattern> exPattern_along( extrapolatedPattern( trk,
Trk::alongMomentum ) );
1177 std::unique_ptr<ExtrapolatedPattern> exPattern_oppos( extrapolatedPattern( trk,
Trk::oppositeMomentum ) );
1179 m_extrapolatedPatternBank.emplace( trk, std::make_pair( std::move(exPattern_along), std::move(exPattern_oppos) ) );
1183 auto& exPattern = m_extrapolatedPatternBank.at( trk );
1185 using LayerCombination = std::vector<int>;
1187 std::map<LayerCombination, unsigned> layerMap;
1188 if( layerMap.empty() ) {
1231 const LayerCombination comb { std::get<detector>( point ), std::get<bec>( point ), std::get<layer>( point ) };
1233 for(
auto& pair : layerMap ) {
1234 if( pair.first == comb ) {
1242 enum { kShouldNotHaveHit, kShouldHaveHit, kMayHaveHit };
1248 auto& exPattern_along = *( exPattern.first );
1250 for(
auto itr = exPattern_along.begin(); itr != exPattern_along.end(); ++itr ) {
1251 if(
std::next( itr ) == exPattern_along.end() )
continue;
1253 const auto& point = *itr;
1254 const auto& nextPoint = *(
std::next( itr ) );
1256 ATH_MSG_VERBOSE(
" > " << __FUNCTION__ <<
": isGood = " << std::get<isGood>( point ) );
1258 const auto& thisPos = std::get<position>( point );
1259 const auto& nextPos = std::get<position>( nextPoint );
1261 auto sectionVector = nextPos - thisPos;
1265 const auto& detectorType = getDetectorType( point );
1267 ATH_MSG_VERBOSE(
" > " << __FUNCTION__ <<
": detType = " << detectorType );
1273 if( vertexVector.Mag() < 10. ) {
1274 expectedHitPattern.at( detectorType ) = kMayHaveHit;
1280 if( !
static_cast<bool>(std::get<isGood>( point )) ) {
1281 expectedHitPattern.at( detectorType ) = kMayHaveHit;
1290 if( sectionVector.Mag() == 0. )
continue;
1293 <<
": hitPos = (" << thisPos.Perp() <<
", " << thisPos.z() <<
", " << thisPos.Phi() <<
")"
1294 <<
", sectionVec = (" << sectionVector.Perp() <<
", " << sectionVector.z() <<
", " << sectionVector.Phi() <<
")"
1295 <<
", vertexVec = (" << vertexVector.Perp() <<
", " << vertexVector.z() <<
", " << vertexVector.Phi() <<
")"
1296 <<
", cos(s,v) = " << sectionVector * vertexVector / ( sectionVector.Mag() * vertexVector.Mag() +
AlgConsts::infinitesimal ) );
1298 if( sectionVector * vertexVector > 0. )
continue;
1300 if( minExpectedRadius > thisPos.Perp() ) minExpectedRadius = thisPos.Perp();
1304 expectedHitPattern.at( detectorType ) = kShouldHaveHit;
1308 auto& exPattern_oppos = *( exPattern.second );
1309 bool oppositeFlag =
false;
1311 for(
auto itr = exPattern_oppos.begin(); itr != exPattern_oppos.end(); ++itr ) {
1312 if(
std::next( itr ) == exPattern_oppos.end() )
continue;
1314 const auto& point = *itr;
1315 const auto& nextPoint = *(
std::next( itr ) );
1317 const auto& thisPos = std::get<position>( point );
1318 const auto& nextPos = std::get<position>( nextPoint );
1320 auto sectionVector = nextPos - thisPos;
1323 const auto& detectorType = getDetectorType( point );
1325 ATH_MSG_VERBOSE(
" > " << __FUNCTION__ <<
": detType = " << detectorType );
1328 <<
": hitPos = (" << thisPos.Perp() <<
", " << thisPos.z() <<
", " << thisPos.Phi() <<
")"
1329 <<
", vertex = (" <<
vertex.perp() <<
", " <<
vertex.z() <<
", " <<
vertex.phi() <<
")"
1330 <<
", cos(s,v) = " << sectionVector * vertexVector / ( sectionVector.Mag() * vertexVector.Mag() +
AlgConsts::infinitesimal ) );
1335 if( sectionVector * vertexVector < 0. ) {
1336 oppositeFlag =
true;
1345 std::string
msg =
"Expected hit pattern: ";
1347 msg += Form(
"%s", expectedHitPattern.at(
i) < kMayHaveHit? Form(
"%u", expectedHitPattern.at(
i)) :
"X" );
1351 msg =
"Recorded hit pattern: ";
1357 std::vector<unsigned> matchedLayers;
1360 const unsigned recordedPattern = ( (trk->
hitPattern()>>
i) & 1 );
1362 if( expectedHitPattern.at(
i) == kMayHaveHit ) {
1363 matchedLayers.emplace_back(
i );
1364 }
else if( expectedHitPattern.at(
i) == kShouldHaveHit ) {
1365 if( expectedHitPattern.at(
i) != recordedPattern ) {
1368 matchedLayers.emplace_back(
i );
1371 if( expectedHitPattern.at(
i) != recordedPattern ) {
1386 while( dphi > TMath::Pi() ) dphi -= TMath::TwoPi();
1387 while( dphi < -TMath::Pi() ) dphi += TMath::TwoPi();
1389 ATH_MSG_DEBUG(
" > " << __FUNCTION__ <<
": vtx phi = " <<
vertex.phi() <<
", track phi = " << trk->
phi() <<
", dphi = " << dphi
1390 <<
", oppositeFlag = " << oppositeFlag
1391 <<
", nPixelHits = " <<
static_cast<int>(PixelHits)
1392 <<
", nSCTHits = " <<
static_cast<int>(SctHits)
1393 <<
", nTRTHits = " <<
static_cast<int>(TRTHits)
1394 <<
", nMatchedLayers = " << matchedLayers.size() );
1396 if( PixelHits == 0 &&
vertex.perp() > 300. ) {
1397 ATH_MSG_DEBUG(
" > " << __FUNCTION__ <<
": vertex r > 300 mm, w/o no pixel hits" );
1402 if( matchedLayers.size() < 2 )
return false;
1406 if( matchedLayers.size() < 4 )
return false;
1412 if( oppositeFlag )
return false;
1418 if( fabs( dphi ) > TMath::Pi()/2.0 )
return false;
1421 TVector3 trkP; trkP.SetPtEtaPhi( trk->
pt(), trk->
eta(), trk->
phi() );
1423 if( trkP.Dot( vtx ) < 0. )
return false;
1443 const double absz = fabs(
vertex.z() );
1452 outsidePixelBarrel0_and_insidePixelBarrel1,
1455 outsidePixelBarrel1_and_insidePixelBarrel2,
1458 outsidePixelBarrel2_and_insidePixelBarrel3,
1461 outsidePixelBarrel3_and_insideSctBarrel0,
1464 outsideSctBarrel0_and_insideSctBarrel1,
1469 int vertex_pattern = 0;
1471 vertex_pattern = insideBeamPipe;
1473 }
else if(
rad < 31.0 && absz < 331.5 ) {
1474 vertex_pattern = insidePixelBarrel0;
1476 }
else if(
rad < 38.4 && absz < 331.5 ) {
1477 vertex_pattern = aroundPixelBarrel0;
1479 }
else if(
rad < 47.7 && absz < 400.5 ) {
1480 vertex_pattern = outsidePixelBarrel0_and_insidePixelBarrel1;
1482 }
else if(
rad < 54.4 && absz < 400.5 ) {
1483 vertex_pattern = aroundPixelBarrel1;
1485 }
else if(
rad < 85.5 && absz < 400.5 ) {
1486 vertex_pattern = outsidePixelBarrel1_and_insidePixelBarrel2;
1488 }
else if(
rad < 92.2 && absz < 400.5 ) {
1489 vertex_pattern = aroundPixelBarrel2;
1491 }
else if(
rad < 119.3 && absz < 400.5 ) {
1492 vertex_pattern = outsidePixelBarrel2_and_insidePixelBarrel3;
1494 }
else if(
rad < 126.1 && absz < 400.5 ) {
1495 vertex_pattern = aroundPixelBarrel3;
1497 }
else if(
rad < 290 && absz < 749.0 ) {
1498 vertex_pattern = outsidePixelBarrel3_and_insideSctBarrel0;
1500 }
else if(
rad < 315 && absz < 749.0 ) {
1501 vertex_pattern = aroundSctBarrel0;
1503 }
else if(
rad < 360 && absz < 749.0 ) {
1504 vertex_pattern = outsideSctBarrel0_and_insideSctBarrel1;
1506 }
else if(
rad < 390 && absz < 749.0 ) {
1507 vertex_pattern = aroundSctBarrel1;
1512 unsigned nPixelLayers { 0 };
1524 if( vertex_pattern == insideBeamPipe ) {
1527 if( nPixelLayers < 3 )
return false;
1530 }
else if( vertex_pattern == insidePixelBarrel0 ) {
1533 if( nPixelLayers < 3 )
return false;
1537 else if( vertex_pattern == aroundPixelBarrel0 ) {
1541 if( nPixelLayers < 2 )
return false;
1545 else if( vertex_pattern == outsidePixelBarrel0_and_insidePixelBarrel1 ) {
1549 if( nPixelLayers < 2 )
return false;
1553 else if( vertex_pattern == aroundPixelBarrel1 ) {
1558 if( nPixelLayers < 2 )
return false;
1562 else if( vertex_pattern == outsidePixelBarrel1_and_insidePixelBarrel2 ) {
1567 if( nPixelLayers < 2 )
return false;
1571 else if( vertex_pattern == aroundPixelBarrel2 ) {
1580 else if( vertex_pattern == outsidePixelBarrel2_and_insidePixelBarrel3 ) {
1588 else if( vertex_pattern == aroundPixelBarrel3 ) {
1598 else if( vertex_pattern == outsidePixelBarrel3_and_insideSctBarrel0 ) {
1608 else if( vertex_pattern == aroundSctBarrel0 ) {
1619 else if( vertex_pattern == outsideSctBarrel0_and_insideSctBarrel1 ) {
1630 else if( vertex_pattern == aroundSctBarrel1 ) {
1657 const double absz = fabs(
vertex.z() );
1666 outsidePixelBarrel0_and_insidePixelBarrel1,
1669 outsidePixelBarrel1_and_insidePixelBarrel2,
1672 outsidePixelBarrel2_and_insidePixelBarrel3,
1675 outsidePixelBarrel3_and_insideSctBarrel0,
1678 outsideSctBarrel0_and_insideSctBarrel1,
1683 int vertex_pattern = 0;
1685 vertex_pattern = insideBeamPipe;
1687 }
else if(
rad < 31.0 && absz < 331.5 ) {
1688 vertex_pattern = insidePixelBarrel0;
1690 }
else if(
rad < 38.4 && absz < 331.5 ) {
1691 vertex_pattern = aroundPixelBarrel0;
1693 }
else if(
rad < 47.7 && absz < 400.5 ) {
1694 vertex_pattern = outsidePixelBarrel0_and_insidePixelBarrel1;
1696 }
else if(
rad < 54.4 && absz < 400.5 ) {
1697 vertex_pattern = aroundPixelBarrel1;
1699 }
else if(
rad < 85.5 && absz < 400.5 ) {
1700 vertex_pattern = outsidePixelBarrel1_and_insidePixelBarrel2;
1702 }
else if(
rad < 92.2 && absz < 400.5 ) {
1703 vertex_pattern = aroundPixelBarrel2;
1705 }
else if(
rad < 119.3 && absz < 400.5 ) {
1706 vertex_pattern = outsidePixelBarrel2_and_insidePixelBarrel3;
1708 }
else if(
rad < 126.1 && absz < 400.5 ) {
1709 vertex_pattern = aroundPixelBarrel3;
1711 }
else if(
rad < 290 && absz < 749.0 ) {
1712 vertex_pattern = outsidePixelBarrel3_and_insideSctBarrel0;
1714 }
else if(
rad < 315 && absz < 749.0 ) {
1715 vertex_pattern = aroundSctBarrel0;
1717 }
else if(
rad < 360 && absz < 749.0 ) {
1718 vertex_pattern = outsideSctBarrel0_and_insideSctBarrel1;
1720 }
else if(
rad < 390 && absz < 749.0 ) {
1721 vertex_pattern = aroundSctBarrel1;
1727 unsigned nPixelLayers { 0 };
1739 if( vertex_pattern == insideBeamPipe ) {
1743 if( nPixelLayers < 3 )
return false;
1745 }
else if( vertex_pattern == insidePixelBarrel0 ) {
1749 if( nPixelLayers < 3 )
return false;
1754 else if( vertex_pattern == aroundPixelBarrel0 ) {
1759 if( nPixelLayers < 3 )
return false;
1763 else if( vertex_pattern == outsidePixelBarrel0_and_insidePixelBarrel1 ) {
1767 if( nPixelLayers < 3 )
return false;
1771 else if( vertex_pattern == aroundPixelBarrel1 ) {
1776 if( nPixelLayers < 2 )
return false;
1780 else if( vertex_pattern == outsidePixelBarrel1_and_insidePixelBarrel2 ) {
1784 if( nPixelLayers < 2 )
return false;
1788 else if( vertex_pattern == aroundPixelBarrel2 ) {
1795 else if( vertex_pattern == outsidePixelBarrel2_and_insidePixelBarrel3 ) {
1800 else if( vertex_pattern == aroundPixelBarrel3 ) {
1808 else if( vertex_pattern == outsidePixelBarrel3_and_insideSctBarrel0 ) {
1815 else if( vertex_pattern == aroundSctBarrel0 ) {
1823 else if( vertex_pattern == outsideSctBarrel0_and_insideSctBarrel1 ) {
1830 else if( vertex_pattern == aroundSctBarrel1 ) {
1851 const double absz = fabs(
vertex.z() );
1860 outsidePixelBarrel1_and_insidePixelBarrel2,
1863 outsidePixelBarrel2_and_insidePixelBarrel3,
1866 outsidePixelBarrel3_and_insideSctBarrel0,
1869 outsideSctBarrel0_and_insideSctBarrel1,
1874 Int_t vertex_pattern = 0;
1876 vertex_pattern = insideBeamPipe;
1878 }
else if(
rad < 47.7 && absz < 400.5 ) {
1879 vertex_pattern = insidePixelBarrel1;
1881 }
else if(
rad < 54.4 && absz < 400.5 ) {
1882 vertex_pattern = aroundPixelBarrel1;
1884 }
else if(
rad < 85.5 && absz < 400.5 ) {
1885 vertex_pattern = outsidePixelBarrel1_and_insidePixelBarrel2;
1887 }
else if(
rad < 92.2 && absz < 400.5 ) {
1888 vertex_pattern = aroundPixelBarrel2;
1890 }
else if(
rad < 119.3 && absz < 400.5 ) {
1891 vertex_pattern = outsidePixelBarrel2_and_insidePixelBarrel3;
1893 }
else if(
rad < 126.1 && absz < 400.5 ) {
1894 vertex_pattern = aroundPixelBarrel3;
1896 }
else if(
rad < 290 && absz < 749.0 ) {
1897 vertex_pattern = outsidePixelBarrel3_and_insideSctBarrel0;
1899 }
else if(
rad < 315 && absz < 749.0 ) {
1900 vertex_pattern = aroundSctBarrel0;
1902 }
else if(
rad < 360 && absz < 749.0 ) {
1903 vertex_pattern = outsideSctBarrel0_and_insideSctBarrel1;
1905 }
else if(
rad < 390 && absz < 749.0 ) {
1906 vertex_pattern = aroundSctBarrel1;
1913 if( vertex_pattern == insideBeamPipe ) {
1920 else if( vertex_pattern == insidePixelBarrel1 ) {
1926 else if( vertex_pattern == aroundPixelBarrel1 ) {
1933 else if( vertex_pattern == outsidePixelBarrel1_and_insidePixelBarrel2 ) {
1940 else if( vertex_pattern == aroundPixelBarrel2 ) {
1948 else if( vertex_pattern == outsidePixelBarrel2_and_insidePixelBarrel3 ) {
1955 else if( vertex_pattern == aroundPixelBarrel3 ) {
1964 else if( vertex_pattern == outsidePixelBarrel3_and_insideSctBarrel0 ) {
1973 else if( vertex_pattern == aroundSctBarrel0 ) {
1983 else if( vertex_pattern == outsideSctBarrel0_and_insideSctBarrel1 ) {
1993 else if( vertex_pattern == aroundSctBarrel1 ) {
2016 const double absz = fabs(
vertex.z() );
2025 outsidePixelBarrel1_and_insidePixelBarrel2,
2028 outsidePixelBarrel2_and_insidePixelBarrel3,
2031 outsidePixelBarrel3_and_insideSctBarrel0,
2034 outsideSctBarrel0_and_insideSctBarrel1,
2039 Int_t vertex_pattern = 0;
2041 vertex_pattern = insideBeamPipe;
2043 }
else if(
rad < 47.7 && absz < 400.5 ) {
2044 vertex_pattern = insidePixelBarrel1;
2046 }
else if(
rad < 54.4 && absz < 400.5 ) {
2047 vertex_pattern = aroundPixelBarrel1;
2049 }
else if(
rad < 85.5 && absz < 400.5 ) {
2050 vertex_pattern = outsidePixelBarrel1_and_insidePixelBarrel2;
2052 }
else if(
rad < 92.2 && absz < 400.5 ) {
2053 vertex_pattern = aroundPixelBarrel2;
2055 }
else if(
rad < 119.3 && absz < 400.5 ) {
2056 vertex_pattern = outsidePixelBarrel2_and_insidePixelBarrel3;
2058 }
else if(
rad < 126.1 && absz < 400.5 ) {
2059 vertex_pattern = aroundPixelBarrel3;
2061 }
else if(
rad < 290 && absz < 749.0 ) {
2062 vertex_pattern = outsidePixelBarrel3_and_insideSctBarrel0;
2064 }
else if(
rad < 315 && absz < 749.0 ) {
2065 vertex_pattern = aroundSctBarrel0;
2067 }
else if(
rad < 360 && absz < 749.0 ) {
2068 vertex_pattern = outsideSctBarrel0_and_insideSctBarrel1;
2070 }
else if(
rad < 390 && absz < 749.0 ) {
2071 vertex_pattern = aroundSctBarrel1;
2078 if( vertex_pattern == insideBeamPipe ) {
2085 else if( vertex_pattern == insidePixelBarrel1 ) {
2091 else if( vertex_pattern == aroundPixelBarrel1 ) {
2098 else if( vertex_pattern == outsidePixelBarrel1_and_insidePixelBarrel2 ) {
2104 else if( vertex_pattern == aroundPixelBarrel2 ) {
2111 else if( vertex_pattern == outsidePixelBarrel2_and_insidePixelBarrel3 ) {
2116 else if( vertex_pattern == aroundPixelBarrel3 ) {
2123 else if( vertex_pattern == outsidePixelBarrel3_and_insideSctBarrel0 ) {
2129 else if( vertex_pattern == aroundSctBarrel0 ) {
2136 else if( vertex_pattern == outsideSctBarrel0_and_insideSctBarrel1 ) {
2142 else if( vertex_pattern == aroundSctBarrel1 ) {
2203 if( m_extrapolatedPatternBank.find( trk ) == m_extrapolatedPatternBank.end() ) {
2205 std::unique_ptr<ExtrapolatedPattern> exPattern_along( extrapolatedPattern( trk,
Trk::alongMomentum ) );
2207 m_extrapolatedPatternBank.emplace( trk, std::make_pair( std::move(exPattern_along),
nullptr ) );
2211 if(
vertex.perp() < 31.0 ) {
2212 double dphi = trk->
phi() -
vertex.phi();
2213 while( dphi > TMath::Pi() ) { dphi -= TMath::TwoPi(); }
2214 while( dphi < -TMath::Pi() ) { dphi += TMath::TwoPi(); }
2215 if(
cos(dphi) < -0.8 )
return false;
2218 auto& exPattern = m_extrapolatedPatternBank.at( trk );
2220 using LayerCombination = std::vector<int>;
2222 std::map<LayerCombination, unsigned> layerMap;
2223 if( layerMap.empty() ) {
2266 const LayerCombination comb { std::get<detector>( point ), std::get<bec>( point ), std::get<layer>( point ) };
2268 for(
auto& pair : layerMap ) {
2269 if( pair.first == comb ) {
2280 auto& exPattern_along = *( exPattern.first );
2282 for(
auto itr = exPattern_along.begin(); itr != exPattern_along.end(); ++itr ) {
2283 if(
std::next( itr ) == exPattern_along.end() )
continue;
2285 const auto& point = *itr;
2287 ATH_MSG_VERBOSE(
" > " << __FUNCTION__ <<
": isGood = " << std::get<isGood>( point ) );
2289 if( !std::get<isGood>( point ) ) {
2290 const auto& detectorType = getDetectorType( point );
2291 disabledPattern += (1 << detectorType);
2296 uint32_t modifiedPattern = disabledPattern | hitPattern;
2298 std::string
msg =
"Disabled hit pattern: ";
2300 msg += Form(
"%u", ( disabledPattern >>
i ) & 1 );
2304 msg =
"Recorded hit pattern: ";
2306 msg += Form(
"%u", ( hitPattern >>
i ) & 1 );
2310 return patternCheck( modifiedPattern,
vertex );
2321 const bool& check_itrk = ( this->*m_patternStrategyFuncs[m_checkPatternStrategy] )( itrk, FitVertex );
2322 const bool& check_jtrk = ( this->*m_patternStrategyFuncs[m_checkPatternStrategy] )( jtrk, FitVertex );
2324 return ( check_itrk && check_jtrk );
2330 void VrtSecInclusive::removeInconsistentTracks(
WrkVrt& wrkvrt )
2335 std::map< std::deque<long int>*, std::vector<const xAOD::TrackParticle*>& > indexMap
2340 for(
auto& pair : indexMap ) {
2343 auto& tracks = pair.second;
2346 [&](
auto&
index ) {
2347 bool isConsistent = (this->*m_patternStrategyFuncs[m_checkPatternStrategy] )( tracks.at(index), vertex );
2348 return !isConsistent;
2359 void VrtSecInclusive::dumpTruthInformation() {
2365 auto sc0 = evtStore()->retrieve( eventInfo,
"EventInfo" );
2366 if( sc0.isFailure() ) {
return; }
2372 auto sc1 = evtStore()->retrieve( truthParticles,
"TruthParticles" );
2373 if( sc1.isFailure() ) {
return; }
2375 auto sc2 = evtStore()->retrieve( truthVertices,
"TruthVertices" );
2376 if( sc2.isFailure() ) {
return; }
2378 if( !truthParticles ) {
return; }
2379 if( !truthVertices ) {
return; }
2381 m_tracingTruthVertices.
clear();
2383 std::vector<const xAOD::TruthParticle*> truthSvTracks;
2390 if( truthVertex->nIncomingParticles() != 1 )
return false;
2391 if( !truthVertex->incomingParticle(0) )
return false;
2392 if( abs(truthVertex->incomingParticle(0)->pdgId()) < 1000000 )
return false;
2393 if( abs(truthVertex->incomingParticle(0)->pdgId()) > 1000000000 )
return false;
2395 bool hasNeutralino =
false;
2396 for(
unsigned ip = 0;
ip < truthVertex->nOutgoingParticles();
ip++ ) {
2397 const auto*
p = truthVertex->outgoingParticle(
ip);
2398 if( abs(
p->pdgId() ) == 1000022 ) {
2399 hasNeutralino =
true;
2403 return hasNeutralino;
2407 if( truthVertex->nIncomingParticles() != 1 )
return false;
2408 if( !truthVertex->incomingParticle(0) )
return false;
2409 if( abs(truthVertex->incomingParticle(0)->pdgId()) != 50 )
return false;
2414 if( truthVertex->nIncomingParticles() != 1 )
return false;
2415 if( !truthVertex->incomingParticle(0) )
return false;
2416 if( abs(truthVertex->incomingParticle(0)->pdgId()) != 36 )
return false;
2421 if( truthVertex->nIncomingParticles() != 1 )
return false;
2422 if( !truthVertex->incomingParticle(0) )
return false;
2423 if( abs(truthVertex->incomingParticle(0)->pdgId()) != 310 )
return false;
2428 if( truthVertex->nIncomingParticles() != 1 )
return false;
2429 if( !truthVertex->incomingParticle(0) )
return false;
2430 if( abs(truthVertex->incomingParticle(0)->pdgId()) <= 500 || abs(truthVertex->incomingParticle(0)->pdgId()) >= 600 )
return false;
2435 if( truthVertex->nIncomingParticles() != 1 )
return false;
2436 if( !truthVertex->incomingParticle(0) )
return false;
2438 const auto*
parent = truthVertex->incomingParticle(0);
2439 if(
parent->isLepton() )
return false;
2441 TLorentzVector p4sum_in;
2442 TLorentzVector p4sum_out;
2443 for(
unsigned ip = 0;
ip < truthVertex->nIncomingParticles();
ip++ ) {
2444 const auto*
particle = truthVertex->incomingParticle(
ip);
2448 for(
unsigned ip = 0;
ip < truthVertex->nOutgoingParticles();
ip++ ) {
2449 const auto*
particle = truthVertex->outgoingParticle(
ip);
2453 return p4sum_out.E() - p4sum_in.E() >= 100.;
2459 std::map<std::string, ParticleSelectFunc> selectFuncs { {
"", selectNone },
2460 {
"Kshort", selectKshort },
2461 {
"Bhadron", selectBhadron },
2462 {
"Rhadron", selectRhadron },
2463 {
"HNL", selectHNL },
2464 {
"Higgs", selectHiggs },
2465 {
"HadInt", selectHadInt } };
2468 if( selectFuncs.find( m_jp.truthParticleFilter ) == selectFuncs.end() ) {
2469 ATH_MSG_WARNING(
" > " << __FUNCTION__ <<
": invalid function specification: " << m_jp.truthParticleFilter );
2473 auto selectFunc = selectFuncs.at( m_jp.truthParticleFilter );
2476 for(
const auto *truthVertex : *truthVertices ) {
2477 if( selectFunc( truthVertex ) ) {
2478 m_tracingTruthVertices.emplace_back( truthVertex );
2480 msg += Form(
"pdgId = %d, (r, z) = (%.2f, %.2f), ", truthVertex->incomingParticle(0)->pdgId(), truthVertex->perp(), truthVertex->z());
2481 msg += Form(
"nOutgoing = %lu, ", truthVertex->nOutgoingParticles() );
2482 msg += Form(
"mass = %.3f GeV, pt = %.3f GeV", truthVertex->incomingParticle(0)->m()/1.e3, truthVertex->incomingParticle(0)->pt()/1.e3 );
2487 if( m_jp.FillHist ) {
2488 for(
const auto* truthVertex : m_tracingTruthVertices ) {
2489 m_hists[
"nMatchedTruths"]->Fill( 0., truthVertex->perp() );