From 10dea0a6004e855538be0a25c836aa7c44705f7f Mon Sep 17 00:00:00 2001
From: Mark Hodgkinson <m.hodgkinson@sheffield.ac.uk>
Date: Thu, 10 Aug 2017 14:45:00 +0100
Subject: [PATCH] Add first version of PFAlgorithm and new cell level
subtraction tool. Include relevant updates in
eflowRec/src/eflowRec_entries.cxx
Former-commit-id: 4cf052947d2ea37824398519b7555be8925b6f48
---
.../eflowRec/eflowRec/PFAlgorithm.h | 42 ++
.../eflowRec/PFCellLevelSubtractionTool.h | 97 ++++
.../eflowRec/eflowRec/PFISubtractionTool.h | 24 +
Reconstruction/eflowRec/src/PFAlgorithm.cxx | 60 +++
.../src/PFCellLevelSubtractionTool.cxx | 465 ++++++++++++++++++
.../src/components/eflowRec_entries.cxx | 6 +
6 files changed, 694 insertions(+)
create mode 100644 Reconstruction/eflowRec/eflowRec/PFAlgorithm.h
create mode 100644 Reconstruction/eflowRec/eflowRec/PFCellLevelSubtractionTool.h
create mode 100644 Reconstruction/eflowRec/eflowRec/PFISubtractionTool.h
create mode 100644 Reconstruction/eflowRec/src/PFAlgorithm.cxx
create mode 100644 Reconstruction/eflowRec/src/PFCellLevelSubtractionTool.cxx
diff --git a/Reconstruction/eflowRec/eflowRec/PFAlgorithm.h b/Reconstruction/eflowRec/eflowRec/PFAlgorithm.h
new file mode 100644
index 00000000000..b72d4700bf9
--- /dev/null
+++ b/Reconstruction/eflowRec/eflowRec/PFAlgorithm.h
@@ -0,0 +1,42 @@
+#ifndef PFALGORITHM_H
+#define PFALGORITHM_H
+
+#include "AthenaBaseComps/AthAlgorithm.h"
+#include "GaudiKernel/ToolHandle.h"
+#include "StoreGate/DataHandle.h"
+
+#include "eflowRec/eflowRecCluster.h"
+#include "eflowRec/eflowRecTrack.h"
+
+#include "xAODCaloEvent/CaloClusterContainer.h"
+
+class PFISubtractionTool;
+
+class PFAlgorithm : public AthAlgorithm {
+
+public:
+ PFAlgorithm(const std::string& name, ISvcLocator* pSvcLocator);
+ ~PFAlgorithm() {};
+
+ StatusCode initialize();
+ StatusCode execute();
+ StatusCode finalize();
+
+private:
+ /** List of PFISubtractionTool, which will be executed by this algorithm */
+ ToolHandleArray<PFISubtractionTool> m_tools;
+
+ /** ReadHandle for the eflowRecTrackContainer to be read in */
+ SG::ReadHandle<eflowRecTrackContainer> m_eflowRecTracksReadHandle;
+
+ /** ReadHandle for the eflowRecClusterContainer to be read in */
+ SG::ReadHandle<eflowRecClusterContainer> m_eflowRecClustersReadHandle;
+
+ /** WriteHandle for CaloClusterContainer to be written out */
+ SG::WriteHandle<xAOD::CaloClusterContainer> m_caloClustersWriteHandle;
+
+ /** Funciton to print out list of tools if in VERBOSE mode */
+ void printTools();
+
+};
+#endif
diff --git a/Reconstruction/eflowRec/eflowRec/PFCellLevelSubtractionTool.h b/Reconstruction/eflowRec/eflowRec/PFCellLevelSubtractionTool.h
new file mode 100644
index 00000000000..960fd5a93ef
--- /dev/null
+++ b/Reconstruction/eflowRec/eflowRec/PFCellLevelSubtractionTool.h
@@ -0,0 +1,97 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef PFCELLLEVELSUBTRACTIONTOOL_H
+#define PFCELLLEVELSUBTRACTIONTOOL_H
+
+#include "AthenaBaseComps/AthAlgTool.h"
+#include "eflowRec/PFISubtractionTool.h"
+#include "GaudiKernel/ToolHandle.h"
+
+#include "eflowRec/eflowCellList.h"
+#include "eflowRec/eflowLayerIntegrator.h"
+#include "eflowRec/eflowEEtaBinnedParameters.h"
+#include "xAODCaloEvent/CaloCluster.h"
+#include "xAODTracking/TrackParticle.h"
+
+#include <vector>
+#include <cassert>
+
+class eflowCaloObjectContainer;
+class eflowRecTrackContainer;
+class eflowRecClusterContainer;
+class IEFlowCellEOverPTool;
+class PFTrackClusterMatchingTool;
+class eflowRecTrack;
+
+static const InterfaceID IID_PFCellLevelSubtractionTool("PFCellLevelSubtractionTool", 1, 0);
+
+class PFCellLevelSubtractionTool : virtual public PFISubtractionTool, public AthAlgTool {
+public:
+
+ PFCellLevelSubtractionTool(const std::string& type,const std::string& name,const IInterface* parent);
+
+ ~PFCellLevelSubtractionTool();
+
+ static const InterfaceID& interfaceID();
+
+ StatusCode initialize();
+ void execute(eflowCaloObjectContainer* theEflowCaloObjectContainer, eflowRecTrackContainer* recTrackContainer, eflowRecClusterContainer* recClusterContainer);
+ StatusCode finalize();
+
+ private:
+
+ void calculateRadialEnergyProfiles();
+ void calculateAverageEnergyDensity();
+ void performSubtraction();
+ bool runInGoldenMode() { return ((m_goldenModeString == "golden1") || (m_goldenModeString == "golden2")); }
+ bool isEOverPFail(double expectedEnergy, double sigma, double clusterEnergy, bool consistencySigmaCut, bool useGoldenMode);
+ bool canAnnihilated(double expectedEnergy, double sigma, double clusterEnergy);
+
+ int matchAndCreateEflowCaloObj(int n);
+ std::string printTrack(const xAOD::TrackParticle* track);
+ std::string printCluster(const xAOD::CaloCluster* cluster);
+ void printAllClusters(const eflowRecClusterContainer& recClusterContainer);
+
+ private:
+
+ eflowCaloObjectContainer* m_eflowCaloObjectContainer;
+ eflowRecTrackContainer* m_eflowTrackContainer;
+ eflowRecClusterContainer* m_eflowClusterContainer;
+
+ ToolHandle<PFTrackClusterMatchingTool> m_matchingTool;
+ /* Matching tools for calculating the pull */
+ ToolHandle<PFTrackClusterMatchingTool> m_matchingToolForPull_015;
+ ToolHandle<PFTrackClusterMatchingTool> m_matchingToolForPull_02;
+
+ /* Tools for "shower simulation" */
+ std::unique_ptr<eflowEEtaBinnedParameters> m_binnedParameters;
+ std::unique_ptr<eflowLayerIntegrator> m_integrator;
+ ToolHandle<IEFlowCellEOverPTool> m_theEOverPTool;
+
+ //double m_rCell;
+
+ double m_subtractionSigmaCut;
+ double m_consistencySigmaCut;
+
+ //flag if we want to caluclate EoVerP from single particles we dont do eflow, so clusters are original unmodified ones.
+ bool m_calcEOverP;
+
+ // string flag to configure for running in golden e/p match mode
+ std::string m_goldenModeString;
+
+ // Number of clusters to match
+ int m_nMatchesInCellLevelSubtraction;
+
+ /** 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& PFCellLevelSubtractionTool::interfaceID()
+{
+ return IID_PFCellLevelSubtractionTool;
+}
+
+#endif
diff --git a/Reconstruction/eflowRec/eflowRec/PFISubtractionTool.h b/Reconstruction/eflowRec/eflowRec/PFISubtractionTool.h
new file mode 100644
index 00000000000..2a5d55f6a06
--- /dev/null
+++ b/Reconstruction/eflowRec/eflowRec/PFISubtractionTool.h
@@ -0,0 +1,24 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef PFISUBTRACTIONALGTOOL_H
+#define PFISUBTRACTIONALGTOOL_H
+
+#include "GaudiKernel/IAlgTool.h"
+
+class eflowCaloObjectContainer;
+class eflowRecTrackContainer;
+class eflowRecClusterContainer;
+
+class PFISubtractionTool : virtual public IAlgTool {
+
+ public:
+
+ virtual StatusCode intialize() {return StatusCode::SUCCESS;}
+ virtual void execute(eflowCaloObjectContainer*, eflowRecTrackContainer*, eflowRecClusterContainer*) = 0;
+ virtual StatusCode finalize() {return StatusCode::SUCCESS;}
+
+};
+
+#endif
diff --git a/Reconstruction/eflowRec/src/PFAlgorithm.cxx b/Reconstruction/eflowRec/src/PFAlgorithm.cxx
new file mode 100644
index 00000000000..cd4e69f2da2
--- /dev/null
+++ b/Reconstruction/eflowRec/src/PFAlgorithm.cxx
@@ -0,0 +1,60 @@
+#include "eflowRec/eflowCaloObject.h"
+#include "eflowRec/PFISubtractionTool.h"
+#include "eflowRec/PFAlgorithm.h"
+
+PFAlgorithm::PFAlgorithm(const std::string& name, ISvcLocator* pSvcLocator) : AthAlgorithm(name, pSvcLocator), m_tools(this), m_eflowRecTracksReadHandle("eflowRecTracks"), m_eflowRecClustersReadHandle("eflowRecCluster"), m_caloClustersWriteHandle("PFCaloCluster")
+{
+ declareProperty( "PrivateToolList", m_tools, "List of Private Subtraction eflowISubtractionAlgTools" );
+ declareProperty("eflowRecTracksInputName", m_eflowRecTracksReadHandle);
+ declareProperty("eflowRecClustersInputName", m_eflowRecClustersReadHandle);
+ declareProperty("PFCaloClustersOutputName", m_caloClustersWriteHandle);
+}
+
+StatusCode PFAlgorithm::initialize(){
+
+ ATH_CHECK(m_tools.retrieve());
+
+ ATH_CHECK(m_eflowRecTracksReadHandle.initialize());
+ ATH_CHECK(m_eflowRecClustersReadHandle.initialize());
+
+ ATH_CHECK(m_caloClustersWriteHandle.initialize());
+
+ printTools();
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode PFAlgorithm::execute(){
+
+ ATH_MSG_DEBUG("Executing");
+
+ eflowCaloObjectContainer theElowCaloObjectContainer;
+
+ /* Run the SubtractionTools
+ * --> CellLevelSubtractionTool, RecoverSplitShowersTool */
+ for (auto thisPFISubtractionTool : m_tools){
+ thisPFISubtractionTool->execute(&theElowCaloObjectContainer,const_cast<eflowRecTrackContainer*>(m_eflowRecTracksReadHandle.ptr()),const_cast<eflowRecClusterContainer*>(m_eflowRecClustersReadHandle.ptr()));
+ }
+
+ return StatusCode::SUCCESS;
+}
+
+StatusCode PFAlgorithm::finalize(){ return StatusCode::SUCCESS;}
+
+void PFAlgorithm::printTools() {
+ ATH_MSG_VERBOSE(" ");
+ ATH_MSG_VERBOSE("List of tools in execution sequence:");
+ ATH_MSG_VERBOSE("------------------------------------");
+ ATH_MSG_VERBOSE(" ");
+ unsigned int toolCtr = 0;
+ for (auto thisPFISubtractionTool : m_tools){
+ toolCtr++;
+ ATH_MSG_VERBOSE(std::setw(2) << std::setiosflags(std::ios_base::right) << toolCtr << ".) "
+ << std::resetiosflags(std::ios_base::right) << std::setw(36) << std::setfill('.')
+ << std::setiosflags(std::ios_base::left) << thisPFISubtractionTool->type() << std::setfill('.')
+ << thisPFISubtractionTool->name() << std::setfill(' '));
+ }
+ ATH_MSG_VERBOSE(" ");
+ ATH_MSG_VERBOSE("------------------------------------");
+}
+
diff --git a/Reconstruction/eflowRec/src/PFCellLevelSubtractionTool.cxx b/Reconstruction/eflowRec/src/PFCellLevelSubtractionTool.cxx
new file mode 100644
index 00000000000..dc561448296
--- /dev/null
+++ b/Reconstruction/eflowRec/src/PFCellLevelSubtractionTool.cxx
@@ -0,0 +1,465 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "eflowRec/PFCellLevelSubtractionTool.h"
+
+#include "eflowRec/eflowRecCluster.h"
+#include "eflowRec/eflowRecTrack.h"
+#include "eflowRec/eflowTrackClusterLink.h"
+#include "eflowRec/eflowDatabase.h"
+#include "eflowRec/eflowDepthCalculator.h"
+#include "eflowRec/eflowTrackCaloPoints.h"
+#include "eflowRec/IEFlowCellEOverPTool.h"
+#include "eflowRec/eflowEEtaBinnedParameters.h"
+#include "eflowRec/eflowLayerIntegrator.h"
+#include "eflowRec/PFTrackClusterMatchingTool.h"
+#include "eflowRec/cycle.h"
+#include "eflowRec/eflowRingSubtractionManager.h"
+#include "eflowRec/eflowCellSubtractionFacilitator.h"
+#include "eflowRec/eflowSubtractor.h"
+#include "eflowRec/eflowRingThicknesses.h"
+
+#include "CaloEvent/CaloCluster.h"
+
+#include "xAODCaloEvent/CaloCluster.h"
+#include "xAODCaloEvent/CaloClusterKineHelper.h"
+
+#include "xAODPFlow/PFO.h"
+
+#include "eflowRec/eflowCaloObject.h"
+#include "eflowRec/eflowCaloObjectMaker.h"
+#include "GaudiKernel/MsgStream.h"
+
+#include "GaudiKernel/ListItem.h"
+
+#include <algorithm>
+#include <iostream>
+#include <cmath>
+#include <vector>
+
+
+using namespace eflowSubtract;
+
+PFCellLevelSubtractionTool::PFCellLevelSubtractionTool(const std::string& type,const std::string& name,const IInterface* parent) :
+ AthAlgTool( type, name, parent),
+ m_matchingTool("PFTrackClusterMatchingTool/CalObjBldMatchingTool", this),
+ m_matchingToolForPull_015("PFTrackClusterMatchingTool/PFPullMatchingTool_015", this),
+ m_matchingToolForPull_02("PFTrackClusterMatchingTool/PFPullMatchingTool_02", this),
+ m_binnedParameters(nullptr),
+ m_integrator(nullptr),
+ m_theEOverPTool("eflowCellEOverPTool",this),
+ //m_rCell(0.75),
+ m_subtractionSigmaCut(1.5),
+ m_consistencySigmaCut(1.0),
+ m_calcEOverP(false),
+ m_goldenModeString(""),
+ m_nMatchesInCellLevelSubtraction(1),
+ m_useUpdated2015ChargedShowerSubtraction(true)
+
+{
+ declareProperty("PFTrackClusterMatchingTool", m_matchingTool, "The track-cluster matching tool");
+ declareProperty("PFTrackClusterMatchingTool_015", m_matchingToolForPull_015, "The 0.15 track-cluster matching tool to calculate the pull");
+ declareProperty("PFTrackClusterMatchingTool_02", m_matchingToolForPull_02, "The 0.2 track-cluster matching tool to calculate the pull");
+ declareProperty("eflowCellEOverPTool", m_theEOverPTool, "Energy Flow E/P Values and Shower Paremeters Tool");
+ declareProperty("SubtractionSigmaCut",m_subtractionSigmaCut);
+ declareProperty("ConsistencySigmaCut",m_consistencySigmaCut);
+ declareProperty("CalcEOverP",m_calcEOverP,"Whether to disable energy flow");
+ declareProperty("goldenModeString",m_goldenModeString,"run in golden match mode only?");
+ declareProperty("nMatchesInCellLevelSubtraction",m_nMatchesInCellLevelSubtraction,"Number of clusters to match");
+ declareProperty("useUpdated2015ChargedShowerSubtraction",m_useUpdated2015ChargedShowerSubtraction,"Toggle whether to use updated 2015 charged shower subtraction, which disables the shower subtraction in high calorimeter energy density region");
+}
+
+PFCellLevelSubtractionTool::~PFCellLevelSubtractionTool() {
+}
+
+StatusCode PFCellLevelSubtractionTool::initialize(){
+ StatusCode sc;
+
+ /* tool service */
+ IToolSvc* myToolSvc;
+ sc = service("ToolSvc", myToolSvc);
+
+ if (m_matchingTool.retrieve().isFailure()) {
+ msg(MSG::WARNING) << "Cannot find PFTrackClusterMatchingTool" << endmsg;
+ }
+
+ if (m_matchingToolForPull_015.retrieve().isFailure()) {
+ msg(MSG::WARNING) << "Cannot find PFTrackClusterMatchingTool_2" << endmsg;
+ }
+
+ if (m_matchingToolForPull_02.retrieve().isFailure()) {
+ msg(MSG::WARNING) << "Cannot find PFTrackClusterMatchingTool_2" << endmsg;
+ }
+
+ m_integrator = std::make_unique<eflowLayerIntegrator>(0.032, 1.0e-3, 3.0);
+ m_binnedParameters = std::make_unique<eflowEEtaBinnedParameters>();
+
+ sc = m_theEOverPTool->execute(m_binnedParameters.get());
+
+ if (sc.isFailure()) {
+ msg(MSG::WARNING) << "Could not execute eflowCellEOverPTool " << endmsg;
+ return StatusCode::SUCCESS;
+ }
+
+ return StatusCode::SUCCESS;
+
+}
+
+void PFCellLevelSubtractionTool::execute(eflowCaloObjectContainer* theEflowCaloObjectContainer, eflowRecTrackContainer* recTrackContainer, eflowRecClusterContainer* recClusterContainer){
+
+ ATH_MSG_VERBOSE("Executing PFCellLevelSubtractionTool");
+
+ m_eflowCaloObjectContainer = theEflowCaloObjectContainer;
+ m_eflowTrackContainer = recTrackContainer;
+ m_eflowClusterContainer = recClusterContainer;
+
+ /* Add each track to its best matching cluster's eflowCaloObject */
+ matchAndCreateEflowCaloObj(m_nMatchesInCellLevelSubtraction);
+
+ if (msgLvl(MSG::DEBUG)) printAllClusters(*recClusterContainer);
+
+ /* Check e/p mode - only perform subtraction if not in this mode */
+ if (!m_calcEOverP) {performSubtraction();}
+
+ /* Check e/p mode - only perform radial profiles calculations if in this mode */
+ if (m_calcEOverP) {calculateRadialEnergyProfiles();}
+
+}
+
+int PFCellLevelSubtractionTool::matchAndCreateEflowCaloObj(int n) {
+ int nMatches(0);
+
+ /* Create eflowTrackClusterLink after matching */
+ unsigned int nTrack = m_eflowTrackContainer->size();
+ for (unsigned int iTrack=0; iTrack<nTrack; ++iTrack) {
+ eflowRecTrack *thisEfRecTrack = static_cast<eflowRecTrack*>(m_eflowTrackContainer->at(iTrack));
+
+ /* Add cluster matches needed for pull calculation*/
+ std::vector<eflowRecCluster*> bestClusters_015 = m_matchingToolForPull_015->doMatches(thisEfRecTrack, m_eflowClusterContainer, -1);
+ std::vector<eflowRecCluster*> bestClusters_02 = m_matchingToolForPull_02->doMatches(thisEfRecTrack, m_eflowClusterContainer, -1);
+
+ for (unsigned int imatch=0; imatch < bestClusters_015.size(); ++imatch) {
+ eflowTrackClusterLink* trackClusterLink = eflowTrackClusterLink::getInstance(thisEfRecTrack, bestClusters_015.at(imatch));
+ thisEfRecTrack->addAlternativeClusterMatch(trackClusterLink,"cone_015");
+ }
+
+ for (unsigned int imatch=0; imatch < bestClusters_02.size(); ++imatch) {
+ eflowTrackClusterLink* trackClusterLink = eflowTrackClusterLink::getInstance(thisEfRecTrack, bestClusters_02.at(imatch));
+ thisEfRecTrack->addAlternativeClusterMatch(trackClusterLink,"cone_02");
+ }
+
+ std::vector<eflowRecCluster*> bestClusters = m_matchingTool->doMatches(thisEfRecTrack, m_eflowClusterContainer, n);
+ if (bestClusters.empty()) { continue; }
+
+ /* Matched cluster: create TrackClusterLink and add it to both the track and the cluster (eflowCaloObject will be created later) */
+ nMatches++;
+ for (unsigned int imatch=0; imatch < bestClusters.size(); ++imatch) {
+ eflowTrackClusterLink* trackClusterLink = eflowTrackClusterLink::getInstance(thisEfRecTrack, bestClusters.at(imatch));
+ thisEfRecTrack->addClusterMatch(trackClusterLink);
+ bestClusters.at(imatch)->addTrackMatch(trackClusterLink);
+
+ }
+
+ }
+
+ /* Create 3 types eflowCaloObjects: track-only, cluster-only, track-cluster-link */
+ eflowCaloObjectMaker makeCaloObject;
+ int nCaloObjects = makeCaloObject.makeTrkCluCaloObjects(m_eflowTrackContainer, m_eflowClusterContainer, m_eflowCaloObjectContainer);
+ msg(MSG::DEBUG) << "PFCellLevelSubtractionTool created total " << nCaloObjects << " CaloObjects." << endmsg;
+
+ /* integrate cells; determine FLI; eoverp */
+ for (unsigned int iCalo=0; iCalo<m_eflowCaloObjectContainer->size(); ++iCalo) {
+ eflowCaloObject *thisEflowCaloObject = static_cast<eflowCaloObject*>(m_eflowCaloObjectContainer->at(iCalo));
+
+ thisEflowCaloObject->simulateShower(m_integrator.get(), m_binnedParameters.get(), m_useUpdated2015ChargedShowerSubtraction);
+ }
+
+ return nMatches;
+}
+
+void PFCellLevelSubtractionTool::calculateRadialEnergyProfiles(){
+
+ msg(MSG::DEBUG) << "Accessed radial energy profile function" << std::endl;
+
+ unsigned int nEFCaloObs = m_eflowCaloObjectContainer->size();
+
+ for (unsigned int iEFCalOb = 0; iEFCalOb < nEFCaloObs; ++iEFCalOb) {
+
+ eflowCaloObject* thisEflowCaloObject = m_eflowCaloObjectContainer->at(iEFCalOb);
+
+ unsigned int nClusters = thisEflowCaloObject->nClusters();
+
+ if (nClusters < 1) {
+ continue;
+ }
+
+ const std::vector<eflowTrackClusterLink*>& matchedTrackList = thisEflowCaloObject->efRecLink();
+
+ int nTrackMatches = thisEflowCaloObject->nTracks();
+
+ for (int iTrack = 0; iTrack < nTrackMatches; ++iTrack) {
+
+ eflowRecTrack* efRecTrack = matchedTrackList[iTrack]->getTrack();
+
+ std::vector<eflowRecCluster*> matchedClusters;
+ matchedClusters.clear();
+ std::vector<eflowTrackClusterLink*> links = efRecTrack->getClusterMatches();
+ std::vector<eflowTrackClusterLink*>::iterator itLink = links.begin();
+ std::vector<eflowTrackClusterLink*>::iterator endLink = links.end();
+
+ for (; itLink != endLink; ++itLink) {
+ matchedClusters.push_back((*itLink)->getCluster());
+ }
+
+ 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(eflowCaloObject::getClusterContainerPtr()));
+ }
+
+ eflowCellList calorimeterCellList;
+ Subtractor::makeOrderedCellList(efRecTrack->getTrackCaloPoints(),clusterSubtractionList,calorimeterCellList);
+
+ eflowRingThicknesses ringThicknessGenerator;
+
+ std::vector<int> layerToStoreVector;
+ std::vector<float> radiusToStoreVector;
+ std::vector<float> avgEdensityToStoreVector;
+
+ for (int i=0; i < eflowCalo::nRegions ;i++){
+
+ eflowCaloENUM layer = (eflowCaloENUM)i;
+ msg(MSG::DEBUG) <<"layer is "<<layer<<std::endl;
+ double ringThickness = ringThicknessGenerator.ringThickness((eflowCaloENUM)i);
+ msg(MSG::DEBUG) <<"ring thickness is "<<ringThickness<<std::endl;
+
+ double eta_extr = calorimeterCellList.etaFF(layer);
+ msg(MSG::DEBUG) <<"extrapolated eta ["<<layer<<"] is "<<eta_extr<<std::endl;
+ double phi_extr = calorimeterCellList.phiFF(layer);
+ msg(MSG::DEBUG) <<"extrapolated phi ["<<layer<<"] is "<<phi_extr<<std::endl;
+
+ if (eta_extr == -999.0){
+ continue;
+ }
+
+ int indexofRing = 0;
+ int layerToStore = -999;
+
+ double radiusToStore = 0;
+ double totalEnergyPerCell = 0;
+
+ double energyDensityPerCell = -666;
+ double totalEnergyPerRing = 0;
+
+ double cellVolume = 0;
+ int totalCells = 0;
+
+ int n;
+ /* 100 is chosen as a number higher than the number of cells found in a normal list */
+ for (n=1; n<100; n++){
+
+ CellIt beginRing = calorimeterCellList.getLowerBound((eflowCaloENUM)i, ringThickness*(n-1));
+
+ if(beginRing == calorimeterCellList.end()){
+ break;
+ }
+
+ indexofRing += 1;
+ std::vector<std::pair<CaloCell*,int> > tempVector = (*beginRing).second;
+ std::vector<std::pair<CaloCell*,int> >::iterator firstPair = tempVector.begin();
+ std::vector<std::pair<CaloCell*,int> >::iterator lastPair = tempVector.end();
+
+ int totalCellsinRing = 0;
+ double energyDensityPerRing = 0;
+ double averageEnergyDensityPerRing = 0;
+
+ for (; firstPair != lastPair; ++firstPair) {
+ const CaloDetDescrElement* DDE = ((*firstPair).first)->caloDDE();
+ CaloCell_ID::CaloSample sampling = DDE->getSampling();
+
+ msg(MSG::DEBUG) << " cell eta and phi are " << ((*firstPair).first)->eta() << " and " << ((*firstPair).first)->phi() << " with index " << (*firstPair).second << " and sampling of " << sampling << std::endl;
+ msg(MSG::DEBUG) << " cell energy is " << ((*firstPair).first)->energy()<<std::endl;
+
+ totalCells += 1;
+ totalCellsinRing += 1;
+
+ totalEnergyPerRing += ((*firstPair).first)->energy();
+ totalEnergyPerCell = ((*firstPair).first)->energy();
+ msg(MSG::DEBUG) << " Total E per Cell is " << totalEnergyPerCell<<std::endl;
+ msg(MSG::DEBUG) << " Total E per Ring is " << totalEnergyPerRing<<std::endl;
+
+ cellVolume = DDE->volume();
+ msg(MSG::DEBUG) << " cell volume is " << cellVolume/1000.<<std::endl;
+
+ energyDensityPerCell = totalEnergyPerCell/(cellVolume/1000.);
+ msg(MSG::DEBUG) << " E density per Cell is " << energyDensityPerCell<<std::endl;
+ msg(MSG::DEBUG) << " Initial added E density per Cell is " << energyDensityPerRing<<std::endl;
+ energyDensityPerRing += energyDensityPerCell;
+ msg(MSG::DEBUG) << " Final added E density per Cell is " << energyDensityPerRing<<std::endl;
+ averageEnergyDensityPerRing = energyDensityPerRing/((totalCellsinRing)*(efRecTrack->getTrack()->e()/1000.));
+ }
+
+ msg(MSG::DEBUG) << " track E is " << efRecTrack->getTrack()->e()/1000.;
+ msg(MSG::DEBUG) << " Average E density per Ring is " << averageEnergyDensityPerRing<<std::endl;
+
+ if (averageEnergyDensityPerRing != 0){
+ avgEdensityToStoreVector.push_back(averageEnergyDensityPerRing);
+ layerToStore = (eflowCaloENUM)i;
+ msg(MSG::DEBUG) <<"layerToStore is "<< layerToStore << std::endl;
+ layerToStoreVector.push_back(layerToStore);
+ radiusToStore = (indexofRing)*ringThickness;
+ msg(MSG::DEBUG) <<"radiusToStore is "<< radiusToStore << std::endl;
+ radiusToStoreVector.push_back(radiusToStore);
+ }
+ else {msg(MSG::DEBUG) <<"averageEnergyDensityPerRing = 0"<<std::endl;}
+ }//loop on 100 cells
+ }//loop on calo regions
+
+ efRecTrack->setLayerCellOrderVector(layerToStoreVector);
+ efRecTrack->setRadiusCellOrderVector(radiusToStoreVector);
+ efRecTrack->setAvgEDensityCellOrderVector(avgEdensityToStoreVector);
+
+ layerToStoreVector.clear();
+ radiusToStoreVector.clear();
+ avgEdensityToStoreVector.clear();
+ }//loop on track matches
+ }//loop on eflowCaloObjects
+}
+
+void PFCellLevelSubtractionTool::performSubtraction() {
+
+ unsigned int nEFCaloObs = m_eflowCaloObjectContainer->size();
+ for (unsigned int iEFCalOb = 0; iEFCalOb < nEFCaloObs; ++iEFCalOb) {
+ eflowCaloObject* thisEflowCaloObject = m_eflowCaloObjectContainer->at(iEFCalOb);
+
+ unsigned int nClusters = thisEflowCaloObject->nClusters();
+ if (nClusters < 1) {
+ continue;
+ }
+
+ /* Get matched cluster via Links */
+
+ double expectedEnergy = thisEflowCaloObject->getExpectedEnergy();
+ double clusterEnergy = thisEflowCaloObject->getClusterEnergy();
+ double expectedSigma = sqrt(thisEflowCaloObject->getExpectedVariance());
+
+ /* Check e/p */
+ if (isEOverPFail(expectedEnergy, expectedSigma, clusterEnergy, m_consistencySigmaCut,
+ runInGoldenMode())
+ || (runInGoldenMode() && thisEflowCaloObject->nTracks() > 1)) {
+ continue;
+ }
+
+ /* Do subtraction */
+
+ int nTrackMatches = thisEflowCaloObject->nTracks();
+ if (nTrackMatches == 0 || runInGoldenMode()) {
+ continue;
+ }
+
+ /* Subtract the track from all matched clusters */
+ const std::vector<eflowTrackClusterLink*>& matchedTrackList = thisEflowCaloObject->efRecLink();
+
+ for (int iTrack = 0; iTrack < nTrackMatches; ++iTrack) {
+ eflowRecTrack* efRecTrack = matchedTrackList[iTrack]->getTrack();
+
+ /* Can't subtract without e/p */
+ if (!efRecTrack->hasBin()) {
+ continue;
+ }
+
+ if (efRecTrack->isInDenseEnvironment()) continue;
+
+ std::vector<eflowRecCluster*> matchedClusters;
+ matchedClusters.clear();
+ std::vector<eflowTrackClusterLink*> links = efRecTrack->getClusterMatches();
+ std::vector<eflowTrackClusterLink*>::iterator itLink = links.begin();
+ std::vector<eflowTrackClusterLink*>::iterator endLink = links.end();
+ for (; itLink != endLink; ++itLink) {
+ matchedClusters.push_back((*itLink)->getCluster());
+ }
+ 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(eflowCaloObject::getClusterContainerPtr()));
+ }
+
+ Subtractor::subtractTracksFromClusters(efRecTrack, clusterSubtractionList);
+
+ /* Annihilate the cluster(s) if the remnant is small (i.e. below k*sigma) */
+ if (canAnnihilated(0, expectedSigma, clusterEnergy)) {
+ Subtractor::annihilateClusters(clusterSubtractionList);
+ }
+ }
+
+ if (canAnnihilated(expectedEnergy, expectedSigma, clusterEnergy)) {
+ /* Check if we can annihilate right away */
+ std::vector<xAOD::CaloCluster*> clusterList;
+ unsigned nCluster = thisEflowCaloObject->nClusters();
+ for (unsigned iCluster = 0; iCluster < nCluster; ++iCluster) {
+ clusterList.push_back(
+ thisEflowCaloObject->efRecCluster(iCluster)->getClusterForModification(
+ eflowCaloObject::getClusterContainerPtr()));
+ }
+ Subtractor::annihilateClusters(clusterList);
+ }
+
+ /* Flag tracks as subtracted */
+ for (unsigned int iTrack = 0; iTrack < thisEflowCaloObject->nTracks(); ++iTrack) {
+ if (!thisEflowCaloObject->efRecTrack(iTrack)->isInDenseEnvironment()) thisEflowCaloObject->efRecTrack(iTrack)->setSubtracted();
+ }
+ }
+}
+
+StatusCode PFCellLevelSubtractionTool::finalize(){
+
+ return StatusCode::SUCCESS;
+}
+
+
+bool PFCellLevelSubtractionTool::isEOverPFail(double expectedEnergy, double sigma, double clusterEnergy, bool consistencySigmaCut, bool useGoldenMode) {
+ if ((expectedEnergy == 0) && (clusterEnergy > 0)) {
+ return false;
+ }
+
+ bool result = useGoldenMode ? fabs(clusterEnergy - expectedEnergy) > consistencySigmaCut*sigma
+ : clusterEnergy < expectedEnergy - consistencySigmaCut*sigma;
+ return result;
+}
+
+bool PFCellLevelSubtractionTool::canAnnihilated(double expectedEnergy, double sigma, double clusterEnergy) {
+ return clusterEnergy - expectedEnergy < m_subtractionSigmaCut * sigma;
+}
+
+std::string PFCellLevelSubtractionTool::printTrack(const xAOD::TrackParticle* track) {
+ std::stringstream result;
+ result << " track with E, eta and phi "<< track->e() << ", " << track->eta() << " and " << track->phi();
+ return result.str();
+}
+
+std::string PFCellLevelSubtractionTool::printCluster(const xAOD::CaloCluster* cluster) {
+ std::stringstream result;
+ result << " cluster with E, eta and phi of " << cluster->e() << ", " << cluster->eta() << " and " << cluster->phi();
+ return result.str();
+}
+
+void PFCellLevelSubtractionTool::printAllClusters(const eflowRecClusterContainer& recClusterContainer) {
+ unsigned int nClusters = recClusterContainer.size();
+ for (unsigned int iCluster = 0; iCluster < nClusters; ++iCluster) {
+ /* Create the eflowCaloObject and put it in the container */
+ const eflowRecCluster* thisEFRecCluster = recClusterContainer.at(iCluster);
+ if (thisEFRecCluster->getTrackMatches().empty()) {
+ ATH_MSG_DEBUG("Isolated" << printCluster(thisEFRecCluster->getCluster()));
+ } else {
+ ATH_MSG_DEBUG("Matched" << printCluster(thisEFRecCluster->getCluster()));
+ int nTrackMatches = thisEFRecCluster->getNTracks();
+ std::vector<eflowTrackClusterLink*> theTrackLinks = thisEFRecCluster->getTrackMatches();
+ for (int iTrackMatch = 0; iTrackMatch < nTrackMatches; ++iTrackMatch) {
+ ATH_MSG_DEBUG("Matched" << printTrack(theTrackLinks[iTrackMatch]->getTrack()->getTrack()));
+ }
+ }
+ }
+}
diff --git a/Reconstruction/eflowRec/src/components/eflowRec_entries.cxx b/Reconstruction/eflowRec/src/components/eflowRec_entries.cxx
index bba1efacb4a..51cceaf9e60 100644
--- a/Reconstruction/eflowRec/src/components/eflowRec_entries.cxx
+++ b/Reconstruction/eflowRec/src/components/eflowRec_entries.cxx
@@ -18,6 +18,8 @@
#include "eflowRec/PFLeptonSelector.h"
#include "eflowRec/PFTrackSelector.h"
#include "eflowRec/PFClusterSelector.h"
+#include "eflowRec/PFAlgorithm.h"
+#include "eflowRec/PFCellLevelSubtractionTool.h"
#include "GaudiKernel/DeclareFactoryEntries.h"
DECLARE_ALGORITHM_FACTORY( eflowBuilder )
@@ -29,6 +31,8 @@ DECLARE_ALGORITHM_FACTORY( eflowVertexInformationSetter )
DECLARE_ALGORITHM_FACTORY( PFLeptonSelector )
DECLARE_ALGORITHM_FACTORY( PFClusterSelector )
DECLARE_ALGORITHM_FACTORY( PFTrackSelector )
+DECLARE_ALGORITHM_FACTORY( PFAlgorithm )
+DECLARE_TOOL_FACTORY( PFCellLevelSubtractionTool )
DECLARE_TOOL_FACTORY( eflowRecoverSplitShowersTool )
DECLARE_TOOL_FACTORY( eflowCellLevelSubtractionTool )
DECLARE_TOOL_FACTORY( eflowLCCalibTool )
@@ -51,6 +55,8 @@ DECLARE_FACTORY_ENTRIES(eflowRec) {
DECLARE_ALGORITHM( PFLeptonSelector )
DECLARE_ALGORITHM( PFClusterSelector )
DECLARE_ALGORITHM( PFTrackSelector )
+ DECLARE_ALGORITHM( PFAlgorithm )
+ DECLARE_TOOL( PFCellLevelSubtractionTool )
DECLARE_TOOL ( eflowRecoverSplitShowersTool )
DECLARE_TOOL ( eflowCellLevelSubtractionTool )
DECLARE_TOOL ( eflowMomentCalculatorTool )
--
GitLab