Merge branch 'mhodkgin_improveDEBUG_June2021' into 'master'

Particle Flow bug fixes to content of Global Particle Flow variables and FPE

See merge request !44488

Reconstruction/eflowRec/eflowRec/PFSubtractionEnergyRatioCalculator.h

+/*                                                                                                                                                                                                                                          
+  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration                                                                                                                                                                   
+*/
+
+#ifndef PFSUBTRACTIONENERGYRATIOCALCULATOR_H
+#define PFSUBTRACTIONENERGYRATIOCALCULATOR_H
+
+#include "AthenaBaseComps/AthMessaging.h"
+
+#include "xAODCaloEvent/CaloCluster.h"
+
+#include <vector>
+#include <map>
+#include <utility>
+
+/** Class to calculate the ratio of new to old energies of CaloClusters after the particle flow
+charged shower subtraction has been run **/
+
+class PFSubtractionEnergyRatioCalculator : public AthMessaging {
+
+  public:
+  PFSubtractionEnergyRatioCalculator();
+  /** For each xAOD::CaloCluster in clusterSubtractionList we calculate the ratio of new to old energy after the charged shower subtraction, which
+  is added in clusterSubtractedEnergyRatios. clusterEnergyMap contains the cluster energy prior to the charged shower energy subtraction */
+  void calculateSubtractedEnergyRatios(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusterSubtractionList, std::map<xAOD::CaloCluster*, double>& clusterEnergyMap, std::vector<float>& clusterSubtractedEnergyRatios);
+  /** If we have decided to annihiliate all clusters in clusterSubtractionList we use this function to set all ratios to zero in clusterEnergySubtractionRatios.
+  If the old energy, in clusterEnergyMap, was zero we instead set the ratio to NAN to denote that no energy was removed from that xAOD::CaloCluster */
+  void calculateSubtractedEnergyRatiosForAnnih(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusterSubtractionList, std::map<xAOD::CaloCluster*, double>& clusterEnergyMap, std::vector<float>& clusterSubtractedEnergyRatios);
+
+};
+#endif

Reconstruction/eflowRec/eflowRec/PFSubtractionStatusSetter.h

@@ -5,16 +5,25 @@
 #ifndef PFSUBTRACTIONSTATUSSETTER_H
 #define PFSUBTRACTIONSTATUSSETTER_H
 
+#include "AthenaBaseComps/AthMessaging.h"
+
 #include "xAODCaloEvent/CaloCluster.h"
 class eflowCaloObject;
 
 #include <vector>
 #include <utility>
 
-class PFSubtractionStatusSetter {
+/** This class contains a few functions to set the amount of energy removed from a xAOD::CaloCluster by a xAOD::Track in eflowRec.
+It can either set this for a specific track attached to an eflowCaloObject or set it for all tracks in an eflowCaloObject **/
+
+class PFSubtractionStatusSetter : public AthMessaging  {
 
 public:
-  void markSubtractionStatus(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusterList, std::vector<float> clusterSubtractedEnergyRatios,eflowCaloObject& thisEflowCaloObject);
-  void markAnnihStatus(eflowCaloObject& thisEflowCaloObject);  
+  PFSubtractionStatusSetter();
+
+  /** Set the ratio of new to old cluster energy for each cluster matched to a track with trackIndex */
+  void markSubtractionStatus(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusterList, std::vector<float>& clusterSubtractedEnergyRatios,eflowCaloObject& thisEflowCaloObject, unsigned int trackIndex);    
+  /** Set the ratio of new to old cluster energy, to zero, for all cluster matched to all tracks attached to the eflowCaloObject */
+  void markAllTracksAnnihStatus(eflowCaloObject& thisEflowCaloObject);  
 };
 #endif

Reconstruction/eflowRec/src/PFCellLevelSubtractionTool.cxx

@@ -19,6 +19,7 @@
 #include "eflowRec/eflowSubtractor.h"
 #include "eflowRec/eflowRingThicknesses.h"
 #include "eflowRec/PFSubtractionStatusSetter.h"
+#include "eflowRec/PFSubtractionEnergyRatioCalculator.h"
 
 #include "CaloEvent/CaloCluster.h"
 
@@ -210,7 +211,7 @@ void PFCellLevelSubtractionTool::calculateRadialEnergyProfiles(eflowData& data)
       std::vector<std::pair<xAOD::CaloCluster*, bool> > clusterSubtractionList;
       clusterSubtractionList.reserve(matchedClusters.size());
 
-for (auto *thisEFlowRecCluster : matchedClusters) clusterSubtractionList.emplace_back(thisEFlowRecCluster->getCluster(),false);
+      for (auto *thisEFlowRecCluster : matchedClusters) clusterSubtractionList.emplace_back(thisEFlowRecCluster->getCluster(),false);
 
       eflowCellList calorimeterCellList;
       Subtractor::makeOrderedCellList(efRecTrack->getTrackCaloPoints(),clusterSubtractionList,calorimeterCellList);
@@ -329,6 +330,11 @@ void PFCellLevelSubtractionTool::performSubtraction(eflowData& data) const {
   ATH_MSG_DEBUG("In performSubtraction");
 
   PFSubtractionStatusSetter pfSubtractionStatusSetter;
+  PFSubtractionEnergyRatioCalculator pfSubtractionEnergyRatioCalculator;
+  if( msgLevel( MSG::DEBUG ) ) {
+    pfSubtractionStatusSetter.setLevel(MSG::DEBUG);
+    pfSubtractionEnergyRatioCalculator.setLevel(MSG::DEBUG);
+  }
 
   unsigned int nEFCaloObs = data.caloObjects->size();
   for (unsigned int iEFCalOb = 0; iEFCalOb < nEFCaloObs; ++iEFCalOb) {
@@ -366,28 +372,42 @@ void PFCellLevelSubtractionTool::performSubtraction(eflowData& data) const {
       continue;
     }
 
-    ATH_MSG_DEBUG("About to perform subtraction for this eflowCaloObject");
-
     const std::vector<std::pair<eflowTrackClusterLink*, float> >& matchedTrackList = thisEflowCaloObject->efRecLink();
 
+    if( msgLevel( MSG::DEBUG ) ){ 
+      for (unsigned int iTrack = 0; iTrack < nTrackMatches; ++iTrack) {
+	      const xAOD::TrackParticle* thisTrack = (matchedTrackList[iTrack].first)->getTrack()->getTrack();
+        ATH_MSG_DEBUG("eflowCaloObject has track with E, pt and eta " << thisTrack->e() << ", " << thisTrack->pt() << " and " << thisTrack->eta());
+      }
+    }
+
+    ATH_MSG_DEBUG("About to perform subtraction for this eflowCaloObject");
+
     /* Check if we can annihilate right away */
     if (canAnnihilated(expectedEnergy, expectedSigma, clusterEnergy)) {
 
       std::vector<std::pair<xAOD::CaloCluster*, bool> > clusterList;
+      std::map<xAOD::CaloCluster*, double> clusterEnergyMap;
       unsigned nCluster = thisEflowCaloObject->nClusters();
       for (unsigned iCluster = 0; iCluster < nCluster; ++iCluster) clusterList.emplace_back(thisEflowCaloObject->efRecCluster(iCluster)->getCluster(),false);
 
+      ATH_MSG_DEBUG("We are going to annihilate. ExpectedEnergy, expectedSigma and clusterEnergy are " << expectedEnergy << ", " << expectedSigma << " and " << clusterEnergy);
+      if( msgLevel( MSG::DEBUG ) ) for (auto thisPair : clusterList) ATH_MSG_DEBUG("Annihilating cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
+
       Subtractor::annihilateClusters(clusterList);
 
       //Now we should mark all of these clusters as being subtracted      
-      pfSubtractionStatusSetter.markAnnihStatus(*thisEflowCaloObject);
+      pfSubtractionStatusSetter.markAllTracksAnnihStatus(*thisEflowCaloObject);
 
     } else {
 
       /* Subtract the track from all matched clusters */
 
+      ATH_MSG_DEBUG("Subtract tracks one by one");
+
       for (unsigned int iTrack = 0; iTrack < nTrackMatches; ++iTrack) {
 	      eflowRecTrack* efRecTrack = (matchedTrackList[iTrack].first)->getTrack();
+        unsigned int trackIndex = efRecTrack->getTrack()->index();
 
 	      ATH_MSG_DEBUG("Have got eflowRecTrack number " << iTrack << " for this eflowCaloObject");
 
@@ -417,25 +437,25 @@ void PFCellLevelSubtractionTool::performSubtraction(eflowData& data) const {
 
 	      ATH_MSG_DEBUG("Have filled clusterSubtractionList for this eflowCaloObject");
 
-	      Subtractor::subtractTracksFromClusters(efRecTrack, clusterSubtractionList);
+        bool debugToggle = false;
+        if( msgLevel( MSG::DEBUG ) ) debugToggle = true;
+	      Subtractor::subtractTracksFromClusters(efRecTrack, clusterSubtractionList,debugToggle);
 
 	      ATH_MSG_DEBUG("Have performed subtraction for this eflowCaloObject");
 
 	      /* Annihilate the cluster(s) if the remnant is small (i.e. below k*sigma) */
 	      if (canAnnihilated(0, expectedSigma, clusterEnergy)) {
+          if( msgLevel( MSG::DEBUG ) ) for (auto thisPair : clusterSubtractionList) ATH_MSG_DEBUG("Annihilating cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
 	        Subtractor::annihilateClusters(clusterSubtractionList);
 	        //Now we should mark all of these clusters as being subtracted
-          pfSubtractionStatusSetter.markAnnihStatus(*thisEflowCaloObject);
+          std::vector<float> clusterSubtractedEnergyRatios;
+          pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatiosForAnnih(clusterSubtractionList,clusterEnergyMap,clusterSubtractedEnergyRatios);           
+          pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject, trackIndex);
 	      }
         else{
           std::vector<float> clusterSubtractedEnergyRatios;
-          for (auto thisCluster: clusterSubtractionList) {
-            if (fabs(thisCluster.first->e() - clusterEnergyMap[thisCluster.first]) > 0.0001) clusterSubtractedEnergyRatios.push_back(thisCluster.first->e()/clusterEnergyMap[thisCluster.first]);
-            else clusterSubtractedEnergyRatios.push_back(NAN);
-          }
-
-	        pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject);
-
+          pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatios(clusterSubtractionList,clusterEnergyMap,clusterSubtractedEnergyRatios);          
+	        pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject, trackIndex);
         }
 
 	      ATH_MSG_DEBUG("Have checked if can annihilate clusters for this eflowCaloOject");

Reconstruction/eflowRec/src/PFChargedFlowElementCreatorAlgorithm.cxx

@@ -110,8 +110,12 @@ void PFChargedFlowElementCreatorAlgorithm::createChargedFlowElements(const eflow
 
       std::vector<std::pair<eflowTrackClusterLink*,float> > trackClusterLinkPairs = energyFlowCaloObject.efRecLink();
 
+      ATH_MSG_DEBUG("Have got " << trackClusterLinkPairs.size() << " trackClusterLinkPairs");
+
       std::vector<eflowTrackClusterLink*> thisTracks_trackClusterLinks = efRecTrack->getClusterMatches();
 
+      ATH_MSG_DEBUG("Have got " << thisTracks_trackClusterLinks.size() << " cluster matches");
+
       /** Each eflowCaloObject has a list of clusters for all the tracks it represents.
        *  We only want the subset of the clusters matched to this track, and collect these in thisTracks_trackClusterLinksSubtracted.
        */
@@ -128,6 +132,7 @@ void PFChargedFlowElementCreatorAlgorithm::createChargedFlowElements(const eflow
             eflowRecCluster* efRecCluster = trackClusterLinkPair.first->getCluster();
             ElementLink<xAOD::CaloClusterContainer> theOriginalClusterLink = efRecCluster->getOriginalClusElementLink();
             ElementLink<xAOD::CaloClusterContainer> theSisterClusterLink = (*theOriginalClusterLink)->getSisterClusterLink();
+	          ATH_MSG_DEBUG("Will add cluster with E and ratio " << (*theOriginalClusterLink)->e() << " and " << trackClusterLinkPair.second);
             if (theSisterClusterLink.isValid()) vectorClusterToSubtractedEnergies.emplace_back(std::pair(theSisterClusterLink,trackClusterLinkPair.second));
             else vectorClusterToSubtractedEnergies.emplace_back(std::pair(theOriginalClusterLink,trackClusterLinkPair.second));
           }

Reconstruction/eflowRec/src/PFRecoverSplitShowersTool.cxx

@@ -20,6 +20,7 @@
 #include "eflowRec/eflowCellSubtractionFacilitator.h"
 #include "eflowRec/eflowSubtractor.h"
 #include "eflowRec/PFSubtractionStatusSetter.h"
+#include "eflowRec/PFSubtractionEnergyRatioCalculator.h"
 
 #include "CaloEvent/CaloClusterContainer.h"
 #include "xAODCaloEvent/CaloClusterKineHelper.h"
@@ -211,9 +212,16 @@ unsigned int PFRecoverSplitShowersTool::matchAndCreateEflowCaloObj(eflowData& da
 void PFRecoverSplitShowersTool::performSubtraction(eflowCaloObject* thisEflowCaloObject) const {
 
   PFSubtractionStatusSetter pfSubtractionStatusSetter;
+  PFSubtractionEnergyRatioCalculator pfSubtractionEnergyRatioCalculator;
+  if( msgLevel( MSG::DEBUG ) ) {
+    pfSubtractionStatusSetter.setLevel(MSG::DEBUG);
+    pfSubtractionEnergyRatioCalculator.setLevel(MSG::DEBUG);
+  }
 
   for (unsigned iTrack = 0; iTrack < thisEflowCaloObject->nTracks(); ++iTrack) {
     eflowRecTrack* thisEfRecTrack = thisEflowCaloObject->efRecTrack(iTrack);
+    unsigned int trackIndex = thisEfRecTrack->getTrack()->index();
+
     ATH_MSG_DEBUG("About to recover track with e, pt, eta and phi of " << thisEfRecTrack->getTrack()->e() << ", " << thisEfRecTrack->getTrack()->pt() << ", " << thisEfRecTrack->getTrack()->eta() << " and "
     << thisEfRecTrack->getTrack()->eta());
     /* Get matched cluster via Links */
@@ -236,42 +244,40 @@ void PFRecoverSplitShowersTool::performSubtraction(eflowCaloObject* thisEflowCal
       clusterEnergyMap[thisCluster] = thisCluster->e();
     }
 
+    ATH_MSG_DEBUG("Have filled clusterSubtractionList for this eflowCaloObject");
+
     if (getSumEnergy(clusterSubtractionList) - thisEfRecTrack->getEExpect() < m_subtractionSigmaCut
         * sqrt(thisEfRecTrack->getVarEExpect())) {
       /* Check if we can annihilate right away */
+      if( msgLevel( MSG::DEBUG ) ) for (auto thisPair : clusterSubtractionList) ATH_MSG_DEBUG("Annihilating cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
+    
       Subtractor::annihilateClusters(clusterSubtractionList);
       //Now we should mark all of these clusters as being subtracted
       //Now need to mark which clusters were modified in the subtraction procedure
-      pfSubtractionStatusSetter.markAnnihStatus(*thisEflowCaloObject);
+      std::vector<float> clusterSubtractedEnergyRatios;
+      pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatiosForAnnih(clusterSubtractionList,clusterEnergyMap,clusterSubtractedEnergyRatios);       
+      pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject,trackIndex);	
     } else {
       /* Subtract the track from all matched clusters */
       const bool debugToggle = msgLvl(MSG::DEBUG);
       Subtractor::subtractTracksFromClusters(thisEfRecTrack, clusterSubtractionList, debugToggle);                
 
+      ATH_MSG_DEBUG("Have performed subtraction for this eflowCaloObject");
+
       /* Annihilate the cluster(s) if the remnant is small (i.e. below k*sigma) */      
       if (getSumEnergy(clusterSubtractionList) < m_subtractionSigmaCut
           * sqrt(thisEfRecTrack->getVarEExpect())) {
+        if( msgLevel( MSG::DEBUG ) ) for (auto thisPair : clusterSubtractionList) ATH_MSG_DEBUG("Annihilating remnant cluster with E and eta " << thisPair.first->e() << " and " << thisPair.first->eta());
         Subtractor::annihilateClusters(clusterSubtractionList);
 	      //Now we should mark all of these clusters as being subtracted
-	      pfSubtractionStatusSetter.markAnnihStatus(*thisEflowCaloObject);
+        std::vector<float> clusterSubtractedEnergyRatios;
+        pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatiosForAnnih(clusterSubtractionList,clusterEnergyMap,clusterSubtractedEnergyRatios); 	      
+        pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject,trackIndex);	
       }
       else{
-        //Else now to mark which clusters were modified in the subtraction procedure
         std::vector<float> clusterSubtractedEnergyRatios;        
-        for (auto thisCluster: clusterSubtractionList) {          
-          //clusterEnergyMap[thisCluster.first can be zero, but this is only a problem if thisCluster.first.e() < 0 - this can happen if a cluster starts with E =0 
-          //from a previous shower subtraction step and then we subtract more such that the new E is < 0. Then the ratio would cause an FPE due to the zero.
-          //If both thisCluster.first->e() and clusterEnergyMap[thisCluster.first] are zero we never enter this step.          
-          if (fabs(thisCluster.first->e() - clusterEnergyMap[thisCluster.first]) > 0.0001){
-            if ( clusterEnergyMap[thisCluster.first] >= 0) clusterSubtractedEnergyRatios.push_back(thisCluster.first->e()/clusterEnergyMap[thisCluster.first]);
-            //approximate zero with 0.0001 to avoid FPE and still give a meaningful ratio (e.g -100/0.0001)
-            else clusterSubtractedEnergyRatios.push_back(thisCluster.first->e()/0.0001);
-          }          
-          else clusterSubtractedEnergyRatios.push_back(NAN);
-        }
-
-	pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject);
-	
+        pfSubtractionEnergyRatioCalculator.calculateSubtractedEnergyRatios(clusterSubtractionList,clusterEnergyMap,clusterSubtractedEnergyRatios);    
+	      pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, clusterSubtractedEnergyRatios, *thisEflowCaloObject,trackIndex);	
       } 
 
     }

Reconstruction/eflowRec/src/PFSubtractionEnergyRatioCalculator.cxx

+#include "eflowRec/PFSubtractionEnergyRatioCalculator.h"
+
+#include "GaudiKernel/IMessageSvc.h"
+#include "GaudiKernel/ISvcLocator.h"
+#include "GaudiKernel/Bootstrap.h"
+
+#include "AthenaKernel/getMessageSvc.h"
+
+PFSubtractionEnergyRatioCalculator::PFSubtractionEnergyRatioCalculator() : AthMessaging(Gaudi::svcLocator()->service<IMessageSvc>("MessageSvc"),
+                  "PFSubtractionEnergyRatioCalculator"){}
+
+void PFSubtractionEnergyRatioCalculator::calculateSubtractedEnergyRatios(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusterSubtractionList, std::map<xAOD::CaloCluster*, double>& clusterEnergyMap, std::vector<float>& clusterSubtractedEnergyRatios){
+
+  ATH_MSG_DEBUG("Setting subtracted energy ratios here");
+  
+  for (auto thisCluster: clusterSubtractionList) {
+    ATH_MSG_DEBUG("Cluster energies are " << thisCluster.first->e() << " and " << clusterEnergyMap[thisCluster.first]);
+
+    //clusterEnergyMap[thisCluster.first can be zero, but this is only a problem if thisCluster.first.e() < 0 - this can happen if a cluster starts with E =0 
+    //from a previous shower subtraction step and then we subtract more such that the new E is < 0. Then the ratio would cause an FPE due to the zero.
+    //If both thisCluster.first->e() and clusterEnergyMap[thisCluster.first] are zero we never enter this step.          
+
+    if (std::abs(thisCluster.first->e() - clusterEnergyMap[thisCluster.first]) > 0.0001) {
+       if ( clusterEnergyMap[thisCluster.first] >= 0) {
+        ATH_MSG_DEBUG("Subtracted energy ratio is " << thisCluster.first->e()/clusterEnergyMap[thisCluster.first]);
+        clusterSubtractedEnergyRatios.push_back(thisCluster.first->e()/clusterEnergyMap[thisCluster.first]);
+       }
+       //approximate zero with 0.0001 to avoid FPE and still give a meaningful ratio (e.g -100/0.0001)
+       else {
+         ATH_MSG_DEBUG("Subtracted energy ratio is " << thisCluster.first->e()/0.0001);
+         clusterSubtractedEnergyRatios.push_back(thisCluster.first->e()/0.0001);
+       }
+    }
+    //else if the cluster enegry did not change then we use NAN to denote that no charged shower subtraction occurred.
+    else {
+      clusterSubtractedEnergyRatios.push_back(NAN);
+      ATH_MSG_DEBUG("Subtracted energy ratio is NAN ");
+    }
+  }//Loop over clusterSubtractionList
+
+}
+
+void PFSubtractionEnergyRatioCalculator::calculateSubtractedEnergyRatiosForAnnih(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusterSubtractionList, std::map<xAOD::CaloCluster*, double>& clusterEnergyMap, std::vector<float>& clusterSubtractedEnergyRatios){
+
+  ATH_MSG_DEBUG("Setting subtracted energy ratios for annih here");
+
+  for (auto thisCluster: clusterSubtractionList) {
+    ATH_MSG_DEBUG("Cluster energies are " << thisCluster.first->e() << " and " << clusterEnergyMap[thisCluster.first]);
+
+    //The energy can be zero if we previously "annihiliated" this cluster because we will have already set the energy to zero.
+    //We don't want to set the energy ratio to zero if we just change the energy from "0" to "0" - in effect in this case
+    //we don't subtract anything and we will set the value of the ratio to NAN which denotes that no subtraction was performed.
+    //So we build a list of indices of clusters which don't need to have the energy set to zero - we compare a floating point
+    //directly to zero, because in a previous annihilation we explicity set the energy to be exactly zero.
+    if (0 != clusterEnergyMap[thisCluster.first]) {
+      ATH_MSG_DEBUG("Setting cluster energy ratio to zero");
+      clusterSubtractedEnergyRatios.push_back(0);
+    }
+    else {
+      ATH_MSG_DEBUG("Setting cluster energy ratio to NAN");
+      clusterSubtractedEnergyRatios.push_back(NAN);
+    }
+
+  }
+
+}
\ No newline at end of file

Reconstruction/eflowRec/src/PFSubtractionStatusSetter.cxx

@@ -3,43 +3,81 @@
 */
 
 
+#include "GaudiKernel/IMessageSvc.h"
+#include "GaudiKernel/ISvcLocator.h"
+#include "GaudiKernel/Bootstrap.h"
+
+#include "AthenaKernel/getMessageSvc.h"
+
 #include "eflowRec/PFSubtractionStatusSetter.h"
 #include "eflowRec/eflowCaloObject.h"
 #include "eflowRec/eflowRecCluster.h"
 #include "eflowRec/eflowTrackClusterLink.h"
 
-void PFSubtractionStatusSetter::markSubtractionStatus(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusterList, std::vector<float> clusterSubtractedEnergyRatios, eflowCaloObject& thisEflowCaloObject){
+PFSubtractionStatusSetter::PFSubtractionStatusSetter() : AthMessaging(Gaudi::svcLocator()->service<IMessageSvc>("MessageSvc"),
+                  "PFSubtractionStatusSetter"){}
+
+void PFSubtractionStatusSetter::markSubtractionStatus(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusterList, std::vector<float>& clusterSubtractedEnergyRatios, eflowCaloObject& thisEflowCaloObject, unsigned int trackIndex){
   //An eflowCaloObject may have one cluster and N tracks, and then one would have N eflowTrackClusterLink* for each track-cluster pair
   //Hence there can be more entries in the track cluster link list due to duplication
 
+  ATH_MSG_DEBUG("Executing markSubtractionStatus and have clusterList of size " << clusterList.size() << ", clusterSubtractedEnergyRatios of size " << clusterSubtractedEnergyRatios.size() << " and trackIndex of " << trackIndex);  
+
   const std::vector<std::pair<eflowTrackClusterLink*, float> >& matchedTrackList = thisEflowCaloObject.efRecLink();
 
   unsigned int clusCounter = 0;
   for (auto thisClusterPair : clusterList){
     xAOD::CaloCluster* thisCluster = thisClusterPair.first;
+    ATH_MSG_DEBUG("Cluster with e and eta " << thisCluster->e() << " and " << thisCluster->eta());
     unsigned int counter = 0;
     for (const auto & thisTrackClusterLinkPair : matchedTrackList){
+
+      if (trackIndex != thisTrackClusterLinkPair.first->getTrack()->getTrack()->index()){
+        counter++;
+        continue;
+      }
+
+      ATH_MSG_DEBUG("Track with e and eta " << thisTrackClusterLinkPair.first->getTrack()->getTrack()->e() << " and " << thisTrackClusterLinkPair.first->getTrack()->getTrack()->eta());
+
       //if the subtraction status is already true, then no need to update it
       if (!std::isnan(thisTrackClusterLinkPair.second)) {
     	  counter++;
 	      continue;
       }
+
+      ATH_MSG_DEBUG("Not already subtracted");
+
       //eflowTrackCluster link returns an eflowRecCluster pointer, which in turn returns an xAOD:;CaloCluster* pointer
       xAOD::CaloCluster* thisMatchedTrackCluster = (thisTrackClusterLinkPair.first)->getCluster()->getCluster();
-      //Now we can do a floating point comparison of the energy to check which cluster we have
-      if (fabs(thisCluster->e() - thisMatchedTrackCluster->e()) < 0.0001){
-	      if (thisClusterPair.second && !std::isnan(clusterSubtractedEnergyRatios[clusCounter])) thisEflowCaloObject.setTrackClusterLinkSubtractionStatus(counter, clusterSubtractedEnergyRatios[clusCounter]);        
+      //Now we can do a pointer comparison of the two clusters to see if we have a match
+      ATH_MSG_DEBUG("Comparing clusters with energies " << thisCluster->e() << " and " << thisMatchedTrackCluster->e());
+      if (thisCluster == thisMatchedTrackCluster){      
+        ATH_MSG_DEBUG("Passed cluster energy check");
+	      if (thisClusterPair.second && !std::isnan(clusterSubtractedEnergyRatios[clusCounter])) {
+          ATH_MSG_DEBUG("Will set ratio to be " << clusterSubtractedEnergyRatios[clusCounter]);
+          thisEflowCaloObject.setTrackClusterLinkSubtractionStatus(counter, clusterSubtractedEnergyRatios[clusCounter]);        
+        }
       }//if have a match of the cluster
       counter++;
     }//loop on track cluster link pairs
     clusCounter++;
   }//loop on cluster pair list
+
 }
 
-void PFSubtractionStatusSetter::markAnnihStatus(eflowCaloObject& thisEflowCaloObject){
+void PFSubtractionStatusSetter::markAllTracksAnnihStatus(eflowCaloObject& thisEflowCaloObject){
+
+  ATH_MSG_DEBUG("Executing markAllTracksAnnihStatus");
+
+  const std::vector<std::pair<eflowTrackClusterLink*, float> >& matchedTrackList = thisEflowCaloObject.efRecLink();
   
-  unsigned int numTrackClusterLinks = thisEflowCaloObject.efRecLink().size();
-  for (unsigned int counter = 0 ; counter < numTrackClusterLinks; ++counter) thisEflowCaloObject.setTrackClusterLinkSubtractionStatus(counter, 0.0);
+  unsigned int counter = 0;
+  for (const auto & thisTrackClusterLinkPair : matchedTrackList){
+    ATH_MSG_DEBUG("Track with e and eta " << thisTrackClusterLinkPair.first->getTrack()->getTrack()->e() << " and " << thisTrackClusterLinkPair.first->getTrack()->getTrack()->eta());
+    ATH_MSG_DEBUG("will annihilate cluster with e and eta " << thisTrackClusterLinkPair.first->getCluster()->getCluster()->e() << " and " << thisTrackClusterLinkPair.first->getCluster()->getCluster()->eta());
+    thisEflowCaloObject.setTrackClusterLinkSubtractionStatus(counter, 0.0);  
+    counter++;
+  }
   
 }