Rich: add index to each PMT pixel by sorting translations inside its pmt volume

core/detectors/rich/CbmRichDetectorData.h

-/* Copyright (C) 2015-2020 GSI/JINR-LIT, Darmstadt/Dubna
+/* Copyright (C) 2015-2023 GSI/JINR-LIT, Darmstadt/Dubna
    SPDX-License-Identifier: GPL-3.0-only
    Authors: Semen Lebedev, Andrey Lebedev [committer], Florian Uhlig */
 
@@ -21,9 +21,10 @@
 class CbmRichPixelData {
 public:
   Int_t fAddress;
-  Double_t fX;
-  Double_t fY;
-  Double_t fZ;
+  Double_t fX;     // global x coordinate
+  Double_t fY;     // global y coordinate
+  Double_t fZ;     // global z coordinate
+  Int_t fPixelId;  // pixel index [0, 63] row major, view from -z to z
   Int_t fPmtId;
 };
 

core/detectors/rich/CbmRichDigiMapManager.cxx

-/* Copyright (C) 2015-2020 GSI/JINR-LIT, Darmstadt/Dubna
+/* Copyright (C) 2015-2023 GSI/JINR-LIT, Darmstadt/Dubna
    SPDX-License-Identifier: GPL-3.0-only
-   Authors: Semen Lebedev, Adrian Amatus Weber, Andrey Lebedev [committer], Florian Uhlig */
+   Authors: Semen Lebedev, Martin Beyer, Adrian Amatus Weber, Andrey Lebedev [committer], Florian Uhlig */
 
 /*
  * CbmRichDigiMapManager.cxx
@@ -22,6 +22,7 @@
 #include <TRandom.h>      // for TRandom, gRandom
 
 #include <string>   // for operator<, stoul
+#include <tuple>    // for tuple, make_tuple, get
 #include <utility>  // for pair
 
 #include <stddef.h>  // for size_t
@@ -43,6 +44,8 @@ CbmRichDigiMapManager::~CbmRichDigiMapManager() {}
 
 void CbmRichDigiMapManager::Init()
 {
+  // temporary map storing local x and y coordinates of pixels in its PMT volume
+  std::map<Int_t, tuple<Double_t, Double_t>> localPixelCoord;
 
   fPixelPathToAddressMap.clear();
   fPixelAddressToDataMap.clear();
@@ -87,6 +90,7 @@ void CbmRichDigiMapManager::Init()
           //std::cout<<"RICH"<<std::endl;
           const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix();
           const Double_t* curNodeTr   = curMatrix->GetTranslation();
+          const TGeoMatrix* localMatrix = curNode->GetMatrix();  // local transformation of pixel in MAPMT volume
           string path                 = string(nodePath.Data());
 
           size_t pmtInd = path.find_last_of("/");
@@ -99,6 +103,8 @@ void CbmRichDigiMapManager::Init()
           pixelData->fY               = curNodeTr[1];
           pixelData->fZ               = curNodeTr[2];
           pixelData->fAddress         = currentPixelAddress;
+          localPixelCoord[currentPixelAddress] =
+            make_tuple(localMatrix->GetTranslation()[0], localMatrix->GetTranslation()[1]);
           fPixelAddressToDataMap.insert(pair<Int_t, CbmRichPixelData*>(pixelData->fAddress, pixelData));
           fPixelAddresses.push_back(currentPixelAddress);
           currentPixelAddress++;
@@ -139,6 +145,7 @@ void CbmRichDigiMapManager::Init()
           //std::cout<<"mRICH"<<std::endl;
           const TGeoMatrix* curMatrix = geoIterator.GetCurrentMatrix();
           const Double_t* curNodeTr   = curMatrix->GetTranslation();
+          const TGeoMatrix* localMatrix = curNode->GetMatrix();  // local transformation of pixel in MAPMT volume
           string path                 = string(nodePath.Data());
 
           size_t pmtInd = path.find_last_of("/");
@@ -157,7 +164,7 @@ void CbmRichDigiMapManager::Init()
           size_t bpInd = path.rfind("/", pmtVolInd - 1);
           if (string::npos == bpInd) continue;
           Int_t posBP = std::stoul(path.substr(bpInd + 14,
-                                               pmtVolInd - bpInd - 14));  // cut away "/pmt_vol_*_" ; position on BP
+                                               pmtVolInd - bpInd - 14));  // cut away "/pmt_cont_vol_" ; position on BP
 
           Int_t x = (posBP / 10) + ((((pmtPosBP - 1) / 3) + 1) % 2);
           Int_t y = (posBP % 10) + (2 - ((pmtPosBP - 1) % 3));
@@ -173,6 +180,7 @@ void CbmRichDigiMapManager::Init()
           pixelData->fY               = curNodeTr[1];
           pixelData->fZ               = curNodeTr[2];
           pixelData->fAddress         = pixelUID;
+          localPixelCoord[pixelUID]   = make_tuple(localMatrix->GetTranslation()[0], localMatrix->GetTranslation()[1]);
           fPixelAddressToDataMap.insert(pair<Int_t, CbmRichPixelData*>(pixelData->fAddress, pixelData));
           fPixelAddresses.push_back(pixelUID);
           //currentPixelAddress++;
@@ -232,6 +240,37 @@ void CbmRichDigiMapManager::Init()
     pmtData->fZ /= pmtData->fPixelAddresses.size();
   }
 
+  // calculate fPixelId for each pixel by sorting x and y translation inside PMT volume
+  for (auto const& pmt : fPmtIdToDataMap) {
+    CbmRichPmtData* pmtData = pmt.second;
+    std::vector<tuple<CbmRichPixelData*, Double_t, Double_t>> pixelsInPmt;
+    for (int pixelAddress : pmtData->fPixelAddresses) {
+      CbmRichPixelData* pixelData = fPixelAddressToDataMap[pixelAddress];
+      if (pixelData == nullptr) continue;
+      pixelsInPmt.push_back(
+        make_tuple(pixelData, get<0>(localPixelCoord[pixelAddress]), get<1>(localPixelCoord[pixelAddress])));
+    }
+
+    std::sort(
+      pixelsInPmt.begin(), pixelsInPmt.end(),
+      [](tuple<CbmRichPixelData*, Double_t, Double_t> a, tuple<CbmRichPixelData*, Double_t, Double_t> b) {
+        return (get<2>(a) > get<2>(b)) || ((abs(get<2>(a) - get<2>(b)) <= 0.3) && (get<1>(a) < get<1>(b)));
+        // first sort by y (up to down) coordinate, then by x (left to right) coordinate
+        // stop sorting by y when y coordinate difference is less than 0.3 (to start sorting by x when y is the "same")
+        // needed because edge/corner pixels are slightly larger
+      });
+
+    if (pixelsInPmt.size() != 64) {
+      LOG(error) << "ERROR: Calculating local pixel indices failed, number of pixels in PMT is not 64. "
+                 << pmtData->ToString();
+    }
+    for (unsigned int i = 0; i < pixelsInPmt.size(); i++) {
+      get<0>(pixelsInPmt[i])->fPixelId = i;
+    }
+  }
+
+  localPixelCoord.clear();
+
   LOG(info) << "CbmRichDigiMapManager is initialized";
   LOG(info) << "fPixelPathToAddressMap.size() = " << fPixelPathToAddressMap.size();
   LOG(info) << "fPixelAddressToDataMap.size() = " << fPixelAddressToDataMap.size();
@@ -283,18 +322,33 @@ CbmRichPmtData* CbmRichDigiMapManager::GetPmtDataById(Int_t id)
   return it->second;
 }
 
-vector<Int_t> CbmRichDigiMapManager::GetDirectNeighbourPixels(Int_t /*address*/)
+vector<Int_t> CbmRichDigiMapManager::GetNeighbourPixels(Int_t address, Int_t N, Bool_t direct, Bool_t diagonal)
 {
-  std::vector<Int_t> v;
-
-  return v;
-}
-
-vector<Int_t> CbmRichDigiMapManager::GetDiagonalNeighbourPixels(Int_t /*address*/)
-{
-  std::vector<Int_t> v;
-
-  return v;
+  vector<Int_t> neighbourPixels;
+  CbmRichPixelData* addressPixelData = GetPixelDataByAddress(address);
+  Int_t indX                         = addressPixelData->fPixelId % 8;
+  Int_t indY                         = addressPixelData->fPixelId / 8;
+  vector<Int_t> pmtPixelAddresses    = GetPmtDataById(addressPixelData->fPmtId)->fPixelAddresses;
+  for (auto const& iAddr : pmtPixelAddresses) {
+    if (iAddr == address) continue;
+    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);
+    }
+    else if (!direct && 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
+      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;
+    }
+  }
+  return neighbourPixels;
 }
 
 int CbmRichDigiMapManager::getDetectorSetup(TString const nodePath)

core/detectors/rich/CbmRichDigiMapManager.h

-/* Copyright (C) 2015-2020 GSI/JINR-LIT, Darmstadt/Dubna
+/* Copyright (C) 2015-2023 GSI/JINR-LIT, Darmstadt/Dubna
    SPDX-License-Identifier: GPL-3.0-only
-   Authors: Semen Lebedev, Andrey Lebedev [committer], Florian Uhlig */
+   Authors: Semen Lebedev, Martin Beyer, Andrey Lebedev [committer], Florian Uhlig */
 
 /*
  * CbmRichDigiMap.h
@@ -28,8 +28,8 @@ private:
 
 public:
   /**
-	 * Return Instance of CbmRichGeoManager.
-	 */
+   * Return Instance of CbmRichGeoManager.
+   */
   static CbmRichDigiMapManager& GetInstance()
   {
     static CbmRichDigiMapManager fInstance;
@@ -37,47 +37,60 @@ public:
   }
 
   /*
-	 * \brief Return digi address by path to node.
-	 */
+   * \brief Return digi address by path to node.
+   */
   Int_t GetPixelAddressByPath(const std::string& path);
 
   /*
-	 * \brief Return CbmRichDataPixel by digi address.
-	 */
+   * \brief Return CbmRichDataPixel by digi address.
+   */
   CbmRichPixelData* GetPixelDataByAddress(Int_t address);
 
   /*
-	 * \brief Return the addresses of the direct neighbour pixels.
-	 * C++11 efficient way to return vector
-	 */
-  std::vector<Int_t> GetDirectNeighbourPixels(Int_t address);
+   * \brief Return random address. Needed for noise digi.
+   */
+  Int_t GetRandomPixelAddress();
 
   /*
-	 * \brief Return the addresses of the diagonal neighbour pixels.
-	 * C++11 efficient way to return vector
-	 */
-  std::vector<Int_t> GetDiagonalNeighbourPixels(Int_t address);
+   * \brief Return addresses of all pixels
+   */
+  std::vector<Int_t> GetPixelAddresses();
 
   /*
-	 * \brief Return random address. Needed for noise digi.
-	 */
-  Int_t GetRandomPixelAddress();
+   * \brief Return ids for all pmts
+   */
+  std::vector<Int_t> GetPmtIds();
 
+  /*
+   * \brief Return CbmRichDataPmt by id.
+   */
+  CbmRichPmtData* GetPmtDataById(Int_t id);
 
   /*
-     * \brief Return addresses of all pixels
-     */
-  std::vector<Int_t> GetPixelAddresses();
+   * \brief Return the addresses of the neighbour pixels.
+   */
+  std::vector<Int_t> GetNeighbourPixels(Int_t address, Int_t N, Bool_t direct = true, Bool_t diagonal = true);
 
   /*
-     * \brief Return ids for all pmts
-     */
-  std::vector<Int_t> GetPmtIds();
+   * \brief Return the addresses of the direct neighbour pixels.
+   */
+  std::vector<Int_t> GetDirectNeighbourPixels(Int_t address) { return GetNeighbourPixels(address, 1, true, false); }
 
   /*
-     * \brief Return CbmRichDataPmt by id.
-     */
-  CbmRichPmtData* GetPmtDataById(Int_t id);
+   * \brief Return the addresses of the diagonal neighbour pixels.
+   */
+  std::vector<Int_t> GetDiagonalNeighbourPixels(Int_t address) { return GetNeighbourPixels(address, 1, 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.
+   * Needed for noise digis caused by charged particles.
+   */
+  std::vector<Int_t> GetNxNNeighbourPixels(Int_t address, Int_t N)
+  {
+    return GetNeighbourPixels(address, N, true, true);
+  }
 
 public:
   virtual ~CbmRichDigiMapManager();
@@ -94,18 +107,18 @@ private:
   int getDetectorSetup(TString const nodePath);
 
   /*
-	 * \brief Initialize maps.
-	 */
+   * \brief Initialize maps.
+   */
   void Init();
 
   /**
-	* \brief Copy constructor.
-	*/
+   * \brief Copy constructor.
+   */
   CbmRichDigiMapManager(const CbmRichDigiMapManager&);
 
   /**
-	* \brief Assignment operator.
-	*/
+   * \brief Assignment operator.
+   */
   CbmRichDigiMapManager& operator=(const CbmRichDigiMapManager&);
 };