Merge branch 'forward_electron_cells' into 'master'

ATLASRECTS-5709 : Forward electron cluster cells

Closes ATLASRECTS-5709

See merge request atlas/athena!39203

Reconstruction/RecExample/RecExCommon/share/RecExCommon_topOptions.py

@@ -1285,6 +1285,8 @@ if ( rec.doAOD() or rec.doWriteAOD()) and not rec.readAOD() :
             if ( rec.readESD() or jobproperties.egammaRecFlags.Enabled ) and not rec.ScopingLevel()==4 and rec.doEgamma :
                 from egammaRec import egammaKeys
                 addClusterToCaloCellAOD(egammaKeys.outputClusterKey())
+                addClusterToCaloCellAOD(egammaKeys.outputFwdClusterKey())
+                addClusterToCaloCellAOD(egammaKeys.outputEgammaLargeFWDClustersKey())
                 if "itemList" in metadata:
                     if ('xAOD::CaloClusterContainer', egammaKeys.EgammaLargeClustersKey()) in metadata["itemList"]:
                         # check first for priority if both keys are in metadata

Reconstruction/egamma/egammaAlgs/python/egammaLargeFWDClusterMakerAlg.py

+# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+
+__doc__ = """ToolFactory to instantiate egammaLargeClusterMaker
+with FWD configuration"""
+__author__ = "Jovan Mitrevski"
+
+from egammaTools.egammaToolsFactories import egammaLargeFWDClusterMakerTool
+from CaloRec import CaloRecConf
+from egammaRec.Factories import AlgFactory, FcnWrapper
+from egammaRec import egammaKeys
+from CaloClusterCorrection.CaloSwCorrections import make_CaloSwCorrections
+
+
+def clusMakerTools():
+    return [egammaLargeFWDClusterMakerTool()]
+
+egammaLargeFWDClusterMakerAlg = AlgFactory(
+    CaloRecConf.CaloClusterMaker,
+    name="egammaLargeFWDClusterMaker",
+    SaveUncalibratedSignalState=False,
+    ClustersOutputName=egammaKeys.EgammaLargeFWDClustersKey(),
+    ClusterMakerTools=FcnWrapper(clusMakerTools),
+    ClusterCorrectionTools=make_CaloSwCorrections("FWDele6_6",
+                                                  suffix="Nocorr",
+                                                  version="none",
+                                                  corrlist=[],
+                                                  cells_name=egammaKeys.caloCellKey()))

Reconstruction/egamma/egammaAlgs/python/egammaLargeFWDClusterMakerAlgConfig.py

+# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+
+__doc__ = "Configure egammaLargeFWDClusterMaker, which chooses cells to store in the AOD" 
+__author__ = "Jovan Mitrevski"
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+from egammaTools.egammaLargeFWDClusterMakerConfig import egammaLargeFWDClusterMakerCfg
+from CaloClusterCorrection.CaloSwCorrections import make_CaloSwCorrections
+CaloClusterMaker=CompFactory.CaloClusterMaker
+
+def egammaLargeFWDClusterMakerAlgCfg(flags, name = "egammaLargeClusterMaker", **kwargs):
+
+    acc = ComponentAccumulator
+    
+    kwargs.setdefault("SaveUncalibratedSignalState", False)
+    kwargs.setdefault("ClustersOutputName", flags.Egamma.Keys.Output.EgammaLargeFWDClusters)
+
+    if "ClusterMakerTools" not in kwargs:
+        toolAcc = egammaLargeFWDClusterMakerCfg(flags)
+        kwargs["ClusterMakerTools"] = [ toolAcc.popPrivateTools() ]
+        acc.merge(toolAcc)
+
+    kwargs.setdefault("ClusterCorrectionTools", make_CaloSwCorrections("FWDele6_6",
+                                                                       suffix="Nocorr",
+                                                                       version="none",
+                                                                       cells_name=flags.Egamma.Keys.Input.CaloCells))
+
+    acc.addEventAlgo(CaloClusterMaker(name, **kwargs))
+    return acc
+

Reconstruction/egamma/egammaConfig/python/egammaConfigFlags.py

@@ -92,7 +92,12 @@ def createEgammaConfigFlags():
     egcf.addFlag("Egamma.Keys.Output.ForwardClusters", 'ForwardElectronClusters')
     egcf.addFlag("Egamma.Keys.Output.ForwardClustersSuppESD", '-SisterCluster')
     egcf.addFlag("Egamma.Keys.Output.ForwardClustersSuppAOD",
-                 '-SisterCluster.-CellLink')
+                 '-SisterCluster')
+
+    # These are the clusters that are used to determine which cells to write out to AOD
+    egcf.addFlag("Egamma.Keys.Output.EgammaLargeFWDClusters", 'egamma66FWDClusters')
+    egcf.addFlag("Egamma.Keys.Output.EgammaLargeFWDClustersSuppESD", '')
+    # don't define SuppAOD because the whole container is suppressed
 
     egcf.addFlag("Egamma.Keys.Output.Photons", 'Photons')
     egcf.addFlag("Egamma.Keys.Output.PhotonsSuppESD", '')

Reconstruction/egamma/egammaRec/python/egammaForwardGetter.py

@@ -47,11 +47,11 @@ egammaForwardBuilder = AlgFactory(
 
 
 class egammaForwardGetter (Configured):
-
+    
     def configure(self):
         mlog = logging.getLogger('egammaForwardGetter.py::configure:')
         mlog.info('entering')
-
+        
         # configure egammaForward here:
         try:
             self._egammaFwdBuilderHandle = egammaForwardBuilder()
@@ -59,7 +59,21 @@ class egammaForwardGetter (Configured):
             mlog.error("could not get handle to egammaForward")
             traceback.print_exc()
             return False
+            
+            # the egammaLargeFWDClusterMaker
+            # (Which chooses the cells to store in the AOD)
+        from egammaAlgs.egammaLargeFWDClusterMakerAlg import (
+            egammaLargeFWDClusterMakerAlg)
+        try:
+            self._egammaLargeClusterMaker = egammaLargeFWDClusterMakerAlg()
+        except Exception:
+            mlog.error("could not get handle to egammaLargeClusterMaker")
+            import traceback
+            traceback.print_exc()
+            return False
+                
         return True
-
+                
     def egammaFwdBuilderHandle(self):
         return self._egammaFwdBuilderHandle
+                    

Reconstruction/egamma/egammaRec/python/egammaKeys.py

@@ -23,6 +23,8 @@ class egammaKeysDict:
         Cluster=['xAOD::CaloClusterContainer', 'egammaClusters', '', ''],
         EgammaLargeClusters=['xAOD::CaloClusterContainer',
                              'egamma711Clusters', '', ''],  # not output to AOD
+        EgammaLargeFWDClusters=['xAOD::CaloClusterContainer',
+                                'egamma66FWDClusters', '', ''],  # not output to AOD
         TopoSeededCluster=['xAOD::CaloClusterContainer',
                            'egammaTopoSeededClusters', '', '-CellLink'],
         Electron=['xAOD::ElectronContainer', 'Electrons',
@@ -35,7 +37,7 @@ class egammaKeysDict:
         FwdElectron=['xAOD::ElectronContainer',
                      'ForwardElectrons', '', FwdElectronisemSupress],
         FwdCluster=['xAOD::CaloClusterContainer',
-                    'ForwardElectronClusters', '-SisterCluster', '.-CellLink'],
+                    'ForwardElectronClusters', '-SisterCluster', ''],
         Photon=['xAOD::PhotonContainer', 'Photons', '',
                 ShowerShapesSuppress+PhotonisemSupress],
         TrackParticle=['xAOD::TrackParticleContainer', 'GSFTrackParticles',
@@ -58,6 +60,10 @@ class egammaKeysDict:
                                               outputs['EgammaLargeClusters'][1] +
                                               '_links',
                                               '', '']
+    outputs['EgammaLargeFWDClustersCellLink'] = ['CaloClusterCellLinkContainer',
+                                              outputs['EgammaLargeFWDClusters'][1] +
+                                              '_links',
+                                              '', '']
     #
 
 

Reconstruction/egamma/egammaRec/share/egammaOutputItemList_jobOptions.py

@@ -44,14 +44,13 @@ def addAuxContainer(outputList, cType, cKey, auxOptionAll='', auxOptionAOD=''):
 AOD_outputs = [i for i, j in egammaKeysDict.outputs.items()
                if i not in ('EgammaRec', 'PhotonSuperRec',
                             'ElectronSuperRec', 'TopoSeededCellLink',
-                            'FwdClusterCellLink', 'EgammaLargeClusters',
-                            'EgammaLargeClustersCellLink')]
+                            'EgammaLargeClusters',
+                            'EgammaLargeClustersCellLink','EgammaLargeFWDClusters', 'EgammaLargeFWDClustersCellLink')]
 
 
 ESD_outputs = [i for i, j in egammaKeysDict.outputs.items()
                if i not in ('EgammaRec', 'PhotonSuperRec',
-                            'ElectronSuperRec', 'TopoSeededCellLink',
-                            'FwdClusterCellLink')]
+                            'ElectronSuperRec', 'TopoSeededCellLink')]
 
 # Define egammaAODList in the proper format (<type>#<key><option>),
 # including aux containers

Reconstruction/egamma/egammaTools/python/egammaLargeFWDClusterMakerConfig.py

+# Copyright (C) 2002-2019 CERN for the benefit of the ATLAS collaboration
+
+__doc__ = "Configure egammaLargeFWDClusterMaker, which chooses cells to store in the AOD" 
+__author__ = "Jovan Mitrevski"
+
+from AthenaConfiguration.ComponentAccumulator import ComponentAccumulator
+from AthenaConfiguration.ComponentFactory import CompFactory
+egammaLargeFWDClusterMaker=CompFactory.egammaLargeFWDClusterMaker
+
+def egammaLargeFWDClusterMakerCfg(flags, **kwargs):
+
+    acc = ComponentAccumulator
+    
+    kwargs.setdefault("CellsName", flags.Egamma.Keys.Input.CaloCells)
+    kwargs.setdefault("InputClusterCollection", flags.Egamma.Keys.Output.FwdCluster)
+    kwargs.setdefault("doFWDelesurraundingWindows", True)
+
+    acc.setPrivateTools(egammaLargeFWDClusterMaker(**kwargs))
+
+    return acc
+

Reconstruction/egamma/egammaTools/python/egammaToolsFactories.py

@@ -124,6 +124,14 @@ egammaLargeClusterMakerTool = ToolFactory(
     CellsName=egammaKeys.caloCellKey()
 )
 
+egammaLargeFWDClusterMakerTool = ToolFactory(
+    egammaToolsConf.egammaLargeClusterMaker,
+    name="egammaLCFWDMakerTool",
+    InputClusterCollection=egammaKeys.FwdClusterKey(),
+    CellsName=egammaKeys.caloCellKey(),
+    doFWDelesurraundingWindows = True
+)
+
 # Electron Selectors
 ElectronPIDBuilder = ToolFactory(
     EMPIDBuilderElectronBase,

Reconstruction/egamma/egammaTools/src/egammaLargeClusterMaker.cxx

@@ -5,6 +5,7 @@
 #include "egammaLargeClusterMaker.h"
 #include "CaloEvent/CaloCellContainer.h"
 #include "CaloEvent/CaloClusterCellLink.h"
+#include "CaloUtils/CaloCellList.h"
 #include "CaloUtils/CaloClusterStoreHelper.h"
 #include "egammaUtils/egammaEnergyPositionAllSamples.h"
 #include "xAODCaloEvent/CaloCluster.h"
@@ -52,83 +53,186 @@ egammaLargeClusterMaker::execute(const EventContext& ctx,
 
   const CaloDetDescrManager* dd_man = nullptr;
   ATH_CHECK(detStore()->retrieve(dd_man, "CaloMgr"));
-
+  
   // The main loop over clusters
   for (const auto *cluster : *inputClusters) {
-
-    // find the center of the cluster, copying the logic of
-    // egammaMiddleShape.cxx
-
-    // check if cluster is in barrel or in the end-cap
-    if (!cluster->inBarrel() && !cluster->inEndcap()) {
-      ATH_MSG_DEBUG(" Cluster is neither in Barrel nor in Endcap; skipping ");
-      continue;
-    }
-
-    // check if cluster is in barrel or end-cap
-    bool in_barrel = egammaEnergyPositionAllSamples::inBarrel(*cluster, 2);
-    CaloSampling::CaloSample sam = CaloSampling::EMB2;
-    if (in_barrel) {
-      sam = CaloSampling::EMB2;
-    } else {
-      sam = CaloSampling::EME2;
-    }
-    // granularity in (eta,phi) in the pre sampler
-    // CaloCellList needs both enums: subCalo and CaloSample
-    auto eta = cluster->etaSample(sam);
-    auto phi = cluster->phiSample(sam);
-
-    if ((eta == 0. && phi == 0.) || fabs(eta) > 100) {
-      return StatusCode::SUCCESS;
-    }
-
-    // Should get overritten
-    bool barrel = false;
-    CaloCell_ID::SUBCALO subcalo = CaloCell_ID::LAREM;
-    int sampling_or_module = 0;
-
-    CaloDetDescrManager::decode_sample(
+    
+    if (!m_isFWD) {
+      // find the center of the cluster, copying the logic of
+      // egammaMiddleShape.cxx
+      
+      // check if cluster is in barrel or in the end-cap
+      if (!cluster->inBarrel() && !cluster->inEndcap()) {
+	ATH_MSG_DEBUG(" Cluster is neither in Barrel nor in Endcap; skipping ");
+	continue;
+      }
+      
+      // check if cluster is in barrel or end-cap
+      bool in_barrel = egammaEnergyPositionAllSamples::inBarrel(*cluster, 2);
+      CaloSampling::CaloSample sam = CaloSampling::EMB2;
+      if (in_barrel) {
+	sam = CaloSampling::EMB2;
+      } else {
+	sam = CaloSampling::EME2;
+      }
+      
+      // granularity in (eta,phi) in the pre sampler
+      // CaloCellList needs both enums: subCalo and CaloSample
+      auto eta = cluster->etaSample(sam);
+      auto phi = cluster->phiSample(sam);
+      
+      if ((eta == 0. && phi == 0.) || fabs(eta) > 100) {
+	return StatusCode::SUCCESS;
+      }
+      
+      // Should get overritten
+      bool barrel = false;
+      CaloCell_ID::SUBCALO subcalo = CaloCell_ID::LAREM;
+      int sampling_or_module = 0;
+      
+      CaloDetDescrManager::decode_sample(
       subcalo, barrel, sampling_or_module, (CaloCell_ID::CaloSample)sam);
-
-    // Get the corresponding grannularities : needs to know where you are
-    //                  the easiest is to look for the CaloDetDescrElement
-
-    const CaloDetDescrElement* dde =
-      dd_man->get_element(CaloCell_ID::LAREM, 0, barrel, eta, phi);
-
-    // if object does not exist then return
-    if (!dde) {
-      return StatusCode::SUCCESS;
-    }
-
-    // local granularity
-    auto deta = dde->deta();
-    auto dphi = dde->dphi();
-
-    // search the hottest cell around the (eta,phi)
-    // (eta,phi) are defined as etaSample() and phiSample
-    // around this position a hot cell is searched for in a window
-    // (m_neta*m_deta,m_nphi*m_dphi), by default (m_neta,m_nphi)=(7,7)
-    CaloLayerCalculator calc;
-    StatusCode sc = calc.fill(*dd_man,
-                              cellcoll.ptr(),
-                              cluster->etaSample(sam),
-                              cluster->phiSample(sam),
-                              m_neta * deta,
-                              m_nphi * dphi,
-                              (CaloSampling::CaloSample)sam);
-    if (sc.isFailure()) {
-      ATH_MSG_WARNING("CaloLayerCalculator failed  fill ");
-    }
-    double etamax = calc.etarmax();
-    double phimax = calc.phirmax();
-
-    // create the new cluster
-    xAOD::CaloCluster* newCluster =
-      CaloClusterStoreHelper::makeCluster(collection, cellcoll.ptr());
-    newCluster->setEta0(etamax);
-    newCluster->setPhi0(phimax);
-    newCluster->setClusterSize(xAOD::CaloCluster::SW_7_11);
-  }
+      
+      // Get the corresponding grannularities : needs to know where you are
+      //                  the easiest is to look for the CaloDetDescrElement
+      
+      const CaloDetDescrElement* dde =
+	dd_man->get_element(CaloCell_ID::LAREM, 0, barrel, eta, phi);
+      
+      // if object does not exist then return
+      if (!dde) {
+	return StatusCode::SUCCESS;
+      }
+      
+      // local granularity
+      auto deta = dde->deta();
+      auto dphi = dde->dphi();
+      
+      // search the hottest cell around the (eta,phi)
+      // (eta,phi) are defined as etaSample() and phiSample
+      // around this position a hot cell is searched for in a window
+      // (m_neta*m_deta,m_nphi*m_dphi), by default (m_neta,m_nphi)=(7,7)
+      CaloLayerCalculator calc;
+      StatusCode sc = calc.fill(*dd_man,
+				cellcoll.ptr(),
+				cluster->etaSample(sam),
+				cluster->phiSample(sam),
+				m_neta * deta,
+				m_nphi * dphi,
+				(CaloSampling::CaloSample)sam);
+      if (sc.isFailure()) {
+	ATH_MSG_WARNING("CaloLayerCalculator failed  fill ");
+      }
+      double etamax = calc.etarmax();
+      double phimax = calc.phirmax();
+      
+      // create the new cluster
+      xAOD::CaloCluster* newCluster =
+	CaloClusterStoreHelper::makeCluster(collection, cellcoll.ptr());
+      newCluster->setEta0(etamax);
+      newCluster->setPhi0(phimax);
+      newCluster->setClusterSize(xAOD::CaloCluster::SW_7_11);
+      
+    } else { 
+      // FWD cluster collection
+      
+      // check if cluster is in EMEC or FCAL 
+      if (cluster->inBarrel()) {
+	ATH_MSG_DEBUG(" Cluster is not in FWD; skip ");
+	continue;
+      } else {
+	ATH_MSG_DEBUG (" CLuster is FWD Cluster ");
+      }
+      
+      // check if cluster is in FCAL or EMEC
+      bool in_EMEC = false;
+      if (cluster->eSample(CaloSampling::EME2) > cluster->eSample(CaloSampling::FCAL0) )
+	in_EMEC = true;
+      
+      CaloSampling::CaloSample sam = CaloSampling::FCAL0;
+      if (in_EMEC) {
+	sam = CaloSampling::EME2;
+      } else {
+	sam = CaloSampling::FCAL0;
+      }
+      
+      // granularity in (eta,phi) in the pre sampler
+      // CaloCellList needs both enums: subCalo and CaloSample
+      auto eta = cluster->etaSample(sam);
+      auto phi = cluster->phiSample(sam);
+       
+      if ((eta == 0. && phi == 0.) || fabs(eta) > 100) {
+	return StatusCode::SUCCESS;
+      }
+      
+      // Should get overritten
+      bool emec = false;
+      CaloCell_ID::SUBCALO subcalo = CaloCell_ID::LARFCAL;
+      int sampling_or_module = 0;
+      
+      CaloDetDescrManager::decode_sample(subcalo, emec, sampling_or_module, (CaloCell_ID::CaloSample)sam);
+      
+      // Get the corresponding grannularities : needs to know where you are
+      //                  the easiest is to look for the CaloDetDescrElement
+      //    const CaloDetDescrElement* get_element_FCAL (const CaloDetDescriptor* reg, double eta, double phi) const;
+      const CaloDetDescrElement* dde =
+	dd_man->get_element(subcalo, sampling_or_module, emec, eta, phi);
+      
+      // if object does not exist then return
+      if (!dde) {
+	return StatusCode::SUCCESS;
+      }
+      
+      // local granularity
+      auto deta = dde->deta();
+      auto dphi = dde->dphi();
+      
+      // search the hottest cell around the (eta,phi)
+      // (eta,phi) are defined as etaSample() and phiSample
+      // around this position a hot cell is searched for in a window
+      // (m_neta*m_deta,m_nphi*m_dphi), by default (m_neta,m_nphi)=(6,6) 
+      // for EMEC-OUTER FWD much bigger cell size 
+      
+      // create the new cluster
+      xAOD::CaloCluster* newCluster =
+	CaloClusterStoreHelper::makeCluster(collection, cellcoll.ptr());
+      
+      // if In EMEC
+      CaloLayerCalculator calc;
+      StatusCode sc = calc.fill(*dd_man,
+				cellcoll.ptr(),
+				cluster->etaSample(sam),
+				cluster->phiSample(sam),
+				m_netaFWD * deta,
+				m_nphiFWD * dphi,
+				(CaloSampling::CaloSample)sam, in_EMEC ? newCluster : nullptr);
+      
+      if (sc.isFailure()) {
+	ATH_MSG_WARNING("CaloLayerCalculator failed  fill for FWD cluster");
+      }
+      double etamax = calc.etarmax();
+      double phimax = calc.phirmax();
+      
+      newCluster->setEta0(etamax);
+      newCluster->setPhi0(phimax);
+      
+      if (!in_EMEC) {
+	// If FCAL need to add cell to cluster in a cone
+	std::vector<const CaloCell*> cells;
+	cells.reserve(300);
+	CaloCellList myList(cellcoll.ptr());
+	myList.select(cluster->etaSample(sam), cluster->phiSample(sam), m_drFWD,(CaloSampling::CaloSample)sam);
+	//  myList.select(dde,newCluster->eta0(), newCluster->phi0(),m_drFWD,(CaloSampling::CaloSample)sam); 
+	cells.insert(cells.end(), myList.begin(), myList.end());
+	
+	for ( const auto *cell : cells ) {
+	  if( !cell || !cell->caloDDE() ) continue;
+	  int index = cellcoll.ptr()->findIndex(cell->caloDDE()->calo_hash());
+	  if( index == -1 ) continue;
+	  newCluster->addCell(index,1.);
+	}
+      }
+    }// end isFWD
+  }// main loop over clusters
   return StatusCode::SUCCESS;
 }

Reconstruction/egamma/egammaTools/src/egammaLargeClusterMaker.h

@@ -51,6 +51,12 @@ private:
   /** @brief Cell container*/
   SG::ReadHandleKey<CaloCellContainer> m_cellsKey {this,
       "CellsName", "AllCalo", "Names of containers which contain cells"};
+
+  /** @brief do FWD cell **/
+  Gaudi::Property<bool> m_isFWD{ this,
+      "doFWDelesurraundingWindows",
+      false,
+      "Save FWD electron surraunding windows" };
   
   Gaudi::Property<double> m_neta {this, "Neta", 7.0,
       "Number of eta cells in each sampling in which to look for hottest cell"};
@@ -58,6 +64,15 @@ private:
   Gaudi::Property<double> m_nphi {this, "Nphi", 7.0,
       "Number of phi cell in each sampling in which to look for hottest cell"};
 
+  Gaudi::Property<double> m_netaFWD {this, "NetaFWD", 6.0,
+      "Number of eta cells in each sampling in which to look for hottest cell"};
+
+  Gaudi::Property<double> m_nphiFWD {this, "NphiFWD", 6.0,
+      "Number of phi cell in each sampling in which to look for hottest cell"};
+
+  Gaudi::Property<double> m_drFWD {this, "deltaR_FCAL", 0.4,
+      "Cone size to collec cell around hottest-cell FCAL"};
+
 
 };