CbmDcaVertexFinder.cxx

Go to the documentation of this file.

/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
   SPDX-License-Identifier: GPL-3.0-only
   Authors: Dario Ramirez [committer] */
 
#include "CbmDcaVertexFinder.h"
 
#include "CbmGlobalTrack.h"
#include "CbmVertex.h"
 
#include <Logger.h>
 
CbmDcaVertexFinder::CbmDcaVertexFinder()
  : fMaxDca(0.1)
  , fInputTracks(std::vector<CbmGlobalTrack*>())
  , fNbPairs(0)
  , fQA(std::nullopt)
{
  LOG(debug) << "Creating CbmDcaVertexFinder ...";
  fCovMatrix.ResizeTo(3, 3);
}
 
CbmDcaVertexFinder::CbmDcaVertexFinder(double max_dca)
  : fMaxDca(max_dca)
  , fInputTracks(std::vector<CbmGlobalTrack*>())
  , fNbPairs(0)
  , fQA(std::nullopt)
{
  LOG(debug) << "Creating CbmDcaVertexFinder ...";
  fCovMatrix.ResizeTo(3, 3);
}

CbmDcaVertexFinder::CbmDcaVertexFinder(const std::vector<CbmGlobalTrack*> tracks)
  : fMaxDca(0.1)
  , fInputTracks(tracks)
  , fNbPairs(0)
  , fQA(std::nullopt)
{
  LOG(debug) << "Creating CbmDcaVertexFinder ...";
  fCovMatrix.ResizeTo(3, 3);
}

CbmDcaVertexFinder::CbmDcaVertexFinder(const std::vector<CbmGlobalTrack*> tracks, const double max_dca)
  : fMaxDca(max_dca)
  , fInputTracks(tracks)
  , fNbPairs(0)
  , fQA(std::nullopt)
{
  LOG(debug) << "Creating CbmDcaVertexFinder ...";
  fCovMatrix.ResizeTo(3, 3);
}
 
void CbmDcaVertexFinder::SetTracks(const std::vector<CbmGlobalTrack*> tracks)
{
  fInputTracks.clear();
  fInputTracks = tracks;
}

std::optional<CbmVertex> CbmDcaVertexFinder::FindVertex()
{
  // Reset the number of track pair used
  fNbPairs     = 0;
  int n_of_trk = fInputTracks.size();
  LOG(debug) << "- PCA - Find event vertex using CbmGlobalTracks: " << n_of_trk;
  TVector3 vtx;
  for (int trk_i_idx = 0; trk_i_idx < n_of_trk - 1; trk_i_idx++) {
    for (int trk_j_idx = trk_i_idx + 1; trk_j_idx < n_of_trk; trk_j_idx++) {
      auto pca = FindPca(fInputTracks[trk_i_idx], fInputTracks[trk_j_idx]);
      if (pca.has_value() && pca->d_trk < fMaxDca) {
        TVector3 pca_i_j = pca->point;
        vtx += pca_i_j;
        fNbPairs++;
 
        if (fQA.has_value()) {
          fQA->pca_y_vs_x->Fill(pca_i_j[0], pca_i_j[1]);
          fQA->pca_x_vs_z->Fill(pca_i_j[2], pca_i_j[0]);
          fQA->pca_y_vs_z->Fill(pca_i_j[2], pca_i_j[1]);
        }
      }
    }
  }
  if (!fNbPairs) return std::nullopt;
 
  vtx *= 1. / fNbPairs;
 
  // WARNING LINE
  return CbmVertex("EventVertex", "EventVertex", vtx[0], vtx[1], vtx[2], 0, 0, fNbPairs, fCovMatrix);
}
 
std::optional<CbmDcaVertexFinder::PCA> CbmDcaVertexFinder::FindPca(CbmGlobalTrack* trk_i, CbmGlobalTrack* trk_j)
{
  TVector3 r1(trk_i->GetParamFirst()->GetX(), trk_i->GetParamFirst()->GetY(), trk_i->GetParamFirst()->GetZ());
  TVector3 e1(trk_i->GetParamFirst()->GetTx(), trk_i->GetParamFirst()->GetTy(), 1);
 
  TVector3 r2(trk_j->GetParamFirst()->GetX(), trk_j->GetParamFirst()->GetY(), trk_j->GetParamFirst()->GetZ());
  TVector3 e2(trk_j->GetParamFirst()->GetTx(), trk_j->GetParamFirst()->GetTy(), 1);
 
  TVector3 n = e1.Cross(e2);  // Directional vector for segment of the pca
 
  float nn = n.Dot(n);
  if (nn != 0) {  // Non-Parallel tracks
 
    TVector3 e1n = e1.Cross(n);
    TVector3 e2n = e1.Cross(n);
    TVector3 r21 = r2 - r1;
 
    float t1 = e2n.Dot(r21) / nn;
    float t2 = e1n.Dot(r21) / nn;
 
    TVector3 p1  = r1 + t1 * e1;  // Closest point on track 1
    TVector3 p2  = r2 + t2 * e2;  // Closest point on track 2
    TVector3 p21 = p2 - p1;
 
    TVector3 point = 0.5 * (p1 + p2);
    double d_trk   = 0.5 * p21.Mag();
 
    return std::make_optional<PCA>({point, d_trk});
  }
  return std::nullopt;
}