From d1c53b7545b11624333d487e5c47fd26733d937c Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Sun, 6 Nov 2022 23:57:50 +0000
Subject: [PATCH] L1: debugger class

---
 reco/L1/CMakeLists.txt              |  1 +
 reco/L1/CbmL1.cxx                   |  1 +
 reco/L1/CbmL1.h                     |  2 +
 reco/L1/L1Algo/L1Algo.cxx           | 13 +++-
 reco/L1/L1Algo/L1Algo.h             | 10 ++-
 reco/L1/L1Algo/L1CATrackFinder.cxx  | 21 +++++-
 reco/L1/catools/CaToolsDebugger.cxx | 91 ++++++++++++++++++++++++++
 reco/L1/catools/CaToolsDebugger.h   | 99 +++++++++++++++++++++++++++++
 8 files changed, 232 insertions(+), 6 deletions(-)
 create mode 100644 reco/L1/catools/CaToolsDebugger.cxx
 create mode 100644 reco/L1/catools/CaToolsDebugger.h

diff --git a/reco/L1/CMakeLists.txt b/reco/L1/CMakeLists.txt
index 4d4deee982..8057a67b64 100644
--- a/reco/L1/CMakeLists.txt
+++ b/reco/L1/CMakeLists.txt
@@ -77,6 +77,7 @@ set(SRCS
   catools/CaToolsMcData.cxx
   catools/CaToolsMcDataManager.cxx
   catools/CaToolsPerformance.cxx
+  catools/CaToolsDebugger.cxx
   ParticleFinder/CbmL1PFFitter.cxx
   ParticleFinder/CbmL1PFMCParticle.cxx
 
diff --git a/reco/L1/CbmL1.cxx b/reco/L1/CbmL1.cxx
index 2c62edb76b..94dac75987 100644
--- a/reco/L1/CbmL1.cxx
+++ b/reco/L1/CbmL1.cxx
@@ -1235,6 +1235,7 @@ void CbmL1::Finish()
 
   gFile      = currentFile;
   gDirectory = curr;
+  fpAlgo->Finish();
 }
 
 
diff --git a/reco/L1/CbmL1.h b/reco/L1/CbmL1.h
index 64b201f2d7..bed97439cd 100644
--- a/reco/L1/CbmL1.h
+++ b/reco/L1/CbmL1.h
@@ -323,6 +323,8 @@ public:
 
   const L1Vector<CbmL1MCPoint>& GetMcPoints() const { return fvMCPoints; }
 
+  const L1Vector<CbmL1MCTrack>& GetMcTracks() const { return fvMCTracks; }
+
   const L1Vector<int>& GetHitMCRefs() const { return fvHitPointIndexes; }
 
   void SetUseMcHit(int MvdUseMcHit = 0, int StsUseMcHit = 0, int MuchUseMcHit = 0, int TrdUseMcHit = 0,
diff --git a/reco/L1/L1Algo/L1Algo.cxx b/reco/L1/L1Algo/L1Algo.cxx
index 55bb0822ca..71fc7a1a9b 100644
--- a/reco/L1/L1Algo/L1Algo.cxx
+++ b/reco/L1/L1Algo/L1Algo.cxx
@@ -6,6 +6,7 @@
 
 #include "CbmL1.h"
 
+#include "CaToolsDebugger.h"
 #include "L1Grid.h"
 #include "L1HitPoint.h"
 
@@ -77,6 +78,8 @@ void L1Algo::Init(const bool UseHitErrors, const TrackingMode mode, const bool M
   fMissingHits  = MissingHits;
 }
 
+void L1Algo::Finish() { ca::tools::Debugger::Instance().Write(); }
+
 // ---------------------------------------------------------------------------------------------------------------------
 //
 void L1Algo::ReceiveInputData(L1InputData&& inputData)
@@ -185,7 +188,7 @@ void L1Algo::CreateHitPoint(const L1Hit& hit, L1HitPoint& point)
   point.Set(hit.z, hit.u, hit.v, hit.du, hit.dv, hit.t, hit.dt);
 }
 
-int L1Algo::GetMcTrackIdForHit(int iHit)
+int L1Algo::GetMcTrackIdForHit(int iHit) const
 {
   int hitId    = iHit;
   int iMcPoint = CbmL1::Instance()->GetHitMCRefs()[hitId];
@@ -193,7 +196,7 @@ int L1Algo::GetMcTrackIdForHit(int iHit)
   return CbmL1::Instance()->GetMcPoints()[iMcPoint].ID;
 }
 
-int L1Algo::GetMcTrackIdForUnusedHit(int iHit)
+int L1Algo::GetMcTrackIdForUnusedHit(int iHit) const
 {
   int hitId    = (*RealIHitP)[iHit];
   int iMcPoint = CbmL1::Instance()->GetHitMCRefs()[hitId];
@@ -201,6 +204,12 @@ int L1Algo::GetMcTrackIdForUnusedHit(int iHit)
   return CbmL1::Instance()->GetMcPoints()[iMcPoint].ID;
 }
 
+const CbmL1MCTrack* L1Algo::GetMcTrackForUnusedHit(int iHit) const
+{
+  int id = GetMcTrackIdForUnusedHit(iHit);
+  if (id < 0) return nullptr;
+  return &CbmL1::Instance()->GetMcTracks()[id];
+}
 
 //   bool L1Algo::SortTrip(TripSort const& a, TripSort const& b) {
 //       return   ( a.trip.GetLevel() >  b.trip.GetLevel() );
diff --git a/reco/L1/L1Algo/L1Algo.h b/reco/L1/L1Algo/L1Algo.h
index f6a561ed1c..f30ce0aaed 100644
--- a/reco/L1/L1Algo/L1Algo.h
+++ b/reco/L1/L1Algo/L1Algo.h
@@ -46,6 +46,8 @@ class L1AlgoDraw;
 #include "L1Utils.h"
 #include "L1Vector.h"
 
+class CbmL1MCTrack;
+
 #ifdef _OPENMP
 #include "omp.h"
 #endif
@@ -366,6 +368,8 @@ public:
 
   void Init(const bool UseHitErrors, const TrackingMode mode, const bool MissingHits);
 
+  void Finish();
+
   void PrintHits();
 
   /// The main procedure - find tracks.
@@ -389,10 +393,12 @@ public:
   int GetNfieldStations() const { return fNfieldStations; }
 
   /// Get mc track ID for a hit (debug tool)
-  int GetMcTrackIdForHit(int iHit);
+  int GetMcTrackIdForHit(int iHit) const;
 
   /// Get mc track ID for a hit (debug tool)
-  int GetMcTrackIdForUnusedHit(int iHit);
+  int GetMcTrackIdForUnusedHit(int iHit) const;
+
+  const CbmL1MCTrack* GetMcTrackForUnusedHit(int iHit) const;
 
 private:
   bool checkTripletMatch(const L1Triplet& l, const L1Triplet& r, fscal& dchi2) const;
diff --git a/reco/L1/L1Algo/L1CATrackFinder.cxx b/reco/L1/L1Algo/L1CATrackFinder.cxx
index d885033dc8..e65c4cbb15 100644
--- a/reco/L1/L1Algo/L1CATrackFinder.cxx
+++ b/reco/L1/L1Algo/L1CATrackFinder.cxx
@@ -19,6 +19,9 @@
  */
 
 
+#include "CbmL1.h"
+#include "CbmL1MCTrack.h"
+
 #include "L1Algo.h"
 #include "L1Assert.h"
 #include "L1Branch.h"
@@ -63,6 +66,8 @@
 
 #include <stdio.h>
 
+#include "CaToolsDebugger.h"
+
 
 // using namespace std;
 using std::cout;
@@ -964,8 +969,7 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
 
   /// Refit Triplets
 
-  L1Fit fit;
-  fit.SetParticleMass(GetDefaultParticleMass());
+  ca::tools::Debugger::Instance().AddNtuple("triplets", "ev:mc:sta:p:vx:vy:vz:chi2");
 
   L1TrackParFit fitNew;
   fitNew.SetParticleMass(GetDefaultParticleMass());
@@ -1121,7 +1125,20 @@ inline void L1Algo::findTripletsStep2(Tindex n3, int istal, int istam, int istar
     //cout << " chi2 " << T3.chi2[i3_4] << " " << T.chi2[0] << endl;
     //T.chi2 = (fscal) T3.chi2[i3_4];
     T3.SetOneEntry(i3_4, T, i3_4);
+
+    {
+      int mc1 = GetMcTrackIdForUnusedHit(ihit[0]);
+      int mc2 = GetMcTrackIdForUnusedHit(ihit[1]);
+      int mc3 = GetMcTrackIdForUnusedHit(ihit[2]);
+      if ((mc1 >= 0) && (mc1 == mc2) && (mc1 == mc3)) {
+        const CbmL1MCTrack& mctr = CbmL1::Instance()->GetMcTracks()[mc1];
+        float ev                 = 0;
+        float chi2               = T.chi2[i3_4];
+        ca::tools::Debugger::Instance().FillNtuple("triplets", ev, mc1, istal, mctr.p, mctr.x, mctr.y, mctr.z, chi2);
+      }
+    }
   }  //i3
+
 }  // findTripletsStep2
 
 
diff --git a/reco/L1/catools/CaToolsDebugger.cxx b/reco/L1/catools/CaToolsDebugger.cxx
new file mode 100644
index 0000000000..1346296d4a
--- /dev/null
+++ b/reco/L1/catools/CaToolsDebugger.cxx
@@ -0,0 +1,91 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov, Sergei Zharko [committer] */
+
+/// \file   CaToolsDebugger.h
+/// \brief  Tracking Debugger class (implementation)
+/// \since  23.09.2022
+/// \author S.Zharko <s.zharko@gsi.de>
+
+/*
+root -l CAdebug.root
+TNtuple *n = (TNtuple*) gDirectory->FindObjectAny("triplets")
+n->Draw("chi2")
+*/
+
+#include "CaToolsDebugger.h"
+
+#include "TFile.h"
+#include "TNtuple.h"
+
+#include <iostream>
+#include <memory>
+#include <string>
+
+#include <stdarg.h>
+
+using namespace ca::tools;
+
+Debugger& Debugger::Instance()
+{
+  static std::unique_ptr<Debugger> instance;
+  if (!instance) { instance.reset(new Debugger); }
+  return *instance.get();
+}
+
+
+Debugger::Debugger(const char* fileName)
+{
+  /// Default constructor
+  fFile = new TFile(fileName, "RECREATE");
+}
+
+void Debugger::Write()
+{
+  /// Write data to file
+  fIsWritten = true;
+  auto currDir {gDirectory};
+  if (fFile) {
+    fFile->cd();
+    for (unsigned int i = 0; i < fNtuples.size(); i++) {
+      fNtuples[i]->Write();
+    }
+  }
+  gDirectory = currDir;
+}
+
+Debugger::~Debugger()
+{
+  /// Destructor
+  if (!fIsWritten) { std::cerr << "CaDebugger: you forgot to Write()!!" << std::endl; }
+}
+
+void Debugger::AddNtuple(const char* name, const char* varlist)
+{
+  /// add ntuple
+  assert(fFile);
+  if (GetNtupleIndex(name) >= 0) return;
+  auto currDir {gDirectory};
+  fFile->cd();
+  fNtuples.push_back(new TNtuple(name, name, varlist));
+  gDirectory = currDir;
+}
+
+int Debugger::GetNtupleIndex(const char* name)
+{
+  int ind = -1;
+  for (unsigned int i = 0; i < fNtuples.size(); i++) {
+    if (strcmp(name, fNtuples[i]->GetName()) == 0) {
+      ind = i;
+      break;
+    }
+  }
+  return ind;
+}
+
+void Debugger::FillNtuple(const char* name, float v[])
+{
+  int iNtuple = GetNtupleIndex(name);
+  if (iNtuple < 0) { std::cerr << "CaDebugger: Ntuple (" << name << ") doesn't exist" << std::endl; }
+  fNtuples[iNtuple]->Fill(v);
+}
\ No newline at end of file
diff --git a/reco/L1/catools/CaToolsDebugger.h b/reco/L1/catools/CaToolsDebugger.h
new file mode 100644
index 0000000000..8c4aa8fffa
--- /dev/null
+++ b/reco/L1/catools/CaToolsDebugger.h
@@ -0,0 +1,99 @@
+/* Copyright (C) 2022 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
+   SPDX-License-Identifier: GPL-3.0-only
+   Authors: Sergey Gorbunov[committer] */
+
+/// \file   CaToolsDebugger.h
+/// \brief  Tracking Debugger class (header)
+/// \since  03.10.2022
+/// \author S.Gorbunov
+
+#ifndef CaToolsDebugger_h
+#define CaToolsDebugger_h 1
+
+#include <iostream>
+#include <vector>
+
+class TFile;
+class TNtuple;
+
+namespace ca
+{
+  namespace tools
+  {
+    /// Class ca::tools::Debugger helps to debug CA tracking
+    ///
+    class Debugger {
+    public:
+      // *****************************************
+      // **     Constructors and destructor     **
+      // *****************************************
+
+
+      /// Default constructor
+      Debugger(const char* fileName = "CAdebug.root");
+
+      /// Destructor
+      ~Debugger();
+
+      /// Copy constructor
+      Debugger(const Debugger& other) = delete;
+
+      /// Move constructor
+      Debugger(Debugger&& other) = delete;
+
+      /// Copy assignment operator
+      Debugger& operator=(const Debugger& other) = delete;
+
+      /// Move assignment operator
+      Debugger& operator=(Debugger&& other) = delete;
+
+      /// Instance
+      static Debugger& Instance();
+
+      /// Write ntuples to the file
+      void Write();
+
+      /// Set new ntuple
+      void AddNtuple(const char* name, const char* varlist);
+
+      /// Add an entry to ntuple
+      void FillNtuple(const char* name, float v[]);
+
+      /// Add an entry to ntuple
+      template<typename... Targs>
+      void FillNtuple(const char* name, Targs... args)
+      {
+        constexpr std::size_t n = sizeof...(Targs);
+        if (n <= 0) return;
+        float v[n];
+        FillFloatArray(v, args...);
+        FillNtuple(name, v);
+      }
+
+      /// Get ntuple index
+      int GetNtupleIndex(const char* name);
+
+    private:
+      template<typename T, typename... Targs>
+      void FillFloatArray(float* v, T val, Targs... args)
+      {
+        v[0] = (float) val;
+        if (sizeof...(args) > 0) { FillFloatArray(v + 1, args...); }
+      }
+
+      template<typename T>
+      void FillFloatArray(float* v, T last)
+      {
+        v[0] = (float) last;
+      }
+
+    private:
+      bool fIsWritten {0};
+      TFile* fFile {nullptr};
+      std::vector<TNtuple*> fNtuples;
+    };
+
+  }  // namespace tools
+}  // namespace ca
+
+#endif  // CaToolsDebugger_h
-- 
GitLab