From 1bc19e77e96e3f1ce2f1c661f5eb558dbc64cc07 Mon Sep 17 00:00:00 2001
From: Mark Hodgkinson <m.hodgkinson@sheffield.ac.uk>
Date: Wed, 16 Aug 2017 13:57:36 +0100
Subject: [PATCH] Add PFRecoverSplitShowersTool class and relevant entries in
 eflowRec_entries.cxx.

---
 .../eflowRec/PFRecoverSplitShowersTool.h      |  93 ++++++
 .../src/PFRecoverSplitShowersTool.cxx         | 292 ++++++++++++++++++
 .../src/components/eflowRec_entries.cxx       |   3 +
 3 files changed, 388 insertions(+)
 create mode 100644 Reconstruction/eflowRec/eflowRec/PFRecoverSplitShowersTool.h
 create mode 100644 Reconstruction/eflowRec/src/PFRecoverSplitShowersTool.cxx

diff --git a/Reconstruction/eflowRec/eflowRec/PFRecoverSplitShowersTool.h b/Reconstruction/eflowRec/eflowRec/PFRecoverSplitShowersTool.h
new file mode 100644
index 00000000000..bbc4af7f717
--- /dev/null
+++ b/Reconstruction/eflowRec/eflowRec/PFRecoverSplitShowersTool.h
@@ -0,0 +1,93 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef PFRECOVERSPLITSHOWERSTOOL_H
+#define PFRECOVERSPLITSHOWERSTOOL_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "xAODCaloEvent/CaloCluster.h"
+#include "xAODCaloEvent/CaloClusterContainer.h"
+#include "eflowRec/IPFSubtractionTool.h"
+
+#include <cassert>
+
+class eflowCaloObjectContainer;
+class eflowRecCluster;
+class eflowRecTrack;
+class eflowTrackCaloPoints;
+class eflowLayerIntegrator;
+class eflowEEtaBinnedParameters;
+class IEFlowCellEOverPTool;
+class PFTrackClusterMatchingTool;
+class eflowRingSubtractionManager;
+class eflowRecTrackContainer;
+class eflowRecClusterContainer;
+class eflowCaloObject;
+
+static const InterfaceID IID_PFRecoverSplitShowersTool("PFRecoverSplitShowersTool", 1, 0);
+
+class PFRecoverSplitShowersTool : virtual public IPFSubtractionTool, public AthAlgTool {
+  
+ public:
+
+  PFRecoverSplitShowersTool(const std::string& type,const std::string& name,const IInterface* parent);
+
+  ~PFRecoverSplitShowersTool();
+
+  static const InterfaceID& interfaceID();
+
+  virtual StatusCode initialize();
+  void execute(eflowCaloObjectContainer* theEflowCaloObjectContainer, eflowRecTrackContainer*, eflowRecClusterContainer*,xAOD::CaloClusterContainer& theCaloClusterContainer);
+  virtual StatusCode finalize();
+
+ private:
+
+  void recycleTracksClusters();
+  void getClustersToConsider();
+  void getTracksToRecover();
+  int matchAndCreateEflowCaloObj();
+  void performRecovery(int const nOriginalObj,xAOD::CaloClusterContainer& theCaloClusterContainer);
+  void subtractTrackFromClusters(const eflowTrackCaloPoints& trackCalo, eflowRingSubtractionManager& ranking, eflowRecTrack* efRecTrack, std::vector<xAOD::CaloCluster*> clusterSubtractionList);
+  double getSumEnergy(const std::vector<xAOD::CaloCluster*>& clusters);
+
+  void printClusterList(std::vector<xAOD::CaloCluster*>& clusters, std::string prefix);
+  void performSubtraction(eflowCaloObject* thisEflowCaloObject,xAOD::CaloClusterContainer& theCaloClusterContainer);
+
+private:
+
+  eflowCaloObjectContainer* m_eflowCaloObjectContainer;
+  std::vector<eflowRecCluster*> m_clustersToConsider;
+  std::vector<eflowRecTrack*> m_tracksToRecover;
+
+  double m_rCell;
+  double m_windowRms;
+
+  /* Tool for getting e/p values and hadronic shower cell ordering principle parameters */
+  ToolHandle<IEFlowCellEOverPTool> m_theEOverPTool;
+
+  /** Track-Cluster matching tool */
+  ToolHandle<PFTrackClusterMatchingTool> m_matchingTool;
+
+  std::unique_ptr<eflowEEtaBinnedParameters> m_binnedParameters;
+  std::unique_ptr<eflowLayerIntegrator> m_integrator;
+
+  double m_subtractionSigmaCut;
+
+  bool m_recoverIsolatedTracks;
+
+  /** Count the number of track-cluster matches -- for the summary in finalize */
+  unsigned int m_nTrackClusterMatches;
+
+  /** Toggle whether to use updated 2015 charged shower subtraction, which disables the shower subtraction in high calorimeter energy density regions  */
+  bool m_useUpdated2015ChargedShowerSubtraction;
+  
+};
+
+inline const InterfaceID& PFRecoverSplitShowersTool::interfaceID()
+{
+  return IID_PFRecoverSplitShowersTool;
+}
+
+#endif
diff --git a/Reconstruction/eflowRec/src/PFRecoverSplitShowersTool.cxx b/Reconstruction/eflowRec/src/PFRecoverSplitShowersTool.cxx
new file mode 100644
index 00000000000..da78251109d
--- /dev/null
+++ b/Reconstruction/eflowRec/src/PFRecoverSplitShowersTool.cxx
@@ -0,0 +1,292 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "eflowRec/PFRecoverSplitShowersTool.h"
+#include "eflowRec/eflowUtil.h"
+
+#include "eflowRec/eflowRecCluster.h"
+#include "eflowRec/eflowRecTrack.h"
+#include "eflowRec/eflowTrackClusterLink.h"
+#include "eflowRec/eflowCaloObject.h"
+#include "eflowRec/eflowCaloObjectMaker.h"
+#include "eflowRec/eflowTrackCaloPoints.h"
+#include "eflowRec/eflowCellList.h"
+#include "eflowRec/IEFlowCellEOverPTool.h"
+#include "eflowRec/PFTrackClusterMatchingTool.h"
+#include "eflowRec/eflowEEtaBinnedParameters.h"
+#include "eflowRec/eflowLayerIntegrator.h"
+#include "eflowRec/eflowRingSubtractionManager.h"
+#include "eflowRec/eflowCellSubtractionFacilitator.h"
+#include "eflowRec/eflowSubtractor.h"
+
+#include "CaloEvent/CaloClusterContainer.h"
+#include "xAODCaloEvent/CaloClusterKineHelper.h"
+
+using namespace eflowSubtract;
+
+PFRecoverSplitShowersTool::PFRecoverSplitShowersTool(const std::string& type,const std::string& name,const IInterface* parent):
+AthAlgTool(type, name, parent),
+m_eflowCaloObjectContainer(0),
+m_rCell(0.75),
+m_windowRms(0.032),
+m_theEOverPTool("eflowCellEOverPTool",this),
+m_matchingTool("PFTrackClusterMatchingTool/RcvrSpltMatchingTool", this),
+m_binnedParameters(std::make_unique<eflowEEtaBinnedParameters>()),
+m_integrator(std::make_unique<eflowLayerIntegrator>(m_windowRms, 1.0e-3, 3.0)),
+m_subtractionSigmaCut(1.5),
+m_recoverIsolatedTracks(false),
+m_nTrackClusterMatches(0),
+m_useUpdated2015ChargedShowerSubtraction(true)
+{
+  declareInterface<PFRecoverSplitShowersTool>(this);
+  declareProperty("SubtractionSigmaCut",m_subtractionSigmaCut);
+  declareProperty("eflowCellEOverPTool", m_theEOverPTool,"Energy Flow E/P Values and Shower Parameters Tool");
+  declareProperty("PFTrackClusterMatchingTool", m_matchingTool, "The track-cluster matching tool");
+  declareProperty("RecoverIsolatedTracks",m_recoverIsolatedTracks,"Whether to recover isolated tracks also");
+  declareProperty("useUpdated2015ChargedShowerSubtraction",m_useUpdated2015ChargedShowerSubtraction,"Toggle whether to use updated 2015 charged shower subtraction, which disables the shower subtraction in high calorimeter energy density region");
+  eflowRingSubtractionManager::setRMaxAndWeightRange(m_rCell, 1.0e6);
+}
+
+PFRecoverSplitShowersTool::~PFRecoverSplitShowersTool() {}
+
+StatusCode PFRecoverSplitShowersTool::initialize(){
+
+  // tool service
+  IToolSvc* myToolSvc;
+  if ( service("ToolSvc",myToolSvc).isFailure() ) {
+    msg(MSG::WARNING) << " Tool Service Not Found" << endmsg;
+    return StatusCode::SUCCESS;
+  }
+
+  if (m_matchingTool.retrieve().isFailure()){
+    msg(MSG::WARNING) << "Couldn't retrieve PFTrackClusterMatchingTool." << endmsg;
+    return StatusCode::SUCCESS;
+  }
+
+  if (m_theEOverPTool.retrieve().isFailure()){
+    msg(MSG::WARNING) << "Cannot find eflowEOverPTool" << endmsg;
+    return StatusCode::SUCCESS;
+  }
+
+  if (m_theEOverPTool->execute(m_binnedParameters.get()).isFailure()){
+    msg(MSG::WARNING) << "Could not execute eflowCellEOverPTool " << endmsg;
+    return StatusCode::SUCCESS;
+  }
+
+  return StatusCode::SUCCESS;
+}
+
+void PFRecoverSplitShowersTool::execute(eflowCaloObjectContainer* theEflowCaloObjectContainer, eflowRecTrackContainer*, eflowRecClusterContainer*,xAOD::CaloClusterContainer& theCaloClusterContainer){
+
+  msg(MSG::DEBUG) << "Executing PFRecoverSplitShowersTool" << endmsg;
+
+  m_eflowCaloObjectContainer = theEflowCaloObjectContainer;
+
+  recycleTracksClusters();
+
+  /* Add each track to its matching cluster's eflowCaloObject */
+  int const nOriginalObj = matchAndCreateEflowCaloObj();
+
+  /* Extrapolate tracks, match clusters, do shower simulation, run cell subtraction */
+  performRecovery(nOriginalObj,theCaloClusterContainer);
+
+}
+
+StatusCode PFRecoverSplitShowersTool::finalize(){
+
+  ATH_MSG_INFO("Produced " << m_nTrackClusterMatches << " track-cluster matches.");
+
+  return StatusCode::SUCCESS;
+
+}
+void PFRecoverSplitShowersTool::recycleTracksClusters() {
+  getTracksToRecover();
+  getClustersToConsider();
+
+}
+
+void PFRecoverSplitShowersTool::getClustersToConsider() {
+
+  m_clustersToConsider.clear();
+
+  eflowCaloObjectContainer::iterator  itEFCalObject = m_eflowCaloObjectContainer->begin();
+  eflowCaloObjectContainer::iterator endEFCalObject = m_eflowCaloObjectContainer->end();
+  for (; itEFCalObject != endEFCalObject; ++itEFCalObject) {
+    eflowCaloObject* thisEflowCaloObject = *itEFCalObject;
+
+    if (thisEflowCaloObject->nClusters() == 0) { continue; }
+
+    for(unsigned i=0; i<thisEflowCaloObject->nClusters(); ++i){
+        /* Skip empty clusters (subtraction remnants) */
+        const CaloClusterCellLink* theCellLink = thisEflowCaloObject->efRecCluster(i)->getCluster()->getCellLinks();
+        if (0 == (int)theCellLink->size()){ continue; }
+
+        thisEflowCaloObject->efRecCluster(i)->clearTrackMatches();
+        m_clustersToConsider.push_back(thisEflowCaloObject->efRecCluster(i));
+        thisEflowCaloObject->clearClusters();
+        thisEflowCaloObject->clearLinks();
+    }
+  }
+
+  std::sort(m_clustersToConsider.begin(),m_clustersToConsider.end(),eflowRecCluster::SortDescendingPt());
+}
+
+void PFRecoverSplitShowersTool::getTracksToRecover() {
+
+  m_tracksToRecover.clear(); int iEFCalOb=0;
+  eflowCaloObjectContainer::iterator  itEFCalObject = m_eflowCaloObjectContainer->begin();
+  eflowCaloObjectContainer::iterator endEFCalObject = m_eflowCaloObjectContainer->end();
+  for (; itEFCalObject != endEFCalObject; ++itEFCalObject, ++iEFCalOb) {
+    eflowCaloObject* thisEflowCaloObject = *itEFCalObject;
+
+    /* Skip isolated tracks if flag set */
+    if (!m_recoverIsolatedTracks && thisEflowCaloObject->nClusters() == 0) {
+      unsigned int nTrk = thisEflowCaloObject->nTracks();
+      /* But make sure we get eflowObjects from them
+       * TODO: replace this mechanism by something better */
+      for (unsigned int iTrk = 0; iTrk < nTrk; ++iTrk) {
+        thisEflowCaloObject->efRecTrack(iTrk)->setSubtracted();
+      }
+      continue;
+    }
+
+    /* Add all tracks on the eflowCaloObject that haven't been subtracted yet*/
+    std::vector<eflowRecTrack*> updatedTracks; updatedTracks.clear();
+    unsigned int nTracks = thisEflowCaloObject->nTracks();
+
+    /* For cluster-only CaloObjects */
+    if(nTracks == 0) continue;
+
+    /* For track-only and track-cluster CaloObjects */
+    for (unsigned int iTrack = 0; iTrack < nTracks; ++iTrack){
+      eflowRecTrack* thisEfRecTrack = thisEflowCaloObject->efRecTrack(iTrack);
+      if (thisEfRecTrack->isSubtracted()){
+        updatedTracks.push_back(thisEfRecTrack);
+        continue;
+      }
+      thisEfRecTrack->clearClusterMatches();
+      m_tracksToRecover.push_back(thisEfRecTrack);
+    }
+    thisEflowCaloObject->clearTracks();
+    if (!updatedTracks.empty()) {
+      thisEflowCaloObject->addTracks(updatedTracks);
+      updatedTracks.clear();
+    } else {
+      thisEflowCaloObject->clearLinks();
+    }
+  }
+
+  std::sort(m_tracksToRecover.begin(),m_tracksToRecover.end(),eflowRecTrack::SortDescendingPt());
+
+}
+
+int PFRecoverSplitShowersTool::matchAndCreateEflowCaloObj() {
+  /* Cache the original number of eflowCaloObjects */
+  const int nCaloObj = m_eflowCaloObjectContainer->size();
+
+  std::vector<eflowRecTrack*>::iterator endEfRecTrack = m_tracksToRecover.end();
+  /* loop tracks in m_tracksToRecover and do matching */
+  for (std::vector<eflowRecTrack*>::iterator itEfRecTrack = m_tracksToRecover.begin();
+      itEfRecTrack != endEfRecTrack; ++itEfRecTrack) {
+    eflowRecTrack* thisEfRecTrack = *itEfRecTrack;
+    /* No point to do anything if bin does not exist */
+    if (!thisEfRecTrack->hasBin()) {
+      std::unique_ptr<eflowCaloObject> thisEflowCaloObject = std::make_unique<eflowCaloObject>();
+      thisEflowCaloObject->addTrack(thisEfRecTrack);
+      m_eflowCaloObjectContainer->push_back(std::move(thisEflowCaloObject));
+      continue;
+    }
+    if (msgLvl(MSG::WARNING)){
+      const xAOD::TrackParticle* track = thisEfRecTrack->getTrack();
+      msg(MSG::DEBUG) << "Recovering charged EFO with e,eta and phi " << track->e() << ", "
+                << track->eta() << " and " << track->phi() << endmsg;
+    }
+    /* Get list of matched clusters */
+    std::vector<eflowRecCluster*> matchedClusters = m_matchingTool->doMatches(thisEfRecTrack, m_clustersToConsider, -1);
+    if (matchedClusters.empty()) { continue; }
+
+    m_nTrackClusterMatches += matchedClusters.size();
+    /* Matched cluster: create TrackClusterLink and add it to both the track and the cluster (eflowCaloObject will be created later) */
+    std::vector<eflowRecCluster*>::const_iterator itCluster = matchedClusters.begin();
+    std::vector<eflowRecCluster*>::const_iterator endCluster = matchedClusters.end();
+    for (; itCluster != endCluster; ++itCluster) {
+      eflowTrackClusterLink* trackClusterLink = eflowTrackClusterLink::getInstance(thisEfRecTrack,
+                                                                                   (*itCluster));
+      thisEfRecTrack->addClusterMatch(trackClusterLink);
+      (*itCluster)->addTrackMatch(trackClusterLink);
+    }
+  }
+
+  /* Create all eflowCaloObjects that contain eflowRecCluster */
+  eflowCaloObjectMaker makeCaloObject;
+  int nCaloObjects = makeCaloObject.makeTrkCluCaloObjects(m_tracksToRecover, m_clustersToConsider,
+                                                          m_eflowCaloObjectContainer);
+  msg(MSG::DEBUG) << "PFRecoverSplitShowersTool created " << nCaloObjects << " CaloObjects." << endmsg;
+
+  /* integrate cells; determine FLI; eoverp */
+  for (unsigned int iCalo = nCaloObj; iCalo < m_eflowCaloObjectContainer->size(); ++iCalo) {
+    eflowCaloObject* thisEflowCaloObject = m_eflowCaloObjectContainer->at(iCalo);
+    thisEflowCaloObject->simulateShower(m_integrator.get(), m_binnedParameters.get(), m_useUpdated2015ChargedShowerSubtraction);
+  }
+  return nCaloObj;
+}
+
+void PFRecoverSplitShowersTool::performSubtraction(eflowCaloObject* thisEflowCaloObject,xAOD::CaloClusterContainer& theCaloClusterContainer) {
+  for (unsigned iTrack = 0; iTrack < thisEflowCaloObject->nTracks(); ++iTrack) {
+    eflowRecTrack* thisEfRecTrack = thisEflowCaloObject->efRecTrack(iTrack);
+    /* Get matched cluster via Links */
+    std::vector<eflowRecCluster*> matchedClusters;
+    matchedClusters.clear();
+    std::vector<eflowTrackClusterLink*> links = thisEfRecTrack->getClusterMatches();
+    std::vector<eflowTrackClusterLink*>::iterator itLink = links.begin();
+    std::vector<eflowTrackClusterLink*>::iterator endLink = links.end();
+    for (; itLink != endLink; ++itLink) {
+      matchedClusters.push_back((*itLink)->getCluster());
+    }
+    /* Do subtraction */
+    std::vector<xAOD::CaloCluster*> clusterSubtractionList;
+    std::vector<eflowRecCluster*>::const_iterator itCluster = matchedClusters.begin();
+    std::vector<eflowRecCluster*>::const_iterator endCluster = matchedClusters.end();
+    for (; itCluster != endCluster; ++itCluster) {
+      clusterSubtractionList.push_back(
+          (*itCluster)->getClusterForModification(&theCaloClusterContainer));
+    }
+    if (getSumEnergy(clusterSubtractionList) - thisEfRecTrack->getEExpect() < m_subtractionSigmaCut
+        * sqrt(thisEfRecTrack->getVarEExpect())) {
+      /* Check if we can annihilate right away */
+      Subtractor::annihilateClusters(clusterSubtractionList);
+    } else {
+      /* Subtract the track from all matched clusters */
+      Subtractor::subtractTracksFromClusters(thisEfRecTrack, clusterSubtractionList);
+      /* Annihilate the cluster(s) if the remnant is small (i.e. below k*sigma) */
+      if (getSumEnergy(clusterSubtractionList) < m_subtractionSigmaCut
+          * sqrt(thisEfRecTrack->getVarEExpect())) {
+        Subtractor::annihilateClusters(clusterSubtractionList);
+      }
+    }
+    /* Flag tracks as subtracted */
+    thisEfRecTrack->setSubtracted();
+  }
+}
+
+void PFRecoverSplitShowersTool::performRecovery(int const nOriginalObj,xAOD::CaloClusterContainer& theCaloClusterContainer) {
+  unsigned int nEFCaloObs = m_eflowCaloObjectContainer->size();
+  for (unsigned int iCalo = nOriginalObj; iCalo < nEFCaloObs; ++iCalo) {
+    eflowCaloObject* thisEflowCaloObject = m_eflowCaloObjectContainer->at(iCalo);
+
+    performSubtraction(thisEflowCaloObject,theCaloClusterContainer);
+  }
+
+}
+
+
+double PFRecoverSplitShowersTool::getSumEnergy(const std::vector<xAOD::CaloCluster*>& clusters) {
+  double result = 0.0;
+  std::vector<xAOD::CaloCluster*>::const_iterator  itCluster = clusters.begin();
+  std::vector<xAOD::CaloCluster*>::const_iterator endCluster = clusters.end();
+  for (; itCluster != endCluster; ++itCluster) {
+    result += (*itCluster)->e();
+  }
+  return result;
+}
diff --git a/Reconstruction/eflowRec/src/components/eflowRec_entries.cxx b/Reconstruction/eflowRec/src/components/eflowRec_entries.cxx
index 51cceaf9e60..e69555956d0 100644
--- a/Reconstruction/eflowRec/src/components/eflowRec_entries.cxx
+++ b/Reconstruction/eflowRec/src/components/eflowRec_entries.cxx
@@ -20,6 +20,7 @@
 #include "eflowRec/PFClusterSelector.h"
 #include "eflowRec/PFAlgorithm.h"
 #include "eflowRec/PFCellLevelSubtractionTool.h"
+#include "eflowRec/PFRecoverSplitShowersTool.h"
 #include "GaudiKernel/DeclareFactoryEntries.h"
 
 DECLARE_ALGORITHM_FACTORY( eflowBuilder )
@@ -33,6 +34,7 @@ DECLARE_ALGORITHM_FACTORY( PFClusterSelector )
 DECLARE_ALGORITHM_FACTORY( PFTrackSelector )
 DECLARE_ALGORITHM_FACTORY( PFAlgorithm )
 DECLARE_TOOL_FACTORY( PFCellLevelSubtractionTool )
+DECLARE_TOOL_FACTORY( PFRecoverSplitShowersTool )
 DECLARE_TOOL_FACTORY( eflowRecoverSplitShowersTool )
 DECLARE_TOOL_FACTORY( eflowCellLevelSubtractionTool )
 DECLARE_TOOL_FACTORY( eflowLCCalibTool )
@@ -57,6 +59,7 @@ DECLARE_FACTORY_ENTRIES(eflowRec) {
     DECLARE_ALGORITHM( PFTrackSelector )
     DECLARE_ALGORITHM( PFAlgorithm )
     DECLARE_TOOL( PFCellLevelSubtractionTool )
+    DECLARE_TOOL( PFRecoverSplitShowersTool )
     DECLARE_TOOL ( eflowRecoverSplitShowersTool )
     DECLARE_TOOL ( eflowCellLevelSubtractionTool )
     DECLARE_TOOL ( eflowMomentCalculatorTool )
-- 
GitLab