From 01c4f2cffa2fce6a5498efaef62cddd5cd81d2ff Mon Sep 17 00:00:00 2001
From: Adam Bailey <adam.bailey@cern.ch>
Date: Thu, 21 May 2020 13:42:15 +0000
Subject: [PATCH] Updated TauSubstructureVariables and CaloClusterVariables to
 use new HelperFunction to navigate jet->clusters. Will now work with PFlow
 seeds

---
 .../tauRec/python/TauAlgorithmsHolder.py      |   9 +-
 Reconstruction/tauRec/python/tauRecFlags.py   |   2 +-
 .../tauRecTools/Root/CaloClusterVariables.cxx |  19 ++-
 .../tauRecTools/Root/HelperFunctions.cxx      | 125 +++++++++++-------
 .../Root/MvaTESVariableDecorator.cxx          |   7 +-
 .../Root/TauSubstructureVariables.cxx         |  26 ++--
 .../tauRecTools/src/TauAxisSetter.cxx         |  16 +--
 .../tauRecTools/src/TauAxisSetter.h           |   1 +
 .../tauRecTools/src/TauPi0ClusterCreator.cxx  |  15 +--
 .../tauRecTools/src/TauPi0ClusterCreator.h    |   2 +
 .../tauRecTools/CaloClusterVariables.h        |  15 ++-
 .../tauRecTools/tauRecTools/HelperFunctions.h |   2 +-
 .../tauRecTools/TauSubstructureVariables.h    |  19 +--
 13 files changed, 146 insertions(+), 112 deletions(-)

diff --git a/Reconstruction/tauRec/python/TauAlgorithmsHolder.py b/Reconstruction/tauRec/python/TauAlgorithmsHolder.py
index b9f73169804..b2df3e55227 100644
--- a/Reconstruction/tauRec/python/TauAlgorithmsHolder.py
+++ b/Reconstruction/tauRec/python/TauAlgorithmsHolder.py
@@ -77,7 +77,8 @@ def getTauAxis():
     from tauRecTools.tauRecToolsConf import TauAxisSetter
     TauAxisSetter = TauAxisSetter(  name = _name, 
                                     ClusterCone = 0.2,
-                                    VertexCorrection = True)
+                                    VertexCorrection = True,
+                                    IncShowerSubtr = tauFlags.useShowerSubClusters() )
                                     
     cached_instances[_name] = TauAxisSetter                
     return TauAxisSetter
@@ -311,7 +312,8 @@ def getTauSubstructure():
                                                           # parameters for CaloIsoCorrected variable
                                                           maxPileUpCorrection = 4000., #MeV
                                                           pileUpAlpha = 1.0,
-                                                          VertexCorrection = True
+                                                          VertexCorrection = True,
+                                                          IncShowerSubtr = tauFlags.useShowerSubClusters()
                                                        )
     
     cached_instances[_name] = TauSubstructureVariables
@@ -388,7 +390,8 @@ def getPi0ClusterCreator():
     
     from tauRecTools.tauRecToolsConf import TauPi0ClusterCreator
     TauPi0ClusterCreator = TauPi0ClusterCreator(name = _name,
-                                                Key_Pi0ClusterContainer="TauPi0SubtractedClusters"
+                                                Key_Pi0ClusterContainer="TauPi0SubtractedClusters",
+                                                IncShowerSubtr = tauFlags.useShowerSubClusters()
                                                 )
     
     cached_instances[_name] = TauPi0ClusterCreator
diff --git a/Reconstruction/tauRec/python/tauRecFlags.py b/Reconstruction/tauRec/python/tauRecFlags.py
index 6fa5ff40ea8..fef5f98b99b 100644
--- a/Reconstruction/tauRec/python/tauRecFlags.py
+++ b/Reconstruction/tauRec/python/tauRecFlags.py
@@ -193,7 +193,7 @@ class useOldVertexFitterAPI(JobProperty):
 class useShowerSubClusters(JobProperty):
     """ switch on use of shower subtracted clusters
     """
-    statusOn=False
+    statusOn=True
     allowedTypes=['bool']
     StoredValue=False
 
diff --git a/Reconstruction/tauRecTools/Root/CaloClusterVariables.cxx b/Reconstruction/tauRecTools/Root/CaloClusterVariables.cxx
index 59e14070850..c78183810d0 100644
--- a/Reconstruction/tauRecTools/Root/CaloClusterVariables.cxx
+++ b/Reconstruction/tauRecTools/Root/CaloClusterVariables.cxx
@@ -3,6 +3,7 @@
 */
 
 #include "tauRecTools/CaloClusterVariables.h"
+#include "tauRecTools/HelperFunctions.h"
 #include <cmath>
 
 const double CaloClusterVariables::DEFAULT = -1111.;
@@ -21,7 +22,8 @@ m_totMass(DEFAULT),
 m_effMass(DEFAULT),
 m_totEnergy(DEFAULT),
 m_effEnergy(DEFAULT),
-m_doVertexCorrection(false) {
+m_doVertexCorrection(false),
+m_incShowerSubtr(true){
 }
 
 //*******************************************
@@ -33,18 +35,15 @@ bool CaloClusterVariables::update(const xAOD::TauJet& pTau) {
     const xAOD::Jet* pSeed = *pTau.jetLink();
     if(!pSeed) return false;
 
-    xAOD::JetConstituentVector::const_iterator nav_it = pSeed->getConstituents().begin();
-    xAOD::JetConstituentVector::const_iterator nav_itE = pSeed->getConstituents().end();
-    const xAOD::CaloCluster* pCluster;
-   
-    std::vector<CaloVertexedClusterType> constituents;
-    for (; nav_it != nav_itE; ++nav_it) {
-      pCluster = dynamic_cast<const xAOD::CaloCluster*> ( (*nav_it)->rawConstituent() );
-      if (!pCluster) continue;
+    // Loop through jets, get links to clusters
+    std::vector<const xAOD::CaloCluster*> clusterList;
+    StatusCode sc = tauRecTools::GetJetClusterList(pSeed, clusterList, m_incShowerSubtr);
 
+    std::vector<CaloVertexedClusterType> constituents;
+    for (auto pCluster : clusterList){
       // correct cluster
       if (pTau.vertexLink() && m_doVertexCorrection)
-        constituents.emplace_back (*pCluster, (*pTau.vertexLink())->position());
+	constituents.emplace_back (*pCluster, (*pTau.vertexLink())->position());
       else
         constituents.emplace_back (*pCluster);
     }
diff --git a/Reconstruction/tauRecTools/Root/HelperFunctions.cxx b/Reconstruction/tauRecTools/Root/HelperFunctions.cxx
index bed43d1062c..890c326a52d 100644
--- a/Reconstruction/tauRecTools/Root/HelperFunctions.cxx
+++ b/Reconstruction/tauRecTools/Root/HelperFunctions.cxx
@@ -195,7 +195,7 @@ float tauRecTools::TRTBDT::GetResponse(){
   return this->bdt->GetResponse();
 }
 
-const StatusCode tauRecTools::GetJetClusterList(const xAOD::Jet* jet, std::vector<const xAOD::CaloCluster*> &clusterList, bool incShowerSubtracted){
+const StatusCode tauRecTools::GetJetClusterList(const xAOD::Jet* jet, std::vector<const xAOD::CaloCluster*> &clusterList, bool incShowerSubtracted, TLorentzVector dRVector, double dRCut){
 
   using namespace tauRecTools::msgHelperFunction;
   // If using subtracted clusters, need to store unmodified to check if charged are duplicates
@@ -206,44 +206,62 @@ const StatusCode tauRecTools::GetJetClusterList(const xAOD::Jet* jet, std::vecto
   for(auto jCon : jVec){
     ANA_MSG_DEBUG("JetConstituent: ");
     ANA_MSG_DEBUG("eta: " << jCon->eta() << " phi: " << jCon->phi() << " e: " << jCon->e()); 
-    if( jCon->type() == xAOD::Type::CaloCluster ) {
-      const xAOD::CaloCluster* cluster = static_cast<const xAOD::CaloCluster*>( jCon->rawConstituent() );
-      ANA_MSG_DEBUG("CaloCluster: ");
-      ANA_MSG_DEBUG("eta: " << cluster->eta() << " phi: " << cluster->phi() << " e: " << cluster->e());
-      ANA_MSG_DEBUG("rawEta: " << cluster->rawEta() << " rawPhi: " << cluster->rawPhi() << " rawE: " << cluster->rawE()); 
-      ANA_MSG_DEBUG("calEta: " << cluster->calEta() << " calPhi: " << cluster->calPhi() << " calE: " << cluster->calE());
-      clusterList.push_back(cluster);
+
+    // do deltaR check against jet constituent
+    bool PassdR = true;
+    if (dRCut > 0){
+      TLorentzVector tempClusterVector;
+      tempClusterVector.SetPtEtaPhiE( jCon->pt(), jCon->eta(), jCon->phi(), jCon->e() );
+      ANA_MSG_DEBUG("Apply dR cut on JetConstituent: " << dRCut );
+      ANA_MSG_DEBUG("JetConstituent Pt: " << tempClusterVector.Pt() << ", Eta: " << tempClusterVector.Eta() << ", Phi: " << tempClusterVector.Phi());
+      ANA_MSG_DEBUG("dR " << dRVector.DeltaR(tempClusterVector));
+      if (dRVector.DeltaR(tempClusterVector) > dRCut){
+	ANA_MSG_DEBUG("Failed dR Cut ");
+	PassdR = false;
+      }
     }
-    else if( jCon->type() == xAOD::Type::ParticleFlow ) {
-      const xAOD::PFO* pfo = static_cast<const xAOD::PFO*>( jCon->rawConstituent() );
-      if( !pfo->isCharged() ){
-	if( pfo->nCaloCluster()==1 ){
-
-	  if (incShowerSubtracted){
-	    ElementLink<xAOD::CaloClusterContainer> subClusLink;
-	    pfo->attribute("PFOShowerSubtractedClusterLink", subClusLink);
-	    if ( !subClusLink.isValid() ){
-	      ANA_MSG_ERROR("Tau HelperFunctions: Found invalid link to shower subtracted cluster");
-	      return StatusCode::FAILURE;
+
+    if (PassdR){
+      if( jCon->type() == xAOD::Type::CaloCluster ) {
+	const xAOD::CaloCluster* cluster = static_cast<const xAOD::CaloCluster*>( jCon->rawConstituent() );
+	ANA_MSG_DEBUG("CaloCluster: ");
+	ANA_MSG_DEBUG("eta: " << cluster->eta() << " phi: " << cluster->phi() << " e: " << cluster->e());
+	ANA_MSG_DEBUG("rawEta: " << cluster->rawEta() << " rawPhi: " << cluster->rawPhi() << " rawE: " << cluster->rawE());
+	ANA_MSG_DEBUG("calEta: " << cluster->calEta() << " calPhi: " << cluster->calPhi() << " calE: " << cluster->calE());
+
+	clusterList.push_back(cluster);
+      }
+      else if( jCon->type() == xAOD::Type::ParticleFlow ) {
+	const xAOD::PFO* pfo = static_cast<const xAOD::PFO*>( jCon->rawConstituent() );
+	if( !pfo->isCharged() ){
+	  if( pfo->nCaloCluster()==1 ){
+
+	    if (incShowerSubtracted){
+	      ElementLink<xAOD::CaloClusterContainer> subClusLink;
+	      pfo->attribute("PFOShowerSubtractedClusterLink", subClusLink);
+	      if ( !subClusLink.isValid() ){
+		ANA_MSG_ERROR("Tau HelperFunctions: Found invalid link to shower subtracted cluster");
+		return StatusCode::FAILURE;
+	      }
+	      else {
+		clusterList.push_back( (*subClusLink) );
+		dupList.push_back( pfo->cluster(0) );
+	      }
 	    }
 	    else {
-	      clusterList.push_back( (*subClusLink) );
-	      dupList.push_back( pfo->cluster(0) );
+	      clusterList.push_back(pfo->cluster(0));
 	    }
-	  }
-	  else {
-	    clusterList.push_back(pfo->cluster(0));
-	  }
 
-	}
-	else ANA_MSG_WARNING("Neutral PFO has " << std::to_string(pfo->nCaloCluster()) << " clusters, expected exactly 1!\n");
+	  }
+	  else ANA_MSG_WARNING("Neutral PFO has " << std::to_string(pfo->nCaloCluster()) << " clusters, expected exactly 1!\n");
 
-      }// neutral PFO check
-    }
-    else{
-      ANA_MSG_ERROR("GetJetConstCluster: Seed jet constituent type not supported!");
-      return StatusCode::FAILURE;
-    }
+	}// neutral PFO check
+      }
+      else{
+	ANA_MSG_ERROR("GetJetConstCluster: Seed jet constituent type not supported!");
+	return StatusCode::FAILURE;
+      }
+    }// dR check
   }
 
   // Get clusters from charged PFOs
@@ -252,21 +270,38 @@ const StatusCode tauRecTools::GetJetClusterList(const xAOD::Jet* jet, std::vecto
   else checkList = clusterList;
 
   for(auto jCon : jVec){
+
     if( jCon->type() == xAOD::Type::ParticleFlow ) {
-      const xAOD::PFO* pfo = static_cast<const xAOD::PFO*>( jCon->rawConstituent() );
-      if( pfo->isCharged() ){
-
-	// loop through clusters linked to charged PFO
-	for (u_int cc=0; cc<pfo->nCaloCluster(); cc++){
-	  const xAOD::CaloCluster* cluster = pfo->cluster(cc);
-	  // check it is not duplicate of one in neutral list
-	  if ( std::find(checkList.begin(), checkList.end(), cluster) == checkList.end() ){
-	    clusterList.push_back(cluster);
-	    checkList.push_back(cluster);
-	  }
-	}
 
+      bool PassdR = true;
+      if (dRCut > 0){
+	TLorentzVector tempClusterVector;
+	tempClusterVector.SetPtEtaPhiE( jCon->pt(), jCon->eta(), jCon->phi(), jCon->e() );
+	ANA_MSG_DEBUG("Apply dR cut on JetConstituent: " << dRCut );
+	ANA_MSG_DEBUG("JetConstituent Pt: " << tempClusterVector.Pt() << ", Eta: " << tempClusterVector.Eta() << ", Phi: " << tempClusterVector.Phi());
+	ANA_MSG_DEBUG("dR " << dRVector.DeltaR(tempClusterVector));
+	if (dRVector.DeltaR(tempClusterVector) > dRCut){
+	  ANA_MSG_DEBUG("Failed dR Cut ");
+	  PassdR = false;
+	}
       }
+
+      if (PassdR){
+	const xAOD::PFO* pfo = static_cast<const xAOD::PFO*>( jCon->rawConstituent() );
+	if( pfo->isCharged() ){
+
+	  // loop through clusters linked to charged PFO
+	  for (u_int cc=0; cc<pfo->nCaloCluster(); cc++){
+	    const xAOD::CaloCluster* cluster = pfo->cluster(cc);
+	    // check it is not duplicate of one in neutral list
+	    if ( std::find(checkList.begin(), checkList.end(), cluster) == checkList.end() ){
+	      clusterList.push_back(cluster);
+	      checkList.push_back(cluster);
+	    }
+	  }
+
+	}
+      }// dR check
     }
   }// loop through jet constituents
 
diff --git a/Reconstruction/tauRecTools/Root/MvaTESVariableDecorator.cxx b/Reconstruction/tauRecTools/Root/MvaTESVariableDecorator.cxx
index f19fbc769a8..313dd1d9e12 100644
--- a/Reconstruction/tauRecTools/Root/MvaTESVariableDecorator.cxx
+++ b/Reconstruction/tauRecTools/Root/MvaTESVariableDecorator.cxx
@@ -84,15 +84,16 @@ StatusCode MvaTESVariableDecorator::execute(xAOD::TauJet& xTau) {
   TLorentzVector clusters_had_P4;
   clusters_had_P4.SetPtEtaPhiM(0,0,0,0);
 
+  TLorentzVector LC_P4;
+  LC_P4.SetPtEtaPhiM(xTau.ptDetectorAxis(), xTau.etaDetectorAxis(), xTau.phiDetectorAxis(), xTau.m());
+
   // Loop through jets, get links to clusters
   std::vector<const xAOD::CaloCluster*> clusterList;
-  ATH_CHECK(tauRecTools::GetJetClusterList(jet_seed, clusterList, m_incShowerSubtr));
+  ATH_CHECK(tauRecTools::GetJetClusterList(jet_seed, clusterList, m_incShowerSubtr, LC_P4, 0.2));
 
   // Loop through clusters and jet constituents
   for (auto cl : clusterList){
 
-    if (xTau.p4(xAOD::TauJetParameters::DetectorAxis).DeltaR(cl->p4()) > 0.2) continue;
-
     clE = cl->calE();
     Etot += clE;
 
diff --git a/Reconstruction/tauRecTools/Root/TauSubstructureVariables.cxx b/Reconstruction/tauRecTools/Root/TauSubstructureVariables.cxx
index 17bfc4df285..9973557f3cf 100644
--- a/Reconstruction/tauRecTools/Root/TauSubstructureVariables.cxx
+++ b/Reconstruction/tauRecTools/Root/TauSubstructureVariables.cxx
@@ -20,6 +20,7 @@
 #include "tauRecTools/TauSubstructureVariables.h"
 
 #include "tauRecTools/KineUtils.h"
+#include "tauRecTools/HelperFunctions.h"
 
 #define GeV 1000
 const double TauSubstructureVariables::DEFAULT = -1111.;
@@ -33,6 +34,7 @@ TauSubstructureVariables::TauSubstructureVariables( const std::string& name )
 	declareProperty("maxPileUpCorrection", m_maxPileUpCorrection = 4 * GeV);
 	declareProperty("pileUpAlpha", m_pileUpAlpha = 1.0);
 	declareProperty("VertexCorrection", m_doVertexCorrection = false);
+	declareProperty("IncShowerSubtr", m_incShowerSubtr = true);
 }
 
 
@@ -79,6 +81,7 @@ StatusCode TauSubstructureVariables::execute(xAOD::TauJet& pTau) {
 
 	CaloClusterVariables CaloClusterVariablesTool;
 	CaloClusterVariablesTool.setVertexCorrection(m_doVertexCorrection);
+	CaloClusterVariablesTool.setIncSub(m_incShowerSubtr);
 
 	bool isFilled = CaloClusterVariablesTool.update(pTau);
 
@@ -132,8 +135,6 @@ StatusCode TauSubstructureVariables::execute(xAOD::TauJet& pTau) {
 	float dr(0.);
 
 	unsigned int num_clusters(0);
-	const xAOD::CaloCluster* incluster;
-	std::vector<const xAOD::CaloCluster*> vClusters;
 
 	TLorentzVector leadClusVec;
 	TLorentzVector subLeadClusVec;
@@ -141,18 +142,12 @@ StatusCode TauSubstructureVariables::execute(xAOD::TauJet& pTau) {
 	double clusELead = -1111.0;
 	double clusESubLead = -1111.0;
 
-	// loop over all clusters of the jet seed
-	xAOD::JetConstituentVector jcv = taujetseed->getConstituents();
-	xAOD::JetConstituentVector::const_iterator nav_it   = jcv.begin();
-	xAOD::JetConstituentVector::const_iterator nav_itE  = jcv.end();
-	for (; nav_it != nav_itE; ++nav_it) {
-		++num_clusters;
-
-		incluster = dynamic_cast<const xAOD::CaloCluster*>( (*nav_it)->rawConstituent() );
-		if (!incluster) continue;
+	// Loop through jets, get links to clusters
+	std::vector<const xAOD::CaloCluster*> vClusters;
+	ATH_CHECK(tauRecTools::GetJetClusterList(taujetseed, vClusters, m_incShowerSubtr));
 
-		// save all clusters of jet seed
-		vClusters.push_back(incluster);
+	for (auto incluster : vClusters){
+	        ++num_clusters;
 
 		// calc total energy
 		totalEnergy += incluster->e();
@@ -200,10 +195,10 @@ StatusCode TauSubstructureVariables::execute(xAOD::TauJet& pTau) {
 	if (clusESubLead > 0.) {
 	  approxSubstructure4Vec += subLeadClusVec;
 	  }
-	
+
 	// now sort cluster by energy
 	std::sort(vClusters.begin(), vClusters.end(), DefCaloClusterCompare());
-	
+
 	// determine energy sum of leading 2 and leading 3 clusters
 	float sum2LeadClusterE(0.);
 	float sum3LeadClusterE(0.);
@@ -276,7 +271,6 @@ StatusCode TauSubstructureVariables::execute(xAOD::TauJet& pTau) {
 	pTau.setDetail(xAOD::TauJetParameters::ChPiEMEOverCaloEME,	static_cast<float>(fChPIEMEOverCaloEME));
 	pTau.setDetail(xAOD::TauJetParameters::EMPOverTrkSysP,		static_cast<float>(fEMPOverTrkSysP));
 
-
 	// jvf and sumPtTrk are now a vector and the old run1-type jvf value is stored in the 0-th element
 	// sumPtTrk is calculated wrt Vertices
 
diff --git a/Reconstruction/tauRecTools/src/TauAxisSetter.cxx b/Reconstruction/tauRecTools/src/TauAxisSetter.cxx
index 726845458ce..e288b4e6fab 100644
--- a/Reconstruction/tauRecTools/src/TauAxisSetter.cxx
+++ b/Reconstruction/tauRecTools/src/TauAxisSetter.cxx
@@ -10,6 +10,7 @@
 #include "CaloUtils/CaloVertexedCluster.h"
 
 #include "TauAxisSetter.h"
+#include "tauRecTools/HelperFunctions.h"
 
 /********************************************************************/
 TauAxisSetter::TauAxisSetter(const std::string& name) :
@@ -89,21 +90,16 @@ StatusCode TauAxisSetter::execute(xAOD::TauJet& pTau)
     if(m_doVertexCorrection)
     {
 	TLorentzVector tauInterAxis;
-	
-	for (cItr = pJetSeed->getConstituents().begin(); cItr != cItrE; ++cItr) {
-	  tempClusterVector.SetPtEtaPhiE( (*cItr)->pt(), (*cItr)->eta(), (*cItr)->phi(), (*cItr)->e() );
-	  if (BaryCenter.DeltaR(tempClusterVector) > m_clusterCone)
-	    continue;
-	  
-	  const xAOD::CaloCluster* cluster = dynamic_cast<const xAOD::CaloCluster*>( (*cItr)->rawConstituent() ); 
-	  if (!cluster) continue;
-	  
+
+	std::vector<const xAOD::CaloCluster*> clusterList;
+	ATH_CHECK(tauRecTools::GetJetClusterList(pJetSeed, clusterList, m_incShowerSubtr, BaryCenter, m_clusterCone));
+	for (auto cluster : clusterList){
 	  if (pTau.vertexLink())
 	    tauInterAxis += xAOD::CaloVertexedCluster(*cluster, (*pTau.vertexLink())->position()).p4();
 	  else
 	    tauInterAxis += xAOD::CaloVertexedCluster(*cluster).p4();
 	}
-	
+
 	// save values for intermediate axis 
 	ATH_MSG_DEBUG("tau axis:" << tauInterAxis.Pt()<< " " << tauInterAxis.Eta() << " " << tauInterAxis.Phi()  << " " << tauInterAxis.E() );
         pTau.setP4(tauInterAxis.Pt(), tauInterAxis.Eta(), tauInterAxis.Phi(), pTau.m());
diff --git a/Reconstruction/tauRecTools/src/TauAxisSetter.h b/Reconstruction/tauRecTools/src/TauAxisSetter.h
index b2d5d9cb055..93ef0de5024 100644
--- a/Reconstruction/tauRecTools/src/TauAxisSetter.h
+++ b/Reconstruction/tauRecTools/src/TauAxisSetter.h
@@ -36,6 +36,7 @@ private:
 
     Gaudi::Property<double> m_clusterCone {this, "ClusterCone", 0.2, "cone of tau candidate"};
     Gaudi::Property<bool> m_doVertexCorrection {this, "VertexCorrection", true, "switch of vertex correction"};
+    Gaudi::Property<bool> m_incShowerSubtr {this, "IncShowerSubtr", true, "use shower subtracted clusters in calo calculations"};
 };
 
 #endif
diff --git a/Reconstruction/tauRecTools/src/TauPi0ClusterCreator.cxx b/Reconstruction/tauRecTools/src/TauPi0ClusterCreator.cxx
index 8cbfcfa25b0..13441c7a25c 100644
--- a/Reconstruction/tauRecTools/src/TauPi0ClusterCreator.cxx
+++ b/Reconstruction/tauRecTools/src/TauPi0ClusterCreator.cxx
@@ -16,7 +16,7 @@
 #include "xAODJet/Jet.h"
 
 #include "TauPi0ClusterCreator.h"
-
+#include "tauRecTools/HelperFunctions.h"
 
 using std::vector;
 using std::string;
@@ -421,18 +421,17 @@ bool TauPi0ClusterCreator::setHadronicClusterPFOs(xAOD::TauJet& pTau, xAOD::PFOC
         ATH_MSG_ERROR("Could not retrieve tau jet seed");
         return false;
     }
-    xAOD::JetConstituentVector::const_iterator clusterItr   = tauJetSeed->getConstituents().begin();
-    xAOD::JetConstituentVector::const_iterator clusterItrE  = tauJetSeed->getConstituents().end();
-    for (; clusterItr != clusterItrE; ++clusterItr){
+    std::vector<const xAOD::CaloCluster*> clusterList;
+
+    StatusCode sc = tauRecTools::GetJetClusterList(tauJetSeed, clusterList, m_incShowerSubtr);
+    if (!sc) return false;
+
+    for (auto cluster : clusterList){
         // Procedure: 
         // - Calculate cluster energy in Hcal. This is to treat -ve energy cells correctly
         // - Then set 4momentum via setP4(E/cosh(eta), eta, phi, m). This forces the PFO to have the correct energy and mass
         // - Ignore clusters outside 0.2 cone and those with overall negative energy or negative energy in Hcal
 
-        // Get xAOD::CaloClusters from jet constituent
-        const xAOD::CaloCluster* cluster = dynamic_cast<const xAOD::CaloCluster*>( (*clusterItr)->rawConstituent() );
-        if (!cluster) continue;
-
         // Don't create PFOs for clusters with overall (Ecal+Hcal) negative energy (noise)
         if(cluster->e()<=0.) continue;
 
diff --git a/Reconstruction/tauRecTools/src/TauPi0ClusterCreator.h b/Reconstruction/tauRecTools/src/TauPi0ClusterCreator.h
index 2c33180d00d..3e42bd24b5c 100644
--- a/Reconstruction/tauRecTools/src/TauPi0ClusterCreator.h
+++ b/Reconstruction/tauRecTools/src/TauPi0ClusterCreator.h
@@ -65,6 +65,8 @@ private:
     /** @brief pt threshold for pi0 candidate clusters */
     Gaudi::Property<double> m_clusterEtCut {this, "ClusterEtCut", 0.5 * Gaudi::Units::GeV, "Et threshould for pi0 candidate clusters"};
     
+    Gaudi::Property<bool> m_incShowerSubtr {this, "IncShowerSubtr", true, "use shower subtracted clusters in calo calculations"};
+
     SG::ReadHandleKey<xAOD::CaloClusterContainer> m_pi0ClusterInputContainer{this,"Key_Pi0ClusterContainer", "TauPi0SubtractedClusters", "input pi0 cluster"};
 
 };
diff --git a/Reconstruction/tauRecTools/tauRecTools/CaloClusterVariables.h b/Reconstruction/tauRecTools/tauRecTools/CaloClusterVariables.h
index c8294413bdc..32fc490e512 100644
--- a/Reconstruction/tauRecTools/tauRecTools/CaloClusterVariables.h
+++ b/Reconstruction/tauRecTools/tauRecTools/CaloClusterVariables.h
@@ -33,6 +33,7 @@ public:
     bool update(const xAOD::TauJet& pTau); //!< update the internal variables for the given tau
 
     void setVertexCorrection(bool flag) {m_doVertexCorrection=flag;}
+    void setIncSub(bool flag) {m_incShowerSubtr=flag;}
 
     // ID Variables
     unsigned int numConstituents() { return (unsigned int) m_numConstit; }
@@ -61,14 +62,16 @@ private:
     double m_totEnergy;
     double m_effEnergy;
 
-    /** Calculate the geometrical center of the tau constituents */
+    // Calculate the geometrical center of the tau constituents
     TLorentzVector calculateTauCentroid(int nConst, const std::vector<CaloVertexedClusterType>& constituents);
-    
-    /** 
-     * Enable cell origin correction.
-     * Eta and phi of the cells are corrected wrt to the origin of the tau vertex
-     */
+
+    // Enable cell origin correction.
+    // Eta and phi of the cells are corrected wrt to the origin of the tau vertex
     bool m_doVertexCorrection;
+
+    // use shower subtracted clusters with PFlow jet seeds
+    bool m_incShowerSubtr;
+
 };
 
 //-------------------------------------------------------------------------
diff --git a/Reconstruction/tauRecTools/tauRecTools/HelperFunctions.h b/Reconstruction/tauRecTools/tauRecTools/HelperFunctions.h
index d394f7c7b97..66121ed93af 100644
--- a/Reconstruction/tauRecTools/tauRecTools/HelperFunctions.h
+++ b/Reconstruction/tauRecTools/tauRecTools/HelperFunctions.h
@@ -22,7 +22,7 @@ namespace tauRecTools
 {
   ANA_MSG_HEADER(msgHelperFunction)
 
-  const StatusCode GetJetClusterList(const xAOD::Jet* jet, std::vector<const xAOD::CaloCluster*> &clusterList, bool incShowerSubtracted);
+    const StatusCode GetJetClusterList(const xAOD::Jet* jet, std::vector<const xAOD::CaloCluster*> &clusterList, bool incShowerSubtracted, TLorentzVector dRVector = TLorentzVector(0.,0.,0.,0.), double dRCut = -1);
 
   xAOD::TauTrack::TrackFlagType isolateClassifiedBits(xAOD::TauTrack::TrackFlagType flag);
   bool sortTracks(const ElementLink<xAOD::TauTrackContainer> &l1, const ElementLink<xAOD::TauTrackContainer> &l2);
diff --git a/Reconstruction/tauRecTools/tauRecTools/TauSubstructureVariables.h b/Reconstruction/tauRecTools/tauRecTools/TauSubstructureVariables.h
index a6038c22f3e..6a24ca76db0 100644
--- a/Reconstruction/tauRecTools/tauRecTools/TauSubstructureVariables.h
+++ b/Reconstruction/tauRecTools/tauRecTools/TauSubstructureVariables.h
@@ -31,17 +31,18 @@ class TauSubstructureVariables : public TauRecToolBase
         virtual StatusCode finalize() override;
 
     private:
-        /** Maximal pile up correction in GeV for a tau candidate.
-         *  Used for the caloIso corrected variable.
-         */
-        double m_maxPileUpCorrection; 
+        // Maximal pile up correction in GeV for a tau candidate.
+        // Used for the caloIso corrected variable.
+	double m_maxPileUpCorrection;
         double m_pileUpAlpha;         //!< slope of the pileup correction
         
-        /** 
-         * enable cell origin correction 
-         * eta and phi of the cells are corrected wrt to the origin of the tau vertex
-         */
-        bool m_doVertexCorrection;
+        // enable cell origin correction
+        // eta and phi of the cells are corrected wrt to the origin of the tau vertex
+	bool m_doVertexCorrection;
+
+	// use shower subtracted clusters with PFlow jet seeds
+	bool m_incShowerSubtr;
+
 };
 
 #endif
-- 
GitLab