egammaSuperClusterBuilderBase.cxx

Go to the documentation of this file.

/*
  Copyright (C) 2002-2023 CERN for the benefit of the ATLAS collaboration
*/
 
#include "egammaSuperClusterBuilderBase.h"
 
#include "CaloUtils/CaloCellList.h"
#include "CaloUtils/CaloClusterStoreHelper.h"
#include "CaloUtils/CaloLayerCalculator.h"
 
#include "egammaRecEvent/egammaRecContainer.h"
#include "xAODCaloEvent/CaloCluster.h"
#include "xAODCaloEvent/CaloClusterAuxContainer.h"
#include "xAODCaloEvent/CaloClusterKineHelper.h"
#include "xAODEgamma/EgammaEnums.h"
#include "xAODEgamma/EgammaxAODHelpers.h"
#include "xAODEgamma/PhotonxAODHelpers.h"
#include "xAODTracking/TrackParticle.h"
#include "xAODTracking/Vertex.h"
 
#include "CaloGeoHelpers/proxim.h"
#include "FourMomUtils/P4Helpers.h"
 
#include <cmath>
#include <optional>
#include <utility>
 
using xAOD::EgammaHelpers::summaryValueInt;
 
namespace {
std::pair<const double, const double>
etaphi_range(const CaloDetDescrManager& dd_man,
             double eta,
             double phi,
             CaloCell_ID::CaloSample sampling,
             const CaloDetDescrElement* elt)
{
  // Should be smaller than the eta half-width of any cell.
  constexpr double eps = 0.001;
 
  // Now look in the negative eta direction, on the low (left) side.
  const CaloDetDescrElement* elt_l = dd_man.get_element_raw(sampling, eta - elt->deta() - eps, phi);
  double deta_l = elt_l ? std::abs(eta - elt_l->eta_raw()) + eps : 0.;
 
  // Now look in the positive eta direction, on the high (right) side.
  const CaloDetDescrElement* elt_r = dd_man.get_element_raw(sampling, eta + elt->deta() + eps, phi);
  double deta_r = elt_r ? std::abs(eta - elt_r->eta_raw()) + eps : 0.; 
 
  // Now for the phi variation.
  // The phi size can change as a function of eta, but not of phi.
  // Thus we have to look again at the adjacent eta cells, and
  // take the largest variation.
 
  // Now look in the negative eta direction, on the low-eta () side.
  elt_l = dd_man.get_element_raw(sampling, eta - elt->deta() - eps, CaloPhiRange::fix(phi - elt->dphi() - eps));
  double dphi_l = elt_l ? std::abs(CaloPhiRange::fix(phi - elt_l->phi_raw())) + eps : 0.; 
 
  // Now look in the positive eta direction, on the positive (down) side.
  elt_r = dd_man.get_element_raw(sampling, eta + elt->deta() + eps, CaloPhiRange::fix(phi - elt->dphi() - eps));
  double dphi_r = elt_r ? std::abs(CaloPhiRange::fix(phi - elt_r->phi_raw())) + eps : 0.; 
  
  // Total is twice the maximum.
  return {2 * std::max(deta_r, deta_l), 2 * std::max(dphi_l, dphi_r)};
}
 
void
fillPositionsInCalo(xAOD::CaloCluster* cluster, const CaloDetDescrManager& mgr)
{
  const bool isBarrel = xAOD::EgammaHelpers::isBarrel(cluster);
  CaloCell_ID::CaloSample sample =
    isBarrel ? CaloCell_ID::EMB2 : CaloCell_ID::EME2;
  // eta and phi of the cluster in the calorimeter frame
  double eta;
  double phi;
  CaloCellDetPos::getDetPosition(
    mgr, sample, cluster->eta(), cluster->phi(), eta, phi);
  cluster->insertMoment(xAOD::CaloCluster::ETACALOFRAME, eta);
  cluster->insertMoment(xAOD::CaloCluster::PHICALOFRAME, phi);
  //  eta in the second sampling
  CaloCellDetPos::getDetPosition(
    mgr, sample, cluster->etaBE(2), cluster->phiBE(2), eta, phi);
  cluster->insertMoment(xAOD::CaloCluster::ETA2CALOFRAME, eta);
  cluster->insertMoment(xAOD::CaloCluster::PHI2CALOFRAME, phi);
  //  eta in the first sampling
  sample = isBarrel ? CaloCell_ID::EMB1 : CaloCell_ID::EME1;
  CaloCellDetPos::getDetPosition(
    mgr, sample, cluster->etaBE(1), cluster->phiBE(1), eta, phi);
  cluster->insertMoment(xAOD::CaloCluster::ETA1CALOFRAME, eta);
  cluster->insertMoment(xAOD::CaloCluster::PHI1CALOFRAME, phi);
}
void
makeCorrection1(xAOD::CaloCluster* cluster,
                const CaloDetDescrManager& mgr,
                const CaloSampling::CaloSample sample)
{
  const double clusterEtaMax = cluster->etamax(sample);
  const double clusterPhiMax = cluster->phimax(sample);
 
  // Protections.
  if (clusterEtaMax == -999. || clusterPhiMax == -999.) {
    return;
  }
  if (std::abs(clusterEtaMax) < 1E-6 && std::abs(clusterPhiMax) < 1E-6) {
    return;
  }
 
  // Get the hottest in raw co-ordinates
  // We have two kinds of enums ...
  const CaloCell_ID::CaloSample xsample = 
    (sample == CaloSampling::EMB1) ? CaloCell_ID::EMB1 : CaloCell_ID::EME1;
  
  const CaloDetDescrElement* dde = mgr.get_element(xsample, 
                                                   clusterEtaMax, 
                                                   clusterPhiMax);
 
  if (!dde) {
    return;
  }
  
  double etamax = dde->eta_raw();
  double phimax = dde->phi_raw();
 
  const CaloDetDescrElement* elt = mgr.get_element_raw(xsample, etamax, phimax);
  if (!elt) {
    return;
  }
 
  // Now Locate the +-1 range, use raw co-ordinates here.
  auto [detastr, dphistr] = etaphi_range(mgr, etamax, phimax, xsample, elt);
 
  // Given the range refine the position employing the smaller window
  if (detastr > 0 && dphistr > 0) {
    CaloLayerCalculator helper;
    const auto* const cellLink = cluster->getCellLinks();
    helper.fill(cellLink->begin(),
                cellLink->end(),
                etamax,
                phimax,
                detastr,
                dphistr,
                sample);
 
    // Here is where we (re-)fill the eta in the 1st sampling
    if (helper.etam() != -999.) {
      // This is "real" atlas co-ordinates
      cluster->setEta(sample, helper.etam());
    }
  }
}
 
void
refineEta1Position(xAOD::CaloCluster* cluster, const CaloDetDescrManager& mgr)
{
  // This only makes sense if we have cells there
  if (!cluster->hasSampling(CaloSampling::EMB1) &&
      !cluster->hasSampling(CaloSampling::EME1)) {
    return;
  }
  // Now calculare the position using cells in barrel or endcap or both
  const double aeta = std::abs(cluster->etaBE(2));
  if (aeta < 1.6 && cluster->hasSampling(CaloSampling::EMB1)) {
    makeCorrection1(cluster, mgr, CaloSampling::EMB1);
  }
  if (aeta > 1.3 && cluster->hasSampling(CaloSampling::EME1)) {
    makeCorrection1(cluster, mgr, CaloSampling::EME1);
  }
}
 
} // end of anonymous namespace
 
/*
 * Gaudi Algorithm implementation
 */
egammaSuperClusterBuilderBase::egammaSuperClusterBuilderBase(
  const std::string& name,
  ISvcLocator* pSvcLocator)
  : AthReentrantAlgorithm(name, pSvcLocator)
{
}
 
StatusCode
egammaSuperClusterBuilderBase::initialize()
{
  ATH_CHECK(m_inputEgammaRecContainerKey.initialize());
  ATH_CHECK(m_precorrClustersKey.initialize(SG::AllowEmpty));
  ATH_CHECK(m_caloDetDescrMgrKey.initialize());
  ATH_CHECK(m_outputEgammaRecContainerKey.initialize());
  ATH_CHECK(m_outputSuperClusterCollectionName.initialize());
 
  m_addCellsWindowEtaBarrel =
    m_addCellsWindowEtaCellsBarrel * s_cellEtaSize * 0.5;
  m_addCellsWindowEtaEndcap =
    m_addCellsWindowEtaCellsEndcap * s_cellEtaSize * 0.5;
  // The +- to account for the different L3 eta granularity
  m_extraL3EtaSize = m_extraL3EtaSizeCells * s_cellEtaSize * 0.5;
  // the + is to account for different L0/L1 phi granularity
  m_extraL0L1PhiSize = m_extraL0L1PhiSizeCells * s_cellPhiSize;
 
  m_searchWindowPhiBarrel = m_searchWindowPhiCellsBarrel * s_cellPhiSize * 0.5;
  m_searchWindowEtaBarrel = m_searchWindowEtaCellsBarrel * s_cellEtaSize * 0.5;
  m_searchWindowPhiEndcap = m_searchWindowPhiCellsEndcap * s_cellPhiSize * 0.5;
  m_searchWindowEtaEndcap = m_searchWindowEtaCellsEndcap * s_cellEtaSize * 0.5;
 
  if (m_addCellsWindowEtaCellsBarrel % 2 == 0 ||
      m_addCellsWindowEtaCellsEndcap % 2 == 0) {
    ATH_MSG_FATAL("For adding cells relative to the hottest cell to be "
                  "symmetric in eta, the AddCells "
                  "window size needs to be odd");
 
    return StatusCode::FAILURE;
  }
  ATH_CHECK(m_clusterCorrectionTool.retrieve());
  ATH_CHECK(m_MVACalibSvc.retrieve());
 
  if (!m_egammaCheckEnergyDepositTool.empty()) {
    ATH_CHECK(m_egammaCheckEnergyDepositTool.retrieve());
  } else {
    m_egammaCheckEnergyDepositTool.disable();
  }
  if (!m_egammaCellRecoveryTool.empty()) {
    ATH_CHECK(m_egammaCellRecoveryTool.retrieve());
  } else {
    m_egammaCellRecoveryTool.disable();
  }
 
  ATH_MSG_INFO("e/gamma super clusters"
               << '\n'
               << "--> Eta Window size for L0/L1/L2 cells : "
               << " Barrel +- " << m_addCellsWindowEtaBarrel << " ,EndCap +- "
               << m_addCellsWindowEtaEndcap << '\n'
               << "--> Eta Window size for L3 cells : "
               << " Barrel +- "
               << (m_addCellsWindowEtaBarrel + m_extraL3EtaSize)
               << " ,EndCap +- "
               << (m_addCellsWindowEtaEndcap + m_extraL3EtaSize) << '\n'
               << " -> Phi window is fully dynamic for L2/L3" << '\n'
               << " -> L0/L1 cells in phi will be collected in a Window : "
               << "(L2 neg extend - " << m_extraL0L1PhiSize << " , "
               << "L2 pos extend + " << m_extraL0L1PhiSize << ")");
 
  return StatusCode::SUCCESS;
}
 
StatusCode
egammaSuperClusterBuilderBase::execute(const EventContext& ctx) const
{
  SG::ReadHandle<EgammaRecContainer> egammaRecs(m_inputEgammaRecContainerKey,
                                                ctx);
  // check is only used for serial running; remove when MT scheduler used
  ATH_CHECK(egammaRecs.isValid());
 
  // Have to register cluster container in order to properly get cluster
  // links.
  SG::WriteHandle<xAOD::CaloClusterContainer> outputClusterContainer(
    m_outputSuperClusterCollectionName, ctx);
 
  ATH_CHECK(CaloClusterStoreHelper::AddContainerWriteHandle(outputClusterContainer));
 
  // Create the new Electron Super Cluster based EgammaRecContainer
  SG::WriteHandle<EgammaRecContainer> newEgammaRecs(
    m_outputEgammaRecContainerKey, ctx);
  ATH_CHECK(newEgammaRecs.record(std::make_unique<EgammaRecContainer>()));
 
  size_t inputSize = egammaRecs->size();
  outputClusterContainer->reserve(inputSize);
  newEgammaRecs->reserve(inputSize);
 
  std::optional<SG::WriteHandle<xAOD::CaloClusterContainer>> precorrClustersH;
  if (!m_precorrClustersKey.empty()) {
    precorrClustersH.emplace(m_precorrClustersKey, ctx);
    ATH_CHECK(CaloClusterStoreHelper::AddContainerWriteHandle(*precorrClustersH));
    precorrClustersH->ptr()->reserve(inputSize);
  }
 
  // The calo Det Descr manager
  SG::ReadCondHandle<CaloDetDescrManager> caloDetDescrMgrHandle{
    m_caloDetDescrMgrKey, ctx
  };
  ATH_CHECK(caloDetDescrMgrHandle.isValid());
  const CaloDetDescrManager* calodetdescrmgr = *caloDetDescrMgrHandle;
 
  // If no input return
  if (egammaRecs->empty()) {
    return StatusCode::SUCCESS;
  }
 
  // Figure the cellCont we need to point to
  const DataLink<CaloCellContainer>& cellCont =
    (*egammaRecs)[0]->caloCluster()->getCellLinks()->getCellContainerLink();
 
  // Loop over input egammaRec objects, build superclusters.
  size_t numInput = egammaRecs->size();
  std::vector<bool> isUsed(numInput, false);
  std::vector<bool> isUsedRevert(numInput, false);
 
  for (std::size_t i = 0; i < egammaRecs->size(); ++i) {
    if (isUsed[i]) {
      continue;
    }
    const auto* egRec = (*egammaRecs)[i];
    // check for good seed cluster
    const xAOD::CaloCluster* clus = egRec->caloCluster();
    if (!seedClusterSelection(clus)) {
      continue;
    }
    
    if (!egammaRecPassesSelection(egRec)) {
      continue; 
    }
 
    // save status in case we fail to create supercluster
    isUsedRevert = isUsed;
    // Mark seed as used,
    isUsed[i] = true;
 
    // Start accumulating the clusters from the seed
    std::vector<const xAOD::CaloCluster*> accumulatedClusters;
    accumulatedClusters.push_back(clus);
 
    // Now we find all the secondary cluster for this seed
    // and we accumulate them
    const std::vector<std::size_t> secondaryIndices =
      searchForSecondaryClusters(i, egammaRecs.cptr(), isUsed);
    for (const auto& secClusIndex : secondaryIndices) {
      const auto* const secRec = (*egammaRecs)[secClusIndex];
      accumulatedClusters.push_back(secRec->caloCluster());
    }
    ATH_MSG_DEBUG("Total clusters " << accumulatedClusters.size());
 
    // Create the new cluster
    std::unique_ptr<xAOD::CaloCluster> newCluster = 
      createNewCluster(ctx,
                       accumulatedClusters,
                       cellCont,
                       *calodetdescrmgr,
                       getEgammaRecType(egRec),
                       precorrClustersH ? precorrClustersH->ptr() : nullptr);
 
    // If we failed to create a cluster revert isUsed for the cluster
    if (newCluster) {
      outputClusterContainer->push_back(std::move(newCluster));
    }
    else {
      isUsed.swap(isUsedRevert);
      continue;
    }
      
    // Add the cluster links to the super cluster
    ElementLink<xAOD::CaloClusterContainer> clusterLink(
      *outputClusterContainer, outputClusterContainer->size() - 1, ctx);
    std::vector<ElementLink<xAOD::CaloClusterContainer>> elClusters{
      clusterLink
    };
  
    // Make egammaRec object, and push it back into output container.
    auto newEgRec = std::make_unique<egammaRec>(*egRec);
    newEgRec->setCaloClusters(elClusters);
    newEgammaRecs->push_back(std::move(newEgRec));
  } // End loop on egammaRecs
 
  ATH_CHECK(redoMatching(ctx, newEgammaRecs));
 
  return StatusCode::SUCCESS;
}
 
bool 
egammaSuperClusterBuilderBase::egammaRecPassesSelection([[maybe_unused]] const egammaRec *egRec) const { 
  return true;
}
 
StatusCode 
egammaSuperClusterBuilderBase::redoMatching(
  [[maybe_unused]] const EventContext &ctx, 
  [[maybe_unused]] SG::WriteHandle<EgammaRecContainer> &newEgammaRecs
) const {
  return StatusCode::SUCCESS;
}
 
bool
egammaSuperClusterBuilderBase::matchesInWindow(
  const xAOD::CaloCluster* ref,
  const xAOD::CaloCluster* clus) const
{
  auto inRange = [](float eta1, float phi1, 
                    float eta2, float phi2,
                    float etaWindow, float phiWindow) {
 
    const float dEta = std::abs(eta1 - eta2);
    const float dPhi = std::abs(P4Helpers::deltaPhi(phi1, phi2));
 
    return dEta < etaWindow && dPhi < phiWindow;
  };
 
  // First the case where the seed is both endcap and barrel, i.e. in the crack
  // Check around both measurements of the seed
  if (ref->hasSampling(CaloSampling::EMB2) &&
      ref->hasSampling(CaloSampling::EME2)) {
    const bool inRangeBarrel = inRange(ref->eta(), 
                                       ref->phi(), 
                                       clus->eta(),
                                       clus->phi(),
                                       m_searchWindowEtaBarrel, 
                                       m_searchWindowPhiBarrel);
 
    const bool inRangeEndcap = inRange(ref->eta(), 
                                       ref->phi(), 
                                       clus->eta(), 
                                       clus->phi(),
                                       m_searchWindowEtaEndcap, 
                                       m_searchWindowPhiEndcap);
 
    const bool inRangeBarrelL2 = inRange(ref->etaSample(CaloSampling::EMB2), 
                                         ref->phiSample(CaloSampling::EMB2),
                                         clus->eta(), 
                                         clus->phi(),
                                         m_searchWindowEtaBarrel, 
                                         m_searchWindowPhiBarrel);
 
    const bool inRangeEndcapL2 = inRange(ref->etaSample(CaloSampling::EME2), 
                                         ref->phiSample(CaloSampling::EME2),
                                         clus->eta(), 
                                         clus->phi(),
                                         m_searchWindowEtaEndcap, 
                                         m_searchWindowPhiEndcap);
    
    // Matches any in case of split
    return inRangeBarrel || inRangeEndcap || inRangeBarrelL2 || inRangeEndcapL2;
  }
 
  if (xAOD::EgammaHelpers::isBarrel(clus)) {
    return inRange(ref->eta(), 
                   ref->phi(), 
                   clus->eta(), 
                   clus->phi(),
                   m_searchWindowEtaBarrel, 
                   m_searchWindowPhiBarrel);
  }
 
  return inRange(ref->eta(), 
                 ref->phi(), 
                 clus->eta(), 
                 clus->phi(),
                 m_searchWindowEtaEndcap, 
                 m_searchWindowPhiEndcap);
}
 
std::unique_ptr<xAOD::CaloCluster>
egammaSuperClusterBuilderBase::createNewCluster(
  const EventContext& ctx,
  const std::vector<const xAOD::CaloCluster*>& clusters,
  const DataLink<CaloCellContainer>& cellCont,
  const CaloDetDescrManager& mgr,
  xAOD::EgammaParameters::EgammaType egType,
  xAOD::CaloClusterContainer* precorrClusters) const
{
  if (clusters.empty()) {
    ATH_MSG_ERROR("Missing the seed cluster! Should not happen.");
    return nullptr;
  }
 
  // create a new empty cluster
  // collection will own it if
  auto newCluster = CaloClusterStoreHelper::makeCluster(cellCont);
 
  if (!newCluster) {
    ATH_MSG_ERROR("CaloClusterStoreHelper::makeCluster failed.");
    return nullptr;
  }
  //
  newCluster->setClusterSize(xAOD::CaloCluster::SuperCluster);
  // Let's try to find the eta and phi of the hottest cell in L2.
  // This will be used as the center for restricting the cluster size.
  CookieCutterHelpers::CentralPosition cp0(clusters, mgr);
 
  // Set the eta0/phi0 based on the references, but in raw coordinates
  if (cp0.emaxB >= cp0.emaxEC) {
    newCluster->setEta0(cp0.etaB);
    newCluster->setPhi0(cp0.phiB);
  } else {
    newCluster->setEta0(cp0.etaEC);
    newCluster->setPhi0(cp0.phiEC);
  }
 
  // Actually fill the cluster here
  if (fillClusterConstrained(*newCluster, clusters, cp0).isFailure()) {
    return nullptr;
  }
  // Apply SW-style summation of TileGap3 cells (if necessary).
  float eta0 = std::abs(newCluster->eta0());
  // In Run2, we did not impose restriction to include TG3 cells at this level.
  // It should have been [1.37,1.63]. It has no impact on performance as TG3 was
  // only used in energy calibration BDT in [1.4,1.6].
  // In Run three we restrict to [1.37,1.75]
  if (!m_useExtendedTG3 ||
      (eta0 > s_ClEtaMinForTG3cell && eta0 < s_ClEtaMaxForTG3cell)) {
    if (addTileGap3CellsinWindow(*newCluster, mgr).isFailure()) {
      ATH_MSG_ERROR("Problem with the input cluster when running "
                    "AddTileGap3CellsinWindow?");
 
      return nullptr;
    }
  }
  CaloClusterKineHelper::calculateKine(newCluster.get(), true, true);
 
  // If adding all EM cells we are somehow below the seed threshold then remove
  if (newCluster->et() < m_EtThresholdCut) {
    return nullptr;
  }
 
  // Check to see if cluster pases basic requirements. If not, kill it.
  if (!m_egammaCheckEnergyDepositTool.empty() &&
      !m_egammaCheckEnergyDepositTool->checkFractioninSamplingCluster(
        newCluster.get())) {
    return nullptr;
  }
 
  // Apply correction calibration
  if (calibrateCluster(ctx, newCluster.get(), mgr, egType, precorrClusters)
        .isFailure()) {
    ATH_MSG_WARNING("There was problem calibrating the object");
    return nullptr;
  }
 
  // Avoid negative energy clusters
  if (newCluster->et() < 0) {
    return nullptr;
  }
 
  if (m_linkToConstituents) {
    // EDM vector to constituent clusters
    std::vector<ElementLink<xAOD::CaloClusterContainer>> constituentLinks;
    for (const xAOD::CaloCluster* cluster : clusters) {
      ElementLink<xAOD::CaloClusterContainer> sisterCluster = 
        cluster->getSisterClusterLink();
 
      // Set the element Link to the constitents
      if (sisterCluster) {
        constituentLinks.push_back(sisterCluster);
      } else {
        ATH_MSG_WARNING("No sister Link available");
      }
    }
    // Set the link from the super cluster to the constituents (accumulated)
    // clusters used.
    static const SG::AuxElement::Accessor<
      std::vector<ElementLink<xAOD::CaloClusterContainer>>>
      caloClusterLinks("constituentClusterLinks");
    caloClusterLinks(*newCluster) = constituentLinks;
  }
  // return the new cluster
  return newCluster;
}
 
bool
egammaSuperClusterBuilderBase::seedClusterSelection(
  const xAOD::CaloCluster* clus) const
{
  // The seed should have 2nd sampling
  if (!clus->hasSampling(CaloSampling::EMB2) &&
      !clus->hasSampling(CaloSampling::EME2)) {
    return false;
  }
  const double eta2 = std::abs(clus->etaBE(2));
  if (eta2 > 10) {
    return false;
  }
  // Accordeon Energy samplings 1 to 3
  const double EMAccEnergy =
    clus->energyBE(1) + clus->energyBE(2) + clus->energyBE(3);
  const double EMAccEt = EMAccEnergy / cosh(eta2);
  // Require minimum energy for supercluster seeding.
  return EMAccEt >= m_EtThresholdCut;
}
 
StatusCode
egammaSuperClusterBuilderBase::fillClusterConstrained(
  xAOD::CaloCluster& tofill,
  const std::vector<const xAOD::CaloCluster*>& clusters,
  const CookieCutterHelpers::CentralPosition& cp0) const
{
  const float addCellsWindowEtaBarrel = m_addCellsWindowEtaBarrel;
  const float addCellsWindowEtaEndcap = m_addCellsWindowEtaEndcap;
  const float addCellsWindowL3EtaBarrel =
    m_addCellsWindowEtaBarrel + m_extraL3EtaSize;
  const float addCellsWindowL3EtaEndcap =
    m_addCellsWindowEtaEndcap + m_extraL3EtaSize;
 
  // Loop for L2/L3
  for (const xAOD::CaloCluster* tocheck : clusters) {
    xAOD::CaloCluster::const_cell_iterator cell_itr = tocheck->begin();
    xAOD::CaloCluster::const_cell_iterator cell_end = tocheck->end();
    // Loop over cells
    for (; cell_itr != cell_end; ++cell_itr) {
      // sanity check on the cell
      const CaloCell* cell = *cell_itr;
      if (!cell) {
        continue;
      }
      const CaloDetDescrElement* dde = cell->caloDDE();
      if (!dde) {
        continue;
      }
      // we want only LAREM
      if (!(dde->getSubCalo() == CaloCell_ID::LAREM)) {
        continue;
      }
      // we want L2 or L3 cells
      const auto sampling = dde->getSampling();
      const bool isL2Cell =
        (CaloCell_ID::EMB2 == sampling || CaloCell_ID::EME2 == sampling);
      const bool isL3Cell =
        (CaloCell_ID::EMB3 == sampling || CaloCell_ID::EME3 == sampling);
 
      if ((!isL2Cell) && (!isL3Cell)) {
        continue;
      }
      // Also exclude the inner wheel Endcap
      if (dde->is_lar_em_endcap_inner()) {
        continue;
      }
 
      bool inEtaRange = false;
      // Check if is inside the eta range wrt to the hottest
      // cell(s) for the cluster we construct
      if (cp0.emaxB > 0) { // barrel
        if (isL2Cell &&
            (std::abs(cp0.etaB - dde->eta_raw()) < addCellsWindowEtaBarrel)) {
          inEtaRange = true;
        }
        if (isL3Cell &&
            (std::abs(cp0.etaB - dde->eta_raw()) < addCellsWindowL3EtaBarrel)) {
          inEtaRange = true;
        }
      }
      if (cp0.emaxEC > 0) { // endcap
        if (isL2Cell &&
            (std::abs(cp0.etaEC - dde->eta_raw()) < addCellsWindowEtaEndcap)) {
          inEtaRange = true;
        }
        if (isL3Cell && (std::abs(cp0.etaEC - dde->eta_raw()) <
                         addCellsWindowL3EtaEndcap)) {
          inEtaRange = true;
        }
      }
      if (!inEtaRange) {
        continue;
      }
      tofill.addCell(cell_itr.index(), cell_itr.weight());
    } // Loop over cells for L2/L3
  }   // Loop over clusters for L2/L3
 
  // We should have a size here
  if (tofill.size() == 0) {
    return StatusCode::FAILURE;
  }
  // Now calculate the cluster size in 2nd layer
  // use that for constraining the L0/L1 cells we add
  const CookieCutterHelpers::PhiSize phiSize(cp0, tofill);
  const float phiPlusB = cp0.phiB + phiSize.plusB + m_extraL0L1PhiSize;
  const float phiMinusB = cp0.phiB - phiSize.minusB - m_extraL0L1PhiSize;
  const float phiPlusEC = cp0.phiEC + phiSize.plusEC + m_extraL0L1PhiSize;
  const float phiMinusEC = cp0.phiEC - phiSize.minusEC - m_extraL0L1PhiSize;
 
  // Loop for L0/L1
  for (const xAOD::CaloCluster* tocheck : clusters) {
    xAOD::CaloCluster::const_cell_iterator cell_itr = tocheck->begin();
    xAOD::CaloCluster::const_cell_iterator cell_end = tocheck->end();
    // Loop over cells
    for (; cell_itr != cell_end; ++cell_itr) {
      // sanity check on the cell
      const CaloCell* cell = *cell_itr;
      if (!cell) {
        continue;
      }
      const CaloDetDescrElement* dde = cell->caloDDE();
      if (!dde) {
        continue;
      }
 
      // only deal with L1 or PS
      const auto sampling = dde->getSampling();
      const bool isL0L1Cell =
        (CaloCell_ID::EMB1 == sampling || CaloCell_ID::EME1 == sampling ||
         CaloCell_ID::PreSamplerB == sampling ||
         CaloCell_ID::PreSamplerE == sampling);
      if (!isL0L1Cell) {
        continue;
      }
 
      bool inEtaRange = false;
      // Check if is inside the eta range wrt to the hottest
      // cell(s) for the cluster we construct
      if (cp0.emaxB > 0) { // barrel
        if (std::abs(cp0.etaB - dde->eta_raw()) < addCellsWindowEtaBarrel) {
          inEtaRange = true;
        }
      }
      if (cp0.emaxEC > 0) { // endcap
        if (std::abs(cp0.etaEC - dde->eta_raw()) < addCellsWindowEtaEndcap) {
          inEtaRange = true;
        }
      }
      if (!inEtaRange) {
        continue;
      }
 
      // Add L0/L1 when we are in the narrow range
      bool inPhiRange = false;
      if (cp0.emaxB > 0) { // barrel
        const double cell_phi = proxim(dde->phi_raw(), cp0.phiB);
        if (cell_phi > phiMinusB && cell_phi < phiPlusB) {
          inPhiRange = true;
        }
      }
      if (cp0.emaxEC > 0) { // endcap
        const double cell_phi = proxim(dde->phi_raw(), cp0.phiEC);
        if (cell_phi > phiMinusEC && cell_phi < phiPlusEC) {
          inPhiRange = true;
        }
      }
      if (!inPhiRange) {
        continue;
      }
 
      tofill.addCell(cell_itr.index(), cell_itr.weight());
    } // Cell Loop for L0/L1
  }   // Cluster loop for L0/L1
 
  if (!m_egammaCellRecoveryTool.empty()) {
    IegammaCellRecoveryTool::Info info;
    if (cp0.emaxB > cp0.emaxEC) {
      info.etamax = cp0.etaB;
      info.phimax = cp0.phiB;
    } else {
      info.etamax = cp0.etaEC;
      info.phimax = cp0.phiEC;
    }
    if (m_egammaCellRecoveryTool->execute(tofill,info).isFailure()) {
      ATH_MSG_WARNING("Issue trying to recover cells");
    }
 
    // And finally add the recovered cells
    const CaloCellContainer* inputcells =
      tofill.getCellLinks()->getCellContainer();
    for (const auto *c : info.addedCells) {
      int index = inputcells->findIndex(c->caloDDE()->calo_hash());
      tofill.addCell(index, 1.);
    }
  }
 
  return StatusCode::SUCCESS;
}
 
StatusCode
egammaSuperClusterBuilderBase::addTileGap3CellsinWindow(
  xAOD::CaloCluster& tofill,
  const CaloDetDescrManager& mgr) const
{
 
  double searchWindowEta = m_useExtendedTG3 ? 0.35 : 0.2;
  constexpr double searchWindowPhi = 2 * M_PI / 64.0 + M_PI / 64; // ~ 0.15 rad
  std::vector<const CaloCell*> cells;
  cells.reserve(16);
  const CaloCellContainer* inputcells =
    tofill.getCellLinks()->getCellContainer();
 
  if (!inputcells) {
    ATH_MSG_ERROR("No cell container in addRemainingCellsToCluster?");
    return StatusCode::FAILURE;
  }
 
  CaloCellList myList(&mgr, inputcells);
 
  const std::vector<CaloSampling::CaloSample> samples = {
    CaloSampling::TileGap3
  };
 
  for (auto samp : samples) {
    // quite slow
    myList.select(
      tofill.eta0(), tofill.phi0(), searchWindowEta, searchWindowPhi, samp);
    cells.insert(cells.end(), myList.begin(), myList.end());
  }
 
  for (const auto* cell : cells) {
    if (!cell) {
      continue;
    }
    const CaloDetDescrElement* dde = cell->caloDDE();
    if (!dde) {
      continue;
    }
 
    float maxEta = s_TG3Run2E4cellEtaMax;
    float minEta = s_TG3Run2E4cellEtaMin;
    if (m_useExtendedTG3) {
      minEta = s_TG3Run3E3cellEtaMin;
      // if |eta2| < 1.56, keep only E3, else keep E3+E4.
      // |eta2| uses as the eta of the highest energy cell in layer 2 as proxy
      if (std::abs(tofill.eta0()) > 1.56) {
        maxEta = s_TG3Run3E4cellEtaMax;
      }
    }
    float cellaEtaRaw = std::abs(dde->eta_raw());
    if (cellaEtaRaw >= minEta && cellaEtaRaw <= maxEta) {
      int index = inputcells->findIndex(dde->calo_hash());
      tofill.addCell(index, 1.);
    }
  }
  return StatusCode::SUCCESS;
}
 
StatusCode
egammaSuperClusterBuilderBase::calibrateCluster(
  const EventContext& ctx,
  xAOD::CaloCluster* newCluster,
  const CaloDetDescrManager& mgr,
  const xAOD::EgammaParameters::EgammaType egType,
  xAOD::CaloClusterContainer* precorrClusters) const
{
 
  refineEta1Position(newCluster, mgr);
  // Save the state before the corrections
  newCluster->setAltE(newCluster->e());
  newCluster->setAltEta(newCluster->eta());
  newCluster->setAltPhi(newCluster->phi());
  // first do the corrections
  if (precorrClusters) {
    precorrClusters->push_back(std::make_unique<xAOD::CaloCluster>());
    *precorrClusters->back() = *newCluster;
  }
  ATH_CHECK(m_clusterCorrectionTool->execute(
    ctx, newCluster, egType, xAOD::EgammaHelpers::isBarrel(newCluster)));
  double aeta = std::abs(newCluster->eta());
  if (aeta > 10) {
    ATH_MSG_DEBUG("Too large eta after S-shape corrections. "
          "SuperCluster rejected");
    return StatusCode::FAILURE;
  }
  newCluster->setRawE(newCluster->e());
  newCluster->setRawEta(newCluster->eta());
  newCluster->setRawPhi(newCluster->phi());
  //
  fillPositionsInCalo(newCluster, mgr);
  ATH_CHECK(m_MVACalibSvc->execute(*newCluster, egType));
 
  return StatusCode::SUCCESS;
}