17 #include "boost/graph/bron_kerbosch_all_cliques.hpp"
34 const double probVrtMergeLimit=0.01;
37 int inpNPart=xAODwrk->
inpTrk.size();
38 std::vector<const xAOD::NeutralParticle*> neutralPartDummy(0);
39 ATH_MSG_DEBUG(
"getVrtSecMulti() called with NPart=" << inpNPart);
41 std::vector<xAOD::Vertex*> finalVertices(0);
43 if( inpNPart < 2 ) {
return finalVertices;}
48 if( nTracks < 2 ) {
return finalVertices;}
54 h.m_hb_ntrksel->Fill( (
double) nTracks,
m_w_1);
65 std::map<long int,std::vector<double>> foundVrt2t;
69 ATH_MSG_DEBUG(
" Defined edges in the graph="<< num_edges(compatibilityGraph));
70 if(num_edges(compatibilityGraph)==0){
return finalVertices;}
80 std::unique_ptr<std::vector<WrkVrt>> wrkVrtSet = std::make_unique<std::vector<WrkVrt>>();
82 std::unique_ptr<Trk::IVKalState> state =
m_fitSvc->makeState();
84 long int NPTR=0, nth=2;
87 std::vector< std::vector<int> > allCliques;
88 bron_kerbosch_all_cliques(compatibilityGraph,
clique_visitor(allCliques));
89 for(
int cq=0; cq<(
int)allCliques.size();cq++){
90 newvrt.selTrk.clear();
91 NPTR=allCliques[cq].size();
92 for(
i=0;
i<NPTR;
i++){ newvrt.selTrk.push_back(allCliques[cq][
i]);}
100 std::vector<double> iniVrtPos=
estimVrtPos(nTracks,newvrt.selTrk,foundVrt2t);
101 m_fitSvc->setApproximateVertex(iniVrtPos[0], iniVrtPos[1], iniVrtPos[2], *state);
103 newvrt.vertex, newvrt.vertexMom, newvrt.vertexCharge, newvrt.vertexCov,
104 newvrt.chi2PerTrk, newvrt.trkAtVrt, newvrt.chi2,
106 if(
sc.isFailure() )
continue;
107 ATH_MSG_VERBOSE(
"Found IniVertex="<<newvrt.vertex[0]<<
", "<<newvrt.vertex[1]<<
", "<<newvrt.vertex[2]);
108 ATH_MSG_VERBOSE(
"with Chi2="<<newvrt.chi2<<
" Ntrk="<<NPTR<<
" trk1,2="<<newvrt.selTrk[0]<<
", "<<newvrt.selTrk[1]);
109 if(NPTR==2 && newvrt.chi2>10.)
continue;
110 if(newvrt.chi2PerTrk.size()==2) newvrt.chi2PerTrk[0]=newvrt.chi2PerTrk[1]=newvrt.chi2/2.;
112 newvrt.projectedVrt=
MomProjDist(newvrt.vertex, primVrt, newvrt.vertexMom);
113 wrkVrtSet->push_back(newvrt);
123 if((*wrkVrtSet).empty())
return finalVertices;
127 std::vector<int> trkNPairs(nTracks,0);
128 for(
auto &vrt : (*wrkVrtSet)){
129 int ntInV=vrt.selTrk.size()-1;
130 for(
auto &trk : vrt.selTrk)trkNPairs.at(trk) += ntInV;
136 std::multimap<double,std::pair<int,int>> vrtWithCommonTrk;
138 int nSoluI=(*wrkVrtSet).size();
139 vrtWithCommonTrk.clear();
140 unsigned int nTComMax=0;
141 for(
int iv=0; iv<nSoluI-1; iv++ ){
142 if(!(*wrkVrtSet)[iv].Good)
continue;
143 if( (*wrkVrtSet)[iv].selTrk.size()<nTComMax)
continue;
144 for(
int jv=iv+1; jv<nSoluI; jv++){
145 if(!(*wrkVrtSet)[jv].Good)
continue;
146 if( (*wrkVrtSet)[jv].selTrk.size()<nTComMax)
continue;
147 unsigned int nTCom=
nTrkCommon( wrkVrtSet.get(), iv, jv);
149 if(nTCom<nTComMax)
continue;
150 double sumChi2=(*wrkVrtSet)[iv].chi2+(*wrkVrtSet)[jv].chi2;
151 sumChi2=
std::min(sumChi2,999.)*1.e-3;
152 vrtWithCommonTrk.emplace(nTCom+sumChi2,std::make_pair(iv,jv));
155 if(vrtWithCommonTrk.empty())
break;
160 unsigned int nTCom=(*vrtWithCommonTrk.rbegin()).
first;
161 WrkVrt & v1 = (*wrkVrtSet)[(*vrtWithCommonTrk.rbegin()).
second.first];
162 WrkVrt &
v2 = (*wrkVrtSet)[(*vrtWithCommonTrk.rbegin()).
second.second];
164 if( nTCom==v1.
selTrk.size() || nTCom==
v2.selTrk.size() ){
165 if(nTCom==v1.
selTrk.size()){v1.
Good =
false;
continue;}
166 if(nTCom==
v2.selTrk.size()){
v2.Good =
false;
continue;}
169 if( nTCom>1 && TMath::Prob( v1.
chi2, 2*v1.
selTrk.size()-3) > probVrtMergeLimit
170 && TMath::Prob(
v2.chi2, 2*
v2.selTrk.size()-3) > probVrtMergeLimit){
172 if(prbV>probVrtMergeLimit){
173 v1.
Good =
false;
v2.Good =
false;
175 newvrt.projectedVrt=
MomProjDist(newvrt.vertex, primVrt, newvrt.vertexMom);
183 int cvgood=0;
for(
const auto& vrt:(*wrkVrtSet))
if(vrt.Good)cvgood++;
185 h.m_hb_rawVrtN->Fill( (
float)cvgood,
m_w_1);
191 for(
auto &v1t : (*wrkVrtSet)){
192 if(v1t.selTrk.size()!=1 || !v1t.Good)
continue;
193 int ind_t=v1t.selTrk[0];
194 if(trkNPairs[ind_t]<2){ v1t.Good=
false;
continue; }
196 for(
auto &vrt :(*wrkVrtSet)){
197 if(!vrt.Good || &v1t==&vrt)
continue;
198 if(
std::find(vrt.selTrk.begin(),vrt.selTrk.end(),ind_t) != vrt.selTrk.end()){ v1t.Good=
false;
break; }
204 int tmpV=0;
while( tmpV<(
int)(*wrkVrtSet).size() )
if( !(*wrkVrtSet)[tmpV].Good ) { (*wrkVrtSet).erase((*wrkVrtSet).begin()+tmpV);}
else {tmpV++;}
205 if((*wrkVrtSet).empty())
return finalVertices;
207 for(
auto &tmpV : (*wrkVrtSet) ) tmpV.projectedVrt=
MomProjDist(tmpV.vertex, primVrt, tmpV.vertexMom );
214 int foundV1=-1, foundV2=-1;
215 std::vector<double> checkedDst(0);
216 double minDistVV=
minVrtVrtDist( wrkVrtSet.get(), foundV1, foundV2, checkedDst);
219 h.m_hb_distVV->Fill( minDistVV,
m_w_1);
222 if(foundV1<foundV2) {
int tmp=foundV1; foundV1=foundV2; foundV2=
tmp;}
224 ATH_MSG_DEBUG(
"Merged vertex prob=" << probV<<
" Vrt1="<<foundV1<<
" Vrt2="<<foundV2<<
" dst="<<minDistVV);
225 if(probV<probVrtMergeLimit){
226 int pos=std::max_element(newvrt.chi2PerTrk.begin(),newvrt.chi2PerTrk.end())-newvrt.chi2PerTrk.begin();
227 newvrt.detachedTrack=newvrt.selTrk[
pos];
228 newvrt.selTrk.erase(newvrt.selTrk.begin()+
pos);
232 if(probV>probVrtMergeLimit && newvrt.vertexMom.M()<
m_vrtMassLimit){
233 newvrt.projectedVrt=
MomProjDist(newvrt.vertex, primVrt, newvrt.vertexMom);
235 (*wrkVrtSet)[foundV2].Good=
false;
236 (*wrkVrtSet)[foundV2].selTrk.clear();
237 }
else checkedDst.push_back(minDistVV);
238 minDistVV=
minVrtVrtDist( wrkVrtSet.get(), foundV1, foundV2, checkedDst);
242 for(
int iv=0; iv<(
int)wrkVrtSet->size(); iv++) {
243 if(!(*wrkVrtSet)[iv].Good )
continue;
244 if( (*wrkVrtSet)[iv].selTrk.size()<3 )
continue;
245 double tmpProb=TMath::Prob( (*wrkVrtSet)[iv].
chi2, 2*(*wrkVrtSet)[iv].selTrk.size()-3 );
249 (*wrkVrtSet)[iv].projectedVrt=
MomProjDist((*wrkVrtSet)[iv].
vertex, primVrt, (*wrkVrtSet)[iv].vertexMom);
254 for(
auto & iv : (*wrkVrtSet)){
256 ATH_MSG_DEBUG(
"Heavy vertex found Mass=" << iv.vertexMom.M());
259 iv.selTrk.erase( iv.selTrk.begin() + it_bad );
261 iv.projectedVrt=
MomProjDist(iv.vertex, primVrt, iv.vertexMom);
266 double signif3D=0., signif2D=0.;
269 for(
int iv=0; iv<(
int)wrkVrtSet->size(); iv++) {
270 WrkVrt & curVrt=(*wrkVrtSet)[iv];
271 nth=(*wrkVrtSet)[iv].
selTrk.size();
272 if(nth == 1)
continue;
273 if(!curVrt.
Good )
continue;
274 (*wrkVrtSet)[iv].Good =
false;
275 if(nth < 1)
continue;
276 if((*wrkVrtSet)[iv].projectedVrt<0.)
continue;
283 h.m_hb_sig3DTot->Fill( signif3D,
m_w_1);
284 if(nth==2)
h.m_hb_sig3D2tr->Fill( signif3D,
m_w_1);
285 if(nth >2)
h.m_hb_sig3DNtr->Fill( signif3D,
m_w_1);
296 h.m_hb_massPiPi->Fill( mass_PiPi,
m_w_1);
297 h.m_hb_massPPi ->Fill( mass_PPi,
m_w_1);
298 if( curVrt.
vertex.perp()>20.)
h.m_hb_massEE ->Fill( mass_EE,
m_w_1);
300 if( std::abs(mass_PiPi-
m_massK0) < 22.)
continue;
301 if( std::abs(mass_PPi-
m_massLam) < 8.)
continue;
302 if( mass_EE < 60. && curVrt.
vertex.perp() > 20.)
continue;
314 std::vector<double> impact,impactError;
315 for(
int iv=0; iv<(
int)wrkVrtSet->size(); iv++) {
316 WrkVrt & curVrt=(*wrkVrtSet)[iv];
318 if(!curVrt.
Good || nth<2)
continue;
319 double minPtT=1.e6, minSig3DT=1.e6, maxSig3DT=0.;
320 int ntrkBC=0,ntrkI=0,sumIBLHits=0,sumBLHits=0;
325 double SigR2 = impact[0]*impact[0]/impactError[0];
326 double SigZ2 = impact[1]*impact[1]/impactError[2];
327 minSig3DT=
std::min( minSig3DT, sqrt( SigR2 + SigZ2) );
328 maxSig3DT=
std::max( maxSig3DT, sqrt( SigR2 + SigZ2) );
336 float vProb=TMath::Prob(curVrt.
chi2, 2*nth-3);
338 float vrtR=curVrt.
vertex.perp();
339 TLorentzVector SVPV(curVrt.
vertex.x()-primVrt.
x(),curVrt.
vertex.y()-primVrt.
y(),curVrt.
vertex.z()-primVrt.
z(), 10.);
345 h.m_curTup->NVrtTrk [
h.m_curTup->nNVrt] = nth;
346 h.m_curTup->NVrtTrkHF [
h.m_curTup->nNVrt] = ntrkBC;
347 h.m_curTup->NVrtTrkI [
h.m_curTup->nNVrt] = ntrkI;
348 h.m_curTup->NVrtProb [
h.m_curTup->nNVrt] = vProb;
349 h.m_curTup->NVrtSig3D [
h.m_curTup->nNVrt] = signif3D;
350 h.m_curTup->NVrtSig2D [
h.m_curTup->nNVrt] = signif2D;
351 h.m_curTup->NVrtDist2D[
h.m_curTup->nNVrt] = vrtR<20. ? Dist2D : vrtR;
352 h.m_curTup->NVrtM [
h.m_curTup->nNVrt] = curVrt.
vertexMom.M();
353 h.m_curTup->NVrtPt [
h.m_curTup->nNVrt] = curVrt.
vertexMom.Pt();
354 h.m_curTup->NVrtEta [
h.m_curTup->nNVrt] = SVPV.Eta();
355 h.m_curTup->NVrtIBL [
h.m_curTup->nNVrt] = sumIBLHits;
356 h.m_curTup->NVrtBL [
h.m_curTup->nNVrt] = sumBLHits;
357 h.m_curTup->NVrtCosSPM[
h.m_curTup->nNVrt] = cosSVPVM;
359 h.m_curTup->NVMinPtT [
h.m_curTup->nNVrt] = minPtT;
360 h.m_curTup->NVMinS3DT [
h.m_curTup->nNVrt] = minSig3DT;
361 h.m_curTup->NVrtBDT [
h.m_curTup->nNVrt] = 1.1;
362 h.m_curTup->NVrtHR1 [
h.m_curTup->nNVrt] = xAODwrk->
listSelTracks[curVrt.
selTrk[0]]->radiusOfFirstHit();
363 h.m_curTup->NVrtHR2 [
h.m_curTup->nNVrt] = xAODwrk->
listSelTracks[curVrt.
selTrk[1]]->radiusOfFirstHit();
372 std::vector<float> VARS(10);
374 VARS[1]=
log(curVrtPt);
376 VARS[3]=
log(vrtR<20. ? SVPV.Perp() : vrtR);
378 VARS[5]=
log(maxSig3DT);
380 VARS[7]=sqrt(std::abs(1.-cosSVPVM*cosSVPVM));
383 float wgtSelect=
m_SV2T_BDT->GetGradBoostMVA(VARS);
384 curVrt.
BDT=wgtSelect;
387 h.m_hb_fakeSVBDT->Fill(wgtSelect,1.);
388 h.m_curTup->NVrtBDT[
h.m_curTup->nNVrt-1] = wgtSelect;
394 for(
auto &vtmp : (*wrkVrtSet)){
395 if(vtmp.selTrk.size()!=1 || (!vtmp.Good))
continue;
396 if(
it==vtmp.detachedTrack)vtmp.Good=
false;
406 for(
auto & vrt : (*wrkVrtSet)) {
407 if( !vrt.Good || vrt.selTrk.size() != 1 )
continue;
409 m_fitSvc->VKalGetImpact(xaodtp, primVrt.
position(), 1, impact, impactError);
410 double SigR2 = std::abs(impact[0]*impact[0]/impactError[0]);
411 double SigZ2 = std::abs(impact[1]*impact[1]/impactError[2]);
412 float dist2D=
vrtVrtDist2D(primVrt,vrt.vertex, vrt.vertexCov, signif2D);
413 h.m_curTup->NVrtTrk [
h.m_curTup->nNVrt] = 1;
414 h.m_curTup->NVrtTrkHF [
h.m_curTup->nNVrt] =
getIdHF(xaodtp);
415 h.m_curTup->NVrtProb [
h.m_curTup->nNVrt] = trkNPairs[vrt.selTrk[0]];
416 h.m_curTup->NVrtSig3D [
h.m_curTup->nNVrt] = 0.;
417 h.m_curTup->NVrtSig2D [
h.m_curTup->nNVrt] = signif2D;
418 h.m_curTup->NVrtDist2D[
h.m_curTup->nNVrt] = dist2D;
419 h.m_curTup->NVrtM [
h.m_curTup->nNVrt] = vrt.vertexMom.M();
420 h.m_curTup->NVrtPt [
h.m_curTup->nNVrt] = vrt.vertexMom.Pt();
421 h.m_curTup->NVrtCosSPM[
h.m_curTup->nNVrt] = 0.;
422 h.m_curTup->NVrtCh [
h.m_curTup->nNVrt] = vrt.vertexCharge;
423 h.m_curTup->NVMinPtT [
h.m_curTup->nNVrt] = xaodtp->pt();
424 h.m_curTup->NVMinS3DT [
h.m_curTup->nNVrt] = sqrt(SigR2 + SigZ2);
425 h.m_curTup->NVrtBDT [
h.m_curTup->nNVrt] = 1.1;
432 std::multimap<double,WrkVrt,std::greater<double>> goodVertexMap;
434 for(
auto & iv : (*wrkVrtSet) ) {
435 nth=iv.selTrk.size();
436 if(nth==1)iv.BDT=-2.;
439 else if(nth>2)
selector=iv.BDT+iv.vertexMom.M()*1.e-5;
440 if( iv.Good && nth>0 ) {
450 return finalVertices;
460 for(
auto & iv : goodVertexMap){
461 WrkVrt & curVrt=iv.second;
476 if(!tmpVertex)
continue;
479 std::vector<float> floatErrMtx(6);
483 std::vector<Trk::VxTrackAtVertex> & tmpVTAV=tmpVertex->
vxTrackAtVertex(); tmpVTAV.clear();
485 CovMtxP.setIdentity();
491 std::move(CovMtxP) );
492 tmpVTAV.emplace_back( 1., tmpMeasPer );
500 wgtBDT (*tmpVertex) =curVrt.
BDT;
501 nTrksDec(*tmpVertex) =curVrt.
selTrk.size();
503 finalVertices.push_back(tmpVertex);
508 h.m_hb_goodvrtN->Fill( finalVertices.size()+0.1,
m_w_1);
509 h.m_hb_goodvrt1N->Fill( n1trVrt+0.1,
m_w_1);
514 return finalVertices;