Implementation of crosstalk for ALLEGRO ECAL Barrel (#82)

* 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

* Add reading crosstalk neighbours

* Add the crosstalk effect to calo cells

* Remove unused header

* Implementation of crosstalk

* ÂMake HitCellsMap a local variable

* Add crosstalk computation test in CMakeLists.txt

* Update the test of cross-talk map for ALLEGRO v03

* Update test macros for crosstalk calculation

---------

Co-authored-by: Zhibo Wu <[email protected]>
Co-authored-by: Zhibo Wu <[email protected]>

RecCalorimeter/src/components/CreateCaloCells.cpp

@@ -21,6 +21,7 @@ GaudiAlgorithm(name, svcLoc), m_geoSvc("GeoSvc", name) {
   declareProperty("hits", m_hits, "Hits from which to create cells (input)");
   declareProperty("cells", m_cells, "The created calorimeter cells (output)");
 
+  declareProperty("crosstalksTool", m_crosstalksTool, "Handle for the cell crosstalk tool");
   declareProperty("calibTool", m_calibTool, "Handle for tool to calibrate Geant4 energy to EM scale tool");
   declareProperty("noiseTool", m_noiseTool, "Handle for the calorimeter cells noise tool");
   declareProperty("geometryTool", m_geoTool, "Handle for the geometry tool");
@@ -37,6 +38,11 @@ StatusCode CreateCaloCells::initialize() {
   info() << "add position information to the cell : " << m_addPosition << endmsg;
 
   // Initialization of tools
+  // Cell crosstalk tool
+  if (!m_crosstalksTool.retrieve()) {
+    error() << "Unable to retrieve the cell crosstalk tool!!!" << endmsg;
+    return StatusCode::FAILURE;
+  }
   // Calibrate Geant4 energy to EM scale tool
   if (m_doCellCalibration) {
     if (!m_calibTool.retrieve()) {
@@ -89,13 +95,41 @@ StatusCode CreateCaloCells::execute() {
     m_cellsMap.clear();
   }
 
+  std::unordered_map<uint64_t, double> HitCellsMap;
+
   // 1. Merge energy deposits into cells
   // If running with noise map already was prepared. Otherwise it is being
   // created below
   for (const auto& hit : *hits) {
     verbose() << "CellID : " << hit.getCellID() << endmsg;
+    if(m_addCrosstalk) { // Crosstalk only applies to energy deposits caused by EM shower hits. Therefore, cell noise needs to be excluded.
+      HitCellsMap[hit.getCellID()] += hit.getEnergy();
+    }
     m_cellsMap[hit.getCellID()] += hit.getEnergy();
   }
+  if(m_addCrosstalk) {
+    m_CrosstalkCellsMap.clear();
+    // loop over cells that get actual EM shower hits
+    for (const auto& this_cell : HitCellsMap) {
+      uint64_t this_cellId=this_cell.first;
+      auto vec_neighbours = m_crosstalksTool->getNeighbours(this_cellId); // a vector of neighbour IDs
+      auto vec_crosstalks = m_crosstalksTool->getCrosstalks(this_cellId); // a vector of crosstalk coefficients
+      // loop over crosstalk neighbrous of the cell under study
+      for (unsigned int i_cell=0; i_cell<vec_neighbours.size(); i_cell++) {
+	// signal transfer = energy deposit brought by EM shower hits * crosstalk coefficient
+        double signal_transfer = this_cell.second * vec_crosstalks[i_cell];
+	// for the cell under study, record the signal transfer that will be subtracted from its final cell energy
+        m_CrosstalkCellsMap[this_cellId] -= signal_transfer;
+	// for the crosstalk neighbour, record the signal transfer that will be added to its final cell energy
+        m_CrosstalkCellsMap[vec_neighbours[i_cell]] += signal_transfer;
+      }
+    }
+
+    for (const auto& this_cell : m_CrosstalkCellsMap) {
+      m_cellsMap[this_cell.first] += this_cell.second;
+    }
+
+  }
   debug() << "Number of calorimeter cells after merging of hits: " << m_cellsMap.size() << endmsg;
 
   // 2. Calibrate simulation energy to EM scale

RecCalorimeter/src/components/CreateCaloCells.h

@@ -7,6 +7,7 @@
 #include "k4Interface/ICalibrateCaloHitsTool.h"
 #include "k4Interface/ICalorimeterTool.h"
 #include "k4Interface/INoiseCaloCellsTool.h"
+#include "k4Interface/ICaloReadCrosstalkMap.h"
 
 // Gaudi
 #include "GaudiAlg/GaudiAlgorithm.h"
@@ -58,13 +59,18 @@ class CreateCaloCells : public GaudiAlgorithm {
   StatusCode finalize();
 
 private:
+
+  /// Handle for the calorimeter cells crosstalk tool
+  ToolHandle<ICaloReadCrosstalkMap> m_crosstalksTool{"ReadCaloCrosstalkMap", this};
   /// Handle for tool to calibrate Geant4 energy to EM scale tool
   ToolHandle<ICalibrateCaloHitsTool> m_calibTool{"CalibrateCaloHitsTool", this};
   /// Handle for the calorimeter cells noise tool
   ToolHandle<INoiseCaloCellsTool> m_noiseTool{"NoiseCaloCellsFlatTool", this};
   /// Handle for the geometry tool
   ToolHandle<ICalorimeterTool> m_geoTool{"TubeLayerPhiEtaCaloTool", this};
 
+  /// Add crosstalk to cells?
+  Gaudi::Property<bool> m_addCrosstalk{this, "addCrosstalk", true, "Add crosstalk effect?"};
   /// Calibrate to EM scale?
   Gaudi::Property<bool> m_doCellCalibration{this, "doCellCalibration", true, "Calibrate to EM scale?"};
   /// Add noise to cells?
@@ -113,8 +119,9 @@ class CreateCaloCells : public GaudiAlgorithm {
   /// Pointer to the geometry service
   ServiceHandle<IGeoSvc> m_geoSvc;
   dd4hep::VolumeManager m_volman;
-  /// Map of cell IDs (corresponding to DD4hep IDs) and energy
+  /// Maps of cell IDs (corresponding to DD4hep IDs) on (1) final energies to be used for clustering, (2) transfer of signals due to crosstalk
   std::unordered_map<uint64_t, double> m_cellsMap;
+  std::unordered_map<uint64_t, double> m_CrosstalkCellsMap;
 };
 
 #endif /* RECCALORIMETER_CREATECALOCELLS_H */

RecCalorimeter/src/components/ReadCaloCrosstalkMap.cpp

@@ -0,0 +1,52 @@
+#include "ReadCaloCrosstalkMap.h"
+
+#include "TFile.h"
+#include "TTree.h"
+#include "TBranch.h"
+
+DECLARE_COMPONENT(ReadCaloCrosstalkMap)
+
+ReadCaloCrosstalkMap::ReadCaloCrosstalkMap(const std::string& type, const std::string& name,
+                                                 const IInterface* parent)
+    : GaudiTool(type, name, parent) {
+  declareInterface<ICaloReadCrosstalkMap>(this);
+}
+
+StatusCode ReadCaloCrosstalkMap::initialize() {
+  StatusCode sc = GaudiTool::initialize();
+  info() <<"Loading crosstalk map..." << endmsg;
+  if (sc.isFailure()) return sc;
+  std::unique_ptr<TFile> file(TFile::Open(m_fileName.value().c_str(),"READ"));
+  TTree* tree = nullptr;
+  file->GetObject("crosstalk_neighbours",tree);
+  ULong64_t read_cellId;
+  std::vector<uint64_t>  *read_neighbours=0;
+  std::vector<double> *read_crosstalks=0;
+
+  tree->SetBranchAddress("cellId",&read_cellId);
+  tree->SetBranchAddress("list_crosstalk_neighbours", &read_neighbours);
+  tree->SetBranchAddress("list_crosstalks", &read_crosstalks);
+  for (uint i = 0; i < tree->GetEntries(); i++) {
+      tree->GetEntry(i);
+      m_mapNeighbours.insert(std::pair<uint64_t, std::vector<uint64_t>>(read_cellId, *read_neighbours));
+      m_mapCrosstalks.insert(std::pair<uint64_t, std::vector<double>>(read_cellId, *read_crosstalks));
+  }
+
+  info() <<"Crosstalk input: " << m_fileName.value().c_str() << endmsg;
+  info() << "Total number of cells = " << tree->GetEntries() << ", Size of crosstalk neighbours = " << m_mapNeighbours.size() << ", Size of coefficients = " << m_mapCrosstalks.size() << endmsg;
+  delete tree;
+  delete read_neighbours;
+  delete read_crosstalks;
+  file->Close();
+  return sc;
+}
+
+StatusCode ReadCaloCrosstalkMap::finalize() { return GaudiTool::finalize(); }
+
+std::vector<uint64_t>& ReadCaloCrosstalkMap::getNeighbours(uint64_t aCellId) {
+  return m_mapNeighbours[aCellId];
+}
+
+std::vector<double>& ReadCaloCrosstalkMap::getCrosstalks(uint64_t aCellId) {
+  return m_mapCrosstalks[aCellId];
+}

RecCalorimeter/src/components/ReadCaloCrosstalkMap.h

@@ -0,0 +1,49 @@
+#ifndef RECCALORIMETER_READCALOXTALKMAP_H
+#define RECCALORIMETER_READCALOXTALKMAP_H
+
+// from Gaudi
+#include "GaudiAlg/GaudiTool.h"
+
+// k4FWCore
+#include "k4Interface/ICaloReadCrosstalkMap.h"
+
+class IGeoSvc;
+
+/** @class ReadCaloCrosstalkMap Reconstruction/RecCalorimeter/src/components/ReadCaloCrosstalkMap.h
+ *TopoCaloNeighbours.h
+ *
+ *  Tool that reads a ROOT file containing the TTree with branches "cellId", "list_crosstalk_neighbours" and "list_crosstalks".
+ *  This tools reads the tree, creates two maps, and allows a lookup of all crosstalk neighbours as well as the corresponding crosstalk coefficients for a given cell.
+ *
+ *  @author Zhibo Wu
+ */
+
+class ReadCaloCrosstalkMap : public GaudiTool, virtual public ICaloReadCrosstalkMap {
+public:
+  ReadCaloCrosstalkMap(const std::string& type, const std::string& name, const IInterface* parent);
+  virtual ~ReadCaloCrosstalkMap() = default;
+
+  virtual StatusCode initialize() final;
+  virtual StatusCode finalize() final;
+
+  /** Function to be called for the crosstalk neighbours of a cell.
+   *   @param[in] aCellId, cellid of the cell of interest.
+   *   @return vector of cellIDs, corresponding to the crosstalk neighbours.
+   */
+  virtual std::vector<uint64_t>& getNeighbours(uint64_t aCellId) final;
+
+  /** Function to be called for the crosstalk coefficients between the input cell and its neighbouring cells.
+   *   @param[in] aCellId, cellid of the cell of interest.
+   *   @return vector of crosstalk coefficients.
+   */
+  virtual std::vector<double>& getCrosstalks(uint64_t aCellId) final;
+
+private:
+  /// Name of input root file that contains the TTree with cellID->vec<list_crosstalk_neighboursCellID> and cellId->vec<list_crosstalksCellID>
+  Gaudi::Property<std::string> m_fileName{this, "fileName", "xtalk_neighbours_map_ecalB_thetamodulemerged.root", "Name of the file that contains the crosstalk map"};
+  /// Output maps to be used for the fast lookup in the creating calo-cells algorithm
+  std::unordered_map<uint64_t, std::vector<uint64_t>> m_mapNeighbours;
+  std::unordered_map<uint64_t, std::vector<double>> m_mapCrosstalks;
+};
+
+#endif /* RECCALORIMETER_READCALOXTALKMAP_H */