Commit 2c110c74 authored by William Keaton Balunas's avatar William Keaton Balunas Committed by Frank Winklmeier
Browse files

Support FlowElements in JetReconstruction

parent 3da23c5c
......@@ -110,6 +110,7 @@ namespace xAOD {
{ "LCPFlow", LCPFlow },
{ "EMPFlow", EMPFlow },
{ "EMCPFlow", EMCPFlow },
{ "EMPFlowFE", EMPFlowFE },
{ "TrackCaloCluster", TrackCaloCluster },
{ "EMTopoOriginSK", EMTopoOriginSK },
{ "EMTopoOriginCS", EMTopoOriginCS },
......@@ -160,6 +161,7 @@ namespace xAOD {
{ LCPFlow, "LCPFlow" },
{ EMPFlow, "EMPFlow" },
{ EMCPFlow, "EMCPFlow" },
{ EMPFlowFE, "EMPFlowFE" },
{ TrackCaloCluster, "TrackCaloCluster" },
{ EMTopoOriginSK, "EMTopoOriginSK" },
{ EMTopoOriginCS, "EMTopoOriginCS" },
......
......@@ -95,6 +95,7 @@ namespace xAOD {
HI,
HIClusters,
Other = 100,
EMPFlowFE = 200, // Temporary, until xAOD::PFO is phased out and replaced with xAOD::FlowElement
Uncategorized= 1000
};
......
......@@ -48,6 +48,9 @@ namespace xAOD {
}
void FlowElement_v1::setP4(FourMom_t p4){
setP4(p4.Pt(), p4.Eta(), p4.Phi(), p4.M());
}
AUXSTORE_PRIMITIVE_SETTER_AND_GETTER(FlowElement_v1, float, charge, setCharge)
bool FlowElement_v1::isCharged() const { return !bool( signalType()& Neutral );}
......
......@@ -76,6 +76,7 @@ namespace xAOD {
virtual Type::ObjectType type() const override ;
void setP4(float pt, float eta, float phi, float m) ;
void setP4(FourMom_t p4);
///@}
// *************************************************
......
......@@ -65,7 +65,8 @@ StatusCode ThinNegativeEnergyNeutralPFOsAlg::initialize()
ATH_MSG_ERROR( "StreamName property has not been initialized." );
return StatusCode::FAILURE;
}
ATH_CHECK( m_neutralPFOsKey.initialize (m_streamName, m_doThinning) );
ATH_CHECK( m_neutralPFOsKey.initialize (m_streamName, m_doThinning && !m_neutralPFOsKey.key().empty()) );
ATH_CHECK( m_neutralPFOsFEKey.initialize (m_streamName, m_doThinning && !m_neutralPFOsFEKey.key().empty()) );
// Initialize the counters to zero
m_nEventsProcessed = 0;
......@@ -101,25 +102,43 @@ StatusCode ThinNegativeEnergyNeutralPFOsAlg::execute()
const EventContext& ctx = Gaudi::Hive::currentContext();
// Retrieve the container
SG::ThinningHandle<xAOD::PFOContainer> neutralPFOs (m_neutralPFOsKey, ctx);
// Set up masks
std::vector<bool> mask;
int nNeutralPFOs = neutralPFOs->size();
m_nNeutralPFOsProcessed += nNeutralPFOs;
mask.assign(nNeutralPFOs,false);
// Loop over NeutralPFOs and update mask
for (int i=0; i<nNeutralPFOs; ++i) {
const xAOD::PFO* neutralPFO = (*neutralPFOs)[i];
// Retain postive energy neutral PFOs
if (neutralPFO->ptEM()>0.0) {mask[i] = true;}
else {++m_nNeutralPFOsThinned;}
if(!m_neutralPFOsFEKey.key().empty()){
// Retrieve the container
SG::ThinningHandle<xAOD::PFOContainer> neutralPFOs (m_neutralPFOsKey, ctx);
// Set up masks
std::vector<bool> mask;
int nNeutralPFOs = neutralPFOs->size();
m_nNeutralPFOsProcessed += nNeutralPFOs;
mask.assign(nNeutralPFOs,false);
// Loop over NeutralPFOs and update mask
for (int i=0; i<nNeutralPFOs; ++i) {
const xAOD::PFO* neutralPFO = (*neutralPFOs)[i];
// Retain postive energy neutral PFOs
if (neutralPFO->ptEM()>0.0) {mask[i] = true;}
else {++m_nNeutralPFOsThinned;}
}
// Apply masks to thinning service
neutralPFOs.keep (mask);
}
// Apply masks to thinning service
neutralPFOs.keep (mask);
if(!m_neutralPFOsFEKey.key().empty()){
SG::ThinningHandle<xAOD::FlowElementContainer> neutralFEs (m_neutralPFOsFEKey, ctx);
std::vector<bool> mask;
int nNeutralFEs = neutralFEs->size();
mask.assign(nNeutralFEs, false);
for(int i=0; i<nNeutralFEs; i++){
const xAOD::FlowElement* neutralFE = (*neutralFEs)[i];
// TODO: Is this OK for LC-scale PFOs?
// Can't access EM-scale momentum without some link to EM-scale FlowElements.
if(neutralFE->pt() > 0.0) mask[i] = true;
}
neutralFEs.keep(mask);
}
return StatusCode::SUCCESS;
}
......
......@@ -23,6 +23,7 @@
#include "StoreGate/ThinningHandleKey.h"
#include "xAODPFlow/PFOContainer.h"
#include "xAODPFlow/FlowElementContainer.h"
class ThinNegativeEnergyNeutralPFOsAlg
: public ::AthAlgorithm
......@@ -54,7 +55,9 @@ private:
/// Names of the containers to thin
SG::ThinningHandleKey<xAOD::PFOContainer> m_neutralPFOsKey
{ this, "NeutralPFOsKey", "JetETMissNeutralParticleFlowObjects", "StoreGate key for the PFOContainer to be thinned" };
{ this, "NeutralPFOsKey", "JetETMissNeutralParticleFlowObjects", "StoreGate key for the PFOContainer to be thinned (if any)" };
SG::ThinningHandleKey<xAOD::FlowElementContainer> m_neutralPFOsFEKey
{ this, "NeutralPFOsFEKey", "", "StoreGate key for the FlowElementContainer to be thinned (if any)" };
/// Counters
unsigned long m_nEventsProcessed;
......
......@@ -68,6 +68,8 @@ ContainersOnTheFly = [
["AntiKt4EMTopoCSSKJetsAux","xAOD::JetAuxContainer"],
["AntiKt4EMPFlowCSSKJets","xAOD::JetContainer"],
["AntiKt4EMPFlowCSSKJetsAux","xAOD::JetAuxContainer"],
["AntiKt4EMPFlowFEJets","xAOD::JetContainer"],
["AntiKt4EMPFlowFEJetsAux","xAOD::JetContainerAux"],
["AntiKt4EMTopoJets_BTagging201810","xAOD::JetContainer"],
["AntiKt4EMTopoJets_BTagging201810Aux","xAOD::ShallowAuxContainer"],
["AntiKt4EMPFlowJets_BTagging201903","xAOD::JetContainer"],
......
......@@ -46,6 +46,7 @@ FullListOfSmartContainers = [
"AntiKt10LCTopoTrimmedPtFrac5SmallR20ExKt2GASubJets",
"AntiKt10LCTopoTrimmedPtFrac5SmallR20ExKt3GASubJets",
"AntiKt4EMPFlowJets",
"AntiKt4EMPFlowFEJets",
"AntiKt4EMPFlowJets_BTagging201903",
"AntiKt4EMPFlowJets_BTagging201810",
"AntiKt4EMTopoJets_BTagging201810",
......
......@@ -499,8 +499,10 @@ class SlimmingHelper:
items.extend(AntiKt10UFOCSSKRecursiveSoftDropBeta100Zcut5NinfJetsCPContent)
elif collectionName=="AntiKt4EMPFlowJets":
from DerivationFrameworkJetEtMiss.AntiKt4EMPFlowJetsCPContent import AntiKt4EMPFlowJetsCPContent
#from DerivationFrameworkCore.AntiKt4EMPFlowJetsCPContent import AntiKt4EMPFlowJetsCPContent
items.extend(AntiKt4EMPFlowJetsCPContent)
elif collectionName=="AntiKt4EMPFlowFEJets":
from DerivationFrameworkJetEtMiss.AntiKt4EMPFlowFEJetsCPContent import AntiKt4EMPFlowFEJetsCPContent
items.extend(AntiKt4EMPFlowFEJetsCPContent)
elif collectionName=="AntiKt4EMPFlowJets_BTagging201810":
if "AntiKt4EMPFlowJets_BTagging201810" not in self.AppendToDictionary:
self.AppendToDictionary["AntiKt4EMPFlowJets_BTagging201810"]='xAOD::JetContainer'
......
# Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
AntiKt4EMPFlowFEJetsCPContent = [
"Kt4EMPFlowFEEventShape",
"Kt4EMPFlowFEEventShapeAux.Density",
"AntiKt4EMPFlowFEJets",
"AntiKt4EMPFlowFEJetsAux.pt.eta.phi.m.JetConstitScaleMomentum_pt.JetConstitScaleMomentum_eta.JetConstitScaleMomentum_phi.JetConstitScaleMomentum_m.NumTrkPt500.SumPtTrkPt500.NumChargedPFOPt500.SumPtChargedPFOPt500.EnergyPerSampling.ActiveArea4vec_eta.ActiveArea4vec_m.ActiveArea4vec_phi.ActiveArea4vec_pt.DetectorEta.DetectorY.FracSamplingMax.FracSamplingMaxIndex.GhostTrack.Jvt.JVFCorr.JvtRpt.NumTrkPt1000.NumChargedPFOPt1000.TrackWidthPt1000.ChargedPFOWidthPt1000.GhostMuonSegmentCount.PartonTruthLabelID.ConeTruthLabelID.HadronConeExclExtendedTruthLabelID.HadronConeExclTruthLabelID.TrueFlavor.Timing",
"PrimaryVertices",
"PrimaryVerticesAux.vertexType"
]
......@@ -122,6 +122,10 @@ def addAntiKt10TruthWZJets(sequence,outputlist):
if DerivationFrameworkIsMonteCarlo:
addStandardJets("AntiKt", 1.0, "TruthWZ", ptmin=40000, mods="truth_ungroomed", algseq=sequence, outputGroup=outputlist)
def addAntiKt4EMPFlowJetsFE(sequence, outputlist):
addCHSPFlowObjectsFE()
addStandardJets("AntiKt", 0.4, "EMPFlowFE", ptmin=10000, ptminFilter=15000, mods="pflow_ungroomed", algseq=sequence, outputGroup=outputlist)
##################################################################
def replaceAODReducedJets(jetlist,sequence,outputlist):
......@@ -395,7 +399,7 @@ def addJetPtAssociation(jetalg, truthjetalg, sequence, algname):
if hasattr(ToolSvc,jetptassociationtoolname):
jetaugtool.JetPtAssociationTool = getattr(ToolSvc,jetptassociationtoolname)
else:
jetptassociationtool = CfgMgr.JetPtAssociationTool(jetptassociationtoolname, InputContainer=truthjetalg, AssociationName="GhostTruth")
jetptassociationtool = CfgMgr.JetPtAssociationTool(jetptassociationtoolname, JetContainer=jetalg+'Jets', MatchingJetContainer=truthjetalg, AssociationName="GhostTruth")
ToolSvc += jetptassociationtool
jetaugtool.JetPtAssociationTool = jetptassociationtool
......@@ -480,7 +484,7 @@ def addQGTaggerTool(jetalg, sequence, algname, truthjetalg=None ):
if hasattr(ToolSvc,jetptassociationtoolname):
jetaugtool.JetPtAssociationTool = getattr(ToolSvc,jetptassociationtoolname)
else:
jetptassociationtool = CfgMgr.JetPtAssociationTool(jetptassociationtoolname, InputContainer=truthjetalg, AssociationName="GhostTruth")
jetptassociationtool = CfgMgr.JetPtAssociationTool(jetptassociationtoolname, JetContainer=jetalg+'Jets', MatchingJetContainer=truthjetalg, AssociationName="GhostTruth")
ToolSvc += jetptassociationtool
jetaugtool.JetPtAssociationTool = jetptassociationtool
......@@ -693,6 +697,8 @@ def addConstModJets(jetalg,radius,inputtype,constmods,sequence,outputlist,custom
constit = JetConstit( xAODType.CaloCluster, ["LC","Origin"])
elif inputtype == "EMPFlow":
constit = JetConstit( xAODType.ParticleFlow )
elif inputtype == "EMPFlowFE":
constit = JetConstit( xAODType.FlowElement )
constit.modifiers += constmods
......@@ -773,6 +779,26 @@ def addCHSPFlowObjects():
job.jetalg.Tools.append(jtm.jetconstitCHSPFlow)
extjetlog.info("Added CHS PFlow sequence to \'jetalg\'")
extjetlog.info(job.jetalg.Tools)
def addCHSPFlowObjectsFE():
# Only act if the collection does not already exist
from RecExConfig.AutoConfiguration import IsInInputFile
if not IsInInputFile("xAOD::FlowElementContainer","CHSFlowElements"):
# Check that an alg doing this has not already been inserted
from AthenaCommon.AlgSequence import AlgSequence
job = AlgSequence()
from JetRec.JetRecStandard import jtm
if not hasattr(job,"jetalgCHSPFlowFE") and not hasattr(jtm,"jetconstitCHSPFlowFE"):
from JetRec.JetRecConf import JetToolRunner
jtm += JetToolRunner("jetconstitCHSPFlowFE",
EventShapeTools=[],
Tools=[jtm.JetConstitSeq_PFlowCHS_FE])
# Add this tool runner to the JetAlgorithm instance "jetalg"
# which runs all preparatory tools
# This was added by JetCommon
job.jetalg.Tools.append(jtm.jetconstitCHSPFlowFE)
extjetlog.info("Added CHS PFlow (FlowElement) sequence to \'jetalg\'")
extjetlog.info(job.jetalg.Tools)
##################################################################
applyJetCalibration_xAODColl("AntiKt4EMTopo")
updateJVT_xAODColl("AntiKt4EMTopo")
......
......@@ -143,11 +143,12 @@ def reCreatePseudoJets(jetalg, rsize, inputtype, variableRMassScale=-1.0, variab
return [jtm.tools[tmpName]]
# no container exist. simply build a new one.
if inputtype=="LCTopo" or inputtype=="EMTopo" or inputtype == "EMPFlow" or inputtype == "EMCPFlow":
if inputtype=="LCTopo" or inputtype=="EMTopo" or inputtype == "EMPFlow" or inputtype == "EMCPFlow" or inputtype == "EMPFlowFE":
defaultmods = {"EMTopo":"emtopo_ungroomed",
"LCTopo":"lctopo_ungroomed",
"EMPFlow":"pflow_ungroomed",
"EMCPFlow":"pflow_ungroomed",
"EMPFlowFE":"pflow_ungroomed",
"Truth":"truth_ungroomed",
"TruthWZ":"truth_ungroomed",
"PV0Track":"track_ungroomed"}
......@@ -191,6 +192,8 @@ def reCreatePseudoJets(jetalg, rsize, inputtype, variableRMassScale=-1.0, variab
constit = JetConstit( xAODType.CaloCluster, ["LC","Origin"])
elif inputtype == "EMPFlow":
constit = JetConstit( xAODType.ParticleFlow )
elif inputtype == "EMPFlowFE":
constit = JetConstit( xAODType.FlowElement )
constit.modifiers += constmods
......@@ -441,6 +444,7 @@ def addStandardJets(jetalg, rsize, inputtype, ptmin=0., ptminFilter=0.,
"LCTopo":"lctopo_ungroomed",
"EMPFlow":"pflow_ungroomed",
"EMCPFlow":"pflow_ungroomed",
"EMPFlowFE":"pflow_ungroomed",
"Truth":"truth_ungroomed",
"TruthWZ":"truth_ungroomed",
"PV0Track":"track_ungroomed",
......@@ -459,7 +463,7 @@ def addStandardJets(jetalg, rsize, inputtype, ptmin=0., ptminFilter=0.,
finderArgs['overwrite']=True
# map the input to the jtm code for PseudoJetGetter
getterMap = dict( LCTopo = 'lctopo', EMTopo = 'emtopo', EMPFlow = 'empflow', EMCPFlow = 'emcpflow',
getterMap = dict( LCTopo = 'lctopo', EMTopo = 'emtopo', EMPFlow = 'empflow', EMCPFlow = 'emcpflow', EMPFlowFE = 'empflowfe',
Truth = 'truth', TruthWZ = 'truthwz', TruthDressedWZ = 'truthdressedwz', TruthCharged = 'truthcharged',
PV0Track='pv0track')
# create the finder for the temporary collection.
......@@ -483,7 +487,7 @@ def addStandardJets(jetalg, rsize, inputtype, ptmin=0., ptminFilter=0.,
alg = JetAlgorithm(algname, Tools = pretools+[finderTool])
dfjetlog.info( "Added "+algname+" to sequence "+algseq.name() )
algseq += alg
DFJetAlgs[algname] = alg;
DFJetAlgs[algname] = alg
##################################################################
# Schedule the adding of BCID info
......
......@@ -30,6 +30,7 @@ public:
protected:
void fillEperSamplingCluster(const xAOD::Jet &jet, std::vector<float> & ePerSampling ) const ;
void fillEperSamplingPFO(const xAOD::Jet &jet, std::vector<float> & ePerSampling ) const ;
void fillEperSamplingFE(const xAOD::Jet &jet, std::vector<float> & ePerSampling ) const ;
private:
Gaudi::Property<std::string> m_jetContainerName{this, "JetContainer", "", "SG key for the input jet container"};
......
......@@ -65,7 +65,8 @@ private: // data
/// Properties.
Gaudi::Property<std::string> m_aname{this, "AssociationName", "", "Key for association vector"};
SG::ReadHandleKey<xAOD::JetContainer> m_jetContainerName{this, "JetContainer", "", "Input jet container"};
Gaudi::Property<std::string> m_jetContainerName{this, "JetContainer", "", "Input jet container"};
SG::ReadHandleKey<xAOD::JetContainer> m_matchingContainerName{this, "MatchingJetContainer", "", "Jet container to match to"};
SG::WriteDecorHandleKey<xAOD::JetContainer> m_assocFracKey{this, "AssociationFractionName", "AssociationFraction", "SG key for output AssociationFraction decoration"};
SG::WriteDecorHandleKey<xAOD::JetContainer> m_assocLinkKey{this, "AssociationLinkName", "AssociationLink", "SG key for output AssociationLink decoration"};
......
......@@ -7,9 +7,11 @@
#include "xAODJet/JetAccessorMap.h"
#include "JetUtils/JetCaloQualityUtils.h"
#include "PFlowUtils/FEHelpers.h"
#include "xAODCaloEvent/CaloCluster.h"
#include "xAODPFlow/PFO.h"
#include "xAODPFlow/FlowElement.h"
#include <vector>
......@@ -62,13 +64,17 @@ StatusCode JetCaloEnergies::decorate(const xAOD::JetContainer& jets) const {
if ( ctype == xAOD::Type::CaloCluster ) {
ATH_MSG_VERBOSE(" Constituents are calo clusters.");
fillEperSamplingCluster(*jet, ePerSampling);
} else if (ctype == xAOD::Type::ParticleFlow) {
ATH_MSG_VERBOSE(" Constituents are pflow objects.");
fillEperSamplingPFO(*jet, ePerSampling);
} else if (ctype == xAOD::Type::FlowElement) {
ATH_MSG_VERBOSE(" Constituents are FlowElements.");
fillEperSamplingFE(*jet, ePerSampling);
}else {
ATH_MSG_VERBOSE("Constituents are not calo clusters nor pflow objects.");
ATH_MSG_VERBOSE("Constituents are not CaloClusters, PFOs, or FlowElements.");
}
}
......@@ -97,10 +103,10 @@ void JetCaloEnergies::fillEperSamplingCluster(const xAOD::Jet& jet, std::vector<
psFracHandle(jet) = jet::JetCaloQualityUtils::presamplerFraction( &jet );
}
#define FillESamplingPFO( LAYERNAME ) \
float E_##LAYERNAME = 0.0; \
#define FillESamplingPFO( LAYERNAME ) \
float E_##LAYERNAME = 0.0; \
if (constit->attribute(xAOD::PFODetails::eflowRec_LAYERENERGY_##LAYERNAME, E_##LAYERNAME)) { \
ePerSampling[CaloSampling::LAYERNAME] += E_##LAYERNAME; \
ePerSampling[CaloSampling::LAYERNAME] += E_##LAYERNAME; \
}
void JetCaloEnergies::fillEperSamplingPFO(const xAOD::Jet & jet, std::vector<float> & ePerSampling ) const {
......@@ -114,50 +120,50 @@ void JetCaloEnergies::fillEperSamplingPFO(const xAOD::Jet & jet, std::vector<flo
if (jet.rawConstituent(i)->type()==xAOD::Type::ParticleFlow){
const xAOD::PFO* constit = static_cast<const xAOD::PFO*>(jet.rawConstituent(i));
if ( fabs(constit->charge())>FLT_MIN ){
eTot += constit->track(0)->e();
eTot += constit->track(0)->e();
} else {
eTot += constit->e();
FillESamplingPFO(PreSamplerB);
FillESamplingPFO(EMB1);
FillESamplingPFO(EMB2);
FillESamplingPFO(EMB3);
FillESamplingPFO(PreSamplerE);
FillESamplingPFO(EME1);
FillESamplingPFO(EME2);
FillESamplingPFO(EME3);
FillESamplingPFO(HEC0);
FillESamplingPFO(HEC1);
FillESamplingPFO(HEC2);
FillESamplingPFO(HEC3);
FillESamplingPFO(TileBar0);
FillESamplingPFO(TileBar1);
FillESamplingPFO(TileBar2);
FillESamplingPFO(TileGap1);
FillESamplingPFO(TileGap2);
FillESamplingPFO(TileGap3);
FillESamplingPFO(TileExt0);
FillESamplingPFO(TileExt1);
FillESamplingPFO(TileExt2);
FillESamplingPFO(FCAL0);
FillESamplingPFO(FCAL1);
FillESamplingPFO(FCAL2);
FillESamplingPFO(MINIFCAL0);
FillESamplingPFO(MINIFCAL1);
FillESamplingPFO(MINIFCAL2);
FillESamplingPFO(MINIFCAL3);
emTot += ( E_PreSamplerB+E_EMB1+E_EMB2+E_EMB3+
E_PreSamplerE+E_EME1+E_EME2+E_EME3+
E_FCAL0 );
hecTot += ( E_HEC0+E_HEC1+E_HEC2+E_HEC3 );
eTot += constit->e();
FillESamplingPFO(PreSamplerB);
FillESamplingPFO(EMB1);
FillESamplingPFO(EMB2);
FillESamplingPFO(EMB3);
FillESamplingPFO(PreSamplerE);
FillESamplingPFO(EME1);
FillESamplingPFO(EME2);
FillESamplingPFO(EME3);
FillESamplingPFO(HEC0);
FillESamplingPFO(HEC1);
FillESamplingPFO(HEC2);
FillESamplingPFO(HEC3);
FillESamplingPFO(TileBar0);
FillESamplingPFO(TileBar1);
FillESamplingPFO(TileBar2);
FillESamplingPFO(TileGap1);
FillESamplingPFO(TileGap2);
FillESamplingPFO(TileGap3);
FillESamplingPFO(TileExt0);
FillESamplingPFO(TileExt1);
FillESamplingPFO(TileExt2);
FillESamplingPFO(FCAL0);
FillESamplingPFO(FCAL1);
FillESamplingPFO(FCAL2);
FillESamplingPFO(MINIFCAL0);
FillESamplingPFO(MINIFCAL1);
FillESamplingPFO(MINIFCAL2);
FillESamplingPFO(MINIFCAL3);
emTot += ( E_PreSamplerB+E_EMB1+E_EMB2+E_EMB3+
E_PreSamplerE+E_EME1+E_EME2+E_EME3+
E_FCAL0 );
hecTot += ( E_HEC0+E_HEC1+E_HEC2+E_HEC3 );
}//only consider neutral PFO
} else {
ATH_MSG_WARNING("Tried to call fillEperSamplingPFlow with a jet constituent that is not a PFO!");
......@@ -186,4 +192,44 @@ void JetCaloEnergies::fillEperSamplingPFO(const xAOD::Jet & jet, std::vector<flo
}
//**********************************************************************
void JetCaloEnergies::fillEperSamplingFE(const xAOD::Jet& jet, std::vector<float> & ePerSampling ) const {
float emTot = 0.;
float hecTot = 0.;
float psTot = 0.;
float eTot = 0.;
size_t numConstit = jet.numConstituents();
for ( size_t i=0; i<numConstit; i++ ) {
if(jet.rawConstituent(i)->type()!=xAOD::Type::FlowElement) {
ATH_MSG_WARNING("Tried to call fillEperSamplingFE with a jet constituent that is not a FlowElement!");
continue;
}
const xAOD::FlowElement* constit = static_cast<const xAOD::FlowElement*>(jet.rawConstituent(i));
if(constit->isCharged()){
eTot += constit->chargedObject(0)->e();
}
else{
eTot += constit->e();
std::vector<float> constitEPerSampling = FEHelpers::getEnergiesPerSampling(*constit);
for ( size_t s = CaloSampling::PreSamplerB; s < CaloSampling::Unknown; s++ ) {
ePerSampling[s] += constitEPerSampling[s];
}
emTot += (constitEPerSampling[CaloSampling::PreSamplerB] + constitEPerSampling[CaloSampling::EMB1]
+ constitEPerSampling[CaloSampling::EMB2] + constitEPerSampling[CaloSampling::EMB3]
+ constitEPerSampling[CaloSampling::PreSamplerE] + constitEPerSampling[CaloSampling::EME1]
+ constitEPerSampling[CaloSampling::EME2] + constitEPerSampling[CaloSampling::EME3]
+ constitEPerSampling[CaloSampling::FCAL0]);
hecTot += (constitEPerSampling[CaloSampling::HEC0] + constitEPerSampling[CaloSampling::HEC1]
+ constitEPerSampling[CaloSampling::HEC2] + constitEPerSampling[CaloSampling::HEC3]);
psTot += (constitEPerSampling[CaloSampling::PreSamplerB] + constitEPerSampling[CaloSampling::PreSamplerE]);
}
}
SG::WriteDecorHandle<xAOD::JetContainer, float> emFracHandle(m_emFracKey);
SG::WriteDecorHandle<xAOD::JetContainer, float> hecFracHandle(m_hecFracKey);
SG::WriteDecorHandle<xAOD::JetContainer, float> psFracHandle(m_psFracKey);
emFracHandle(jet) = eTot != 0. ? emTot/eTot : 0.;
hecFracHandle(jet) = eTot != 0. ? hecTot/eTot : 0.;
psFracHandle(jet) = eTot != 0. ? psTot/eTot : 0.;
}
\ No newline at end of file
......@@ -8,11 +8,11 @@
#include "JetMomentTools/JetECPSFractionTool.h"
#include "xAODCaloEvent/CaloCluster.h"
#include "xAODPFlow/PFO.h"
#include "xAODPFlow/FlowElement.h"
#include "CaloGeoHelpers/CaloSampling.h"
using xAOD::IParticle;
using xAOD::CaloCluster;
using xAOD::PFO;
using xAOD::JetConstituent;
using xAOD::JetFourMom_t;
......@@ -62,11 +62,16 @@ double JetECPSFractionTool::energyFraction(const xAOD::Jet& jet) const {
continue;
}
const IParticle* ppar = pcon->rawConstituent();
if(ppar == nullptr){
ATH_MSG_WARNING("Constituent has no associated raw constituent.");
continue;
}
double conEnergy = pcon->e();
// If constituent is pflow, then look to see if it is cluster-based.
const PFO* ppfl = dynamic_cast<const PFO*>(ppar);
const CaloCluster* pclu = nullptr;
if ( ppfl != nullptr ) {
if ( ppar->type() == xAOD::Type::ParticleFlow){
const xAOD::PFO* ppfl = static_cast<const xAOD::PFO*>(ppar);
if ( !ppfl->isCharged() ) {
pclu = ppfl->cluster(0); // Assume PFO has exactly one cluster.
if ( pclu != nullptr ) ATH_MSG_VERBOSE(" Constituent is a PFO pointing to a cluster");
......@@ -74,8 +79,22 @@ double JetECPSFractionTool::energyFraction(const xAOD::Jet& jet) const {
} else {
ATH_MSG_VERBOSE(" Constituent is a PFO pointing to a track");
}
// Otherwise check if constituent is cluster.
} else {
}
else if ( ppar->type() == xAOD::Type::FlowElement){
const xAOD::FlowElement* pfe = static_cast<const xAOD::FlowElement*>(ppar);
if(!(pfe->signalType() & xAOD::FlowElement::PFlow)){
ATH_MSG_WARNING(" Constituent is a FlowElement but not a PFO. This isn't supported; skipping.");
continue;
}
if ( !pfe->isCharged() ) {
pclu = dynamic_cast<const CaloCluster*>(pfe->otherObject(0)); // Assume PFO has exactly one cluster.
if ( pclu != nullptr ) ATH_MSG_VERBOSE(" Constituent is a PFO pointing to a cluster");
else ATH_MSG_WARNING(" Constituent is a non-charge PFO but no cluster is found.");
}
else ATH_MSG_VERBOSE(" Constituent is a charged PFO");
}
else {
// Not a PFO. Check if constituent is a cluster.
pclu = dynamic_cast<const CaloCluster*>(ppar);
if ( pclu != nullptr ) ATH_MSG_VERBOSE(" Constituent is a cluster");
}
......
......@@ -31,12 +31,14 @@ StatusCode JetPtAssociationTool::initialize() {
ATH_MSG_ERROR("JetPtAssociationTool needs to have its input jet container configured!");
return StatusCode::FAILURE;
}
m_assocFracKey = m_jetContainerName.key() + "." + m_aname + m_assocFracKey.key();
m_assocLinkKey = m_jetContainerName.key() + "." + m_aname + m_assocLinkKey.key();
m_assocFracKey = m_jetContainerName + "." + m_aname + m_assocFracKey.key();
m_assocLinkKey = m_jetContainerName + "." + m_aname + m_assocLinkKey.key();
ATH_CHECK(m_assocFracKey.initialize());
ATH_CHECK(m_assocLinkKey.initialize());
ATH_CHECK(m_matchingContainerName.initialize());
return StatusCode::SUCCESS;