308 {
309 int nTrkToVertex(0);
311
313 for (const Tracklet &trk : tracklets) aveTrkPos += trk.globalPosition();
314 aveTrkPos /= tracklets.size();
315
316
317 double avePhi = aveTrkPos.phi();
318 double LoF = std::atan2(aveTrkPos.perp(), aveTrkPos.z());
320
321
322 std::vector<double> Rpos;
324 double LoFdist = std::abs(RadialDist / std::sin(LoF));
326 double PlaneSpacing = std::abs(
m_VxPlaneDist / std::cos(LoF));
328
329
330 std::array< std::vector<std::unique_ptr<Trk::TrackParameters>>, MAXPLANES> TracksForVertexing{};
331 std::array< std::vector<std::unique_ptr<Trk::TrackParameters>>, MAXPLANES> TracksForErrors{};
332 std::array< std::vector<bool>, MAXPLANES> isNeutralTrack{};
333
334 for (const Tracklet &trk : tracklets) {
336 ++nTrkToVertex;
337
338 Amg::Vector3D trkgpos(trk.globalPosition().perp() * std::cos(avePhi),
339 trk.globalPosition().perp() * std::sin(avePhi),
340 trk.globalPosition().z());
341 double x0 = trkgpos.x();
342 double y0 = trkgpos.y();
343 double r0 = trkgpos.perp();
344
345
347 double NominalTrkAng = anglesign * NominalAngle;
348 double MaxTrkAng = anglesign * RotationAngle;
349
350
351 for (
int k = 0;
k < nplanes; ++
k) {
352
353 if (Rpos[k] > trk.globalPosition().perp()) break;
354
355
356 double Xp = Rpos[
k] * std::cos(avePhi);
357 double Yp = Rpos[
k] * std::sin(avePhi);
358
359
360 double DelR = std::hypot(x0 - Xp, y0 - Yp) / std::cos(NominalAngle);
361 double X1 = DelR * std::cos(NominalTrkAng + avePhi) + Xp;
362 double Y1 = DelR * std::sin(NominalTrkAng + avePhi) + Yp;
363 double R1 = std::hypot(X1, Y1);
367 double Dirmag = std::hypot(X1 - Xp, Y1 - Yp);
368 double Xdir = (X1 - Xp) / Dirmag;
369 double Ydir = (Y1 - Yp) / Dirmag;
370 double trkpx = Xdir * trk.momentum().perp();
371 double trkpy = Ydir * trk.momentum().perp();
372 double trkpz = trk.momentum().z();
373
374
375 double charge = trk.charge();
376 if (std::abs(
charge) < 0.1) {
378 isNeutralTrack[
k].push_back(
true);
379 } else
380 isNeutralTrack[
k].push_back(
false);
381
382
386 TracksForVertexing[k].push_back(std::make_unique<Trk::
Perigee>(0., 0., trkmomentum.
phi(), trkmomentum.
theta(),
charge / trkmomentum.
mag(), Trk::PerigeeSurface(trkgpos), covariance));
387
388
389 double xp = Rpos[k] * std::cos(avePhi);
390 double yp = Rpos[k] * std::sin(avePhi);
391 double delR = std::hypot(x0 - xp, y0 - yp) / std::cos(RotationAngle);
392 double x1 = delR * std::cos(MaxTrkAng + avePhi) + xp;
393 double y1 = delR * std::sin(MaxTrkAng + avePhi) + yp;
394 double r1 = std::hypot(x1, y1);
395 double norm = r0 / r1;
396 x1 = x1 * norm;
397 y1 = y1 * norm;
398 double dirmag = std::hypot(x1 - xp, y1 - yp);
399 double xdir = (x1 - xp) / dirmag;
400 double ydir = (y1 - yp) / dirmag;
401 double errpx = xdir * trk.momentum().
perp();
402 double errpy = ydir * trk.momentum().
perp();
403 double errpz = trk.momentum().
z();
404
405
407 Amg::
Vector3D trkerrmom(errpx, errpy, errpz);
408 Amg::
Vector3D trkerrpos(x1, y1, trk.globalPosition().
z());
409 TracksForErrors[k].push_back(std::make_unique<Trk::
Perigee>(0., 0., trkerrmom.
phi(), trkerrmom.
theta(),
charge / trkerrmom.
mag(), Trk::PerigeeSurface(trkerrpos), covariance2));
410 }
411 }
412
413
414 if (nTrkToVertex < 3) return;
415
416
417 bool boundaryCheck = true;
418 std::array<std::vector<double>, MAXPLANES> ExtrapZ{};
419 std::array<std::vector<double>, MAXPLANES> dlength{};
420 std::array<std::vector<std::pair<unsigned int, unsigned int>>, MAXPLANES> UsedTracks{};
421 std::array<std::vector<bool>, MAXPLANES> ExtrapSuc{};
422 std::vector<std::unique_ptr<MSVertex>> vertices; vertices.reserve(nplanes);
423 std::array<std::vector<double>, MAXPLANES>
sigmaZ{};
424 std::array<std::vector<Amg::Vector3D>, MAXPLANES> pAtVx{};
425
426
427 for (
int k = 0;
k < nplanes; ++
k) {
428 double rpos = Rpos[
k];
429 for (
unsigned int i = 0;
i < TracksForVertexing[
k].size(); ++
i) {
430
431 if (TracksForVertexing[k].size() < 3) break;
432
434 surfaceTransformMatrix.setIdentity();
435 Trk::CylinderSurface cyl(surfaceTransformMatrix, rpos, 10000.);
436
437 std::unique_ptr<const Trk::TrackParameters> extrap_par(
440
442
443 if (extrap) {
444
445 if (isNeutralTrack[k].at(i)) {
446 double pTot = std::hypot(TracksForVertexing[k].at(i)->
momentum().
perp(), TracksForVertexing[k].at(i)->
momentum().
z());
448 double extrapRdist = TracksForVertexing[
k].at(i)->position().perp() - Rpos[
k];
449 double sz = std::abs(20 * dirErr * extrapRdist * std::pow(pTot,2) / std::pow(TracksForVertexing[k].at(i)->
momentum().
perp(), 2));
450 double ExtrapErr =
sz;
452 ExtrapSuc[
k].push_back(
false);
453 else {
454 ExtrapSuc[
k].push_back(
true);
455 std::pair<unsigned int, unsigned int> trkmap(ExtrapZ[k].size(), i);
456 UsedTracks[
k].push_back(trkmap);
460 dlength[
k].push_back(0);
461 }
462 }
463
464 else {
465
467 srfTransMat2.setIdentity();
468 Trk::CylinderSurface cyl2(srfTransMat2, rpos, 10000.);
469 std::unique_ptr<const Trk::TrackParameters> extrap_par2(
473
474 if (extrap2) {
477 double ExtrapErr = std::hypot(
sz, zdiff);
479 ExtrapSuc[
k].push_back(
false);
480 else {
481
482 ExtrapSuc[
k].push_back(
true);
483 std::pair<unsigned int, unsigned int> trkmap(ExtrapZ[k].size(), i);
484 UsedTracks[
k].push_back(trkmap);
488 dlength[
k].push_back(zdiff);
489 }
490 } else
491 ExtrapSuc[
k].push_back(
false);
492 }
493 }
494
495 else
496 ExtrapSuc[
k].push_back(
false);
497 }
498 }
499
500
501
502 std::array<std::vector<Amg::Vector3D>, MAXPLANES> trkp{};
503
504 for (
int k = 0;
k < nplanes; ++
k) {
505 if (ExtrapZ[k].size() < 3) continue;
506
507 double zLoF = Rpos[
k] / std::tan(LoF);
508 double dzLoF(10);
509 double aveZpos(0), posWeight(0);
510 for (
unsigned int i = 0;
i < ExtrapZ[
k].size(); ++
i) {
511 double ExtrapErr = std::hypot(sigmaZ[k][i], dlength[k][i], dzLoF);
512 if (isNeutralTrack[k][i]) ExtrapErr = std::hypot(sigmaZ[k][i], dzLoF);
513 aveZpos += ExtrapZ[
k][
i] / std::pow(ExtrapErr,2);
514 posWeight += 1. / std::pow(ExtrapErr,2);
515 }
516
517 zLoF = aveZpos / posWeight;
518 double zpossigma(dzLoF), Chi2(0), Chi2Prob(-1);
519 unsigned int Nitr(0);
520 std::vector<unsigned int> vxtracks;
521 std::vector<bool> blocklist(ExtrapZ[k].size(), false);
522
523
524 while (true) {
525 vxtracks.clear(); trkp[
k].clear();
526 int tmpnTrks(0);
527 double tmpzLoF(0), tmpzpossigma(0), tmpchi2(0), posWeight(0), worstdelz(0);
528 unsigned int iworst(0);
529
530 for (
unsigned int i = 0;
i < ExtrapZ[
k].size(); ++
i) {
531 if (blocklist[i]) continue;
532 trkp[
k].push_back(pAtVx[k][i]);
533 double delz = zLoF - ExtrapZ[
k][
i];
534 double ExtrapErr = std::hypot(sigmaZ[k][i], dlength[k][i], dzLoF);
535 double trkchi2 = std::pow(delz,2) / std::pow(ExtrapErr,2);
536 if (trkchi2 > worstdelz) {
538 worstdelz = trkchi2;
539 }
540 tmpzLoF += ExtrapZ[
k][
i] / std::pow(ExtrapErr,2);
541 posWeight += 1. / std::pow(ExtrapErr,2);
542 tmpzpossigma += std::pow(delz,2);
543 tmpchi2 += trkchi2;
544 ++tmpnTrks;
545 }
546
547 if (tmpnTrks < 3) break;
548 tmpzpossigma = std::sqrt(tmpzpossigma / (double)tmpnTrks);
549 zLoF = tmpzLoF / posWeight;
550 zpossigma = tmpzpossigma;
551 double testChi2 = TMath::Prob(tmpchi2, tmpnTrks - 1);
553 blocklist[iworst] = true;
554 else {
555 Chi2 = tmpchi2;
556 Chi2Prob = testChi2;
557
558 for (
unsigned int i = 0;
i < ExtrapZ[
k].size(); ++
i) {
559 double delz = zLoF - ExtrapZ[
k][
i];
560 double ExtrapErr = std::hypot(sigmaZ[k][i], dlength[k][i], dzLoF);
561 double trkErr = std::hypot(ExtrapErr, zpossigma) + 0.001;
562 double trkNsigma = std::abs(delz / trkErr);
563 if (trkNsigma < 3) vxtracks.push_back(i);
564 }
565 break;
566 }
567 if (Nitr >= (ExtrapZ[k].size() - 3)) break;
568 ++Nitr;
569 }
570
571 if (vxtracks.size() < 3) continue;
572
573
574 std::vector<const xAOD::TrackParticle*> vxTrackParticles;
575 vxTrackParticles.reserve(vxtracks.size());
576 for (std::vector<unsigned int>::iterator vxtrk = vxtracks.begin(); vxtrk != vxtracks.end(); ++vxtrk) {
577 for (
unsigned int i = 0;
i < UsedTracks[
k].size(); ++
i) {
578 if ((*vxtrk) != UsedTracks[k].at(i).first) continue;
579 const Tracklet& trklt = tracklets.at(UsedTracks[k].at(i).second);
581 break;
582 }
583 }
584 Amg::Vector3D position(Rpos[k] * std::cos(avePhi), Rpos[k] * std::sin(avePhi), zLoF);
585 vertices.push_back(std::make_unique<MSVertex>(1, position, vxTrackParticles, Chi2Prob, Chi2, 0, 0, 0));
586 }
587
588 if (vertices.empty()) return;
589
590
591
592 unsigned int bestVx(0);
593 for (
unsigned int k = 1;
k < vertices.size(); ++
k) {
594 if (vertices[k]->getChi2Probability() <
m_VxChi2ProbCUT || vertices[k]->getNTracks() < 3)
continue;
595 if (vertices[k]->getNTracks() < vertices[bestVx]->getNTracks()) continue;
596 if (vertices[k]->getNTracks() == vertices[bestVx]->getNTracks() &&
597 vertices[k]->getChi2Probability() < vertices[bestVx]->getChi2Probability())
598 continue;
600 }
601 vtx = std::make_unique<MSVertex>(*vertices[bestVx]);
602 vertices.clear();
603 }
Scalar perp() const
perp method - perpendicular length
Scalar phi() const
phi method
Scalar theta() const
theta method
Scalar mag() const
mag method
double charge(const T &p)
#define AmgSymMatrix(dim)
const xAOD::TrackParticle * getTrackParticle() const
const Amg::Vector3D & momentum() const
Access method for the momentum.
Amg::Vector2D localPosition() const
Access method for the local coordinates, local parameter definitions differ for each surface type.
void Norm(TH1 *h, double scale)
double error(const Amg::MatrixX &mat, int index)
return diagonal error of the matrix caller should ensure the matrix is symmetric and the index is in ...
Eigen::Affine3d Transform3D
const double r0
electron radius{cm}
ParametersT< TrackParametersDim, Charged, CylinderSurface > AtaCylinder
1.14.0