Rich: make Digitizer crosstalk probabilities more flexible

core/detectors/rich/CbmRichDigiMapManager.cxx

-/* Copyright (C) 2015-2023 GSI/JINR-LIT, Darmstadt/Dubna
+/* Copyright (C) 2015-2024 GSI/JINR-LIT, Darmstadt/Dubna
    SPDX-License-Identifier: GPL-3.0-only
    Authors: Semen Lebedev, Martin Beyer, Adrian Amatus Weber, Andrey Lebedev [committer], Florian Uhlig */
 
@@ -322,7 +322,8 @@ CbmRichPmtData* CbmRichDigiMapManager::GetPmtDataById(Int_t id)
   return it->second;
 }
 
-vector<Int_t> CbmRichDigiMapManager::GetNeighbourPixels(Int_t address, Int_t n, Bool_t direct, Bool_t diagonal)
+vector<Int_t> CbmRichDigiMapManager::GetNeighbourPixels(Int_t address, Int_t n, Bool_t horizontal, Bool_t vertical,
+                                                        Bool_t diagonal)
 {
   vector<Int_t> neighbourPixels;
   if (n == 0) return neighbourPixels;
@@ -335,18 +336,25 @@ vector<Int_t> CbmRichDigiMapManager::GetNeighbourPixels(Int_t address, Int_t n,
     CbmRichPixelData* iPixelData = GetPixelDataByAddress(iAddr);
     Int_t iIndX                  = iPixelData->fPixelId % 8;
     Int_t iIndY                  = iPixelData->fPixelId / 8;
-    if (direct && !diagonal && n == 1) {  // direct neighbours
-      if ((abs(iIndX - indX) + abs(iIndY - indY)) == 1) neighbourPixels.push_back(iAddr);
+    if (horizontal && !vertical && !diagonal && n == 1) {  // direct horizontal neighbours
+      if (abs(iIndX - indX) == 1 && abs(iIndY - indY) == 0) neighbourPixels.push_back(iAddr);
+    }
+    else if (vertical && !horizontal && !diagonal && n == 1) {  // direct vertical neighbours
+      if (abs(iIndX - indX) == 0 && abs(iIndY - indY) == 1) neighbourPixels.push_back(iAddr);
     }
-    else if (!direct && diagonal && n == 1) {  // diagonal neighbours
+    else if (!horizontal && !vertical && diagonal && n == 1) {  // diagonal neighbours
       if ((abs(iIndX - indX) == 1) && (abs(iIndY - indY) == 1)) neighbourPixels.push_back(iAddr);
     }
-    else if (direct && diagonal) {  // all neighbours in (2N+1)*(2N+1) grid
+    else if (horizontal && vertical && !diagonal && n == 1) {  // direct horizontal and vertical neighbours
+      if ((abs(iIndX - indX) + abs(iIndY - indY)) == 1) neighbourPixels.push_back(iAddr);
+    }
+    else if (horizontal && vertical && diagonal) {  // all neighbours in (2N+1)*(2N+1) grid
       if ((abs(iIndX - indX) <= n) && (abs(iIndY - indY) <= n)) neighbourPixels.push_back(iAddr);
     }
     else {
       LOG(error) << "ERROR: Unrecogniced option in CbmRichDigiMapManager::GetNeighbourPixels " << endl
-                 << " n = " << n << " direct = " << direct << " diagonal = " << diagonal;
+                 << " n = " << n << " horizontal = " << horizontal << " vertical = " << vertical
+                 << " diagonal = " << diagonal;
     }
   }
   return neighbourPixels;

core/detectors/rich/CbmRichDigiMapManager.h

-/* Copyright (C) 2015-2023 GSI/JINR-LIT, Darmstadt/Dubna
+/* Copyright (C) 2015-2024 GSI/JINR-LIT, Darmstadt/Dubna
    SPDX-License-Identifier: GPL-3.0-only
    Authors: Semen Lebedev, Martin Beyer, Andrey Lebedev [committer], Florian Uhlig */
 
@@ -23,12 +23,12 @@ class CbmRichPixelData;
 class CbmRichPmtData;
 
 class CbmRichDigiMapManager {
-private:
+ private:
   CbmRichDigiMapManager();
 
-public:
+ public:
   /**
-   * Return Instance of CbmRichGeoManager.
+   * \brief Return Instance of CbmRichGeoManager.
    */
   static CbmRichDigiMapManager& GetInstance()
   {
@@ -36,58 +36,66 @@ public:
     return fInstance;
   }
 
-  /*
+  /**
    * \brief Return digi address by path to node.
    */
   Int_t GetPixelAddressByPath(const std::string& path);
 
-  /*
+  /**
    * \brief Return CbmRichDataPixel by digi address.
    */
   CbmRichPixelData* GetPixelDataByAddress(Int_t address);
 
-  /*
+  /**
    * \brief Return random address. Needed for noise digi.
    */
   Int_t GetRandomPixelAddress();
 
-  /*
+  /**
    * \brief Return addresses of all pixels
    */
   std::vector<Int_t> GetPixelAddresses();
 
-  /*
+  /**
    * \brief Return ids for all pmts
    */
   std::vector<Int_t> GetPmtIds();
 
-  /*
+  /**
    * \brief Return CbmRichDataPmt by id.
    */
   CbmRichPmtData* GetPmtDataById(Int_t id);
 
-  /*
+  /**
    * \brief Return the addresses of the neighbour pixels.
    * \param address Pixel address
    * \param n Size of the grid (2n+1)*(2n+1)
-   * \param direct return direct neighbours
+   * \param horizontal return horizontal neighbours
+   * \param vertical return vertical neighbours
    * \param diagonal return diagonal neighbours
    */
-  std::vector<Int_t> GetNeighbourPixels(Int_t address, Int_t N, Bool_t direct = true, Bool_t diagonal = true);
+  std::vector<Int_t> GetNeighbourPixels(Int_t address, Int_t N, Bool_t horizontal = true, Bool_t vertical = true,
+                                        Bool_t diagonal = true);
 
-  /*
+  /**
    * \brief Return the addresses of the direct neighbour pixels.
    * \param address Pixel address
    */
-  std::vector<Int_t> GetDirectNeighbourPixels(Int_t address) { return GetNeighbourPixels(address, 1, true, false); }
+  std::vector<Int_t> GetDirectNeighbourPixels(Int_t address, Bool_t horizontal = true, Bool_t vertical = true)
+  {
+    return GetNeighbourPixels(address, 1, horizontal, vertical, false);
+  }
 
-  /*
+  /**
    * \brief Return the addresses of the diagonal neighbour pixels.
    * \param address Pixel address
    */
-  std::vector<Int_t> GetDiagonalNeighbourPixels(Int_t address) { return GetNeighbourPixels(address, 1, false, true); }
+  std::vector<Int_t> GetDiagonalNeighbourPixels(Int_t address)
+  {
+    return GetNeighbourPixels(address, 1, false, false, true);
+  }
 
-  /*
+  /**
    * \brief Return the addresses of pixels in a (2n+1)*(2n+1) grid,
    * with the address pixel in the center of the grid.
    * Addresses are limited to the same MAPMT as the input address.
@@ -97,13 +105,13 @@ public:
    */
   std::vector<Int_t> GetNxNNeighbourPixels(Int_t address, Int_t n)
   {
-    return GetNeighbourPixels(address, n, true, true);
+    return GetNeighbourPixels(address, n, true, true, true);
   }
 
-public:
+ public:
   virtual ~CbmRichDigiMapManager();
 
-private:
+ private:
   std::map<std::string, Int_t> fPixelPathToAddressMap;
   std::map<Int_t, CbmRichPixelData*> fPixelAddressToDataMap;
   std::vector<Int_t> fPixelAddresses;  // vector of all  pixel addresses
@@ -114,7 +122,7 @@ private:
 
   int getDetectorSetup(TString const nodePath);
 
-  /*
+  /**
    * \brief Initialize maps.
    */
   void Init();

reco/detectors/rich/qa/CbmRichDigiQa.cxx

@@ -83,14 +83,14 @@ void CbmRichDigiQa::Exec(Option_t* /*option*/)
     }
   }
 
-  // Fill crosstalk digis
+  // Fill neighbour digis
   for (auto& [pmt, digis] : fPmtDigisTimeAddress) {
     std::sort(digis.begin(), digis.end());
     for (auto it0 = digis.begin(); it0 != digis.end(); ++it0) {
       Int_t indX0 = CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(it0->second)->fPixelId % 8;
       Int_t indY0 = CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(it0->second)->fPixelId / 8;
       for (auto it1 = it0 + 1; it1 != digis.end(); ++it1) {
-        if (abs(it0->first - it1->first) < fCrosstalkTimeLimit) {
+        if (abs(it0->first - it1->first) < fNeighbourTimeLimit) {
           Int_t indX1 = CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(it1->second)->fPixelId % 8;
           Int_t indY1 = CbmRichDigiMapManager::GetInstance().GetPixelDataByAddress(it1->second)->fPixelId / 8;
           if (indX0 == indX1 && indY0 == indY1) continue;
@@ -99,7 +99,7 @@ void CbmRichDigiQa::Exec(Option_t* /*option*/)
           if (abs(indX0 - indX1) > 1 || abs(indY0 - indY1) > 1) continue;
           fHM->H2("fhDigiNeighbours")->Fill(indX1 - indX0, indY0 - indY1);
         }
-        else if (it1->first - it0->first > fCrosstalkTimeLimit) {
+        else if (it1->first - it0->first > fNeighbourTimeLimit) {
           break;
         }
       }

reco/detectors/rich/qa/CbmRichDigiQa.h

-/* Copyright (C) 2023-2023 UGiessen, Giessen
+/* Copyright (C) 2023-2024 UGiessen, Giessen
    SPDX-License-Identifier: GPL-3.0-only
    Authors: Martin Beyer [committer] */
 
@@ -67,12 +67,12 @@ public:
   }
 
   /**
-    * \brief Set time limit in which a neighbour digi is considered crosstalk
+    * \brief Set time limit in which a digi is considered neighbour.
     * \param[in] limit absolue time limit
     */
-  void SetCrosstalkTimeLimit(Double_t limit) { fCrosstalkTimeLimit = limit; }
+  void SetNeighbourTimeLimit(Double_t limit) { fNeighbourTimeLimit = limit; }
 
-private:
+ private:
   int fEventNum {};
 
   CbmDigiManager* fDigiMan {nullptr};
@@ -84,7 +84,7 @@ private:
   std::map<Int_t, std::vector<std::pair<Double_t, Int_t>>>
     fPmtDigisTimeAddress {};  // Store times and addresses of fired digis for each pmt in Ts
 
-  Double_t fCrosstalkTimeLimit {5.};  // Time in which a neighbour digi is considered crosstalk
+  Double_t fNeighbourTimeLimit{5.};   // Time in which digi is considered a neighbour
   Double_t fToTLimitLow {-1.};        // Lower limit for ToT cut
   Double_t fToTLimitHigh {-1.};       // Upper limit for ToT cut
 

sim/detectors/rich/CbmRichDigitizer.cxx

-/* Copyright (C) 2015-2023 GSI/JINR-LIT, Darmstadt/Dubna
+/* Copyright (C) 2015-2024 GSI/JINR-LIT, Darmstadt/Dubna
    SPDX-License-Identifier: GPL-3.0-only
    Authors: Martin Beyer, Andrey Lebedev [committer], Volker Friese, Florian Uhlig, Semen Lebedev */
 
@@ -188,7 +188,9 @@ void CbmRichDigitizer::ProcessPoint(CbmRichPoint* point, Int_t pointId, Int_t ev
     CbmLink link(1., pointId, eventNum, inputNum);
     AddSignalToBuffer(address, time, link);
     if (gcode == 50000050 && address > 0) AddCrossTalk(address, time, link);
-    if (mcTrack->GetCharge() != 0. && address > 0) AddChargedParticleCluster(address, time, eventNum, inputNum);
+    if (fClusterSize > 0 && mcTrack->GetCharge() != 0. && address > 0) {
+      AddChargedParticleCluster(address, time, eventNum, inputNum);
+    }
   }
 }
 
@@ -199,18 +201,30 @@ void CbmRichDigitizer::AddSignalToBuffer(Int_t address, Double_t time, const Cbm
 
 void CbmRichDigitizer::AddCrossTalk(Int_t address, Double_t time, const CbmLink& link)
 {
-  std::vector<Int_t> directPixels   = CbmRichDigiMapManager::GetInstance().GetDirectNeighbourPixels(address);
+  std::vector<Int_t> directHorizontalPixels =
+    CbmRichDigiMapManager::GetInstance().GetDirectNeighbourPixels(address, true, false);
+  std::vector<Int_t> directVerticalPixels =
+    CbmRichDigiMapManager::GetInstance().GetDirectNeighbourPixels(address, false, true);
   std::vector<Int_t> diagonalPixels = CbmRichDigiMapManager::GetInstance().GetDiagonalNeighbourPixels(address);
   /* clang-format off */
-  Double_t crosstalkCreatedProbability = 1 - pow(1 - fCrossTalkProbability, directPixels.size()) *
-                                             pow(1 - fCrossTalkProbability / 4., diagonalPixels.size());
+  Double_t crosstalkProbability = 1 - pow(1 - fCrossTalkProbability[0], directHorizontalPixels.size()) *
+                                      pow(1 - fCrossTalkProbability[1], directVerticalPixels.size()) *
+                                      pow(1 - fCrossTalkProbability[2], diagonalPixels.size());
   Int_t crosstalkAddress = -1;
-  if (gRandom->Uniform() < crosstalkCreatedProbability) {
-    Double_t thresholdDirectDiagonal =
-      static_cast<Double_t>(directPixels.size())
-      / (static_cast<Double_t>(directPixels.size()) + static_cast<Double_t>(diagonalPixels.size()) / 4.);
-    if (gRandom->Uniform() < thresholdDirectDiagonal) {
-      crosstalkAddress = directPixels[gRandom->Integer(directPixels.size())];
+  if (gRandom->Uniform() < crosstalkProbability) {
+    // Split into 3 intervals based on crosstalk probability and number of pixels
+    Double_t denom = directHorizontalPixels.size() * fCrossTalkProbability[0] + 
+                     directVerticalPixels.size()   * fCrossTalkProbability[1] + 
+                     diagonalPixels.size()         * fCrossTalkProbability[2];
+    // Threshold values for end of each interval
+    Double_t thHorizontal = directHorizontalPixels.size() * fCrossTalkProbability[0] / denom;
+    Double_t thVertical = directVerticalPixels.size()     * fCrossTalkProbability[1] / denom + thHorizontal;    
+    Double_t rnd = gRandom->Uniform();
+    if (rnd < thHorizontal) {
+      crosstalkAddress = directHorizontalPixels[gRandom->Integer(directHorizontalPixels.size())];
+    }
+    else if (rnd < thVertical) {
+      crosstalkAddress = directVerticalPixels[gRandom->Integer(directVerticalPixels.size())];
     } else {
       crosstalkAddress = diagonalPixels[gRandom->Integer(diagonalPixels.size())];
     }

sim/detectors/rich/CbmRichDigitizer.h

-/* Copyright (C) 2015-2023 GSI/JINR-LIT, Darmstadt/Dubna
+/* Copyright (C) 2015-2024 GSI/JINR-LIT, Darmstadt/Dubna
    SPDX-License-Identifier: GPL-3.0-only
    Authors: Semen Lebedev, Andrey Lebedev [committer], Martin Beyer */
 
@@ -76,8 +76,18 @@ public:
     */
   void SetPixelDeadTime(Double_t dt) { fPixelDeadTime = dt; }
 
-  /** Set crosstalk probability */
-  void SetCrossTalkProbability(Double_t crosstalk) { fCrossTalkProbability = crosstalk; }
+  /** 
+   * \brief Set crosstalk probability for horizontal, vertical and diagonal directions.
+   * \param horizontal probability of crosstalk for direct horizontal pixels
+   * \param vertical probability of crosstalk for direct vertical pixels
+   * \param diagonal probability of crosstalk for direct diagonal pixels
+   */
+  void SetCrossTalkProbability(Double_t horizontal, Double_t vertical, Double_t diagonal)
+  {
+    fCrossTalkProbability[0] = horizontal;
+    fCrossTalkProbability[1] = vertical;
+    fCrossTalkProbability[2] = diagonal;
+  }
 
   /**
     * \brief Set event noise rate per McRichPoint / per  pixel / per second :
@@ -130,7 +140,9 @@ private:
 
   Double_t fTimeResolution {1.};              // in ns, time resolution for digis
   Double_t fPixelDeadTime {50.};              // in ns, deadtime for one pixel
-  Double_t fCrossTalkProbability {0.02};      // probability of the crosstalk for direct neighbor for one pixel
+  // probability of crosstalk depending on direction of neighbour pixel
+  // index 0: horizontal, 1: vertical, 2: diagonal
+  std::array<Double_t, 3> fCrossTalkProbability{0.02, 0.02, 0.005};
   Double_t fEventNoiseRate {5.};              // noise rate per McRichPoint / per  pixel / per second
   Double_t fEventNoiseInterval {50.};         // in ns, interval to generate event noise signals
   Double_t fDarkRatePerPixel {1000.};         // dark rate per pixel in Hz