48 const double probVrtMergeLimit = 0.01;
52 inpNPart = xAODwrk->InpTrk.size();
53 xAODwrk->FoundSecondTracks.clear();
56 ATH_MSG_DEBUG(
"InDet getVrtSecMulti() called with NPart=" << inpNPart );
58 std::vector<xAOD::Vertex*> finalVertices(0);
65 const EventContext & ctx = Gaudi::Hive::currentContext();
67 if (eventInfo.isValid()) {
68 if(eventInfo->hasBeamSpotWeight()) evtWgt *= eventInfo->beamSpotWeight();
73 TLorentzVector momentumJet;
76 nRefPVTrk =
selGoodTrkParticle(xAODwrk->InpTrk, primVrt, jetDir, xAODwrk->listJetTracks,evtWgt);
77 while (!xAODwrk->listJetTracks.empty() &&
78 xAODwrk->listJetTracks[0]->pt() / jetDir.Pt() > 1.){
79 xAODwrk->listJetTracks.erase(xAODwrk->listJetTracks.begin());
81 nTracks = xAODwrk->listJetTracks.size();
82 momentumJet =
totalMom(xAODwrk->listJetTracks);
90 ATH_MSG_DEBUG(
"Number of selected tracks inside jet= " << nTracks);
94 h.m_hb_jmom->Fill(momentumJet.Perp(), evtWgt);
95 h.m_hb_ntrkjet->Fill(
static_cast<double>(nTracks), evtWgt);
106 std::vector<double> inpMass(nTracks,
m_massPi);
107 double vrt2TrackNumber = 0;
111 xAODwrk->TracksForFit,
117 xAODwrk->listSecondTracks,
119 if (xAODwrk->TrkFromV0.size() > 1) {
123 vrt2TrackNumber = xAODwrk->listSecondTracks.size()/2.0;
129 ATH_MSG_DEBUG(
" Found different tracks in pairs=" << vrt2TrackNumber);
141 std::unique_ptr<std::vector<WrkVrt>> wrkVrtSet = std::make_unique<std::vector<WrkVrt>>();
149 std::vector<std::vector<int>> allCliques;
150 bron_kerbosch_all_cliques(compatibilityGraph, clique_visitor(allCliques));
152 for (
const auto& clique : allCliques) {
154 newvrt.selTrk.clear();
155 if (xAODwrk) xAODwrk->tmpListTracks.clear();
157 for (
int i_trk : clique){
158 newvrt.selTrk.push_back(i_trk);
159 if(xAODwrk) xAODwrk->tmpListTracks.push_back(xAODwrk->listJetTracks.at(i_trk));
163 if (
sc.isFailure() ||
174 double jetVrtDir = jetDir.Px() * vDist.x() + jetDir.Py() * vDist.y() +
175 jetDir.Pz() * vDist.z();
177 if (jetVrtDir > 0.) {
188 sc = StatusCode::FAILURE;
201 if (
sc.isFailure())
continue;
203 if (clique.size() == 2 && newvrt.chi2 > 10.)
continue;
205 if (newvrt.chi2PerTrk.size() == 2){
206 newvrt.chi2PerTrk[0] = newvrt.chi2PerTrk[1] = newvrt.chi2 / 2.;
209 newvrt.nCloseVrt = 0;
210 newvrt.dCloseVrt = 1000000.;
211 newvrt.projectedVrt =
213 wrkVrtSet->push_back(newvrt);
219 bool disassembled =
false;
222 disassembled =
false;
223 int NSoluI = (*wrkVrtSet).size();
224 for (
int iv = 0; iv < NSoluI; iv++) {
225 WrkVrt vrt = (*wrkVrtSet)[iv];
226 if (!vrt.Good || vrt.selTrk.size() == 2)
continue;
227 if (TMath::Prob(vrt.chi2, 2 * vrt.selTrk.size() - 3) < 1.e-3) {
228 if (xAODwrk)
disassembleVertex(wrkVrtSet.get(), iv, xAODwrk->listJetTracks, *state);
232 }
while(disassembled);
235 while( (*wrkVrtSet).size()>1 ){
236 int tmpN = (*wrkVrtSet).size();
239 for(; iv<tmpN-1; iv++){
240 WrkVrt vert_i = (*wrkVrtSet)[iv];
241 if(!vert_i.Good )
continue;
242 int ntrk_i = (*wrkVrtSet)[iv].selTrk.size();
245 for(; jv<tmpN; jv++){
246 WrkVrt vert_j = (*wrkVrtSet)[jv];
247 if(!vert_j.Good )
continue;
248 int ntrk_j = (*wrkVrtSet)[jv].selTrk.size();
250 int nTCom =
nTrkCommon(wrkVrtSet.get(), iv, jv);
252 (*wrkVrtSet).erase((*wrkVrtSet).begin()+iv);
255 else if(nTCom==ntrk_j){
256 (*wrkVrtSet).erase((*wrkVrtSet).begin()+jv);
263 if(iv==tmpN-1)
break;
270 std::multimap<int,std::pair<int,int>> vrtWithCommonTrk;
271 std::multimap<int,std::pair<int,int>>::reverse_iterator icvrt;
273 int NSoluI = (*wrkVrtSet).size();
274 vrtWithCommonTrk.clear();
275 for(
int iv = 0; iv<NSoluI-1; iv++ ){
276 for(
int jv = iv+1; jv<NSoluI; jv++){
277 if(!(*wrkVrtSet)[iv].Good || !(*wrkVrtSet)[jv].Good)
continue;
278 int nTCom=
nTrkCommon(wrkVrtSet.get(), iv, jv);
280 vrtWithCommonTrk.emplace(nTCom,std::make_pair(iv,jv));
284 for(icvrt = vrtWithCommonTrk.rbegin(); icvrt!=vrtWithCommonTrk.rend(); ++icvrt){
285 int nTCom = (*icvrt).first;
286 int iv = (*icvrt).second.first;
287 int jv = (*icvrt).second.second;
288 int nTrkI = (*wrkVrtSet)[iv].selTrk.size();
289 int nTrkJ = (*wrkVrtSet)[jv].selTrk.size();
295 if(probV<probVrtMergeLimit){
296 if(nTrkI==2 || nTrkJ==2 || nTCom<2) {
299 if(nTCom>nTrkI-nTCom || nTCom>nTrkJ-nTCom){
308 (*wrkVrtSet).push_back(newvrt);
309 (*wrkVrtSet)[iv].Good =
false;
310 (*wrkVrtSet)[jv].Good =
false;
314 }
while( icvrt != vrtWithCommonTrk.rend() );
319 for(
const auto& vrt : (*wrkVrtSet)){
320 if(vrt.Good) cvgood++;
323 h.m_hb_rawVrtN->Fill(
static_cast<float>(cvgood), evtWgt);
327 for(
auto &tmpV : (*wrkVrtSet) ) {
328 if(tmpV.vertex.perp()>
m_rLayer3+10.) tmpV.Good =
false;
329 TLorentzVector SVPV(tmpV.vertex.x()-primVrt.
x(),
330 tmpV.vertex.y()-primVrt.
y(),
331 tmpV.vertex.z()-primVrt.
z(), 1.);
335 unsigned int tmpV = 0;
336 while( tmpV<(*wrkVrtSet).size() ){
337 if( !(*wrkVrtSet)[tmpV].Good ) (*wrkVrtSet).erase((*wrkVrtSet).begin()+tmpV);
341 if((*wrkVrtSet).empty()){
342 return finalVertices;
345 std::vector< std::vector<float> > trkScore(0);
347 for(
auto &trk : xAODwrk->listJetTracks) trkScore.push_back(
m_trackClassificator->trkTypeWgts(trk, primVrt, jetDir));
350 for(
auto &tmpV : (*wrkVrtSet) ) tmpV.projectedVrt=
jetProjDist(tmpV.vertex, primVrt, jetDir);
357 std::unique_ptr<std::vector< std::deque<long int> > > trkInVrt = std::make_unique<std::vector<std::deque<long int>>>(nTracks);
361 long int selectedTrack, selectedVertex;
362 int foundV1, foundV2;
365 while( ( foundMaxT =
MaxOfShared(wrkVrtSet.get(), trkInVrt.get(), selectedTrack, selectedVertex) ) > 0) {
367 double foundMinVrtDst = 1000000.;
372 &&
nTrkCommon(wrkVrtSet.get(), foundV1, foundV2) ){
374 bool vrtMerged=
false;
376 if(foundV1<foundV2)
std::swap(foundV1, foundV2);
379 if (xAODwrk) probV =
mergeAndRefitVertices(wrkVrtSet.get(), foundV1, foundV2, newvrt, xAODwrk->listJetTracks, *state);
383 double tstDst =
jetProjDist(newvrt.vertex, primVrt, jetDir);
386 std::swap((*wrkVrtSet)[foundV1], newvrt);
387 (*wrkVrtSet)[foundV1].projectedVrt = tstDst;
388 (*wrkVrtSet)[foundV2].Good =
false;
389 (*wrkVrtSet)[foundV2].selTrk.clear();
394 (*wrkVrtSet)[foundV1].nCloseVrt = -1;
395 (*wrkVrtSet)[foundV1].dCloseVrt = 1000000.;
396 (*wrkVrtSet)[foundV2].nCloseVrt = -1;
397 (*wrkVrtSet)[foundV2].dCloseVrt = 1000000.;
403 trkInVrt->resize(nTracks);
412 sc = StatusCode::FAILURE;
413 if(xAODwrk)
sc =
refitVertex( wrkVrtSet.get(), selectedVertex, xAODwrk->listJetTracks, *state,
false);
415 (*wrkVrtSet)[selectedVertex].projectedVrt =
jetProjDist((*wrkVrtSet)[selectedVertex].
vertex, primVrt, jetDir);
417 if(
sc.isFailure() ) (*wrkVrtSet)[selectedVertex].Good =
false;
418 if( (*wrkVrtSet)[selectedVertex].projectedVrt<0.
419 && (*wrkVrtSet)[selectedVertex].selTrk.size()==2 ){
420 (*wrkVrtSet)[selectedVertex].Good =
false;
429 double minDistVV =
minVrtVrtDist( wrkVrtSet.get(), foundV1, foundV2);
431 if(foundV1<foundV2)
std::swap(foundV1, foundV2);
434 if(xAODwrk) probV =
mergeAndRefitVertices(wrkVrtSet.get(), foundV1, foundV2, newvrt, xAODwrk->listJetTracks, *state);
438 double tstDst =
jetProjDist(newvrt.vertex, primVrt, jetDir);
442 (*wrkVrtSet)[foundV1].projectedVrt = tstDst;
443 (*wrkVrtSet)[foundV2].Good =
false;
444 (*wrkVrtSet)[foundV2].selTrk.clear();
448 (*wrkVrtSet)[foundV1].nCloseVrt = -1;
449 (*wrkVrtSet)[foundV1].dCloseVrt = 1000000.;
450 (*wrkVrtSet)[foundV2].nCloseVrt = -1;
451 (*wrkVrtSet)[foundV2].dCloseVrt = 1000000.;
457 for(
unsigned int iv=0; iv<wrkVrtSet->size(); iv++) {
458 WrkVrt vert_i = (*wrkVrtSet)[iv];
459 if(!vert_i.Good)
continue;
460 if( vert_i.selTrk.size()<3 )
continue;
462 double tmpProb = TMath::Prob( vert_i.chi2, 2*vert_i.selTrk.size()-3 );
464 if(xAODwrk) tmpProb =
improveVertexChi2(wrkVrtSet.get(), iv, xAODwrk->listJetTracks, *state,
false);
465 if(tmpProb<0.001) vert_i.Good =
false;
466 vert_i.projectedVrt =
jetProjDist(vert_i.vertex, primVrt, jetDir);
475 for(
auto &ntrVrt : (*wrkVrtSet)){
476 if(!ntrVrt.Good || ntrVrt.selTrk.size()<=1)
continue;
478 for(
auto &onetVrt : (*wrkVrtSet)){
479 if(!onetVrt.Good || onetVrt.selTrk.size()!=1)
continue;
481 if(ntrVrt.detachedTrack==onetVrt.selTrk[0]){
483 newV.selTrk.push_back(ntrVrt.detachedTrack);
485 if(xAODwrk) vProb =
refitVertex(newV, xAODwrk->listJetTracks, *state,
true);
486 if(vProb>probVrtMergeLimit){
489 ntrVrt.detachedTrack=-1;
497 for(
auto &iVrt : (*wrkVrtSet)){
498 if(!iVrt.Good || iVrt.selTrk.size()!=1)
continue;
500 for(
auto &jVrt : (*wrkVrtSet)){
501 if(!jVrt.Good || jVrt.selTrk.size()!=1)
continue;
503 if(iVrt.detachedTrack==jVrt.selTrk[0]){
505 newV.selTrk.push_back(iVrt.detachedTrack);
507 if(xAODwrk) vProb =
refitVertex(newV, xAODwrk->listJetTracks, *state,
true);
508 if(vProb>probVrtMergeLimit){
511 iVrt.detachedTrack = -1;
520 for(
auto &ntrVrt : (*wrkVrtSet)){
521 if(!ntrVrt.Good || ntrVrt.selTrk.size()<=1)
continue;
523 for(
auto tr : ntrVrt.selTrk){
524 for(
auto &onetVrt : (*wrkVrtSet)){
525 if(!onetVrt.Good || onetVrt.selTrk.size()!=1)
continue;
527 if(onetVrt.detachedTrack==tr){
529 newV.selTrk.push_back(onetVrt.selTrk[0]);
531 if(xAODwrk) vProb =
refitVertex( newV, xAODwrk->listJetTracks, *state,
true);
532 if(vProb>probVrtMergeLimit){
533 onetVrt.Good =
false;
544 for(
auto & curVrt : (*wrkVrtSet) ) {
545 if(!curVrt.Good )
continue;
546 if( std::abs(curVrt.vertex.z())>650. ){
550 if(curVrt.selTrk.size() != 1)
continue;
556 double Signif3D = 0.;
557 vrtVrtDist(primVrt,curVrt.vertex, curVrt.vertexCov, Signif3D);
558 double tmpProb = TMath::Prob(curVrt.chi2, 1);
560 bool trkGood =
false;
561 if(xAODwrk) trkGood =
check1TrVertexInPixel(xAODwrk->listJetTracks[curVrt.selTrk[0]], curVrt.vertex, curVrt.vertexCov);
563 if(trkGood && tmpProb>0.01){
565 std::vector<double> Impact,ImpactError;
566 double Signif3DP = 0;
567 if(xAODwrk) Signif3DP =
m_fitSvc->
VKalGetImpact(xAODwrk->listJetTracks[curVrt.selTrk[0]], primVrt.
position(), 1, Impact, ImpactError, *state);
571 h.m_hb_diffPS->Fill(Signif3DP, evtWgt);
582 for(
auto& curVrt : (*wrkVrtSet)){
583 if(!curVrt.Good )
continue;
584 int nth = curVrt.selTrk.size();
587 double Signif3D = 0.;
588 vrtVrtDist(primVrt,curVrt.vertex, curVrt.vertexCov, Signif3D);
589 if(xAODwrk) xAODwrk->tmpListTracks.resize(nth);
591 for(
int i = 0;
i < nth;
i++) {
592 if(xAODwrk) xAODwrk->tmpListTracks[
i] = xAODwrk->listJetTracks[curVrt.selTrk[
i]];
596 if(nth <= 1)
continue;
597 if(curVrt.projectedVrt<0.)
continue;
599 if(TMath::Prob( curVrt.chi2, 2*nth-3)<0.001)
continue;
601 if(nth==2 && xAODwrk){
604 if(!
check2TrVertexInPixel(xAODwrk->tmpListTracks[0],xAODwrk->tmpListTracks[1],curVrt.vertex,curVrt.vertexCov))
continue;
609 float ihitR = xAODwrk->tmpListTracks[0]->radiusOfFirstHit();
610 float jhitR = xAODwrk->tmpListTracks[1]->radiusOfFirstHit();
612 if(std::abs(ihitR-jhitR)>25.)
continue;
613 if(curVrt.vertex.perp()-
std::min(ihitR,jhitR) > 2.*vrErr)
continue;
621 h.m_hb_r2d->Fill(curVrt.vertex.perp(), evtWgt);
625 double xvt = curVrt.vertex.x();
626 double yvt = curVrt.vertex.y();
627 double zvt = curVrt.vertex.z();
628 double Rvt = std::hypot(xvt,yvt);
635 h.m_hb_r2d->Fill(curVrt.vertex.perp(), evtWgt);
639 if(nth==2 && curVrt.vertexCharge==0 && curVrt.detachedTrack<0){
640 double mass_PiPi = curVrt.vertexMom.M();
645 h.m_hb_massPiPi->Fill(mass_PiPi, evtWgt);
646 h.m_hb_massPPi ->Fill(mass_PPi, evtWgt);
647 if(curVrt.vertex.perp()>20.)
h.m_hb_massEE->Fill(mass_EE, evtWgt);
650 if(std::abs(mass_PiPi-
m_massK0) < 22.)
continue;
651 if(std::abs(mass_PPi-
m_massLam) < 8.)
continue;
652 if(mass_EE < 60. && curVrt.vertex.perp() > 20.)
continue;
657 h.m_hb_sig3DTot->Fill(Signif3D, evtWgt);
671 int ngoodVertices = 0;
672 std::vector<WrkVrt> goodVertices(0);
674 for(
auto & iv : (*wrkVrtSet)) {
675 int nth = iv.selTrk.size();
676 if(nth == 0)
continue;
680 goodVertices.emplace_back(iv);
681 if(nth==2) n2trVrt++;
682 if(nth>=3) nNtrVrt++;
686 if(ngoodVertices == 0 || (n2trVrt+nNtrVrt)==0){
687 return finalVertices;
692 bool swapDone =
false;
693 for(
unsigned int iv = 0; iv<goodVertices.size()-1; iv++) {
694 if(goodVertices[iv].
vertex.perp() > goodVertices[iv+1].vertex.perp()){
695 std::swap( goodVertices[iv], goodVertices[iv+1] );
703 bool swapDone =
false;
704 for(
unsigned int iv = 0; iv<goodVertices.size()-1; iv++) {
705 if(goodVertices[iv+1].
vertex.perp() > 400.)
continue;
706 if(goodVertices[iv].projectedVrt > goodVertices[iv+1].projectedVrt){
707 std::swap( goodVertices[iv], goodVertices[iv+1] );
715 if( goodVertices[1].vertexMom.M()-goodVertices[0].vertexMom.M() > 5000.){
716 std::swap( goodVertices[0], goodVertices[1] );
722 h.m_hb_distVV->Fill(
minVrtVrtDist(wrkVrtSet.get(), foundV1, foundV2), evtWgt);
731 std::vector<int> nonusedTrk(0);
732 for(
int trk=0; trk<nTracks; trk++){
733 bool present =
false;
735 for(
const auto& iv : (*wrkVrtSet)){
737 for(
auto t_in_V : iv.selTrk){
747 if(!present) nonusedTrk.push_back(trk);
750 struct MatchedSV{
int indVrt;
double Signif3D;};
751 std::vector<MatchedSV> matchSV(0);
753 for(
const auto& i_trk : nonusedTrk){
754 MatchedSV tmpV = {0, 1.e9};
756 for(
unsigned int iv = 0; iv<goodVertices.size(); iv++){
758 if(!goodVertices[iv].Good)
continue;
759 if(goodVertices[iv].selTrk.size()<2)
continue;
760 if(
vrtVrtDist(primVrt, goodVertices[iv].
vertex, goodVertices[iv].vertexCov, jetDir) < 10.)
continue;
763 std::vector<double> Impact, ImpactError;
764 if(xAODwrk) Signif =
m_fitSvc->
VKalGetImpact(xAODwrk->listJetTracks[i_trk], goodVertices[iv].vertex, 1, Impact, ImpactError, *state);
765 if(Signif<tmpV.Signif3D){
766 tmpV.Signif3D = Signif;
771 matchSV.push_back(tmpV);
774 for(
int iv = 0; iv<static_cast<int>(goodVertices.size()); iv++){
775 WrkVrt newV(goodVertices[iv]);
777 std::map<double,int> addTrk;
779 for(
unsigned int it = 0;
it<nonusedTrk.size();
it++){
780 if(matchSV[
it].indVrt==iv) addTrk[matchSV[
it].Signif3D]=
it;
786 newV.selTrk.push_back(nonusedTrk[atrk->second]);
794 newV.selTrk.push_back(nonusedTrk[atrk->second]);
800 if(xAODwrk) vProb =
refitVertex(newV, xAODwrk->listJetTracks, *state,
true);
801 if(vProb>0.01) goodVertices[iv] = newV;
803 std::vector<WrkVrt> TestVertices(1,newV);
804 if(xAODwrk) vProb =
improveVertexChi2(&TestVertices, 0, xAODwrk->listJetTracks, *state,
true);
805 if(vProb>0.01) goodVertices[iv] = TestVertices[0];
815 int ngoodVertices_final = 0;
817 TLorentzVector vertexMom;
818 bool isPrimaryVertex =
true;
820 for(
auto& vrt : goodVertices){
821 int nth = vrt.selTrk.size();
822 if(xAODwrk) xAODwrk->tmpListTracks.clear();
824 for(
int i = 0;
i < nth;
i++) {
825 int j = vrt.selTrk[
i];
826 if(xAODwrk) xAODwrk->tmpListTracks.push_back( xAODwrk->listJetTracks[j] );
831 if(nth==1)
h.m_hb_r1dc->Fill(vrt.vertex.perp(), evtWgt);
832 if(nth==2)
h.m_hb_r2dc->Fill(vrt.vertex.perp(), evtWgt);
833 if(nth==3)
h.m_hb_r3dc->Fill(vrt.vertex.perp(), evtWgt);
834 if(nth> 3)
h.m_hb_rNdc->Fill(vrt.vertex.perp(), evtWgt);
835 double Signif3D =
vrtVrtDist(primVrt, vrt.vertex, vrt.vertexCov, jetDir);
838 if(vrt.vertexCharge==0)
h.m_hb_totmass2T1->Fill(vrt.vertexMom.M(), evtWgt);
839 else h.m_hb_totmass2T2->Fill(vrt.vertexMom.M(), evtWgt);
840 h.m_hb_sig3D2tr->Fill(Signif3D , evtWgt);
844 else if(vrt.vertexMom.M()>6000.)
h.m_hb_sig3D1tr->Fill(Signif3D, evtWgt);
845 else h.m_hb_sig3DNtr->Fill(Signif3D, evtWgt);
851 return finalVertices;
856 std::vector<float> floatErrMtx;
857 floatErrMtx.resize(6);
858 for(
int i = 0;
i<6;
i++) floatErrMtx[
i] = vrt.vertexCov[
i];
863 std::vector<Trk::VxTrackAtVertex> & tmpVTAV = tmpVertex->
vxTrackAtVertex();
866 for(
int ii = 0; ii<nth; ii++) {
868 CovMtxP.setIdentity();
869 Trk::Perigee * tmpMeasPer =
new Trk::Perigee( 0., 0.,
874 std::move(CovMtxP) );
875 tmpVTAV.emplace_back( 1., tmpMeasPer );
886 finalVertices.push_back(tmpVertex);
887 ngoodVertices_final++;
888 if(nth==1) n1trVrt++;
891 isPrimaryVertex =
false;
895 vertexMom += vrt.vertexMom;
901 h.m_hb_goodvrtN->Fill(ngoodVertices_final+0.1, evtWgt);
902 h.m_curTup->ewgt = evtWgt;
903 if(n1trVrt)
h.m_hb_goodvrtN->Fill(n1trVrt+15., evtWgt);
907 if(ngoodVertices_final==0){
908 return finalVertices;
915 double totMass = vertexMom.M();
917 double eRatio = vertexMom.E()/momentumJet.E();
918 results.push_back(eRatio<1. ? eRatio : 1.);
919 results.push_back(vrt2TrackNumber);
921 if(xAODwrk)
results.push_back(xAODwrk->listSecondTracks.size());
923 results.push_back(momentumJet.E());
927 h.m_hb_ratio->Fill(
results[1], evtWgt);
928 h.m_hb_totmass->Fill(
results[0], evtWgt);
929 h.m_hb_nvrt2->Fill(
results[2], evtWgt);
930 h.m_hb_mom->Fill( momentumJet.Perp(), evtWgt);
933 return finalVertices;
1.8.18