From be184d4b1ad9659c93a6e7fa3919686b3bff33e8 Mon Sep 17 00:00:00 2001
From: "se.gorbunov" <se.gorbunov@gsi.de>
Date: Thu, 21 Sep 2023 16:20:34 +0000
Subject: [PATCH] CA: tune cuts for mcbm tracker

---
 macro/L1/configs/ca_params_mcbm.yaml    |  2 +-
 reco/L1/CbmL1ReadEvent.cxx              | 16 ++++++++++++++
 reco/L1/L1Algo/L1CaTrackFinderSlice.cxx | 22 ++++++++++++-------
 reco/L1/L1Algo/L1TripletConstructor.cxx | 29 ++++++++++++++++---------
 4 files changed, 50 insertions(+), 19 deletions(-)

diff --git a/macro/L1/configs/ca_params_mcbm.yaml b/macro/L1/configs/ca_params_mcbm.yaml
index 260f7219e8..4cb3acbb00 100644
--- a/macro/L1/configs/ca_params_mcbm.yaml
+++ b/macro/L1/configs/ca_params_mcbm.yaml
@@ -70,7 +70,7 @@ ca:
       triplet_final_chi2_cut:   23.4450
       doublet_chi2_cut:         7.56327
       pick_gather:              4.0
-      triplet_link_chi2:        25.
+      triplet_link_chi2:        125.
       max_qp:                   2.   # 1 / 0.5 [GeV/c]^-1
       max_slope_pv:             9.5
       max_slope:                2.748
diff --git a/reco/L1/CbmL1ReadEvent.cxx b/reco/L1/CbmL1ReadEvent.cxx
index 36f6276a45..e6cb66f771 100644
--- a/reco/L1/CbmL1ReadEvent.cxx
+++ b/reco/L1/CbmL1ReadEvent.cxx
@@ -48,6 +48,7 @@
 
 #include <iostream>
 
+#include "CaToolsDebugger.h"
 #include "L1Algo/L1Algo.h"
 
 using std::cout;
@@ -778,6 +779,21 @@ void CbmL1::ReadEvent(CbmEvent* event)
   fpIODataManager->ResetInputData(nHits);
   fpIODataManager->SetNhitKeys(NStrips);
 
+  if (0) {
+    ca::tools::Debugger::Instance().AddNtuple("hits", "ev:mcpoint:mc:mcMother:nmc:sta:x:y:z");
+    for (int iHit = 0; iHit < nHits; ++iHit) {
+      TmpHit& h    = tmpHits[iHit];
+      int mc       = -1;
+      int mcMother = -1;
+      if (h.iBestMc >= 0) {
+        mc       = fvMCPoints[h.iBestMc].ID;
+        mcMother = fvMCPoints[h.iBestMc].mother_ID;
+      }
+      ca::tools::Debugger::Instance().FillNtuple("hits", nCalls, h.iBestMc, mc, mcMother, h.vMc.size(), h.iStation, h.x,
+                                                 h.y, h.z);
+    }
+  }
+
   // ----- Fill
   for (int iHit = 0; iHit < nHits; ++iHit) {
     TmpHit& th = tmpHits[iHit];
diff --git a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
index 3b5e1574f2..504e7d8f65 100644
--- a/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
+++ b/reco/L1/L1Algo/L1CaTrackFinderSlice.cxx
@@ -87,7 +87,6 @@ bool L1Algo::checkTripletMatch(const L1Triplet& l, const L1Triplet& r, fscal& dc
     dchi2 = dqp * dqp / Cqp;
   }
   else {
-
     fscal dtx = l.GetTx() - r.GetTx();
     fscal Ctx = l.GetCtx() + r.GetCtx();
 
@@ -104,7 +103,8 @@ bool L1Algo::checkTripletMatch(const L1Triplet& l, const L1Triplet& r, fscal& dc
     if (dty * dty > fTripletLinkChi2 * Cty) return false;
     if (dtx * dtx > fTripletLinkChi2 * Ctx) return false;
 
-    dchi2 = 0.5f * (dtx * dtx / Ctx + dty * dty / Cty);
+    //dchi2 = 0.5f * (dtx * dtx / Ctx + dty * dty / Cty);
+    dchi2 = 0.;
   }
 
   if (!std::isfinite(dchi2)) return false;
@@ -514,8 +514,8 @@ void L1Algo::CaTrackFinderSlice()
 
           CAFindTrack(istaF, best_tr, best_L, best_chi2, &first_trip, (curr_tr), curr_L, curr_chi2, min_best_l,
                       new_tr);  /// reqursive func to build a tree of possible track-candidates and choose the best
-          //              if ( best_L < min_best_l ) continue;
-          if (best_L < firstTripletLevel + 2) continue;  // lose maximum one hit
+
+          if (best_L < firstTripletLevel + 2) continue;  // loose maximum one hit
 
           if (best_L < min_level + 3) continue;  // should find all hits for min_level
 
@@ -841,7 +841,7 @@ void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best_L, fsc
 
     if (kGlobal == fTrackingMode || kMcbm == fTrackingMode) { ndf = curr_L * 2 - 4; }
 
-    if (curr_chi2 > fTrackChi2Cut * ndf) return;
+    if (curr_chi2 > fTrackChi2Cut * ndf) { return; }
 
     // TODO: SZh 04.10.2022: Does one should delete this code?
     //       // try to find more hits
@@ -913,14 +913,20 @@ void L1Algo::CAFindTrack(int ista, L1Branch& best_tr, unsigned char& best_L, fsc
           }
         }
 
+        int ndf = 2 * (new_L + 2) - 5;
+
         if (kGlobal == fTrackingMode) {  //SGtrd2d!!!
-          // nHits in chi2 calculation == new_L+2, ndf == 4
-          if (new_chi2 > fTrackChi2Cut * (2 * (new_L + 2) - 4)) continue;
+          ndf = 2 * (new_L + 2) - 4;
+        }
+        else if (kMcbm == fTrackingMode) {
+          ndf = 2 * (new_L + 2) - 4;
         }
         else {
-          if (new_chi2 > fTrackChi2Cut * new_L) continue;  // TODO: SG: it must be  ( 2 * new_L )
+          ndf = new_L;  // TODO: SG:  2 * (new_L + 2) - 5
         }
 
+        if (new_chi2 > fTrackChi2Cut * ndf) continue;
+
         // add new hit
         new_tr[ista] = curr_tr;
         new_tr[ista].fHits.push_back(fGridHitIds[new_trip.GetLHit()]);
diff --git a/reco/L1/L1Algo/L1TripletConstructor.cxx b/reco/L1/L1Algo/L1TripletConstructor.cxx
index ad90a88bf5..d07a1552a3 100644
--- a/reco/L1/L1Algo/L1TripletConstructor.cxx
+++ b/reco/L1/L1Algo/L1TripletConstructor.cxx
@@ -375,7 +375,10 @@ void L1TripletConstructor::FitTriplets()
   fTracks_3.reset(n3, L1TrackPar());
 
   /// Refit Triplets
-  if (dumpTriplets) { ca::tools::Debugger::Instance().AddNtuple("triplets", "ev:mc:sta:p:vx:vy:vz:chi2:ndf"); }
+  if (dumpTriplets) {
+    ca::tools::Debugger::Instance().AddNtuple("triplets",
+                                              "ev:iter:i0:x0:y0:z0:i1:x1:y1:z1:i2:x2:y2:z2:mc:sta:p:vx:vy:vz:chi2:ndf");
+  }
 
   L1Fit fit;
   fit.SetMask(fmask::One());
@@ -531,21 +534,27 @@ void L1TripletConstructor::FitTriplets()
     fTracks_3[i3] = T;
 
     if (dumpTriplets) {
-      int mc1 = fAlgo->GetMcTrackIdForCaHit(ihit[0]);
-      int mc2 = fAlgo->GetMcTrackIdForCaHit(ihit[1]);
-      int mc3 = fAlgo->GetMcTrackIdForCaHit(ihit[2]);
+      int ih0 = ihit[0];
+      int ih1 = ihit[1];
+      int ih2 = ihit[2];
+      int mc1 = fAlgo->GetMcTrackIdForCaHit(ih0);
+      int mc2 = fAlgo->GetMcTrackIdForCaHit(ih1);
+      int mc3 = fAlgo->GetMcTrackIdForCaHit(ih2);
+
+      const L1Hit& h0 = fAlgo->fInputData.GetHit(ih0);
+      const L1Hit& h1 = fAlgo->fInputData.GetHit(ih1);
+      const L1Hit& h2 = fAlgo->fInputData.GetHit(ih2);
+
       if ((mc1 >= 0) && (mc1 == mc2) && (mc1 == mc3)) {
         const CbmL1MCTrack& mctr = CbmL1::Instance()->GetMcTracks()[mc1];
-        float ev                 = 0;
-        float chi2               = T.chi2[0];
-        float nd                 = T.NDF[0];
-        ca::tools::Debugger::Instance().FillNtuple("triplets", ev, mc1, fIstaL, mctr.p, mctr.x, mctr.y, mctr.z, chi2,
-                                                   nd);
+        ca::tools::Debugger::Instance().FillNtuple("triplets", mctr.iEvent, fAlgo->isec, ih0, h0.x, h0.y, h0.z, ih1,
+                                                   h1.x, h1.y, h1.z, ih2, h2.x, h2.y, h2.z, mc1, fIstaL, mctr.p, mctr.x,
+                                                   mctr.y, mctr.z, (float) T.chi2[0], (float) T.NDF[0]);
       }
     }
   }  //i3
 
-}  // findTripletsStep2
+}  // FitTriplets
 
 
 void L1TripletConstructor::StoreTriplets()
-- 
GitLab