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:
  template <class AddDeadMaterialEnergy>
  void accumulateDeadMaterialEnergy(
      const std::vector<const CaloCalibrationHitContainer*>& v_dmcchc,
      const xAOD::CaloClusterContainer& theClusColl,
      const ClusInfo_t& clusInfoVec,
      const ClusList& clusList,
      bool useParticleID,
      AddDeadMaterialEnergy&& addDeadMaterialEnergy) const
  {
    for (const CaloCalibrationHitContainer* dmcchc : v_dmcchc) {
      for (const CaloCalibrationHit* hit : *dmcchc) {
        const Identifier myId = hit->cellID();
        if (!m_calo_id->is_lar_dm(myId) && !m_calo_id->is_tile_dm(myId)) {
          continue;
        }
 
        const CaloDmDescrElement* myCDDE = m_caloDmDescrManager->get_element(myId);
        if (!myCDDE) {
          continue;
        }
 
        int uniqueID = HepMC::UNDEFINED_ID;
        if (useParticleID) {
          uniqueID = HepMC::uniqueID(hit);
        }
 
        const int jeO = static_cast<int>(std::floor(m_n_eta_out * (myCDDE->eta() / m_out_eta_max)));
        if (jeO < -m_n_eta_out || jeO >= m_n_eta_out) {
          continue;
        }
 
        int jpO = static_cast<int>(std::floor(m_n_phi_out * (myCDDE->phi() / m_out_phi_max)));
        if (jpO < -m_n_phi_out) {
          jpO += 2 * m_n_phi_out;
        }
        if (jpO >= m_n_phi_out) {
          jpO -= 2 * m_n_phi_out;
        }
 
        const int nDmArea = m_caloDmDescrManager->get_dm_area(myId);
        const CaloDmRegion* dmRegion = m_caloDmDescrManager->get_dm_region(myId);
        if (!dmRegion) {
          continue;
        }
 
        std::vector<int> hitClusIndex;
        std::vector<double> hitClusEffEnergy;
        hitClusIndex.reserve(theClusColl.size());
        hitClusEffEnergy.reserve(theClusColl.size());
        double hitClusNorm = 0.0;
 
        const std::vector<int>& matchingClusters =
            clusList[(jpO + m_n_phi_out) * (2 * m_n_eta_out + 1) + jeO + m_n_eta_out];
 
        for (int iClus : matchingClusters) {
          const xAOD::CaloCluster* theCluster = theClusColl.at(iClus);
          const MyClusInfo& clusInfo = clusInfoVec[iClus];
          auto pos = clusInfo.engCalibParticle.find(uniqueID);
          if (pos == clusInfo.engCalibParticle.end()) {
            continue;
          }
 
          const double engClusTruthUniqueIDCalib = pos->second.engTot;
          if (engClusTruthUniqueIDCalib <= m_energyMinCalib || theCluster->e() <= m_energyMin) {
            continue;
          }
 
          double sum_smp_energy = 0.0;
          for (unsigned int i_smp = 0; i_smp < dmRegion->m_CaloSampleNeighbours.size(); ++i_smp) {
            const CaloSampling::CaloSample nsmp =
                static_cast<CaloSampling::CaloSample>(dmRegion->m_CaloSampleNeighbours[i_smp]);
            if ((dmRegion->m_CaloSampleEtaMin[i_smp] - 0.5) <= theCluster->eta() &&
                theCluster->eta() <= (dmRegion->m_CaloSampleEtaMax[i_smp] + 0.5)) {
              sum_smp_energy += pos->second.engSmp[nsmp];
            }
          }
          if (sum_smp_energy <= 0.0) {
            continue;
          }
 
          double phi_diff = myCDDE->phi() - theCluster->phi();
          if (phi_diff <= -M_PI) {
            phi_diff += 2. * M_PI;
          }
          else if (phi_diff > M_PI) {
            phi_diff -= 2. * M_PI;
          }
          const double eta_diff = myCDDE->eta() - theCluster->eta();
          const float distance = std::sqrt(eta_diff * eta_diff + phi_diff * phi_diff);
          const double effEner = std::pow(sum_smp_energy, m_apars_alpha) * std::exp(-distance / m_apars_r0);
 
          hitClusIndex.push_back(iClus);
          hitClusEffEnergy.push_back(effEner);
          hitClusNorm += effEner;
        }
 
        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 dm_weight = hitClusEffEnergy[i] * inv_hitClusNorm;
          addDeadMaterialEnergy(iClus, uniqueID, nDmArea, hit->energyTotal() * dm_weight);
        }
      }
    }
  }
 
  std::vector<std::string> m_momentsNames;
  moment_name_vector m_validNames;
  moment_name_set m_validMoments;
  std::vector<std::string> m_momentsNamesAOD;
  std::set<xAOD::CaloCluster::MomentType> m_momentsAOD;
  SG::ReadHandleKeyArray<CaloCalibrationHitContainer> m_CalibrationHitContainerNames;
  SG::ReadHandleKeyArray<CaloCalibrationHitContainer> m_DMCalibrationHitContainerNames;
  SG::ReadHandleKey<xAOD::TruthParticleContainer> m_truthParticleContainerKey{this,"TruthParticles","TruthParticles","ReadHandleKey for truth particle container"};
  SG::ReadCondHandleKey<CaloDetDescrManager> m_caloDetDescrMgrKey{this,"CaloDetDescrManager", "CaloDetDescrManager"};
  const CaloCell_ID* m_calo_id{};
  const CaloDM_ID* m_caloDM_ID{};
  const CaloDmDescrManager* m_caloDmDescrManager{};
  int m_n_phi_out{};
  int m_n_eta_out{};
  double m_out_phi_max{};
  double m_out_eta_max{};
  double m_rmaxOut[3]{};
  std::array<std::vector<std::vector<CalibHitIPhiIEtaRange> >, 3> m_i_phi_eta;
  mutable std::atomic<bool> m_foundAllContainers{};
  enum keys_dm_energy_sharing {kMatchDmOff, kMatchDmLoose, kMatchDmMedium, kMatchDmTight};
  enum keys_calib_frac_origin {kCalibFracEM, kCalibFracHAD, kCalibFracREST, kCalibFracMax};
  bool m_doDeadEnergySharing{};
  bool m_doOutOfClusterL{};
  bool m_doOutOfClusterM{};
  bool m_doOutOfClusterT{};
  bool m_doDeadL{};
  bool m_doDeadM{};
  bool m_doDeadT{};
  bool m_useParticleID{};
  bool m_doCalibFrac{};
  float m_energyMin{};
  float m_energyMinCalib{};
  float m_apars_alpha{};
  float m_apars_r0{};
  int m_MatchDmType{};
 
  static double angle_mollier_factor(double x);
};
 
#endif // CALOCALIBCLUSTERMOMENTSMAKER2_H