46{
47
48 const double probVrtMergeLimit = 0.01;
49
50 int inpNPart = 0;
51 if (xAODwrk) {
52 inpNPart = xAODwrk->InpTrk.size();
53 xAODwrk->FoundSecondTracks.clear();
55 }
56 ATH_MSG_DEBUG(
"InDet getVrtSecMulti() called with NPart=" << inpNPart );
57
58 std::vector<xAOD::Vertex*> finalVertices(0);
59
60 if (inpNPart < 2) {
61 return finalVertices;
62 }
63
64 float evtWgt = 1.;
65 const EventContext & ctx = Gaudi::Hive::currentContext();
68 if(eventInfo->hasBeamSpotWeight()) evtWgt *= eventInfo->beamSpotWeight();
70
71
72 long int nTracks = 0;
73 TLorentzVector momentumJet;
74 int nRefPVTrk = 0;
75 if (xAODwrk) {
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());
80 }
81 nTracks = xAODwrk->listJetTracks.size();
82 momentumJet =
totalMom(xAODwrk->listJetTracks);
83 }
84
85 if (nTracks < 2) {
86 return finalVertices;
87 }
88
89
90 ATH_MSG_DEBUG(
"Number of selected tracks inside jet= " << nTracks);
91
94 h.m_hb_jmom->Fill(momentumJet.Perp(), evtWgt);
95 h.m_hb_ntrkjet->Fill(
static_cast<double>(nTracks), evtWgt);
96 }
97
98
99
100
101
102
103
104
105
106 std::vector<double> inpMass(nTracks,
m_massPi);
107 double vrt2TrackNumber = 0;
108
109 if (xAODwrk) {
111 xAODwrk->listJetTracks,
112 xAODwrk->TracksForFit,
113 primVrt,
114 jetDir,
115 inpMass,
116 nRefPVTrk,
117 xAODwrk->TrkFromV0,
118 xAODwrk->listSecondTracks,
119 compatibilityGraph);
120 if (xAODwrk->TrkFromV0.size() > 1) {
123 }
124 vrt2TrackNumber = xAODwrk->listSecondTracks.size()/2.0;
127 }
128
129
130 ATH_MSG_DEBUG(
" Found different tracks in pairs=" << vrt2TrackNumber);
131
132
133
134
135
136
137
138
139
140
141
142 std::unique_ptr<std::vector<WrkVrt>> wrkVrtSet = std::make_unique<std::vector<WrkVrt>>();
145 std::unique_ptr<Trk::IVKalState> state =
m_fitSvc->makeState(ctx);
147
148
149
150 std::vector<std::vector<int>> allCliques;
151 bron_kerbosch_all_cliques(compatibilityGraph, clique_visitor(allCliques));
152
153 for (const auto& clique : allCliques) {
154
155 newvrt.selTrk.clear();
156 if (xAODwrk) xAODwrk->tmpListTracks.clear();
157
158 for (int i_trk : clique){
159 newvrt.selTrk.push_back(i_trk);
160 if(xAODwrk) xAODwrk->tmpListTracks.push_back(xAODwrk->listJetTracks.at(i_trk));
161 }
162
164 if (
sc.isFailure() ||
166
167 m_fitSvc->setApproximateVertex(primVrt.
x(),
170 *state);
172 m_fitSvc->setApproximateVertex(0., 0., 0., *state);
173 } else {
175 double jetVrtDir = jetDir.Px() * vDist.x() + jetDir.Py() * vDist.y() +
176 jetDir.Pz() * vDist.z();
178 if (jetVrtDir > 0.) {
179
180 m_fitSvc->setApproximateVertex(newvrt.vertex.x(),
181 newvrt.vertex.y(),
182 newvrt.vertex.z(),
183 *state);
184 } else {
185 m_fitSvc->setApproximateVertex(primVrt.
x(), primVrt.
y(), primVrt.
z(), *state);
186 }
187 }
188
189 sc = StatusCode::FAILURE;
190 if (xAODwrk) {
192 newvrt.vertex,
193 newvrt.vertexMom,
194 newvrt.vertexCharge,
195 newvrt.vertexCov,
196 newvrt.chi2PerTrk,
197 newvrt.trkAtVrt,
198 newvrt.chi2,
199 *state,
200 false);
201 }
202 if (
sc.isFailure())
continue;
203
204 if (clique.size() == 2 && newvrt.chi2 > 10.) continue;
205
206 if (newvrt.chi2PerTrk.size() == 2){
207 newvrt.chi2PerTrk[0] = newvrt.chi2PerTrk[1] = newvrt.chi2 / 2.;
208 }
209 newvrt.Good = true;
210 newvrt.nCloseVrt = 0;
211 newvrt.dCloseVrt = 1000000.;
212 newvrt.projectedVrt =
214 wrkVrtSet->push_back(newvrt);
215 }
216
217
218
219
220 bool disassembled = false;
221
222 do {
223 disassembled = false;
224 int NSoluI = (*wrkVrtSet).size();
225 for (int iv = 0; iv < NSoluI; iv++) {
226 WrkVrt vrt = (*wrkVrtSet)[iv];
227 if (!vrt.Good || vrt.selTrk.size() == 2) continue;
228 if (TMath::Prob(vrt.chi2, 2 * vrt.selTrk.size() - 3) < 1.e-3) {
229 if (xAODwrk)
disassembleVertex(wrkVrtSet.get(), iv, xAODwrk->listJetTracks, *state);
230 disassembled = true;
231 }
232 }
233 }while(disassembled);
234
235
236 while( (*wrkVrtSet).size()>1 ){
237 int tmpN = (*wrkVrtSet).size();
238
239 int iv = 0;
240 for(; iv<tmpN-1; iv++){
241 WrkVrt vert_i = (*wrkVrtSet)[iv];
242 if(!vert_i.Good ) continue;
243 int ntrk_i = (*wrkVrtSet)[iv].
selTrk.size();
244
245 int jv = iv+1;
246 for(; jv<tmpN; jv++){
247 WrkVrt vert_j = (*wrkVrtSet)[jv];
248 if(!vert_j.Good ) continue;
249 int ntrk_j = (*wrkVrtSet)[jv].
selTrk.size();
250
251 int nTCom =
nTrkCommon(wrkVrtSet.get(), iv, jv);
252 if(nTCom==ntrk_i){
253 (*wrkVrtSet).erase((*wrkVrtSet).begin()+iv);
254 break;
255 }
256 else if(nTCom==ntrk_j){
257 (*wrkVrtSet).erase((*wrkVrtSet).begin()+jv);
258 break;
259 }
260 }
261 if(jv!=tmpN) break;
262
263 }
264 if(iv==tmpN-1) break;
265
266 }
267
268
269
270
271 std::multimap<int,std::pair<int,int>> vrtWithCommonTrk;
272 std::multimap<int,std::pair<int,int>>::reverse_iterator icvrt;
273 do{
274 int NSoluI = (*wrkVrtSet).size();
275 vrtWithCommonTrk.clear();
276 for(int iv = 0; iv<NSoluI-1; iv++ ){
277 for(int jv = iv+1; jv<NSoluI; jv++){
278 if(!(*wrkVrtSet)[iv].Good || !(*wrkVrtSet)[jv].Good) continue;
279 int nTCom=
nTrkCommon(wrkVrtSet.get(), iv, jv);
280 if(!nTCom)continue;
281 vrtWithCommonTrk.emplace(nTCom,std::make_pair(iv,jv));
282 }
283 }
284
285 for(icvrt = vrtWithCommonTrk.rbegin(); icvrt!=vrtWithCommonTrk.rend(); ++icvrt){
286 int nTCom = (*icvrt).first;
287 int iv = (*icvrt).second.first;
288 int jv = (*icvrt).second.second;
289 int nTrkI = (*wrkVrtSet)[iv].selTrk.size();
290 int nTrkJ = (*wrkVrtSet)[jv].selTrk.size();
291 double probV = -1.;
292 if (xAODwrk) {
294 }
295
296 if(probV<probVrtMergeLimit){
297 if(nTrkI==2 || nTrkJ==2 || nTCom<2) {
298 continue;
299 }
300 if(nTCom>nTrkI-nTCom || nTCom>nTrkJ-nTCom){
301
303 break;
304 }
305 continue;
306 }
307
308 newvrt.Good = true;
309 (*wrkVrtSet).push_back(newvrt);
310 (*wrkVrtSet)[iv].Good = false;
311 (*wrkVrtSet)[jv].Good = false;
312 break;
313 }
314
315 } while( icvrt != vrtWithCommonTrk.rend() );
316
317
319 int cvgood=0;
320 for(const auto& vrt : (*wrkVrtSet)){
321 if(vrt.Good) cvgood++;
322 }
324 h.m_hb_rawVrtN->Fill(
static_cast<float>(cvgood), evtWgt);
325 }
326
327
328 for(auto &tmpV : (*wrkVrtSet) ) {
329 if(tmpV.vertex.perp()>
m_rLayer3+10.) tmpV.Good =
false;
330 TLorentzVector SVPV(tmpV.vertex.x()-primVrt.
x(),
331 tmpV.vertex.y()-primVrt.
y(),
332 tmpV.vertex.z()-primVrt.
z(), 1.);
334 }
335
336 unsigned int tmpV = 0;
337 while( tmpV<(*wrkVrtSet).size() ){
338 if( !(*wrkVrtSet)[tmpV].Good ) (*wrkVrtSet).erase((*wrkVrtSet).begin()+tmpV);
339 else tmpV++;
340 }
341
342 if((*wrkVrtSet).empty()){
343 return finalVertices;
344 }
345
346 std::vector< std::vector<float> > trkScore(0);
348 for(
auto &trk : xAODwrk->listJetTracks) trkScore.push_back(
m_trackClassificator->trkTypeWgts(ctx, trk, primVrt, jetDir));
349 }
350
351 for(
auto &tmpV : (*wrkVrtSet) ) tmpV.projectedVrt=
jetProjDist(tmpV.vertex, primVrt, jetDir);
352
353
354
355
356
357
358 std::unique_ptr<std::vector< std::deque<long int> > > trkInVrt = std::make_unique<std::vector<std::deque<long int>>>(nTracks);
360
361 double foundMaxT;
362 long int selectedTrack, selectedVertex;
363 int foundV1, foundV2;
364
366 while( ( foundMaxT =
MaxOfShared(wrkVrtSet.get(), trkInVrt.get(), selectedTrack, selectedVertex) ) > 0) {
367
368 double foundMinVrtDst = 1000000.;
370
371
373 &&
nTrkCommon(wrkVrtSet.get(), foundV1, foundV2) ){
374
375 bool vrtMerged=false;
377 if(foundV1<foundV2)
std::swap(foundV1, foundV2);
378
379 double probV = 0.;
380 if (xAODwrk) probV =
mergeAndRefitVertices(wrkVrtSet.get(), foundV1, foundV2, newvrt, xAODwrk->listJetTracks, *state);
381
382 if(probV>probVrtMergeLimit && newvrt.vertexMom.M()<c_vrtBCMassLimit){
383
384 double tstDst =
jetProjDist(newvrt.vertex, primVrt, jetDir);
385 if(tstDst>0.){
386
387 std::swap((*wrkVrtSet)[foundV1], newvrt);
388 (*wrkVrtSet)[foundV1].projectedVrt = tstDst;
389 (*wrkVrtSet)[foundV2].Good = false;
390 (*wrkVrtSet)[foundV2].selTrk.clear();
391 vrtMerged=true;
392 }
393 }
394
395 (*wrkVrtSet)[foundV1].nCloseVrt = -1;
396 (*wrkVrtSet)[foundV1].dCloseVrt = 1000000.;
397 (*wrkVrtSet)[foundV2].nCloseVrt = -1;
398 (*wrkVrtSet)[foundV2].dCloseVrt = 1000000.;
400
401 }
402
403 if(vrtMerged){
404 trkInVrt->resize(nTracks);
406 continue;
407 }
408
409 }
410
412
413 sc = StatusCode::FAILURE;
414 if(xAODwrk)
sc =
refitVertex( wrkVrtSet.get(), selectedVertex, xAODwrk->listJetTracks, *state,
false);
415
416 (*wrkVrtSet)[selectedVertex].projectedVrt =
jetProjDist((*wrkVrtSet)[selectedVertex].vertex, primVrt, jetDir);
417
418 if(
sc.isFailure() ) (*wrkVrtSet)[selectedVertex].Good =
false;
419 if( (*wrkVrtSet)[selectedVertex].projectedVrt<0.
420 && (*wrkVrtSet)[selectedVertex].selTrk.size()==2 ){
421 (*wrkVrtSet)[selectedVertex].Good = false;
422 }
423
424 }
425
426
427
428
429
430 double minDistVV =
minVrtVrtDist( wrkVrtSet.get(), foundV1, foundV2);
432 if(foundV1<foundV2)
std::swap(foundV1, foundV2);
433
434 double probV = 0.;
435 if(xAODwrk) probV =
mergeAndRefitVertices(wrkVrtSet.get(), foundV1, foundV2, newvrt, xAODwrk->listJetTracks, *state);
436
437 if(probV>probVrtMergeLimit && newvrt.vertexMom.M()<c_vrtBCMassLimit){
438
439 double tstDst =
jetProjDist(newvrt.vertex, primVrt, jetDir);
440 if(tstDst>0.){
441
443 (*wrkVrtSet)[foundV1].projectedVrt = tstDst;
444 (*wrkVrtSet)[foundV2].Good = false;
445 (*wrkVrtSet)[foundV2].selTrk.clear();
446 }
447 }
448
449 (*wrkVrtSet)[foundV1].nCloseVrt = -1;
450 (*wrkVrtSet)[foundV1].dCloseVrt = 1000000.;
451 (*wrkVrtSet)[foundV2].nCloseVrt = -1;
452 (*wrkVrtSet)[foundV2].dCloseVrt = 1000000.;
454 }
455
456
457
458 for(unsigned int iv=0; iv<wrkVrtSet->size(); iv++) {
459 WrkVrt vert_i = (*wrkVrtSet)[iv];
460 if(!vert_i.Good) continue;
461 if( vert_i.selTrk.size()<3 ) continue;
462
463 double tmpProb = TMath::Prob( vert_i.chi2, 2*vert_i.selTrk.size()-3 );
464 if(tmpProb<0.001){
465 if(xAODwrk) tmpProb =
improveVertexChi2(wrkVrtSet.get(), iv, xAODwrk->listJetTracks, *state,
false);
466 if(tmpProb<0.001) vert_i.Good = false;
467 vert_i.projectedVrt =
jetProjDist(vert_i.vertex, primVrt, jetDir);
468 }
469 }
470
471
473
474
475
476 for(auto &ntrVrt : (*wrkVrtSet)){
477 if(!ntrVrt.Good || ntrVrt.selTrk.size()<=1) continue;
478
479 for(auto &onetVrt : (*wrkVrtSet)){
480 if(!onetVrt.Good || onetVrt.selTrk.size()!=1) continue;
481
482 if(ntrVrt.detachedTrack==onetVrt.selTrk[0]){
484 newV.selTrk.push_back(ntrVrt.detachedTrack);
485 double vProb = 0.;
486 if(xAODwrk) vProb =
refitVertex(newV, xAODwrk->listJetTracks, *state,
true);
487 if(vProb>probVrtMergeLimit){
488 onetVrt.Good=false;
489 ntrVrt=newV;
490 ntrVrt.detachedTrack=-1;
491 }
492 break;
493 }
494 }
495 }
496
497
498 for(auto &iVrt : (*wrkVrtSet)){
499 if(!iVrt.Good || iVrt.selTrk.size()!=1) continue;
500
501 for(auto &jVrt : (*wrkVrtSet)){
502 if(!jVrt.Good || jVrt.selTrk.size()!=1) continue;
503
504 if(iVrt.detachedTrack==jVrt.selTrk[0]){
506 newV.selTrk.push_back(iVrt.detachedTrack);
507 double vProb = 0.;
508 if(xAODwrk) vProb =
refitVertex(newV, xAODwrk->listJetTracks, *state,
true);
509 if(vProb>probVrtMergeLimit){
510 jVrt.Good = false;
511 iVrt = newV;
512 iVrt.detachedTrack = -1;
513 }
514 break;
515 }
516 }
517
518 }
519
520
521 for(auto &ntrVrt : (*wrkVrtSet)){
522 if(!ntrVrt.Good || ntrVrt.selTrk.size()<=1) continue;
523
524 for(auto tr : ntrVrt.selTrk){
525 for(auto &onetVrt : (*wrkVrtSet)){
526 if(!onetVrt.Good || onetVrt.selTrk.size()!=1) continue;
527
528 if(onetVrt.detachedTrack==tr){
530 newV.selTrk.push_back(onetVrt.selTrk[0]);
531 double vProb = 0.;
532 if(xAODwrk) vProb =
refitVertex( newV, xAODwrk->listJetTracks, *state,
true);
533 if(vProb>probVrtMergeLimit){
534 onetVrt.Good = false;
535 ntrVrt = newV;
536 }
537 }
538
539 }
540 }
541
542 }
543
544
545 for(auto & curVrt : (*wrkVrtSet) ) {
546 if(!curVrt.Good ) continue;
547 if( std::abs(curVrt.vertex.z())>650. ){
548 curVrt.Good = false;
549 continue;
550 }
551 if(curVrt.selTrk.size() != 1) continue;
552
553 curVrt.Good = false;
554
556
557 double Signif3D = 0.;
558 vrtVrtDist(primVrt,curVrt.vertex, curVrt.vertexCov, Signif3D);
559 double tmpProb = TMath::Prob(curVrt.chi2, 1);
560
561 bool trkGood = false;
562 if(xAODwrk) trkGood =
check1TrVertexInPixel(xAODwrk->listJetTracks[curVrt.selTrk[0]], curVrt.vertex, curVrt.vertexCov);
563
564 if(trkGood && tmpProb>0.01){
565
566 std::vector<double> Impact,ImpactError;
567 double Signif3DP = 0;
568 if(xAODwrk) Signif3DP =
m_fitSvc->VKalGetImpact(xAODwrk->listJetTracks[curVrt.selTrk[0]], primVrt.
position(), 1, Impact, ImpactError, *state);
569
572 h.m_hb_diffPS->Fill(Signif3DP, evtWgt);
573 }
574
576 }
577
578 }
579
580 }
581
582
583 for(auto& curVrt : (*wrkVrtSet)){
584 if(!curVrt.Good ) continue;
585 int nth = curVrt.selTrk.size();
586 if(nth==1) continue;
587
588 double Signif3D = 0.;
589 vrtVrtDist(primVrt,curVrt.vertex, curVrt.vertexCov, Signif3D);
590 if(xAODwrk) xAODwrk->tmpListTracks.resize(nth);
591
592 for(
int i = 0;
i < nth;
i++) {
593 if(xAODwrk) xAODwrk->tmpListTracks[
i] = xAODwrk->listJetTracks[curVrt.selTrk[
i]];
594 }
595 curVrt.Good = false;
596
597 if(nth <= 1) continue;
598 if(curVrt.projectedVrt<0.) continue;
599
600 if(TMath::Prob( curVrt.chi2, 2*nth-3)<0.001) continue;
601
602 if(nth==2 && xAODwrk){
603
605 if(!
check2TrVertexInPixel(xAODwrk->tmpListTracks[0],xAODwrk->tmpListTracks[1],curVrt.vertex,curVrt.vertexCov))
continue;
606 }
607
608
610 float ihitR = xAODwrk->tmpListTracks[0]->radiusOfFirstHit();
611 float jhitR = xAODwrk->tmpListTracks[1]->radiusOfFirstHit();
613 if(std::abs(ihitR-jhitR)>25.) continue;
614 if(curVrt.vertex.perp()-std::min(ihitR,jhitR) > 2.*vrErr) continue;
615 }
616
617 }
618
619
622 h.m_hb_r2d->Fill(curVrt.vertex.perp(), evtWgt);
623 }
624
626 double xvt = curVrt.vertex.x();
627 double yvt = curVrt.vertex.y();
628 double zvt = curVrt.vertex.z();
629 double Rvt = std::hypot(xvt,yvt);
632 }
633
636 h.m_hb_r2d->Fill(curVrt.vertex.perp(), evtWgt);
637 }
638
639
640 if(nth==2 && curVrt.vertexCharge==0 && curVrt.detachedTrack<0){
641 double mass_PiPi = curVrt.vertexMom.M();
646 h.m_hb_massPiPi->Fill(mass_PiPi, evtWgt);
647 h.m_hb_massPPi ->Fill(mass_PPi, evtWgt);
648 if(curVrt.vertex.perp()>20.)
h.m_hb_massEE->Fill(mass_EE, evtWgt);
649 }
650
651 if(std::abs(mass_PiPi-
m_massK0) < 22.)
continue;
652 if(std::abs(mass_PPi-
m_massLam) < 8.)
continue;
653 if(mass_EE < 60. && curVrt.vertex.perp() > 20.) continue;
654 }
655
658 h.m_hb_sig3DTot->Fill(Signif3D, evtWgt);
659 }
661
662 curVrt.Good = true;
663
664 }
665
666
668
669
670 int n2trVrt = 0;
671 int nNtrVrt = 0;
672 int ngoodVertices = 0;
673 std::vector<WrkVrt> goodVertices(0);
674
675 for(auto & iv : (*wrkVrtSet)) {
676 int nth = iv.selTrk.size();
677 if(nth == 0) continue;
678
679 if(iv.Good){
680 ngoodVertices++;
681 goodVertices.emplace_back(iv);
682 if(nth==2) n2trVrt++;
683 if(nth>=3) nNtrVrt++;
684 }
685 }
686
687 if(ngoodVertices == 0 || (n2trVrt+nNtrVrt)==0){
688 return finalVertices;
689 }
690
691
692 while(1){
693 bool swapDone = false;
694 for(unsigned int iv = 0; iv<goodVertices.size()-1; iv++) {
695 if(goodVertices[iv].
vertex.perp() > goodVertices[iv+1].vertex.perp()){
696 std::swap( goodVertices[iv], goodVertices[iv+1] );
697 swapDone = true;
698 }
699 }
700 if(!swapDone) break;
701 }
702
703 while(1){
704 bool swapDone = false;
705 for(unsigned int iv = 0; iv<goodVertices.size()-1; iv++) {
706 if(goodVertices[iv+1].
vertex.perp() > 400.)
continue;
707 if(goodVertices[iv].projectedVrt > goodVertices[iv+1].projectedVrt){
708 std::swap( goodVertices[iv], goodVertices[iv+1] );
709 swapDone = true;
710 }
711 }
712 if(!swapDone) break;
713 }
714
715 if(ngoodVertices>1){
716 if( goodVertices[1].vertexMom.M()-goodVertices[0].vertexMom.M() > 5000.){
717 std::swap( goodVertices[0], goodVertices[1] );
718 }
719 }
720
723 h.m_hb_distVV->Fill(
minVrtVrtDist(wrkVrtSet.get(), foundV1, foundV2), evtWgt);
724 }
725
726
727
728
729
730
731
732 std::vector<int> nonusedTrk(0);
733 for(int trk=0; trk<nTracks; trk++){
734 bool present = false;
735
736 for(const auto& iv : (*wrkVrtSet)){
737 if(iv.Good){
738 for(auto t_in_V : iv.selTrk){
739 if(trk==t_in_V){
740 present = true;
741 break;
742 }
743 }
744 }
745 if(present) break;
746 }
747
748 if(!present) nonusedTrk.push_back(trk);
749 }
750
751 struct MatchedSV{int indVrt; double Signif3D;};
752 std::vector<MatchedSV> matchSV(0);
753
754 for(const auto& i_trk : nonusedTrk){
755 MatchedSV tmpV = {0, 1.e9};
756
757 for(unsigned int iv = 0; iv<goodVertices.size(); iv++){
758
759 if(!goodVertices[iv].Good) continue;
760 if(goodVertices[iv].selTrk.size()<2) continue;
761 if(
vrtVrtDist(primVrt, goodVertices[iv].vertex, goodVertices[iv].vertexCov, jetDir) < 10.)
continue;
762
763 double Signif = 0.;
764 std::vector<double> Impact, ImpactError;
765 if(xAODwrk) Signif =
m_fitSvc->VKalGetImpact(xAODwrk->listJetTracks[i_trk], goodVertices[iv].vertex, 1, Impact, ImpactError, *state);
766 if(Signif<tmpV.Signif3D){
767 tmpV.Signif3D = Signif;
768 tmpV.indVrt=iv;
769 }
770 }
771
772 matchSV.push_back(tmpV);
773 }
774
775 for(int iv = 0; iv<static_cast<int>(goodVertices.size()); iv++){
776 WrkVrt newV(goodVertices[iv]);
777 bool addedT = false;
778 std::map<double,int> addTrk;
779
780 for(
unsigned int it = 0;
it<nonusedTrk.size();
it++){
781 if(matchSV[it].indVrt==iv) addTrk[matchSV[
it].Signif3D]=
it;
782 }
783
784 std::map<double,int>::iterator atrk=addTrk.begin();
785 if(!addTrk.empty()){
786 if(atrk->first<4.){
787 newV.selTrk.push_back(nonusedTrk[atrk->second]);
788 addedT=true;
789 }
790 }
791
792 if(addTrk.size()>1){
793 ++atrk;
794 if(atrk->first<4.){
795 newV.selTrk.push_back(nonusedTrk[atrk->second]);
796 }
797 }
798
799 if(addedT){
800 double vProb = 0.;
801 if(xAODwrk) vProb =
refitVertex(newV, xAODwrk->listJetTracks, *state,
true);
802 if(vProb>0.01) goodVertices[iv] = newV;
803 else{
804 std::vector<WrkVrt> TestVertices(1,newV);
805 if(xAODwrk) vProb =
improveVertexChi2(&TestVertices, 0, xAODwrk->listJetTracks, *state,
true);
806 if(vProb>0.01) goodVertices[iv] = TestVertices[0];
807 }
808 }
809
810 }
811
812
813
814
815
816 int ngoodVertices_final = 0;
817 int n1trVrt = 0;
818 TLorentzVector vertexMom;
819 bool isPrimaryVertex = true;
820
821 for(auto& vrt : goodVertices){
822 int nth = vrt.selTrk.size();
823 if(xAODwrk) xAODwrk->tmpListTracks.clear();
824
825 for(
int i = 0;
i < nth;
i++) {
826 int j = vrt.selTrk[
i];
827 if(xAODwrk) xAODwrk->tmpListTracks.push_back( xAODwrk->listJetTracks[j] );
828 }
829
832 if(nth==1)
h.m_hb_r1dc->Fill(vrt.vertex.perp(), evtWgt);
833 if(nth==2)
h.m_hb_r2dc->Fill(vrt.vertex.perp(), evtWgt);
834 if(nth==3)
h.m_hb_r3dc->Fill(vrt.vertex.perp(), evtWgt);
835 if(nth> 3)
h.m_hb_rNdc->Fill(vrt.vertex.perp(), evtWgt);
836 double Signif3D =
vrtVrtDist(primVrt, vrt.vertex, vrt.vertexCov, jetDir);
837
838 if(nth==2){
839 if(vrt.vertexCharge==0)
h.m_hb_totmass2T1->Fill(vrt.vertexMom.M(), evtWgt);
840 else h.m_hb_totmass2T2->Fill(vrt.vertexMom.M(), evtWgt);
841 h.m_hb_sig3D2tr->Fill(Signif3D , evtWgt);
843
844 }
845 else if(vrt.vertexMom.M()>6000.)
h.m_hb_sig3D1tr->Fill(Signif3D, evtWgt);
846 else h.m_hb_sig3DNtr->Fill(Signif3D, evtWgt);
847 }
848
849
851 if(!tmpVertex){
852 return finalVertices;
853 }
854 tmpVertex->makePrivateStore();
856
857 std::vector<float> floatErrMtx;
858 floatErrMtx.resize(6);
859 for(
int i = 0;
i<6;
i++) floatErrMtx[i] = vrt.vertexCov[i];
861
863
864 std::vector<Trk::VxTrackAtVertex> & tmpVTAV = tmpVertex->
vxTrackAtVertex();
865 tmpVTAV.clear();
866
867 for(int ii = 0; ii<nth; ii++) {
869 CovMtxP.setIdentity();
871 vrt.trkAtVrt[ii][0],
872 vrt.trkAtVrt[ii][1],
873 vrt.trkAtVrt[ii][2],
874 Trk::PerigeeSurface(vrt.vertex),
875 std::move(CovMtxP) );
876 tmpVTAV.emplace_back( 1., tmpMeasPer );
877
879 ElementLink<xAOD::TrackParticleContainer> TEL;
884 }
885 }
886
887 finalVertices.push_back(tmpVertex);
888 ngoodVertices_final++;
889 if(nth==1) n1trVrt++;
890
892 isPrimaryVertex = false;
893 continue;
894 }
895
896 vertexMom += vrt.vertexMom;
897
898 }
899
902 h.m_hb_goodvrtN->Fill(ngoodVertices_final+0.1, evtWgt);
903 h.m_curTup->ewgt = evtWgt;
904 if(n1trVrt)
h.m_hb_goodvrtN->Fill(n1trVrt+15., evtWgt);
906 }
907
908 if(ngoodVertices_final==0){
909 return finalVertices;
910 }
911
912
913
914
915
916 double totMass = vertexMom.M();
918 double eRatio = vertexMom.E()/momentumJet.E();
919 results.push_back(eRatio<1. ? eRatio : 1.);
920 results.push_back(vrt2TrackNumber);
922 if(xAODwrk)
results.push_back(xAODwrk->listSecondTracks.size());
924 results.push_back(momentumJet.E());
925
928 h.m_hb_ratio->Fill( results[1], evtWgt);
929 h.m_hb_totmass->Fill( results[0], evtWgt);
930 h.m_hb_nvrt2->Fill( results[2], evtWgt);
931 h.m_hb_mom->Fill( momentumJet.Perp(), evtWgt);
932 }
933
934 return finalVertices;
935
936}
#define AmgSymMatrix(dim)
std::vector< Trk::VxTrackAtVertex > & vxTrackAtVertex()
Non-const access to the VxTrackAtVertex vector.
float j(const xAOD::IParticle &, const xAOD::TrackMeasurementValidation &hit, const Eigen::Matrix3d &jab_inv)
void swap(ElementLinkVector< DOBJ > &lhs, ElementLinkVector< DOBJ > &rhs)
1.16.1