From 61655656faa7e58d239e2ff548a7b0704ee1a907 Mon Sep 17 00:00:00 2001
From: "s.zharko@gsi.de" <s.zharko@gsi.de>
Date: Mon, 9 Oct 2023 02:30:47 +0200
Subject: [PATCH] KF:   introduction of GlobalTracks as source for PV finding
 KF:   introduction of CbmPVFinderKFGlobal mCBM: introduction of the primary
 vertex finder based on the global tracks

UPD: clang-format refresh
---
 macro/mcbm/mcbm_reco_event.C              |   6 +-
 reco/KF/CbmKFPrimaryVertexFinder.cxx      |   2 +-
 reco/KF/Interface/CbmKFTrack.cxx          |  18 ++++
 reco/KF/Interface/CbmKFTrack.h            |   3 +
 reco/KF/Interface/CbmPVFinderKF.cxx       |  13 +--
 reco/KF/Interface/CbmPVFinderKF.h         |  19 ++--
 reco/KF/Interface/CbmPVFinderKFGlobal.cxx | 123 ++++++++++++++++++++++
 reco/KF/Interface/CbmPVFinderKFGlobal.h   |  40 +++++++
 reco/KF/KF.cmake                          |   1 +
 reco/KF/KFLinkDef.h                       |   1 +
 reco/global/CbmFindPrimaryVertex.cxx      |  10 +-
 reco/global/CbmFindPrimaryVertex.h        |  11 +-
 12 files changed, 216 insertions(+), 31 deletions(-)
 create mode 100644 reco/KF/Interface/CbmPVFinderKFGlobal.cxx
 create mode 100644 reco/KF/Interface/CbmPVFinderKFGlobal.h

diff --git a/macro/mcbm/mcbm_reco_event.C b/macro/mcbm/mcbm_reco_event.C
index b5cf9e132d..74e4c5bed3 100644
--- a/macro/mcbm/mcbm_reco_event.C
+++ b/macro/mcbm/mcbm_reco_event.C
@@ -394,12 +394,12 @@ void mcbm_reco_event(Int_t nEvents = 10, TString dataset = "data/test",
 
   //run->AddTask(new CbmKfFitTracksTask(CbmKfFitTracksTask::FitMode::kMcbm));
   // Primary vertex finder from global tracks
-  auto* pvFinder = new CbmPVFinderKF();
-  pvFinder->SetSourceTrackType(CbmPVFinderKF::ESourceTrackType::kGlobalTrack);
+  auto* pvFinder                   = new CbmPVFinderKFGlobal();
   CbmFindPrimaryVertex* findVertex = new CbmFindPrimaryVertex(pvFinder);
+  findVertex->SetTrackType(ECbmDataType::kGlobalTrack);
   run->AddTask(findVertex);
 
-  if (debugWithMC) {  // match tracks 
+  if (debugWithMC) {  // match tracks
     CbmMatchRecoToMC* match2 = new CbmMatchRecoToMC();
     run->AddTask(match2);
   }
diff --git a/reco/KF/CbmKFPrimaryVertexFinder.cxx b/reco/KF/CbmKFPrimaryVertexFinder.cxx
index 8015e993c9..b5dc0fb8ff 100644
--- a/reco/KF/CbmKFPrimaryVertexFinder.cxx
+++ b/reco/KF/CbmKFPrimaryVertexFinder.cxx
@@ -139,7 +139,7 @@ void CbmKFPrimaryVertexFinder::Fit(CbmKFVertexInterface& vtx)
       //* Invert S
       {
         Double_t w = S[0] * S[2] - S[1] * S[1];
-        if (w < 1.E-20) continue;
+        if (w < 1.E-200) continue;
         w           = 1. / w;
         Double_t S0 = S[0];
         S[0]        = w * S[2];
diff --git a/reco/KF/Interface/CbmKFTrack.cxx b/reco/KF/Interface/CbmKFTrack.cxx
index 8bcba6dabf..793dc9ea25 100644
--- a/reco/KF/Interface/CbmKFTrack.cxx
+++ b/reco/KF/Interface/CbmKFTrack.cxx
@@ -4,6 +4,7 @@
 
 #include "CbmKFTrack.h"
 
+#include "CbmGlobalTrack.h"
 #include "CbmKFMath.h"
 #include "CbmStsTrack.h"
 #include "FairTrackParam.h"
@@ -53,6 +54,14 @@ void CbmKFTrack::SetStsTrack(CbmStsTrack& track, bool first)
   GetRefNDF()  = track.GetNDF();
 }
 
+void CbmKFTrack::SetGlobalTrack(CbmGlobalTrack& track, bool first)
+{
+  SetPID(track.GetPidHypo());
+  SetTrackParam(first ? *track.GetParamFirst() : *track.GetParamLast());
+  GetRefChi2() = track.GetChi2();
+  GetRefNDF()  = track.GetNDF();
+}
+
 void CbmKFTrack::GetTrackParam(FairTrackParam& track) { CbmKFMath::CopyTC2TrackParam(&track, fT, fC); }
 
 void CbmKFTrack::GetStsTrack(CbmStsTrack& track, bool first)
@@ -64,6 +73,15 @@ void CbmKFTrack::GetStsTrack(CbmStsTrack& track, bool first)
   track.SetNDF(GetRefNDF());
 }
 
+void CbmKFTrack::GetGlobalTrack(CbmGlobalTrack& track, bool first)
+{
+  FairTrackParam par(first ? *track.GetParamFirst() : *track.GetParamLast());
+  GetTrackParam(par);  //first? *track.GetParamFirst() : *track.GetParamLast() );
+  first ? track.SetParamFirst(&par) : track.SetParamLast(&par);
+  track.SetChi2(GetRefChi2());
+  track.SetNDF(GetRefNDF());
+}
+
 void CbmKFTrack::SetPID(Int_t pidHypo)
 {
   TParticlePDG* particlePDG = TDatabasePDG::Instance()->GetParticle(pidHypo);
diff --git a/reco/KF/Interface/CbmKFTrack.h b/reco/KF/Interface/CbmKFTrack.h
index 02da1d12c8..5c2c9d22a5 100644
--- a/reco/KF/Interface/CbmKFTrack.h
+++ b/reco/KF/Interface/CbmKFTrack.h
@@ -21,6 +21,7 @@
 class CbmKFHit;
 class FairTrackParam;
 class CbmStsTrack;
+class CbmGlobalTrack;
 
 class CbmKFTrack : public CbmKFTrackInterface {
 
@@ -65,9 +66,11 @@ class CbmKFTrack : public CbmKFTrackInterface {
   void SetTrack(CbmKFTrackInterface& track);
   void SetTrackParam(const FairTrackParam& track);
   void SetStsTrack(CbmStsTrack& track, bool first = 1);
+  void SetGlobalTrack(CbmGlobalTrack& track, bool first = 1);
 
   void GetTrackParam(FairTrackParam& track);
   void GetStsTrack(CbmStsTrack& track, bool first = 1);
+  void GetGlobalTrack(CbmGlobalTrack& track, bool first = 1);
 
   void SetPID(Int_t pidHypo);
 
diff --git a/reco/KF/Interface/CbmPVFinderKF.cxx b/reco/KF/Interface/CbmPVFinderKF.cxx
index 958c8f7b71..be7acad497 100644
--- a/reco/KF/Interface/CbmPVFinderKF.cxx
+++ b/reco/KF/Interface/CbmPVFinderKF.cxx
@@ -13,9 +13,9 @@
 
 #include <cmath>
 
-ClassImp(CbmPVFinderKF)
-
-  Int_t CbmPVFinderKF::FindPrimaryVertex(TClonesArray* tracks, CbmVertex* vertex)
+// ---------------------------------------------------------------------------------------------------------------------
+//
+Int_t CbmPVFinderKF::FindPrimaryVertex(TClonesArray* tracks, CbmVertex* vertex)
 {
 
   Int_t NTracks = tracks->GetEntriesFast();
@@ -48,8 +48,8 @@ ClassImp(CbmPVFinderKF)
   return vertex->GetNTracks();
 }
 
-
-// ----   Find vertex for one event   ---------------------------------------
+// ---------------------------------------------------------------------------------------------------------------------
+//
 Int_t CbmPVFinderKF::FindEventVertex(CbmEvent* event, TClonesArray* tracks)
 {
 
@@ -97,4 +97,5 @@ Int_t CbmPVFinderKF::FindEventVertex(CbmEvent* event, TClonesArray* tracks)
   delete[] trackArray;
   return vertex->GetNTracks();
 }
-// --------------------------------------------------------------------------
+
+ClassImp(CbmPVFinderKF);
diff --git a/reco/KF/Interface/CbmPVFinderKF.h b/reco/KF/Interface/CbmPVFinderKF.h
index ab9fea4e40..fc32b5db67 100644
--- a/reco/KF/Interface/CbmPVFinderKF.h
+++ b/reco/KF/Interface/CbmPVFinderKF.h
@@ -10,24 +10,23 @@
 
 #include "CbmPrimaryVertexFinder.h"
 
-/// \class CbmPVFinderKF
-/// \brief Implementation of the primary vertex finder using KF utility
+/// \class  CbmPVFinderKF
+/// \brief  Implementation of the primary vertex finder using KF utility
 ///
 class CbmPVFinderKF : public CbmPrimaryVertexFinder {
-
  public:
   /// \brief Track type for PV recnostruction
-  enum ESourceTrackType {
+  enum ESourceTrackType
+  {
     kStsTrack    = 0,
     kGlobalTrack = 1
   };
 
   /// \brief Default constructor
-  CbmPVFinderKF() {};
+  CbmPVFinderKF(){};
 
   /// \brief Destructior
-  ~CbmPVFinderKF() {};
-
+  ~CbmPVFinderKF(){};
 
   /// \brief Execution of PV finding.
   /// \param tracks   TClonesArray of CbmStsTracks
@@ -41,12 +40,6 @@ class CbmPVFinderKF : public CbmPrimaryVertexFinder {
   /// \param tracks   TClonesArray of CbmStsTracks
   virtual Int_t FindEventVertex(CbmEvent* event, TClonesArray* tracks);
 
-  /// \brief Sets source track type for the primary vertex reconstruction
-  void SetSourceTrackType(ESourceTrackType type) { fTrackType = type; }
-
-private:
-  ESourceTrackType fTrackType = kStsTrack;
-
   ClassDef(CbmPVFinderKF, 1);
 };
 
diff --git a/reco/KF/Interface/CbmPVFinderKFGlobal.cxx b/reco/KF/Interface/CbmPVFinderKFGlobal.cxx
new file mode 100644
index 0000000000..c707f04dc6
--- /dev/null
+++ b/reco/KF/Interface/CbmPVFinderKFGlobal.cxx
@@ -0,0 +1,123 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file    CbmPVFinderKFGlobal.h
+/// \brief   Primary vertex finder from the global tracks (implementation)
+/// \since   09.10.2023
+/// \authors Anna Senger, Sergei Zharko
+
+#include "CbmPVFinderKFGlobal.h"
+
+#include "CbmEvent.h"
+#include "CbmGlobalTrack.h"
+#include "CbmKFPrimaryVertexFinder.h"
+#include "CbmKFTrack.h"
+#include "CbmKFVertex.h"
+#include "TClonesArray.h"
+
+#include <cmath>
+#include <vector>
+
+ClassImp(CbmPVFinderKFGlobal);
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+Int_t CbmPVFinderKFGlobal::FindPrimaryVertex(TClonesArray* tracks, CbmVertex* vertex)
+{
+  Int_t nTracks = tracks->GetEntriesFast();
+
+  CbmKFPrimaryVertexFinder finder;
+  std::vector<CbmKFTrack> vKFTracks(nTracks);
+  for (int iT = 0; iT < nTracks; ++iT) {
+    auto* globalTrack = dynamic_cast<CbmGlobalTrack*>(tracks->At(iT));
+
+    if (globalTrack->GetStsTrackIndex() < 0) {
+      continue;
+    }
+    if (globalTrack->GetTrdTrackIndex() < 0) {
+      continue;
+    }
+    if (globalTrack->GetTofTrackIndex() < 0) {
+      continue;
+    }
+
+    if (globalTrack->GetChi2() < 0.) continue;
+    CbmKFTrack& kfTrack = vKFTracks[iT];
+    kfTrack.SetGlobalTrack(*globalTrack);
+    if (!isfinite(kfTrack.GetTrack()[0]) || !isfinite(kfTrack.GetCovMatrix()[0])) continue;
+    finder.AddTrack(&kfTrack);
+  }
+  CbmKFVertex kfVertex;
+  finder.Fit(kfVertex);
+  kfVertex.GetVertex(*vertex);
+
+  // Re-fit vertices of the global tracks
+  for (int iT = 0; iT < nTracks; ++iT) {
+    auto* globalTrack = dynamic_cast<CbmGlobalTrack*>(tracks->At(iT));
+    auto& kfTrack     = vKFTracks[iT];
+    kfTrack.Fit2Vertex(kfVertex);
+    FairTrackParam par;
+    kfTrack.GetTrackParam(par);
+    globalTrack->SetParamPrimaryVertex(&par);
+  }
+
+  return vertex->GetNTracks();
+}
+
+// ---------------------------------------------------------------------------------------------------------------------
+//
+Int_t CbmPVFinderKFGlobal::FindEventVertex(CbmEvent* event, TClonesArray* tracks)
+{
+  assert(event);
+  CbmKFPrimaryVertexFinder finder;
+
+  // Get vertex object
+  CbmVertex* vertex = event->GetVertex();
+
+  // Copy input tracks to KF tracks
+  Int_t nTracks = event->GetNofData(ECbmDataType::kStsTrack);
+  if (nTracks <= 0) return 0;
+
+  std::vector<CbmKFTrack> vKFTracks(nTracks);
+  for (Int_t iT = 0; iT < nTracks; ++iT) {
+    auto iTrkEvent    = event->GetIndex(ECbmDataType::kGlobalTrack, iT);
+    auto* globalTrack = dynamic_cast<CbmGlobalTrack*>(tracks->At(iTrkEvent));
+
+    if (globalTrack->GetStsTrackIndex() < 0) {
+      continue;
+    }
+    if (globalTrack->GetTrdTrackIndex() < 0) {
+      continue;
+    }
+    if (globalTrack->GetTofTrackIndex() < 0) {
+      continue;
+    }
+
+    if (globalTrack->GetChi2() < 0.) continue;
+    CbmKFTrack& kfTrack = vKFTracks[iTrkEvent];
+    kfTrack.SetGlobalTrack(*globalTrack);
+    if (!isfinite(kfTrack.GetTrack()[0]) || !isfinite(kfTrack.GetCovMatrix()[0])) continue;
+    finder.AddTrack(&kfTrack);
+  }
+
+  // Do the vertex finding
+  CbmKFVertex kfVertex;
+  finder.Fit(kfVertex);
+
+  // Copy KFVertex into CbmVertex
+  kfVertex.GetVertex(*vertex);
+
+  // Re-fit vertices of the global tracks
+  for (int iT = 0; iT < nTracks; ++iT) {
+    auto iTrkEvent    = event->GetIndex(ECbmDataType::kGlobalTrack, iT);
+    auto* globalTrack = dynamic_cast<CbmGlobalTrack*>(tracks->At(iTrkEvent));
+    auto& kfTrack     = vKFTracks[iTrkEvent];
+    kfTrack.Fit2Vertex(kfVertex);
+    FairTrackParam par;
+    kfTrack.GetTrackParam(par);
+    globalTrack->SetParamPrimaryVertex(&par);
+  }
+
+  return vertex->GetNTracks();
+}
diff --git a/reco/KF/Interface/CbmPVFinderKFGlobal.h b/reco/KF/Interface/CbmPVFinderKFGlobal.h
new file mode 100644
index 0000000000..685f506790
--- /dev/null
+++ b/reco/KF/Interface/CbmPVFinderKFGlobal.h
@@ -0,0 +1,40 @@
+/* Copyright (C) 2023 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergei Zharko [committer] */
+
+/// \file    CbmPVFinderKFGlobal.h
+/// \brief   Primary vertex finder from the global tracks (header)
+/// \since   09.10.2023
+/// \authors Anna Senger, Sergei Zharko
+
+#ifndef CbmPVFinderKFGlobal_h
+#define CbmPVFinderKFGlobal_h 1
+
+#include "CbmPrimaryVertexFinder.h"
+
+/// \class  CbmPVFinderKFGlobal
+/// \brief  Implementation of the primary vertex finder using KF utility
+///
+class CbmPVFinderKFGlobal : public CbmPrimaryVertexFinder {
+ public:
+  /// \brief Default constructor
+  CbmPVFinderKFGlobal() = default;
+
+  /// \brief Destructior
+  ~CbmPVFinderKFGlobal() = default;
+
+  /// \brief Execution of PV finding.
+  /// \param tracks   TClonesArray of CbmStsTracks
+  /// \param vertex   Primary vertex (output)
+  /// \param event    Pointer to event object
+  virtual Int_t FindPrimaryVertex(TClonesArray* tracks, CbmVertex* vertex);
+
+  /// \brief Execution of PV finding.
+  /// \param event    Pointer to event object
+  /// \param tracks   TClonesArray of CbmStsTracks
+  virtual Int_t FindEventVertex(CbmEvent* event, TClonesArray* tracks);
+
+  ClassDef(CbmPVFinderKFGlobal, 1);
+};
+
+#endif
diff --git a/reco/KF/KF.cmake b/reco/KF/KF.cmake
index 6e17f6ec7d..f84ecc56e2 100644
--- a/reco/KF/KF.cmake
+++ b/reco/KF/KF.cmake
@@ -30,6 +30,7 @@ set(SRCS
   Interface/CbmKFTrdHit.cxx 
   Interface/CbmKFTofHit.cxx 
   Interface/CbmPVFinderKF.cxx 
+  Interface/CbmPVFinderKFGlobal.cxx 
   Interface/CbmStsFitPerformanceTask.cxx 
   Interface/CbmStsKFTrackFitter.cxx 
   Interface/CbmStsKFSecondaryVertexFinder.cxx 
diff --git a/reco/KF/KFLinkDef.h b/reco/KF/KFLinkDef.h
index d2fe9ead38..5afbe7be46 100644
--- a/reco/KF/KFLinkDef.h
+++ b/reco/KF/KFLinkDef.h
@@ -31,6 +31,7 @@
 #pragma link C++ class CbmStsKFTrackFitter + ;
 #pragma link C++ class CbmStsKFSecondaryVertexFinder + ;
 #pragma link C++ class CbmPVFinderKF + ;
+#pragma link C++ class CbmPVFinderKFGlobal + ;
 #pragma link C++ class CbmStsFitPerformanceTask + ;
 //#pragma link C++ class  CbmEcalTrackExtrapolationKF+;
 #pragma link C++ class CbmTrdTrackFitterKF + ;
diff --git a/reco/global/CbmFindPrimaryVertex.cxx b/reco/global/CbmFindPrimaryVertex.cxx
index eac1c4eb52..d480793b8b 100644
--- a/reco/global/CbmFindPrimaryVertex.cxx
+++ b/reco/global/CbmFindPrimaryVertex.cxx
@@ -8,15 +8,15 @@
 // -------------------------------------------------------------------------
 #include "CbmFindPrimaryVertex.h"
 
+#include "CbmDefs.h"
 #include "CbmEvent.h"
 #include "CbmPrimaryVertexFinder.h"
 #include "CbmVertex.h"
-
 #include "FairRootManager.h"
-#include <Logger.h>
-
 #include "TClonesArray.h"
 
+#include <Logger.h>
+
 #include <iomanip>
 #include <iostream>
 
@@ -90,8 +90,8 @@ InitStatus CbmFindPrimaryVertex::Init()
     return kFATAL;
   }
 
-  // Get CbmStsTrack array
-  fTracks = (TClonesArray*) ioman->GetObject("StsTrack");
+  // Get CbmStsTrack or CbmGlobalTrack array
+  fTracks = (TClonesArray*) ioman->GetObject((ECbmDataType::kGlobalTrack == fTrackType) ? "GlobalTrack" : "StsTrack");
   if (!fTracks) {
     cout << "-W- CbmFindPrimaryVertex::Init: No STSTrack array!" << endl;
     return kERROR;
diff --git a/reco/global/CbmFindPrimaryVertex.h b/reco/global/CbmFindPrimaryVertex.h
index 4f3fef5a5e..b9e37498fd 100644
--- a/reco/global/CbmFindPrimaryVertex.h
+++ b/reco/global/CbmFindPrimaryVertex.h
@@ -24,10 +24,9 @@
 #define CBMFINDPRIMARYVERTEX_H 1
 
 
+#include "CbmDefs.h"
 #include "CbmVertex.h"
-
 #include "FairTask.h"
-
 #include "TStopwatch.h"
 
 class TClonesArray;
@@ -71,13 +70,19 @@ public:
   /** Finish **/
   virtual void Finish();
 
+  /** 
+   * @brief Sets type of tracks used for building the primary vertex 
+   * @note  the ECbmDataType::kStsTrack type use as default
+   **/
+  void SetTrackType(ECbmDataType trackType) { fTrackType = trackType; }
 
-private:
+ private:
   TStopwatch fTimer;
   CbmPrimaryVertexFinder* fFinder;
   TClonesArray* fEvents = nullptr;
   TClonesArray* fTracks;
   CbmVertex* fPrimVert;
+  ECbmDataType fTrackType = ECbmDataType::kStsTrack;
 
   Int_t fNofTs            = 0;   ///< Number of processed timeslices
   Int_t fNofEvents        = 0;   ///< Number of processed events
-- 
GitLab