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->TracksForFit,
112 primVrt,
113 jetDir,
114 inpMass,
115 nRefPVTrk,
116 xAODwrk->TrkFromV0,
117 xAODwrk->listSecondTracks,
118 compatibilityGraph);
119 if (xAODwrk->TrkFromV0.size() > 1) {
122 }
123 vrt2TrackNumber = xAODwrk->listSecondTracks.size()/2.0;
126 }
127
128
129 ATH_MSG_DEBUG(
" Found different tracks in pairs=" << vrt2TrackNumber);
130
131
132
133
134
135
136
137
138
139
140
141 std::unique_ptr<std::vector<WrkVrt>> wrkVrtSet = std::make_unique<std::vector<WrkVrt>>();
144 std::unique_ptr<Trk::IVKalState> state =
m_fitSvc->makeState();
146
147
148
149 std::vector<std::vector<int>> allCliques;
150 bron_kerbosch_all_cliques(compatibilityGraph, clique_visitor(allCliques));
151
152 for (const auto& clique : allCliques) {
153
154 newvrt.selTrk.clear();
155 if (xAODwrk) xAODwrk->tmpListTracks.clear();
156
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));
160 }
161
163 if (
sc.isFailure() ||
165
166 m_fitSvc->setApproximateVertex(primVrt.
x(),
169 *state);
171 m_fitSvc->setApproximateVertex(0., 0., 0., *state);
172 } else {
174 double jetVrtDir = jetDir.Px() * vDist.x() + jetDir.Py() * vDist.y() +
175 jetDir.Pz() * vDist.z();
177 if (jetVrtDir > 0.) {
178
179 m_fitSvc->setApproximateVertex(newvrt.vertex.x(),
180 newvrt.vertex.y(),
181 newvrt.vertex.z(),
182 *state);
183 } else {
184 m_fitSvc->setApproximateVertex(primVrt.
x(), primVrt.
y(), primVrt.
z(), *state);
185 }
186 }
187
188 sc = StatusCode::FAILURE;
189 if (xAODwrk) {
191 newvrt.vertex,
192 newvrt.vertexMom,
193 newvrt.vertexCharge,
194 newvrt.vertexCov,
195 newvrt.chi2PerTrk,
196 newvrt.trkAtVrt,
197 newvrt.chi2,
198 *state,
199 false);
200 }
201 if (
sc.isFailure())
continue;
202
203 if (clique.size() == 2 && newvrt.chi2 > 10.) continue;
204
205 if (newvrt.chi2PerTrk.size() == 2){
206 newvrt.chi2PerTrk[0] = newvrt.chi2PerTrk[1] = newvrt.chi2 / 2.;
207 }
208 newvrt.Good = true;
209 newvrt.nCloseVrt = 0;
210 newvrt.dCloseVrt = 1000000.;
211 newvrt.projectedVrt =
213 wrkVrtSet->push_back(newvrt);
214 }
215
216
217
218
219 bool disassembled = false;
220
221 do {
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);
229 disassembled = true;
230 }
231 }
232 }while(disassembled);
233
234
235 while( (*wrkVrtSet).size()>1 ){
236 int tmpN = (*wrkVrtSet).size();
237
238 int iv = 0;
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();
243
244 int jv = iv+1;
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();
249
250 int nTCom =
nTrkCommon(wrkVrtSet.get(), iv, jv);
251 if(nTCom==ntrk_i){
252 (*wrkVrtSet).erase((*wrkVrtSet).begin()+iv);
253 break;
254 }
255 else if(nTCom==ntrk_j){
256 (*wrkVrtSet).erase((*wrkVrtSet).begin()+jv);
257 break;
258 }
259 }
260 if(jv!=tmpN) break;
261
262 }
263 if(iv==tmpN-1) break;
264
265 }
266
267
268
269
270 std::multimap<int,std::pair<int,int>> vrtWithCommonTrk;
271 std::multimap<int,std::pair<int,int>>::reverse_iterator icvrt;
272 do{
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);
279 if(!nTCom)continue;
280 vrtWithCommonTrk.emplace(nTCom,std::make_pair(iv,jv));
281 }
282 }
283
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();
290 double probV = -1.;
291 if (xAODwrk) {
293 }
294
295 if(probV<probVrtMergeLimit){
296 if(nTrkI==2 || nTrkJ==2 || nTCom<2) {
297 continue;
298 }
299 if(nTCom>nTrkI-nTCom || nTCom>nTrkJ-nTCom){
300
302 break;
303 }
304 continue;
305 }
306
307 newvrt.Good = true;
308 (*wrkVrtSet).push_back(newvrt);
309 (*wrkVrtSet)[iv].Good = false;
310 (*wrkVrtSet)[jv].Good = false;
311 break;
312 }
313
314 } while( icvrt != vrtWithCommonTrk.rend() );
315
316
318 int cvgood=0;
319 for(const auto& vrt : (*wrkVrtSet)){
320 if(vrt.Good) cvgood++;
321 }
323 h.m_hb_rawVrtN->Fill(
static_cast<float>(cvgood), evtWgt);
324 }
325
326
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.);
333 }
334
335 unsigned int tmpV = 0;
336 while( tmpV<(*wrkVrtSet).size() ){
337 if( !(*wrkVrtSet)[tmpV].Good ) (*wrkVrtSet).erase((*wrkVrtSet).begin()+tmpV);
338 else tmpV++;
339 }
340
341 if((*wrkVrtSet).empty()){
342 return finalVertices;
343 }
344
345 std::vector< std::vector<float> > trkScore(0);
347 for(
auto &trk : xAODwrk->listJetTracks) trkScore.push_back(
m_trackClassificator->trkTypeWgts(trk, primVrt, jetDir));
348 }
349
350 for(
auto &tmpV : (*wrkVrtSet) ) tmpV.projectedVrt=
jetProjDist(tmpV.vertex, primVrt, jetDir);
351
352
353
354
355
356
357 std::unique_ptr<std::vector< std::deque<long int> > > trkInVrt = std::make_unique<std::vector<std::deque<long int>>>(nTracks);
359
360 double foundMaxT;
361 long int selectedTrack, selectedVertex;
362 int foundV1, foundV2;
363
365 while( ( foundMaxT =
MaxOfShared(wrkVrtSet.get(), trkInVrt.get(), selectedTrack, selectedVertex) ) > 0) {
366
367 double foundMinVrtDst = 1000000.;
369
370
372 &&
nTrkCommon(wrkVrtSet.get(), foundV1, foundV2) ){
373
374 bool vrtMerged=false;
376 if(foundV1<foundV2)
std::swap(foundV1, foundV2);
377
378 double probV = 0.;
379 if (xAODwrk) probV =
mergeAndRefitVertices(wrkVrtSet.get(), foundV1, foundV2, newvrt, xAODwrk->listJetTracks, *state);
380
381 if(probV>probVrtMergeLimit && newvrt.vertexMom.M()<c_vrtBCMassLimit){
382
383 double tstDst =
jetProjDist(newvrt.vertex, primVrt, jetDir);
384 if(tstDst>0.){
385
386 std::swap((*wrkVrtSet)[foundV1], newvrt);
387 (*wrkVrtSet)[foundV1].projectedVrt = tstDst;
388 (*wrkVrtSet)[foundV2].Good = false;
389 (*wrkVrtSet)[foundV2].selTrk.clear();
390 vrtMerged=true;
391 }
392 }
393
394 (*wrkVrtSet)[foundV1].nCloseVrt = -1;
395 (*wrkVrtSet)[foundV1].dCloseVrt = 1000000.;
396 (*wrkVrtSet)[foundV2].nCloseVrt = -1;
397 (*wrkVrtSet)[foundV2].dCloseVrt = 1000000.;
399
400 }
401
402 if(vrtMerged){
403 trkInVrt->resize(nTracks);
405 continue;
406 }
407
408 }
409
411
412 sc = StatusCode::FAILURE;
413 if(xAODwrk)
sc =
refitVertex( wrkVrtSet.get(), selectedVertex, xAODwrk->listJetTracks, *state,
false);
414
415 (*wrkVrtSet)[selectedVertex].projectedVrt =
jetProjDist((*wrkVrtSet)[selectedVertex].vertex, primVrt, jetDir);
416
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;
421 }
422
423 }
424
425
426
427
428
429 double minDistVV =
minVrtVrtDist( wrkVrtSet.get(), foundV1, foundV2);
431 if(foundV1<foundV2)
std::swap(foundV1, foundV2);
432
433 double probV = 0.;
434 if(xAODwrk) probV =
mergeAndRefitVertices(wrkVrtSet.get(), foundV1, foundV2, newvrt, xAODwrk->listJetTracks, *state);
435
436 if(probV>probVrtMergeLimit && newvrt.vertexMom.M()<c_vrtBCMassLimit){
437
438 double tstDst =
jetProjDist(newvrt.vertex, primVrt, jetDir);
439 if(tstDst>0.){
440
442 (*wrkVrtSet)[foundV1].projectedVrt = tstDst;
443 (*wrkVrtSet)[foundV2].Good = false;
444 (*wrkVrtSet)[foundV2].selTrk.clear();
445 }
446 }
447
448 (*wrkVrtSet)[foundV1].nCloseVrt = -1;
449 (*wrkVrtSet)[foundV1].dCloseVrt = 1000000.;
450 (*wrkVrtSet)[foundV2].nCloseVrt = -1;
451 (*wrkVrtSet)[foundV2].dCloseVrt = 1000000.;
453 }
454
455
456
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;
461
462 double tmpProb = TMath::Prob( vert_i.chi2, 2*vert_i.selTrk.size()-3 );
463 if(tmpProb<0.001){
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);
467 }
468 }
469
470
472
473
474
475 for(auto &ntrVrt : (*wrkVrtSet)){
476 if(!ntrVrt.Good || ntrVrt.selTrk.size()<=1) continue;
477
478 for(auto &onetVrt : (*wrkVrtSet)){
479 if(!onetVrt.Good || onetVrt.selTrk.size()!=1) continue;
480
481 if(ntrVrt.detachedTrack==onetVrt.selTrk[0]){
483 newV.selTrk.push_back(ntrVrt.detachedTrack);
484 double vProb = 0.;
485 if(xAODwrk) vProb =
refitVertex(newV, xAODwrk->listJetTracks, *state,
true);
486 if(vProb>probVrtMergeLimit){
487 onetVrt.Good=false;
488 ntrVrt=newV;
489 ntrVrt.detachedTrack=-1;
490 }
491 break;
492 }
493 }
494 }
495
496
497 for(auto &iVrt : (*wrkVrtSet)){
498 if(!iVrt.Good || iVrt.selTrk.size()!=1) continue;
499
500 for(auto &jVrt : (*wrkVrtSet)){
501 if(!jVrt.Good || jVrt.selTrk.size()!=1) continue;
502
503 if(iVrt.detachedTrack==jVrt.selTrk[0]){
505 newV.selTrk.push_back(iVrt.detachedTrack);
506 double vProb = 0.;
507 if(xAODwrk) vProb =
refitVertex(newV, xAODwrk->listJetTracks, *state,
true);
508 if(vProb>probVrtMergeLimit){
509 jVrt.Good = false;
510 iVrt = newV;
511 iVrt.detachedTrack = -1;
512 }
513 break;
514 }
515 }
516
517 }
518
519
520 for(auto &ntrVrt : (*wrkVrtSet)){
521 if(!ntrVrt.Good || ntrVrt.selTrk.size()<=1) continue;
522
523 for(auto tr : ntrVrt.selTrk){
524 for(auto &onetVrt : (*wrkVrtSet)){
525 if(!onetVrt.Good || onetVrt.selTrk.size()!=1) continue;
526
527 if(onetVrt.detachedTrack==tr){
529 newV.selTrk.push_back(onetVrt.selTrk[0]);
530 double vProb = 0.;
531 if(xAODwrk) vProb =
refitVertex( newV, xAODwrk->listJetTracks, *state,
true);
532 if(vProb>probVrtMergeLimit){
533 onetVrt.Good = false;
534 ntrVrt = newV;
535 }
536 }
537
538 }
539 }
540
541 }
542
543
544 for(auto & curVrt : (*wrkVrtSet) ) {
545 if(!curVrt.Good ) continue;
546 if( std::abs(curVrt.vertex.z())>650. ){
547 curVrt.Good = false;
548 continue;
549 }
550 if(curVrt.selTrk.size() != 1) continue;
551
552 curVrt.Good = false;
553
555
556 double Signif3D = 0.;
557 vrtVrtDist(primVrt,curVrt.vertex, curVrt.vertexCov, Signif3D);
558 double tmpProb = TMath::Prob(curVrt.chi2, 1);
559
560 bool trkGood = false;
561 if(xAODwrk) trkGood =
check1TrVertexInPixel(xAODwrk->listJetTracks[curVrt.selTrk[0]], curVrt.vertex, curVrt.vertexCov);
562
563 if(trkGood && tmpProb>0.01){
564
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);
568
571 h.m_hb_diffPS->Fill(Signif3DP, evtWgt);
572 }
573
575 }
576
577 }
578
579 }
580
581
582 for(auto& curVrt : (*wrkVrtSet)){
583 if(!curVrt.Good ) continue;
584 int nth = curVrt.selTrk.size();
585 if(nth==1) continue;
586
587 double Signif3D = 0.;
588 vrtVrtDist(primVrt,curVrt.vertex, curVrt.vertexCov, Signif3D);
589 if(xAODwrk) xAODwrk->tmpListTracks.resize(nth);
590
591 for(
int i = 0;
i < nth;
i++) {
592 if(xAODwrk) xAODwrk->tmpListTracks[
i] = xAODwrk->listJetTracks[curVrt.selTrk[
i]];
593 }
594 curVrt.Good = false;
595
596 if(nth <= 1) continue;
597 if(curVrt.projectedVrt<0.) continue;
598
599 if(TMath::Prob( curVrt.chi2, 2*nth-3)<0.001) continue;
600
601 if(nth==2 && xAODwrk){
602
604 if(!
check2TrVertexInPixel(xAODwrk->tmpListTracks[0],xAODwrk->tmpListTracks[1],curVrt.vertex,curVrt.vertexCov))
continue;
605 }
606
607
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;
614 }
615
616 }
617
618
621 h.m_hb_r2d->Fill(curVrt.vertex.perp(), evtWgt);
622 }
623
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);
631 }
632
635 h.m_hb_r2d->Fill(curVrt.vertex.perp(), evtWgt);
636 }
637
638
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);
648 }
649
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;
653 }
654
657 h.m_hb_sig3DTot->Fill(Signif3D, evtWgt);
658 }
660
661 curVrt.Good = true;
662
663 }
664
665
667
668
669 int n2trVrt = 0;
670 int nNtrVrt = 0;
671 int ngoodVertices = 0;
672 std::vector<WrkVrt> goodVertices(0);
673
674 for(auto & iv : (*wrkVrtSet)) {
675 int nth = iv.selTrk.size();
676 if(nth == 0) continue;
677
678 if(iv.Good){
679 ngoodVertices++;
680 goodVertices.emplace_back(iv);
681 if(nth==2) n2trVrt++;
682 if(nth>=3) nNtrVrt++;
683 }
684 }
685
686 if(ngoodVertices == 0 || (n2trVrt+nNtrVrt)==0){
687 return finalVertices;
688 }
689
690
691 while(1){
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] );
696 swapDone = true;
697 }
698 }
699 if(!swapDone) break;
700 }
701
702 while(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] );
708 swapDone = true;
709 }
710 }
711 if(!swapDone) break;
712 }
713
714 if(ngoodVertices>1){
715 if( goodVertices[1].vertexMom.M()-goodVertices[0].vertexMom.M() > 5000.){
716 std::swap( goodVertices[0], goodVertices[1] );
717 }
718 }
719
722 h.m_hb_distVV->Fill(
minVrtVrtDist(wrkVrtSet.get(), foundV1, foundV2), evtWgt);
723 }
724
725
726
727
728
729
730
731 std::vector<int> nonusedTrk(0);
732 for(int trk=0; trk<nTracks; trk++){
733 bool present = false;
734
735 for(const auto& iv : (*wrkVrtSet)){
736 if(iv.Good){
737 for(auto t_in_V : iv.selTrk){
738 if(trk==t_in_V){
739 present = true;
740 break;
741 }
742 }
743 }
744 if(present) break;
745 }
746
747 if(!present) nonusedTrk.push_back(trk);
748 }
749
750 struct MatchedSV{int indVrt; double Signif3D;};
751 std::vector<MatchedSV> matchSV(0);
752
753 for(const auto& i_trk : nonusedTrk){
754 MatchedSV tmpV = {0, 1.e9};
755
756 for(unsigned int iv = 0; iv<goodVertices.size(); iv++){
757
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;
761
762 double Signif = 0.;
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;
767 tmpV.indVrt=iv;
768 }
769 }
770
771 matchSV.push_back(tmpV);
772 }
773
774 for(int iv = 0; iv<static_cast<int>(goodVertices.size()); iv++){
775 WrkVrt newV(goodVertices[iv]);
776 bool addedT = false;
777 std::map<double,int> addTrk;
778
779 for(
unsigned int it = 0;
it<nonusedTrk.size();
it++){
780 if(matchSV[it].indVrt==iv) addTrk[matchSV[
it].Signif3D]=
it;
781 }
782
783 std::map<double,int>::iterator atrk=addTrk.begin();
784 if(!addTrk.empty()){
785 if(atrk->first<4.){
786 newV.selTrk.push_back(nonusedTrk[atrk->second]);
787 addedT=true;
788 }
789 }
790
791 if(addTrk.size()>1){
792 ++atrk;
793 if(atrk->first<4.){
794 newV.selTrk.push_back(nonusedTrk[atrk->second]);
795 }
796 }
797
798 if(addedT){
799 double vProb = 0.;
800 if(xAODwrk) vProb =
refitVertex(newV, xAODwrk->listJetTracks, *state,
true);
801 if(vProb>0.01) goodVertices[iv] = newV;
802 else{
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];
806 }
807 }
808
809 }
810
811
812
813
814
815 int ngoodVertices_final = 0;
816 int n1trVrt = 0;
817 TLorentzVector vertexMom;
818 bool isPrimaryVertex = true;
819
820 for(auto& vrt : goodVertices){
821 int nth = vrt.selTrk.size();
822 if(xAODwrk) xAODwrk->tmpListTracks.clear();
823
824 for(
int i = 0;
i < nth;
i++) {
825 int j = vrt.selTrk[
i];
826 if(xAODwrk) xAODwrk->tmpListTracks.push_back( xAODwrk->listJetTracks[j] );
827 }
828
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);
836
837 if(nth==2){
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);
842
843 }
844 else if(vrt.vertexMom.M()>6000.)
h.m_hb_sig3D1tr->Fill(Signif3D, evtWgt);
845 else h.m_hb_sig3DNtr->Fill(Signif3D, evtWgt);
846 }
847
848
850 if(!tmpVertex){
851 return finalVertices;
852 }
855
856 std::vector<float> floatErrMtx;
857 floatErrMtx.resize(6);
858 for(
int i = 0;
i<6;
i++) floatErrMtx[i] = vrt.vertexCov[i];
860
862
863 std::vector<Trk::VxTrackAtVertex> & tmpVTAV = tmpVertex->
vxTrackAtVertex();
864 tmpVTAV.clear();
865
866 for(int ii = 0; ii<nth; ii++) {
868 CovMtxP.setIdentity();
870 vrt.trkAtVrt[ii][0],
871 vrt.trkAtVrt[ii][1],
872 vrt.trkAtVrt[ii][2],
873 Trk::PerigeeSurface(vrt.vertex),
874 std::move(CovMtxP) );
875 tmpVTAV.emplace_back( 1., tmpMeasPer );
876
878 ElementLink<xAOD::TrackParticleContainer> TEL;
883 }
884 }
885
886 finalVertices.push_back(tmpVertex);
887 ngoodVertices_final++;
888 if(nth==1) n1trVrt++;
889
891 isPrimaryVertex = false;
892 continue;
893 }
894
895 vertexMom += vrt.vertexMom;
896
897 }
898
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);
905 }
906
907 if(ngoodVertices_final==0){
908 return finalVertices;
909 }
910
911
912
913
914
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());
924
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);
931 }
932
933 return finalVertices;
934
935}
#define AmgSymMatrix(dim)
std::vector< Trk::VxTrackAtVertex > & vxTrackAtVertex()
Non-const access to the VxTrackAtVertex vector.
void swap(ElementLinkVector< DOBJ > &lhs, ElementLinkVector< DOBJ > &rhs)
1.14.0