Commit 28ecbd31 authored by Teng Jian Khoo's avatar Teng Jian Khoo Committed by Graeme Stewart
Browse files

update of WeightPFOTool interface (JetRecTools-00-03-50)

	* Root/PFlowPseudoJetGetter: Update to WeightPFOTool interface
	* Tag as JetRecTools-00-03-50

2016-09-10    <khoo@cern.ch>
	* Root/PFlowPseudoJetGetter: initial change to retain CPFOs with weight 0 as jet constituents
	* Tag as JetRecTools-00-03-49

2016-09-10    <delsart@lpsc1120x.in2p3.fr>

	* Root/LinkDef.h: added
	* Tag as JetRecTools-00-03-48

2016-08-31  scott snyder  <snyder@bnl.gov>

	* Tagging JetRecTools-00-03-47.
	* Comply with ATLAS naming conventions.

2016-08-23  scott snyder  <snyder@bnl.gov>

	* Tagging JetRecTools-00-03-46.
...
(Long ChangeLog diff - truncated)
parent 8c1d238c
......@@ -45,3 +45,4 @@ atlas_add_component( JetRecTools
# Install files from the package:
atlas_install_headers( JetRecTools )
atlas_install_python_modules( python/*.py )
......@@ -24,6 +24,7 @@
#include "JetRec/PseudoJetGetter.h"
#include "AsgTools/ToolHandle.h"
#include "PFlowUtils/IRetrievePFOTool.h"
#include "PFlowUtils/IWeightPFOTool.h"
namespace jet{
class TrackVertexAssociation;
......@@ -48,6 +49,7 @@ protected:
const LabelIndex* pli) const {return PseudoJetGetter::buildCUI(ppar,idx,pli);}
ToolHandle<CP::IRetrievePFOTool> m_retrievePFOTool; /// Retrieval tool
ToolHandle<CP::IWeightPFOTool> m_weightPFOTool; /// Retrieval tool
bool m_inputIsEM; /// If true EM clusters are used for neutral PFOs.
bool m_calibrate; /// If true, EM PFOs are calibrated to LC.
bool m_useneutral; // IF true, neutral pflow is included.
......
......@@ -26,8 +26,10 @@
// InputCollection: Name of the input cluster collection.
// GridSize: The grid size that should be applied for the SK
// algorithm. Suggested values between 0.3 and 0.6
// SKRapMin: The minimum rapidity over which to calculate and apply SK
// SKRapMax: The maximum rapidity over which to calculate and apply SK
// SKRapMin: The minimum (absolute) rapidity over which to calculate SK
// SKRapMax: The maximum (absolute) rapidity over which to calculate SK
// SKRapMinApplied: The minimum (absolute) rapidity over which to apply SK
// SKRapMaxApplied: The maximum (absolute) rapidity over which to apply SK
// isCaloSplit: If false, SK is run the same on all clusters. If
// true, SK is run separately for clusters in the ECal and the
// HCal.
......@@ -47,8 +49,6 @@
#include <string>
//#include "AsgTools/ToolHandle.h"
//#include "AsgTools/AsgTool.h"
#include "JetRecTools/JetConstituentModifierBase.h"
#include "xAODBase/IParticleContainer.h"
......@@ -95,6 +95,8 @@ private:
float m_hCalGrid;
float m_rapmin;
float m_rapmax;
float m_rapminApplied;
float m_rapmaxApplied;
mutable double m_minPt;
mutable double m_minPtECal;
mutable double m_minPtHCal;
......
......@@ -32,7 +32,7 @@ StatusCode CaloClusterConstituentsOrigin::process(xAOD::IParticleContainer* cont
const xAOD::VertexContainer* vertexContainer=0;
//get vertexcontainer from eventstore
CHECK( evtStore()->retrieve(vertexContainer, m_vertexContName) );
ATH_CHECK( evtStore()->retrieve(vertexContainer, m_vertexContName) );
return process(clust, vertexContainer->at(0));
}
......@@ -55,9 +55,11 @@ StatusCode CaloClusterConstituentsOrigin::processLC(xAOD::CaloClusterContainer*
for(xAOD::CaloCluster* cl : *cont) {
xAOD::CaloVertexedTopoCluster corrCL( *cl,vert->position());
cl->setEta(corrCL.eta());
cl->setPhi(corrCL.phi());
if(cl->calE()>1e-9) {
xAOD::CaloVertexedTopoCluster corrCL( *cl,vert->position());
cl->setEta(corrCL.eta());
cl->setPhi(corrCL.phi());
}
}
return StatusCode::SUCCESS;
}
......@@ -66,10 +68,12 @@ StatusCode CaloClusterConstituentsOrigin::processEM(xAOD::CaloClusterContainer*
for(xAOD::CaloCluster* cl : *cont) {
xAOD::CaloVertexedTopoCluster corrCL( *cl,xAOD::CaloCluster::UNCALIBRATED, vert->position());
cl->setE(corrCL.e());
cl->setEta(corrCL.eta());
cl->setPhi(corrCL.phi());
if(cl->rawE()>1e-9) {
xAOD::CaloVertexedTopoCluster corrCL( *cl,xAOD::CaloCluster::UNCALIBRATED, vert->position());
cl->setE(corrCL.e());
cl->setEta(corrCL.eta());
cl->setPhi(corrCL.phi());
}
}
return StatusCode::SUCCESS;
}
......
......@@ -30,9 +30,8 @@ JetConstituentModSequence::JetConstituentModSequence(const std::string &name): a
}
#ifdef ASGTOOL_ATHENA
StatusCode JetConstituentModSequence::initialize() {
CHECK( m_modifiers.retrieve() );
ATH_CHECK( m_modifiers.retrieve() );
if( m_modifiers.empty() ) {
ATH_MSG_ERROR(" empty container !!" );
return StatusCode::FAILURE;
......@@ -52,7 +51,7 @@ StatusCode JetConstituentModSequence::initialize() {
int JetConstituentModSequence::execute() const {
const xAOD::IParticleContainer* cont = 0;
CHECK( evtStore()->retrieve(cont, m_inputContainer) );
ATH_CHECK( evtStore()->retrieve(cont, m_inputContainer) );
xAOD::IParticleContainer* modifiedCont = 0;
......@@ -96,5 +95,4 @@ int JetConstituentModSequence::execute() const {
return 0;
}
#endif
/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#include <JetRecTools/JetTrackSelectionTool.h>
#include <JetRecTools/TrackVertexAssociationTool.h>
#include <JetRecTools/TrackPseudoJetGetter.h>
#ifdef __CINT__
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ nestedclass;
#pragma link C++ class JetTrackSelectionTool+;
#pragma link C++ class TrackVertexAssociationTool+;
#pragma link C++ class TrackPseudoJetGetter+;
#endif
......@@ -53,11 +53,12 @@ namespace PFlowPJHelper{
}
PFlowPseudoJetGetter::PFlowPseudoJetGetter(const std::string &name)
: PseudoJetGetter(name), m_retrievePFOTool("RetrievePFOTool"), m_trkVtxAssocName("JetTrackVtxAssoc") {
: PseudoJetGetter(name), m_retrievePFOTool("RetrievePFOTool"), m_weightPFOTool("WeightPFOTool"), m_trkVtxAssocName("JetTrackVtxAssoc") {
declareProperty("RetrievePFOTool", m_retrievePFOTool, "Name of tool that builds the PFO collection.");
declareProperty("WeightPFOTool", m_weightPFOTool, "Name of tool that extracts the cPFO weights.");
declareProperty("InputIsEM", m_inputIsEM =false, "True if neutral PFOs are EM scale clusters.");
declareProperty("CalibratePFO", m_calibrate =true, "True if LC calibration should be applied to EM PFOs.");
declareProperty("UseNeutral", m_useneutral =true, "True to use the nuetral component of PFlow.");
declareProperty("UseNeutral", m_useneutral =true, "True to use the neutral component of PFlow.");
declareProperty("UseCharged", m_usecharged =true, "True if use the charged component of PFlow.");
declareProperty("UseVertices", m_usevertices = true, "True if we make use of the primary vertex information.");
declareProperty("UseChargedWeights",m_useChargedWeights = true, "True if we make use of weighting scheme for charged PFO");
......@@ -131,6 +132,7 @@ int PFlowPseudoJetGetter::appendTo(PseudoJetVector& psjs, const LabelIndex* pli)
// Get the charged pflow.
if ( m_usecharged ) {
const xAOD::PFOContainer* pcpfs = m_retrievePFOTool->retrievePFO(CP::EM,CP::charged);
TLorentzVector cpfo_p4;
for ( const xAOD::PFO* pcpf : *pcpfs ) {
if ( pcpf == 0 ) {
ATH_MSG_WARNING("Have NULL pointer to charged PFO");
......@@ -141,6 +143,7 @@ int PFlowPseudoJetGetter::appendTo(PseudoJetVector& psjs, const LabelIndex* pli)
ATH_MSG_WARNING("Skipping charged PFO with null track pointer.");
continue;
}
cpfo_p4.Clear();
bool matchedToPrimaryVertex = false;
if (true == m_useTrackToVertexTool && true == m_usevertices){
......@@ -166,47 +169,15 @@ int PFlowPseudoJetGetter::appendTo(PseudoJetVector& psjs, const LabelIndex* pli)
}
if ( true == matchedToPrimaryVertex){
float weight = 0.0;
if (true == m_useChargedWeights) {
if (ptrk){
//This weight allows us to linearly de-weight higher pt tracks as we move towards the calo only regime
if (ptrk->pt() < 30000) weight = 1.0;
else if (ptrk->pt() >= 30000 && ptrk->pt() <= 60000) weight = (1.0 - (ptrk->pt()-30000)/30000);
if (0.0 != weight){
xAOD::PFODetails::PFOAttributes myAttribute_isInDenseEnvironment = xAOD::PFODetails::PFOAttributes::eflowRec_isInDenseEnvironment;
int isInDenseEnvironment = false;
bool gotVariable = pcpf->attribute(myAttribute_isInDenseEnvironment,isInDenseEnvironment);
if (false == gotVariable) ATH_MSG_WARNING("This charged PFO did not have eflowRec_isInDenseEnvironment set");
float expectedEnergy = 0.0;
xAOD::PFODetails::PFOAttributes myAttribute_tracksExpectedEnergyDeposit = xAOD::PFODetails::PFOAttributes::eflowRec_tracksExpectedEnergyDeposit;
gotVariable = pcpf->attribute(myAttribute_tracksExpectedEnergyDeposit,expectedEnergy);
if (false == gotVariable) {
ATH_MSG_WARNING("This charged PFO did not have eflowRec_tracksExpectedEnergyDeposit set");
weight = 1.0;
filler.fill(pcpf, pcpf->p4()*weight);
}
else{
if (true == isInDenseEnvironment) {
TLorentzVector modifiedFourVector;//we cannot directly modify the xAOD::PFO because it is const
modifiedFourVector.SetPtEtaPhiM((ptrk->e()-expectedEnergy)/cosh(pcpf->eta()),pcpf->eta(),pcpf->phi(),pcpf->m());
filler.fill(pcpf, modifiedFourVector*weight );
}
else{
float expectedPt = expectedEnergy/cosh(pcpf->eta());
if (1.0 == weight) filler.fill(pcpf, pcpf->p4()*weight );
else{
float secondWeight = (expectedPt + weight*(ptrk->pt()-expectedPt))/ptrk->pt();
filler.fill(pcpf, pcpf->p4()*secondWeight );
}
}
}
}
}
else ATH_MSG_WARNING("This charged PFO has no track attached to it");
float weight = 0.0;
ATH_CHECK( m_weightPFOTool->fillWeight( *pcpf, weight ) );
//if (weight>FLT_MIN){ // check against float precision
ATH_MSG_VERBOSE("Fill pseudojet for CPFO with weighted pt: " << pcpf->pt()*weight);
filler.fill(pcpf, pcpf->p4()*weight);
//} else {
// ATH_MSG_VERBOSE("CPFO had a weight of 0, do not fill.");
//} // received a weight
}//if should use charged PFO weighting scheme
else filler.fill(pcpf, pcpf->p4());
}
......
......@@ -27,12 +27,14 @@ using namespace std;
SoftKillerWeightTool::SoftKillerWeightTool(const std::string& name) : JetConstituentModifierBase(name)
, m_lambdaCalDivide(317)
, m_gridSpacing(0.6)
, m_rapmin(0.0)
, m_rapmax(2.5)
, m_isCaloSplit(false)
, m_gridSpacing(0.6)
, m_eCalGrid(0.5)
, m_hCalGrid(0.5)
, m_rapmin(0.0)
, m_rapmax(2.5)
, m_rapminApplied(0)
, m_rapmaxApplied(10)
{
......@@ -43,6 +45,8 @@ SoftKillerWeightTool::SoftKillerWeightTool(const std::string& name) : JetConstit
declareProperty("SKGridSize", m_gridSpacing = 0.6);
declareProperty("SKRapMin", m_rapmin = 0);
declareProperty("SKRapMax", m_rapmax = 2.5);
declareProperty("SKRapMinApplied", m_rapminApplied = 0);
declareProperty("SKRapMaxApplied", m_rapmaxApplied = 10);
declareProperty("isCaloSplit", m_isCaloSplit = false);
declareProperty("ECalGridSize", m_eCalGrid = 0.5);
declareProperty("HCalGridSize", m_hCalGrid = 0.5);
......@@ -75,74 +79,74 @@ StatusCode SoftKillerWeightTool::process(xAOD::IParticleContainer* cont) const {
// Finds the pT cut for this event based on the SK results
// The clustSK collection contains all clusters that aren't cut, so clusters below
// its min pT are cut
double SoftKillerWeightTool::findMinPt(vector<fastjet::PseudoJet> *m_clustSK) const {
double m_minPt = 999999999999;
double SoftKillerWeightTool::findMinPt(vector<fastjet::PseudoJet> *clustSK) const {
double minPt = 999999999999;
for(unsigned int i=0; i < m_clustSK->size(); i++){
if( (*m_clustSK)[i].pt() < m_minPt) m_minPt = (*m_clustSK)[i].pt();
for(unsigned int i=0; i < clustSK->size(); i++){
if( (*clustSK)[i].pt() < minPt) minPt = (*clustSK)[i].pt();
}
// There is a small rounding error which I account for this way
return (m_minPt - 1e-12);
return (minPt - 1e-12);
}
// Reweights clusters (when calo isn't split)
void SoftKillerWeightTool::RunClusters(xAOD::CaloClusterContainer* m_clust) const {
vector<fastjet::PseudoJet> m_clustPJ;
void SoftKillerWeightTool::RunClusters(xAOD::CaloClusterContainer* clust) const {
vector<fastjet::PseudoJet> clustPJ;
for(xAOD::CaloCluster* cl : *m_clust){
if((cl)->e() > 0) m_clustPJ.push_back( fastjet::PseudoJet( (cl)->p4() ));
for(xAOD::CaloCluster* cl : *clust){
if((cl)->e() > 0) clustPJ.push_back( fastjet::PseudoJet( (cl)->p4() ));
}
fastjet::Selector m_selector = fastjet::SelectorAbsRapRange(m_rapmin, m_rapmax);
fastjet::RectangularGrid m_SKgrid(-m_rapmax, m_rapmax, m_gridSpacing, m_gridSpacing, m_selector);
fastjet::contrib::SoftKiller m_softkiller(m_SKgrid);
std::vector<fastjet::PseudoJet> m_clustSK = m_softkiller(m_selector(m_clustPJ));
fastjet::Selector selector = fastjet::SelectorAbsRapRange(m_rapmin, m_rapmax);
fastjet::RectangularGrid SKgrid(-m_rapmax, m_rapmax, m_gridSpacing, m_gridSpacing, selector);
fastjet::contrib::SoftKiller softkiller(SKgrid);
std::vector<fastjet::PseudoJet> clustSK = softkiller(selector(clustPJ));
m_minPt = findMinPt(&m_clustSK);
m_minPt = findMinPt(&clustSK);
}
void SoftKillerWeightTool::RunSplitClusters(xAOD::CaloClusterContainer* m_clust) const {
vector<fastjet::PseudoJet> m_clustPJ_ECal;
vector<fastjet::PseudoJet> m_clustPJ_HCal;
void SoftKillerWeightTool::RunSplitClusters(xAOD::CaloClusterContainer* clust) const {
vector<fastjet::PseudoJet> clustPJ_ECal;
vector<fastjet::PseudoJet> clustPJ_HCal;
for(xAOD::CaloCluster* cl : *m_clust){
double m_center_lambda;
(cl)->retrieveMoment(xAOD::CaloCluster::CENTER_LAMBDA,m_center_lambda);
if( m_center_lambda < m_lambdaCalDivide && (cl)->e() > 0) m_clustPJ_ECal.push_back( fastjet::PseudoJet( (cl)->p4() ));
if( m_center_lambda >= m_lambdaCalDivide && (cl)->e() > 0) m_clustPJ_HCal.push_back( fastjet::PseudoJet( (cl)->p4() ));
for(xAOD::CaloCluster* cl : *clust){
double center_lambda;
(cl)->retrieveMoment(xAOD::CaloCluster::CENTER_LAMBDA,center_lambda);
if( center_lambda < m_lambdaCalDivide && (cl)->e() > 0) clustPJ_ECal.push_back( fastjet::PseudoJet( (cl)->p4() ));
if( center_lambda >= m_lambdaCalDivide && (cl)->e() > 0) clustPJ_HCal.push_back( fastjet::PseudoJet( (cl)->p4() ));
}
fastjet::Selector m_selector = fastjet::SelectorAbsRapRange(m_rapmin, m_rapmax);
fastjet::RectangularGrid m_SKgridECal(-m_rapmax, m_rapmax, m_eCalGrid, m_eCalGrid, m_selector);
fastjet::contrib::SoftKiller m_softkillerECal(m_SKgridECal);
std::vector<fastjet::PseudoJet> m_clustSK_ECal = m_softkillerECal(m_selector(m_clustPJ_ECal));
m_minPtECal = findMinPt(&m_clustSK_ECal);
fastjet::Selector selector = fastjet::SelectorAbsRapRange(m_rapmin, m_rapmax);
fastjet::RectangularGrid SKgridECal(-m_rapmax, m_rapmax, m_eCalGrid, m_eCalGrid, selector);
fastjet::contrib::SoftKiller softkillerECal(SKgridECal);
std::vector<fastjet::PseudoJet> clustSK_ECal = softkillerECal(selector(clustPJ_ECal));
m_minPtECal = findMinPt(&clustSK_ECal);
fastjet::RectangularGrid m_SKgridHCal(-m_rapmax, m_rapmax, m_hCalGrid, m_hCalGrid, m_selector);
fastjet::contrib::SoftKiller m_softkillerHCal(m_SKgridHCal);
std::vector<fastjet::PseudoJet> m_clustSK_HCal = m_softkillerHCal(m_selector(m_clustPJ_HCal));
m_minPtHCal = findMinPt(&m_clustSK_HCal);
fastjet::RectangularGrid SKgridHCal(-m_rapmax, m_rapmax, m_hCalGrid, m_hCalGrid, selector);
fastjet::contrib::SoftKiller softkillerHCal(SKgridHCal);
std::vector<fastjet::PseudoJet> clustSK_HCal = softkillerHCal(selector(clustPJ_HCal));
m_minPtHCal = findMinPt(&clustSK_HCal);
}
float SoftKillerWeightTool::calculateWeight(xAOD::CaloCluster* cl) const{
// If the cluster pT is below the SoftKiller pT cut, rescale 4-momentum to 0
if( abs(cl->eta()) < m_rapmin || abs(cl->eta()) > m_rapmax) return 1;
if( abs(cl->eta()) < m_rapminApplied || abs(cl->eta()) > m_rapmaxApplied) return 1;
if( cl->pt() < m_minPt) return 0;
return 1;
}
float SoftKillerWeightTool::calculateSplitWeight(xAOD::CaloCluster* cl) const{
if( abs(cl->eta()) < m_rapmin || abs(cl->eta()) > m_rapmax) return 1;
double m_center_lambda;
if( abs(cl->eta()) < m_rapminApplied || abs(cl->eta()) > m_rapmaxApplied) return 1;
double center_lambda;
if(!cl->retrieveMoment(xAOD::CaloCluster::CENTER_LAMBDA,m_center_lambda)) m_center_lambda = 0;
if(!cl->retrieveMoment(xAOD::CaloCluster::CENTER_LAMBDA,center_lambda)) center_lambda = 0;
//Make a separate pT cut for the ECal and HCal
if( m_center_lambda < m_lambdaCalDivide && cl->pt() < m_minPtECal) return 0;
if( center_lambda < m_lambdaCalDivide && cl->pt() < m_minPtECal) return 0;
if( cl->pt() < m_minPtHCal) return 0;
return 1;
}
......
# this file is -*- C++ -*-
# this makefile also gets parsed by shell scripts
# therefore it does not support full make syntax and features
# edit with care
......@@ -9,10 +10,10 @@ PACKAGE = JetRecTools
PACKAGE_PRELOAD =
PACKAGE_CXXFLAGS =
PACKAGE_OBJFLAGS =
PACKAGE_LDFLAGS =
PACKAGE_LDFLAGS = -lConstituentSubtractor -lSoftKiller
PACKAGE_BINFLAGS =
PACKAGE_LIBFLAGS =
PACKAGE_DEP = xAODEventInfo JetRec PFlowUtils InDetTrackSelectionTool TrackVertexAssociationTool
PACKAGE_DEP = xAODEventInfo xAODTruth JetRec PFlowUtils Asg_FastJet InDetTrackSelectionTool TrackVertexAssociationTool
PACKAGE_TRYDEP =
PACKAGE_CLEAN =
PACKAGE_NOGRID =
......
......@@ -39,10 +39,8 @@ use xAODBase xAODBase-* Event/xAOD
use xAODCaloEvent xAODCaloEvent-* Event/xAOD
use xAODTracking xAODTracking-* Event/xAOD
use AsgTools AsgTools-* Control/AthToolSupport
use PFlowUtils PFlowUtils-* Reconstruction/PFlow
private
# Add explicit dependency against xAODPFlow to avoid indirect dependency from PFlowUtils
# which creates a component library and cmake (correctly) does not propagate such dependencies
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment