CaloCalibClusterMomentsMaker2.h

Go to the documentation of this file.

/*
  Copyright (C) 2002-2026 CERN for the benefit of the ATLAS collaboration
*/
 
#ifndef CALOCALIBHITREC_CALOCALIBCLUSTERMOMENTSMAKER2_H
#define CALOCALIBHITREC_CALOCALIBCLUSTERMOMENTSMAKER2_H
 
class CaloCell_ID;
class CaloDM_ID;
class CaloDmDescrManager;
class McEventCollection;
class TruthParticleContainer;
 
#include "GaudiKernel/ToolHandle.h"
 
#include "CaloUtils/CaloClusterCollectionProcessor.h"
#include "CaloGeoHelpers/CaloSampling.h"
#include "CaloIdentifier/CaloCell_ID.h"
#include "CaloDmDetDescr/CaloDmDescrManager.h"
#include "CaloSimEvent/CaloCalibrationHitContainer.h"
#include "StoreGate/ReadHandleKeyArray.h"
#include "CaloDetDescr/CaloDetDescrManager.h"
#include "xAODTruth/TruthParticleContainer.h"
#include "CLHEP/Units/SystemOfUnits.h"
#include "TruthUtils/MagicNumbers.h"
 
#include <string>
#include <vector>
#include <set>
#include <map>
#include <atomic>
#include <array>
#include <cmath>
 
class CaloCalibClusterMomentsMaker2: public AthAlgTool, virtual public CaloClusterCollectionProcessor
{
 public:
 
  class CalibHitIPhiIEtaRange {
   public:
    char iPhi, iEtaMin, iEtaMax;
  };
 
  class MyCellInfo : public std::vector<std::pair<int, double> > {
   public:
    MyCellInfo(int iClus, double w) { this->emplace_back(iClus, w); }
    void Add(const MyCellInfo& other) { this->insert(this->end(), other.begin(), other.end()); }
  };
 
  typedef std::map<Identifier, MyCellInfo> CellInfoSet_t;
 
  class MyClusInfo {
   public:
    class ClusCalibEnergy {
     public:
      double engTot = 0.0;
      std::array<double, CaloSampling::Unknown + 1> engSmp{};
      void Add(double eng, int nsmp)
      {
        engTot += eng;
        engSmp[nsmp] += eng;
      }
    };
 
    void Add(double eng, int nsmp, int pid = 0)
    {
      engCalibIn.Add(eng, nsmp);
      engCalibParticle[pid].Add(eng, nsmp);
    }
 
    ClusCalibEnergy engCalibIn;
    double engCalibOut = 0.0;
    double engCalibDead = 0.0;
    std::array<double, CaloDmDescrArea::DMA_MAX> engCalibDeadInArea{};
    std::map<int, ClusCalibEnergy> engCalibParticle{};
  };
  typedef std::vector<MyClusInfo> ClusInfo_t;
  using ClusList = std::vector<std::vector<int> >;
 
  typedef std::pair<std::string, xAOD::CaloCluster::MomentType> moment_name_pair;
  typedef std::vector<moment_name_pair> moment_name_vector;
  typedef std::set<moment_name_pair> moment_name_set;
 
  CaloCalibClusterMomentsMaker2(const std::string& type, const std::string& name,
                                const IInterface* parent);
 
  using CaloClusterCollectionProcessor::execute;
  virtual StatusCode execute(const EventContext& ctx,
                             xAOD::CaloClusterContainer* theClusColl) const override;
  virtual StatusCode initialize() override;

  static void initializeOutOfClusterDistanceTables(
      int n_phi_out,
      int n_eta_out,
      double out_phi_max,
      double out_eta_max,
      const double (&rmaxOut)[3],
      std::array<std::vector<std::vector<CalibHitIPhiIEtaRange> >, 3>& i_phi_eta);
  static void buildCellInfoMap(const xAOD::CaloClusterContainer& theClusColl,
                               CellInfoSet_t& cellInfo);
  
 
  static void buildOutOfClusterClusterLists(
      const xAOD::CaloClusterContainer& theClusColl,
      const ClusInfo_t& clusInfoVec,
      int n_phi_out,
      int n_eta_out,
      double out_phi_max,
      double out_eta_max,
      const std::array<std::vector<std::vector<CalibHitIPhiIEtaRange> >, 3>& i_phi_eta,
      const std::array<bool, 3>& doOutOfCluster,
      const std::array<ClusList*, 3>& clusLists);
  
 
  template <class InvalidUniqueIdHandler>
  static void accumulateClusterCalibHits(
      const std::vector<const CaloCalibrationHitContainer*>& v_cchc,
      const CellInfoSet_t& cellInfo,
      const CaloCell_ID& calo_id,
      ClusInfo_t& clusInfoVec,
      unsigned int& nHitsTotal,
      unsigned int& nHitsWithoutParticleUID,
      bool useParticleID,
      InvalidUniqueIdHandler&& invalidUniqueIdHandler)
  {
    for (const CaloCalibrationHitContainer* cchc : v_cchc) {
      for (const CaloCalibrationHit* hit : *cchc) {
        const Identifier myId = hit->cellID();
        typename CellInfoSet_t::const_iterator pos = cellInfo.find(myId);
 
        if (pos != cellInfo.end()) {
          const CaloSampling::CaloSample nsmp =
              CaloSampling::CaloSample(calo_id.calo_sample(myId));
 
          for (const std::pair<int, double>& p : pos->second) {
            const int iClus = p.first;
            const double weight = p.second;
 
            if (useParticleID) {
              const int uniqueID = HepMC::uniqueID(hit);
              if (uniqueID == HepMC::INVALID_PARTICLE_ID) {
                invalidUniqueIdHandler();
                break;
              }
              clusInfoVec[iClus].Add(weight * hit->energyTotal(),
                                     nsmp,
                                     static_cast<unsigned int>(uniqueID));
            }
            else {
              clusInfoVec[iClus].Add(weight * hit->energyTotal(), nsmp);
            }
          }
        }
 
        if (useParticleID && HepMC::uniqueID(hit) == HepMC::UNDEFINED_ID) {
          ++nHitsWithoutParticleUID;
        }
        ++nHitsTotal;
      }
    }
  }
 
  template <class AddOutOfClusterEnergy>
  static void accumulateOutOfClusterEnergy(
      const std::vector<const CaloCalibrationHitContainer*>& v_cchc,
      const CellInfoSet_t& cellInfo,
      const CaloDetDescrManager& calo_dd_man,
      const ClusInfo_t& clusInfoVec,
      int n_phi_out,
      int n_eta_out,
      double out_phi_max,
      double out_eta_max,
      const std::array<bool, 3>& doOutOfCluster,
      const std::array<const ClusList*, 3>& clusLists,
      AddOutOfClusterEnergy&& addOutOfClusterEnergy)
  {
    for (const CaloCalibrationHitContainer* cchc : v_cchc) {
      for (const CaloCalibrationHit* hit : *cchc) {
        const Identifier myId = hit->cellID();
        typename CellInfoSet_t::const_iterator pos = cellInfo.find(myId);
        if (pos != cellInfo.end()) {
          continue;
        }
 
        const CaloDetDescrElement* myCDDE = calo_dd_man.get_element(myId);
        const int uniqueID = HepMC::uniqueID(hit);
        if (!myCDDE) {
          continue;
        }
 
        const int jeO = static_cast<int>(std::floor(n_eta_out * (myCDDE->eta() / out_eta_max)));
        if (jeO < -n_eta_out || jeO >= n_eta_out) {
          continue;
        }
 
        int jpO = static_cast<int>(std::floor(n_phi_out * (myCDDE->phi() / out_phi_max)));
        if (jpO < -n_phi_out) jpO += 2 * n_phi_out;
        if (jpO >= n_phi_out) jpO -= 2 * n_phi_out;
 
        for (unsigned int ii = 0; ii < 3; ++ii) {
          const ClusList* pClusList = clusLists[ii];
          if (!doOutOfCluster[ii] || pClusList == nullptr) {
            continue;
          }
 
          double hitClusNorm = 0.0;
          std::vector<int> hitClusIndex;
          std::vector<double> hitClusEffEnergy;
          const std::vector<int>& matchingClusters =
              (*pClusList)[(jpO + n_phi_out) * (2 * n_eta_out + 1) + jeO + n_eta_out];
 
          for (int iClus : matchingClusters) {
            const MyClusInfo& clusInfo = clusInfoVec[iClus];
            auto partPos = clusInfo.engCalibParticle.find(uniqueID);
            if (partPos == clusInfo.engCalibParticle.end()) {
              continue;
            }
            hitClusNorm += partPos->second.engTot;
            hitClusEffEnergy.push_back(partPos->second.engTot);
            hitClusIndex.push_back(iClus);
          }
 
          if (hitClusNorm <= 0.0) {
            continue;
          }
 
          const double inv_hitClusNorm = 1.0 / hitClusNorm;
          for (std::size_t i = 0; i < hitClusIndex.size(); ++i) {
            const int iClus = hitClusIndex[i];
            const double w = hitClusEffEnergy[i] * inv_hitClusNorm;
            addOutOfClusterEnergy(ii, iClus, uniqueID, w * hit->energyTotal());
          }
        }
      }
    }
  }
 
 private:
  // Properties
  SG::ReadHandleKeyArray<CaloCalibrationHitContainer> m_CalibrationHitContainerNames{this, "CalibrationHitContainerNames", {}};
  SG::ReadHandleKeyArray<CaloCalibrationHitContainer> m_DMCalibrationHitContainerNames{this, "DMCalibrationHitContainerNames", {}};
  SG::ReadHandleKey<xAOD::TruthParticleContainer> m_truthParticleContainerKey{this,"TruthParticles","TruthParticles","ReadHandleKey for truth particle container"};
  SG::ReadCondHandleKey<CaloDetDescrManager> m_caloDetDescrMgrKey{this,"CaloDetDescrManager", "CaloDetDescrManager"};
  
  Gaudi::Property<std::vector<std::string>> m_momentsNames{this, "MomentsNames", {}};
  Gaudi::Property<int> m_MatchDmType{this, "MatchDmType", kMatchDmLoose};
  Gaudi::Property<bool> m_useParticleID{this, "UseParticleID", true};
 
  // Name(s) of Moments which can be stored on the AOD - all others go to ESD
  Gaudi::Property<std::vector<std::string>> m_momentsNamesAOD {this, "AODMomentsNames",
                                   {"ENG_CALIB_TOT"
                                , "ENG_CALIB_OUT_L"
                                , "ENG_CALIB_OUT_M"
                                , "ENG_CALIB_OUT_T"
                                , "ENG_CALIB_EMB0"
                                , "ENG_CALIB_EME0"
                                , "ENG_CALIB_TILEG3"
                                , "ENG_CALIB_DEAD_TOT"
                                , "ENG_CALIB_DEAD_EMB0"
                                , "ENG_CALIB_DEAD_TILE0"
                                , "ENG_CALIB_DEAD_TILEG3"
                                , "ENG_CALIB_DEAD_EME0"
                                , "ENG_CALIB_DEAD_HEC0"
                                , "ENG_CALIB_DEAD_FCAL"
                                , "ENG_CALIB_DEAD_LEAKAGE"
                                , "ENG_CALIB_DEAD_UNCLASS"
                                , "ENG_CALIB_FRAC_EM"
                                , "ENG_CALIB_FRAC_HAD"
                                , "ENG_CALIB_FRAC_REST"
                                   }
  };
 
 
  // Name(s) of Moments to calculate
  moment_name_vector m_validNames{std::make_pair(std::string("ENG_CALIB_TOT"), xAOD::CaloCluster::ENG_CALIB_TOT)
    , std::make_pair(std::string("ENG_CALIB_OUT_L"), xAOD::CaloCluster::ENG_CALIB_OUT_L)
    , std::make_pair(std::string("ENG_CALIB_OUT_M"), xAOD::CaloCluster::ENG_CALIB_OUT_M)
    , std::make_pair(std::string("ENG_CALIB_OUT_T"), xAOD::CaloCluster::ENG_CALIB_OUT_T)
    , std::make_pair(std::string("ENG_CALIB_EMB0"), xAOD::CaloCluster::ENG_CALIB_EMB0)
    , std::make_pair(std::string("ENG_CALIB_EME0"), xAOD::CaloCluster::ENG_CALIB_EME0)
    , std::make_pair(std::string("ENG_CALIB_TILEG3"), xAOD::CaloCluster::ENG_CALIB_TILEG3)
    , std::make_pair(std::string("ENG_CALIB_DEAD_TOT"), xAOD::CaloCluster::ENG_CALIB_DEAD_TOT)
    , std::make_pair(std::string("ENG_CALIB_DEAD_EMB0"), xAOD::CaloCluster::ENG_CALIB_DEAD_EMB0)
    , std::make_pair(std::string("ENG_CALIB_DEAD_TILE0"), xAOD::CaloCluster::ENG_CALIB_DEAD_TILE0)
    , std::make_pair(std::string("ENG_CALIB_DEAD_TILEG3"), xAOD::CaloCluster::ENG_CALIB_DEAD_TILEG3)
    , std::make_pair(std::string("ENG_CALIB_DEAD_EME0"), xAOD::CaloCluster::ENG_CALIB_DEAD_EME0)
    , std::make_pair(std::string("ENG_CALIB_DEAD_HEC0"), xAOD::CaloCluster::ENG_CALIB_DEAD_HEC0)
    , std::make_pair(std::string("ENG_CALIB_DEAD_FCAL"), xAOD::CaloCluster::ENG_CALIB_DEAD_FCAL)
    , std::make_pair(std::string("ENG_CALIB_DEAD_LEAKAGE"), xAOD::CaloCluster::ENG_CALIB_DEAD_LEAKAGE)
    , std::make_pair(std::string("ENG_CALIB_DEAD_UNCLASS"), xAOD::CaloCluster::ENG_CALIB_DEAD_UNCLASS)
    , std::make_pair(std::string("ENG_CALIB_FRAC_EM"), xAOD::CaloCluster::ENG_CALIB_FRAC_EM)
    , std::make_pair(std::string("ENG_CALIB_FRAC_HAD"), xAOD::CaloCluster::ENG_CALIB_FRAC_HAD)
    , std::make_pair(std::string("ENG_CALIB_FRAC_REST"), xAOD::CaloCluster::ENG_CALIB_FRAC_REST)
  };
 
  moment_name_set m_validMoments;
 
  std::set<xAOD::CaloCluster::MomentType> m_momentsAOD;
  const CaloCell_ID* m_calo_id{};
  const CaloDM_ID* m_caloDM_ID{};
  const CaloDmDescrManager* m_caloDmDescrManager{};
  int m_n_phi_out{127};
  int m_n_eta_out{127};
  double m_out_phi_max{M_PI};
  double m_out_eta_max{6};
  double m_rmaxOut[3] = {1.0, 0.5, 0.3};
  std::array<std::vector<std::vector<CalibHitIPhiIEtaRange> >, 3> m_i_phi_eta;
  mutable std::atomic<bool> m_foundAllContainers{false};
  enum keys_dm_energy_sharing {kMatchDmOff, kMatchDmLoose, kMatchDmMedium, kMatchDmTight};
  enum keys_calib_frac_origin {kCalibFracEM, kCalibFracHAD, kCalibFracREST, kCalibFracMax};
  bool m_doDeadEnergySharing{false};
  bool m_doOutOfClusterL{false};
  bool m_doOutOfClusterM{false};
  bool m_doOutOfClusterT{false};
  bool m_doDeadL{false};
  bool m_doDeadM{false};
  bool m_doDeadT{false};
 
  bool m_doCalibFrac{false};
  float m_energyMin{200*CLHEP::MeV};
  float m_energyMinCalib{20*CLHEP::MeV};
  float m_apars_alpha{0.5};
  float m_apars_r0{0.2};
 
  static double angle_mollier_factor(double x);
};
 
#endif // CALOCALIBCLUSTERMOMENTSMAKER2_H