From 9e4939f9a327840ee85721db73fbede7e665593a Mon Sep 17 00:00:00 2001
From: Nora Bluhme <bluhme@fias.uni-frankfurt.de>
Date: Fri, 31 Jan 2025 21:58:39 +0100
Subject: [PATCH] Reduced BBA output and changed log level for chi square
 mismatch message in CbmKfTrackFitter.

---
 reco/alignment/CbmBbaAlignmentTask.cxx | 79 ++++++++++++++++++++++----
 reco/kfnew/CbmKfTrackFitter.cxx        |  2 +-
 2 files changed, 70 insertions(+), 11 deletions(-)

diff --git a/reco/alignment/CbmBbaAlignmentTask.cxx b/reco/alignment/CbmBbaAlignmentTask.cxx
index 22a7a29f31..1f4aa751ba 100644
--- a/reco/alignment/CbmBbaAlignmentTask.cxx
+++ b/reco/alignment/CbmBbaAlignmentTask.cxx
@@ -195,6 +195,8 @@ InitStatus CbmBbaAlignmentTask::Init()
   ca::Framework* l1                       = CbmL1::Instance()->fpAlgo;
   const ca::Parameters<ca::fvec>& l1Param = l1->GetParameters();
   fNtrackingStations                      = l1Param.GetNstationsActive();
+  // TO get rid of weird station counting:
+  // fNtrackingStations                      = l1Param.GetNstationsActive() - l1Param. GetNstationsActive(ca::EDetectorID::kMvd);
 
   TDirectory* curDirectory = gDirectory;
 
@@ -301,6 +303,7 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
   fFitter.SetDoSmooth(true);
 
   // select tracks for alignment and store them
+  unsigned num_rejected_fit = 0;
   if (fTrackingMode == kSts && fInputStsTracks) {
 
     fFitter.SetDefaultMomentumForMs(0.1);
@@ -317,10 +320,15 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
 
       for (auto& n : t.fUnalignedTrack.fNodes) {
         n.fIsTimeSet = false;
+        // attempt to switch off MVD
+        // if (n.fHitSystemId == ECbmModuleId::kMvd) {
+        //   n.fIsXySet = false;
+        // }
       }
 
       if (!fFitter.FitTrajectory(t.fUnalignedTrack)) {
-        LOG(warning) << "failed to fit the track! ";
+        LOG(debug) << "failed to fit the track! ";
+        num_rejected_fit++;
         continue;
       }
 
@@ -371,7 +379,8 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
       if (t.fNtofHits < 2) continue;
       if (t.fNtrd1dHits + t.fNtrd2dHits < 2) continue;
       if (!fFitter.FitTrajectory(t.fUnalignedTrack)) {
-        LOG(warning) << "failed to fit the track! ";
+        LOG(debug) << "failed to fit the track! ";
+        num_rejected_fit++;
         continue;
       }
       t.fAlignedTrack = t.fUnalignedTrack;
@@ -379,31 +388,50 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
     }
   }  // end of mcbm tracking mode
 
+  if (num_rejected_fit != 0) {
+    LOG(warn) << "failed to fit " << num_rejected_fit << " tracks";
+  }
 
   // fix the material radiation thickness to avoid the chi2 rattling at the material map bin boundaries
   // (to be improved)
+  unsigned num_fix_rad_max   = 0;
+  unsigned num_fix_rad_min   = 0;
+  unsigned num_rejected_fit2 = 0;
   for (TrackContainer& t : fTracks) {
     if (!t.fIsActive) continue;
     for (unsigned int in = 0; in < t.fUnalignedTrack.fNodes.size(); in++) {
       CbmKfTrackFitter::TrajectoryNode& node = t.fUnalignedTrack.fNodes[in];
       if (!node.fIsFitted) {
         t.fIsActive = false;
+        num_rejected_fit2++;
         break;
       }
       node.fIsRadThickFixed = true;
       if (node.fRadThick > 0.01) {
-        LOG(warning) << "CbmBbaAlignmentTask: radiation thickness is too large: " << node.fRadThick;
+        LOG(debug) << "CbmBbaAlignmentTask: radiation thickness is too large: " << node.fRadThick;
+        num_fix_rad_max++;
         node.fRadThick = 0.01;
       }
       if (node.fRadThick < 0.001) {
-        LOG(warning) << "CbmBbaAlignmentTask: radiation thickness is too small: " << node.fRadThick;
+        LOG(debug) << "CbmBbaAlignmentTask: radiation thickness is too small: " << node.fRadThick;
+        num_fix_rad_min++;
         node.fRadThick = 0.001;
       }
     }
   }
+  if (num_fix_rad_max != 0) {
+    LOG(warn) << "CbmBbaAlignmentTask: radiation thickness is too large " << num_fix_rad_max << " times";
+  }
+  if (num_fix_rad_min != 0) {
+    LOG(warn) << "CbmBbaAlignmentTask: radiation thickness is too small " << num_fix_rad_min << " times";
+  }
+  if (num_rejected_fit2 != 0) {
+    LOG(warn) << "Rejected fit 2  " << num_rejected_fit2 << " tracks";
+  }
 
   // ensure that all the hits are not too much deviated from the trajectory
   // (to be improved)
+  unsigned num_rejected_traj = 0;
   for (TrackContainer& t : fTracks) {
     if (!t.fIsActive) continue;
     const double cutX = 8.;
@@ -417,6 +445,7 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
       node.fIsXySet = false;  // exclude the node from the fit
       if (!fFitter.FitTrajectory(tr)) {
         t.fIsActive = false;
+        num_rejected_traj++;
         break;
       }
       double x_residual = node.fMxy.X() - node.fParamDn.X();  // this should be the difference vector
@@ -426,10 +455,14 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
       if (!node.fIsFitted || (node.fMxy.NdfX() > 0. && fabs(x_pull) > cutX)
           || (node.fMxy.NdfY() > 0. && fabs(y_pull) > cutY)) {
         t.fIsActive = false;
+        num_rejected_traj++;
         break;
       }
     }
   }
+  if (num_rejected_traj != 0) {
+    LOG(warn) << "Rejected  " << num_rejected_traj << " tracks for hits deviating from the trajectory";
+  }
 
   TClonesArray* inputTracks = nullptr;
   if (fTrackingMode == kSts) {
@@ -441,8 +474,13 @@ void CbmBbaAlignmentTask::Exec(Option_t* /*opt*/)
   static int statTracks = 0;
   statTracks += inputTracks->GetEntriesFast();
 
+  unsigned active_tracks = 0;
+  for (auto& t : fTracks) {
+    if (t.fIsActive) active_tracks++;
+  }
+
   LOG(info) << "Bba: " << inputTracks->GetEntriesFast() << " tracks in event, " << statTracks << " tracks in total, "
-            << fTracks.size() << " tracks selected for alignment";
+            << fTracks.size() << " tracks selected for alignment, " << active_tracks << " tracks active";
 }
 
 
@@ -614,15 +652,16 @@ double CbmBbaAlignmentTask::CostFunction(const std::vector<double>& par)
 
 void CbmBbaAlignmentTask::Finish()
 {
-  //
   // perform the alignment
-  //
 
   LOG(info) << "BBA: start the alignment procedure with " << fTracks.size() << " tracks ...";
 
   // collect the sensors from the hits
   // TODO: read it from the detector setup
 
+  if (fTrackingMode == kSts) {
+    fNtrackingStations = CbmStsTrackingInterface::Instance()->GetNtrackingStations();
+  }
 
   std::set<Sensor> sensorSet;
   for (auto& t : fTracks) {
@@ -635,6 +674,11 @@ void CbmBbaAlignmentTask::Finish()
       s.fSensorId = n.fHitAddress;
       // TODO: get the station index from n.fHitSystemId, n.fHitAddress
       s.fTrackingStation = n.fMaterialLayer;
+
+      if (fTrackingMode == kSts && s.fSystemId == ECbmModuleId::kSts) {
+        s.fTrackingStation = CbmStsTrackingInterface::Instance()->GetTrackingStationIndex(n.fHitAddress);
+      }
+
       if (s.fSystemId == ECbmModuleId::kTrd || s.fSystemId == ECbmModuleId::kTrd2d) {
         s.fSensorId = s.fTrackingStation;
         // hardcode path to physical node in geometry
@@ -651,6 +695,7 @@ void CbmBbaAlignmentTask::Finish()
       else if (s.fSystemId == ECbmModuleId::kTof) {
         s.fSensorId = CbmTofAddress::GetRpcFullId(n.fHitAddress);
       }
+
       sensorSet.insert(s);
     }
   }
@@ -670,6 +715,11 @@ void CbmBbaAlignmentTask::Finish()
       s.fSystemId        = n.fHitSystemId;
       s.fSensorId        = n.fHitAddress;
       s.fTrackingStation = n.fMaterialLayer;
+
+      if (fTrackingMode == kSts && s.fSystemId == ECbmModuleId::kSts) {
+        s.fTrackingStation = CbmStsTrackingInterface::Instance()->GetTrackingStationIndex(n.fHitAddress);
+      }
+
       if (s.fSystemId == ECbmModuleId::kTrd || s.fSystemId == ECbmModuleId::kTrd2d) {
         s.fSensorId = s.fTrackingStation;
       }
@@ -722,7 +772,7 @@ void CbmBbaAlignmentTask::Finish()
   {
     for (int is = 0; is < fNalignmentBodies; is++) {
       auto& b = fAlignmentBodies[is];
-      LOG(info) << "Body " << is << " sta " << b.fTrackingStation;
+      LOG(info) << "Body: " << is << " station: " << b.fTrackingStation;
     }
   }
 
@@ -844,14 +894,17 @@ void CbmBbaAlignmentTask::Finish()
 
     fCostInitial = CostFunction(parMisaligned);
 
+    unsigned tracks_rejected = 0;
     for (auto& t : fTracks) {
       double ndf  = t.fAlignedTrack.fNodes.back().fParamDn.GetNdf();
       double chi2 = t.fAlignedTrack.fNodes.back().fParamDn.GetChiSq();
       if (ndf < 0. || chi2 < 0. || !std::isfinite(ndf) || !std::isfinite(chi2)) {
         t.fIsActive = false;
+        tracks_rejected++;
       }
     }
-    std::cout << "reject bad tracks and recalculate the initial cost function..." << std::endl;
+    std::cout << "reject " << tracks_rejected << " bad tracks and recalculate the initial cost function..."
+              << std::endl;
     fCostInitial = CostFunction(parMisaligned);
     fFixedNdf    = -1;  //fNdfTotal;
 
@@ -1008,10 +1061,15 @@ void CbmBbaAlignmentTask::Finish()
     }
   }
 
+  unsigned active_tracks = 0;
+  for (auto& t : fTracks) {
+    if (t.fIsActive) active_tracks++;
+  }
 
   LOG(info) << "\n";
 
   LOG(info) << "Bba: " << fTracks.size() << " tracks used for alignment\n";
+  LOG(info) << "     " << active_tracks << " tracks active\n";
 
   LOG(info) << " Cost function for the true parameters: " << fCostIdeal;
   LOG(info) << " Cost function for the initial parameters: " << fCostInitial;
@@ -1019,7 +1077,8 @@ void CbmBbaAlignmentTask::Finish()
 
   for (auto& t : fTracks) {
     if (!t.fIsActive) continue;
-    // calculate the residuals for all tracks at each TrajectoryNode before alignment
+    // calculate the residuals for all tracks at each TrajectoryNode after alignment
+    // TODO: put into function
     for (unsigned int in = 0; in < t.fAlignedTrack.fNodes.size(); in++) {
       CbmKfTrackFitter::Trajectory tr        = t.fAlignedTrack;
       CbmKfTrackFitter::TrajectoryNode& node = tr.fNodes[in];
diff --git a/reco/kfnew/CbmKfTrackFitter.cxx b/reco/kfnew/CbmKfTrackFitter.cxx
index 4905e7b2d4..25674692b7 100644
--- a/reco/kfnew/CbmKfTrackFitter.cxx
+++ b/reco/kfnew/CbmKfTrackFitter.cxx
@@ -781,7 +781,7 @@ bool CbmKfTrackFitter::FitTrajectory(CbmKfTrackFitter::Trajectory& t)
     double upChi2 = fFit.Tr().GetChiSq();
     if (fabs(upChi2 - dnChi2) > 1.e-1) {
       //if (upChi2 > dnChi2 + 1.e-2) {
-      LOG(error) << "CbmKfTrackFitter: " << fDebugInfo << ": chi2 mismatch: dn " << dnChi2 << " != up " << upChi2
+      LOG(debug) << "CbmKfTrackFitter: " << fDebugInfo << ": chi2 mismatch: dn " << dnChi2 << " != up " << upChi2
                  << " first node " << firstHitNode << " last node " << lastHitNode << " of " << nNodes;
       //if (fVerbosityLevel > 0) {
       //exit(1);
-- 
GitLab