First round of updates to set subtracted energy, instead of bool....

First round of updates to set subtracted energy, instead of bool. PFCellLevelSubtractionTool migrated.

Event/xAOD/xAODPFlow/xAODPFlow/selection.xml

@@ -40,7 +40,9 @@
   <class name="ElementLink<xAOD::TrackCaloClusterContainer_v1>" />
   <class name="std::vector<ElementLink<xAOD::TrackCaloClusterContainer_v1> >" />
   <class name="std::vector<std::vector<ElementLink<xAOD::TrackCaloClusterContainer_v1> > >" />
-      
+
+  <!-- Extra types needed -->
+  <class name="std::vector<std::vector<std::pair<ElementLink<DataVector<xAOD::CaloCluster_v1> >,double> > >" />
   <!-- Suppress the unwanted classes found by ROOT 6. -->
   <!-- Hopefully we can remove these extra lines at one point... -->
   <exclusion>

Event/xAOD/xAODPFlow/xAODPFlow/xAODPFlowDict.h

@@ -50,6 +50,8 @@ namespace {
     ElementLink< xAOD::TrackCaloClusterContainer_v1 >                               l15;
     std::vector< ElementLink< xAOD::TrackCaloClusterContainer_v1 > >                l16;
     std::vector< std::vector< ElementLink< xAOD::TrackCaloClusterContainer_v1 > > > l17;
+
+    std::vector<std::vector<std::pair<ElementLink<DataVector<xAOD::CaloCluster_v1> >,double> > > l18;
     
    };
 }

Reconstruction/eflowRec/eflowRec/PFSubtractionStatusSetter.h

@@ -14,7 +14,7 @@ class eflowCaloObject;
 class PFSubtractionStatusSetter {
 
 public:
-  void markSubtractionStatus(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusterList, eflowCaloObject& thisEflowCaloObject);
-    
+  void markSubtractionStatus(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusterList, std::vector<float> clusterSubtractedEnergyRatios,eflowCaloObject& thisEflowCaloObject);
+  void markAnnihStatus(eflowCaloObject& thisEflowCaloObject);  
 };
 #endif

Reconstruction/eflowRec/eflowRec/eflowCaloObject.h

@@ -54,9 +54,10 @@ public:
     }
   }
 
-  /* For a specific eflowTrackClusterLink indicate whether or not it has been fully/partially subtracted via the bool */
-  /* True indicates it has been fully or partially subtracted, whilst false indicates it has not been modified at all */
-  void setTrackClusterLinkSubtractionStatus(unsigned int index, bool status) { m_trackClusterLinks[index].second = status; }
+  /* For a specific eflowTrackClusterLink indicate whether or not it has been fully/partially subtracted by setting the energy ratio
+  ** of subtracted cluster energy to original cluster enegry.
+  A value other than nan indicates it has been fully or partially subtracted */
+  void setTrackClusterLinkSubtractionStatus(unsigned int index, float energyRatio) { m_trackClusterLinks[index].second = energyRatio; }
   
   /* Track accessor methods */
   eflowRecTrack* efRecTrack(int i) const { return m_eflowRecTracks[i]; }
@@ -70,7 +71,7 @@ public:
 
 
   /* Link accessor methods */
-  std::vector<std::pair<eflowTrackClusterLink*,bool> > efRecLink() const { return m_trackClusterLinks; }
+  std::vector<std::pair<eflowTrackClusterLink*,float> > efRecLink() const { return m_trackClusterLinks; }
   void clearLinks() { m_trackClusterLinks.clear(); }
 
   /* Calculate total tracks energy, total tracks energy variance, total cluster energy for subtraction */
@@ -82,14 +83,17 @@ public:
 
 private:
 
-  void addTrackClusterLink(eflowTrackClusterLink* trackClusterLink) { m_trackClusterLinks.push_back(std::pair(trackClusterLink,false)); }
+  void addTrackClusterLink(eflowTrackClusterLink* trackClusterLink) { m_trackClusterLinks.push_back(std::pair(trackClusterLink,NAN)); }
 
  private:
   /* Vector of clusters */
   std::vector<eflowRecCluster*> m_eflowRecClusters;
 
-  /* Vector of track-cluster matches - the bool is to be used to indicate whether a cluster was subtracted or not */
-  std::vector<std::pair<eflowTrackClusterLink*,bool> > m_trackClusterLinks;
+  /* Vector of track-cluster matches - the float is to be used to indicate the ratio of subtracted cluster energy
+  ** to original unsubtracted cluster energy. It is initialiased to nan, and you hence one can use isNan to verify
+  ** whether the cluster was subtracted or not (we leave it as nan it not).
+  */
+  std::vector<std::pair<eflowTrackClusterLink*,float> > m_trackClusterLinks;
 
   /* Vector of tracks */
   std::vector<eflowRecTrack*> m_eflowRecTracks;

Reconstruction/eflowRec/src/PFCellLevelSubtractionTool.cxx

@@ -36,6 +36,7 @@
 #include <iostream>
 #include <cmath>
 #include <vector>
+#include <map>
 
 
 using namespace eflowSubtract;
@@ -175,7 +176,7 @@ void PFCellLevelSubtractionTool::calculateRadialEnergyProfiles(){
     continue;
     }
     
-    const std::vector<std::pair<eflowTrackClusterLink*,bool> >& matchedTrackList = thisEflowCaloObject->efRecLink();
+    const std::vector<std::pair<eflowTrackClusterLink*,float> >& matchedTrackList = thisEflowCaloObject->efRecLink();
 
     int nTrackMatches = thisEflowCaloObject->nTracks();
     
@@ -351,65 +352,77 @@ void PFCellLevelSubtractionTool::performSubtraction() {
 
     ATH_MSG_DEBUG("About to perform subtraction for this eflowCaloObject");
 
-    const std::vector<std::pair<eflowTrackClusterLink*, bool> >& matchedTrackList = thisEflowCaloObject->efRecLink();
+    const std::vector<std::pair<eflowTrackClusterLink*, float> >& matchedTrackList = thisEflowCaloObject->efRecLink();
 
+    /* Check if we can annihilate right away */
     if (canAnnihilated(expectedEnergy, expectedSigma, clusterEnergy)) {
-
-      /* Check if we can annihilate right away */
+      
       std::vector<std::pair<xAOD::CaloCluster*, bool> > clusterList;
       unsigned nCluster = thisEflowCaloObject->nClusters();
-      for (unsigned iCluster = 0; iCluster < nCluster; ++iCluster) {
-        clusterList.push_back(std::pair(thisEflowCaloObject->efRecCluster(iCluster)->getCluster(),false));
-      }
-      Subtractor::annihilateClusters(clusterList);
+      for (unsigned iCluster = 0; iCluster < nCluster; ++iCluster) clusterList.push_back(std::pair(thisEflowCaloObject->efRecCluster(iCluster)->getCluster(),false));
+
+      Subtractor::annihilateClusters(clusterList);      
 
       //Now we should mark all of these clusters as being subtracted
-      pfSubtractionStatusSetter.markSubtractionStatus(clusterList, *thisEflowCaloObject);
+      pfSubtractionStatusSetter.markAnnihStatus(*thisEflowCaloObject);
 
     } else {
 
       /* Subtract the track from all matched clusters */
       
       for (int iTrack = 0; iTrack < nTrackMatches; ++iTrack) {
-	eflowRecTrack* efRecTrack = (matchedTrackList[iTrack].first)->getTrack();
+	      eflowRecTrack* efRecTrack = (matchedTrackList[iTrack].first)->getTrack();
 	
-	ATH_MSG_DEBUG("Have got eflowRecTrack number " << iTrack << " for this eflowCaloObject");
+	      ATH_MSG_DEBUG("Have got eflowRecTrack number " << iTrack << " for this eflowCaloObject");
 	
-	/* Can't subtract without e/p */
-	if (!efRecTrack->hasBin()) {
-	  continue;
-	}
+	      /* Can't subtract without e/p */
+	      if (!efRecTrack->hasBin()) {
+	        continue;
+	      }
      
-	if (efRecTrack->isInDenseEnvironment()) continue;
+	      if (efRecTrack->isInDenseEnvironment()) continue;
 
-	ATH_MSG_DEBUG("Have bin and am not in dense environment for this eflowCaloObject");
+	      ATH_MSG_DEBUG("Have bin and am not in dense environment for this eflowCaloObject");
       
-	std::vector<eflowRecCluster*> matchedClusters;
-	matchedClusters.clear();
-	std::vector<eflowTrackClusterLink*> links = efRecTrack->getClusterMatches();
-	for (auto thisEFlowTrackClusterLink : links) matchedClusters.push_back(thisEFlowTrackClusterLink->getCluster());
+	      std::vector<eflowRecCluster*> matchedClusters;
+	      matchedClusters.clear();
+	      std::vector<eflowTrackClusterLink*> links = efRecTrack->getClusterMatches();
+	      for (auto thisEFlowTrackClusterLink : links) matchedClusters.push_back(thisEFlowTrackClusterLink->getCluster());
 
-	ATH_MSG_DEBUG("Have filled matchedClusters list for this eflowCaloObject");
+	      ATH_MSG_DEBUG("Have filled matchedClusters list for this eflowCaloObject");
 	
-	std::vector<std::pair<xAOD::CaloCluster*, bool> > clusterSubtractionList;
-	for (auto thisEFlowRecCluster : matchedClusters) clusterSubtractionList.push_back(std::pair(thisEFlowRecCluster->getCluster(),false));
-
-	ATH_MSG_DEBUG("Have filled clusterSubtractionList for this eflowCaloObject");
+	      std::vector<std::pair<xAOD::CaloCluster*, bool> > clusterSubtractionList;
+        std::map<xAOD::CaloCluster*, double> clusterEnergyMap;
+	      for (auto thisEFlowRecCluster : matchedClusters) {
+          xAOD::CaloCluster* thisCluster = thisEFlowRecCluster->getCluster();
+          clusterSubtractionList.push_back(std::pair(thisCluster,false));
+          clusterEnergyMap[thisCluster] = thisCluster->e();
+        }
+
+	      ATH_MSG_DEBUG("Have filled clusterSubtractionList for this eflowCaloObject");
       
-	Subtractor::subtractTracksFromClusters(efRecTrack, clusterSubtractionList);
-	//Now need to mark which clusters were modified in the subtraction procedure
-	pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, *thisEflowCaloObject);
+	      Subtractor::subtractTracksFromClusters(efRecTrack, clusterSubtractionList);
 	
-	ATH_MSG_DEBUG("Have performed subtraction for this eflowCaloObject");
+	      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)) {
-	  Subtractor::annihilateClusters(clusterSubtractionList);
-	  //Now we should mark all of these clusters as being subtracted
-	  pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, *thisEflowCaloObject);
-	}
+	      /* Annihilate the cluster(s) if the remnant is small (i.e. below k*sigma) */
+	      if (canAnnihilated(0, expectedSigma, clusterEnergy)) {
+	        Subtractor::annihilateClusters(clusterSubtractionList);
+	        //Now we should mark all of these clusters as being subtracted
+          pfSubtractionStatusSetter.markAnnihStatus(*thisEflowCaloObject);	       
+	      }
+        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);
+	
+        } 
 	
-	ATH_MSG_DEBUG("Have checked if can annihilate clusters for this eflowCaloOject");
+	      ATH_MSG_DEBUG("Have checked if can annihilate clusters for this eflowCaloOject");
       
       }
     }

Reconstruction/eflowRec/src/PFOChargedCreatorAlgorithm.cxx

@@ -104,7 +104,7 @@ void PFOChargedCreatorAlgorithm::createChargedPFO(const eflowCaloObject& energyF
     /* Optionally we add the links to clusters to the xAOD::PFO */
     if (true == addClusters){
 
-      std::vector<std::pair<eflowTrackClusterLink*,bool> > trackClusterLinkPairs = energyFlowCaloObject.efRecLink();
+      std::vector<std::pair<eflowTrackClusterLink*,float> > trackClusterLinkPairs = energyFlowCaloObject.efRecLink();
 
       std::vector<eflowTrackClusterLink*> thisTracks_trackClusterLinks = efRecTrack->getClusterMatches();
 
@@ -114,15 +114,27 @@ void PFOChargedCreatorAlgorithm::createChargedPFO(const eflowCaloObject& energyF
 
       std::vector<eflowTrackClusterLink*> thisTracks_trackClusterLinksSubtracted;
 
+      //Create vector of pairs which map each CaloCluster to the ratio of its new energy to unstracted energy
+      std::vector<std::pair<ElementLink<xAOD::CaloClusterContainer>, double> > vectorClusterToSubtractedEnergies;
+
       for (auto trackClusterLink : thisTracks_trackClusterLinks){
         for (auto trackClusterLinkPair : trackClusterLinkPairs){
-          if (!m_eOverPMode && trackClusterLinkPair.first == trackClusterLink && true == trackClusterLinkPair.second) {
+          if (!m_eOverPMode && trackClusterLinkPair.first == trackClusterLink && !std::isnan(trackClusterLinkPair.second)) {
             thisTracks_trackClusterLinksSubtracted.push_back(trackClusterLink);
+            eflowRecCluster* efRecCluster = trackClusterLinkPair.first->getCluster();
+            ElementLink<xAOD::CaloClusterContainer> theOriginalClusterLink = efRecCluster->getOriginalClusElementLink();
+    	      ElementLink<xAOD::CaloClusterContainer> theSisterClusterLink = (*theOriginalClusterLink)->getSisterClusterLink();
+            if (theSisterClusterLink.isValid()) vectorClusterToSubtractedEnergies.push_back(std::pair(theSisterClusterLink,trackClusterLinkPair.second));
+            else vectorClusterToSubtractedEnergies.push_back(std::pair(theOriginalClusterLink,trackClusterLinkPair.second));
           }
           else if (m_eOverPMode) thisTracks_trackClusterLinksSubtracted.push_back(trackClusterLink);
         }
       }
 
+      //sort the vectorClusterToSubtractedEnergies in order of subtracted energy ratio from low (most subtracted) to high (least subtracted)
+      std::sort(vectorClusterToSubtractedEnergies.begin(),vectorClusterToSubtractedEnergies.end(), [](auto const& a, auto const&b){return a.second < b.second;});
+      thisPFO->setAttribute("PF_vectorClusterToSubtractedEnergies",vectorClusterToSubtractedEnergies);
+
       //Now loop over the list of eflowTrackClusterLink which correspond to subtracted clusters matched to this track.
       bool isFirstCluster = true;
       for (auto trackClusterLink : thisTracks_trackClusterLinksSubtracted){	    

Reconstruction/eflowRec/src/PFRecoverSplitShowersTool.cxx

@@ -224,29 +224,43 @@ void PFRecoverSplitShowersTool::performSubtraction(eflowCaloObject* thisEflowCal
     /* Do subtraction */
     std::vector<std::pair<xAOD::CaloCluster*, bool> > clusterSubtractionList;
     clusterSubtractionList.reserve(matchedClusters.size());
-    for (auto thisEFlowRecCluster : matchedClusters) clusterSubtractionList.push_back(std::pair(thisEFlowRecCluster->getCluster(),false));
-    
+    std::map<xAOD::CaloCluster*, double> clusterEnergyMap;
+    for (auto thisEFlowRecCluster : matchedClusters) {
+      xAOD::CaloCluster* thisCluster = thisEFlowRecCluster->getCluster();
+      clusterSubtractionList.push_back(std::pair(thisCluster,false));
+      clusterEnergyMap[thisCluster] = thisCluster->e();
+    }
+
     if (getSumEnergy(clusterSubtractionList) - thisEfRecTrack->getEExpect() < m_subtractionSigmaCut
         * sqrt(thisEfRecTrack->getVarEExpect())) {
       /* Check if we can annihilate right away */
       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.markSubtractionStatus(clusterSubtractionList, *thisEflowCaloObject);
-
+      pfSubtractionStatusSetter.markAnnihStatus(*thisEflowCaloObject);
     } else {
       /* Subtract the track from all matched clusters */
       Subtractor::subtractTracksFromClusters(thisEfRecTrack, clusterSubtractionList);
-      //Now need to mark which clusters were modified in the subtraction procedure
-      pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, *thisEflowCaloObject);
-      
+          
       /* 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);
-	//Now we should mark all of these clusters as being subtracted
-	pfSubtractionStatusSetter.markSubtractionStatus(clusterSubtractionList, *thisEflowCaloObject);
+	      //Now we should mark all of these clusters as being subtracted
+	      pfSubtractionStatusSetter.markAnnihStatus(*thisEflowCaloObject);
       }
+      else{
+        //Else now to mark which clusters were modified in the subtraction procedure
+        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);
+	
+        } 
+
     }
     /* Flag tracks as subtracted */
     thisEfRecTrack->setSubtracted();

Reconstruction/eflowRec/src/PFSubtractionStatusSetter.cxx

@@ -8,30 +8,38 @@
 #include "eflowRec/eflowRecCluster.h"
 #include "eflowRec/eflowTrackClusterLink.h"
 
-void PFSubtractionStatusSetter::markSubtractionStatus(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusterList, eflowCaloObject& thisEflowCaloObject){
+void PFSubtractionStatusSetter::markSubtractionStatus(const std::vector<std::pair<xAOD::CaloCluster*, bool> >& clusterList, std::vector<float> clusterSubtractedEnergyRatios, eflowCaloObject& thisEflowCaloObject){
   //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
 
-  const std::vector<std::pair<eflowTrackClusterLink*, bool> >& matchedTrackList = thisEflowCaloObject.efRecLink();
+  const std::vector<std::pair<eflowTrackClusterLink*, float> >& matchedTrackList = thisEflowCaloObject.efRecLink();
 
+  unsigned int clusCounter = 0;
   for (auto thisClusterPair : clusterList){
     xAOD::CaloCluster* thisCluster = thisClusterPair.first;
     unsigned int counter = 0;
     for (auto& thisTrackClusterLinkPair : matchedTrackList){
       //if the subtraction status is already true, then no need to update it
-      if (true == thisTrackClusterLinkPair.second) {
-	counter++;
-	continue;
+      if (!std::isnan(thisTrackClusterLinkPair.second)) {
+    	  counter++;
+	      continue;
       }
       //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 (true == thisClusterPair.second) thisEflowCaloObject.setTrackClusterLinkSubtractionStatus(counter, true);
+	      if (true == thisClusterPair.second && !std::isnan(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){
+
+  unsigned int numTrackClusterLinks = thisEflowCaloObject.efRecLink().size();
+  for (unsigned int counter = 0 ; counter < numTrackClusterLinks; ++counter) thisEflowCaloObject.setTrackClusterLinkSubtractionStatus(counter, 0.0);
+
+}
+