diff --git a/Calorimeter/CaloEvent/CaloEvent/CaloCellClusterWeights.h b/Calorimeter/CaloEvent/CaloEvent/CaloCellClusterWeights.h index b8b0ec8688918b2649a62c32b224a146b1a4b118..956378604ca1d511aa0d93d731faf097deeb7bf6 100644 --- a/Calorimeter/CaloEvent/CaloEvent/CaloCellClusterWeights.h +++ b/Calorimeter/CaloEvent/CaloEvent/CaloCellClusterWeights.h @@ -32,8 +32,6 @@ public: typedef store_t::iterator iterator; ///< @brief Iterator type /// @} - /// @name Constructors - /// @{ /// @brief Default constructor /// /// The constructed data object provides a reserved and sized store appropriate for the total number of cells in the calorimeter. @@ -45,8 +43,8 @@ public: /// /// @param size requested store size CaloCellClusterWeights(size_t size); - /// @} - + /// @brief Copy constructor + CaloCellClusterWeights(const CaloCellClusterWeights& cellClusterWeights); /// @brief Destructor virtual ~CaloCellClusterWeights(); diff --git a/Calorimeter/CaloEvent/src/CaloCellClusterWeights.cxx b/Calorimeter/CaloEvent/src/CaloCellClusterWeights.cxx index 4efe29d829c934c75ae0ada0c873adc514750155..8ea3633c3efb8015809b691da3194cb42b1f877e 100644 --- a/Calorimeter/CaloEvent/src/CaloCellClusterWeights.cxx +++ b/Calorimeter/CaloEvent/src/CaloCellClusterWeights.cxx @@ -21,6 +21,10 @@ CaloCellClusterWeights::CaloCellClusterWeights() : CaloCellClusterWeights( CELLCLUSTERLOOKUP ) { } +CaloCellClusterWeights::CaloCellClusterWeights(const CaloCellClusterWeights& cellClusterWeights) + : _hashTable(cellClusterWeights._hashTable) +{ } + CaloCellClusterWeights::~CaloCellClusterWeights() { } diff --git a/Calorimeter/CaloExample/CaloRecEx/share/CaloRecOutputItemList_jobOptions.py b/Calorimeter/CaloExample/CaloRecEx/share/CaloRecOutputItemList_jobOptions.py index 6d7041f7d903aa99217a21eb90999f15029e6d73..ac1be5c7246b363a033773abcddadbefeea7ea0a 100644 --- a/Calorimeter/CaloExample/CaloRecEx/share/CaloRecOutputItemList_jobOptions.py +++ b/Calorimeter/CaloExample/CaloRecEx/share/CaloRecOutputItemList_jobOptions.py @@ -7,10 +7,6 @@ CaloESDList = [] CaloESDList += [ "CaloCellContainer#AllCalo" ] -# compactify calo cell -from CaloTools.CaloToolsConf import CaloCompactCellTool -svcMgr.ToolSvc += CaloCompactCellTool() - # add explicitly E4', MBTS cells and trigger output to ESD CaloESDList += [ "TileCellContainer#E4prContainer" ] CaloESDList += [ "TileCellContainer#MBTSContainer" ] @@ -21,6 +17,11 @@ CaloClusterItemList=[] CaloClusterKeys=[] CaloClusterKeys+=["CaloCalTopoClusters"] +if jobproperties.CaloRecFlags.doCaloTopoTower.get_Value(): + CaloClusterKeys+=["CaloCalTopoTowers"] +if jobproperties.CaloRecFlags.doCaloTopoSignal.get_Value(): + CaloClusterKeys+=["CaloCalTopoSignals"] +##CaloClusterKeys+=["CaloCalFwdTopoTowers"] CaloClusterKeys+=["CombinedCluster"] #CaloClusterKeys+=["EMTopoCluster430"] CaloClusterKeys+=["EMTopoSW35"] @@ -84,9 +85,9 @@ if jobproperties.Beam.beamType() == 'cosmics' or jobproperties.Beam.beamType() = #List of AOD moments: (copied from CaloClusterTopoGetter) -AODMoments=[#"LATERAL" +AODMoments=["SECOND_R" + #,"LATERAL" #,"LONGITUDINAL" - "SECOND_R" ,"SECOND_LAMBDA" ,"CENTER_MAG" ,"CENTER_LAMBDA" @@ -107,36 +108,55 @@ AODMoments=[#"LATERAL" ,"EM_PROBABILITY" #,"PTD" ,"BadChannelList" - ,#"LATERAL" + #,"MASS" ] + +if jobproperties.CaloRecFlags.doExtendedClusterMoments.get_Value(): + AODMoments += ["LATERAL" + ,"ENG_BAD_HV_CELLS" + ,"N_BAD_HV_CELLS" + ,"SIGNIFICANCE" + ,"CELL_SIGNIFICANCE" + ,"CELL_SIG_SAMPLING" + ,"PTD" + ,"MASS" + ] try: from Digitization.DigitizationFlags import digitizationFlags if digitizationFlags.doDigiTruth(): - AODMoments+=["SECOND_R_DigiHSTruth" - ,"SECOND_LAMBDA_DigiHSTruth" - ,"CENTER_MAG_DigiHSTruth" - ,"CENTER_LAMBDA_DigiHSTruth" - ,"FIRST_ENG_DENS_DigiHSTruth" - ,"ENG_FRAC_MAX_DigiHSTruth" - ,"ISOLATION_DigiHSTruth" - ,"ENG_BAD_CELLS_DigiHSTruth" - ,"N_BAD_CELLS_DigiHSTruth" - ,"BADLARQ_FRAC_DigiHSTruth" - #,"ENG_BAD_HV_CELLS_Truth" - #,"N_BAD_HV_CELLS_Truth" - ,"ENG_POS_DigiHSTruth" - #,"SIGNIFICANCE_Truth" - #,"CELL_SIGNIFICANCE_Truth" - #,"CELL_SIG_SAMPLING_Truth" - ,"AVG_LAR_Q_DigiHSTruth" - ,"AVG_TILE_Q_DigiHSTruth" - ,"EM_PROBABILITY_DigiHSTruth" - #,"PTD_Truth" - ,"ENERGY_DigiHSTruth" - ,"ETA_DigiHSTruth" - ,"PHI_DigiHSTruth" - ] + AODMoments+=["SECOND_R_DigiHSTruth" + ,"SECOND_LAMBDA_DigiHSTruth" + ,"CENTER_MAG_DigiHSTruth" + ,"CENTER_LAMBDA_DigiHSTruth" + ,"FIRST_ENG_DENS_DigiHSTruth" + ,"ENG_FRAC_MAX_DigiHSTruth" + ,"ISOLATION_DigiHSTruth" + ,"ENG_BAD_CELLS_DigiHSTruth" + ,"N_BAD_CELLS_DigiHSTruth" + ,"BADLARQ_FRAC_DigiHSTruth" + #,"ENG_BAD_HV_CELLS_Truth" + #,"N_BAD_HV_CELLS_Truth" + ,"ENG_POS_DigiHSTruth" + #,"SIGNIFICANCE_Truth" + #,"CELL_SIGNIFICANCE_Truth" + #,"CELL_SIG_SAMPLING_Truth" + ,"AVG_LAR_Q_DigiHSTruth" + ,"AVG_TILE_Q_DigiHSTruth" + ,"EM_PROBABILITY_DigiHSTruth" + #,"PTD_Truth" + ,"ENERGY_DigiHSTruth" + ,"ETA_DigiHSTruth" + ,"PHI_DigiHSTruth" + ] + if jobproperties.CaloRecFlags.doExtendedClusterMoments.get_Value(): + AODMoments+=["ENG_BAD_HV_CELLS_Truth" + ,"N_BAD_HV_CELLS_Truth" + ,"SIGNIFICANCE_Truth" + ,"CELL_SIGNIFICANCE_Truth" + ,"CELL_SIG_SAMPLING_Truth" + ,"PTD_Truth" + ] except: log = logging.getLogger('CaloRecOutputItemList') log.info('Unable to import DigitizationFlags in CaloRecOutputItemList_jobOptions. Expected in AthenaP1') @@ -153,6 +173,11 @@ CaloClusterKeys=[] CaloClusterKeys+=["CaloCalTopoClusters"] +if jobproperties.CaloRecFlags.doCaloTopoTower.get_Value(): + CaloClusterKeys+=["CaloCalTopoTowers"] +if jobproperties.CaloRecFlags.doCaloTopoSignal.get_Value(): + CaloClusterKeys+=["CaloCalTopoSignals"] +##CaloClusterKeys+=["CaloCalFwdTopoTowers"] CaloClusterKeys+=["CombinedCluster"] #CaloClusterKeys+=["EMTopoCluster430"] CaloClusterKeys+=["EMTopoSW35"] diff --git a/Calorimeter/CaloRec/CaloRec/CaloTopoClusterFromTowerCalibrator.h b/Calorimeter/CaloRec/CaloRec/CaloTopoClusterFromTowerCalibrator.h index 52c5a21cb3ef285d77b32495b26ed130d68d722b..dd6b13bbedfea55d726448bea8c313e4da91f955 100644 --- a/Calorimeter/CaloRec/CaloRec/CaloTopoClusterFromTowerCalibrator.h +++ b/Calorimeter/CaloRec/CaloRec/CaloTopoClusterFromTowerCalibrator.h @@ -3,6 +3,8 @@ #ifndef CALOREC_CALOTOPOCLUSTERFROMTOWERCALIBRATOR_H #define CALOREC_CALOTOPOCLUSTERFROMTOWERCALIBRATOR_H +#include "StoreGate/ReadHandleKey.h" + #include "AthenaBaseComps/AthAlgTool.h" #include "CaloRec/CaloClusterCollectionProcessor.h" @@ -18,6 +20,9 @@ public: /// @brief Tool constructor CaloTopoClusterFromTowerCalibrator(const std::string& type,const std::string& name,const IInterface* pParent); + /// @brief Tool initialization + StatusCode initialize(); + /// @brief Tool execution StatusCode execute(xAOD::CaloClusterContainer* pClusCont); @@ -25,16 +30,20 @@ private: /// @name Tool properties /// @{ - std::string m_cellClusterWeightsKey; ///> @brief Key for cell-cluster weight look up object - bool m_orderByPt; ///> @brief Turn on pT ordering + SG::ReadHandleKey m_cellClusterWeightKey; ///< Key for handle to cell-cluster weight look up object + bool m_orderByPt; ///< Turn on pT ordering if @c true /// @} - - /// @brief Default container key - static std::string m_defaultKey; }; /// @class CaloTopoClusterFromTowerCalibrator /// -/// @brief A cluster builder tool calibrated topo-clusters formed from calorimeter towers to the LCW scale. +/// @brief A cluster builder tool to calibrate topo-clusters formed from (EM) calorimeter towers to the LCW scale. /// +/// This module applies LCW weights to cells contributing to towers represented by @c xAOD::CaloCluster objects. +/// The overall energy contribution of a given cell contributing to a given tower is then @f$ w_{\rm geo} \times w_{\rm LCW} \times E_{\rm cell} @f$, +/// where @f$ w_{\rm geo} @f$ is the geometrical weight, @f$ w_{\rm LCW} @f$ is +/// the calibration weight the cell received from the LCW calibration in the context of the @c xAOD::CaloCluster objects it contributes to (at most two), +/// and @f$ E_{\rm cell} @f$ is the cell energy on EM scale. More details on the weights are given on +/// this page. +/// /// @author Peter Loch #endif diff --git a/Calorimeter/CaloRec/CaloRec/CaloTopoClusterFromTowerHelpers.h b/Calorimeter/CaloRec/CaloRec/CaloTopoClusterFromTowerHelpers.h index bf56e0146012d557a9dacfc7078d4dd5d3f3df8a..797de7a586b8ec3570d9c301d7204e6b3b413193 100644 --- a/Calorimeter/CaloRec/CaloRec/CaloTopoClusterFromTowerHelpers.h +++ b/Calorimeter/CaloRec/CaloRec/CaloTopoClusterFromTowerHelpers.h @@ -8,6 +8,14 @@ #include "xAODCaloEvent/CaloCluster.h" #include +#include +#include + +#define _SNAMELU( NAME ) { CaloSampling::NAME , std::string( #NAME ) } +#define _SIDLU( ID ) { std::string( #ID ) , CaloSampling::ID } + +#define _MNAMELU( NAME ) { xAOD::CaloCluster::NAME, std::string( #NAME ) } +#define _MIDLU( ID ) { std::string( #ID ) , xAOD::CaloCluster::ID } class CaloCell; @@ -147,5 +155,117 @@ namespace CaloRec { bool calculateKine(xAOD::CaloCluster* pClus,bool onlyKine=false); /// @} } // Helpers + + ///@brief Stores + namespace Lookup { + ///@name Sampling names and identifier maps + ///@{ + static const std::map samplingNames = { + // EM + _SNAMELU( PreSamplerB ), _SNAMELU( EMB1 ), _SNAMELU( EMB2 ), _SNAMELU( EMB3 ), + _SNAMELU( PreSamplerE ), _SNAMELU( EME1 ), _SNAMELU( EME2 ), _SNAMELU( EME3 ), + // HAD + _SNAMELU( HEC0 ), _SNAMELU( HEC1 ), _SNAMELU( HEC2 ), _SNAMELU( HEC3 ), + _SNAMELU( TileBar0 ), _SNAMELU( TileBar1 ), _SNAMELU( TileBar2 ), + // Tile + _SNAMELU( TileBar0 ), _SNAMELU( TileBar1 ), _SNAMELU( TileBar2 ), + _SNAMELU( TileGap1 ), _SNAMELU( TileGap2 ), _SNAMELU( TileGap3 ), + _SNAMELU( TileExt0 ), _SNAMELU( TileExt1 ), _SNAMELU( TileExt2 ), + // FCal + _SNAMELU( FCAL0 ), _SNAMELU( FCAL1 ), _SNAMELU( FCAL2 ), + // Mini-FCal + _SNAMELU( MINIFCAL0 ), _SNAMELU( MINIFCAL1 ), _SNAMELU( MINIFCAL2 ), _SNAMELU( MINIFCAL3 ), + // unknown + _SNAMELU( Unknown ) + }; ///< Name lookup by sampling identifier (data) + static const std::map samplingIds = { + // EM + _SIDLU( PreSamplerB ), _SIDLU( EMB1 ), _SIDLU( EMB2 ), _SIDLU( EMB3 ), + _SIDLU( PreSamplerE ), _SIDLU( EME1 ), _SIDLU( EME2 ), _SIDLU( EME3 ), + // HAD + _SIDLU( HEC0 ), _SIDLU( HEC1 ), _SIDLU( HEC2 ), _SIDLU( HEC3 ), + _SIDLU( TileBar0 ), _SIDLU( TileBar1 ), _SIDLU( TileBar2 ), + // Tile + _SIDLU( TileBar0 ), _SIDLU( TileBar1 ), _SIDLU( TileBar2 ), + _SIDLU( TileGap1 ), _SIDLU( TileGap2 ), _SIDLU( TileGap3 ), + _SIDLU( TileExt0 ), _SIDLU( TileExt1 ), _SIDLU( TileExt2 ), + // FCal + _SIDLU( FCAL0 ), _SIDLU( FCAL1 ), _SIDLU( FCAL2 ), + // Mini-FCal + _SIDLU( MINIFCAL0 ), _SIDLU( MINIFCAL1 ), _SIDLU( MINIFCAL2 ), _SIDLU( MINIFCAL3 ), + // unknown + _SIDLU( Unknown ) + }; ///< Identifier lookup by sampling name (data) + static const std::string& getSamplingName(CaloSampling::CaloSample sid); ///< Lookup sampling name by identifier (function) + static CaloSampling::CaloSample getSamplingId(const std::string& sname); ///< Lookup sampling identifier by name (function) + ///@} + ///@name Cluster moment names and identifier maps + ///@{ + static const std::map clusterMomentNames = { + _MNAMELU( FIRST_PHI ),_MNAMELU( FIRST_ETA ),_MNAMELU( SECOND_R ),_MNAMELU( SECOND_LAMBDA ),_MNAMELU( DELTA_PHI ),_MNAMELU( DELTA_THETA ),_MNAMELU( DELTA_ALPHA ), + _MNAMELU( CENTER_X ),_MNAMELU( CENTER_Y ),_MNAMELU( CENTER_Z ),_MNAMELU( CENTER_MAG ),_MNAMELU( CENTER_LAMBDA ),_MNAMELU( LATERAL ),_MNAMELU( LONGITUDINAL ), + _MNAMELU( ENG_FRAC_EM ),_MNAMELU( ENG_FRAC_MAX ),_MNAMELU( ENG_FRAC_CORE ),_MNAMELU( FIRST_ENG_DENS ),_MNAMELU( SECOND_ENG_DENS ),_MNAMELU( ENG_POS ), + _MNAMELU( ISOLATION ),_MNAMELU( ENG_BAD_CELLS ),_MNAMELU( N_BAD_CELLS ),_MNAMELU( N_BAD_CELLS_CORR ),_MNAMELU( BAD_CELLS_CORR_E ),_MNAMELU( BADLARQ_FRAC ), + _MNAMELU( SIGNIFICANCE ),_MNAMELU( CELL_SIGNIFICANCE ),_MNAMELU( CELL_SIG_SAMPLING ),_MNAMELU( AVG_LAR_Q ),_MNAMELU( AVG_TILE_Q ),_MNAMELU( ENG_BAD_HV_CELLS ), + _MNAMELU( N_BAD_HV_CELLS ),_MNAMELU( PTD ),_MNAMELU( MASS ),_MNAMELU( EM_PROBABILITY ),_MNAMELU( HAD_WEIGHT ),_MNAMELU( OOC_WEIGHT ),_MNAMELU( DM_WEIGHT ), + _MNAMELU( TILE_CONFIDENCE_LEVEL ),_MNAMELU( VERTEX_FRACTION ),_MNAMELU( NVERTEX_FRACTION ),_MNAMELU( ETACALOFRAME ),_MNAMELU( PHICALOFRAME ),_MNAMELU( ETA1CALOFRAME ), + _MNAMELU( PHI1CALOFRAME ),_MNAMELU( ETA2CALOFRAME ),_MNAMELU( PHI2CALOFRAME ),_MNAMELU( ENG_CALIB_TOT ),_MNAMELU( ENG_CALIB_OUT_L ),_MNAMELU( ENG_CALIB_OUT_M ), + _MNAMELU( ENG_CALIB_OUT_T ),_MNAMELU( ENG_CALIB_DEAD_L ),_MNAMELU( ENG_CALIB_DEAD_M ),_MNAMELU( ENG_CALIB_DEAD_T ),_MNAMELU( ENG_CALIB_EMB0 ),_MNAMELU( ENG_CALIB_EME0 ), + _MNAMELU( ENG_CALIB_TILEG3 ),_MNAMELU( ENG_CALIB_DEAD_TOT ),_MNAMELU( ENG_CALIB_DEAD_EMB0 ),_MNAMELU( ENG_CALIB_DEAD_TILE0 ),_MNAMELU( ENG_CALIB_DEAD_TILEG3 ), + _MNAMELU( ENG_CALIB_DEAD_EME0 ),_MNAMELU( ENG_CALIB_DEAD_HEC0 ),_MNAMELU( ENG_CALIB_DEAD_FCAL ),_MNAMELU( ENG_CALIB_DEAD_LEAKAGE ),_MNAMELU( ENG_CALIB_DEAD_UNCLASS ), + _MNAMELU( ENG_CALIB_FRAC_EM ),_MNAMELU( ENG_CALIB_FRAC_HAD ),_MNAMELU( ENG_CALIB_FRAC_REST ),_MNAMELU( ENERGY_DigiHSTruth ),_MNAMELU( ETA_DigiHSTruth ),_MNAMELU( PHI_DigiHSTruth ), + _MNAMELU( TIME_DigiHSTruth ),_MNAMELU( ENERGY_CALIB_DigiHSTruth ),_MNAMELU( ETA_CALIB_DigiHSTruth ),_MNAMELU( PHI_CALIB_DigiHSTruth ),_MNAMELU( TIME_CALIB_DigiHSTruth ), + _MNAMELU( FIRST_PHI_DigiHSTruth ),_MNAMELU( FIRST_ETA_DigiHSTruth ),_MNAMELU( SECOND_R_DigiHSTruth ),_MNAMELU( SECOND_LAMBDA_DigiHSTruth ),_MNAMELU( DELTA_PHI_DigiHSTruth ), + _MNAMELU( DELTA_THETA_DigiHSTruth ),_MNAMELU( DELTA_ALPHA_DigiHSTruth ),_MNAMELU( CENTER_X_DigiHSTruth ),_MNAMELU( CENTER_Y_DigiHSTruth ),_MNAMELU( CENTER_Z_DigiHSTruth ), + _MNAMELU( CENTER_MAG_DigiHSTruth ),_MNAMELU( CENTER_LAMBDA_DigiHSTruth ),_MNAMELU( LATERAL_DigiHSTruth ),_MNAMELU( LONGITUDINAL_DigiHSTruth ),_MNAMELU( ENG_FRAC_EM_DigiHSTruth ), + _MNAMELU( ENG_FRAC_MAX_DigiHSTruth ),_MNAMELU( ENG_FRAC_CORE_DigiHSTruth ),_MNAMELU( FIRST_ENG_DENS_DigiHSTruth ),_MNAMELU( SECOND_ENG_DENS_DigiHSTruth ), + _MNAMELU( ISOLATION_DigiHSTruth ),_MNAMELU( ENG_BAD_CELLS_DigiHSTruth ),_MNAMELU( N_BAD_CELLS_DigiHSTruth ),_MNAMELU( N_BAD_CELLS_CORR_DigiHSTruth ), + _MNAMELU( BAD_CELLS_CORR_E_DigiHSTruth ),_MNAMELU( BADLARQ_FRAC_DigiHSTruth ),_MNAMELU( ENG_POS_DigiHSTruth ),_MNAMELU( SIGNIFICANCE_DigiHSTruth ), + _MNAMELU( CELL_SIGNIFICANCE_DigiHSTruth ),_MNAMELU( CELL_SIG_SAMPLING_DigiHSTruth ),_MNAMELU( AVG_LAR_Q_DigiHSTruth ),_MNAMELU( AVG_TILE_Q_DigiHSTruth ), + _MNAMELU( ENG_BAD_HV_CELLS_DigiHSTruth ),_MNAMELU( N_BAD_HV_CELLS_DigiHSTruth ),_MNAMELU( EM_PROBABILITY_DigiHSTruth ),_MNAMELU( HAD_WEIGHT_DigiHSTruth ), + _MNAMELU( OOC_WEIGHT_DigiHSTruth ),_MNAMELU( DM_WEIGHT_DigiHSTruth ) + }; ///< Moment names by moment types + static const std::map clusterMomentTypes = { + _MIDLU( FIRST_PHI ),_MIDLU( FIRST_ETA ),_MIDLU( SECOND_R ),_MIDLU( SECOND_LAMBDA ),_MIDLU( DELTA_PHI ),_MIDLU( DELTA_THETA ),_MIDLU( DELTA_ALPHA ), + _MIDLU( CENTER_X ),_MIDLU( CENTER_Y ),_MIDLU( CENTER_Z ),_MIDLU( CENTER_MAG ),_MIDLU( CENTER_LAMBDA ),_MIDLU( LATERAL ),_MIDLU( LONGITUDINAL ), + _MIDLU( ENG_FRAC_EM ),_MIDLU( ENG_FRAC_MAX ),_MIDLU( ENG_FRAC_CORE ),_MIDLU( FIRST_ENG_DENS ),_MIDLU( SECOND_ENG_DENS ),_MIDLU( ENG_POS ), + _MIDLU( ISOLATION ),_MIDLU( ENG_BAD_CELLS ),_MIDLU( N_BAD_CELLS ),_MIDLU( N_BAD_CELLS_CORR ),_MIDLU( BAD_CELLS_CORR_E ),_MIDLU( BADLARQ_FRAC ), + _MIDLU( SIGNIFICANCE ),_MIDLU( CELL_SIGNIFICANCE ),_MIDLU( CELL_SIG_SAMPLING ),_MIDLU( AVG_LAR_Q ),_MIDLU( AVG_TILE_Q ),_MIDLU( ENG_BAD_HV_CELLS ), + _MIDLU( N_BAD_HV_CELLS ),_MIDLU( PTD ),_MIDLU( MASS ),_MIDLU( EM_PROBABILITY ),_MIDLU( HAD_WEIGHT ),_MIDLU( OOC_WEIGHT ),_MIDLU( DM_WEIGHT ), + _MIDLU( TILE_CONFIDENCE_LEVEL ),_MIDLU( VERTEX_FRACTION ),_MIDLU( NVERTEX_FRACTION ),_MIDLU( ETACALOFRAME ),_MIDLU( PHICALOFRAME ),_MIDLU( ETA1CALOFRAME ), + _MIDLU( PHI1CALOFRAME ),_MIDLU( ETA2CALOFRAME ),_MIDLU( PHI2CALOFRAME ),_MIDLU( ENG_CALIB_TOT ),_MIDLU( ENG_CALIB_OUT_L ),_MIDLU( ENG_CALIB_OUT_M ), + _MIDLU( ENG_CALIB_OUT_T ),_MIDLU( ENG_CALIB_DEAD_L ),_MIDLU( ENG_CALIB_DEAD_M ),_MIDLU( ENG_CALIB_DEAD_T ),_MIDLU( ENG_CALIB_EMB0 ),_MIDLU( ENG_CALIB_EME0 ), + _MIDLU( ENG_CALIB_TILEG3 ),_MIDLU( ENG_CALIB_DEAD_TOT ),_MIDLU( ENG_CALIB_DEAD_EMB0 ),_MIDLU( ENG_CALIB_DEAD_TILE0 ),_MIDLU( ENG_CALIB_DEAD_TILEG3 ), + _MIDLU( ENG_CALIB_DEAD_EME0 ),_MIDLU( ENG_CALIB_DEAD_HEC0 ),_MIDLU( ENG_CALIB_DEAD_FCAL ),_MIDLU( ENG_CALIB_DEAD_LEAKAGE ),_MIDLU( ENG_CALIB_DEAD_UNCLASS ), + _MIDLU( ENG_CALIB_FRAC_EM ),_MIDLU( ENG_CALIB_FRAC_HAD ),_MIDLU( ENG_CALIB_FRAC_REST ),_MIDLU( ENERGY_DigiHSTruth ),_MIDLU( ETA_DigiHSTruth ),_MIDLU( PHI_DigiHSTruth ), + _MIDLU( TIME_DigiHSTruth ),_MIDLU( ENERGY_CALIB_DigiHSTruth ),_MIDLU( ETA_CALIB_DigiHSTruth ),_MIDLU( PHI_CALIB_DigiHSTruth ),_MIDLU( TIME_CALIB_DigiHSTruth ), + _MIDLU( FIRST_PHI_DigiHSTruth ),_MIDLU( FIRST_ETA_DigiHSTruth ),_MIDLU( SECOND_R_DigiHSTruth ),_MIDLU( SECOND_LAMBDA_DigiHSTruth ),_MIDLU( DELTA_PHI_DigiHSTruth ), + _MIDLU( DELTA_THETA_DigiHSTruth ),_MIDLU( DELTA_ALPHA_DigiHSTruth ),_MIDLU( CENTER_X_DigiHSTruth ),_MIDLU( CENTER_Y_DigiHSTruth ),_MIDLU( CENTER_Z_DigiHSTruth ), + _MIDLU( CENTER_MAG_DigiHSTruth ),_MIDLU( CENTER_LAMBDA_DigiHSTruth ),_MIDLU( LATERAL_DigiHSTruth ),_MIDLU( LONGITUDINAL_DigiHSTruth ),_MIDLU( ENG_FRAC_EM_DigiHSTruth ), + _MIDLU( ENG_FRAC_MAX_DigiHSTruth ),_MIDLU( ENG_FRAC_CORE_DigiHSTruth ),_MIDLU( FIRST_ENG_DENS_DigiHSTruth ),_MIDLU( SECOND_ENG_DENS_DigiHSTruth ), + _MIDLU( ISOLATION_DigiHSTruth ),_MIDLU( ENG_BAD_CELLS_DigiHSTruth ),_MIDLU( N_BAD_CELLS_DigiHSTruth ),_MIDLU( N_BAD_CELLS_CORR_DigiHSTruth ), + _MIDLU( BAD_CELLS_CORR_E_DigiHSTruth ),_MIDLU( BADLARQ_FRAC_DigiHSTruth ),_MIDLU( ENG_POS_DigiHSTruth ),_MIDLU( SIGNIFICANCE_DigiHSTruth ), + _MIDLU( CELL_SIGNIFICANCE_DigiHSTruth ),_MIDLU( CELL_SIG_SAMPLING_DigiHSTruth ),_MIDLU( AVG_LAR_Q_DigiHSTruth ),_MIDLU( AVG_TILE_Q_DigiHSTruth ), + _MIDLU( ENG_BAD_HV_CELLS_DigiHSTruth ),_MIDLU( N_BAD_HV_CELLS_DigiHSTruth ),_MIDLU( EM_PROBABILITY_DigiHSTruth ),_MIDLU( HAD_WEIGHT_DigiHSTruth ), + _MIDLU( OOC_WEIGHT_DigiHSTruth ),_MIDLU( DM_WEIGHT_DigiHSTruth ) + }; ///< Moment names buy moment identifiers + static const std::string& getMomentName(xAOD::CaloCluster::MomentType momentType); ///< Get moment name associated with a moment type + static bool getMomentType(const std::string& momentName,xAOD::CaloCluster::MomentType& momentType); ///< Get moment type associated with a moment name + static bool haveMomentType(const std::string& momentName); ///< Returns @c true in case moment is known + ///@} + } // Lookup } // CaloRec + +inline const std::string& CaloRec::Lookup::getSamplingName(CaloSampling::CaloSample sid) { return samplingNames.find(sid)->second; } +inline CaloSampling::CaloSample CaloRec::Lookup::getSamplingId(const std::string& sname) { + auto fid(samplingIds.find(sname)); return fid != samplingIds.end() ? fid->second : samplingIds.find("Unknown")->second; +} + +inline bool CaloRec::Lookup::haveMomentType(const std::string& momentName) { return clusterMomentTypes.find(momentName) != clusterMomentTypes.end(); } +inline const std::string& CaloRec::Lookup::getMomentName(xAOD::CaloCluster::MomentType momentType) { return clusterMomentNames.at(momentType); } +inline bool CaloRec::Lookup::getMomentType(const std::string& momentName,xAOD::CaloCluster::MomentType& momentType) { + bool isOk(haveMomentType(momentName)); if ( isOk ) { momentType = clusterMomentTypes.at(momentName); } return isOk; +} + + #endif diff --git a/Calorimeter/CaloRec/CaloRec/CaloTopoClusterFromTowerMaker.h b/Calorimeter/CaloRec/CaloRec/CaloTopoClusterFromTowerMaker.h index c9f493dc982eecdde5c04799d4b8c76f71f05c1b..858bb1dfbcda02d2b2e01266c6e828fc3aa0fdab 100644 --- a/Calorimeter/CaloRec/CaloRec/CaloTopoClusterFromTowerMaker.h +++ b/Calorimeter/CaloRec/CaloRec/CaloTopoClusterFromTowerMaker.h @@ -1,248 +1,254 @@ // -*- c++ -*- -/* Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ +/* Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration */ #ifndef CALOREC_CALOTOPOCLUSTERFROMTOWERMAKER_H #define CALOREC_CALOTOPOCLUSTERFROMTOWERMAKER_H -#include "AthenaBaseComps/AthAlgTool.h" +#include "StoreGate/ReadHandleKey.h" +#include "StoreGate/WriteHandleKey.h" + +#include "GaudiKernel/ServiceHandle.h" -#include "CaloEvent/CaloCellTowerLink.h" +#include "AthenaBaseComps/AthAlgTool.h" #include "CaloRec/CaloClusterCollectionProcessor.h" +#include "CaloRec/CaloTowerGeometrySvc.h" #include "CaloRec/CaloProtoCluster.h" +#include "CaloEvent/CaloCell.h" +#include "CaloEvent/CaloClusterCellLink.h" +#include "CaloEvent/CaloCellClusterWeights.h" +#include "CaloEvent/CaloCellContainer.h" + +#include "CaloGeoHelpers/CaloSampling.h" + #include "xAODCaloEvent/CaloClusterContainer.h" -#include "xAODCaloEvent/CaloTowerContainer.h" #include -#include - #include -#include +#include +#include -class CaloClusterCellLink; -class CaloCell; -class CaloCellClusterWeights; +#ifndef _CALOTOPOCLUSTERFROMTOWERMAKER_BITSET_SIZE +#define _CALOTOPOCLUSTERFROMTOWERMAKER_BITSET_SIZE 28 +#endif class CaloTopoClusterFromTowerMaker : public AthAlgTool, virtual public CaloClusterCollectionProcessor { public: - /// @brief Tool constructor + ///@brief Tool constructor CaloTopoClusterFromTowerMaker(const std::string& type,const std::string& name,const IInterface* pParent); - /// @brief Tool initialization - StatusCode initialize(); - - /// @brief Tool execution - StatusCode execute(xAOD::CaloClusterContainer* pClusCont); - - /// @brief Tool finalize - StatusCode finalize(); + ///@name @c AthAlgTool and @c CaloClusterCellProcessor interface implementations + ///@{ + StatusCode initialize(); ///< Setting up the operational mode and corresponding parameters + StatusCode execute(xAOD::CaloClusterContainer* pClusCont); ///< Execute the tool and fill the @c xAOD::CaloClusterContainer pointed to by @c pClusCont + StatusCode finalize(); ///< Finalize the tool (no action) + ///@} private: - /// @name Tool properties - /// @{ - std::string m_towerContainerKey; ///> @brief Key (name) of input tower container - std::string m_cellTowerLinksKey; ///> @brief Derived property (cannot be set by user) - std::string m_cellContainerKey; ///> @brief Key (name) of input cell container - std::string m_clusterContainerKey; ///> @brief Key of topo-cluster container to be used as source for cells - std::string m_cellClusterWeightKey; ///> @brief Key for cell-cluster weight look up object - bool m_orderByPt; ///> @brief Orders cluster container by @f$ p_{\text{T}} @f$, default @c true - double m_energyThreshold; ///> @brief Cell energy threshold - bool m_applyLCW; ///> @brief Prepare LCW calibration - /// @} + ///@name Internally used types + ///@{ + typedef std::vector protocont_t; ///< Container for @c CaloProtoCluster objects + typedef std::size_t uint_t; ///< Unsigned integral type + ///@} - /// @name Processing flags + /// @name Tool properties /// @{ - /// @brief Energy threshold flag - bool m_applyEnergyThreshold; - /// @brief Use cells from topo-clusters flag - bool m_useCellsFromClusters; - /// @brief Prepare LCW calibration - bool m_prepareLCW; + ServiceHandle m_towerGeometrySvc; ///< Tower geometry service + SG::ReadHandleKey m_clusterContainerKey; ///< Topo-cluster container key + SG::ReadHandleKey m_cellContainerKey; ///< Calorimeter cell container + SG::WriteHandleKey m_cellClusterWeightKey; ///< Cell signal weights in clusters key + bool m_orderByPt; ///< Orders cluster container by @f$ p_{\text{T}} @f$, default @c true + bool m_prepareLCW; ///< Prepare LCW calibration, default is @c false + bool m_useCellsFromClusters; ///< Use cells from topo-clusters if @c true, else use all cells, default is @c true + bool m_applyCellEnergyThreshold; ///< Apply cell energy threshold, default is @c false + bool m_doCellIndexCheck; ///< Check cell hash index consistency if @c true (default @c false) + bool m_buildCombinedSignal; ///< Build topo-clusters within given @f$ y @f$ range, else topo-towers + double m_energyThreshold; ///< Cell energy threshold, default is set in @c m_energyThresholdDef + double m_clusterRange; ///< Range where topo-clusters are used when m_buildCombinedSignal = true /// @} /// @name Constants and parameters /// @{ - /// @brief Number of cells - int m_numberOfCells; - /// @brief Number of samplings - int m_numberOfSamplings; - /// @brief Default energy threshold - static double m_energyThresholdDef; - /// @brief Default container key - static std::string m_defaultKey; + uint_t m_numberOfCells; ///< Number of cells (highest cell index + 1) + uint_t m_maxCellHash; ///< Maximum hash index of cell ( number of cells - 1) + uint_t m_numberOfSamplings; ///< Number of samplings + uint_t m_numberOfTowers; ///< Number of towers + static double m_energyThresholdDef; ///< Default energy threshold + static double m_clusterRangeDef; ///< Default cluster @f$ y @f$ range + static std::string m_defaultKey; ///< Default container key + static uint_t m_errorValueUINT; ///< Error value for @c uint_t type values /// @} - /// @brief Access to cell-to-tower links - const CaloCellTowerLink::Map* m_cellTowerLinksHandle; - - /// @name Cache for cell weights - /// @{ - /// @brief Data structure tagging a cell as a cluster cell and storing its weight - typedef boost::tuples::tuple celltag_t; - /// @brief Random access lookup of cell tags - typedef std::vector celltagstore_t; - /// @brief Cell tag store - celltagstore_t m_cellTagStore; - /// @brief Previously used cell indices - std::vector m_cellIdxStore; - /// @} + ///@name Internally used helpers + ///@{ + xAOD::CaloCluster::ClusterSize getClusterSize(uint_t etaBins,uint_t phiBins); ///< Returns a cluster size tag from number of eta and phi bins in tower grid + xAOD::CaloCluster::ClusterSize getClusterSize(uint_t towerBins); ///< Returns a cluster size tag from number of towers (bins) in tower grid + int cleanupCells(CaloClusterCellLink* clk,uint_t nclus); ///< Checks @c CaloClusterCellLink for consistency + ///@} - /// @brief Pathologies - /// - /// Internally used object collects pathologies observed for @c CaloCell objects in execution. - /// Each pathology is identified by a pre-defined message, and the number of pathologies is produced by - /// counting the number of occurancies of a given message. - class CellPathology - { - public: - /// @name Internally used data types - /// @{ - typedef std::string pathology_t; - typedef int count_t; - typedef std::map cache_t; - /// @} - - /// @brief Constructor - CellPathology(); - /// @name Set pathologies - /// @{ - /// @brief Invalid @c CaloCell pointer - /// - /// An invalid @c NULL pointer was found at a given valid index in the @c CaloCellContainer. - /// - /// @param cellhash hash identifier (index in @c CaloCellContainer) for a given @c CaloCell. - /// @param cptr cell pointer (should be @c NULL in this case) - void setInvalidPointer(int cellhash,const CaloCell* cptr); - /// @brief Invalid @c CaloCell hash - /// - /// A cell identifier is found to be out of range. - /// - /// @param cellhash hash identifier (index in @c CaloCellContainer - here invalid!) for a given @c CaloCell. - /// @param cptr cell pointer (should be @c NULL in this case) - void setInvalidHash(int cellhash,const CaloCell* cptr); - /// @brief Invalid sampling identifier for given @c CaloCell - /// - /// The sampling identifier associated with a @c CaloCell is not valid (not in the enumerator list). - /// - /// @param cellhash hash identifier (index in @c CaloCellContainer - here invalid!) for a given @c CaloCell. - /// @param cptr cell pointer (should be @c NULL in this case) - void setInvalidSampling(int cellhash,const CaloCell* cptr); - /// @} - - /// @name Accessors to pathologies stats - /// @{ - const cache_t& invalidPointer() const; ///< @brief Occurances of invalid pointers - const cache_t& invalidHash() const; ///< @brief Occurances of invalid hash identifier - const cache_t& invalidSampling() const; ///< @brief Occurances of invalid sampling identifier - /// @} - - private: - - /// @name Caches - /// @{ - cache_t _invalidPointer; ///< @brief Invalid pointer - cache_t _invalidHash; ///< @brief Invalid hash - cache_t _invalidSampling; ///< @brief Invalid sampling - /// @} - - /// @name Helpers - /// @{ - /// @brief Add a pathology to the cache - /// - /// @param msg message describing pathology - /// @param map cache - void _addMsg(const pathology_t& msg,cache_t& map); - /// @} - }; - - /// @brief Cluster size tag - xAOD::CaloCluster::ClusterSize getClusterSize(const xAOD::CaloTowerContainer& towerCont); - - /// @brief Clean up cell links - remove null pointers - int cleanupCells(CaloClusterCellLink* clk); - - // /// @brief Calculate kinemactics - // bool calculateKine(xAOD::CaloCluster* pClus); - - /// @brief Book keeping - void addWarning(const std::string& msg); - std::map m_warning; - void addMessage(const std::string& msg); - std::map m_message; - CellPathology m_cellProblems; - void printProblems(const std::map& map); - - typedef std::vector protocont_t; - - /// @name Helpers - /// @{ - /// @brief Build inclusive tower clusters - /// - /// Fills all cells into proto-clusters. - /// - /// @return @c false in case of problems with data access or inconsistent data structures. - /// - /// @param pCellCont reference to non-modifiable @c CaloCellContainer - /// @param pProtoCont reference to @c CaloProtoCluster container filled on output. - /// - bool buildInclClusters(const CaloCellContainer& pCellCont,protocont_t& pProtoCont); - /// @brief Build exclusive tower clusters + ///@name Tower builders /// - /// Fills cells above energy threshold into proto-clusters. + ///@return @c false in case of problems with data access or inconsistent data structures /// - /// @return @c false in case of problems with data access or inconsistent data structures. + ///@param pCellCont reference to non-modifiable @c CaloCellContainer + ///@param pProtoCont reference to @c CaloProtoCluster container filled on output. + ///@param clusCont reference to non-modifiable @c xAOD::CaloClusterContainer + ///@param protoCont reference to modifiable proto-cluster container /// - /// @param pCellCont reference to non-modifiable @c CaloCellContainer - /// @param pProtoCont reference to @c CaloProtoCluster container filled on output. - bool buildExclClusters(const CaloCellContainer& pCellCont,protocont_t& pProtoCont); - // bool buildFilteredClusters(const CaloCellContainer& pCellCont,protocont_t& pProtoCont); ///< @brief Build tower clusters from topo-cluster cells - /// @brief Add cells to proto-clusters + ///@return + ///@{ + uint_t buildInclTowers(const CaloCellContainer& pCellCont,protocont_t& pProtoCont); ///< Inclusive towers + uint_t buildExclTowers(const CaloCellContainer& pCellCont,protocont_t& pProtoCont); ///< Exclusive towers + uint_t buildEMTopoTowers(const xAOD::CaloClusterContainer& clusCont,protocont_t& protoCont); ///< EM topo-towers + uint_t buildLCWTopoTowers(const xAOD::CaloClusterContainer& clusCont,protocont_t& protoCont); ///< LCW topo-towers + ///@} + /// @brief Adding cells to proto-clusters /// /// @return @c true if cell successfully added to one (or more) proto-clusters /// /// @param cptr pointer ton non-modifiable @c CaloCell object /// @param pProtoCont reference to proto-cluster container - bool addCellToProtoCluster(const CaloCell* cptr,protocont_t& pProtoCont); - /// @} - - /// @brief Cache - CaloCellClusterWeights* m_cellClusterWeights; + /// @param weight additional (global) weight of cell (e.g. for geometrical weight for combined EM-scale signals) + bool addCellToProtoCluster(const CaloCell* cptr,protocont_t& pProtoCont,double weight=1.); + + ///@name Helpers + ///@{ + bool filterProtoCluster(const CaloClusterCellLink& clnk) const; ///< Checks for and removes invalid cell links + bool checkCellIndices(const CaloCellContainer* pCellCont) const; ///< Checks consistency between cell indices and hash identifiers + bool isValidIndex(uint_t idx) const; ///< Checks if argument is a valid index value + uint_t badIndexValue() const; ///< Returns value indicating a bad index + ///@} + + ///@name Excluded samplings + ///@{ + std::vector m_excludedSamplings; ///< List of excluded samplings (@c CaloSampling::CaloSample enumerators) + std::vector m_excludedSamplingsName; ///< List of excluded samplings (human-readable names) + std::bitset< _CALOTOPOCLUSTERFROMTOWERMAKER_BITSET_SIZE > m_excludedSamplingsPattern; ///< Bit pattern indicates if sampling is excluded + ///@} + + ///@name Monitoring + ///@{ + ///@} }; -inline void CaloTopoClusterFromTowerMaker::CellPathology::_addMsg(const std::string& msg,std::map& map) -{ if ( map.find(msg) == map.end() ) { map[msg] = 0; } map[msg] += 1; } +inline CaloTopoClusterFromTowerMaker::uint_t CaloTopoClusterFromTowerMaker::badIndexValue() const { return m_errorValueUINT; } +inline bool CaloTopoClusterFromTowerMaker::isValidIndex(uint_t idx) const { return idx != badIndexValue(); } -inline const std::map& CaloTopoClusterFromTowerMaker::CellPathology::invalidPointer() const -{ return _invalidPointer; } -inline const std::map& CaloTopoClusterFromTowerMaker::CellPathology::invalidHash() const -{ return _invalidHash; } -inline const std::map& CaloTopoClusterFromTowerMaker::CellPathology::invalidSampling() const -{ return _invalidSampling; } - -/// @class CaloTopoClusterFromTowerMaker -/// @author Peter Loch +///@class CaloTopoClusterFromTowerMaker /// -/// @brief A cluster builder tool forming topo-clusters representing calorimeter tower signals on a regular grid in @f$ (\eta,\phi) @f$ space. +/// @brief A cluster builder tool forming topo-clusters representing calorimeter tower signals on a regular grid in @f$ (\eta,\phi) @f$ space. By default, +/// EM-scale topo-towers are created from cells in topo-clusters. /// -/// This tool provides functionality allowing the creation of topo-cluster collections representing towers. The output data objects are -/// @c xAOD::CaloCluster , collected into @c xAOD::CaloClusterContainer . Several configuration options are supported: +/// This tool fills EM-scale towers and stores them as @c xAOD::CaloCluster. It supports several operational modes, which are +/// controlled by tool properties. It fills a container of type @c xAOD::CaloClusterContainer. The properties controlling its +/// specific behavior are: /// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +/// +///
Properties defining tool behavior
Property nameProperty typeDefault valueComment
OrderClusterByPtboolfalseif @c true, the @c xAOD::CaloClusterContainer is ordered by @f$ p_{\rm T}^{\rm clus} @f$. See further comments below.
PrepareLCWboolfalseif @c true, the tool fills a @c CaloCellClusterWeights object and records it into the event store to be used by @c CaloTopoClusterFromTowerCalibrator
UseCellsFromClustersbooltrueif @c true, only cells from topo-clusters are used to fill the towers (topo-towers); else, @a inclusive @a towers are filled with all cells.
Properties setting variables for operational modes
Property nameProperty typeDefault valueComment
CellEnergyThresholddoublem_energyThresholdDefcell energy threshold used in exclusive mode only. See further comments below.
CellContainerKeySG::ReadHandleKey"AllCalo"cell container key is needed to pick up @c CaloCellContainer for all operational modes.
ClusterContainerKeySG::ReadHandleKey"CaloTopoCluster"cluster container key is needed to pick up @c xAOD::CaloClusterContainer for filtered mode (UseCellsFromCluster = true) +///
CellClusterWeightKeySG::WriteHandleKey−N/A−key for @c CaloCellClusterWeights object is needed if PrepareLCW = true. Default is empty key. +///
BuildCombinedTopoSignalboolfalseturns on combined topo-cluster/topo-tower output, with topo-clusters used within the rapidity range defined by TopoClusterRange and topo-towers elsewhere.
TopoClusterRangedouble5.sets the range @f$ y_{\rm topo-cluster}^{\rm max} @f$ for using topo-clusters when BuildCombinedTopoSignal = true; +/// topo-clusters with @f$ \left|y_{\rm topo-cluster}\right| < y_{\rm topo-cluster}^{\rm max} @f$ are used. +///
+/// +/// The towers can be classified as: /// -# inclusive cell towers -/// All cells are collected into towers, independent of their signal. This is the default behaviour. +/// All cells are collected into inclusive towers, independent of their signal. Requires properties UseCellsFromClusters = false and UseCellEnergyThreshold = false. Only EM +/// towers are possible, as cells not collected into topo-clustersdo not have a calibration applied. /// -# exclusive cell towers -/// Cells with @f$ E > E_{\rm min} @f$ are collected into towers. This behaviour is turned on by poviding a property @c CellEnergyThreshold setting @f$ E_{\rm min} @f$. +/// Cells with @f$ E > E_{\rm min} @f$ are collected into exclusive towers. This behaviour is turned on by UseCellsFromClusters = false and UseCellEnergyThreshold = true. A +/// meaningful CellEnergyThreshold value needs to be provided in addition. /// -# filtered mode -/// Cells contributing to standard topo-clusters are collected into towers. This behaviour is triggered by providing the @c CaloTopoClusterContainerKey property with a -/// valid key for a topo-cluster collection in the event store. +/// Cells contributing to standard topo-clusters are collected into topo-towers. This behaviour is triggered by UseCellsFromClusters = true. Optionally, LCW calibration can be applied +/// to these towers by setting PrepareLCW = true and scheduling a @c CaloTopoClusterFromTowerCalibrator tool after the cluster moment calculators. The values of the UseEnergyThreshold +/// and CellEnergyThreshold properties are ignored in this mode. A valid event store key needs to be provided in the to pick up the topo-cluster container. Note that building EM +/// topo-towers requires topo-clusters on EM scale (no LCW applied) to get the correct geometrical cell weights only. LCW topo-towers require LCW scale topo-clusters to get the correct full geometrical +/// and calibration weights. +/// -# mixed mode +/// Cells contributing to standard topo-clusters are collected into towers if these topo-clusters are outside of a give rapidity range. The rapidity range is defined by the TopoClusterRange +/// property. This mode is turned on by setting the property BuildCombinedTopoSignal = true. It is turned off by default (BuildCombinedTopoSignal = false). +/// EM scale and LCW scale is possible, as in the filtered mode. /// /// Configuration 2 and 3 are exclusive, with 3 overwriting 2. The output topo-clusters represent calorimeter towers on the EM scale. The can be handed to cluster moment /// tools (needs EM scale) and, if desired, to a dedicated cluster calibration tool of type @c xAOD::CaloTowerClusterFromTowerCalibrator . /// -/// To avoid mulitple retrievals of the same weights by searching for cells in (many) topo-clusters, the overall weight of the cell signal is stored in a random access +/// To avoid multiple retrievals of the same weights by searching for cells in (many) topo-clusters, the overall weight of the cell signal is stored in a random access /// look-up table stored in a @c CaloCellClusterWeights object in the detector store. This object is created by this tool, if needed. If the tool property -/// @c CellWeightLookupKey is set, this object will be generated, filled, and recorded. -/// +/// @c CellWeightLookupKey is set, this object will be generated, filled, and recorded. This is essential for calibrated topo-towers! +/// +///@note The @c OrderByPt property, which orders the container by descending transverse momentum, is only useful for EM towers. Applyin LCW may lead to a different +/// order - if a container with LCW towers should be ordered, the corresponding property @c OrderByPt of the @c CaloTopoClusterFromTowerCalibrator tool should +/// be set to @c true. +/// +///@note Many more details on the towers are available on +/// this page. +/// +/// @author Peter Loch #endif diff --git a/Calorimeter/CaloRec/CaloRec/CaloTopoClusterFromTowerMonitor.h b/Calorimeter/CaloRec/CaloRec/CaloTopoClusterFromTowerMonitor.h new file mode 100644 index 0000000000000000000000000000000000000000..52abb14cae72c62adc049bf768c2d524b454455d --- /dev/null +++ b/Calorimeter/CaloRec/CaloRec/CaloTopoClusterFromTowerMonitor.h @@ -0,0 +1,342 @@ +// -*- c++ -*- +#ifndef CALOREC_CALOTOPOCLUSTERFROMTOWERMONITOR_H +#define CALOREC_CALOTOPOCLUSTERFROMTOWERMONITOR_H + +#include "StoreGate/ReadHandleKey.h" + +#include "GaudiKernel/ServiceHandle.h" + +#include "AthenaBaseComps/AthHistogramAlgorithm.h" + +#include "xAODCaloEvent/CaloClusterContainer.h" + +#include "CaloRec/CaloTowerGeometrySvc.h" +#include "CaloRec/CaloTopoClusterFromTowerHelpers.h" + +#include "xAODCaloEvent/CaloCluster.h" +#include "CaloEvent/CaloCell.h" + +#include +#include +#include +#include + +#include "TH1D.h" +#include "TH2D.h" + +///@brief Algorithm to monitor kinematics and composition of towers (and clusters) +/// +/// This algorithm fills a few histograms relevant for monitoring tower-related information, As both topo-clusters and towers +/// use the same representation in data, both signal definitions can be monitored - even though most histograms are designed +/// for distributions of tower-specific signal features. +/// +/// The algorithm has several process control properties and configured values: +/// +/// | Property | Default | Controled feature or value +/// :------------------ | :----------: |:------------------------------------------------------------------------------------------- +/// | @c DoGeoAutoBins | @c true | if @c true the binning of @f$ y @f$ and @f$ \phi @f$ is constructed from the tower geometry database; +/// | ^ | ^ | else it is user-defined by the corresponding value property +/// | Hotspot analysis is on if @f$ y_{\rm min}^{\rm hotspot} < y_{\rm max}^{\rm hotspot} @f$ and @f$ \phi_{\rm min}^{\rm hotspot} < \phi_{\rm max}^[\rm hotspot} @f$ (on by defaults) || +/// | @c EtaMinHotspot | 0. | defines lower limit @f$ y_{\rm min}^{\rm hotspot} @f$ of rapidity for the hotspot analysis +/// | @c EtaMaxHotspot | 0.1 | defines upper limit @f$ y_{\rm max}^{\rm hotspot} @f$ of rapidity for the hotspot analysis +/// | @c PhiMinHotspot | 0. | defines lower limit @f$ \phi_{\rm min}^{\rm hotspot} @f$ of azimuth for the hotspot analysis +/// | @c PhiMaxHotspot | 0.1 | defines upper limit @f$ \phi_{\rm max}^{\rm hotspot} @f$ of azimuth for the hotspot analysis +/// | Histogram binning: (pseudo)-rapidity (@f$ y_{\rm min} \geq y_{\rm max} @f$ returns error) || +/// | @c EtaTowersBins | 100 | number of rapidity bins +/// | @c EtaTowersMin | -5. | lower limit of rapidity range @f$ y_{\rm min} @f$ +/// | @c EtaTowersMax | 5. | upper limit of rapidity range @f$ y_{\rm max} @f$ +/// | Histogram binning: azimuth (@f$ \phi_{\rm min \geq \phi_[\rm max} @f$ returns error) || +/// | @c PhiTowersBins | 64 | number of azimuth bins +/// | @c PhiTowersMin | @f$ -\pi @f$ | lower limit of azimuth range @f$ \phi_{\rm min} @f$ +/// | @c PhiTowersMax | @f$ +\pi @f$ | upper limit of azimuth range @f$ \phi_{\rm max} @f$ +/// | Histogram binning: transverse momentum (@f$ p_{\rm T,min} < p_{\rm T,max} @f$ returns error) || +/// | @c PtTowersBins | 220 | number of transverse momentum bins +/// | @c PtTowersMin | -10 GeV | lower limit of transverse momentum range @f$ p_{\rm T,min} @f$ +/// | @c PtTowersMax | 100 GeV | upper limit of transverse momentum range @f$ p_{\rm T,max} @f$ +class CaloTopoClusterFromTowerMonitor : public AthHistogramAlgorithm +{ +public: + + /// Default algorithm constructor + CaloTopoClusterFromTowerMonitor(const std::string& name,ISvcLocator* pSvcLocator); + /// Base-class destructor + virtual ~CaloTopoClusterFromTowerMonitor(); + + ///@brief Initialization + /// + /// This method configures the algorithm, allocates the tower geometry service and books the histograms by + /// invoking the @c book() method. In addition, histograms are filled with geometry information retrieved + /// from the tower geometry service (static information not depending on event variables). + virtual StatusCode initialize(); + ///@brief Execution + /// + /// This method allocates the input data container and fills all histograms. If configured, it also provides + /// + virtual StatusCode execute(); ///< Execution fills histograms. + + +private: + + ///@name Data access properties + ///@{ + SG::ReadHandleKey m_towerContainerKey; ///< Allocator for input @c xAOD::CaloClusterContainer + ServiceHandle m_towerGeometrySvc; ///< Allocator for tower geometry services + ///@} + + ///@name Histogram binning properties + ///@{ + int m_ncBins; ///< Number of cells in towers - number of bins + double m_ncMin; ///< Number of cells in towers - lower limit of value range + double m_ncMax; ///< Number of cells in towers - upper limit of value range + int m_nBins; ///< Tower multiplicity - number of bins + double m_nMin; ///< Tower multiplicity - lower limit of value range + double m_nMax; ///< Tower multiplicity - upper limit of value range + int m_ptBins; ///< Tower @f$ p_{\rm T} @f$ - number of bins + double m_ptMin; ///< Tower @f$ p_{\rm T} @f$ - lower limit of value range (in GeV) + double m_ptMax; ///< Tower @f$ p_{\rm T} @f$ - upper limit of value range (in GeV) + int m_etaBins; ///< Tower rapidity - number of bins + double m_etaMin; ///< Tower rapidity - lower limit of value range + double m_etaMax; ///< Tower rapidity - upper limit of value range + int m_phiBins; ///< Tower azimuth - number of bins + double m_phiMin; ///< Tower azimuth - lower limit of value range + double m_phiMax; ///< Tower azimuth - upper limit of value range + + double m_hsEtaMin; double m_hsEtaMax; + double m_hsPhiMin; double m_hsPhiMax; + + bool m_doGeoAutoBins; + bool m_doHotspot; + ///@} + + TH1D* h_n; + TH1D* h_pt; + TH1D* h_eta; + TH1D* h_phi; + TH1D* h_nc; + TH1D* h_samp; + + TH2D* d_n_eta_phi; + TH2D* d_nc_eta_phi; + TH2D* d_pt_eta; + TH2D* d_nc_eta; + + TH2D* d_wgt_samp; + TH2D* d_ntt_samp; + TH2D* d_geo_samp; + TH2D* d_maxtowers_samp; + TH2D* d_wgttowers_samp; + + TH2D* d_maxcells_eta; + TH2D* d_allwghts_eta; + + TH2D* d_deta_eta; + TH2D* d_dphi_eta; + TH2D* d_dphi_deta; + + TH2D* d_detac_eta; + TH2D* d_dphic_eta; + TH2D* d_dphic_detac; + + TH2D* d_detac_samp; + TH2D* d_dphic_samp; + + // hot spot + TH1D* h_nc_hs; + TH1D* h_n_hs; + TH1D* h_pt_hs; + TH1D* h_eta_hs; + TH1D* h_phi_hs; + TH1D* h_samp_hs; + + TH2D* d_n_eta_phi_hs; + TH2D* d_nc_eta_phi_hs; + + TH2D* d_deta_eta_hs; + TH2D* d_dphi_eta_hs; + TH2D* d_dphi_deta_hs; + + TH2D* d_detac_eta_hs; + TH2D* d_dphic_eta_hs; + TH2D* d_dphic_detac_hs; + + TH2D* d_detac_samp_hs; + TH2D* d_dphic_samp_hs; + + std::vector d_maxcells_phi_eta_slice; + std::vector d_allwghts_phi_eta_slice; + +protected: + + bool isInHotspot(const xAOD::CaloCluster& ptow) const; + bool fillComposition(const xAOD::CaloCluster& ptow,std::vector& deta,std::vector& dphi,std::vector& csam) const; + bool setAxisTitle(TH1* h,const std::string& title,const std::string& axis="x"); + + virtual StatusCode book(); + std::bitset<200000> m_cellTags; + + /////////////////////// + // BookAny Templates // + /////////////////////// + + template + H* bookAny(const std::string& hname,const std::string& htitle,const std::string& xtitle,int nxbins,double xmin,double xmax) { + H* hptr = (H*)bookGetPointer( H(hname.c_str(),htitle.c_str(),nxbins,xmin,xmax) ); + if ( hptr == 0 ) { + ATH_MSG_WARNING( CaloRec::Helpers::fmtMsg("Cannot book distribution \042%s\042 with title \042%s\042",hname.c_str(),htitle.c_str()) ); + } else { + hptr->Sumw2(); + if ( !xtitle.empty() && xtitle != "" ) { hptr->GetXaxis()->SetTitle(xtitle.c_str()); } + } + return hptr; + } + + template + H* bookAny(const std::string& hname,const std::string& htitle,const std::string& xtitle,int nxbins,double xmin,double xmax,int nybins,double ymin,double ymax) { + H* hptr = (H*)bookGetPointer( H(hname.c_str(),htitle.c_str(),nxbins,xmin,xmax,nybins,ymin,ymax) ); + if ( hptr == 0 ) { + ATH_MSG_WARNING( CaloRec::Helpers::fmtMsg("Cannot book distribution \042%s\042 with title \042%s\042",hname.c_str(),htitle.c_str()) ); + } else { + hptr->Sumw2(); + if ( !xtitle.empty() && xtitle != "" ) { hptr->GetXaxis()->SetTitle(xtitle.c_str()); } + } + return hptr; + } + + template + H* bookAny(const std::string& hname,const std::string& htitle,const std::string& xtitle, + int nxbins,double xmin,double xmax,int nybins,double ymin,double ymax,int nzbins,double zmin,double zmax) { + H* hptr = (H*)bookGetPointer( H(hname.c_str(),htitle.c_str(),nxbins,xmin,xmax,nybins,ymin,ymax,nzbins,zmin,zmax) ); + if ( hptr == 0 ) { + ATH_MSG_WARNING( CaloRec::Helpers::fmtMsg("Cannot book distribution \042%s\042 with title \042%s\042",hname.c_str(),htitle.c_str()) ); + } else { + hptr->Sumw2(); + if ( !xtitle.empty() && xtitle != "" ) { hptr->GetXaxis()->SetTitle(xtitle.c_str()); } + } + return hptr; + } + + template + H* bookAny(const std::string& hname,const std::string& htitle,int nxbins,double xmin,double xmax) { + H* hptr = (H*)bookGetPointer( H(hname.c_str(),htitle.c_str(),nxbins,xmin,xmax) ); + if ( hptr == 0 ) { + ATH_MSG_WARNING( CaloRec::Helpers::fmtMsg("Cannot book distribution \042%s\042 with title \042%s\042",hname.c_str(),htitle.c_str()) ); + } else { + hptr->Sumw2(); + hptr->GetXaxis()->SetTitle(hptr->GetTitle()); + } + return hptr; + } + + template + H* bookAny(const std::string& hname,const std::string& htitle,int nxbins,double xmin,double xmax,int nybins,double ymin,double ymax) { + H* hptr = (H*)bookGetPointer( H(hname.c_str(),htitle.c_str(),nxbins,xmin,xmax,nybins,ymin,ymax) ); + if ( hptr == 0 ) { + ATH_MSG_WARNING( CaloRec::Helpers::fmtMsg("Cannot book distribution \042%s\042 with title \042%s\042",hname.c_str(),htitle.c_str()) ); + } else { + hptr->Sumw2(); + hptr->GetXaxis()->SetTitle(hptr->GetTitle()); + } + return hptr; + } + + template + H* bookAny(const std::string& hname,const std::string& htitle,int nxbins,double xmin,double xmax,int nybins,double ymin,double ymax,int nzbins,double zmin,double zmax) { + H* hptr = (H*)bookGetPointer( H(hname.c_str(),htitle.c_str(),nxbins,xmin,xmax,nybins,ymin,ymax,nzbins,zmin,zmax) ); + if ( hptr == 0 ) { + ATH_MSG_WARNING( CaloRec::Helpers::fmtMsg("Cannot book distribution \042%s\042 with title \042%s\042",hname.c_str(),htitle.c_str()) ); + } else { + hptr->Sumw2(); + hptr->GetXaxis()->SetTitle(hptr->GetTitle()); + } + return hptr; + } + + //////////////////////// + // Book for Samplings // + //////////////////////// + + template + H* bookForSamplings(const std::string& hname,const std::string& htitle) { + int nsamp((int)CaloRec::Lookup::getSamplingId("MINIFCAL0")); double xmin(-0.5); double xmax(xmin+1.*(nsamp+1)); + H* hptr = bookAny(hname,htitle,"",nsamp+1,xmin,xmax); + if ( hptr != 0 ) { + for ( int isamp(0); isamp < nsamp; ++isamp ) { + hptr->GetXaxis()->SetBinLabel(isamp+1,CaloRec::Lookup::getSamplingName((CaloSampling::CaloSample)isamp).c_str()); + } + hptr->GetXaxis()->SetBinLabel(hptr->GetNbinsX(),CaloRec::Lookup::getSamplingName(CaloSampling::Unknown).c_str()); + } + return hptr; + } + + template + H* bookForSamplings(const std::string& hname,const std::string& htitle,int nybins,double ymin,double ymax) { + int nsamp((int)CaloRec::Lookup::getSamplingId("MINIFCAL0")); double xmin(-0.5); double xmax(xmin+1.*(nsamp+1)); + H* hptr = bookAny(hname,htitle,"",nsamp+1,xmin,xmax,nybins,ymin,ymax); + if ( hptr != 0 ) { + for ( int isamp(0); isamp < nsamp; ++isamp ) { + hptr->GetXaxis()->SetBinLabel(isamp+1,CaloRec::Lookup::getSamplingName((CaloSampling::CaloSample)isamp).c_str()); + } + hptr->GetXaxis()->SetBinLabel(hptr->GetNbinsX(),CaloRec::Lookup::getSamplingName(CaloSampling::Unknown).c_str()); + } + return hptr; + } + + /////////////////////// + // Book for Rapidity // + /////////////////////// + + template + H* bookForEta(const std::string& hname,const std::string& htitle) { return bookAny(hname,htitle,"y_{tower}",m_etaBins,m_etaMin,m_etaMax); } + + template + H* bookForEta(const std::string& hname,const std::string& htitle,int nybins,double ymin,double ymax) { + return bookAny(hname,htitle,"y_{tower}",m_etaBins,m_etaMin,m_etaMax,nybins,ymin,ymax); + } + + ////////////////////// + // Book for Azimuth // + ////////////////////// + + template + H* bookForPhi(const std::string& hname,const std::string& htitle) { return bookAny(hname,htitle,"#phi_{tower} [rad]",m_phiBins,m_phiMin,m_phiMax); } + + template + H* bookForPhi(const std::string& hname,const std::string& htitle,int nybins,double ymin,double ymax) { + return bookAny(hname,htitle,"#phi_{tower} [rad]",m_phiBins,m_phiMin,m_phiMax,nybins,ymin,ymax); + } + + /////////////////// + // Fill Sampling // + /////////////////// + + template + void fillSampling(H* hptr,CaloSampling::CaloSample csamp) { + int isamp(std::min(static_cast(csamp),static_cast(CaloRec::Lookup::getSamplingId("MINIFCAL0")))); + hptr->Fill(1.*isamp); + } + + template + void fillSampling(H* hptr,CaloSampling::CaloSample csamp,double yval) { + int isamp(std::min(static_cast(csamp),static_cast(CaloRec::Lookup::getSamplingId("MINIFCAL0")))); + hptr->Fill(1.*isamp,yval); + } + + template + void fillSampling(H* hptr,CaloSampling::CaloSample csamp,double yval,double zval) { + int isamp(std::min(static_cast(csamp),static_cast(CaloRec::Lookup::getSamplingId("MINIFCAL0")))); + hptr->Fill(1.*isamp,yval,zval); + } +}; // CaloTopoClusterFromTowerMonitor + +inline bool CaloTopoClusterFromTowerMonitor::isInHotspot(const xAOD::CaloCluster& ptow) const +{ return m_doHotspot && ( ptow.eta() >= m_hsEtaMin && ptow.eta() < m_hsEtaMax ) && ( ptow.phi() >= m_hsPhiMin && ptow.phi() < m_hsPhiMax ); } + +inline bool CaloTopoClusterFromTowerMonitor::setAxisTitle(TH1* hptr,const std::string& title,const std::string& axis) { + if ( axis == "x" ) { hptr->GetXaxis()->SetTitle(title.c_str()); return true; } + if ( axis == "y" ) { hptr->GetYaxis()->SetTitle(title.c_str()); return true; } + if ( axis == "z" ) { hptr->GetZaxis()->SetTitle(title.c_str()); return true; } + return false; +} +#endif + + diff --git a/Calorimeter/CaloRec/CaloRec/CaloTopoClusterTowerMerger.h b/Calorimeter/CaloRec/CaloRec/CaloTopoClusterTowerMerger.h new file mode 100644 index 0000000000000000000000000000000000000000..792ae2cf9f7fb14c16004a9d1b33c8f687c24206 --- /dev/null +++ b/Calorimeter/CaloRec/CaloRec/CaloTopoClusterTowerMerger.h @@ -0,0 +1,86 @@ +// -*- c++ -*- + +/* + Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +*/ + + +#ifndef CALOREC_CALOTOPOCLUSTERTOWERMERGER_H +#define CALOREC_CALOTOPOCLUSTERTOWERMERGER_H + +#include "StoreGate/ReadHandleKey.h" +#include "StoreGate/WriteHandleKey.h" + +#include "AthenaBaseComps/AthAlgorithm.h" + +#include "xAODCaloEvent/CaloCluster.h" +#include "xAODCaloEvent/CaloClusterContainer.h" + +#include +//#include +#include +#include + +class CaloTopoClusterTowerMerger : public AthAlgorithm +{ +public: + ///@brief Algorithm constructor + CaloTopoClusterTowerMerger(const std::string& name,ISvcLocator* pSvcLocator); + ///@brief Baseclass destructor + virtual ~CaloTopoClusterTowerMerger(); + + ///@name Algorithm interface + ///@{ + virtual StatusCode initialize(); ///< Initialization sets up read and write handle keys + virtual StatusCode execute(); ///< Execution merges the container contents + ///@} + +private: + + ///@name Internally used types + ///@{ + typedef SG::ReadHandleKey rhandlekey_t; ///< Input data handle key type + typedef SG::WriteHandleKey whandlekey_t; ///< Output data handle key type + typedef SG::ReadHandle rhandle_t; ///< Input data handle type + typedef SG::WriteHandle whandle_t; ///< Output data handle type + ///@} + + ///@name Algorithm properties + ///@{ + rhandlekey_t m_clusterContainerKey; ///< Input topo-cluster container + rhandlekey_t m_towerContainerKey; ///< Input topo-tower container + whandlekey_t m_topoSignalContainerKey; ///< Output merged (view) container + double m_clusterRange; ///< Rapidity range for topo-clusters + // bool m_copyMoments; ///< Explicitely copy cluster moments + ///@} + + ///@name Helpers + ///@{ + bool makeDeepCopy(const xAOD::CaloCluster& rClus,xAOD::CaloClusterContainer* pClusCont) const; ///< Attaches a deep copy to container, returns @c true if successful. + // bool fillMoments(const xAOD::CaloCluster& rClus); ///< Copies list of filled moments into lookup + bool clusterFilter(const xAOD::CaloCluster& rClus) const; ///< Filter topo-cluster + bool towerFilter(const xAOD::CaloCluster& rTowr) const; ///< Filter topo-tower + StatusCode addContainerWriteHandle(whandle_t& signalHandle); ///< Add a write handle for a container (in CaloClusterStoreHelper from r21.9) + ///@} + + ///@name Moment lookup + ///@{ + // static std::vector > m_momentMap; ///< Map moment types to human readable names + // static std::vector m_momentList; ///< List of used moments for given tower collection + ///@} +}; + + +inline bool CaloTopoClusterTowerMerger::clusterFilter(const xAOD::CaloCluster& rClus) const { return std::abs(rClus.eta()) <= m_clusterRange; } +inline bool CaloTopoClusterTowerMerger::towerFilter(const xAOD::CaloCluster& /*rTowr*/) const { return true; } + +///@class CaloTopoClusterTowerMerger +/// +/// This algorithm merges objects from two @c xAOD::CaloClusterContainer. In the context of the +/// mixed topo-tower/topo-cluster output, the objects in the topo-cluster container are taken +/// up to a client-defined (symmetric) rapidity range. The rest of the phase space is then filled +/// with the objects from the topo-tower container. It is assumed that the overlap resolution +/// is performed when the topo-tower container is filled. The mixed object container contains deep copies +/// of the objects in the input containers. +/// +#endif diff --git a/Calorimeter/CaloRec/CaloRec/CaloTowerGeometrySvc.h b/Calorimeter/CaloRec/CaloRec/CaloTowerGeometrySvc.h new file mode 100644 index 0000000000000000000000000000000000000000..6d2f445548ea860896ee62dad476ca57bcd72b70 --- /dev/null +++ b/Calorimeter/CaloRec/CaloRec/CaloTowerGeometrySvc.h @@ -0,0 +1,414 @@ +// -*- c++ -*- +#ifndef CALOREC_CALOTOWERGEOMETRYSVC_H +#define CALOREC_CALOTOWERGEOMETRYSVC_H + +/* Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration */ + +#include "AthenaBaseComps/AthService.h" + +#include "Identifier/IdentifierHash.h" + +#include "CaloGeoHelpers/CaloPhiRange.h" + +#include "CaloDetDescr/CaloDetDescrElement.h" +#include "CaloDetDescr/CaloDetDescrManager.h" + +#include "CaloEvent/CaloCell.h" + +#include +#include +#include +#include +#include +#include + +///@name Forward declarations +///@{ +///template class SvcFactory; +///@} + +///@brief Interface accessor for @c CaloTowerGeometrySvc +static const InterfaceID IID_CaloTowerGeometrySvc( "CaloTowerGeometrySvc", 1, 0 ); + +class CaloTowerGeometrySvc : public AthService +{ +protected: + + ///@brief To allow interface queries by the service factory (?) + friend class SvcFactory; + +public: + + ///@name Gaudi interfaces and implementations + ///@{ + static const InterfaceID& interfaceID() { return IID_CaloTowerGeometrySvc; } ///< Interface indentifier needed by Gaudi + virtual StatusCode queryInterface(const InterfaceID& riid, void** ppcInterface); ///< Interface query with fallbacks + ///@} + + ///@{ + typedef std::size_t uint_t; ///< Type for unsigned integer + typedef IdentifierHash::value_type index_t; ///< Type for scalar (pseudorapidity,azimuth) index (is an unsigned int type) + typedef std::tuple towerindex_t; ///< Type for composite (tower) index + typedef std::tuple element_t; ///< Type storing tower index and geometrical weight + typedef std::vector elementvector_t; ///< Type for list of elements holding tower index and list of weights + typedef std::vector elementmap_t; ///< Type for lists of lists of elements (lookup table type) + ///@} + + ///@brief Standard constructor + /// + /// Constructs towergrid as defined by properties. + /// + ///@param[in] name name of service (implementation of @c Gaudi::Service ) + ///@param[in] pSvc pointer to service locator (from configuration framework) + CaloTowerGeometrySvc(const std::string& name,ISvcLocator* pSvc); + ///@brief Base class destructor + virtual ~CaloTowerGeometrySvc() { }; + + ///@name Implementation of (Ath)Service interfaces + ///@{ + virtual StatusCode initialize() override; ///< Initialize service + virtual StatusCode finalize() override; ///< Finalize service + // virtual StatusCode queryInterface(const InterfaceID& iid,void** ppvInterface); ///< Need to fix compiler issue (FIXME) + ///@} + + // --- Full documentation of this block after end of class definition! + ///@name Retrieve towers for cells + ///@{ + StatusCode access(const CaloCell* pCell,std::vector& towerIdx,std::vector& towerWghts) const; + StatusCode access(IdentifierHash cellHash,std::vector& towerIdx,std::vector& towerWghts) const; + elementvector_t getTowers(const CaloCell* pCell) const; + elementvector_t getTowers(IdentifierHash cellHash) const; + ///@} + + ///@name Tower bin descriptors and other size information + ///@{ + uint_t maxCellHash() const; ///< Maximum cell hash value + uint_t totalNumberCells() const; ///< Total number of cells + uint_t etaBins() const; ///< Number of pseudorapidity bins + double etaMin() const; ///< Lower boundary of pseudorapidity range + double etaMax() const; ///< Upper boundary of pseudorapidity range + double etaWidth() const; ///< Width of pseudorapidity bin + uint_t phiBins() const; ///< Number of azimuth bins + double phiMin() const; ///< Lower boundary of azimuth + double phiMax() const; ///< Upper boundary of azimuth + double phiWidth() const; ///< Width of azimuth bin + uint_t towerBins() const; ///< Total number of towers + double towerArea() const; ///< Area of individual tower + ///@} + + ///@name Tower index calculators, translaters, manipulators and converters + ///@{ + index_t etaIndex(const CaloCell* pCell) const; ///< Get tower @f$ \eta @f$ bin index for a calorimeter cell referenced by a pointer + index_t etaIndex(IdentifierHash cellHash) const; ///< Get tower @f$ \eta @f$ bin index for a calorimeter cell referenced by its hash identifier + index_t etaIndex(double eta) const; ///< Get tower @f$ \eta @f$ bin index for a given value of @f$ \eta @f$ + index_t etaIndexFromTowerIndex(index_t towerIdx) const; ///< Get tower @f$ \eta @f$ bin index for a given global tower index + index_t phiIndex(const CaloCell* pCell) const; ///< Get tower @f$ \phi @f$ bin index for a calorimeter cell referenced by a pointer + index_t phiIndex(IdentifierHash cellHash) const; ///< Get tower @f$ \phi @f$ bin index for a calorimeter cell referenced by its hash identifier + index_t phiIndex(double phi) const; ///< Get tower @f$ \phi @f$ bin index for a given value of @f$ \phi @f$ + index_t phiIndexFromTowerIndex(index_t towerIdx) const; ///< Get tower @f$ \phi @f$ bin index for a given global tower index + index_t towerIndex(const CaloCell* pCell) const; ///< Get global tower index for a calorimeter cell referenced by a pointer + index_t towerIndex(IdentifierHash cellHash) const; ///< Get global tower index for a calorimeter cell referenced by its hash identifier + index_t towerIndex(double eta,double phi) const; ///< Get global tower index for a pair of @f$ (\eta,\phi) @f$ values + index_t towerIndex(index_t etaIdx,index_t phiIdx) const; ///< Get global tower index for a pair of @f$ (\eta,\phi) @f$ indices + index_t towerIndex(const element_t& elm) const; ///< Get global tower index from payload data + index_t invalidIndex() const; ///< Returns value of invalid index + bool isInvalidIndex(index_t idx) const; ///< Returns @c true if argument is equal to the value provided by @c invalidIndex() + bool isInvalidIndex(index_t idx,index_t maxIdx) const; ///< Returns @c true if first argument is equal to the value provided by @c invalidIndex() or if first argument is not smaller than second argument + bool isInvalidEtaIndex(index_t idx) const; ///< Returns @c true if argument is not a valid pseudorapidity index + bool isInvalidPhiIndex(index_t idx) const; ///< Returns @c true if argumant is not a valid azimuth index + bool isInvalidTowerIndex(index_t idx) const; ///< Returns @c true if argument is not a valid tower index + ///@} + + ///@name Variable generators using tower indices + ///@{ + double towerEtaLocal(index_t etaIndex) const; ///< Return pseudorapdity from local index (bin center) + double towerPhiLocal(index_t phiIndex) const; ///< Return azimuth from local index (bin center) + double towerEta(index_t towerIndex) const; ///< Return pseudorapidity from global tower index (bin center) + double towerPhi(index_t towerIndex) const; ///< Return azimuth from global tower index (bin center) + double invalidValue() const; ///< Return invalid value + bool isInvalidValue(double val) const; ///< Return @c true if given value is invalid + ///@} + + ///@name Helper functions + ///@{ + double cellWeight(const element_t& elm) const; ///< Retrieve cell signal weight from lookup table entry + double cellWeight(IdentifierHash cellHash,index_t towerIdx) const; ///< Retrieve cell signal weight from cell identifier and tower index + double cellWeight(const CaloCell* pCell, index_t towerIdx) const; ///< Retrieve cell signal weight from pointer to cell object and tower index + ///@} + + ///@name Access to storage + ///@{ + elementmap_t::const_iterator begin() const; ///< Iterator points to first entry in internal look-up table (only @c const access!) + elementmap_t::const_iterator end() const; ///< Iterator marks end of internal look-up table (only @c const access) + size_t size() const; ///< Size of internal look-up table + bool empty() const; ///< Internal look-up table is empty if @c true + ///@} + +private: + + ///@name Helpers + ///@{ + StatusCode f_setupSvc(); ///< Internally used function setting up other services needed by this service + StatusCode f_setupTowerGrid(); ///< Internally used function setting up the lookup store + StatusCode f_setupTowerGridFCal(const CaloDetDescrElement* pCaloDDE,std::ofstream& logger); ///< Internally used function mapping an FCal cell onto the tower grid + StatusCode f_setupTowerGridProj(const CaloDetDescrElement* pCaloDDE,std::ofstream& logger); ///< Internally used function mapping a projective cell onto the tower grid + double f_assign(IdentifierHash cellHash,index_t towerIdx,double wgt); ///< Internally used function assigning tower to cell with update of weight if same tower is already assigned + ///@} + + ///@name Access to detector store and other services and stores + ///@{ + const CaloDetDescrManager* m_caloDDM; ///< Pointer to calorimeter detector description + std::string m_logFileName; ///< Name of log file + ///@} + +protected: + + ///@name Internal stores and derived parameters + ///@{ + elementmap_t m_towerLookup; ///< Cell-to-tower mapping lookup store + double m_towerEtaWidth; ///< Width of tower bin in pseudorapidity + double m_towerPhiWidth; ///< Width of tower bin in azimuth + double m_towerArea; ///< Area of individual tower + uint_t m_towerBins; ///< Maximum number of towers + uint_t m_maxCellHash; ///< Maximum cell hash value + uint_t m_numberOfCells; ///< Total number of cells + ///@} + + ///@name Properties + ///@{ + ///@brief Internally stored tower grid descriptors + uint_t m_towerEtaBins; ///< Number of @f$ \eta @f$ bins + double m_towerEtaMin; ///< Lower boundary @f$ \eta_{\rm min} @f$ + double m_towerEtaMax; ///< Upper boundary @f$ \eta_{\rm max} @f$ + bool m_adjustEta; ///< Adjust FCal cells to eta boundary (default @c true ) + uint_t m_towerPhiBins; ///< Number of @f$ \phi @f$ bins + double m_towerPhiMin; ///< Lower boundary @f$ \phi_{\rm min} @f$ + double m_towerPhiMax; ///< Upper boundary @f$ \phi_{\rm max} @f$ + double m_fcal1Xslice; ///< Number of x slices for cells in FCal1 + double m_fcal1Yslice; ///< Number of y slices for cells in FCal1 + double m_fcal2Xslice; ///< Number of x slices for cells in FCal2 + double m_fcal2Yslice; ///< Number of y slices for cells in FCal2 + double m_fcal3Xslice; ///< Number of x slices for cells in FCal3 + double m_fcal3Yslice; ///< Number of y slices for cells in FCal3 + ///@} + + ///@name Process flags, helpers and numerical constants + ///@{ + static index_t m_invalidIndex; ///< Invalid index indicator + static double m_invalidValue; ///< Return value for out-of-range indices andother invalid conversions to a physical quantity + const CaloDetDescrManager* f_caloDDM() const; ///< Pointer to calorimeter detector description manager + const CaloDetDescrElement* f_caloDDE(const CaloCell* pCell) const; ///< Retrieve calorimeter detector description element for a cell object referenced by a pointer + const CaloDetDescrElement* f_caloDDE(IdentifierHash cellHash) const; ///< Retrieve calorimeter detector description element for a given cell hash identifier + double f_cellEta(const CaloCell* pCell) const; ///< Retrieve calorimeter cell pseudorapidity for a cell object referenced by a pointer + double f_cellEta(IdentifierHash cellHash) const; ///< Retrieve calorimeter cell pseudorapidity for a given cell hash identifier + double f_cellPhi(const CaloCell* pCell) const; ///< Retrieve calorimeter cell azimuth for a cell object referenced by a pointer + double f_cellPhi(IdentifierHash cellHash) const; ///< Retrieve calorimeter cell azimuth for a given cell hash identifier + ///@} + + ///@name Stores + ///@{ + std::array m_ndxFCal; ///< Stores number of fragments along x for each FCal module + std::array m_ndyFCal; ///< Stores number of fragments along y for each FCal module + std::array m_wgtFCal; ///< Stores geometrical weights + ///@} +}; + +//-------------------------------------------------// +// Documentation for grouped methods and functions // +// (removed from before/after method for better // +// formatting by doxygen in html). // +//-------------------------------------------------// + +///@fn StatusCode CaloTowerGeometrySvc::access(const CaloCell* pCell,std::vector& towerIdx,std::vector& towerWghts) const +/// +///@brief Retrieve the list of towers associated with a calorimeter cell referenced by a pointer +/// +/// The tower indices and weights are returned in two index-parallel vectors. +/// Previous content of these two vectors is removed if this method finds towers for the cell. +/// +///@return Returns @c StatusCode::SUCCESS if list of towers found, else @s StatusCode::FAILURE. +/// +///@param[in] pCell pointer to non-modifiable @c CaloCell object. +///@param[in] towerIdx reference to modifiable vector of indices (payload type @c index_t ); vector is filled if cell is successfully mapped. +///@param[in] towerWghts reference to modifiable vector of weights (payload type @c double ); vector is filled if cell is successfully mapped. + +///@fn StatusCode CaloTowerGeometrySvc::access(IdentifierHash cellHash,std::vector& towerIdx,std::vector& towerWghts) const; +/// +///@brief Retrieve the list of towers associated with a calorimeter cell referenced its hash identifier +/// +/// The tower indices and weights are returned in two index-parallel vectors. +/// Previous content of these two vectors is removed if this method finds towers for the cell. +/// +///@return Returns @c StatusCode::SUCCESS if list of towers found, else @s StatusCode::FAILURE. +/// +///@param[in] cellHash hash identifier referencing a calorimeter cell. +///@param[in] towerIdx reference to modifiable vector of indices (payload type @c index_t ); vector is filled if cell is successfully mapped. +///@param[in] towerWghts reference to modifiable vector of weights (payload type @c double ); vector is filled if cell is successfully mapped. + +///@fn CaloTowerGeometrySvc::elementvector_t CaloTowerGeometrySvc::getTowers(const CaloCell* pCell) const; +/// +///@brief Retrieve the list of towers associated with a calorimeter cell referenced by a pointer +/// +///@return Returns a vector of (index,weight) pairs as a @c elementvector_t container. The container is empty +/// if the cell does not have any overlap with a tower. +/// +///@param[in] pCell pointer to non-modifiable @c CaloCell object. + +///@fn CaloTowerGeometrySvc::elementvector_t CaloTowerGeometrySvc::getTowers(IdentifierHash cellHash) const; +/// +///@brief Retrieve the list of towers associated with a calorimeter cell referenced by its hash identifier +/// +///@return Returns a vector of (index,weight) pairs as a @c elementvector_t container. The container is empty +/// if the cell does not have any overlap with a tower. +/// +///@param[in] cellHash hash identifier referencing a calorimeter cell. + +//---------------------// +// Class documentation // +//---------------------// + +/// @class CaloTowerGeometrySvc +/// +/// @brief Tower geometry store and description provider +/// +/// This service sets up a lookup table storing the geometrical area overlap fraction of a calorimeter cell +/// with towers in a given grid. This lookup table is set up at instantiation of the service. It can only be +/// defined at that time. The default setup is a @f$ \Delta\eta\times\Delta\phi = 0.1 \times \pi/32 @f$ grid. +/// Any regular grid can be constructed. The grid definition can be provided as property. +/// +/// The cell-to-tower information is stored internally as a (random access) lookup table. For a given cell, +/// the hash index is used to retrieve a list of towers this cell overlaps with, and the overlap paramater +/// (area fraction used as a geometrical weight). This indices and geometrical weights are represented by +/// a list of pairs of @c int and @c double numbers. Each cell can potential overlap with more than one +/// tower. A more detailed description of towers and the geometrical overlap is available on the +/// calorimeter tower project page. +/// +/// The lookup table is implemented for random access and using the cell hash identifier to retrieve the +/// requested list of twoer indices and weights. Several retrieval mechanisms are supported (see documentation +/// of the corresponding methods). +/// +/// To map the azimuth of a cell to a tower, @f$ -\pi < \phi < \pi @f$ is used (ATLAS standard). For +/// consistency, all @f$ \phi @f$ values are mapped into this range. +/// +/// The service inherits from @c AthService and thus from the @c Service base class in Gaudi. The managed tower grid +/// is defined by service properties, with the following naming convention: +/// - pseudorapidity range +/// - number of bins TowerEtaBins (default 100) +/// - lower boundary of pseudorapidity range TowerEtaMin (default -5.0) +/// - upper boundary of pseudorapidity range TowerEtaMax (default 5.0) +/// - azimuth range +/// - number of bins TowerPhiBins (default 64) +/// - lower boundary of azimuthal range TowerPhiMin (default -π) +/// - upper boundary of azimuthal range TowerPhiMax (default π) +/// +/// Addtional properties of this service define the granularity of the cell splitting in the ATLAS FCal. This +/// is used to map the FCal readout cells (rectangular slabs) onto the tower grid and calculate the geometrical +/// (area) overlap fraction, which is used to distribute the cell energy to the towers. +/// - horizontal FCal cell splitting (along @a x axis) +/// - number of @a x slices in FCal1 FCal1NSlicesX (default 4) +/// - number of @a x slices in FCal2 FCal2NSlicesX (default 4) +/// - number of @a x slices in FCal3 FCal3NSlicesX (default 6) +/// - vertical FCal cell splitting (along @a y axis) +/// - number of @a y slices in FCal1 FCal1NSlicesY (default 4) +/// - number of @a y slices in FCal2 FCal2NSlicesY (default 6) +/// - number of @a y slices in FCal3 FCal3NSlicesY (default 6) +/// +/// @warning It is recommended to @b not change the parameters for the FCal cell slicing. This configuration option is provided for expert use for R & D purposes only. +/// +/// @todo Allow regional grids (varying segmentation as function of @f$ \eta @f$ . This requires additional interfaces (or interface changes) and +/// and modifications of the index construction. +/// +/// @author Peter Loch +/// + +//------------------// +// Inline Functions // +//------------------// + +//---------------------------// +// Control and configuration // +//---------------------------// +inline CaloTowerGeometrySvc::index_t CaloTowerGeometrySvc::invalidIndex() const { return m_invalidIndex; } +inline StatusCode CaloTowerGeometrySvc::f_setupSvc() { + m_caloDDM = CaloDetDescrManager::instance(); + return f_caloDDM() != 0 ? StatusCode::SUCCESS : StatusCode::FAILURE; +} + +//------------------------------------// +// Public access to tower descriptors // +//------------------------------------// +inline CaloTowerGeometrySvc::uint_t CaloTowerGeometrySvc::maxCellHash() const { return m_maxCellHash; } +inline CaloTowerGeometrySvc::uint_t CaloTowerGeometrySvc::totalNumberCells() const { return m_numberOfCells; } + +inline CaloTowerGeometrySvc::uint_t CaloTowerGeometrySvc::etaBins() const { return m_towerEtaBins; } +inline double CaloTowerGeometrySvc::etaMin() const { return m_towerEtaMin; } +inline double CaloTowerGeometrySvc::etaMax() const { return m_towerEtaMax; } +inline double CaloTowerGeometrySvc::etaWidth() const { return m_towerEtaWidth; } + +inline CaloTowerGeometrySvc::uint_t CaloTowerGeometrySvc::phiBins() const { return m_towerPhiBins; } +inline double CaloTowerGeometrySvc::phiMin() const { return m_towerPhiMin; } +inline double CaloTowerGeometrySvc::phiMax() const { return m_towerPhiMax; } +inline double CaloTowerGeometrySvc::phiWidth() const { return m_towerPhiWidth; } + +inline CaloTowerGeometrySvc::uint_t CaloTowerGeometrySvc::towerBins() const { return m_towerBins; } +inline double CaloTowerGeometrySvc::towerArea() const { return m_towerArea; } + +//----------------// +// Index checking // +//----------------// +inline bool CaloTowerGeometrySvc::isInvalidIndex(index_t idx) const { return idx == invalidIndex(); } +inline bool CaloTowerGeometrySvc::isInvalidIndex(index_t idx,index_t maxIdx) const { return idx == invalidIndex() || idx >= maxIdx; } +inline bool CaloTowerGeometrySvc::isInvalidEtaIndex(index_t idx) const { return isInvalidIndex(idx,m_towerEtaBins); } +inline bool CaloTowerGeometrySvc::isInvalidPhiIndex(index_t idx) const { return isInvalidIndex(idx,m_towerPhiBins); } +inline bool CaloTowerGeometrySvc::isInvalidTowerIndex(index_t idx) const { return isInvalidIndex(idx,m_towerBins); } + +//------------------------------// +// Index retrieval/construction // +//------------------------------// +inline CaloTowerGeometrySvc::index_t CaloTowerGeometrySvc::etaIndex(const CaloCell* pCell) const { return etaIndex(pCell->eta()); } +inline CaloTowerGeometrySvc::index_t CaloTowerGeometrySvc::phiIndex(const CaloCell* pCell) const { return phiIndex(pCell->phi()); } + +inline CaloTowerGeometrySvc::index_t CaloTowerGeometrySvc::towerIndex(IdentifierHash cellHash) const { return towerIndex(etaIndex(cellHash),phiIndex(cellHash)); } +inline CaloTowerGeometrySvc::index_t CaloTowerGeometrySvc::towerIndex(double eta,double phi) const { return towerIndex(etaIndex(eta),phiIndex(phi)); } +inline CaloTowerGeometrySvc::index_t CaloTowerGeometrySvc::towerIndex(index_t etaIdx,index_t phiIdx) const { return !isInvalidEtaIndex(etaIdx) && !isInvalidPhiIndex(phiIdx) ? phiIdx+etaIdx*m_towerPhiBins : invalidIndex(); } +inline CaloTowerGeometrySvc::index_t CaloTowerGeometrySvc::towerIndex(const CaloCell* pCell) const { return towerIndex(pCell->eta(),pCell->phi()); } +inline CaloTowerGeometrySvc::index_t CaloTowerGeometrySvc::towerIndex(const element_t& elm) const { return std::get<0>(elm); } + +inline CaloTowerGeometrySvc::index_t CaloTowerGeometrySvc::etaIndexFromTowerIndex(index_t towerIdx) const { return (index_t)(towerIdx/phiBins()); } +inline CaloTowerGeometrySvc::index_t CaloTowerGeometrySvc::phiIndexFromTowerIndex(index_t towerIdx) const { return (index_t)(towerIdx%phiBins()); } + +//-----------------------------// +// Access to tower description // +//-----------------------------// + +inline double CaloTowerGeometrySvc::invalidValue() const { return m_invalidValue; } +inline bool CaloTowerGeometrySvc::isInvalidValue(double val) const { return val == invalidValue(); } +inline double CaloTowerGeometrySvc::towerEtaLocal(index_t etaIdx) const { return !isInvalidEtaIndex(etaIdx) ? etaMin()+(static_cast(etaIdx)+0.5)*etaWidth() : invalidValue(); } +inline double CaloTowerGeometrySvc::towerPhiLocal(index_t phiIdx) const { return !isInvalidPhiIndex(phiIdx) ? phiMin()+(static_cast(phiIdx)+0.5)*phiWidth() : invalidValue(); } +inline double CaloTowerGeometrySvc::towerEta(index_t towerIdx) const { return towerEtaLocal(etaIndexFromTowerIndex(towerIdx)); } +inline double CaloTowerGeometrySvc::towerPhi(index_t towerIdx) const { return towerPhiLocal(phiIndexFromTowerIndex(towerIdx)); } + +inline CaloTowerGeometrySvc::elementmap_t::const_iterator CaloTowerGeometrySvc::begin() const { return m_towerLookup.begin(); } +inline CaloTowerGeometrySvc::elementmap_t::const_iterator CaloTowerGeometrySvc::end() const { return m_towerLookup.end(); } +inline size_t CaloTowerGeometrySvc::size() const { return m_towerLookup.size(); } +inline bool CaloTowerGeometrySvc::empty() const { return m_towerLookup.empty(); } + +//-------------------// +// Other data access // +//-------------------// +inline double CaloTowerGeometrySvc::cellWeight(const CaloCell* pCell,index_t towerIdx) const { return cellWeight(pCell->caloDDE()->calo_hash(),towerIdx); } +inline double CaloTowerGeometrySvc::cellWeight(const element_t& elm) const { return std::get<1>(elm); } + + +//----------------------------------// +// Internal functions and accessors // +//----------------------------------// +inline const CaloDetDescrManager* CaloTowerGeometrySvc::f_caloDDM() const { return m_caloDDM; } +inline const CaloDetDescrElement* CaloTowerGeometrySvc::f_caloDDE(const CaloCell* pCell) const { return pCell->caloDDE(); } +inline const CaloDetDescrElement* CaloTowerGeometrySvc::f_caloDDE(IdentifierHash cellHash) const { return f_caloDDM()->get_element(cellHash); } + +inline double CaloTowerGeometrySvc::f_cellEta(IdentifierHash cellHash) const { return f_caloDDE(cellHash)->eta(); } +inline double CaloTowerGeometrySvc::f_cellEta(const CaloCell* pCell) const { return pCell->eta(); } +inline double CaloTowerGeometrySvc::f_cellPhi(IdentifierHash cellHash) const { return CaloPhiRange::fix(f_caloDDE(cellHash)->phi()); } +inline double CaloTowerGeometrySvc::f_cellPhi(const CaloCell* pCell) const { return pCell->phi(); } +#endif diff --git a/Calorimeter/CaloRec/python/CaloRecFlags.py b/Calorimeter/CaloRec/python/CaloRecFlags.py index bb6f10ecd1022ec247009a4d03861282596929ee..5f94d15e10d139c3611e4e85cd4d0becd458ec86 100644 --- a/Calorimeter/CaloRec/python/CaloRecFlags.py +++ b/Calorimeter/CaloRec/python/CaloRecFlags.py @@ -87,7 +87,6 @@ class doCaloTowerFromCluster(CaloRecFlagsJobProperty): allowedTypes=['bool'] StoredValue=False - class doCaloTopoTower(CaloRecFlagsJobProperty): """ switch noise suppressed towers based on standard tower + topo clusters """ @@ -95,6 +94,20 @@ class doCaloTopoTower(CaloRecFlagsJobProperty): allowedTypes=['bool'] StoredValue=False +class doCaloTopoSignal(CaloRecFlagsJobProperty): + """ produce mixed topo-cluster and topo-tower container + """ + statusOn=True + allowedTypes=['bool'] + storedValue=False + +class doExtendedClusterMoments(CaloRecFlagsJobProperty): + """ add more cluster moments for R&D + """ + statusOn=True + allowedTypes=['bool'] + storedValue=False + class doTileMuId(CaloRecFlagsJobProperty): """ switch for TileMuId - muon finding algorighm """ @@ -170,7 +183,7 @@ jobproperties.add_Container(CaloRecFlags) # I want always the following flags in the Rec container -_list_Calo=[Enabled,doCaloTopoCluster,doEmCluster,doCaloEMTopoCluster,emTopoClusterThreshold,doCaloCluster,doCaloTopoTower,doTileMuId,doTileCellCorrection,doLArAffectedRegion,doLArAutoConfiguration,doLArNoisyRO,doEMDigits,doFillMBTSBackgroundBit,doLArNoiseBurstVeto,clusterCellGetterName,doCaloTowerFromCells,doCaloTowerFromCluster] +_list_Calo=[Enabled,doCaloTopoCluster,doEmCluster,doCaloEMTopoCluster,emTopoClusterThreshold,doCaloCluster,doCaloTopoTower,doCaloTopoSignal,doExtendedClusterMoments,doTileMuId,doTileCellCorrection,doLArAffectedRegion,doLArAutoConfiguration,doLArNoisyRO,doEMDigits,doFillMBTSBackgroundBit,doLArNoiseBurstVeto,clusterCellGetterName,doCaloTowerFromCells,doCaloTowerFromCluster] for j in _list_Calo: jobproperties.CaloRecFlags.add_JobProperty(j) del _list_Calo diff --git a/Calorimeter/CaloRec/python/MakeClustersFromTowers.py b/Calorimeter/CaloRec/python/MakeClustersFromTowers.py index 6bc75b051bf3c1f72f4a94ab13a76ad389879347..afa1f7e3cc89f884175c9d573e5b8d8c840c280b 100644 --- a/Calorimeter/CaloRec/python/MakeClustersFromTowers.py +++ b/Calorimeter/CaloRec/python/MakeClustersFromTowers.py @@ -1,164 +1,261 @@ # Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -from CaloRec.CaloRecConf import CaloTowerxAODFromCells -from CaloRec.CaloRecConf import CaloTopoClusterFromTowerMaker +from AthenaCommon.AppMgr import ServiceMgr as svcMgr -from AthenaCommon.Logging import logging +from CaloRec.CaloRecConf import CaloTopoClusterFromTowerMaker +from CaloRec.CaloRecConf import CaloTowerGeometrySvc +from AthenaCommon.Logging import logging import AthenaCommon.Constants as Lvl -def ClustersFromTowersDict(clusterBuilderName='TowerFromClusterMaker', - towerBuilderAlgo=CaloTowerxAODFromCells('CmbTowerAlgo'), - orderedClusterByPt=False, - doLCW=False, #### don't apply LCW by default - caloTopoClusterKey='NONE', #### all cells is default - caloCellClusterWeightKey='NONE'): #### no cell weight data object needed +#################################### +## Tower configuration dictionary ## +#################################### + +def ClustersFromTowersDict(clusterBuilderName = 'TowerFromClusterTool', + towerGeometrySvc = CaloTowerGeometrySvc('CaloTowerGeometrySvc'), + cellContainerKey = 'AllCalo', + buildTopoTowers = True, + topoClusterContainerKey = 'CaloTopoCluster', + cellClusterWeightKey = 'CaloCellClusterWeightKey', + orderClusterByPt = False, + applyCellEnergyThreshold = False, + doCellIndexCheck = False, + cellEnergyThreshold = 0., + applyLCW = False, + buildCombinedSignal = False, + clusterRange = 5.): ''' Configuration dictionary for tower-to-cluster converter ''' - configDict = { 'ClusterBuilderName' : clusterBuilderName, - 'CaloTowerBuilder' : towerBuilderAlgo, - 'OrderClusterByPt' : orderedClusterByPt, - 'CaloTopoClusterContainerKey' : caloTopoClusterKey, - 'CellClusterWeightKey' : caloCellClusterWeightKey, - 'ApplyLCW' : doLCW + configDict = { 'ClusterBuilderName' : clusterBuilderName, ### name of the tower builder tool + 'CaloTowerGeometrySvc' : towerGeometrySvc, ### tower geometry provider + 'CaloCellContainerKey' : cellContainerKey, ### (input) cell container key + 'CaloTopoClusterContainerKey' : topoClusterContainerKey, ### (input) topo-cluster container key + 'CellClusterWeightKey' : cellClusterWeightKey, ### (output) cell weights in topo-clusters + 'BuildTopoTowers' : buildTopoTowers, ### (control) form topo-towers + 'OrderClusterByPt' : orderClusterByPt, ### (control) order final clusters by Pt + 'ApplyCellEnergyThreshold' : applyCellEnergyThreshold, ### (control) apply energy thresholds to cells + 'CellEnergyThreshold' : cellEnergyThreshold, ### (control) value of energy threshold + 'PrepareLCW' : applyLCW, ### (control) prepare (and apply) LCW + 'DoCellIndexCheck' : doCellIndexCheck, ### (control) check cell hash indices + 'BuildCombinedTopoSignal' : buildCombinedSignal, ### (control) build combined topo-cluster/topo-tower container + 'TopoClusterRange' : clusterRange, ### (control) range for topo-cluster in combined mode } return configDict -def MakeClustersFromTowers(clusterMakerName='CaloClusterMaker',clusterContainerKey='TowerTopoCluster',configDict=ClustersFromTowersDict(),applyEnergyThreshold=False,debugOn=False): - ''' This function generates an instance of a cluster algorithm producting clusters trom towers with or without moments +################################### +## Tower algorithm configuration ## +################################### + +def MakeClustersFromTowers(towerMakerName = 'CaloTowerBuilderAlg', ### name of tower builder algorithm + towerContainerKey = 'CaloTowerTopoCluster', ### output container key + configDict = ClustersFromTowersDict(), ### tower builder tool configuration + debugOn = False): + ''' This function generates an instance of a cluster algorithm configuration producting clusters trom towers with or without moments. ''' - # collect inputs mlog = logging.getLogger('MakeClustersFromTowers.py:: ') - mlog.info('ClusterMakerName = "'+clusterMakerName+'"') - mlog.info('ClusterContainerKey = <'+clusterContainerKey+'>') - mlog.info('Converter parameters: ',configDict) - - # configure cluster builder - cnvname = configDict['ClusterBuilderName'] - twralgo = configDict['CaloTowerBuilder'] - towerkey = twralgo.CaloTowerContainer - mlog.info('(input) CaloTowerContainer <'+towerkey+'>') - cellkey = 'AllCalo' ### twralgo.InputCellContainer --> this does not work, why? - mlog.info('(input) CaloCellContainer <'+cellkey+'>') - ptorder = configDict['OrderClusterByPt'] - tcluskey = configDict['CaloTopoClusterContainerKey'] - tcwgtkey = configDict['CellClusterWeightKey'] - #### buildtt = ( tcluskey != 'NONE' and tcwgtkey != 'NONE' ) - buildtt = configDict['ApplyLCW'] - ''' Configuration module for the tower converter + mlog.info('TowerMakerName = "'+towerMakerName+'"') + mlog.info('TowerContainerKey = <'+towerContainerKey+'>') + + ######################################## + ## Configuring the tower builder tool ## + ######################################## + + ''' collect properties from dictionary and set correct dependencies + ''' + mlog.info('Converter properties: ',configDict) + excludedKeys = [ 'ClusterBuilderName' ] + if configDict['PrepareLCW']: + towerBuilder = CaloTopoClusterFromTowerMaker(configDict['ClusterBuilderName'],OrderClusterByPt=False) ### order by pt after LCW calibration! + excludedKeys += [ 'OrderClusterByPt' ] + else: + towerBuilder = CaloTopoClusterFromTowerMaker(configDict['ClusterBuilderName']) + + ''' Copy properties from dictionary ''' - towerConverter = CaloTopoClusterFromTowerMaker(cnvname,CaloTowerContainerKey=towerkey,CaloCellContainerKey=cellkey,OrderClusterByPt=ptorder) - ''' Refinement of converter configuration + for key,value in configDict.items(): + if key not in excludedKeys: + setattr(towerBuilder,key,value) + + ''' Check basic consistency of configuration + ''' + mlog.info('Consistency check') + if towerBuilder.PrepareLCW and not towerBuilder.BuildTopoTowers: + raise RuntimeError('{0}[inconsistent configuration] applying LCW requires to build topo-towers'.format(towerBuilder.name())) + if towerBuilder.BuildCombinedTopoSignal and not towerBuilder.BuildTopoTowers: + raise RuntimeError('{0}[inconsistent configuration] building combined topo-cluster/topo-tower signals requires to build topo-towers'.format(towerBuilder.name())) + if towerBuilder.ApplyCellEnergyThreshold and towerBuilder.BuildTopoTowers: + raise RuntimeError('{0}[inconsistent configuration] applying cell energy thresholds for topo-towres not yet implemented'.format(towerBuilder.name())) + + ''' Tower converter configuration summary + ''' + if towerBuilder.BuildTopoTowers: + if towerBuilder.PrepareLCW: + ''' LCW topo-towers + ''' + mlog.info('################################################') + mlog.info('## Produce LCW calibrated topo-tower clusters ##') + mlog.info('################################################') + mlog.info('CaloTopoClusterContainerKey .. {0}'.format(towerBuilder.CaloTopoClusterContainerKey)) + mlog.info('CellClusterWeightKey ......... {0}'.format(towerBuilder.CellClusterWeightKey)) + else: + ''' EM topo-towers + ''' + mlog.info('###############################################') + mlog.info('## Produce EM calibrated topo-tower clusters ##') + mlog.info('###############################################') + mlog.info('CaloTopoClusterContainerKey .. {0}'.format(towerBuilder.CaloTopoClusterContainerKey)) + + if towerBuilder.BuildCombinedTopoSignal: + mlog.info(' ') + mlog.info('Combined topo-cluster/topo-towermode with y_max = {0}'.format(towerBuilder.TopoClusterRange)) + else: + ''' EM towers + ''' + mlog.info('########################################') + mlog.info('## Produce EM standard tower clusters ##') + mlog.info('########################################') + + ''' Set debug flag (not a property in the dictionary) ''' - mlog.info(' ') - if buildtt: - mlog.info('################################################') - mlog.info('## Produce LCW calibrated topo-tower clusters ##') - mlog.info('################################################') - mlog.info('CaloTopoClusterContainerKey .. {0}'.format(tcluskey)) - mlog.info('CellClusterWeightKey ......... {0}'.format(tcwgtkey)) - towerConverter.CaloTopoClusterContainerKey = tcluskey - towerConverter.CellClusterWeightKey = tcwgtkey - towerConverter.ApplyLCW = True - else: - mlog.info('####################################################') - mlog.info('## Produce EM calibrated inclusive tower clusters ##') - mlog.info('####################################################') - mlog.info(' ') - if applyEnergyThreshold: - towerConverter.CellEnergyThreshold = twralgo.CellEnergyThreshold if debugOn: - towerConverter.OutputLevel = Lvl.DEBUG - # setting up the moments: external tools + towerBuilder.OutputLevel = Lvl.DEBUG + + towerCoreName = towerMakerName + if towerCoreName.find('Builder') > 0: + (towerCoreName.replace('Builder',' ')).rstrip(' ') + elif towerCoreName.find('Maker') > 0: + (towerCoreName.replace('Maker',' ')).rstrip(' ') + + ############################ + ## Setting up the moments ## + ############################ + + ''' External tools for moment calculation + ''' from CaloTools.CaloNoiseToolDefault import CaloNoiseToolDefault - caloNoiseTool = CaloNoiseToolDefault() from AthenaCommon.AppMgr import ToolSvc - ToolSvc += caloNoiseTool + caloNoiseTool = CaloNoiseToolDefault() + ToolSvc += caloNoiseTool - # moment maker + ''' Cluster moment maker (all general moments) + ''' from CaloRec.CaloTopoClusterFlags import jobproperties - from AthenaCommon.SystemOfUnits import deg, GeV, MeV - from CaloRec.CaloRecConf import CaloClusterMomentsMaker - clusterMoments = CaloClusterMomentsMaker (clusterMakerName+'MomentMaker') - clusterMoments.MaxAxisAngle = 20*deg - clusterMoments.CaloNoiseTool = caloNoiseTool - clusterMoments.UsePileUpNoise = True + from AthenaCommon.SystemOfUnits import deg, GeV, MeV + from CaloRec.CaloRecConf import CaloClusterMomentsMaker + clusterMoments = CaloClusterMomentsMaker (towerMakerName+'MomentMaker') + clusterMoments.MaxAxisAngle = 20*deg + clusterMoments.CaloNoiseTool = caloNoiseTool + clusterMoments.UsePileUpNoise = True clusterMoments.TwoGaussianNoise = jobproperties.CaloTopoClusterFlags.doTwoGaussianNoise() clusterMoments.MinBadLArQuality = 4000 - clusterMoments.MomentsNames = ["FIRST_PHI" - ,"FIRST_ETA" - ,"SECOND_R" - ,"SECOND_LAMBDA" - ,"DELTA_PHI" - ,"DELTA_THETA" - ,"DELTA_ALPHA" - ,"CENTER_X" - ,"CENTER_Y" - ,"CENTER_Z" - ,"CENTER_MAG" - ,"CENTER_LAMBDA" - ,"LATERAL" - ,"LONGITUDINAL" - ,"FIRST_ENG_DENS" - ,"ENG_FRAC_EM" - ,"ENG_FRAC_MAX" - ,"ENG_FRAC_CORE" - ,"SECOND_ENG_DENS" - ,"ISOLATION" - ,"ENG_BAD_CELLS" - ,"N_BAD_CELLS" - ,"N_BAD_CELLS_CORR" - ,"BAD_CELLS_CORR_E" - ,"BADLARQ_FRAC" - ,"ENG_POS" - ,"SIGNIFICANCE" - ,"CELL_SIGNIFICANCE" - ,"CELL_SIG_SAMPLING" - ,"AVG_LAR_Q" - ,"AVG_TILE_Q" - ,"PTD" - ,"MASS" - ] - - # only add HV related moments if it is offline. + clusterMoments.MomentsNames = [ + "FIRST_PHI" + ,"FIRST_ETA" + ,"SECOND_R" + ,"SECOND_LAMBDA" + ,"DELTA_PHI" + ,"DELTA_THETA" + ,"DELTA_ALPHA" + ,"CENTER_X" + ,"CENTER_Y" + ,"CENTER_Z" + ,"CENTER_MAG" + ,"CENTER_LAMBDA" + ,"LATERAL" + ,"LONGITUDINAL" + ,"FIRST_ENG_DENS" + ,"ENG_FRAC_EM" + ,"ENG_FRAC_MAX" + ,"ENG_FRAC_CORE" + ,"SECOND_ENG_DENS" + ,"ISOLATION" + ,"ENG_BAD_CELLS" + ,"N_BAD_CELLS" + ,"N_BAD_CELLS_CORR" + ,"BAD_CELLS_CORR_E" + ,"BADLARQ_FRAC" + ,"ENG_POS" + ,"SIGNIFICANCE" + ,"CELL_SIGNIFICANCE" + ,"CELL_SIG_SAMPLING" + ,"AVG_LAR_Q" + ,"AVG_TILE_Q" + ,"PTD" + ,"MASS" + ] + + ''' HV related moments for offline data + ''' from IOVDbSvc.CondDB import conddb if not conddb.isOnline: from LArRecUtils.LArHVScaleRetrieverDefault import LArHVScaleRetrieverDefault - clusterMoments.LArHVScaleRetriever=LArHVScaleRetrieverDefault() - clusterMoments.MomentsNames += ["ENG_BAD_HV_CELLS" - ,"N_BAD_HV_CELLS" - ] + clusterMoments.LArHVScaleRetriever = LArHVScaleRetrieverDefault() + clusterMoments.MomentsNames += ["ENG_BAD_HV_CELLS","N_BAD_HV_CELLS"] - # cluster maker + ############################################################### + ## Set up the tower builder algorithm - as a cluster builder ## + ############################################################### + + ''' Basic algorithm properties + ''' from CaloRec.CaloRecConf import CaloClusterMaker - clusterMaker = CaloClusterMaker(clusterMakerName) - clusterMaker.ClustersOutputName = clusterContainerKey - clusterMaker.ClusterMakerTools = [ towerConverter ] - mlog.info('instantiated CaloClusterMaker "{0}"'.format(clusterMaker.name())) + towerMaker = CaloClusterMaker(towerMakerName) + towerMaker.ClustersOutputName = towerContainerKey + towerMaker.ClusterMakerTools = [ towerBuilder ] + towerMaker += towerBuilder + mlog.info('instantiated CaloClusterMaker "{0}" configuration'.format(towerMaker.name())) - # bad cell corrections -## from CaloClusterCorrection.CaloClusterBadChannelListCorr import CaloClusterBadChannelListCorr -## badChannelCorr = CaloClusterBadChannelListCorr() + ''' Set up bad cell corrections + ''' + from CaloClusterCorrection.CaloClusterBadChannelListCorr import CaloClusterBadChannelListCorr + badChannelCorr = CaloClusterBadChannelListCorr() - # Correction tools -## clusterMaker.ClusterCorrectionTools += [ badChannelCorr ] - clusterMaker.ClusterCorrectionTools += [ clusterMoments ] + ''' Register correction and moment tools + ''' + towerMaker.ClusterCorrectionTools += [ badChannelCorr ] + towerMaker.ClusterCorrectionTools += [ clusterMoments ] + towerMaker += clusterMoments - # configuring the algorithm - clusterMaker += towerConverter -## clusterMaker += badChannelCorr - clusterMaker += clusterMoments + #################################### + ## Configure LCW calibration tool ## + #################################### - if buildtt: + if towerBuilder.PrepareLCW: from CaloRec.CaloRecConf import CaloTopoClusterFromTowerCalibrator - calgname = clusterMakerName+'Calibrator' - mlog.info('TopoTowers: add LCW calibration tool <'+calgname+'>') - clusterCalibrator = CaloTopoClusterFromTowerCalibrator(calgname) - mlog.info('TopoTowers: '+calgname+'.CellClusterWeightKey = "'+tcwgtkey+'"') - clusterCalibrator.CellClusterWeightKey = tcwgtkey - clusterCalibrator.OrderClusterByPt = ptorder - clusterMaker.ClusterCorrectionTools += [ clusterCalibrator ] - clusterMaker += clusterCalibrator - - # done - return clusterMaker + ''' Configure name for calibration tool + ''' + towerCalName = towerCoreName+'Calibrator' + towerCalibrator = CaloTopoClusterFromTowerCalibrator(towerCalName) + mlog.info('add LCW calibration tool <'+towerCalName+'>') + mlog.info('TopoTowers: '+towerCalName+'.CellClusterWeightKey = "'+towerBuilder.CellClusterWeightKey+'"') + towerCalibrator.CellClusterWeightKey = towerBuilder.CellClusterWeightKey + towerCalibrator.OrderClusterByPt = configDict['OrderClusterByPt'] + if debugOn: + towerCalibrator.OutputLevel = Lvl.DEBUG + ''' Schedule calibration tool + ''' + towerMaker.ClusterCorrectionTools += [ towerCalibrator ] + towerMaker += towerCalibrator + + ####################### + # Configuration done ## + ####################### + + return towerMaker + +## +## toolname = configDict['ClusterBuilderName'] ### name of the tower builder tool +## cellkey = configDict['CaloCellContainerKey'] ### cell container key +## buildtopotower = configDict['BuildTopoTowers'] ### controls if topo-towers or inclusive towers are built +## towergeosvc = configDict['CaloTowerGeometrySvc'] ### tower geometry provider +## if ( buildtopotower ): +## topoclusterkey = configDict['CaloTopoClusterContainerKey'] +## else: +## topoclusterkey = 'N/A' +## +## cellweightkey = configDict['CellClusterWeightKey'] +## +## mlog.info('(input) CaloCellContainer <'+cellkey+'>') +## mlog.info('(input) CaloTopoClusterContainer <'+topoclusterkey+'>') +## mlog.info('(input) CellClusterWeightKey <'+cellweightkey+'>') diff --git a/Calorimeter/CaloRec/share/CaloRec_jobOptions.py b/Calorimeter/CaloRec/share/CaloRec_jobOptions.py index de9da525e1346511f3fee547907118882a33ca95..100a5278e0f7c6607d8037b36ee26899478bc003 100644 --- a/Calorimeter/CaloRec/share/CaloRec_jobOptions.py +++ b/Calorimeter/CaloRec/share/CaloRec_jobOptions.py @@ -248,14 +248,14 @@ else: # # functionality : Noise suppressed tower # -if jobproperties.CaloRecFlags.doCaloTopoTower() and DetFlags.haveRIO.Calo_on(): - try: - include ("CaloRec/CaloTopoTower_jobOptions.py") - except Exception: - treatException("Problem with CaloTopoTower. Switched off.") - jobproperties.CaloRecFlags.doCaloTopoTower=False -else: - jobproperties.CaloRecFlags.doCaloTopoTower=False +#if jobproperties.CaloRecFlags.doCaloTopoTower() and DetFlags.haveRIO.Calo_on(): +# try: +# include ("CaloRec/CaloTopoTower_jobOptions.py") +# except Exception: +# treatException("Problem with CaloTopoTower. Switched off.") +# jobproperties.CaloRecFlags.doCaloTopoTower=False +#else: +# jobproperties.CaloRecFlags.doCaloTopoTower=False # # functionality : muon candidates in Tile @@ -361,3 +361,10 @@ if rec.doWritexAOD(): #L1Calo Trigger tower decoration if globalflags.DataSource()=='data' and rec.doESD() and rec.doCalo() and rec.doTrigger(): include("TrigT1CaloCalibTools/DecorateL1CaloTriggerTowers_prodJobOFragment.py") + +#new style CaloTopoTowers +if jobproperties.CaloRecFlags.doCaloTopoTower(): + include ( "CaloRec/CaloTopoTowerFragment.py" ) +#mixed topo-cluster/topo-tower +if jobproperties.CaloRecFlags.doCaloTopoSignal(): + include ("CaloRec/CaloTopoSignalFragment.py" ) diff --git a/Calorimeter/CaloRec/share/CaloTopoSignalFragment.py b/Calorimeter/CaloRec/share/CaloTopoSignalFragment.py new file mode 100644 index 0000000000000000000000000000000000000000..5e936582044478363aab955eff6ba93da427ad18 --- /dev/null +++ b/Calorimeter/CaloRec/share/CaloTopoSignalFragment.py @@ -0,0 +1,65 @@ +###################################### +## Create standard 0.1 x 0.1 towers ## +###################################### + + +from AthenaCommon.Logging import logging +mlog = logging.getLogger('CaloTopoSognalFragment.py:: ') + +import AthenaCommon.Constants as Lvl +from AthenaCommon.AppMgr import ServiceMgr as svcMgr + +from CaloRec.MakeClustersFromTowers import ClustersFromTowersDict, MakeClustersFromTowers +from CaloRec.CaloRecConf import CaloTowerGeometrySvc, CaloTopoClusterTowerMerger + +mlog.info(' ') +mlog.info('##################################') +mlog.info('## Topological Signal Formation ##') +mlog.info('##################################') +mlog.info(' ') + +############################# +## Tower Geometry Provider ## +############################# + +if not hasattr(svcMgr,'CaloTowerGeometryProvider'): + mlog.info("setting up tower geometry provider") + caloTowerGeoSvc = CaloTowerGeometrySvc('CaloTowerGeometryProvider') + caloTowerGeoSvc.TowerEtaBins = 100 + caloTowerGeoSvc.TowerEtaMin = -5. + caloTowerGeoSvc.TowerEtaMax = 5. + svcMgr += caloTowerGeoSvc + + +############################# +## CaloTopoTower Formation ## +############################# + +caloTowerDict = ClustersFromTowersDict(clusterBuilderName='CaloFwdTopoTowerBuilder', + towerGeometrySvc=svcMgr.CaloTowerGeometryProvider, + cellContainerKey='AllCalo', + buildTopoTowers=True, + topoClusterContainerKey='CaloCalTopoClusters', + cellClusterWeightKey='CaloCalFwdTopoTowerCellWeights', + orderClusterByPt=False, + applyCellEnergyThreshold=False, + doCellIndexCheck=False, + cellEnergyThreshold=0., + applyLCW=True, + buildCombinedSignal=True, + clusterRange=2.5) + +caloTowerAlgo = MakeClustersFromTowers(towerMakerName = 'CaloFwdTopoTowerMaker', + towerContainerKey = 'CaloCalFwdTopoTowers', + configDict = caloTowerDict, + debugOn = False) +#merging +caloTowerMerger = CaloTopoClusterTowerMerger("CaloTopoSignalMaker") +caloTowerMerger.TopoClusterRange = caloTowerAlgo.CaloFwdTopoTowerBuilder.TopoClusterRange +caloTowerMerger.TopoClusterContainerKey = caloTowerAlgo.CaloFwdTopoTowerBuilder.CaloTopoClusterContainerKey +caloTowerMerger.TopoTowerContainerKey = caloTowerAlgo.ClustersOutputName +caloTowerMerger.TopoSignalContainerKey = 'CaloCalTopoSignals' +caloTowerMerger.OutputLevel = Lvl.DEBUG + +topSequence+=caloTowerAlgo +topSequence+=caloTowerMerger diff --git a/Calorimeter/CaloRec/share/CaloTopoTowerFragment.py b/Calorimeter/CaloRec/share/CaloTopoTowerFragment.py new file mode 100644 index 0000000000000000000000000000000000000000..56abf8c621187e77f2187bc65d939d3c759e8082 --- /dev/null +++ b/Calorimeter/CaloRec/share/CaloTopoTowerFragment.py @@ -0,0 +1,56 @@ +###################################### +## Create standard 0.1 x 0.1 towers ## +###################################### + + +from AthenaCommon.Logging import logging +mlog = logging.getLogger('CaloTopoTowerFragment.py:: ') + +import AthenaCommon.Constants as Lvl +from AthenaCommon.AppMgr import ServiceMgr as svcMgr + +from CaloRec.MakeClustersFromTowers import ClustersFromTowersDict, MakeClustersFromTowers +from CaloRec.CaloRecConf import CaloTowerGeometrySvc + +mlog.info(' ') +mlog.info('##################################') +mlog.info('## Standard Tower Configuration ##') +mlog.info('##################################') +mlog.info(' ') + +############################# +## Tower Geometry Provider ## +############################# + +if not hasattr(svcMgr,'CaloTowerGeometryProvider'): + mlog.info("setting up tower geometry provider") + caloTowerGeoSvc = CaloTowerGeometrySvc('CaloTowerGeometryProvider') + caloTowerGeoSvc.TowerEtaBins = 100 + caloTowerGeoSvc.TowerEtaMin = -5. + caloTowerGeoSvc.TowerEtaMax = 5. + svcMgr += caloTowerGeoSvc + +############################# +## CaloTopoTower Formation ## +############################# + +caloTowerDict = ClustersFromTowersDict(clusterBuilderName='CaloTopoTowerBuilder', + towerGeometrySvc=svcMgr.CaloTowerGeometryProvider, + cellContainerKey='AllCalo', + buildTopoTowers=True, + topoClusterContainerKey='CaloCalTopoClusters', + cellClusterWeightKey='CaloCalTopoTowerCellWeights', + orderClusterByPt=False, + applyCellEnergyThreshold=False, + doCellIndexCheck=False, + cellEnergyThreshold=0., + applyLCW=True, + buildCombinedSignal=False, + clusterRange=5.) + +caloTowerAlgo = MakeClustersFromTowers(towerMakerName = 'CaloTopoTowerMaker', + towerContainerKey = 'CaloCalTopoTowers', + configDict = caloTowerDict, + debugOn = False) + +topSequence+=caloTowerAlgo diff --git a/Calorimeter/CaloRec/src/CaloClusterMomentsMaker.cxx b/Calorimeter/CaloRec/src/CaloClusterMomentsMaker.cxx index 7b80840b591317c65e1903bf01bdcd92a6cfcedb..bdf1d6d8c7fea57daa37ba28b38bf744f80b4ad9 100644 --- a/Calorimeter/CaloRec/src/CaloClusterMomentsMaker.cxx +++ b/Calorimeter/CaloRec/src/CaloClusterMomentsMaker.cxx @@ -40,6 +40,7 @@ #include #include #include +//#include using HepGeom::Vector3D; @@ -375,20 +376,22 @@ StatusCode CaloClusterMomentsMaker::execute(xAOD::CaloClusterContainer *theClusC // loop over all cell members and fill cell vector for used cells xAOD::CaloCluster::cell_iterator cellIter = theCluster->cell_begin(); xAOD::CaloCluster::cell_iterator cellIterEnd = theCluster->cell_end(); + // printf("CaloClusterMomentMaker::execute(...) - number of cells cluster %4i is %zu\n",iClus,theCluster->size()); + // int iCell(0); for(; cellIter != cellIterEnd; cellIter++ ){ - CxxUtils::prefetchNext(cellIter, cellIterEnd); + CxxUtils::prefetchNext(cellIter, cellIterEnd); const CaloCell* pCell = *cellIter; - - Identifier myId = pCell->ID(); - IdentifierHash myHashId = m_calo_id->calo_cell_hash(myId); - if ( clusterIdx[(unsigned int)myHashId].first != noCluster) { - // check weight and assign to current cluster if weight is > 0.5 - double weight = cellIter.weight(); - if ( weight > 0.5 ) + // printf("CaloClusterMomentMaker::execute(...) - cell %4i/%4zu reference %p\n",++iCell,theCluster->size(),(void*)pCell); + if ( pCell != 0 ) { + Identifier myId = pCell->ID(); + IdentifierHash myHashId = m_calo_id->calo_cell_hash(myId); + if ( clusterIdx[(unsigned int)myHashId].first != noCluster) { + // check weight and assign to current cluster if weight is > 0.5 + double weight = cellIter.weight(); + if ( weight > 0.5 ) clusterIdx[(unsigned int)myHashId].first = iClus; + } else { clusterIdx[(unsigned int)myHashId].first = iClus; - } - else { - clusterIdx[(unsigned int)myHashId].first = iClus; + } } } ++iClus; diff --git a/Calorimeter/CaloRec/src/CaloTopoClusterFromTowerCalibrator.cxx b/Calorimeter/CaloRec/src/CaloTopoClusterFromTowerCalibrator.cxx index 4a7c082f7420d3140da8abb728c82cd6522f3213..c75e4053aabb4e92aa9b3dacef043305f183d58c 100644 --- a/Calorimeter/CaloRec/src/CaloTopoClusterFromTowerCalibrator.cxx +++ b/Calorimeter/CaloRec/src/CaloTopoClusterFromTowerCalibrator.cxx @@ -1,9 +1,12 @@ /* Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ + #include "CaloRec/CaloTopoClusterFromTowerCalibrator.h" #include "CaloRec/CaloTopoClusterFromTowerHelpers.h" #include "CaloEvent/CaloCellClusterWeights.h" +#include "CaloGeoHelpers/CaloSampling.h" + #include namespace { ///////////////////////////////////////////////////////////// @@ -16,24 +19,30 @@ namespace { } } -std::string CaloTopoClusterFromTowerCalibrator::m_defaultKey = "NONE"; - CaloTopoClusterFromTowerCalibrator::CaloTopoClusterFromTowerCalibrator(const std::string& type,const std::string& name,const IInterface* pParent) : AthAlgTool(type,name,pParent) - , m_cellClusterWeightsKey(m_defaultKey) + , m_cellClusterWeightKey("CaloTopoClusterCellWeights") , m_orderByPt(false) { declareInterface(this); - declareProperty("CellClusterWeightKey",m_cellClusterWeightsKey,"SG Key for CellClusterWeights (input)"); - declareProperty("OrderClusterByPt", m_orderByPt, "Order clusters by calibrated Pt (input)"); + declareProperty("CellClusterWeightKey",m_cellClusterWeightKey,"SG Key for CellClusterWeights (input)"); + declareProperty("OrderClusterByPt", m_orderByPt, "Order clusters by calibrated Pt (input)"); } +StatusCode CaloTopoClusterFromTowerCalibrator::initialize() +{ + ATH_CHECK(m_cellClusterWeightKey.initialize()); + return StatusCode::SUCCESS; +} StatusCode CaloTopoClusterFromTowerCalibrator::execute(xAOD::CaloClusterContainer* pClusCont) { // retrieve weights - const CaloCellClusterWeights* pCellWeights = 0; - CHECK(evtStore()->retrieve(pCellWeights,CaloCellClusterWeights::key(m_cellClusterWeightsKey) )); // PA change + SG::ReadHandle pCellWeights(m_cellClusterWeightKey); + if ( !pCellWeights.isValid() ) { + ATH_MSG_ERROR( "Cannot allocate CaloCellClusterWeights with key <" << m_cellClusterWeightKey.key() << ">" ); + return StatusCode::FAILURE; + } ///////////////////////// // Calibrated clusters // @@ -46,7 +55,15 @@ StatusCode CaloTopoClusterFromTowerCalibrator::execute(xAOD::CaloClusterContaine const CaloCellClusterWeights::weight_t& wght(pCellWeights->at(*fCell)); // Retrieve list of weights associated with this cluster. // double weight(fCell.weight()); // Get cell-in-tower weight. // weight *= accumulateWeight(wght); // Combine with calibration weights. // - pClus->reweightCell(fCell,weight); // Set new weight. // + if ( weight == 0. ) { + ATH_MSG_DEBUG( CaloRec::Helpers::fmtMsg("[NO_LCW_REWEIGHT] Tower (%6.3f,%6.3f) cell [%6zu] weight = %6.3f [# LCW weights %zu geo %6.3f LCW %6.3f] SamplingID %2u Name \042%s\042", + pClus->eta(),pClus->phi(),(size_t)fCell->caloDDE()->calo_hash(),weight,wght.size(),fCell.weight(),weight/std::max(fCell.weight(),1e-08), + (unsigned int)fCell->caloDDE()->getSampling(),CaloSampling::getSamplingName(fCell->caloDDE()->getSampling()).c_str()) ); + } else { + ATH_MSG_DEBUG( CaloRec::Helpers::fmtMsg("[DO_LCW_REWEIGHT] Tower (%6.3f,%6.3f) cell [%6zu] weight = %6.3f [# LCW weights %zu geo %6.3f LCW %6.3f]", + pClus->eta(),pClus->phi(),(size_t)fCell->caloDDE()->calo_hash(),weight,wght.size(),fCell.weight(),weight/fCell.weight()) ); + pClus->reweightCell(fCell,weight); // Set new weight. // + } } //////////////////////////////////////////////////////////// // preserve raw (EM) kinematics double rawE(pClus->e()); //////////////////////////////////////////////////////////////////////////////////////////// diff --git a/Calorimeter/CaloRec/src/CaloTopoClusterFromTowerMaker.cxx b/Calorimeter/CaloRec/src/CaloTopoClusterFromTowerMaker.cxx index 2ef79a197b00eda28935230512d108ad6ea85212..bb4fa41e25257d1de7b3db657191bcbbb1f17616 100644 --- a/Calorimeter/CaloRec/src/CaloTopoClusterFromTowerMaker.cxx +++ b/Calorimeter/CaloRec/src/CaloTopoClusterFromTowerMaker.cxx @@ -1,14 +1,15 @@ -/* Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration */ +/* Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration */ +#include "GaudiKernel/MsgStream.h" +#include "GaudiKernel/SystemOfUnits.h" + +#include "AthenaKernel/Units.h" + #include "CaloRec/CaloTopoClusterFromTowerMaker.h" #include "CaloRec/CaloTopoClusterFromTowerHelpers.h" - -#include "xAODCaloEvent/CaloTower.h" -#include "xAODCaloEvent/CaloTower.h" #include "xAODCaloEvent/CaloClusterKineHelper.h" #include "CaloEvent/CaloClusterCellLink.h" -#include "CaloEvent/CaloCellTowerLink.h" #include "CaloEvent/CaloCellClusterWeights.h" #include "CaloGeoHelpers/CaloSampling.h" @@ -16,277 +17,280 @@ #include "CaloGeoHelpers/proxim.h" #include "CaloRec/CaloProtoCluster.h" +#include "CaloRec/CaloTowerGeometrySvc.h" + +#include "CaloDetDescr/CaloDetDescrElement.h" +#include "CaloDetDescr/CaloDetDescrManager.h" #include #include +#include #include #include -#include #include #include +#include +#include -////////////////////////////////////////////////// -// CaloTopoClusterFromTowerMaker::CellPathology // -////////////////////////////////////////////////// - -CaloTopoClusterFromTowerMaker::CellPathology::CellPathology() -{ } - -void CaloTopoClusterFromTowerMaker::CellPathology::setInvalidPointer(int cellhash,const CaloCell* cptr) -{ - std::string msg(CaloRec::Helpers::fmtMsg("CaloCell [%06i] has invalid pointer %p",cellhash,(void*)cptr)); - _addMsg(msg,_invalidPointer); +namespace { + MsgStream& operator<<(MsgStream& mstr,const SG::ReadHandleKey& ckey) { mstr << ckey.key(); return mstr; } + MsgStream& operator<<(MsgStream& mstr,const SG::ReadHandleKey& ckey) { mstr << ckey.key(); return mstr; } + MsgStream& operator<<(MsgStream& mstr,const SG::WriteHandleKey& ckey) { mstr << ckey.key(); return mstr; } } -void CaloTopoClusterFromTowerMaker::CellPathology::setInvalidHash(int cellhash,const CaloCell* cptr) -{ - std::string msg(CaloRec::Helpers::fmtMsg("CaloCell [%06i] at %p has invalid hash %i",cellhash,(void*)cptr,cellhash)); - _addMsg(msg,_invalidHash); -} - -void CaloTopoClusterFromTowerMaker::CellPathology::setInvalidSampling(int cellhash,const CaloCell* cptr) -{ - int hash(CaloSampling::Unknown); - if ( cptr != 0 ) { hash = (int)cptr->caloDDE()->getSampling(); } - std::string msg(CaloRec::Helpers::fmtMsg("CaloCell [%06i] at %p has invalid sampling %2i (not in [0,%2i]",cellhash,(void*)cptr,hash,(int)CaloSampling::Unknown)); - _addMsg(msg,_invalidSampling); -} +std::atomic CaloTopoClusterFromTowerMaker_checkCellIndices(false); /////////////////////////////////// // CaloTopoClusterFromTowerMaker // /////////////////////////////////// -double CaloTopoClusterFromTowerMaker::m_energyThresholdDef = -100000000.; // in MeV -std::string CaloTopoClusterFromTowerMaker::m_defaultKey = "NONE"; +double CaloTopoClusterFromTowerMaker::m_energyThresholdDef = -100000000.; // in MeV +double CaloTopoClusterFromTowerMaker::m_clusterRangeDef = 5.; +std::string CaloTopoClusterFromTowerMaker::m_defaultKey = "NONE"; +CaloTopoClusterFromTowerMaker::uint_t CaloTopoClusterFromTowerMaker::m_errorValueUINT = uint_t(-1); CaloTopoClusterFromTowerMaker::CaloTopoClusterFromTowerMaker(const std::string& type, const std::string& name, const IInterface* pParent) : AthAlgTool(type,name,pParent) - , m_towerContainerKey(m_defaultKey) - , m_cellTowerLinksKey(m_defaultKey) + , m_towerGeometrySvc("CaloTowerGeometrySvc",name) + , m_clusterContainerKey("CaloTopoCluster") , m_cellContainerKey("AllCalo") - , m_clusterContainerKey(m_defaultKey) - , m_cellClusterWeightKey(m_defaultKey) + , m_cellClusterWeightKey("CaloTopoClusterCellWeights") , m_orderByPt(false) - , m_energyThreshold(m_energyThresholdDef-1.) - , m_applyLCW(false) - , m_applyEnergyThreshold(false) - , m_useCellsFromClusters(false) , m_prepareLCW(false) + , m_useCellsFromClusters(true) + , m_applyCellEnergyThreshold(false) + , m_doCellIndexCheck(false) + , m_buildCombinedSignal(false) + , m_energyThreshold(m_energyThresholdDef-1.) + , m_clusterRange(m_clusterRangeDef) , m_numberOfCells(0) - , m_numberOfSamplings(CaloSampling::Unknown) - , m_cellTowerLinksHandle(0) - , m_cellClusterWeights(0) + , m_maxCellHash(0) + , m_numberOfSamplings(static_cast(CaloSampling::Unknown)) + , m_numberOfTowers(0) { declareInterface(this); - declareProperty("CaloTowerContainerKey", m_towerContainerKey, "SG Key for CaloTowerContainer (input)"); - declareProperty("CaloCellContainerKey", m_cellContainerKey, "SG Key for CaloCellContainer (input)"); - declareProperty("CaloTopoClusterContainerKey", m_clusterContainerKey, "SG Key for CaloClusterContainer (input)"); - declareProperty("CellClusterWeightKey", m_cellClusterWeightKey,"SG Key for CellClusterWeights (output)"); - declareProperty("OrderClusterByPt", m_orderByPt, "Turn on/off pT-ordering of CaloClusterContainer (output)"); - declareProperty("CellEnergyThreshold", m_energyThreshold, "Energy threshold for cells filled in clusters"); - declareProperty("ApplyLCW", m_applyLCW, "Prepare data structure to apply LCW"); + declareProperty("CaloTowerGeometrySvc", m_towerGeometrySvc=ServiceHandle("CaloTowerGeometrySvc",name), "Service providing tower geometry"); + declareProperty("CaloCellContainerKey", m_cellContainerKey, "SG Key for CaloCellContainer (input)"); + declareProperty("BuildTopoTowers", m_useCellsFromClusters, "Turn on/off topo-tower formation"); + declareProperty("CaloTopoClusterContainerKey", m_clusterContainerKey, "SG Key for CaloClusterContainer (input)"); + declareProperty("CellClusterWeightKey", m_cellClusterWeightKey, "SG Key for CellClusterWeights (output)"); + declareProperty("OrderClusterByPt", m_orderByPt, "Turn on/off pT-ordering of CaloClusterContainer (output)"); + declareProperty("ApplyCellEnergyThreshold", m_applyCellEnergyThreshold, "Turn on/off cell energy thresholds"); + declareProperty("CellEnergyThreshold", m_energyThreshold, "Energy threshold for cells filled in clusters"); + declareProperty("PrepareLCW", m_prepareLCW, "Prepare data structure to apply LCW"); + declareProperty("ExcludedSamplings", m_excludedSamplingsName, "Excluded samplings by name"); + declareProperty("DoCellIndexCheck", m_doCellIndexCheck, "Check cell hash indices for consistency"); + declareProperty("BuildCombinedTopoSignal", m_buildCombinedSignal, "Build topo-clusters and topo-towers"); + declareProperty("TopoClusterRange", m_clusterRange, "Rapidity range for using topo-clusters in combined signal mode"); } StatusCode CaloTopoClusterFromTowerMaker::initialize() { - // build SG key for cell-tower link map object - m_cellTowerLinksKey = m_towerContainerKey + CaloCellTowerLink::extTag; - ATH_MSG_INFO("Input retrieved from CaloTowerxAODContainer <" << m_towerContainerKey << "> and CaloCellTowerLink::Map <" << m_cellTowerLinksKey << ">"); - // check if energy threshold is to be applied - m_applyEnergyThreshold = m_energyThreshold > m_energyThresholdDef; - // check if filtered cells to be used - m_useCellsFromClusters = m_clusterContainerKey != m_defaultKey; - // analyze configuration + //--------------------// + // Set up handle keys // + //--------------------// + + ATH_CHECK(m_cellContainerKey.initialize()); + + //---------------------// + // Check configuration // + //---------------------// + + // tower geometry service + ATH_MSG_INFO("Allocate tower geometry service:"); + if ( !m_towerGeometrySvc.isValid() ) { + ATH_MSG_ERROR("[reject] cannot allocate tower geometry service - fatal"); + return StatusCode::FAILURE; + } else { + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("Tower geometry service is allocated, describes %6zu towers in grid:", m_towerGeometrySvc->towerBins()) ); + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("[accept] %3zu eta bins in [%5.2f,%5.2f]",m_towerGeometrySvc->etaBins(),m_towerGeometrySvc->etaMin(),m_towerGeometrySvc->etaMax()) ); + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("[accept] %3zu phi bins in [%5.2f,%5.2f]",m_towerGeometrySvc->phiBins(),m_towerGeometrySvc->phiMin(),m_towerGeometrySvc->phiMax()) ); + } + + // tower builder configurations if ( m_useCellsFromClusters ) { - if ( m_applyEnergyThreshold ) { - ATH_MSG_WARNING("Misconfiguration - cannot apply energy threshold to cells from TopoClusters, will accept all clustered cells"); - m_applyEnergyThreshold = false; - } else { - ATH_MSG_INFO("Filtered mode, cluster from towers filled with cells from TopoClusters"); + // topo-tower + ATH_MSG_INFO("Configure for building topo-towers (filtered mode):"); + // energy threshold not (yet) implemented for topo-towers + if ( m_applyCellEnergyThreshold ) { + ATH_MSG_WARNING("[ignore] cannot apply energy thresholds to topo-towers!"); + m_applyCellEnergyThreshold = false; } + ATH_CHECK(m_clusterContainerKey.initialize()); + // check on request for LCW + if ( m_prepareLCW ) { + ATH_CHECK(m_cellClusterWeightKey.initialize()); + ATH_MSG_INFO("[accept] prepare for LCW calibration - initialize CaloCellClusterWeights key object <" << m_cellClusterWeightKey << ">"); + } else { + ATH_MSG_INFO("[accept] use EM scale"); + } + } else { + // inclusive/exclusive towers + ATH_MSG_INFO("Configure for building cell towers:"); + if ( m_applyCellEnergyThreshold ) { + ATH_MSG_INFO("[accept] configure exclusive towers: use cell energy threshold"); + if ( m_energyThreshold < m_energyThresholdDef ) { + ATH_MSG_ERROR("######## [reject] invalid cell energy threshold " << m_energyThreshold/Athena::Units::GeV + << " GeV is smaller than default (no-op) " << m_energyThresholdDef/Athena::Units::GeV << " GeV - fatal"); + return StatusCode::FAILURE; + } + ATH_MSG_INFO("######## [accept] energy threshold for cells to contribute to towers is " << m_energyThreshold/Athena::Units::GeV << " GeV"); + } else { + ATH_MSG_INFO("[accept] configure inclusive towers"); + } // end inclusive/exclusive tower configuration + } // end tower builder configuration + + // local data (constant parameters) + m_numberOfCells = m_towerGeometrySvc->totalNumberCells(); + m_maxCellHash = m_towerGeometrySvc->maxCellHash(); + m_numberOfTowers = m_towerGeometrySvc->towerBins(); + ATH_MSG_INFO("Additional tool parameters:"); + if ( m_numberOfCells > 0 ) { + ATH_MSG_INFO("[accept] maximum cell hash index is " << m_maxCellHash); + ATH_MSG_INFO("[accept] maximum number of cells is " << m_numberOfCells); } else { - if ( m_applyEnergyThreshold ) { - ATH_MSG_INFO(CaloRec::Helpers::fmtMsg("Exclusive mode, cluster from towers filled with cells with energy > %.3f MeV",m_energyThreshold)); - } else { - ATH_MSG_INFO("Inclusive mode, all cells are used"); + ATH_MSG_ERROR("[reject] invalid maximum cell hash index/total number of cells " << m_maxCellHash << "/" << m_numberOfCells << " - fatal"); + return StatusCode::FAILURE; + } + if ( m_numberOfTowers > 0 ) { + ATH_MSG_INFO("[accept] maximum number of towers is " << m_numberOfTowers); + } else { + ATH_MSG_ERROR("[reject] invalid maximum number of towers " << m_numberOfTowers << " - fatal"); + return StatusCode::FAILURE; + } + + if ( m_excludedSamplingsName.empty() ) { + m_excludedSamplings.clear(); + m_excludedSamplingsPattern.reset(); + ATH_MSG_INFO("Cells from all samplings used for topo-cluster included"); + } else { + size_t nex(std::min(m_excludedSamplingsName.size(), m_excludedSamplingsPattern.size())); + if ( m_excludedSamplingsName.size() > m_excludedSamplingsPattern.size() ) { + ATH_MSG_WARNING( CaloRec::Helpers::fmtMsg("Configuration problem: number of excluded sampling names %zu exceeds expected maximum %zu - ignore last %zu name(s)", + m_excludedSamplingsName.size(), m_excludedSamplingsPattern.size(),m_excludedSamplingsName.size()-m_excludedSamplingsPattern.size()) ); + } + m_excludedSamplings.resize(nex); + m_excludedSamplingsPattern.reset(); + for ( size_t i(0); i blu { { true, "true" }, { false, "false" } }; + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("PrepareLCW ................. %s", blu[m_prepareLCW].c_str()) ); + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("BuildTopoTowers ............ %s", blu[m_useCellsFromClusters].c_str()) ); + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("ApplyCellEnergyThreshold ... %s", blu[m_applyCellEnergyThreshold].c_str()) ); + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("OrderClusterByPt ........... %s", blu[m_orderByPt].c_str()) ); + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("DoCellIndexCheck ........... %s", blu[m_doCellIndexCheck].c_str()) ); + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("BuildCombinedTopoSignal .... %s", blu[m_buildCombinedSignal].c_str()) ); + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("TopoClusterRange ........... %.2f", m_clusterRange) ); + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("ExcludedSamplings .......... %zu (number of)",m_excludedSamplingsName.size()) ); return StatusCode::SUCCESS; } StatusCode CaloTopoClusterFromTowerMaker::finalize() +{ return StatusCode::SUCCESS; } + +StatusCode CaloTopoClusterFromTowerMaker::execute(xAOD::CaloClusterContainer* pClusCont) { - ATH_MSG_INFO(" "); - ATH_MSG_INFO("Summary of warnings in execute():"); - ATH_MSG_INFO("Count: Warning:"); - for ( auto fiter(m_warning.begin()); fiter!=m_warning.end(); ++fiter ) { - ATH_MSG_INFO(CaloRec::Helpers::fmtMsg("%6i \042%s\042",(int)fiter->second,fiter->first.c_str())); + ///////////////// + // Check input // + ///////////////// + + // CaloCellContainer is needed to construct CaloProtoCluster + SG::ReadHandle pCellCont(m_cellContainerKey); + if ( !pCellCont.isValid() ) { + ATH_MSG_ERROR("Cannot allocate CaloCellContainer with key <" << m_cellContainerKey << ">"); + return StatusCode::FAILURE; } - ATH_MSG_INFO(" "); - ATH_MSG_INFO("Summary of messages in execute():"); - ATH_MSG_INFO("Count: Message:"); - for ( auto fiter(m_message.begin()); fiter!=m_message.end(); ++fiter ) { - ATH_MSG_INFO(CaloRec::Helpers::fmtMsg("%6i \042%s\042",(int)fiter->second,fiter->first.c_str())); + if ( msgLvl(MSG::DEBUG) && m_numberOfCells != pCellCont->size() ) { + ATH_MSG_DEBUG( CaloRec::Helpers::fmtMsg("[mismatch] number of cells in CaloCellContainer %6zu, total number of cell descriptors %6zu", + pCellCont->size(),m_towerGeometrySvc->totalNumberCells()) ); } - if ( !m_cellProblems.invalidPointer().empty() ) { - ATH_MSG_INFO(" "); - ATH_MSG_INFO("Observed cell pathologies [invalid pointer]:"); - printProblems(m_cellProblems.invalidPointer()); - } - if ( !m_cellProblems.invalidHash().empty() ) { - ATH_MSG_INFO(" "); - ATH_MSG_INFO("Observed cell pathologies [invalid hash]:"); - printProblems(m_cellProblems.invalidHash()); - } - if ( !m_cellProblems.invalidSampling().empty() ) { - ATH_MSG_INFO(" "); - ATH_MSG_INFO("Observed cell pathologies [invalid sampling]:"); - printProblems(m_cellProblems.invalidSampling()); - } - ATH_MSG_INFO(" "); - - return StatusCode::SUCCESS; -} -StatusCode CaloTopoClusterFromTowerMaker::execute(xAOD::CaloClusterContainer* pClusCont) -{ - // find cell-tower links - CHECK(detStore()->retrieve(m_cellTowerLinksHandle,m_cellTowerLinksKey)); - - // find input: CaloCell container - const CaloCellContainer* pCellCont = 0; - CHECK(evtStore()->retrieve(pCellCont,m_cellContainerKey)); - m_numberOfCells = (int)pCellCont->size(); - - // cross check - this->addMessage(CaloRec::Helpers::fmtMsg("CaloCellContainer::%s size %i, CaloCellTowerlink::Map size %i", - m_cellContainerKey.c_str(),(int)pCellCont->size(),(int)m_cellTowerLinksHandle->size())); - for ( size_t itl(0); itlsize(); ++itl ) { - if ( itl >= pCellCont->size() ) { - this->addWarning(CaloRec::Helpers::fmtMsg("CaloCell [%06i] has hash out of reach (max hash is %i)",(int)itl,(int)pCellCont->size()-1)); - } else { - if ( pCellCont->at(itl) == 0 ) { this->addWarning(CaloRec::Helpers::fmtMsg("CaloCell [%06i] has invalid pointer",(int)itl)); } - } - } - // int icell(0); - // int ncell(0); - // for ( auto cptr : *pCellCont ) { printf("CaloCell [%06i/%06i] ptr = %p\n",icell++,m_numberOfCells,(void*)cptr); if ( cptr == 0 ) { ++ncell; } } - // this->addMessage(CaloRec::Helpers::fmtMsg("Allocated CaloCellContainer with key <%s> at %p",m_cellContainerKey.c_str(),(void*)pCellCont)); - - // find input: tower container - const xAOD::CaloTowerContainer* pTowCont = 0; - CHECK(evtStore()->retrieve(pTowCont,m_towerContainerKey)); - // this->addMessage(CaloRec::Helpers::fmtMsg("Allocated CaloTowerContainer with key <%s> at %p",m_towerContainerKey.c_str(),(void*)pTowCont)); - - ////////////////////////////////////////////////////////////////////////////////////////// - // Tricky bit: // - // Especially for fine tower grids some towers may have 0 energy - in particular if the // - // grid spans -5.0 < eta < 5.0. The proto-cluster container has the same size as the // - // tower container, so it will have empty proto-clusters. These are removed when the // - // proto-clusters are converted to full CaloCluster objects. The tower center in // - // (eta,phi) is preserved in eta0(), phi0() of the CaloCluster. The cluster container // - // index has no relation to the tower position anymore. // - ////////////////////////////////////////////////////////////////////////////////////////// + if ( m_doCellIndexCheck ) { this->checkCellIndices(pCellCont.cptr()); } ///////////////////////// // Set up ProtoCluster // ///////////////////////// - protocont_t pProtoCont(pTowCont->size(),(CaloProtoCluster*)0); - for ( size_t itow(0); itowsize(); ++itow) { pProtoCont[itow] = new CaloProtoCluster(pCellCont); } + // index of CaloProtoCluster in container relates to tower position! DO NOT sort, shuffle, or remove elements! + protocont_t pProtoCont; pProtoCont.reserve(m_numberOfTowers); + for ( uint_t i(0); iretrieve(pTopoClusCont,m_clusterContainerKey)); // of the geometrical weight in cell splitting // - // instantiate cell-cluster link object // between clusters and the calibration weights. // - m_cellClusterWeights = new CaloCellClusterWeights(pCellCont->size()); // The geometrical weight of a cell in the towers is // - // fill the weights: cluster loop // used in the addCellToProtoCluster method. // - for ( auto pClus : *pTopoClusCont ) { /////////////////////////////////////////////////////// - // cell loop - for ( auto fCell(pClus->cell_begin()); fCell != pClus->cell_end(); ++fCell ) { - if ( !m_cellClusterWeights->check(*fCell) ) { this->addCellToProtoCluster(*fCell,pProtoCont); } // fill protocluster only once/cell - m_cellClusterWeights->set(*fCell,fCell.weight()); // prep for calibration tool - } // cell in cluster loop - } // cluster loop - // store LCW weights for calibrator - if ( m_prepareLCW ) { - CHECK(evtStore()->record(m_cellClusterWeights,CaloCellClusterWeights::key(m_cellClusterWeightKey))); - } else { - delete m_cellClusterWeights; m_cellClusterWeights = 0; - } + // The weights extracted for cells from clusters are LCW weights (typically). The total + // contribution of a LCW-weighted cell to towers is Ecell*Weight_LCW*Weight_tower. - ///////////////////////////////////////////////////////////// - // Generate and fill inclusive or exclusive proto-clusters // - ///////////////////////////////////////////////////////////// - - } else { - if ( !m_applyEnergyThreshold ) { buildInclClusters(*pCellCont,pProtoCont); } else { buildExclClusters(*pCellCont,pProtoCont); } - } + // If EM clusters are used, the weights of a clustered cell are completely defined + // by the tower grid. As cells can shared between clusters, each cell can only be + // projected onto the towergrid once, with Ecell*Weight_tower - //////////////////////////// - // Convert proto-clusters // - //////////////////////////// + // The CaloCellClusterWeights object is used to store the combined LCW weight. + // for each clustered cell. In case of EM, the LCW weights are ignored and this + // object is not used - a simple vector tags cells already put into towers. - pClusCont->reserve(pProtoCont.size()); - - xAOD::CaloCluster::ClusterSize csize = this->getClusterSize(*pTowCont); - int ictr(0); - for ( size_t ipc(0); ipcreleaseCellLinks(); // take over CaloClusterCellLink object - // convert proto-cluster to cluster - if ( lptr->size() > 0 ) { // remove empty proto-clusters - xAOD::CaloCluster* clptr = new xAOD::CaloCluster(); // new empty cluster - pClusCont->push_back(clptr); // put into container - clptr->addCellLink(lptr); // transfer cell links to CaloCluster - clptr->setClusterSize(csize); // set the cluster size spec - CaloRec::Helpers::calculateKine(clptr,false); // calculate kinematics and other signals from cells - clptr->setEta0(pTowCont->eta(ipc)); // save the tower center eta - clptr->setPhi0(pTowCont->phi(ipc)); // save the tower center phi - ++ictr; // counts the number of CaloCluster produced in the end - } // tower has cells - delete pProto; // delete the proto-cluster object - pProtoCont[ipc] = (CaloProtoCluster*)0; // FIXME (needed?) - } // valid proto-cluster + uint_t cCtr(0); + if ( m_useCellsFromClusters ) { + // retrieve topo-cluster container for topo-towers + SG::ReadHandle pTopoClusCont(m_clusterContainerKey); + if ( !pTopoClusCont.isValid() ) { + ATH_MSG_ERROR("Cannot allocate xAOD::CaloClusterContainer with key <" << m_clusterContainerKey << ">"); + return StatusCode::FAILURE; + } // check on ReadHandle validity + cCtr = m_prepareLCW ? this->buildLCWTopoTowers(*pTopoClusCont,pProtoCont) : this->buildEMTopoTowers(*pTopoClusCont,pProtoCont); + if ( !isValidIndex(cCtr) ) { ATH_MSG_WARNING("problems building EM or LCW topo-towers"); return StatusCode::SUCCESS; } + } else { + // fill inclusive/exclusive towers + cCtr = m_applyCellEnergyThreshold ? this->buildExclTowers(*pCellCont,pProtoCont) : this->buildInclTowers(*pCellCont,pProtoCont); + if ( !isValidIndex(cCtr) ) { ATH_MSG_WARNING("problems building EM inclusive or exclusive towers"); return StatusCode::SUCCESS; } + } // end topo-towers/inclusive-exclusive towers + + // allocate sufficient space in vector + pClusCont->reserve(cCtr); + // pick up cluster size tag and set up counter + xAOD::CaloCluster::ClusterSize csize = this->getClusterSize(m_numberOfTowers); + // loop proto-clusters + for ( uint_t ipc(0); ipccleanupCells(lptr,ipc); // clean up cell links + if ( this->filterProtoCluster(*lptr) ) { // ignore empty proto-clusters (no cells assigned) + xAOD::CaloCluster* clptr = new xAOD::CaloCluster(); // new empty cluster + pClusCont->push_back(clptr); // put into container + clptr->addCellLink(lptr); // transfer cell links to CaloCluster + clptr->setClusterSize(csize); // set the cluster size spec + CaloRec::Helpers::calculateKine(clptr,false); // calculate kinematics and other signals from cells + clptr->setEta0(m_towerGeometrySvc->towerEta(ipc)); // save the tower center eta + clptr->setPhi0(m_towerGeometrySvc->towerPhi(ipc)); // save the tower center phi + } else { + delete lptr; + } } // proto-cluster loop - pProtoCont.clear(); // clean up (objects in this container have already been destroyed) - - //////////////////// - // Check clusters // - //////////////////// - - size_t idc(0); - for ( auto pClus : *pClusCont ) { - if ( pClus->getCellLinks() == 0 || pClus->getCellLinks()->size() == 0 ) { - this->addWarning(CaloRec::Helpers::fmtMsg("Cluster #%05i has invalid pointer (%p) or no cells",idc,(void*)pClus)); - } - ++idc; - } + // clean up proto-cluster container + pProtoCont.clear(); - //////////////////////// - // Sort protoclusters // - //////////////////////// + ///////////// + // Sorting // + ///////////// - // definitively removes the index-to-tower center relation in all cases! + // All towers/clusters at this point are on EM scale. Sorting LCW towers by pT should be done in the + // CaloTopoClusterFromTowerCalibrator tool to assure desired ordering on the final scale. + // The link between tower location and index of tower representation (CaloCluster) in its + // container is definitively broken after sorting (was never ok in mixed cluster/tower mode). if ( m_orderByPt ) { std::sort(pClusCont->begin(),pClusCont->end(),[](xAOD::CaloCluster* pc1,xAOD::CaloCluster* pc2) { volatile double pt1(pc1->pt()); // FIXME needed? (this was just copied) @@ -294,226 +298,317 @@ StatusCode CaloTopoClusterFromTowerMaker::execute(xAOD::CaloClusterContainer* pC return ( pt1 > pt2 ); } ); + } // end ordered by pT + + return StatusCode::SUCCESS; +} // end execute + +////////////////////// +// Fill topo-towers // +////////////////////// + +// EM +CaloTopoClusterFromTowerMaker::uint_t CaloTopoClusterFromTowerMaker::buildEMTopoTowers(const xAOD::CaloClusterContainer& pClusCont,protocont_t& pProtoCont) +{ + // presets + uint_t cCtr(0); + std::vector cellTags(m_numberOfCells,false); + + // -- EM scale clusters + if ( !m_buildCombinedSignal ) { + // topo-towers + for ( auto pClus : pClusCont ) { + for ( auto fCell(pClus->cell_begin()); fCell != pClus->cell_end(); ++fCell ) { + uint_t cidx(static_cast((*fCell)->caloDDE()->calo_hash())); + if ( cidx < cellTags.size() ) { + if ( !cellTags.at(cidx) ) { cellTags[cidx] = this->addCellToProtoCluster(*fCell,pProtoCont); } + } else { + ATH_MSG_ERROR( CaloRec::Helpers::fmtMsg("Invalid cell hash index %6zu >= maximum index %6zu for cell in %s at (eta,phi) = (%6.3,%f6.3)", + cidx,cellTags.size(),CaloSampling::getSamplingName((*fCell)->caloDDE()->getSampling()).c_str(),(*fCell)->eta(),(*fCell)->phi()) ); + return m_errorValueUINT; + } + } // end cells-in-cluster loop + } // end cluster loop + } else { + // selected topo-towers for combined signal + std::vector > cellList(m_numberOfCells,std::tuple(0,0.)); + for ( auto pClus : pClusCont ) { + if ( std::abs(pClus->eta()) > m_clusterRange ) { + for ( auto fCell(pClus->cell_begin()); fCell != pClus->cell_end(); ++fCell ) { + uint_t cidx(static_cast((*fCell)->caloDDE()->calo_hash())); + if ( cellTags.at(cidx) ) { + std::get<1>(cellList[cidx]) += fCell.weight(); + } else { + cellList[cidx] = std::tuple(*fCell,fCell.weight()); + cellTags[cidx] = true; + } + } // cell in cluster loop + } else { + ++cCtr; + } // cluster range check + } // cluster loop + // fill proto-cluster + for ( auto tpl : cellList ) { this->addCellToProtoCluster(std::get<0>(tpl),pProtoCont,std::get<1>(tpl)); } + } // end of fill mode + + // + return cCtr+pProtoCont.size(); +} + +// LCW +CaloTopoClusterFromTowerMaker::uint_t CaloTopoClusterFromTowerMaker::buildLCWTopoTowers(const xAOD::CaloClusterContainer& pClusCont,protocont_t& pProtoCont) +{ + // Need to keep track of LCW weights (up to two per cell) available from the topo-cluster(s) the cell is assigned to. + // Each cell in a topo-cluster is, at first occurance, added to the CaloProtoCluster(s) representing the tower(s) and its + // LCW calibration weight is stored in a lookup table indexed by the calorimeter hash id of the cell. The second + // time the same cell is found in another topo-cluster, only its LCW weight is added to the lookup table (stored in + // CaloCellClusterWeights for use in the downstream tower calibration tool) - the assignment to tower(s) has already + // happened. + + uint_t cCtr(0); + // write handle object on the stack + SG::WriteHandle cellClusterWeightHandle(m_cellClusterWeightKey); + // record output container + if ( cellClusterWeightHandle.record(std::make_unique(m_towerGeometrySvc->maxCellHash())).isFailure() ) { + ATH_MSG_WARNING( "Cannot record cluster cell weights with key <" << m_cellClusterWeightKey << ">" ); + return m_errorValueUINT; } + // project cells on tower grid + if ( !m_buildCombinedSignal ) { + for ( auto pClus : pClusCont ) { + for ( auto fCell(pClus->cell_begin()); fCell != pClus->cell_end(); ++fCell ) { + // map to towers only once + if ( !cellClusterWeightHandle.ptr()->check(*fCell) ) { this->addCellToProtoCluster(*fCell,pProtoCont); } + // store all associated LCW weights + cellClusterWeightHandle.ptr()->set(*fCell,fCell.weight()); + } // end cells-in-cluster loop + } // end cluster loop + } else { + for ( auto pClus : pClusCont ) { + if ( std::abs(pClus->eta()) > m_clusterRange ) { + for ( auto fCell(pClus->cell_begin()); fCell != pClus->cell_end(); ++fCell ) { + // map to towers only once + if ( !cellClusterWeightHandle.ptr()->check(*fCell) ) { this->addCellToProtoCluster(*fCell,pProtoCont); } + // store all associated LCW weights + cellClusterWeightHandle.ptr()->set(*fCell,fCell.weight()); + } // end cells-in-cluster loop + } else { + ++cCtr; + } // end range check + } // end cluster loop + } // end combined signal check + + // + return cCtr+pProtoCont.size(); +} + +///////////////// +// Fill towers // +///////////////// +// inclusive +CaloTopoClusterFromTowerMaker::uint_t CaloTopoClusterFromTowerMaker::buildInclTowers(const CaloCellContainer& pCellCont,protocont_t& pProtoCont) +{ + // loop cell container - counter icl replaces cell hash index for NULL pointers in cell container + uint_t icl(0); + for ( auto cptr : pCellCont ) { + if ( cptr == 0 ) { + ATH_MSG_ERROR( CaloRec::Helpers::fmtMsg("CaloCellContainer[%6zu] contains invalid cell object pointer %p",icl,(void*)cptr) ); + return m_errorValueUINT; + } else { + // existing cell with non-zero energy (negative or positive) + if ( std::fabs(cptr->e()) > 0. ) { this->addCellToProtoCluster(cptr,pProtoCont); } + } // end pointer check + ++icl; + } // end cell loop + return pProtoCont.size(); +} + +// exclusive +CaloTopoClusterFromTowerMaker::uint_t CaloTopoClusterFromTowerMaker::buildExclTowers(const CaloCellContainer& pCellCont,protocont_t& pProtoCont) +{ + // loop cell container + uint_t icl(0); + for ( auto cptr : pCellCont ) { + if ( cptr == 0 ) { + ATH_MSG_ERROR( CaloRec::Helpers::fmtMsg("CaloCellContainer[%6zu] contains invalid cell object pointer %p",icl,(void*)cptr) ); + return m_errorValueUINT; + } else { + // existing cell with energy above threshold + if ( cptr->e() > m_energyThreshold ) { this->addCellToProtoCluster(cptr,pProtoCont); } + } // end pointer check + ++icl; + } // end cell loop return StatusCode::SUCCESS; } -xAOD::CaloCluster::ClusterSize CaloTopoClusterFromTowerMaker::getClusterSize(const xAOD::CaloTowerContainer& towerCont) +bool CaloTopoClusterFromTowerMaker::addCellToProtoCluster(const CaloCell* cptr,protocont_t& pProtoCont,double weight) { - return towerCont.nTowers() == 6400 - ? xAOD::CaloCluster::Tower_01_01 : towerCont.nTowers() == 25600 - ? xAOD::CaloCluster::Tower_005_005 : xAOD::CaloCluster::CSize_Unknown; + // invalid input + if ( cptr == 0 ) { return false; } + + // get towers for cell from geometry service + uint_t nctr(0); + for ( auto elm : m_towerGeometrySvc->getTowers(cptr) ) { + auto towerIdx(m_towerGeometrySvc->towerIndex(elm)); + if ( !m_towerGeometrySvc->isInvalidIndex(towerIdx) ) { + if ( !m_excludedSamplingsPattern[(size_t)cptr->caloDDE()->getSampling()] ) { + uint_t cellIdx(pProtoCont.at(towerIdx).getCellLinks()->getCellContainer()->findIndex(cptr->caloDDE()->calo_hash())); + pProtoCont[towerIdx].addCell(cellIdx,m_towerGeometrySvc->cellWeight(elm)*weight); ++nctr; + } + } + } + return nctr > 0; +} + +///////////// +// Helpers // +///////////// + +xAOD::CaloCluster::ClusterSize CaloTopoClusterFromTowerMaker::getClusterSize(uint_t etaBins,uint_t phiBins) +{ return this->getClusterSize(etaBins*phiBins); } + +xAOD::CaloCluster::ClusterSize CaloTopoClusterFromTowerMaker::getClusterSize(uint_t nTowers) +{ + // check for tower sizes + return nTowers == 6400 // known "standard" towers 0,1 x 0.1 + ? xAOD::CaloCluster::Tower_01_01 + : nTowers == 25600 // known "fine" towers 0.05 x 0.05 + ? xAOD::CaloCluster::Tower_005_005 + : xAOD::CaloCluster::Tower_fixed_area; // unspecified towers } -int CaloTopoClusterFromTowerMaker::cleanupCells(CaloClusterCellLink* clk) +int CaloTopoClusterFromTowerMaker::cleanupCells(CaloClusterCellLink* clk,uint_t nclus) { + // Any pathology here probably indicates a configuration problem with the conditions (geometry) + // database (wrong tag for data?) + // check on null pointers in cell links int nrc(0); int hid(0); auto fcell(clk->begin()); while ( fcell != clk->end() ) { const CaloCell* pCell = *fcell; - this->addMessage(CaloRec::Helpers::fmtMsg("CaloCell* = %p in CaloClusterCellLink* = %p")); - bool removeCell(false); - if ( pCell == 0 ) { - m_cellProblems.setInvalidPointer(hid,pCell); - removeCell = true; + auto nc(clk->getCellContainer()->size()); + const CaloCell* aCell = fcell.index() < nc ? clk->getCellContainer()->at(fcell.index()) : (const CaloCell*)0; + if ( pCell == 0 ) { + ATH_MSG_WARNING( CaloRec::Helpers::fmtMsg("CaloCellContainer[%6zu/%6zu] - tower %5zu at (%6.3f,%6.3f) - cell pointer invalid (%p/%p) [removed %3i of %3zu cells]", + fcell.index(),nc-1,nclus,m_towerGeometrySvc->towerEta(nclus),m_towerGeometrySvc->towerPhi(nclus), + (void*)pCell,(void*)aCell,++nrc,clk->size()) ); + fcell = clk->removeCell(fcell); } else { - int hash((int)pCell->caloDDE()->calo_hash()); - if ( hash >= m_numberOfCells ) { m_cellProblems.setInvalidHash(hash,pCell); removeCell = true; } - int samp((int)pCell->caloDDE()->getSampling()); - if ( samp >= m_numberOfSamplings ) { m_cellProblems.setInvalidSampling(hash,pCell); removeCell = true; } - } - if ( removeCell ) { fcell = clk->removeCell(fcell); ++nrc; } else { ++fcell; } + uint_t chash(static_cast(pCell->caloDDE()->calo_hash())); + uint_t csamp(static_cast(pCell->caloDDE()->getSampling())); + if (chash > m_maxCellHash ) { + // check cell hash + ATH_MSG_WARNING( CaloRec::Helpers::fmtMsg("Tower %5zu at (%6.3f,%6.3f) linked cell %3i - cell hash index (%6zu/%6zu) invalid", + nclus,m_towerGeometrySvc->towerEta(nclus),m_towerGeometrySvc->towerPhi(nclus),hid,chash,m_maxCellHash) ); + fcell = clk->removeCell(fcell); ++nrc; + } else if ( csamp >= m_numberOfSamplings ) { + // check sampling id + ATH_MSG_WARNING( CaloRec::Helpers::fmtMsg("Tower %5zu at (%6.3f,%6.3f) linked cell %3i -cell sampling id (%3zu/%3zu) invalid", + nclus,m_towerGeometrySvc->towerEta(nclus),m_towerGeometrySvc->towerPhi(nclus),hid,csamp,m_numberOfSamplings) ); + fcell = clk->removeCell(fcell); ++nrc; + } else if ( fcell.weight() <= 0.0000001 ) { + // remove cells with 0 weight + fcell = clk->removeCell(fcell); ++nrc; + } else { + // next cell + ++fcell; + } + } // end remove cell due to pointer invalid ++hid; - } + } // end loop on cells in cell link object return nrc; } -// bool CaloTopoClusterFromTowerMaker::calculateKine(xAOD::CaloCluster* pClus) -// { return CaloRec::Helpers::updateKine -// // input -// if ( pClus == 0 ) { -// ATH_MSG_WARNING("Invalid (NULL) pointer to xAOD::CaloCluster object"); -// return false; -// } -// // copy from CaloKineHelper::calculateKine(...) - but without the prefetch! -// const CaloClusterCellLink* clk = pClus->getCellLinks(); -// if ( clk == 0 ) { -// this->addWarning(CaloRec::Helpers::fmtMsg("Invalid reference (pointer %p) to CaloClusterCellLink object",(void*)clk)); -// return false; -// } -// if ( clk->size() == 0 ) { -// this->addWarning(CaloRec::Helpers::fmtMsg("Number of linked cells in cluster is %i",(int)clk->size())); -// return false; -// } - -// // accumulator object -// CaloRec::Helpers::CaloClusterSignalAccumulator accum; -// // accumulate cells -// for ( auto cIter(clk->begin()); cIter!=clk->end(); ++cIter) { -// const CaloCell* pCell = *cIter; -// if ( pCell != 0 ) { -// if ( !CaloRec::Helpers::caloCellAccumulator(*pCell,accum,cIter.weight(),false) ) { -// this->addWarning(CaloRec::Helpers::fmtMsg("CaloCell [%06i] has invalid sampling id %2i (not in [0,%2i])",(int)(*cIter)->caloDDE()->calo_hash(),(int)(*cIter)->caloDDE()->getSampling(),(int)CaloSampling::Unknown)); -// } -// } else { -// this->addWarning(CaloRec::Helpers::fmtMsg("Unexpected CaloCell pointer %p detected",(void*)pCell)); -// } -// } - -// // set cluster kinematics: energy & mass -// pClus->setE(accum.cluster.accumE); -// pClus->setM(0.); - -// // set cluster kinematics: directions -// if ( accum.cluster.accumAbsE != 0. ) { -// double invPosNorm(1./accum.cluster.accumAbsE); -// pClus->setEta(accum.cluster.accumEta*invPosNorm); -// pClus->setPhi(CaloPhiRange::fix(accum.cluster.accumPhi*invPosNorm)); -// } else { -// pClus->setEta(0.); -// pClus->setPhi(0.); -// } - -// // set cluster kinematics: time -// if ( accum.cluster.accumTimeNorm != 0. ) { -// pClus->setTime(accum.cluster.accumTime/accum.cluster.accumTimeNorm); -// } else { -// pClus->setTime(0.); -// } - -// // set sampling pattern -// uint32_t samplingPattern(0); -// for ( int i(0); i<(int)CaloSampling::Unknown; ++i ) { -// if ( accum.sampling.presenceInSample[i] ) { samplingPattern |= (0x1U<samplingPattern() ) { -// pClus->clearSamplingData(); -// pClus->setSamplingPattern(samplingPattern,true); -// } -// } - -// // set sampling variables -// for ( int i(0); i<(int)CaloSampling::Unknown; ++i ) { -// if ( accum.presenceInSample[i] ) { -// CaloSampling::CaloSample sam = (CaloSampling::CaloSample)i; -// pClus->setEnergy(sam,accum.energyInSample[i]); -// if ( accum.posNormInSample[i] != 0. ) { -// pClus->setEta(sam,accum.etaInSample[i]/accum.posNormInSample[i]); -// pClus->setPhi(sam,accum.phiInSample[i]/accum.posNormInSample[i]); -// } else { -// pClus->setEta(sam,accum.etaInSample[i]); -// pClus->setPhi(sam,accum.phiInSample[i]); -// } -// pClus->setEmax(sam,accum.maxEnergyInSample[i]); -// pClus->setEtamax(sam,accum.etaMaxEnergyInSample[i]); -// pClus->setPhimax(sam,accum.phiMaxEnergyInSample[i]); -// } // check if sampling is in cluster -// } // loop on samplings - -// return true; -// } - -void CaloTopoClusterFromTowerMaker::addWarning(const std::string& msg) -{ - ATH_MSG_WARNING(msg); - auto fiter(m_warning.find(msg)); - if ( fiter == m_warning.end() ) { m_warning[msg] = 0; } - ++m_warning[msg]; -} +bool CaloTopoClusterFromTowerMaker::filterProtoCluster(const CaloClusterCellLink& clnk) const +{ return clnk.size() > 0; } -void CaloTopoClusterFromTowerMaker::addMessage(const std::string& msg) +bool CaloTopoClusterFromTowerMaker::checkCellIndices(const CaloCellContainer* pCellCont) const { - ATH_MSG_DEBUG(msg); - auto fiter(m_message.find(msg)); - if ( fiter == m_message.end() ) { m_message[msg] = 0; } - ++m_message[msg]; -} + //////////////////////////// + // input and setup checks // + //////////////////////////// -void CaloTopoClusterFromTowerMaker::printProblems(const std::map& map) -{ - static std::string _dstr = "...................................................................."; - static size_t _lds = _dstr.length(); - static size_t _occ = _lds - 3; - // - for ( auto fmap(map.begin()); fmap!=map.end(); ++fmap ) { - std::string msg = fmap->first.length() > _occ ? fmap->first.substr(0,_occ) : fmap->first; - std::string prn(std::string(_dstr).replace(0,msg.length(),msg)); - ATH_MSG_INFO(CaloRec::Helpers::fmtMsg("%s %i",prn.c_str(),fmap->second)); + // check argument + if ( pCellCont == 0 ) { + ATH_MSG_WARNING( CaloRec::Helpers::fmtMsg("Invalid pointer to CaloCellContainer (%p)",(void*)pCellCont) ); return false; + } else if ( pCellCont->empty() ) { + ATH_MSG_WARNING( CaloRec::Helpers::fmtMsg("CaloCellContainer at %p is empty (size %zu)",(void*)pCellCont,pCellCont->size()) ); return false; } -} + // check the atomic state + if ( CaloTopoClusterFromTowerMaker_checkCellIndices ) { return true; } + // set the atomic flag + ATH_MSG_INFO( "Cell hash index check requested" ); + CaloTopoClusterFromTowerMaker_checkCellIndices = true; + // assign output file + std::string algname(this->name()); + if ( algname.find_last_of('.') != std::string::npos ) { algname = algname.substr(algname.find_last_of('.')+1); } + std::string logname(CaloRec::Helpers::fmtMsg("%s.cellhash_index_check.dat",this->name().c_str())); + std::ofstream logstream; logstream.open(logname); + if ( !logstream.is_open() ) { + ATH_MSG_WARNING( CaloRec::Helpers::fmtMsg("Cannot open log file \042%s\042 - no hash index checking",logname.c_str()) ); + return false; + } + logstream << "##########################################################################" << std::endl; + logstream << "### This file contains a list of CaloCell indices in CaloCellContainer ###" << std::endl; + logstream << "### for which this index is not the same as the calorimeter cell hash ###" << std::endl; + logstream << "### identifier. An empty list indicates full consistency between this ###" << std::endl; + logstream << "### index and the hash identifier for all cells. ###" << std::endl; + logstream << "##########################################################################" << std::endl; + logstream << "--------------------------------------------------------------" << std::endl; -bool CaloTopoClusterFromTowerMaker::buildInclClusters(const CaloCellContainer& pCellCont,protocont_t& pProtoCont) -{ - // inclusive mode - this->addMessage(CaloRec::Helpers::fmtMsg("Found %6i cell links in CaloCellTowerLink::Map::%s",(int)m_cellTowerLinksHandle->size(),m_cellTowerLinksKey.c_str())); + ///////////////////////// + // loop cell container // + ///////////////////////// - // loop cell container - size_t icl(0); - for ( auto cptr : pCellCont ) { - // FIXME check pointer value - if ( cptr == 0 ) { m_cellProblems.setInvalidPointer(icl,cptr); } - // existing cell - else { - if ( std::fabs(cptr->e()) > 0. ) { this->addCellToProtoCluster(cptr,pProtoCont); } + // prepare tag store + size_t ifc(0); std::bitset<200000> chkflg; chkflg.reset(); + for ( size_t i(0); isize(); ++i ) { + if ( pCellCont->at(i) != 0 ) { + size_t chash((size_t)pCellCont->at(i)->caloDDE()->calo_hash()); + if ( chash != i ) { + std::string cni("UKNOWN"); + double etai(0.); double phii(0.); + const CaloDetDescrElement* iel = i < CaloDetDescrManager::instance()->element_size() ? CaloDetDescrManager::instance()->get_element(i) : 0; + if ( iel != 0 ) { + cni = CaloRec::Lookup::getSamplingName(iel->getSampling()); + etai = iel->eta_raw(); + phii = iel->phi_raw(); + } + std::string cnc("UNKNOWN"); + double etac(0.); double phic(0.); + const CaloDetDescrElement* cel = chash < CaloDetDescrManager::instance()->element_size() ? CaloDetDescrManager::instance()->get_element(chash) : 0; + if ( cel != 0 ) { + cnc = CaloRec::Lookup::getSamplingName(cel->getSampling()); + etac = cel->eta_raw(); + phic = cel->phi_raw(); + } + size_t cidx(pCellCont->findIndex(chash)); + logstream << CaloRec::Helpers::fmtMsg("[%06zu] Cell %6zu [%12.12s %5.3f %5.3f] non-matching id %6zu [%12.12s %5.3f %5.3f] findCell() index %6zu", + ++ifc,i,cni.c_str(),etai,phii,chash,cnc.c_str(),etac,phic,cidx) << std::endl; + } + chkflg.set(chash); } - ++icl; - } // cell loop - return true; -} - -bool CaloTopoClusterFromTowerMaker::buildExclClusters(const CaloCellContainer& pCellCont,protocont_t& pProtoCont) -{ - // inclusive mode - this->addMessage(CaloRec::Helpers::fmtMsg("Found %6i cell links in CaloCellTowerLink::Map::%s",(int)m_cellTowerLinksHandle->size(),m_cellTowerLinksKey.c_str())); + } + logstream << "----------------------------------------------------------------" << std::endl; + logstream.close(); - // loop cell container - for ( auto cptr : pCellCont ) { - // FIXME check pointer value - if ( cptr == 0 ) { m_cellProblems.setInvalidPointer(999999,cptr); } - // existing cell - else { - if ( cptr->e() > m_energyThreshold ) { this->addCellToProtoCluster(cptr,pProtoCont); } - } - } // cell loop - return true; -} + ///////////////////////// + // check missed hashes // + ///////////////////////// -// bool CaloTopoClusterFromTowerMaker::buildFilteredClusters(const CaloCellContainer& pCellCont,protocont_t& pProtoCont) -// { -// // inclusive mode -// this->addMessage(CaloRec::Helpers::fmtMsg("Found %6i cell links in CaloCellTowerLink::Map::%s",(int)m_cellTowerLinksHandle->size(),m_cellTowerLinksKey.c_str())); - -// // loop cell container -// size_t icl(0); -// for ( auto cptr : pCellCont ) { -// // FIXME check pointer value -// if ( cptr == 0 ) { m_cellProblems.setInvalidPointer(icl,cptr); } -// // existing cell -// else { -// size_t idx(static_cast(cptr->caloDDE()->calo_hash_id())); -// if ( m_cellTagStore[idx].get<0>() ) { if ( std::fabs(cptr->e()) > 0. ) { this->addCellToProtoCluster(cptr,pProtoCont); } } -// } -// ++icl; -// } // cell loop -// return true; -// } - -bool CaloTopoClusterFromTowerMaker::addCellToProtoCluster(const CaloCell* cptr,protocont_t& pProtoCont) -{ - if ( cptr == 0 ) { return false; } - // cell hash - int cellidx(cptr->caloDDE()->calo_hash()); - if ( cellidx >= (int)m_cellTowerLinksHandle->size() ) { - this->addWarning(CaloRec::Helpers::fmtMsg("CaloCell [%06i] has hash outside of link lookup (max hash is %i)",cellidx,(int)m_cellTowerLinksHandle->size())); - return false; + // number of non-matched hashes + if ( ifc > 0 ) { + ATH_MSG_DEBUG( CaloRec::Helpers::fmtMsg("Found %zu non-matching cell hashes",ifc) ); + } + // list of non-matched hashes + std::vector chl; chl.reserve(m_towerGeometrySvc->totalNumberCells()); + for ( size_t i(0); ielementList(cellidx) ) { - int towidx(CaloCellTowerLink::towerIndex(link)); - if ( towidx < (int)pProtoCont.size() ) { pProtoCont[towidx]->addCell(cellidx,CaloCellTowerLink::cellWeight(link)); } - else { this->addWarning(CaloRec::Helpers::fmtMsg("Invalid tower index %i (not in [0,%i])",(int)towidx,(int)pProtoCont.size()-1)); } - } // loop all cell->tower links - } // valid hash index + return true; } + diff --git a/Calorimeter/CaloRec/src/CaloTopoClusterFromTowerMonitor.cxx b/Calorimeter/CaloRec/src/CaloTopoClusterFromTowerMonitor.cxx new file mode 100644 index 0000000000000000000000000000000000000000..ad12d4b149b8279ef12c6b55620e99888c6a3c4b --- /dev/null +++ b/Calorimeter/CaloRec/src/CaloTopoClusterFromTowerMonitor.cxx @@ -0,0 +1,421 @@ + +#include "GaudiKernel/MsgStream.h" +#include "GaudiKernel/SystemOfUnits.h" +#include "GaudiKernel/PhysicalConstants.h" + +#include "CaloRec/CaloTopoClusterFromTowerMonitor.h" + +#include "CaloGeoHelpers/CaloPhiRange.h" + +#include "CaloDetDescr/CaloDetDescrElement.h" +#include "CaloDetDescr/CaloDetDescrManager.h" + +#include "TH1D.h" +#include "TH2D.h" + +#include + +namespace { + MsgStream& operator<<(MsgStream& mstr,const SG::ReadHandleKey& ckey) { mstr << ckey.key(); return mstr; } + MsgStream& operator<<(MsgStream& mstr,const ServiceHandle& shdl) { mstr << shdl->name(); return mstr; } +} + +CaloTopoClusterFromTowerMonitor::CaloTopoClusterFromTowerMonitor(const std::string& name,ISvcLocator* pSvcLocator) + : AthHistogramAlgorithm(name,pSvcLocator) + , m_towerContainerKey("LCWTowerTopoClusterStd") + , m_towerGeometrySvc("CaloTowerGeometrySvc",name) + , m_ncBins(100) + , m_ncMin(-0.5) + , m_ncMax(m_ncMin+m_ncBins) + , m_nBins(320) + , m_nMin(-0.05) + , m_nMax(m_nMin+m_nBins*20.) + , m_ptBins(220) + , m_ptMin(-10.) + , m_ptMax(100.) + , m_etaBins(100) + , m_etaMin(-5.) + , m_etaMax(5.) + , m_phiBins(64) + , m_phiMin(-Gaudi::Units::pi) + , m_phiMax(Gaudi::Units::pi) + , m_hsEtaMin(0.) + , m_hsEtaMax(0.1) + , m_hsPhiMin(0.) + , m_hsPhiMax(0.1) + , m_doGeoAutoBins(true) + , m_doHotspot(false) + , h_n((TH1D*)0), h_pt((TH1D*)0), h_eta((TH1D*)0), h_phi((TH1D*)0), h_nc((TH1D*)0), h_samp((TH1D*)0) + , d_n_eta_phi((TH2D*)0), d_nc_eta_phi((TH2D*)0) + , d_pt_eta((TH2D*)0), d_nc_eta((TH2D*)0) + , d_wgt_samp((TH2D*)0), d_ntt_samp((TH2D*)0), d_geo_samp((TH2D*)0) + , d_maxtowers_samp((TH2D*)0), d_wgttowers_samp((TH2D*)0) + , d_maxcells_eta((TH2D*)0), d_allwghts_eta((TH2D*)0) + , d_deta_eta((TH2D*)0), d_dphi_eta((TH2D*)0), d_dphi_deta((TH2D*)0) + , d_detac_eta((TH2D*)0), d_dphic_eta((TH2D*)0), d_dphic_detac((TH2D*)0), d_detac_samp((TH2D*)0), d_dphic_samp((TH2D*)0) + , h_nc_hs((TH1D*)0), h_n_hs((TH1D*)0), h_pt_hs((TH1D*)0), h_eta_hs((TH1D*)0), h_phi_hs((TH1D*)0), h_samp_hs((TH1D*)0) + , d_n_eta_phi_hs((TH2D*)0), d_nc_eta_phi_hs((TH2D*)0) + , d_deta_eta_hs((TH2D*)0), d_dphi_eta_hs((TH2D*)0), d_dphi_deta_hs((TH2D*)0) + , d_detac_eta_hs((TH2D*)0), d_dphic_eta_hs((TH2D*)0), d_dphic_detac_hs((TH2D*)0), d_detac_samp_hs((TH2D*)0), d_dphic_samp_hs((TH2D*)0) +{ + declareProperty("CaloTowerContainerKey", m_towerContainerKey, "Input container key" ); + declareProperty("CaloTowerGeometrySvc", m_towerGeometrySvc, "Tower geometry provider" ); + declareProperty("NTowerCellsBins", m_ncBins, "Number of bins in cell-in-tower multiplicity"); + declareProperty("NTowerCellsMin", m_ncMin, "Lower limit in cell-in-tower multiplicity"); + declareProperty("NTowerCellsMax", m_ncMax, "Upper limit in cell-in-tower multiplicity"); + declareProperty("NTowersBins", m_nBins, "Number of bins in tower multiplicity binning"); + declareProperty("NTowersMin", m_nMin, "Lower limit in tower multiplicity binning"); + declareProperty("NTowersMax", m_nMax, "Upper limit in tower multiplicity binning"); + declareProperty("PtTowersBins", m_ptBins, "Number of bins in tower pT binning"); + declareProperty("PtTowersMin", m_ptMin, "Lower limit in tower pT binning"); + declareProperty("PtTowersMax", m_ptMax, "Upper limit in tower pT binning"); + declareProperty("EtaTowersBins", m_etaBins, "Number of bins in tower rapidity binning"); + declareProperty("EtaTowersMin", m_etaMin, "Lower limit in tower rapidity binning"); + declareProperty("EtaTowersMax", m_etaMax, "Upper limit in tower rapidity binning"); + declareProperty("PhiTowersBins", m_phiBins, "Number of bins in tower azimuth binning"); + declareProperty("PhiTowersMin", m_phiMin, "Lower limit in tower azimuth binning"); + declareProperty("PhiTowersMax", m_phiMax, "Upper limit in tower azimuth binning"); + declareProperty("EtaMinHotspot", m_hsEtaMin, "lower limit in tower rapidity for hotspot"); + declareProperty("EtaMaxHotspot", m_hsEtaMax, "Upper limit in tower rapidity for hotspot"); + declareProperty("PhiMinHotspot", m_hsPhiMin, "lower limit in tower azimuth for hotspot"); + declareProperty("PhiMaxHotspot", m_hsPhiMax, "Upper limit in tower azimuth for hotspot"); + declareProperty("DoGeoAutoBins", m_doGeoAutoBins, "Flag controls automatic binning for rapidity and azimuth"); +} + +CaloTopoClusterFromTowerMonitor::~CaloTopoClusterFromTowerMonitor() +{ } + +StatusCode CaloTopoClusterFromTowerMonitor::initialize() +{ + // initialize read handle key + ATH_CHECK(m_towerContainerKey.initialize()); + + // tower geometry service + ATH_MSG_INFO( "Allocate tower geometry service:" ); + if ( !m_towerGeometrySvc.isValid() ) { + ATH_MSG_ERROR("[reject] - cannot allocate tower geometry service - fatal"); + return StatusCode::FAILURE; + } else { + ATH_MSG_INFO( "[accept] - allocated tower geometry provider \042" << m_towerGeometrySvc << "\042"); + } + + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("Tower geometry service is allocated, describes %6zu towers in grid:", m_towerGeometrySvc->towerBins()) ); + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("[accept] %3zu eta bins in [%5.2f,%5.2f]",m_towerGeometrySvc->etaBins(),m_towerGeometrySvc->etaMin(),m_towerGeometrySvc->etaMax()) ); + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("[accept] %3zu phi bins in [%5.2f,%5.2f]",m_towerGeometrySvc->phiBins(),m_towerGeometrySvc->phiMin(),m_towerGeometrySvc->phiMax()) ); + + std::string cfgStr; + if ( m_doGeoAutoBins ) { + ATH_MSG_INFO("Eta and phi binning taken from tower geometry provider"); + m_etaBins = m_towerGeometrySvc->etaBins(); + m_etaMin = m_towerGeometrySvc->etaMin(); + m_etaMax = m_towerGeometrySvc->etaMax(); + m_phiBins = m_towerGeometrySvc->phiBins(); + m_phiMin = m_towerGeometrySvc->phiMin(); + m_phiMax = m_towerGeometrySvc->phiMax(); + m_hsPhiMax = m_hsPhiMin+m_towerGeometrySvc->phiWidth(); + cfgStr = "autoconfig"; + } else { + ATH_MSG_INFO("Eta and phi binning taken from specific configuration"); + cfgStr = "userconfig"; + } + + m_doHotspot = m_hsEtaMin < m_hsEtaMax && m_hsPhiMin < m_hsPhiMax; + + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("Eta and phi binning for histograms (total %zu towers)",m_etaBins*m_phiBins) ); + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("[%s] %3zu eta bins in [%5.2f,%5.2f]",cfgStr.c_str(),m_etaBins,m_etaMin,m_etaMax) ); + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("[%s] %3zu phi bins in [%5.2f,%5.2f]",cfgStr.c_str(),m_phiBins,m_phiMin,m_phiMax) ); + + ATH_CHECK( book() ); + + // database plots + std::vector > cWghts(m_towerGeometrySvc->towerBins(),std::vector()); + size_t chash(0); + for ( auto fLnk(m_towerGeometrySvc->begin()); fLnk != m_towerGeometrySvc->end(); ++fLnk, ++chash ) { + CaloSampling::CaloSample samp = CaloDetDescrManager::instance()->get_element(chash)->getSampling(); + double cs(static_cast(samp)); + double nt(static_cast(fLnk->size())); + d_maxtowers_samp->Fill(cs,nt); + for ( auto elm : *fLnk ) { + // collect for towers + size_t tidx(static_cast(m_towerGeometrySvc->towerIndex(elm))); + double weight(m_towerGeometrySvc->cellWeight(elm)); + cWghts[tidx].push_back(weight); + // plot for smaplings + d_wgttowers_samp->Fill(cs,weight); + } + } + + for ( size_t tidx(0); tidxtowerEta(tidx)); + double phi(m_towerGeometrySvc->towerPhi(tidx)); + double nc(1.0*cWghts.at(tidx).size()); + d_maxcells_eta->Fill(eta,nc); + size_t etaIdx(m_towerGeometrySvc->etaIndexFromTowerIndex(tidx)); + d_maxcells_phi_eta_slice[etaIdx]->Fill(phi,nc); + for ( auto w : cWghts.at(tidx) ) { d_allwghts_eta->Fill(eta,w); d_allwghts_phi_eta_slice[etaIdx]->Fill(phi,w); } + } + + if ( msgLvl(MSG::DEBUG) ) { + ATH_MSG_DEBUG( CaloRec::Helpers::fmtMsg("+------------+--------------+") ); + ATH_MSG_DEBUG( CaloRec::Helpers::fmtMsg("| SamplingId | SamplingName |") ); + ATH_MSG_DEBUG( CaloRec::Helpers::fmtMsg("+------------+--------------+") ); + const auto& smap = CaloRec::Lookup::samplingNames; + for ( auto fsamp(smap.begin()); fsamp != smap.end(); ++fsamp ) { + ATH_MSG_DEBUG( CaloRec::Helpers::fmtMsg("| %02i | %11.11s |",fsamp->first,fsamp->second.c_str()) ); + } + ATH_MSG_DEBUG( CaloRec::Helpers::fmtMsg("+------------+--------------+") ); + } + + return StatusCode::SUCCESS; +} + +StatusCode CaloTopoClusterFromTowerMonitor::book() +{ + //////////////////////////////////////////////////////////// + // composition, multiplicity, kinematics (event by event) // + //////////////////////////////////////////////////////////// + + int detaBins(105); int dphiBins(128); + double detaMin(-1.05); double dphiMin(-Gaudi::Units::pi/2.); + double detaMax(1.05); double dphiMax(Gaudi::Units::pi/2.); + + // re-center delta phi + double dphi((dphiMax-dphiMin)/(1.*dphiBins)); double dphiMinOld(dphiMin); double dphiMaxOld(dphiMax); //int dphiBinsOld(dphiBins); + int dphiBinsOld = (dphiMax-dphiMin)/dphi; + dphiMin -= (dphi/2.); dphiMax += (dphi/2.); + double dphim(dphi); + dphiBins = (dphiMax-dphiMin)/dphim; + + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("re-center delta_phi distributions, old/new binning [%6.3f,%6.3f]/[%6.3f,%6.3f] with %3i/%3i bins %5.3f/%5.3f rad wide", + dphiMinOld,dphiMaxOld,dphiMin,dphiMax,dphiBinsOld,dphiBins,dphi,dphim) ); + + + h_n = bookAny("TowerMultiplicity", "N_{tower}", m_nBins, m_nMin, m_nMax); + h_pt = bookAny("TowerPt", "p_{T}^{tower} [GeV]",m_ptBins,m_ptMin,m_ptMax); + h_nc = bookAny("TowerCellMultiplicity","N_{cell}^{tower}", m_ncBins,m_ncMin,m_ncMax); + h_eta = bookForEta("TowerRapidity", "y_{tower}"); + h_phi = bookForPhi("TowerAzimuth", "#phi_{tower} [rad]"); + + h_samp = bookForSamplings("TowerCellSamplings","TowerCellSamplings"); setAxisTitle(h_samp,"N_{cells} in sampling","y"); + // relations (event by event) + d_n_eta_phi = bookAny("TowerMultiplicityVsEtaVsPhi", "N_{tower}(y_{tower},#phi_{tower})", "y_{tower}",m_etaBins,m_etaMin,m_etaMax,m_phiBins,m_phiMin,m_phiMax); + setAxisTitle(d_n_eta_phi,"#phi_{tower} [rad]","y"); setAxisTitle(d_n_eta_phi,"N_{tower}","z"); + d_nc_eta_phi = bookAny("TowerCellMultiplicityVsEtaBsPhi","N_{cell}^{tower}(y_{tower},#phi_{tower})","y_{tower}",m_etaBins,m_etaMin,m_etaMax,m_phiBins,m_phiMin,m_phiMax); + setAxisTitle(d_nc_eta_phi,"#phi_{tower} [rad]","y"); setAxisTitle(d_nc_eta_phi,"N_{cell}^{tower}","z"); + d_pt_eta = bookAny("TowerPtVsEta", "p_{T}^{tower}(y_{tower})", "y_{tower}",m_etaBins,m_etaMin,m_etaMax,m_ptBins, m_ptMin, m_ptMax ); + setAxisTitle(d_pt_eta,"p_{T}^{tower} [GeV]","y"); + d_nc_eta = bookAny("TowerCellMultiplicityVsEta", "N_{cell}^{tower}(y_{tower})", "y_{tower}",m_etaBins,m_etaMin,m_etaMax,m_ncBins, m_ncMin, m_ncMax ); + setAxisTitle(d_pt_eta,"N_{cell}^{tower} [GeV]","y"); + // distance to nominal + d_deta_eta = bookAny("TowerDeltaEtaVsEta", "(y_{tower}-y_{tower,0})(y_{tower})", "y_{tower}", m_etaBins,m_etaMin,m_etaMax,detaBins,detaMin,detaMax); + setAxisTitle(d_deta_eta,"#Deltay_{tower}","y"); + d_dphi_eta = bookAny("TowerDeltaPhiVsEta", "(#phi_{tower}-#phi_{tower,0})(y_{tower})","y_{tower}", m_etaBins,m_etaMin,m_etaMax,dphiBins,dphiMin,dphiMax); + setAxisTitle(d_dphi_eta,"#Delta#phi_{tower} [rad]","y"); + d_dphi_deta = bookAny("TowerDeltaPhiVsDeltaEta", "#Delta#phi_{tower}(#Deltay_{tower})", "#Deltay_{tower}",detaBins, detaMin, detaMax, dphiBins,dphiMin,dphiMax); + setAxisTitle(d_dphi_deta,"#Delta#phi_{tower} [rad]","y"); + // cell distance to nominal + d_detac_eta = bookAny("CellDeltaEtaVsEta", "(y_{cell #in tower}#minusy_{tower,0})(y_{tower})", "y_{tower}", + m_etaBins,m_etaMin,m_etaMax,detaBins,detaMin,detaMax); + setAxisTitle(d_detac_eta,"#Deltay_{cell #in tower}","y"); + d_dphic_eta = bookAny("CellDeltaPhiVsEta", "(#phi_{cell #in tower}^{cell}#minus#phi_{tower,0})(y_{tower})","y_{tower}", + m_etaBins,m_etaMin,m_etaMax,dphiBins,dphiMin,dphiMax); + setAxisTitle(d_dphic_eta,"#Delta#phi_{cell #in tower}","y"); + d_dphic_detac = bookAny("CellDeltaPhiVsDeltaEta", "#Delta#phi_{cell #in tower}^{cell}(#Deltay_{tower}^{cell})", "#Deltay_{cell #in tower}", + detaBins,detaMin,detaMax,dphiBins,dphiMin,dphiMax); + setAxisTitle(d_dphic_detac,"#Delta#phi_{cell #in tower}","y"); + + /////////////////////////// + // features in samplings // + /////////////////////////// + + int fwbins(250); + double fwmin(-0.1); + double fwmax(4.9); + + int gwbins(60); + double gwmin(-0.1); + double gwmax(1.1); + + int ctbins(20); + double ctmin(-0.5); + double ctmax(19.5); + + // per event + d_wgt_samp = bookForSamplings("EvtFullWeightvsSampling", "[Event] Full weight vs sampling", fwbins,fwmin,fwmax); + setAxisTitle(d_wgt_samp,"w_{geo#plusLCW}","y"); + d_ntt_samp = bookForSamplings("EvtNTowersPerCellvsSampling","[Event] Number of towers/cell vs sampling", ctbins,ctmin,ctmax); + setAxisTitle(d_ntt_samp,"N_{towers}/cell","y"); + d_geo_samp = bookForSamplings("EvtGeoWeightvsSampling", "[Event] Geometrical weight vs sampling", gwbins,gwmin,gwmax); + setAxisTitle(d_geo_samp,"w_{geo}","y"); + d_detac_samp = bookForSamplings("EvtDetaCellvsSampling", "[Event] y_{cell #in tower}#minusy_{0}^{tower}", detaBins,detaMin,detaMax); + setAxisTitle(d_detac_samp,"#Deltay_{cell #in tower}","y"); + d_dphic_samp = bookForSamplings("EvtDphiCellvsSampling", "[Event] #phi_{cell #in tower}#minus#phi_{0}^{tower}",dphiBins,dphiMin,dphiMax); + setAxisTitle(d_dphic_samp,"#Delta#phi_{cell #in tower}","y"); + + // from database + d_maxtowers_samp = bookForSamplings("DBNTowersPerCellvsSampling", "[DB] Number of towers/cell vs sampling",ctbins,ctmin,ctmax); + setAxisTitle(d_maxtowers_samp,"N_{towers}/cell (maximum sharing)","y"); + d_wgttowers_samp = bookForSamplings("DBGeoWeightvsSampling", "[DB] Geometrical weight vs sampling", gwbins,gwmin,gwmax); + setAxisTitle(d_wgttowers_samp,"w_{geo} (from database)","y"); + + ///////////////////// + // features in eta // + ///////////////////// + + int tcbins(100); + double tcmin(-0.5); + double tcmax(99.5); + + d_maxcells_eta = bookForEta("DBNTowerCellsvsEta","[DB] Cells in towers vs y_{tower}", tcbins,tcmin,tcmax); setAxisTitle(d_maxcells_eta,"N_{cells #in tower} (from database)","y"); + d_allwghts_eta = bookForEta("DBGeoWeightvsEta", "[DB] Geometrical weight vs y_{tower}",gwbins,gwmin,gwmax); setAxisTitle(d_allwghts_eta,"w_{geo} (from database)","y"); + + // analyzing the database (in tower eta slices) + double deta((m_etaMax-m_etaMin)/(1.*m_etaBins)); + double eta0(m_etaMin); + + d_maxcells_phi_eta_slice.resize(m_etaBins,(TH2D*)0); + d_allwghts_phi_eta_slice.resize(m_etaBins,(TH2D*)0); + + std::string hname; std::string htitle; + for ( int i(0); i(hname,htitle,tcbins,tcmin,tcmax); setAxisTitle(d_maxcells_phi_eta_slice.at(i),"N_{cells #in tower} (from database)","y"); + hname = CaloRec::Helpers::fmtMsg("DBGeoWeightvsPhiEta%02i",i); + htitle = CaloRec::Helpers::fmtMsg("[DB] Geometrical weight vs #phi_{tower} for y #in [%6.3f,%6.3f]",eta0,eta0+deta); + d_allwghts_phi_eta_slice[i] = bookForPhi(hname,htitle,gwbins,gwmin,gwmax); setAxisTitle(d_allwghts_phi_eta_slice.at(i),"w_{geo} (from database)","y"); + } + + ///////////////////////// + // debugging y=0,phi=0 // + ///////////////////////// + + if ( m_doHotspot ) { + h_n_hs = bookAny("HSTowerMultiplicity", "[Hotspot] N_{tower}", "N_{tower}", m_nBins, m_nMin, m_nMax); + h_pt_hs = bookAny("HSTowerPt", "[Hotspot] p_{T}^{tower}", "p_{T}^{tower} [GeV]",m_ptBins,m_ptMin,m_ptMax); + h_nc_hs = bookAny("HSTowerCellMultiplicity","[Hotspot] N_{cell}^{tower}","N_{cell}^{tower}", m_ncBins,m_ncMin,m_ncMax); + h_eta_hs = bookForEta("HSTowerEta", "[Hotspot] y_{tower}"); + h_phi_hs = bookForPhi("HSTowerPhi", "[Hotspot] #phi_{tower}"); + + h_samp_hs = bookForSamplings("HSTowerCellSampling","HSTowerCellSampling"); + + d_n_eta_phi_hs = bookAny("HSTowerMultiplicityVsEtaVsPhi", "[Hotspot] N_{tower}(y_{tower},#phi_{tower})" "y_{tower}",m_etaBins,m_etaMin,m_etaMax,m_phiBins,m_phiMin,m_phiMax); + setAxisTitle(d_n_eta_phi_hs,"#phi_{tower} [rad]","y"); setAxisTitle(d_n_eta_phi_hs,"N_{tower}","z"); + d_nc_eta_phi_hs = bookAny("HSTowerCellMultiplicityVsEtaVsPhi","[Hotspot] N_{cell}^{tower}(y_{tower},#phi_{tower})","y_{tower}",m_etaBins,m_etaMin,m_etaMax,m_phiBins,m_phiMin,m_phiMax); + setAxisTitle(d_nc_eta_phi_hs,"#phi_{tower} [rad]","y"); setAxisTitle(d_nc_eta_phi_hs,"N_{cell}^{tower}","z"); + d_deta_eta_hs = bookForEta("HSTowerDeltaEtaVsEta", "[Hotspot] #Deltay(y_{tower})",detaBins,detaMin,detaMax); setAxisTitle(d_deta_eta_hs,"#Deltay_{tower}","y"); + setAxisTitle(d_deta_eta_hs,"#Deltay_{tower}","y"); + d_dphi_eta_hs = bookForEta("HSTowerDeltaPhiVsEta", "[Hotspot] #Delta#phi(y_{tower})",dphiBins,dphiMin,dphiMax); + setAxisTitle(d_dphi_eta_hs,"#Delta#phi_{tower} [rad]","y"); + d_dphi_deta_hs = bookAny("HSTowerDeltaPhiVsDeltaEta", "[Hotspot] #Delta#phi_{tower}(#Deltay_{tower})","#Deltay_{tower}",detaBins,detaMin,detaMax,dphiBins,dphiMin,dphiMax); + setAxisTitle(d_dphi_deta_hs,"#Delta#phi_{tower} [rad]","y"); + d_detac_eta_hs = bookForEta("HSCellDeltaEtaVsEta", "[Hotspot] #Deltay_{cell #in tower}(y_{tower})",detaBins,detaMin,detaMax); + setAxisTitle(d_detac_eta_hs,"#Deltay_{cell #in tower}","y"); + d_dphic_eta_hs = bookForEta("HSCellDeltaPhiVsEta", "[Hotspot] #Delta#phi_{cell #in tower}(y_{tower})",dphiBins,dphiMin,dphiMax); + setAxisTitle(d_dphic_eta_hs,"#Delta#phi_{cell #in tower}","y"); + d_dphic_detac_hs = bookAny("HSCellDeltaPhiVsCellDeltaEta", "[Hotspot] #Delta#phi_{cell #in tower}(#Deltay_{cell #in tower})","#Deltay_{cell #in tower}", + detaBins,detaMin,detaMax,dphiBins,dphiMin,dphiMax); + setAxisTitle(d_dphic_detac_hs,"#Delta#phi_{cell #in tower}","y"); + d_detac_samp_hs = bookForSamplings("HSCellDeltaEtavsSampling", "[Hotspot] #Deltay_{cell #in tower} in samplings", detaBins,detaMin,detaMax); + setAxisTitle(d_detac_samp_hs,"#Deltay_{cell #in tower}","y"); + d_dphic_samp_hs = bookForSamplings("HSCelLDeltaPhivsSampling", "[Hotspot] #Delta#phi_{cell #in tower} in samplings",dphiBins,dphiMin,dphiMax); + setAxisTitle(d_dphic_samp_hs,"#Delta#phi_{cell #in tower}","y"); + } + + return StatusCode::SUCCESS; +} + +StatusCode CaloTopoClusterFromTowerMonitor::execute() +{ + SG::ReadHandle towerHandle(m_towerContainerKey); + if ( !towerHandle.isValid() ) { + ATH_MSG_ERROR( "cannot allocate cluster container with key <" << m_towerContainerKey << ">" ); + return StatusCode::FAILURE; + } + + // fill plots + size_t nhs(0); + m_cellTags.reset(); + std::vector deta; deta.reserve(1000); + std::vector dphi; dphi.reserve(deta.capacity()); + std::vector csam; csam.reserve(deta.capacity()); + h_n->Fill(1.0*towerHandle->size()); + for ( auto ftow(towerHandle->begin()); ftow!=towerHandle->end(); ++ftow ) { + const xAOD::CaloCluster* ptow = *ftow; + // tower kinematics + double pt(ptow->pt()/Gaudi::Units::GeV); + double eta(ptow->eta()); + double phi(CaloPhiRange::fix(ptow->phi())); + double nc(1.0*ptow->size()); + // tower composition and get distances + this->fillComposition(*ptow,deta,dphi,csam); + double deltaEta(ptow->eta()-ptow->eta0()); + double deltaPhi(CaloPhiRange::fix(ptow->phi()-ptow->phi0())); + // inclusive plots + h_pt->Fill(pt); h_eta->Fill(eta); h_phi->Fill(phi); h_nc->Fill(nc); + d_n_eta_phi->Fill(eta,phi); d_nc_eta_phi->Fill(eta,phi,nc); + d_pt_eta->Fill(eta,pt); + d_nc_eta->Fill(eta,nc); + d_deta_eta->Fill(eta,deltaEta); + d_dphi_eta->Fill(eta,deltaPhi); + d_dphi_deta->Fill(deltaEta,deltaPhi); + for ( size_t ic(0); icFill(eta,deta.at(ic)); + d_dphic_eta->Fill(eta,dphi.at(ic)); + d_dphic_detac->Fill(deta.at(ic),dphi.at(ic)); + fillSampling(d_detac_samp,csam.at(ic),deta.at(ic)); + fillSampling(d_dphic_samp,csam.at(ic),dphi.at(ic)); + fillSampling(h_samp,csam.at(ic)); + } + // hot spot + if ( this->isInHotspot(*ptow) ) { + ++nhs; + h_pt_hs->Fill(pt); + h_eta_hs->Fill(eta); + h_phi_hs->Fill(phi); + h_nc_hs->Fill(nc); + d_n_eta_phi_hs->Fill(eta,phi); + d_nc_eta_phi_hs->Fill(eta,phi,nc); + d_deta_eta_hs->Fill(eta,deltaEta); + d_dphi_eta_hs->Fill(eta,deltaPhi); + d_dphi_deta_hs->Fill(deltaEta,deltaPhi); + for ( size_t ic(0); icFill(eta,deta.at(ic)); + d_dphic_eta_hs->Fill(eta,dphi.at(ic)); + d_dphic_detac_hs->Fill(deta.at(ic),dphi.at(ic)); + fillSampling(d_detac_samp_hs,csam.at(ic),deta.at(ic)); + fillSampling(d_dphic_samp_hs,csam.at(ic),dphi.at(ic)); + fillSampling(h_samp_hs,csam.at(ic)); + } + } + // fill in samplings + for ( auto fCell(ptow->getCellLinks()->begin()); fCell!=ptow->getCellLinks()->end(); ++fCell ) { + CaloSampling::CaloSample csamp = (*fCell)->caloDDE()->getSampling(); + fillSampling(d_wgt_samp,csamp,fCell.weight()); + // check for cell which other tower it contributes to and plot geometrical weights (only 1/cell) + size_t chash(static_cast((*fCell)->caloDDE()->calo_hash())); + if ( !m_cellTags.test(chash) ) { + CaloTowerGeometrySvc::elementvector_t lOfTowers = m_towerGeometrySvc->getTowers(chash); + double ntowers(1.0*lOfTowers.size()); + fillSampling(d_ntt_samp,csamp,ntowers); + for ( auto elm : lOfTowers ) { fillSampling(d_geo_samp,csamp,m_towerGeometrySvc->cellWeight(elm)); } + m_cellTags.set(chash); + } + } + } + // number of towers in hot spot + h_n_hs->Fill(1.0*nhs); + + return StatusCode::SUCCESS; +} + +bool CaloTopoClusterFromTowerMonitor::fillComposition(const xAOD::CaloCluster& ptow,std::vector& deta,std::vector& dphi,std::vector& csam) const +{ + deta.clear(); dphi.clear(); csam.clear(); + for ( auto fCell(ptow.getCellLinks()->begin()); fCell != ptow.getCellLinks()->end(); ++fCell ) { + deta.push_back((*fCell)->eta()-ptow.eta0()); + dphi.push_back(CaloPhiRange::fix((*fCell)->phi()-ptow.phi0())); + csam.push_back((*fCell)->caloDDE()->getSampling()); + } + return !deta.empty(); +} diff --git a/Calorimeter/CaloRec/src/CaloTopoClusterTowerMerger.cxx b/Calorimeter/CaloRec/src/CaloTopoClusterTowerMerger.cxx new file mode 100644 index 0000000000000000000000000000000000000000..19e7b6085afa6613fa2c444f19d95b6c718beba4 --- /dev/null +++ b/Calorimeter/CaloRec/src/CaloTopoClusterTowerMerger.cxx @@ -0,0 +1,214 @@ +/* + Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration +*/ + +#include "StoreGate/ReadHandle.h" +#include "StoreGate/WriteHandle.h" + +#include "xAODCaloEvent/CaloClusterAuxContainer.h" + +#include "CaloRec/CaloTopoClusterTowerMerger.h" +#include "CaloRec/CaloTopoClusterFromTowerHelpers.h" + +#include "CaloUtils/CaloClusterStoreHelper.h" + +// #include + +// std::atomic CaloTopoClusterTowerMerger_checkMoments(true); + +// #define CL_MNAME( NAME ) xAOD::CaloCluster::NAME +// #define CL_ENTRY( NAME ) std::tuple< CL_MNAME( MomentType ),std::string>( CL_MNAME( NAME ), #NAME ) + +#define CL_RHMSG( NAME ) MsgStream& operator<<(MsgStream& mstr,const SG::ReadHandleKey< NAME >& ckey ) { mstr << ckey.key(); return mstr; } +#define CL_WHMSG( NAME ) MsgStream& operator<<(MsgStream& mstr,const SG::WriteHandleKey< NAME >& ckey ) { mstr << ckey.key(); return mstr; } + +namespace { + CL_RHMSG( xAOD::CaloClusterContainer ) +} + +// std::vector CaloTopoClusterTowerMerger::m_momentList = std::vector< CL_MNAME( MomentType ) >(); +// std::vector > CaloTopoClusterTowerMerger::m_momentMap { +// CL_ENTRY( FIRST_PHI ), +// CL_ENTRY( FIRST_ETA ), +// CL_ENTRY( SECOND_R ), +// CL_ENTRY( SECOND_LAMBDA ), +// CL_ENTRY( DELTA_PHI ), +// CL_ENTRY( DELTA_THETA ), +// CL_ENTRY( DELTA_ALPHA ), +// CL_ENTRY( CENTER_X ), +// CL_ENTRY( CENTER_Y ), +// CL_ENTRY( CENTER_Z ), +// CL_ENTRY( CENTER_MAG ), +// CL_ENTRY( CENTER_LAMBDA ), +// CL_ENTRY( LATERAL ), +// CL_ENTRY( LONGITUDINAL ), +// CL_ENTRY( ENG_FRAC_EM ), +// CL_ENTRY( ENG_FRAC_MAX ), +// CL_ENTRY( ENG_FRAC_CORE ), +// CL_ENTRY( FIRST_ENG_DENS ), +// CL_ENTRY( SECOND_ENG_DENS ), +// CL_ENTRY( ISOLATION ), +// CL_ENTRY( ENG_BAD_CELLS ), +// CL_ENTRY( N_BAD_CELLS ), +// CL_ENTRY( N_BAD_CELLS_CORR ), +// CL_ENTRY( BAD_CELLS_CORR_E ), +// CL_ENTRY( BADLARQ_FRAC ), +// CL_ENTRY( ENG_POS ), +// CL_ENTRY( SIGNIFICANCE ), +// CL_ENTRY( CELL_SIGNIFICANCE ), +// CL_ENTRY( CELL_SIG_SAMPLING ), +// CL_ENTRY( AVG_LAR_Q ), +// CL_ENTRY( AVG_TILE_Q ), +// CL_ENTRY( ENG_BAD_HV_CELLS ), +// CL_ENTRY( N_BAD_HV_CELLS ), +// CL_ENTRY( PTD ), +// CL_ENTRY( EM_PROBABILITY ), +// CL_ENTRY( HAD_WEIGHT ), +// CL_ENTRY( OOC_WEIGHT ), +// CL_ENTRY( DM_WEIGHT ), +// CL_ENTRY( TILE_CONFIDENCE_LEVEL ), +// CL_ENTRY( VERTEX_FRACTION ), +// CL_ENTRY( NVERTEX_FRACTION ), +// CL_ENTRY( ETACALOFRAME ), +// CL_ENTRY( PHICALOFRAME ), +// CL_ENTRY( ETA1CALOFRAME ), +// CL_ENTRY( PHI1CALOFRAME ), +// CL_ENTRY( ETA2CALOFRAME ), +// CL_ENTRY( PHI2CALOFRAME ), +// CL_ENTRY( ENG_CALIB_TOT ), +// CL_ENTRY( ENG_CALIB_OUT_L ), +// CL_ENTRY( ENG_CALIB_OUT_M ), +// CL_ENTRY( ENG_CALIB_OUT_T ), +// CL_ENTRY( ENG_CALIB_DEAD_L ), +// CL_ENTRY( ENG_CALIB_DEAD_M ), +// CL_ENTRY( ENG_CALIB_DEAD_T ), +// CL_ENTRY( ENG_CALIB_EMB0 ), +// CL_ENTRY( ENG_CALIB_EME0 ), +// CL_ENTRY( ENG_CALIB_TILEG3 ), +// CL_ENTRY( ENG_CALIB_DEAD_TOT ), +// CL_ENTRY( ENG_CALIB_DEAD_EMB0 ), +// CL_ENTRY( ENG_CALIB_DEAD_TILE0 ), +// CL_ENTRY( ENG_CALIB_DEAD_TILEG3 ), +// CL_ENTRY( ENG_CALIB_DEAD_EME0 ), +// CL_ENTRY( ENG_CALIB_DEAD_HEC0 ), +// CL_ENTRY( ENG_CALIB_DEAD_FCAL ), +// // CL_ENTRY( ENG_CALIB_DEAD_LEAKAG ), // not in r21 +// // CL_ENTRY( ENG_CALIB_DEAD_UNCLAS ), // not in r21 +// CL_ENTRY( ENG_CALIB_FRAC_EM ), +// CL_ENTRY( ENG_CALIB_FRAC_HAD ), +// CL_ENTRY( ENG_CALIB_FRAC_REST ) +// }; + +CaloTopoClusterTowerMerger::CaloTopoClusterTowerMerger(const std::string& name,ISvcLocator* pSvcLocator) + : AthAlgorithm(name,pSvcLocator) + , m_clusterContainerKey("CaloCalTopoCluster") + , m_towerContainerKey("CaloCalFwdTopoTower") + , m_topoSignalContainerKey("CaloCalTopoSignal") + // ** not in r21 **, m_cellLinkContainerKey("") + , m_clusterRange(3.2) +{ + declareProperty("TopoClusterContainerKey",m_clusterContainerKey, "Topo-cluster container key" ); + declareProperty("TopoTowerContainerKey", m_towerContainerKey, "Topo-tower container key" ); + declareProperty("TopoSignalContainerKey", m_topoSignalContainerKey,"Topo-signal container key" ); + declareProperty("TopoClusterRange", m_clusterRange, "Rapidity range for using topo-clusters in combined signal mode"); +} + +CaloTopoClusterTowerMerger::~CaloTopoClusterTowerMerger() +{ } + +StatusCode CaloTopoClusterTowerMerger::initialize() +{ + if ( m_clusterRange <= 0. ) { + ATH_MSG_ERROR( CaloRec::Helpers::fmtMsg("Invalid topo-cluster range |y| < %6.3f - algorithm non-functional",m_clusterRange) ); + return StatusCode::FAILURE; + } + + // ** not in r21 ** if ( m_cellLinkContainerKey.key().empty() ) { m_cellLinkKey = m_topoSignalContainerKey.key() + std::string("_links"); } + + ATH_CHECK( m_clusterContainerKey.initialize() ); + ATH_CHECK( m_towerContainerKey.initialize() ); + ATH_CHECK( m_topoSignalContainerKey.initialize() ); + // ** not in r21 ** ATH_CHECK( m_cellLinkContainerKey.initialize() ); + + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("Topo_cluster with |y| < %.2f will be merged with topo-towers with |y| > %.2f",m_clusterRange,m_clusterRange) ); + // if ( m_copyMoments ) { ATH_MSG_INFO( "Cluster moments are explicitly copied." ); } + + return StatusCode::SUCCESS; +} + +StatusCode CaloTopoClusterTowerMerger::execute() +{ + + // collect input + rhandle_t clusterHandle(m_clusterContainerKey); + if ( !clusterHandle.isValid() ) { + ATH_MSG_WARNING( "Topo-cluster container with key <" << m_clusterContainerKey << "> not found" ); + return StatusCode::SUCCESS; + } + rhandle_t towerHandle(m_towerContainerKey); + if ( !towerHandle.isValid() ) { + ATH_MSG_WARNING( "Topo-tower container with key <" << m_towerContainerKey << "> not found" ); + return StatusCode::SUCCESS; + } + + // prepare output + whandle_t signalHandle(m_topoSignalContainerKey); + ATH_CHECK(this->addContainerWriteHandle(signalHandle)); + + // check cluster moments + // if ( m_copyMoments ) { if ( CaloTopoClusterTowerMerger_checkMoments ) { this->fillMoments(*(towerHandle.ptr()->front())); CaloTopoClusterTowerMerger_checkMoments = false; } } + + // fill output from topo-clusters + for ( auto pClus : *clusterHandle ) { if ( clusterFilter(*pClus) ) { this->makeDeepCopy(*pClus,signalHandle.ptr()); } } + // fill output from topo-towers + for ( auto pTowr : *towerHandle ) { if ( towerFilter(*pTowr) ) { this->makeDeepCopy(*pTowr,signalHandle.ptr()); } } + + // take care of cell links + ATH_CHECK(CaloClusterStoreHelper::finalizeClusters(&(*evtStore()),signalHandle.ptr(),m_topoSignalContainerKey.key(),msg())); + + return StatusCode::SUCCESS; +} + +// bool CaloTopoClusterTowerMerger::fillMoments(const xAOD::CaloCluster& rClus) +// { +// m_momentList.reserve(m_momentMap.size()); +// double val(0); +// if ( !msgLvl(MSG::DEBUG) ) { +// for ( auto fmom(m_momentMap.begin()); fmom != m_momentMap.end(); ++fmom ) { +// if ( rClus.retrieveMoment(std::get<0>(*fmom),val) ) { m_momentList.push_back(std::get<0>(*fmom)); } +// } +// } else { +// ATH_MSG_DEBUG("###############################"); +// ATH_MSG_DEBUG("## Moments (topo-towers) ##"); +// ATH_MSG_DEBUG("##---------------------------##"); +// for ( auto fmom(m_momentMap.begin()); fmom != m_momentMap.end(); ++fmom ) { +// if ( rClus.retrieveMoment( std::get<0>(*fmom), val ) ) { +// m_momentList.push_back(std::get<0>(*fmom)); +// ATH_MSG_DEBUG( CaloRec::Helpers::fmtMsg("## %5.5x | \042%16.16s\042 ##", std::get<0>(*fmom), std::get<1>(*fmom).c_str()) ); +// } +// } +// ATH_MSG_DEBUG("###############################"); +// } // debug mode +// return true; +// } + +bool CaloTopoClusterTowerMerger::makeDeepCopy(const xAOD::CaloCluster& rClus,xAOD::CaloClusterContainer* pClusCont) const +{ pClusCont->push_back(new xAOD::CaloCluster(rClus)); return true; } + +StatusCode CaloTopoClusterTowerMerger::addContainerWriteHandle(whandle_t& signalHandle) +{ + // get a new signal handle + signalHandle = std::unique_ptr(new xAOD::CaloClusterContainer()); + if ( !signalHandle.isValid() ) { return StatusCode::FAILURE; } + // get AUX container + xAOD::CaloClusterAuxContainer* auxData = new xAOD::CaloClusterAuxContainer(); + std::string auxName(m_topoSignalContainerKey.key()+"Aux."); + if ( evtStore()->overwrite(auxData,auxName).isFailure() ) { + ATH_MSG_ERROR("Failed to record xAOD::CaloClusterAuxContainer with key <" << auxName << ">"); + delete auxData; + return StatusCode::FAILURE; + } + // connect store with object container + signalHandle.ptr()->setStore(auxData); + return StatusCode::SUCCESS; +} diff --git a/Calorimeter/CaloRec/src/CaloTowerGeometrySvc.cxx b/Calorimeter/CaloRec/src/CaloTowerGeometrySvc.cxx new file mode 100644 index 0000000000000000000000000000000000000000..e1dfdcabc7b382581853a7a7dc9fd908d4864310 --- /dev/null +++ b/Calorimeter/CaloRec/src/CaloTowerGeometrySvc.cxx @@ -0,0 +1,410 @@ +/* Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration */ +#include "CaloRec/CaloTowerGeometrySvc.h" +#include "CaloRec/CaloTopoClusterFromTowerHelpers.h" + +#include + +namespace { constexpr auto _pi = 3.14159265358979323846; } + +CaloTowerGeometrySvc::index_t CaloTowerGeometrySvc::m_invalidIndex = index_t(-1); +double CaloTowerGeometrySvc::m_invalidValue = -999.; + +CaloTowerGeometrySvc::CaloTowerGeometrySvc(const std::string& name,ISvcLocator* pSvcLocator) + : AthService(name,pSvcLocator) + , m_caloDDM((const CaloDetDescrManager*)0) + , m_logFileName("logfile") + , m_towerEtaWidth(0.) + , m_towerPhiWidth(0.) + , m_towerArea(0.) + , m_towerBins(0) + , m_maxCellHash(0) //----------------------------------------// + , m_towerEtaBins(100) // Default tower definition is "hadronic" // + , m_towerEtaMin(-5.0) // towers of size 0.1 x pi/32. // + , m_towerEtaMax(5.0) //----------------------------------------// + , m_adjustEta(true) + , m_towerPhiBins(64) + , m_towerPhiMin(-_pi) //----------------------------------------// + , m_towerPhiMax(_pi) // FCal vertical and horizontal cell // + , m_fcal1Xslice(8.) // slicing creates "mini-cells" which are // + , m_fcal1Yslice(8.) // then projected onto towers. The mini- // + , m_fcal2Xslice(8.) // cell signal is 1/(Nx x Ny) x Ecell, // + , m_fcal2Yslice(12.) // where Nx(y) are the number of x(y) // + , m_fcal3Xslice(12.) // slices. // + , m_fcal3Yslice(12.) //----------------------------------------// +{ + // Properties + declareProperty("TowerEtaBins", m_towerEtaBins, "Number of pseudorapidity bins in tower grid"); + declareProperty("TowerEtaMin", m_towerEtaMin, "Lower boundary of pseudorapidity range"); + declareProperty("TowerEtaMax", m_towerEtaMax, "Upper boundary of pseudorapidity range"); + declareProperty("AdjustFCalToTowerEtaBounds", m_adjustEta, "Adjust FCal cells to eta boundaries"); + declareProperty("TowerPhiBins", m_towerPhiBins, "Number of azimuthal bins in tower grid"); + declareProperty("TowerPhiMin", m_towerPhiMin, "Lower boundary of azimuthal range"); + declareProperty("TowerPhiMax", m_towerPhiMax, "Upper boundary of azimuthal range"); + // change only for R&D + declareProperty("FCal1NSlicesX", m_fcal1Xslice, "Number of X slices for FCal1 cells"); + declareProperty("FCal1NSlicesY", m_fcal1Yslice, "Number of Y slices for FCal1 cells"); + declareProperty("FCal2NSlicesX", m_fcal2Xslice, "Number of X slices for FCal2 cells"); + declareProperty("FCal2NSlicesY", m_fcal2Yslice, "Number of Y slices for FCal2 cells"); + declareProperty("FCal3NSlicesX", m_fcal3Xslice, "Number of X slices for FCal3 cells"); + declareProperty("FCal3NSlicesY", m_fcal3Yslice, "Number of Y slices for FCal3 cells"); +} + +//-------------// +// Gaudi stuff // +//-------------// + +StatusCode CaloTowerGeometrySvc::queryInterface(const InterfaceID& riid,void** ppvInterface) +{ + if ( CaloTowerGeometrySvc::interfaceID().versionMatch(riid) ) { + *ppvInterface = this; + return StatusCode::SUCCESS; + } else { + return AthService::queryInterface(riid,ppvInterface); + } +} + +//-------------------------// +// Initialize and Finalize // +//-------------------------// + +StatusCode CaloTowerGeometrySvc::initialize() +{ + // set up services and managers + if ( f_setupSvc().isFailure() ) { + ATH_MSG_ERROR("Cannot allocate the calorimeter detector description manager"); + return StatusCode::FAILURE; + } + + // prepare FCal segmentation + m_ndxFCal[0] = m_fcal1Xslice; m_ndxFCal[1] = m_fcal2Xslice; m_ndxFCal[2] = m_fcal3Xslice; + m_ndyFCal[0] = m_fcal1Yslice; m_ndyFCal[1] = m_fcal2Yslice; m_ndyFCal[2] = m_fcal3Yslice; + for ( uint_t i(0); i 0 ) { + m_towerEtaWidth = (m_towerEtaMax-m_towerEtaMin)/((double)m_towerEtaBins); + } else { + ATH_MSG_ERROR("Number of tower eta bins is invalid (" << m_towerEtaBins << " bins)"); + return StatusCode::FAILURE; + } + if ( m_towerPhiBins > 0 ) { + m_towerPhiWidth = (m_towerPhiMax-m_towerPhiMin)/((double)m_towerPhiBins); + } else { + ATH_MSG_ERROR("Number of tower phi bins is invalid (" << m_towerPhiBins << " bins)"); + return StatusCode::FAILURE; + } + + m_towerBins = m_towerEtaBins*m_towerPhiBins; + m_towerArea = m_towerEtaWidth*m_towerPhiWidth; + m_maxCellHash = f_caloDDM()->element_size()-1; + m_numberOfCells = f_caloDDM()->element_size(); + + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("Tower description - Eta(bins,min,max,width) = (%3zu,%6.3f,%6.3f,%6.4f)",m_towerEtaBins,m_towerEtaMin,m_towerEtaMax,m_towerEtaWidth) ); + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("Tower description - Phi(bins,min,max,width) = (%3zu,%6.3f,%6.3f,%6.4f)",m_towerPhiBins,m_towerPhiMin,m_towerPhiMax,m_towerPhiWidth) ); + ATH_MSG_INFO( CaloRec::Helpers::fmtMsg("Tower description - maximum number of towers %5zu; (eta,phi)-space area %6.4f; maximum cell hash index %6zu",m_towerBins,m_towerArea,m_maxCellHash) ); + + if ( this->msgLvl(MSG::VERBOSE) ) { + std::vector logger; logger.reserve(m_towerBins+10); + logger.push_back(CaloRec::Helpers::fmtMsg("[CaloTowerGeometrySvc::%s] +-------- Tower Index Mapping ---------+------------+------------+",this->name().c_str())); + logger.push_back(CaloRec::Helpers::fmtMsg("[CaloTowerGeometrySvc::%s] | %10.10s | %10.10s | %10.10s | %10.10s | %10.10s |",this->name().c_str(),"TowerIndex", " EtaIndex ", " PhiIndex "," Eta "," Phi ")); + logger.push_back(CaloRec::Helpers::fmtMsg("[CaloTowerGeometrySvc::%s] +------------+------------+------------+------------+------------+",this->name().c_str())); + for ( size_t i(0); ietaIndexFromTowerIndex(i)); + size_t phiIndex(this->phiIndexFromTowerIndex(i)); + double eta(this->towerEta(i)); double phi(this->towerPhi(i)); + logger.push_back(CaloRec::Helpers::fmtMsg("[CaloTowerGeometrySvc::%s] | %5zu | %5zu | %5zu | %6.3f | %6.3f |",this->name().c_str(),i,etaIndex,phiIndex,eta,phi)); + } + logger.push_back(CaloRec::Helpers::fmtMsg("[CaloTowerGeometrySvc::%s] +------------+------------+------------+------------+------------+",this->name().c_str())); + // + std::string logFileName = m_logFileName+"."+this->name()+".dat"; + std::ofstream logfile; + logfile.open(logFileName); for ( auto mlog : logger ) { logfile << mlog << std::endl; } logfile.close(); + } + return f_setupTowerGrid(); +} + +StatusCode CaloTowerGeometrySvc::finalize() +{ return StatusCode::SUCCESS; } + +//-------// +// Setup // +//-------// + +StatusCode CaloTowerGeometrySvc::f_setupTowerGrid() +{ + // initialized + if ( m_maxCellHash == 0 ) { + ATH_MSG_ERROR("Service not initialized? Maximum cell hash is " << m_maxCellHash ); + return StatusCode::FAILURE; + } + + // payload template + elementvector_t ev; + + // set up lookup table + ATH_MSG_INFO( "Setting up cell-to-tower lookup for " << m_numberOfCells << " calorimeter cells" ); + m_towerLookup.resize(m_numberOfCells,ev); + + // dump projections + std::ofstream logger; + if ( msgLvl(MSG::DEBUG) ) { + logger.open(CaloRec::Helpers::fmtMsg("%s.calocellprojections_dump.dat",this->name().c_str())); + logger << CaloRec::Helpers::fmtMsg("############################################################################") << std::endl; + logger << CaloRec::Helpers::fmtMsg("## CaloCell projections onto %3.1f x %3.1f tower grid (projective cells only) ##",m_towerEtaWidth,m_towerPhiWidth) << std::endl; + logger << CaloRec::Helpers::fmtMsg("############################################################################") << std::endl; + logger << std::endl; + } + + // loop cells + for ( auto fcell(m_caloDDM->element_begin()); fcell != m_caloDDM->element_end(); ++fcell ) { + // reference cell descriptor + const CaloDetDescrElement* pCaloDDE = *fcell; + // check hash id validity + index_t cidx(pCaloDDE->calo_hash()); + if ( cidx >= m_towerLookup.size() ) { + ATH_MSG_WARNING( "Cell hash identifier out of range " << cidx << "/" << m_towerLookup.size() << ", ignore cell" ); + } else { + if ( pCaloDDE->is_lar_fcal() ) { + if ( this->f_setupTowerGridFCal(pCaloDDE,logger).isFailure() ) { return StatusCode::FAILURE; } + } else { + if ( this->f_setupTowerGridProj(pCaloDDE,logger).isFailure() ) { return StatusCode::FAILURE; } + } + } // cell hash in range? + } // loop cell descriptors + if ( logger.is_open() ) { logger.close(); } + + return StatusCode::SUCCESS; +} + +StatusCode CaloTowerGeometrySvc::f_setupTowerGridFCal(const CaloDetDescrElement* pCaloDDE,std::ofstream& /*logger*/) +{ + //-----------------------------------------------------------------------------------------// + // FCal special - the rectangular (in linear space) calorimeter cells are sub-divided into // + // small cells and then shared across as many towers as the small cells cover. // + //-----------------------------------------------------------------------------------------// + + // collect geometrical variables + int cLayer(pCaloDDE->getLayer()-1); // FCal layer number 1..3 -> array indices 0..2 + + double cXpos(pCaloDDE->x()); // FCal cell x position (cell center) + double cYpos(pCaloDDE->y()); // FCal cell y position (cell center) + double cZpos(pCaloDDE->z()); // FCal cell z position (cell center) + double cZpos2(cZpos*cZpos); + + double cXwid(pCaloDDE->dx()); // FCal cell x full width + double cYwid(pCaloDDE->dy()); // FCal cell y full width + + // double xSlice(cXwid/m_ndxFCal[cLayer]); // FCal cell x slize width + // double ySlice(cYwid/m_ndyFCal[cLayer]); // FCal cell y slice width + double cWght(m_wgtFCal[cLayer]); // FCal cell geometrical (signal) weight + + int nXslice((int)m_ndxFCal[cLayer]); // FCal number of x slices + int nYslice((int)m_ndyFCal[cLayer]); // FCal number of y slices + double cXstp(cXwid/((double)nXslice)); // FCal slice x width + double cYstp(cYwid/((double)nYslice)); // FCal slice y width + + // fill cell fragments + // double xoff(cXpos-cXwid/2.+cXstp/2.); double yoff(cYpos-cYwid/2.+cYstp/2.); + double x(cXpos-(cXwid-cXstp)/2.); + double xlim(cXpos+cXwid/2.); double ylim(cYpos+cYwid/2.); + double etaOrig(0.); + // for ( int ix(0); ix < nXslice; ++ix ) { + // double x(xoff+ix*cXstp); + while ( x < xlim ) { + // for ( int iy(0); iy < nYslice; ++iy ) { + // double y(yoff+iy*cYstp); + double y(cYpos-(cYwid-cYstp)/2.); + while ( y < ylim ) { + double r(std::sqrt(x*x+y*y+cZpos2)); + double eta(-0.5*std::log((r-cZpos)/(r+cZpos))); + bool etaAdjusted(false); + if ( m_adjustEta ) { + if ( eta < m_towerEtaMin ) { + etaAdjusted = true; + etaOrig = eta; + eta = m_towerEtaMin+m_towerEtaWidth/2.; + } else if ( eta > m_towerEtaMax ) { + etaAdjusted = true; + etaOrig = eta; + eta = m_towerEtaMax-m_towerEtaWidth/2.; + } + } + double phi(CaloPhiRange::fix(std::atan2(y,x))); + index_t towerIdx(this->towerIndex(eta,phi)); + // tower index not valid + if ( isInvalidIndex(towerIdx) ) { + ATH_MSG_WARNING("Found invalid tower index for FCal cell (eta,phi) = (" << eta << "," << phi << ") at (x,y,z) = (" << x << "," << y << "," << cZpos << ") [cell ignored]"); + } else { + // add tower to lookup + if ( etaAdjusted ) { + ATH_MSG_WARNING("FCal cell direction (eta,phi) = (" << etaOrig << "," << phi << ") for cell at (x,y,z) = (" + << x << "," << y << "," << cZpos << ") adjusted to (eta,phi) = (" << eta << "," << phi << ") [cell adjusted]"); + } + f_assign(pCaloDDE->calo_hash(),towerIdx,cWght); + } // tower index ok + y += cYstp; + } // loop on y fragments + x += cXstp; + } // loop on x fragments + return StatusCode::SUCCESS; +} + +StatusCode CaloTowerGeometrySvc::f_setupTowerGridProj(const CaloDetDescrElement* pCaloDDE,std::ofstream& logger) +{ + // projective readout calos - collect geometrical variables + double cEtaPos(pCaloDDE->eta_raw()); // projective cell center in pseudorapidity + double cEtaWid(pCaloDDE->deta()); // projective cell width in pseudorapidity + double cPhiPos(pCaloDDE->phi_raw()); // projective cell center in azimuth + double cPhiWid(pCaloDDE->dphi()); // projective cell width in azimuth + + // check cell-tower overlap area fractions + uint_t kEta(static_cast(cEtaWid/m_towerEtaWidth+0.5)); kEta = kEta == 0 ? 1 : kEta; // fully contained cell may have 0 fragments (protection) + uint_t kPhi(static_cast(cPhiWid/m_towerPhiWidth+0.5)); kPhi = kPhi == 0 ? 1 : kPhi; + + // print out + if ( kEta > 1 || kPhi > 1 ) { + ATH_MSG_VERBOSE("Found cell [" << pCaloDDE->calo_hash() << "/0x" << std::hex << pCaloDDE->identify().get_compact() << std::dec << "] spawning several towers." + << " Neta = " << kEta << ", Nphi = " << kPhi ); + } + + // share cells + double cWght(1./((double)kEta*kPhi)); // area weight + double sEta(cEtaWid/((double)kEta)); // step size (pseudorapidity) + double sPhi(cPhiWid/((double)kPhi)); // step size (azimuth) + double oEta(cEtaPos-sEta/2.); // offset (pseudorapidity) + double oPhi(cPhiPos-sPhi/2.); // offset (azimuth) + + // loop over cell fragments + for ( uint_t ie(1); ie<=kEta; ++ie ) { + double ceta(oEta+((double)ie-0.5)*sEta); // eta of fragment + for ( uint_t ip(1); ip<=kPhi; ++ip ) { + double cphi(oPhi+((double)ip-0.5)*sPhi); // phi fragment + // tower index + index_t towerIdx(this->towerIndex(ceta,cphi)); + if ( isInvalidIndex(towerIdx) ) { + ATH_MSG_ERROR("Found invalid tower index for non-FCal cell (id,eta,phi) = (" << pCaloDDE->calo_hash() << "," << ceta << "," << cphi << ")"); + return StatusCode::FAILURE; + } // invalid tower index + // m_towerLookup[pCaloDDE->calo_hash()].emplace_back(towerIdx,cWght); // add to tower lookup + f_assign(pCaloDDE->calo_hash(),towerIdx,cWght); + if ( logger.is_open() ) { + double el(this->towerEta(towerIdx)-this->etaWidth()/2.); double eh(this->towerEta(towerIdx)+this->etaWidth()/2.); + double pl(this->towerPhi(towerIdx)-this->phiWidth()/2.); double ph(this->towerPhi(towerIdx)+this->phiWidth()/2.); + if ( cEtaPos >= el && cEtaPos <= eh && cPhiPos >= pl && cPhiPos <= ph ) { + logger << CaloRec::Helpers::fmtMsg("Cell [%6zu] at (%6.3f,%6.3f) in Sampling %10.10s contributes to tower [%5zu] at ([%6.3f,%6.3f],[%6.3f,%6.3f]) with weight %6.4f [IN_BOUNDS]", + (size_t)pCaloDDE->calo_hash(),cEtaPos,cPhiPos,CaloRec::Lookup::getSamplingName(pCaloDDE->getSampling()).c_str(),towerIdx,el,eh,pl,ph,cWght) << std::endl; + } else { + logger << CaloRec::Helpers::fmtMsg("Cell [%6zu] at (%6.3f,%6.3f) in Sampling %10.10s contributes to tower [%5zu] at ([%6.3f,%6.3f],[%6.3f,%6.3f]) with weight %6.4f [OUT_OF_BOUNDS]", + (size_t)pCaloDDE->calo_hash(),cEtaPos,cPhiPos,CaloRec::Lookup::getSamplingName(pCaloDDE->getSampling()).c_str(),towerIdx,el,eh,pl,ph,cWght) << std::endl; + } + } // debugging central region + } // phi fragment loop + } // eta fragment loop + return StatusCode::SUCCESS; +} + +//------// +// Fill // +//------// + +double CaloTowerGeometrySvc::f_assign(IdentifierHash cellHash,index_t towerIdx,double wght) +{ + // check if cell-tower already related + uint_t cidx(static_cast(cellHash)); + for ( element_t& elm : m_towerLookup.at(cidx) ) { + if ( towerIndex(elm) == towerIdx ) { std::get<1>(elm) += wght; return cellWeight(elm); } + } + // not yet related + m_towerLookup[cidx].emplace_back(towerIdx,wght); + return cellWeight(m_towerLookup.at(cidx).back()); +} + +//--------// +// Access // +//--------// + +StatusCode CaloTowerGeometrySvc::access(const CaloCell* pCell,std::vector& towerIdx,std::vector& towerWghts) const +{ return this->access(f_caloDDE(pCell)->calo_hash(),towerIdx,towerWghts); } + +StatusCode CaloTowerGeometrySvc::access(IdentifierHash cellHash,std::vector& towerIdx,std::vector& towerWghts) const +{ + towerIdx.clear(); + towerWghts.clear(); + + uint_t cidx(static_cast(cellHash)); + + if ( cidx >= m_towerLookup.size() ) { + ATH_MSG_WARNING("Invalid cell hash index " << cellHash << ", corresponding index " << cidx << " not found in tower lookup"); + return StatusCode::SUCCESS; + } + + if ( towerIdx.capacity() < m_towerLookup.at(cidx).size() ) { towerIdx.reserve(m_towerLookup.at(cidx).size()); } + if ( towerWghts.capacity() < m_towerLookup.at(cidx).size() ) { towerWghts.reserve(m_towerLookup.at(cidx).size()); } + + for ( const auto& elm : m_towerLookup.at(cidx) ) { towerIdx.push_back(towerIndex(elm)); towerWghts.push_back(cellWeight(elm)); } + + return StatusCode::SUCCESS; +} + +CaloTowerGeometrySvc::elementvector_t CaloTowerGeometrySvc::getTowers(const CaloCell* pCell) const +{ return pCell != 0 ? this->getTowers(f_caloDDE(pCell)->calo_hash()) : elementvector_t(); } + +CaloTowerGeometrySvc::elementvector_t CaloTowerGeometrySvc::getTowers(IdentifierHash cellHash) const +{ + // check input + uint_t cidx(static_cast(cellHash)); + if ( cidx >= m_towerLookup.size() ) { + ATH_MSG_WARNING( CaloRec::Helpers::fmtMsg("invalid cell hash %6zu beyond range (max hash is %6zu)",cidx,m_maxCellHash) ); + return elementvector_t(); + } else { + return m_towerLookup.at(cidx); + } +} + +//-----------------------// +// Tower Geometry Helper // +//-----------------------// + +double CaloTowerGeometrySvc::cellWeight(IdentifierHash cellHash,index_t towerIdx) const +{ + index_t cidx(static_cast(cellHash)); + double cwght(0.); + + if ( cidx < m_towerLookup.size() ) { + for ( auto elm : m_towerLookup.at(cidx) ) { + if ( towerIndex(elm) == towerIdx ) { cwght = cellWeight(elm); break; } + } + } + return cwght; +} + +//---------------// +// Index Helpers // +//---------------// + +CaloTowerGeometrySvc::index_t CaloTowerGeometrySvc::etaIndex(IdentifierHash cellHash) const +{ + const auto cdde = f_caloDDE(cellHash); + return cdde != 0 ? etaIndex(cdde->eta()) : invalidIndex(); +} + +CaloTowerGeometrySvc::index_t CaloTowerGeometrySvc::etaIndex(double eta) const +{ + return eta >= m_towerEtaMin && eta <= m_towerEtaMax + ? index_t(std::min(static_cast((eta-m_towerEtaMin)/m_towerEtaWidth),m_towerEtaBins-1)) + : invalidIndex(); +} + +CaloTowerGeometrySvc::index_t CaloTowerGeometrySvc::phiIndex(IdentifierHash cellHash) const +{ + const auto cdde = f_caloDDE(cellHash); + return cdde != 0 ? phiIndex(cdde->phi()) : invalidIndex(); +} + +CaloTowerGeometrySvc::index_t CaloTowerGeometrySvc::phiIndex(double phi) const +{ + double dphi(CaloPhiRange::diff(phi,m_towerPhiMin)); + return dphi >= m_towerPhiMin && dphi <= m_towerPhiMax + ? index_t(std::min(static_cast((phi-m_towerPhiMin)/m_towerPhiWidth),m_towerPhiBins-1)) + : invalidIndex(); +} diff --git a/Calorimeter/CaloRec/src/components/CaloRec_entries.cxx b/Calorimeter/CaloRec/src/components/CaloRec_entries.cxx index 4ca13346c76443b425ba96669d377841c2776dcb..90657db39b7455bd8fc9202fc757f3ac8dbfa3f1 100644 --- a/Calorimeter/CaloRec/src/components/CaloRec_entries.cxx +++ b/Calorimeter/CaloRec/src/components/CaloRec_entries.cxx @@ -33,7 +33,9 @@ #include "CaloRec/CaloTowerxAODAlgoBase.h" #include "CaloRec/CaloTopoClusterFromTowerMaker.h" #include "CaloRec/CaloTopoClusterFromTowerCalibrator.h" - +#include "CaloRec/CaloTopoClusterFromTowerMonitor.h" +#include "CaloRec/CaloTowerGeometrySvc.h" +#include "CaloRec/CaloTopoClusterTowerMerger.h" #include "GaudiKernel/DeclareFactoryEntries.h" @@ -57,6 +59,8 @@ DECLARE_ALGORITHM_FACTORY( CaloClusterCellLinksUpdater ) DECLARE_ALGORITHM_FACTORY( CaloTowerxAODFromCells ) DECLARE_ALGORITHM_FACTORY( CaloTowerxAODFromClusters ) DECLARE_ALGORITHM_FACTORY( CaloTowerxAODAlgoBase ) +DECLARE_ALGORITHM_FACTORY( CaloTopoClusterFromTowerMonitor ) + DECLARE_TOOL_FACTORY ( CaloTopoClusterFromTowerMaker ) DECLARE_TOOL_FACTORY ( CaloTopoClusterFromTowerCalibrator ) @@ -96,6 +100,7 @@ DECLARE_FACTORY_ENTRIES(CaloRec) { DECLARE_ALGORITHM( CaloTowerxAODAlgoBase ) DECLARE_ALGORITHM( CaloTowerxAODFromCells ) DECLARE_ALGORITHM( CaloTowerxAODFromClusters ) + DECLARE_ALGORITHM( CaloTopoClusterFromTowerMonitor ) DECLARE_TOOL( CaloTopoClusterFromTowerMaker ) DECLARE_TOOL( CaloTopoClusterFromTowerCalibrator ) @@ -116,3 +121,7 @@ DECLARE_FACTORY_ENTRIES(CaloRec) { DECLARE_TOOL( CaloCellFastCopyTool ) DECLARE_TOOL( CaloClusterSnapshot ) } + +DECLARE_COMPONENT( CaloTowerGeometrySvc ) +DECLARE_COMPONENT( CaloTopoClusterTowerMerger ) + diff --git a/Event/xAOD/xAODCaloEvent/xAODCaloEvent/versions/CaloCluster_v1.h b/Event/xAOD/xAODCaloEvent/xAODCaloEvent/versions/CaloCluster_v1.h index 8e50de066ee8daa7fb4343b9d3ff0377110692fd..02dbf5504a5c379748176d5b55543d688375a083 100644 --- a/Event/xAOD/xAODCaloEvent/xAODCaloEvent/versions/CaloCluster_v1.h +++ b/Event/xAOD/xAODCaloEvent/xAODCaloEvent/versions/CaloCluster_v1.h @@ -84,14 +84,12 @@ namespace xAOD { Topo_633 = 12, // transient cluster for AODCellContainer SW_7_11 = 13, - // cluster representation of towers - Tower_01_01 = 14, - Tower_005_005 = 15, - - - //New (2016) egamma cluster SuperCluster=14, + //New (2019) cluster representation of towers + Tower_01_01 = 15, + Tower_005_005 = 16, + Tower_fixed_area = 17, CSize_Unknown = 99 };