CbmRichGeoHandler.cxx

Go to the documentation of this file.

/* Copyright (C) 2015-2025 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 */

 
#include "CbmRichGeoHandler.h"
 
#include "CbmRichAddress.h"
#include "CbmRichDetectorData.h"  // for CbmRichPmtData, CbmRichPixelData
 
#include <Logger.h>
 
#include <TGeoBBox.h>
#include <TGeoManager.h>
#include <TGeoMatrix.h>
#include <TGeoNode.h>
#include <TGeoVolume.h>
#include <TRandom.h>
 
#include <string>
#include <tuple>
 
void CbmRichGeoHandler::Init()
{
  // Temporary map storing local x and y coordinates of pixels in its PMT volume
  std::map<Int_t, std::tuple<Double_t, Double_t>> localPixelCoord{};
 
  fPixelPathToAddress.clear();
  fPixelAddressToData.clear();
  fPixelAddresses.clear();
 
  fPmtPathToId.clear();
  fPmtIdToData.clear();
  fPmtIds.clear();
 
  Int_t currentPixelAddress{0};
  Int_t currentPmtId{0};
 
  // Get MAPMT dimensions
  Double_t pmtHeight{5.2};
  Double_t pmtWidth{5.2};
  TGeoVolume* pmtVolume = gGeoManager->FindVolumeFast("pmt");
  if (!pmtVolume) pmtVolume = gGeoManager->FindVolumeFast("pmt_vol_0");
  if (pmtVolume) {
    auto shape = static_cast<const TGeoBBox*>(pmtVolume->GetShape());
    if (shape) {
      pmtHeight = 2. * shape->GetDY();
      pmtWidth  = 2. * shape->GetDX();
    }
  }
 
  const TString pixelNameStr{"pmt_pixel"};  // MAPMT pixel search string
 
  TGeoIterator geoIterator{gGeoManager->GetTopNode()->GetVolume()};
  geoIterator.SetTopName("/cave_1");
  geoIterator.Reset();
  TGeoNode* curNode{nullptr};
  while ((curNode = geoIterator())) {
    if (!TString(curNode->GetVolume()->GetName()).Contains(pixelNameStr)) continue;
 
    const TGeoMatrix* curMatrix   = geoIterator.GetCurrentMatrix();
    const Double_t* curNodeTr     = curMatrix->GetTranslation();
    const TGeoMatrix* localMatrix = curNode->GetMatrix();  // local transformation of pixel in MAPMT volume
 
    TString nodePath{};
    geoIterator.GetPath(nodePath);
 
    // Check if this is a fallback version < v25
    // for deriving the pixel address from the node path
    auto vPos = nodePath.Index("rich_v");
    if (vPos != TString::kNPOS) {
      auto versionStr = std::string(nodePath.Data() + vPos + 6, 2);
      int version     = std::stoi(versionStr);
      fFallback       = version < 25;
    }
 
    const std::string path = std::string(nodePath.Data());
    auto pmtInd            = path.find_last_of('/');
    if (std::string::npos == pmtInd) continue;
    std::string pmtPath = path.substr(0, pmtInd + 1);
 
    // Derive pixel address and pmt Id from node path
    auto setup = GetDetectorSetup(nodePath);
    switch (setup) {
      case CbmRichGeoSetup::Rich:  // RICH detector
        if (fFallback) {
          currentPixelAddress++;
          break;
        }
        if (CreateAddressRich(path, currentPixelAddress, currentPmtId)) {
          LOG(error) << "CbmRichGeoHandler::Init: CbmRichGeoHandler::CreateAddressRich failed!";
        }
        break;
      case CbmRichGeoSetup::MiniRich:  // mRICH detector
        if (CreateAddressMiniRich(path, currentPixelAddress, currentPmtId)) {
          LOG(error) << "CbmRichGeoHandler::Init: CbmRichGeoHandler::CreateAddressMiniRich failed!";
          continue;
        }
        break;
      default: LOG(error) << "ERROR: Could not identify Detector setup!"; break;
    }
 
    // Create mappings, PixelData and PmtData
    fPixelPathToAddress.insert(std::pair<std::string, Int_t>(path, currentPixelAddress));
    auto pixelData      = std::make_unique<CbmRichPixelData>();
    pixelData->fX       = curNodeTr[0];
    pixelData->fY       = curNodeTr[1];
    pixelData->fZ       = curNodeTr[2];
    pixelData->fAddress = currentPixelAddress;
    localPixelCoord[currentPixelAddress] =
      std::make_tuple(localMatrix->GetTranslation()[0], localMatrix->GetTranslation()[1]);
    fPixelAddresses.push_back(currentPixelAddress);
 
    if (fPmtPathToId.count(pmtPath) == 0) {
      if (fFallback && setup == CbmRichGeoSetup::Rich) currentPmtId++;
      fPmtPathToId.insert(std::pair<std::string, Int_t>(pmtPath, currentPmtId));
      auto pmtData = std::make_unique<CbmRichPmtData>();
      pmtData->fId = currentPmtId;
      pmtData->fPixelAddresses.push_back(pixelData->fAddress);
      pmtData->fHeight  = pmtHeight;
      pmtData->fWidth   = pmtWidth;
      pixelData->fPmtId = pmtData->fId;
 
      fPmtIds.push_back(pmtData->fId);
      fPmtIdToData.insert(std::make_pair(pmtData->fId, std::move(pmtData)));
    }
    else {
      Int_t pmtId  = fPmtPathToId[pmtPath];
      auto pmtData = fPmtIdToData[pmtId].get();
      if (!pmtData || pmtId != pmtData->fId) LOG(error) << "(!pmtData || pmtId != pmtData->fId) ";
      pmtData->fPixelAddresses.push_back(pixelData->fAddress);
      pixelData->fPmtId = pmtData->fId;
      if (pmtData->fPixelAddresses.size() > 64) {
        LOG(info) << "size:" << pmtData->fPixelAddresses.size() << " pmtData->fId:" << pmtData->fId
                  << " pmtPath:" << pmtPath << std::endl
                  << " path:" << path;
      }
    }
 
    fPixelAddressToData.insert(std::make_pair(pixelData->fAddress, std::move(pixelData)));
  }
 
  // Calculate Pmt center as center of gravity of 64 pixels
  for (const auto& pmt : fPmtIdToData) {
    auto pmtData = pmt.second.get();
    pmtData->fX  = 0.;
    pmtData->fY  = 0.;
    pmtData->fZ  = 0.;
    for (int pixelId : pmtData->fPixelAddresses) {
      auto pixelData = fPixelAddressToData[pixelId].get();
      if (!pixelData) continue;
      pmtData->fX += pixelData->fX;
      pmtData->fY += pixelData->fY;
      pmtData->fZ += pixelData->fZ;
    }
    Double_t nPixels = static_cast<Double_t>(pmtData->fPixelAddresses.size());
    pmtData->fX /= nPixels;
    pmtData->fY /= nPixels;
    pmtData->fZ /= nPixels;
  }
 
  // Calculate CbmRichPixelData.fPixelId for each pixel by sorting x and y translation inside PMT volume
  for (const auto& pmt : fPmtIdToData) {
    auto pmtData = pmt.second.get();
    std::vector<std::tuple<CbmRichPixelData*, Double_t, Double_t>> pixelsInPmt{};
    for (int pixelAddress : pmtData->fPixelAddresses) {
      auto pixelData = fPixelAddressToData[pixelAddress].get();
      if (!pixelData) continue;
      pixelsInPmt.emplace_back(std::make_tuple(pixelData, std::get<0>(localPixelCoord[pixelAddress]),
                                               std::get<1>(localPixelCoord[pixelAddress])));
    }
 
    std::sort(
      pixelsInPmt.begin(), pixelsInPmt.end(),
      [](std::tuple<CbmRichPixelData*, Double_t, Double_t> a, std::tuple<CbmRichPixelData*, Double_t, Double_t> b) {
        return (std::get<2>(a) > std::get<2>(b))
               || ((std::abs(std::get<2>(a) - std::get<2>(b)) <= 0.3) && (std::get<1>(a) < std::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 (std::size_t i = 0; i < pixelsInPmt.size(); i++) {
      std::get<0>(pixelsInPmt[i])->fPixelId = static_cast<Int_t>(i);
    }
  }
 
  localPixelCoord.clear();
 
  LOG(info) << "CbmRichGeoHandler is initialized";
  LOG(info) << "fPixelPathToAddress.size() = " << fPixelPathToAddress.size();
  LOG(info) << "fPixelAddressToData.size() = " << fPixelAddressToData.size();
 
  LOG(info) << "fPmtPathToId.size() = " << fPmtPathToId.size();
  LOG(info) << "fPmtIdToData.size() = " << fPmtIdToData.size();
}
 
Int_t CbmRichGeoHandler::GetPixelAddressByPath(const std::string& path)
{
  return fPixelPathToAddress.count(path) ? fPixelPathToAddress[path] : -1;
}
 
CbmRichPixelData* CbmRichGeoHandler::GetPixelDataByAddress(Int_t address)
{
  if (fPixelAddressToData.count(address) == 0) {
    LOG(fatal) << "CbmRichGeoHandler::GetPixelDataByAddress: Pixel with address " << address << " not found.";
    return nullptr;
  }
  return fPixelAddressToData[address].get();
}
 
Int_t CbmRichGeoHandler::GetRandomPixelAddress() { return fPixelAddresses[gRandom->Integer(fPixelAddresses.size())]; }
 
CbmRichPmtData* CbmRichGeoHandler::GetPmtDataById(Int_t id)
{
  return fPmtIdToData.count(id) ? fPmtIdToData[id].get() : nullptr;
}
 
std::vector<Int_t> CbmRichGeoHandler::GetNeighbourPixels(Int_t address, Int_t n, Bool_t horizontal, Bool_t vertical,
                                                         Bool_t diagonal)
{
  std::vector<Int_t> neighbourPixels{};
  if (n == 0) return neighbourPixels;
  CbmRichPixelData* addressPixelData   = GetPixelDataByAddress(address);
  Int_t indX                           = addressPixelData->fPixelId % 8;
  Int_t indY                           = addressPixelData->fPixelId / 8;
  std::vector<Int_t> pmtPixelAddresses = GetPmtDataById(addressPixelData->fPmtId)->fPixelAddresses;
  for (const auto& iAddr : pmtPixelAddresses) {
    if (iAddr == address) continue;
    CbmRichPixelData* iPixelData = GetPixelDataByAddress(iAddr);
    Int_t iIndX                  = iPixelData->fPixelId % 8;
    Int_t iIndY                  = iPixelData->fPixelId / 8;
    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 (!horizontal && !vertical && diagonal && n == 1) {  // diagonal neighbours
      if ((abs(iIndX - indX) == 1) && (abs(iIndY - indY) == 1)) neighbourPixels.push_back(iAddr);
    }
    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 CbmRichGeoHandler::GetNeighbourPixels " << std::endl
                 << " n = " << n << " horizontal = " << horizontal << " vertical = " << vertical
                 << " diagonal = " << diagonal;
    }
  }
  return neighbourPixels;
}
 
Int_t CbmRichGeoHandler::CreateAddressRich(const std::string& nodePath, Int_t& channelAddr, Int_t& pmtId)
{
  // Example path for RICH v25
  // /cave_1/rich_v25a_0/rich_cont_1/rich_gas_1/camera_cont_1/camera_strip_1/camera_module_7/pmt_1/pmt_pixel_corner_1
  // Indexing starts from 1 in geo paths, decrease by 1 for creating the RichAddress
 
  auto pmtPixelInd = nodePath.find_last_of('/');
  if (std::string::npos == pmtPixelInd) return 1;
 
  // "/pmt_pixel_*_[1-64]" -> Pixel number [0, 63]
  auto pixelUnderscore = nodePath.find_last_of('_');
  if (std::string::npos == pixelUnderscore) return 1;
  uint32_t pixelId = std::stoul(nodePath.substr(pixelUnderscore + 1)) - 1;
 
  // "/pmt_[1-6]" -> Pmt number [0, 5] in backplane
  auto pmtContInd = nodePath.rfind('/', pmtPixelInd - 1);
  uint32_t pmt    = std::stoul(nodePath.substr(pmtContInd + 5, pmtPixelInd - pmtContInd - 5)) - 1;
  if (std::string::npos == pmtContInd) return 1;
 
  // "/camera_module_[1-7]" -> Backplane number [0, 6]
  auto backplaneInd  = nodePath.rfind('/', pmtContInd - 1);
  uint32_t backplane = std::stoul(nodePath.substr(backplaneInd + 15, pmtContInd - backplaneInd - 15)) - 1;
  if (std::string::npos == backplaneInd) return 1;
 
  // "/camera_strip_[1-14]" -> Strip number [0, 13]
  auto stripInd  = nodePath.rfind('/', backplaneInd - 1);
  uint32_t strip = std::stoul(nodePath.substr(stripInd + 14, backplaneInd - stripInd - 14)) - 1;
  if (std::string::npos == stripInd) return 1;
 
  // "/camera_cont_[1-2]" -> Camera number [0: upper, 1: lower camera]
  auto cameraInd  = nodePath.rfind('/', stripInd - 1);
  uint32_t camera = std::stoul(nodePath.substr(cameraInd + 13, stripInd - cameraInd - 13)) - 1;
  if (std::string::npos == cameraInd) return 1;
 
  channelAddr = CbmRichAddress::GetAddress(camera, strip, backplane, pmt, pixelId);
  pmtId       = CbmRichAddress::GetMotherAddress(channelAddr, 4);
 
  return 0;
}
 
Int_t CbmRichGeoHandler::CreateAddressMiniRich(const std::string& nodePath, Int_t& channelAddr, Int_t& pmtId)
{
  auto pmtInd = nodePath.find_last_of('/');
  if (std::string::npos == pmtInd) return 1;
 
  Int_t channel  = static_cast<Int_t>(std::stoul(nodePath.substr(pmtInd + 11)));  // cut away "pmt_pixel_"
  Int_t posAtPmt = channel / 100;
  channel        = channel % 100;
 
  auto pmtVolInd = nodePath.rfind('/', pmtInd - 1);
  if (std::string::npos == pmtVolInd) return 1;
  Int_t pmtPosBP = static_cast<Int_t>(
    std::stoul(nodePath.substr(pmtVolInd + 11,
                               pmtInd - pmtVolInd - 11)));  // cut away "/pmt_vol_*_" ; position on BP
 
  auto bpInd = nodePath.rfind('/', pmtVolInd - 1);
  if (std::string::npos == bpInd) return 1;
  Int_t posBP = static_cast<Int_t>(
    std::stoul(nodePath.substr(bpInd + 14,
                               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));
 
  Int_t DiRICH_Add = ((7 & 0xF) << 12) + ((x & 0xF) << 8) + ((y & 0xF) << 4) + (posAtPmt & 0xF);
  Int_t pixelUID   = ((DiRICH_Add << 16) | (channel & 0x00FF));
  Int_t pmtUID     = ((x & 0xF) << 4) + (y & 0xF);
 
  channelAddr = pixelUID;
  pmtId       = pmtUID;
 
  return 0;
}