From d0197b5f8c51b7c844b9b2a038cdd64d8645f4e3 Mon Sep 17 00:00:00 2001
From: Axel Puntke <axel.puntke@uni-muenster.de>
Date: Mon, 28 Aug 2023 12:01:26 +0000
Subject: [PATCH] TRD Incomplete cluster recognition

---
 reco/detectors/trd/CbmTrdModuleRecR.cxx | 54 ++++++++++++++++++++++++-
 reco/detectors/trd/CbmTrdModuleRecR.h   |  1 +
 2 files changed, 54 insertions(+), 1 deletion(-)

diff --git a/reco/detectors/trd/CbmTrdModuleRecR.cxx b/reco/detectors/trd/CbmTrdModuleRecR.cxx
index b99feaec98..875891595a 100644
--- a/reco/detectors/trd/CbmTrdModuleRecR.cxx
+++ b/reco/detectors/trd/CbmTrdModuleRecR.cxx
@@ -440,7 +440,7 @@ void CbmTrdModuleRecR::addClusters(std::deque<std::pair<Int_t, const CbmTrdDigi*
 Bool_t CbmTrdModuleRecR::MakeHits() { return kTRUE; }
 
 //_______________________________________________________________________________
-CbmTrdHit* CbmTrdModuleRecR::MakeHit(Int_t clusterId, const CbmTrdCluster* /*cluster*/,
+CbmTrdHit* CbmTrdModuleRecR::MakeHit(Int_t clusterId, const CbmTrdCluster* cluster,
                                      std::vector<const CbmTrdDigi*>* digis)
 {
 
@@ -541,6 +541,9 @@ CbmTrdHit* CbmTrdModuleRecR::MakeHit(Int_t clusterId, const CbmTrdCluster* /*clu
     cluster_pad_dposV[1] = sqrt(fDigiPar->GetPadSizeY(1));
   }
 
+  // Set charge of incomplete clusters (missing NTs) to -1 (not deleting them because they are still relevant for tracking)
+  if (!IsClusterComplete(cluster)) totalCharge = -1.0;
+
   Int_t nofHits = fHits->GetEntriesFast();
 
   //  return new ((*fHits)[nofHits]) CbmTrdHit(fModAddress, global,
@@ -579,4 +582,53 @@ Double_t CbmTrdModuleRecR::GetSpaceResolution(Double_t val)
   return selval;
 }
 
+bool CbmTrdModuleRecR::IsClusterComplete(const CbmTrdCluster* cluster)
+{
+  uint32_t colMin = fDigiPar->GetNofColumns();
+  uint32_t rowMin = fDigiPar->GetNofRows();
+
+  for (int i = 0; i < cluster->GetNofDigis(); ++i) {
+    const CbmTrdDigi* digi = CbmDigiManager::Instance()->Get<CbmTrdDigi>(cluster->GetDigi(i));
+    int digiCol            = fDigiPar->GetPadColumn(digi->GetAddressChannel());
+    int digiRow            = fDigiPar->GetPadRow(digi->GetAddressChannel());
+
+    if (digiCol < colMin) colMin = digiCol;
+    if (digiRow < rowMin) rowMin = digiRow;
+  }
+
+  const UShort_t nCols = cluster->GetNCols();
+  const UShort_t nRows = cluster->GetNRows();
+
+  CbmTrdDigi* digiMap[nRows][nCols];                        //create array on stack for optimal performance
+  memset(digiMap, 0, sizeof(CbmTrdDigi*) * nCols * nRows);  //init with nullpointers
+
+  for (int i = 0; i < cluster->GetNofDigis(); ++i) {
+    const CbmTrdDigi* digi = CbmDigiManager::Instance()->Get<CbmTrdDigi>(cluster->GetDigi(i));
+    int digiCol            = fDigiPar->GetPadColumn(digi->GetAddressChannel());
+    int digiRow            = fDigiPar->GetPadRow(digi->GetAddressChannel());
+
+    if (digiMap[digiRow - rowMin][digiCol - colMin])
+      return false;  // To be investigated why this sometimes happens (Redmin Issue 2914)
+
+    digiMap[digiRow - rowMin][digiCol - colMin] = const_cast<CbmTrdDigi*>(digi);
+  }
+
+  // check if each row of the cluster starts and ends with a kNeighbor digi
+  for (int iRow = 0; iRow < nRows; ++iRow) {
+    int colStart = 0;
+    while (digiMap[iRow][colStart] == nullptr)
+      ++colStart;
+    if (digiMap[iRow][colStart]->GetTriggerType() != static_cast<Int_t>(CbmTrdDigi::eTriggerType::kNeighbor))
+      return false;
+
+    int colStop = nCols - 1;
+    while (digiMap[iRow][colStop] == nullptr)
+      --colStop;
+    if (digiMap[iRow][colStop]->GetTriggerType() != static_cast<Int_t>(CbmTrdDigi::eTriggerType::kNeighbor))
+      return false;
+  }
+
+  return true;
+}
+
 ClassImp(CbmTrdModuleRecR)
diff --git a/reco/detectors/trd/CbmTrdModuleRecR.h b/reco/detectors/trd/CbmTrdModuleRecR.h
index c54c0edc5e..f04e2cb11c 100644
--- a/reco/detectors/trd/CbmTrdModuleRecR.h
+++ b/reco/detectors/trd/CbmTrdModuleRecR.h
@@ -41,6 +41,7 @@ public:
 
   Int_t GetOverThreshold() const { return fDigiCounter; }
   Double_t GetSpaceResolution(Double_t val = 3.0);
+  bool IsClusterComplete(const CbmTrdCluster* cluster);
   /**
    * \brief Steering routine for building hits
    **/
-- 
GitLab