Add crosstalk neighbour class (#69)

* Add crosstalk neighbour class

* Use c++ static_cast

* Attempt to add automatic test

* Adapt to the new headfile name of the k4geo update

* revert neighbours to ExtSvc in the test

* change function name: xtalk_get_cell_indices

* Adapt to new k4geo function names for the crosstalk neighbours

RecFCCeeCalorimeter/CMakeLists.txt

@@ -64,3 +64,8 @@ add_test(NAME FCCeeLAr_benchmarkCalibration
 )
 set_test_env(FCCeeLAr_benchmarkCalibration)
 
+add_test(NAME FCCeeLAr_xtalkNeighbours
+  COMMAND k4run RecFCCeeCalorimeter/tests/options/runCaloXTalkNeighbours.py
+         WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
+)
+set_test_env(FCCeeLAr_xtalkNeighbours)

RecFCCeeCalorimeter/src/components/CreateFCCeeCaloXTalkNeighbours.cpp

@@ -0,0 +1,256 @@
+#include "CreateFCCeeCaloXTalkNeighbours.h"
+
+// DD4hep
+#include "DD4hep/Detector.h"
+
+// k4FWCore
+#include "k4Interface/IGeoSvc.h"
+
+// k4geo
+#include "detectorCommon/DetUtils_k4geo.h"
+#include "detectorSegmentations/GridTheta_k4geo.h"
+#include "detectorSegmentations/FCCSWGridPhiTheta_k4geo.h"
+#include "detectorSegmentations/FCCSWGridModuleThetaMerged_k4geo.h"
+
+// ROOT
+#include "TFile.h"
+#include "TTree.h"
+
+#define DEBUGCELLINFO 1
+
+DECLARE_COMPONENT(CreateFCCeeCaloXTalkNeighbours)
+
+CreateFCCeeCaloXTalkNeighbours::CreateFCCeeCaloXTalkNeighbours(const std::string &aName, ISvcLocator *aSL)
+    : base_class(aName, aSL)
+{
+  declareProperty("outputFileName", m_outputFileName, "Name of the output file");
+}
+
+CreateFCCeeCaloXTalkNeighbours::~CreateFCCeeCaloXTalkNeighbours() {}
+
+StatusCode CreateFCCeeCaloXTalkNeighbours::initialize()
+{
+  // Initialize necessary Gaudi components
+  if (Service::initialize().isFailure())
+  {
+    error() << "Unable to initialize Service()" << endmsg;
+    return StatusCode::FAILURE;
+  }
+  m_geoSvc = service("GeoSvc");
+  if (!m_geoSvc)
+  {
+    error() << "Unable to locate Geometry Service. "
+            << "Make sure you have GeoSvc and SimSvc in the right order in the configuration." << endmsg;
+    return StatusCode::FAILURE;
+  }
+  std::unordered_map<uint64_t, std::vector<std::pair<uint64_t,double>>> map;
+#if DEBUGCELLINFO ==1
+  std::unordered_map<uint64_t, std::vector<int>> CellInfo;
+#endif
+
+  for (uint iSys = 0; iSys < m_readoutNamesSegmented.size(); iSys++)
+  {
+    // Check if readout exists
+    info() << "Readout: " << m_readoutNamesSegmented[iSys] << endmsg;
+    if (m_geoSvc->getDetector()->readouts().find(m_readoutNamesSegmented[iSys]) == m_geoSvc->getDetector()->readouts().end())
+    {
+      error() << "Readout <<" << m_readoutNamesSegmented[iSys] << ">> does not exist." << endmsg;
+      return StatusCode::FAILURE;
+    }
+
+    // get theta-based segmentation
+    dd4hep::DDSegmentation::Segmentation *aSegmentation = m_geoSvc->getDetector()->readout(m_readoutNamesSegmented[iSys]).segmentation().segmentation();
+    if (aSegmentation == nullptr)
+    {
+      error() << "Segmentation does not exist." << endmsg;
+      return StatusCode::FAILURE;
+    }
+
+    std::string segmentationType = aSegmentation->type();
+    info() << "Segmentation type : " << segmentationType << endmsg;
+
+    dd4hep::DDSegmentation::GridTheta_k4geo *segmentation = nullptr;
+    dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo *moduleThetaSegmentation = nullptr;
+    if (segmentationType == "FCCSWGridModuleThetaMerged_k4geo")
+    {
+      segmentation = dynamic_cast<dd4hep::DDSegmentation::GridTheta_k4geo *>(aSegmentation);
+      moduleThetaSegmentation = dynamic_cast<dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo *>(aSegmentation);
+    }
+    else
+    {
+      error() << "Segmentation type not handled." << endmsg;
+      return StatusCode::FAILURE;
+    }
+
+    if (segmentation == nullptr || moduleThetaSegmentation == nullptr)
+    {
+      error() << "Unable to cast segmentation pointer!!!!" << endmsg;
+      return StatusCode::FAILURE;
+    }
+
+    info() << "Segmentation: size in Theta " << segmentation->gridSizeTheta() << endmsg;
+    info() << "Segmentation: offset in Theta " << segmentation->offsetTheta() << endmsg;
+    if (segmentationType == "FCCSWGridModuleThetaMerged_k4geo")
+    {
+      info() << "Segmentation: bins in Module " << moduleThetaSegmentation->nModules() << endmsg;
+    }
+
+    // retrieve decoders for the crosstalk neighbour finding
+    auto decoder = m_geoSvc->getDetector()->readout(m_readoutNamesSegmented[iSys]).idSpec().decoder();
+
+    // Loop over all cells in the calorimeter and retrieve existing cellIDs and find neihbours
+    if (segmentationType == "FCCSWGridModuleThetaMerged_k4geo")
+    {
+      std::vector<std::vector<std::pair<int, int>>> extrema_layer;
+      for (unsigned int ilayer = 0; ilayer < m_activeVolumesNumbersSegmented[iSys]; ilayer++)
+      {
+        dd4hep::DDSegmentation::CellID volumeId = 0;
+        static_cast<dd4hep::DDSegmentation::BitFieldCoder>(*decoder)[m_fieldNamesSegmented[iSys]].set(volumeId, m_fieldValuesSegmented[iSys]);
+        static_cast<dd4hep::DDSegmentation::BitFieldCoder>(*decoder)[m_activeFieldNamesSegmented[iSys]].set(volumeId, ilayer);
+        static_cast<dd4hep::DDSegmentation::BitFieldCoder>(*decoder)["theta"].set(volumeId, 0);
+        static_cast<dd4hep::DDSegmentation::BitFieldCoder>(*decoder)["module"].set(volumeId, 0);
+        auto numCells = det::utils::numberOfCells(volumeId, *moduleThetaSegmentation);
+
+	std::vector<std::pair<int, int>> extrema;
+        // extrema[0]: min layer, n layers
+        extrema.emplace_back(std::make_pair(0, m_activeVolumesNumbersSegmented[iSys] - 1));
+        // extrema[1]: min module number (0), max module number
+        extrema.emplace_back(std::make_pair(0, (numCells[0] - 1) * moduleThetaSegmentation->mergedModules(ilayer)));
+        // extrema[2]: min theta ID, n (merged) theta cells
+        extrema.emplace_back(std::make_pair(numCells[2], numCells[2] + (numCells[1] - 1) * moduleThetaSegmentation->mergedThetaCells(ilayer)));
+        debug() << "Layer: " << ilayer << endmsg;
+        debug() << "Extrema[0]: " << extrema[0].first << " , " << extrema[0].second << endmsg;
+        debug() << "Extrema[1]: " << extrema[1].first << " , " << extrema[1].second << endmsg;
+        debug() << "Extrema[2]: " << extrema[2].first << " , " << extrema[2].second << endmsg;
+        debug() << "Number of segmentation cells in (module,theta): " << numCells << endmsg;
+	extrema_layer.emplace_back(extrema);
+      }
+
+      for (unsigned int ilayer = 0; ilayer < m_activeVolumesNumbersSegmented[iSys]; ilayer++)
+      {
+        dd4hep::DDSegmentation::CellID volumeId = 0;
+        // Get VolumeID
+        static_cast<dd4hep::DDSegmentation::BitFieldCoder>(*decoder)[m_fieldNamesSegmented[iSys]].set(volumeId, m_fieldValuesSegmented[iSys]);
+        static_cast<dd4hep::DDSegmentation::BitFieldCoder>(*decoder)[m_activeFieldNamesSegmented[iSys]].set(volumeId, ilayer);
+        static_cast<dd4hep::DDSegmentation::BitFieldCoder>(*decoder)["theta"].set(volumeId, 0);
+        static_cast<dd4hep::DDSegmentation::BitFieldCoder>(*decoder)["module"].set(volumeId, 0);
+        // Loop over segmentation cells to find crosstalk neighbours in ECAL
+        for (int imodule = extrema_layer[ilayer][1].first; imodule <= extrema_layer[ilayer][1].second; imodule += moduleThetaSegmentation->mergedModules(ilayer))
+        {
+          for (int itheta = extrema_layer[ilayer][2].first; itheta <= extrema_layer[ilayer][2].second; itheta += moduleThetaSegmentation->mergedThetaCells(ilayer))
+          {
+            dd4hep::DDSegmentation::CellID cellId = volumeId;
+            decoder->set(cellId, "module", imodule);
+            decoder->set(cellId, "theta", itheta); // start from the minimum existing theta cell in this layer
+            uint64_t id = cellId;
+            map.insert(std::pair<uint64_t, std::vector<std::pair<uint64_t,double>>>(
+                id,
+                det::crosstalk::getNeighboursModuleThetaMerged(
+                    *moduleThetaSegmentation,
+                    *decoder,
+                    {m_activeFieldNamesSegmented[iSys],
+                     "module", "theta"},
+                    extrema_layer,
+                    id)));
+          }
+        }
+      }
+    }
+
+    if (msgLevel() <= MSG::DEBUG)
+    {
+      std::vector<int> counter;
+      counter.assign(40, 0);
+      for (const auto &item : map)
+      {
+        counter[item.second.size()]++;
+      }
+      for (uint iCount = 0; iCount < counter.size(); iCount++)
+      {
+        if (counter[iCount] != 0)
+        {
+          info() << counter[iCount] << " cells have " << iCount << " neighbours" << endmsg;
+        }
+      }
+    }
+    info() << "total number of cells:  " << map.size() << endmsg;
+  }
+
+  if (msgLevel() <= MSG::DEBUG)
+  {
+    std::vector<int> counter;
+    counter.assign(40, 0);
+    for (const auto &item : map)
+    {
+      counter[item.second.size()]++;
+    }
+    for (uint iCount = 0; iCount < counter.size(); iCount++)
+    {
+      if (counter[iCount] != 0)
+      {
+        debug() << counter[iCount] << " cells have " << iCount << " neighbours" << endmsg;
+      }
+    }
+  }
+  //debug() << "cells with neighbours across Calo boundaries: " << count << endmsg;
+
+  std::unique_ptr<TFile> file(TFile::Open(m_outputFileName.c_str(), "RECREATE"));
+  file->cd();
+  TTree tree("crosstalk_neighbours", "Tree with map of neighbours");
+  uint64_t saveCellId;
+  std::vector<uint64_t> saveNeighbours;
+  std::vector<double> saveCrosstalks;
+  tree.Branch("cellId", &saveCellId, "cellId/l");
+  tree.Branch("list_crosstalk_neighbours", &saveNeighbours);
+  tree.Branch("list_crosstalks", &saveCrosstalks);
+
+  // Debug: save cell position
+#if DEBUGCELLINFO==1
+  std::vector<int> saveCellInfo;
+  tree.Branch("CellInfo", &saveCellInfo);
+#endif
+
+  int count_map=0;
+  for (const auto &item : map)
+  {
+    saveCellId = item.first;
+    std::vector<std::pair<uint64_t, double>> temp_pair=item.second;
+    std::vector<uint64_t> temp_neighbours;
+    std::vector<double> temp_crosstalks;
+    for (const auto &this_temp : temp_pair)
+    {
+      temp_neighbours.push_back(this_temp.first);
+      temp_crosstalks.push_back(this_temp.second);
+    }
+    saveNeighbours = temp_neighbours;
+    saveCrosstalks = temp_crosstalks;
+    // Debug: save cell position
+#if DEBUGCELLINFO==1
+    for (uint iSys = 0; iSys < m_readoutNamesSegmented.size(); iSys++)
+    {
+      dd4hep::DDSegmentation::Segmentation *aSegmentation = m_geoSvc->getDetector()->readout(m_readoutNamesSegmented[iSys]).segmentation().segmentation();
+      std::string segmentationType = aSegmentation->type();
+      dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo *moduleThetaSegmentation = nullptr;
+      auto decoder = m_geoSvc->getDetector()->readout(m_readoutNamesSegmented[iSys]).idSpec().decoder();
+      if (segmentationType == "FCCSWGridModuleThetaMerged_k4geo")
+      {
+        moduleThetaSegmentation = dynamic_cast<dd4hep::DDSegmentation::FCCSWGridModuleThetaMerged_k4geo *>(aSegmentation);
+      }
+      saveCellInfo=det::crosstalk::getCellIndices(
+          *moduleThetaSegmentation,
+          *decoder,
+          {m_activeFieldNamesSegmented[iSys],"module", "theta"},
+          saveCellId);
+    }
+#endif
+    tree.Fill();
+    count_map++;
+    if(!count_map%1000) std::cout<<"Number of cells: "<<count_map<<std::endl;
+  }
+  file->Write();
+  file->Close();
+
+  return StatusCode::SUCCESS;
+}
+
+StatusCode CreateFCCeeCaloXTalkNeighbours::finalize() { return Service::finalize(); }

RecFCCeeCalorimeter/src/components/CreateFCCeeCaloXTalkNeighbours.h

@@ -0,0 +1,86 @@
+#ifndef RECFCCEECALORIMETER_CREATEFCCEECALOXTALKNEIGHBOURS_H
+#define RECFCCEECALORIMETER_CREATEFCCEECALOXTALKNEIGHBOURS_H
+
+// Gaudi
+#include "GaudiKernel/Service.h"
+
+// k4FWCore
+#include "k4Interface/ICaloCreateMap.h"
+#include "k4Interface/IGeoSvc.h"
+#include "k4FWCore/DataHandle.h"
+#include "k4Interface/ICellPositionsTool.h"
+
+// k4geo
+#include "detectorCommon/DetUtils_k4geo.h"
+#include "detectorCommon/xtalk_neighbors_moduleThetaMergedSegmentation.h"
+
+// DD4hep
+#include "DD4hep/Readout.h"
+#include "DD4hep/Volumes.h"
+#include "DDSegmentation/Segmentation.h"
+
+// ROOT
+#include "TGeoManager.h"
+
+class IGeoSvc;
+
+/** @class CreateFCCeeCaloXTalkNeighbours
+ *
+ *  Service building a map of crosstalk neighbours for all existing cells in the geometry.
+ *  Only applicable to the ALLEGRO ECAL barrel with the theta-merged segmentation
+ *
+ *  @author Zhibo Wu
+ */
+
+class CreateFCCeeCaloXTalkNeighbours : public extends1<Service, ICaloCreateMap> {
+public:
+  /// Standard constructor
+  explicit CreateFCCeeCaloXTalkNeighbours(const std::string& aName, ISvcLocator* aSL);
+  /// Standard destructor
+  virtual ~CreateFCCeeCaloXTalkNeighbours();
+  /**  Initialize the map creator service.
+   *   @return status code
+   */
+  virtual StatusCode initialize() final;
+  /**  Finalize the map creator service.
+   *   @return status code
+   */
+  virtual StatusCode finalize() final;
+
+private:
+  /// Pointer to the geometry service
+  SmartIF<IGeoSvc> m_geoSvc;
+
+  /// Names of the detector readout for the volumes
+  Gaudi::Property<std::vector<std::string>> m_readoutNamesSegmented{this, "readoutNames", {"ECalBarrelModuleThetaMerged", "BarHCal_Readout_phitheta"}};
+  /// Name of the fields describing the segmented volume
+  Gaudi::Property<std::vector<std::string>> m_fieldNamesSegmented{this, "systemNames", {"system", "system"}};
+  /// Values of the fields describing the segmented volume
+  Gaudi::Property<std::vector<int>> m_fieldValuesSegmented{this, "systemValues", {4, 8}};
+  /// Names of the active volume in geometry along radial axis (e.g. layer), the others are "module" or "phi", "theta"
+  Gaudi::Property<std::vector<std::string>> m_activeFieldNamesSegmented{this, "activeFieldNames", {"layer", "layer"}};
+  /// Number of layers in the segmented volumes
+  Gaudi::Property<std::vector<unsigned int>> m_activeVolumesNumbersSegmented{this, "activeVolumesNumbers", {12, 13}};
+  // Theta ranges of layers in the segmented volumes
+  Gaudi::Property<std::vector<std::vector<double>>> m_activeVolumesTheta{this, "activeVolumesTheta"};
+
+  // System ID of ECAL and HCAL barrels
+  Gaudi::Property<uint> m_ecalBarrelSysId{this, "ecalBarrelSysId", 4};
+  Gaudi::Property<uint> m_hcalBarrelSysId{this, "hcalBarrelSysId", 8};
+
+  // For combination of barrels: flag if ECal and HCal barrels should be merged
+  Gaudi::Property<bool> m_connectBarrels{this, "connectBarrels", true};
+  // For combination of barrels: size of HCal cell along z axis
+  Gaudi::Property<double> m_hCalZSize{this, "hCalZsize", 18};
+  // For combination of barrels: offset of HCal detector in z (lower edge)
+  Gaudi::Property<double> m_hCalZOffset{this, "hCalZoffset", -4590};
+  // For combination of barrels: HCal inner radius for calculation of eta from z ???
+  Gaudi::Property<double> m_hCalRinner{this, "hCalRinner", 2850};
+  // For combination of barrels: offset of HCal modules in phi (lower edge)
+  Gaudi::Property<double> m_hCalPhiOffset{this, "hCalPhiOffset"};
+
+  /// Name of output file
+  std::string m_outputFileName;
+};
+
+#endif /* RECFCCEECALORIMETER_CREATEFCCEECALOXTALKNEIGHBOURS_H */

RecFCCeeCalorimeter/tests/options/runCaloXTalkNeighbours.py

@@ -0,0 +1,37 @@
+from Configurables import GeoSvc
+from Configurables import ApplicationMgr
+from Configurables import CreateFCCeeCaloXTalkNeighbours
+import os
+from Gaudi.Configuration import INFO, DEBUG
+
+# Detector geometry
+geoservice = GeoSvc("GeoSvc")
+# if K4GEO is empty, this should use relative path to working directory
+path_to_detector = os.environ.get("K4GEO", "")
+print(path_to_detector)
+detectors_to_use = [
+    'FCCee/ALLEGRO/compact/ALLEGRO_o1_v02/ALLEGRO_o1_v02.xml'
+]
+
+# prefix all xmls with path_to_detector
+geoservice.detectors = [os.path.join(path_to_detector, _det) for _det in detectors_to_use]
+geoservice.OutputLevel = INFO
+
+# create the crosstalk neighbour file for ECAL barrel cells
+neighbours = CreateFCCeeCaloXTalkNeighbours("xtalk_neighbours",
+                                       outputFileName="xtalk_neighbours_map_ecalB_thetamodulemerged.root",
+                                       readoutNames=["ECalBarrelModuleThetaMerged"],
+                                       systemNames=["system"],
+                                       systemValues=[4],
+                                       activeFieldNames=["layer"],
+                                       activeVolumesNumbers=[12],
+                                       OutputLevel=DEBUG)
+
+# ApplicationMgr
+ApplicationMgr(TopAlg=[],
+               EvtSel='NONE',
+               EvtMax=1,
+               # order is important, as GeoSvc is needed by G4SimSvc
+               ExtSvc=[geoservice, neighbours],
+               OutputLevel=INFO
+               )