diff --git a/RecFCCeeCalorimeter/src/components/CreatePositionedCaloCells.cpp b/RecFCCeeCalorimeter/src/components/CreatePositionedCaloCells.cpp new file mode 100644 index 00000000..7dc89b1a --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/CreatePositionedCaloCells.cpp @@ -0,0 +1,197 @@ +#include "CreatePositionedCaloCells.h" + +// k4geo +#include "detectorCommon/DetUtils_k4geo.h" + +// k4FWCore +#include "k4Interface/IGeoSvc.h" + +// DD4hep +//#include "DD4hep/Detector.h" +//#include "DD4hep/Volumes.h" +//#include "TGeoManager.h" + +// edm4hep +#include "edm4hep/CalorimeterHit.h" + +DECLARE_COMPONENT(CreatePositionedCaloCells) + +CreatePositionedCaloCells::CreatePositionedCaloCells(const std::string& name, ISvcLocator* svcLoc) : +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("positionsTool", m_cellPositionsTool, "Handle for cell positions tool"); + declareProperty("crosstalkTool", m_crosstalkTool, "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"); +} + +StatusCode CreatePositionedCaloCells::initialize() { + StatusCode sc = GaudiAlgorithm::initialize(); + if (sc.isFailure()) return sc; + + info() << "CreatePositionedCaloCells initialized" << endmsg; + info() << "do calibration : " << m_doCellCalibration << endmsg; + info() << "add cell noise : " << m_addCellNoise << endmsg; + info() << "remove cells below threshold : " << m_filterCellNoise << endmsg; + info() << "emulate crosstalk : " << m_addCrosstalk << endmsg; + + // Initialization of tools + + // Cell position tool + if (!m_cellPositionsTool.retrieve()) { + error() << "Unable to retrieve the cell positions tool!!!" << endmsg; + return StatusCode::FAILURE; + } + // Cell crosstalk tool + if (m_addCrosstalk) { + if (!m_crosstalkTool.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()) { + error() << "Unable to retrieve the calo cells calibration tool!!!" << endmsg; + return StatusCode::FAILURE; + } + } + // Cell noise tool + if (m_addCellNoise || m_filterCellNoise) { + if (!m_noiseTool.retrieve()) { + error() << "Unable to retrieve the calo cells noise tool!!!" << endmsg; + return StatusCode::FAILURE; + } + // Geometry settings + if (!m_geoTool.retrieve()) { + error() << "Unable to retrieve the geometry tool!!!" << endmsg; + return StatusCode::FAILURE; + } + // Prepare map of all existing cells in calorimeter to add noise to all + StatusCode sc_prepareCells = m_geoTool->prepareEmptyCells(m_cellsMap); + if (sc_prepareCells.isFailure()) { + error() << "Unable to create empty cells!" << endmsg; + return StatusCode::FAILURE; + } + } + + // Copy over the CellIDEncoding string from the input collection to the output collection + auto hitsEncoding = m_hitsCellIDEncoding.get_optional(); + if (!hitsEncoding.has_value()) { + error () << "Missing cellID encoding for input collection" << endmsg; + return StatusCode::FAILURE; + } + m_cellsCellIDEncoding.put(hitsEncoding.value()); + + return StatusCode::SUCCESS; +} + +StatusCode CreatePositionedCaloCells::execute() { + // Get the input collection with Geant4 hits + const edm4hep::SimCalorimeterHitCollection* hits = m_hits.get(); + debug() << "Input Hit collection size: " << hits->size() << endmsg; + + // 0. Clear all cells + if (m_addCellNoise) { + std::for_each(m_cellsMap.begin(), m_cellsMap.end(), [](std::pair& p) { p.second = 0; }); + } else { + m_cellsMap.clear(); + } + + + // 1. Merge energy deposits into cells + // If running with noise, map was already prepared in initialize(). + // Otherwise it is being created below + for (const auto& hit : *hits) { + verbose() << "CellID : " << hit.getCellID() << endmsg; + m_cellsMap[hit.getCellID()] += hit.getEnergy(); + } + debug() << "Number of calorimeter cells after merging of hits: " << m_cellsMap.size() << endmsg; + + // 2. Emulate cross-talk (if asked) + if (m_addCrosstalk) { + // Derive the cross-talk contributions without affecting yet the nominal energy + // (one has to emulate crosstalk based on cells free from any cross-talk contributions) + m_crosstalkCellsMap.clear(); // this is a temporary map to hold energy exchange due to cross-talk, without affecting yet the nominal energy + // loop over cells with nominal energies + for (const auto& this_cell : m_cellsMap) { + uint64_t this_cellId = this_cell.first; + auto vec_neighbours = m_crosstalkTool->getNeighbours(this_cellId); // a vector of neighbour IDs + auto vec_crosstalks = m_crosstalkTool->getCrosstalks(this_cellId); // a vector of crosstalk coefficients + // loop over crosstalk neighbours of the cell under study + for (unsigned int i_cell=0; i_cellcalibrate(m_cellsMap); + } + + // 4. Add noise to all cells + if (m_addCellNoise) { + m_noiseTool->addRandomCellNoise(m_cellsMap); + } + + // 5. Filter cells + if (m_filterCellNoise) { + m_noiseTool->filterCellNoise(m_cellsMap); + } + + // 6. Copy information to CaloHitCollection + edm4hep::CalorimeterHitCollection* edmCellsCollection = new edm4hep::CalorimeterHitCollection(); + for (const auto& cell : m_cellsMap) { + if (m_addCellNoise || (!m_addCellNoise && cell.second != 0)) { + auto newCell = edmCellsCollection->create(); + newCell.setEnergy(cell.second); + uint64_t cellid = cell.first; + newCell.setCellID(cellid); + + // add cell position + auto cached_pos = m_positions_cache.find(cellid); + if(cached_pos == m_positions_cache.end()) { + // retrieve position from tool + dd4hep::Position posCell = m_cellPositionsTool->xyzPosition(cellid); + edm4hep::Vector3f edmPos; + edmPos.x = posCell.x() / dd4hep::mm; + edmPos.y = posCell.y() / dd4hep::mm; + edmPos.z = posCell.z() / dd4hep::mm; + m_positions_cache[cellid] = edmPos; + newCell.setPosition(edmPos); + } + else { + newCell.setPosition(cached_pos->second); + } + + debug() << "Cell energy (GeV) : " << newCell.getEnergy() << "\tcellID " << newCell.getCellID() << endmsg; + debug() << "Position of cell (mm) : \t" << newCell.getPosition().x/dd4hep::mm + << "\t" << newCell.getPosition().y/dd4hep::mm + << "\t" << newCell.getPosition().z/dd4hep::mm << endmsg; + } + } + + // push the CaloHitCollection to event store + m_cells.put(edmCellsCollection); + + debug() << "Output Cell collection size: " << edmCellsCollection->size() << endmsg; + + return StatusCode::SUCCESS; +} + +StatusCode CreatePositionedCaloCells::finalize() { return GaudiAlgorithm::finalize(); } diff --git a/RecFCCeeCalorimeter/src/components/CreatePositionedCaloCells.h b/RecFCCeeCalorimeter/src/components/CreatePositionedCaloCells.h new file mode 100644 index 00000000..16b771f0 --- /dev/null +++ b/RecFCCeeCalorimeter/src/components/CreatePositionedCaloCells.h @@ -0,0 +1,108 @@ +#ifndef RECFCCEECALORIMETER_CREATEPOSITIONEDCALOCELLS_H +#define RECFCCEECALORIMETER_CREATEPOSITIONEDCALOCELLS_H + +// k4FWCore +#include "k4FWCore/DataHandle.h" +#include "k4FWCore/MetaDataHandle.h" +#include "k4Interface/ICellPositionsTool.h" +#include "k4Interface/ICalibrateCaloHitsTool.h" +#include "k4Interface/ICalorimeterTool.h" +#include "k4Interface/INoiseCaloCellsTool.h" +#include "k4Interface/ICaloReadCrosstalkMap.h" + +// Gaudi +#include "GaudiAlg/GaudiAlgorithm.h" +#include "GaudiKernel/ToolHandle.h" + +// edm4hep +#include "edm4hep/CalorimeterHitCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/Constants.h" + +// DD4hep +// #include "DD4hep/Detector.h" +// #include "DD4hep/Volumes.h" +// #include "TGeoManager.h" + +class IGeoSvc; + +/** @class CreatePositionedCaloCells + * + * Algorithm for creating positioned calorimeter cells from Geant4 hits. + * Based on CreateCaloCells, adding a positioning tool so that output + * cells from the digitisation contain the correct position in 3D space. + * + * Flow of the program: + * 1/ Merge Geant4 energy deposits with same cellID + * 2/ Emulate cross-talk (if switched on) + * 3/ Calibrate to electromagnetic scale (if calibration switched on) + * 4/ Add random noise to each cell (if noise switched on) + * 5/ Filter cells and remove those with energy below threshold (if noise + + * filtering switched on) + * 6/ Add cell positions + * + * Tools called: + * - CalibrateCaloHitsTool + * - NoiseCaloCellsTool + * - CaloReadCrosstalkMap + * - CalorimeterTool (for geometry) + * - CellPositionsTool + * + * @author Giovanni Marchiori + * @date 2024-07 + * + */ + +class CreatePositionedCaloCells : public GaudiAlgorithm { + +public: + CreatePositionedCaloCells(const std::string& name, ISvcLocator* svcLoc); + + StatusCode initialize(); + + StatusCode execute(); + + StatusCode finalize(); + +private: + /// Handle for tool to get cells positions + ToolHandle m_cellPositionsTool{"CellPositionsTool", this}; + /// Handle for the calorimeter cells crosstalk tool + ToolHandle m_crosstalkTool{"ReadCaloCrosstalkMap", this}; + /// Handle for tool to calibrate Geant4 energy to EM scale tool + ToolHandle m_calibTool{"CalibrateCaloHitsTool", this}; + /// Handle for the calorimeter cells noise tool + ToolHandle m_noiseTool{"NoiseCaloCellsFlatTool", this}; + /// Handle for the geometry tool + ToolHandle m_geoTool{"TubeLayerPhiEtaCaloTool", this}; + + /// Add crosstalk to cells? + Gaudi::Property m_addCrosstalk{this, "addCrosstalk", false, "Add crosstalk effect?"}; + /// Calibrate to EM scale? + Gaudi::Property m_doCellCalibration{this, "doCellCalibration", true, "Calibrate to EM scale?"}; + /// Add noise to cells? + Gaudi::Property m_addCellNoise{this, "addCellNoise", true, "Add noise to cells?"}; + /// Save only cells with energy above threshold? + Gaudi::Property m_filterCellNoise{this, "filterCellNoise", false, + "Save only cells with energy above threshold?"}; + + /// Handle for calo hits (input collection) + DataHandle m_hits{"hits", Gaudi::DataHandle::Reader, this}; + /// Handle for the cellID encoding string of the input collection + MetaDataHandle m_hitsCellIDEncoding{m_hits, edm4hep::labels::CellIDEncoding, Gaudi::DataHandle::Reader}; + /// Handle for calo cells (output collection) + DataHandle m_cells{"cells", Gaudi::DataHandle::Writer, this}; + MetaDataHandle m_cellsCellIDEncoding{m_cells, edm4hep::labels::CellIDEncoding, Gaudi::DataHandle::Writer}; + /// Pointer to the geometry service (is this needed?) + ServiceHandle m_geoSvc; + /// dd4hep::VolumeManager m_volman; (not needed) + /// Maps of cell IDs (corresponding to DD4hep IDs) on final energies to be used for clustering + std::unordered_map m_cellsMap; + /// Maps of cell IDs (corresponding to DD4hep IDs) on transfer of signals due to crosstalk + std::unordered_map m_crosstalkCellsMap; + + /// Cache position vs cellID + std::unordered_map m_positions_cache{}; +}; + +#endif /* RECFCCEECALORIMETER_CREATEPOSITIONEDCALOCELLS_H */