diff --git a/PhysicsAnalysis/HeavyIonPhys/HIEventUtils/Root/HITowerWeightTool.cxx b/PhysicsAnalysis/HeavyIonPhys/HIEventUtils/Root/HITowerWeightTool.cxx
index 1159f05437de616cf56207a48263a842429b9773..d61455cbb1b77f3d281cfa18c45e22298a918a55 100644
--- a/PhysicsAnalysis/HeavyIonPhys/HIEventUtils/Root/HITowerWeightTool.cxx
+++ b/PhysicsAnalysis/HeavyIonPhys/HIEventUtils/Root/HITowerWeightTool.cxx
@@ -85,8 +85,17 @@ StatusCode HITowerWeightTool::configureEvent()
auto itr=m_runMap.find(run_number);
if(itr==m_runMap.end())
{
- m_runIndex=0;
- ATH_MSG_WARNING("No calibration avaliable for " << run_number << ". Doing no eta-phi correction.");
+ //trying generic run number <=> no run dependence
+ run_number = 226000;
+ auto itrg=m_runMap.find(run_number);
+ if(itrg==m_runMap.end()){
+ m_runIndex=0;
+ ATH_MSG_WARNING("No generic calibration or calibration for " << run_number << " is avaliable. Doing no eta-phi correction.");
+ }
+ else {
+ m_runIndex=itrg->second;
+ ATH_MSG_DEBUG("Using generic calibration for eta-phi correction.");
+ }
}
else m_runIndex=itr->second;
m_runNumber=run_number;
@@ -161,7 +170,9 @@ StatusCode HITowerWeightTool::initialize()
ATH_MSG_FATAL("Cannot find TH3F h1_run_index in config file " << full_path );
return StatusCode::FAILURE;
}
- for(int xbin=1; xbin<h1_run_index->GetNbinsX(); xbin++) m_runMap.emplace_hint(m_runMap.end(),std::make_pair(h1_run_index->GetBinContent(xbin),xbin));
+ for(int xbin=1; xbin<=h1_run_index->GetNbinsX(); xbin++) {
+ m_runMap.emplace_hint(m_runMap.end(),std::make_pair(h1_run_index->GetBinContent(xbin),xbin));
+ }
f->Close();
m_init=true;
return StatusCode::SUCCESS;
diff --git a/Reconstruction/HeavyIonRec/HIGlobal/src/HIEventShapeFillerTool.cxx b/Reconstruction/HeavyIonRec/HIGlobal/src/HIEventShapeFillerTool.cxx
index 226f59b3aae88193b75b5c149aa812e6621f56d0..afb7b5a186e819f4bdca8166c5a56686ca686fcb 100644
--- a/Reconstruction/HeavyIonRec/HIGlobal/src/HIEventShapeFillerTool.cxx
+++ b/Reconstruction/HeavyIonRec/HIGlobal/src/HIEventShapeFillerTool.cxx
@@ -29,7 +29,6 @@ StatusCode HIEventShapeFillerTool::initializeCollection(xAOD::HIEventShapeContai
{
//change m_evtShape to m_evtShape
m_evtShape=evtShape_;
-
//tool is initialized only once
if(!m_index)
{
@@ -105,6 +104,10 @@ StatusCode HIEventShapeFillerTool::fillCollectionFromClusterContainer(const xAOD
weight_vector->reserve(theClusters->size());
SG::AuxElement::Decorator< float > decorator("HIEtaPhiWeight");
+ std::unique_ptr<std::vector<float> > cm_vector(new std::vector<float>());
+ cm_vector->reserve(theClusters->size());
+ SG::AuxElement::Decorator< float > cm_decorator("HIMag");
+
constexpr float area_cluster=HI::TowerBins::getBinArea();
if(m_towerWeightTool) CHECK(m_towerWeightTool->configureEvent());
@@ -125,6 +128,25 @@ StatusCode HIEventShapeFillerTool::fillCollectionFromClusterContainer(const xAOD
weight_vector->push_back(weight);
decorator(*cl)=weight;
+ //HIMag back in rel 22 (removed by mistake in 21)
+ float etot2=0;
+ float er2=0;
+
+ for(unsigned int sample=0; sample<24; sample++)
+ {
+ CaloSampling::CaloSample s=static_cast<CaloSampling::CaloSample>(sample);
+ if(!cl->hasSampling(s)) continue;
+ float esamp=std::abs(cl->eSample(s));
+ float w1=m_towerWeightTool->getWeight(eta,phi,s);
+ float wr=m_towerWeightTool->getWeightMag(eta,phi,s);
+ etot2+=esamp*w1;
+ er2+=esamp*wr;
+
+ }
+ float cm=er2/etot2;
+ cm_vector->push_back(cm);
+ cm_decorator(*cl)=cm;
+
//update members
slice->setEt(slice->et()+weight*ET);
slice->setRho(slice->rho() + weight*ET/area_cluster);
@@ -158,7 +180,6 @@ StatusCode HIEventShapeFillerTool::fillCollectionFromClusterContainer(const xAOD
StatusCode HIEventShapeFillerTool::fillCollectionFromCells(const SG::ReadHandleKey<CaloCellContainer> &cell_container_key)
{
- ATH_MSG_INFO("INSIDE FillCollectionFromCells");
//retrieve the cell container from store
SG::ReadHandle<CaloCellContainer> read_handle_caloCell ( cell_container_key );
const CaloCellContainer* CellContainer = read_handle_caloCell.get();
@@ -166,7 +187,6 @@ StatusCode HIEventShapeFillerTool::fillCollectionFromCells(const SG::ReadHandleK
}
StatusCode HIEventShapeFillerTool::fillCollectionFromCellContainer(const CaloCellContainer* CellContainer)
{
- ATH_MSG_INFO("INSIDE FillCollectionFromCellContainer");
//loop on Cells
for(const auto cellItr : *CellContainer) updateShape(m_evtShape,m_index,cellItr,1.,cellItr->eta(),cellItr->phi());
return StatusCode::SUCCESS;
diff --git a/Reconstruction/HeavyIonRec/HIJetRec/HIJetRec/HIJetConstituentModifierTool.h b/Reconstruction/HeavyIonRec/HIJetRec/HIJetRec/HIJetConstituentModifierTool.h
new file mode 100644
index 0000000000000000000000000000000000000000..f9e5f03e740d0e0ae11d25aa046edfe0c91572b5
--- /dev/null
+++ b/Reconstruction/HeavyIonRec/HIJetRec/HIJetRec/HIJetConstituentModifierTool.h
@@ -0,0 +1,47 @@
+// This file is -*- C++ -*-
+/*
+ Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+ */
+
+#ifndef HIJETREC_HIJETCONSTITUENTMODIFIERTOOL_H
+#define HIJETREC_HIJETCONSTITUENTMODIFIERTOOL_H
+ ///////////////////////////////////////////////////////
+ /// \class HIJetConstituentModifier
+ /// \author R.Longo - riccardo.longo@cern.ch
+ /// \date May 2020
+ ///
+ /// Tool to subtract final HIEventShape from constituents
+ /// on the fly (w/o modifying clusters)
+ ///////////////////////////////////////////////////////
+
+#include "JetRec/JetModifierBase.h"
+#include "xAODCaloEvent/CaloClusterContainer.h"
+#include "HIJetRec/IHISubtractorTool.h"
+
+#include "StoreGate/ReadHandleKey.h"
+#include "StoreGate/WriteHandleKey.h"
+
+class HIJetConstituentModifierTool : public JetModifierBase {
+
+ ASG_TOOL_CLASS(HIJetConstituentModifierTool, IJetModifier);
+
+ public:
+
+ // Constructor from tool name.
+ HIJetConstituentModifierTool(const std::string& myname);
+
+ virtual StatusCode initialize() override;
+ // Inherited method to modify a jet.
+ virtual int modifyJet(xAOD::Jet& jet) const;
+
+ private:
+ /// \brief Name of input cluster container
+ SG::ReadHandleKey< xAOD::CaloClusterContainer > m_clusterKey { this, "ClusterKey", "ClusterKey", "Name of the input Cluster Container"};
+ /// \brief handle to IHISubtractorTool that determines the subtracted kinematics for each constituent
+ ToolHandle<IHISubtractorTool> m_subtractorTool { this, "Subtractor", "HIJetSubtractorToolBase", "" };
+ /// |brief boolean switch to drive the JetScale settings after constituents are added
+ Gaudi::Property< bool > m_originCorrection { this, "ApplyOriginCorrection", false, "Apply Origin Correction boolean switch"};
+
+};
+
+#endif
diff --git a/Reconstruction/HeavyIonRec/HIJetRec/HIJetRec/HIJetConstituentSubtractionTool.h b/Reconstruction/HeavyIonRec/HIJetRec/HIJetRec/HIJetConstituentSubtractionTool.h
index 0e346b76c3e072c8dfef596e0832e7ac93db8f6b..30699de280eae07e8df4b351c698860e20c2ca5f 100644
--- a/Reconstruction/HeavyIonRec/HIJetRec/HIJetRec/HIJetConstituentSubtractionTool.h
+++ b/Reconstruction/HeavyIonRec/HIJetRec/HIJetRec/HIJetConstituentSubtractionTool.h
@@ -62,8 +62,6 @@ private:
/// \brief Name of HIEventShapeContainer
SG::ReadHandleKey< xAOD::HIEventShapeContainer > m_eventShapeKey { this, "EventShapeKey", "", "The input HI Event Shape"};
- SG::ReadHandleKey< xAOD::VertexContainer > m_vertexContainer { this, "VertexContainer", "PrimaryVertices", "Vertex container for primary vertices"};
-
//That looks useless. commented out for the moment
//std::string m_modulation_key;
/// \brief Subtracted jet kinematics are stored
diff --git a/Reconstruction/HeavyIonRec/HIJetRec/HIJetRec/HIJetRecDefs.h b/Reconstruction/HeavyIonRec/HIJetRec/HIJetRec/HIJetRecDefs.h
index 81a9cd5d92a82f10b7e22db4d1441fb740f7e38c..21d109e6484328eedcf64b719ad7ba72cc83c658 100644
--- a/Reconstruction/HeavyIonRec/HIJetRec/HIJetRec/HIJetRecDefs.h
+++ b/Reconstruction/HeavyIonRec/HIJetRec/HIJetRec/HIJetRecDefs.h
@@ -9,17 +9,22 @@
#include "xAODJet/JetTypes.h"
#include "HIEventUtils/HIEventDefs.h"
#include "FourMomUtils/xAODP4Helpers.h"
+#include "TLorentzVector.h"
namespace HIJetRec{
//define conventions for HIJets in terms of various signal states
- constexpr const char* unsubtractedJetState() { return "JetUnsubtractedScaleMomentum";}
- constexpr const char* subtractedJetState() { return "JetSubtractedScaleMomentum";}
- constexpr xAOD::JetConstitScale subtractedConstitState() {return xAOD::UncalibratedJetConstituent;}
+ constexpr const char* unsubtractedJetState() { return "JetUnsubtractedScaleMomentum";}
+ constexpr const char* subtractedJetState() { return "JetSubtractedScaleMomentum";}
+ constexpr const char* subtractedOriginCorrectedJetState() { return "JetSubtractedOriginCorrectedScaleMomentum";}
+
+ constexpr xAOD::JetConstitScale subtractedConstitState() { return xAOD::UncalibratedJetConstituent;}
+ constexpr xAOD::JetConstitScale subtractedOriginCorrectedConstitState() { return xAOD::CalibratedJetConstituent;}
constexpr xAOD::CaloCluster::State unsubtractedClusterState() {return xAOD::CaloCluster::ALTCALIBRATED;}
constexpr xAOD::CaloCluster::State subtractedClusterState() {return xAOD::CaloCluster::UNCALIBRATED;}
+ constexpr xAOD::CaloCluster::State subtractedOriginCorrectedClusterState() {return xAOD::CaloCluster::CALIBRATED;}
inline bool inTowerBoundary(float eta0, float phi0, float eta, float phi)
{
@@ -27,13 +32,13 @@ namespace HIJetRec{
if( 2.*std::abs(xAOD::P4Helpers::deltaPhi(phi0,phi)) > HI::TowerBins::getBinSizePhi() ) return false;
return true;
}
-
+
inline void setClusterP4(const xAOD::CaloCluster::FourMom_t& p, xAOD::CaloCluster* cl, xAOD::CaloCluster::State s)
{
float energy=p.Energy();
float eta=p.Eta();
float phi=p.Phi();
- if(energy < 0.)
+ if(energy < 0.)
{
eta*=-1.;
if(phi > 0.) phi-=M_PI;
@@ -62,5 +67,23 @@ namespace HIJetRec{
}
}
+/// \brief : getClusterP4 Added during Rel22 migration. We should never use CaloCluster_v1::p4() because this four vector can
+/// never have a negative energy. Therefore, if one tries to use setClusterP4 and then retrieve the cluster four momenta,
+/// must use this new method that keeps the cluster information as set in setClusterP4 (negative E and flipped eta phi)
+ inline TLorentzVector getClusterP4( const xAOD::CaloCluster* cl, xAOD::CaloCluster::State s ){
+ TLorentzVector correct_clusP4;
+ float energy=cl->e(s);
+ float eta=cl->eta(s);
+ float phi=cl->phi(s);
+ if(energy < 0.)
+ {
+ eta*=-1.;
+ if(phi > 0.) phi-=M_PI;
+ else phi+=M_PI;
+ }
+ float pt=energy/std::cosh(eta);
+ correct_clusP4.SetPtEtaPhiE(pt, eta,phi ,energy);
+ return correct_clusP4;
+ }
}
#endif
diff --git a/Reconstruction/HeavyIonRec/HIJetRec/Root/HIEventShapeJetIteration.cxx b/Reconstruction/HeavyIonRec/HIJetRec/Root/HIEventShapeJetIteration.cxx
index 1dcb08a83f76920d6c3096c37d54dd19d9696442..555e6036eff2215d5ccea43cc13bfcfa447ef053 100644
--- a/Reconstruction/HeavyIonRec/HIJetRec/Root/HIEventShapeJetIteration.cxx
+++ b/Reconstruction/HeavyIonRec/HIJetRec/Root/HIEventShapeJetIteration.cxx
@@ -55,6 +55,8 @@ int HIEventShapeJetIteration::execute() const
const HIEventShapeIndex* es_index = HIEventShapeMap::getIndex(m_inputEventShapeKey.key());
//New implementation after moving away from mutable
+ ATH_MSG_INFO("HIEventShapeJetIteration: found index for " << m_inputEventShapeKey.key());
+
if(es_index==nullptr)
{
ATH_MSG_INFO("No HIEventShapeIndex w/ name " << m_inputEventShapeKey.key() << " adding it to the map");
@@ -64,7 +66,9 @@ int HIEventShapeJetIteration::execute() const
}
const HIEventShapeIndex* other_index = HIEventShapeMap::getIndex(m_outputEventShapeKey.key());
- if(!other_index) HIEventShapeMap::insert( m_inputEventShapeKey.key(), *es_index );
+ if(!other_index) {
+ HIEventShapeMap::getMap()->insert( m_outputEventShapeKey.key(), *es_index );
+ }
//End of new implementation
const xAOD::JetContainer* theCaloJets=0;
diff --git a/Reconstruction/HeavyIonRec/HIJetRec/Root/HIJetConstituentModifierTool.cxx b/Reconstruction/HeavyIonRec/HIJetRec/Root/HIJetConstituentModifierTool.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..7230177af05af31871d3a574b77ba357eab772ca
--- /dev/null
+++ b/Reconstruction/HeavyIonRec/HIJetRec/Root/HIJetConstituentModifierTool.cxx
@@ -0,0 +1,112 @@
+/*
+ Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "HIJetRec/HIJetConstituentModifierTool.h"
+#include "HIJetRec/HIJetRecDefs.h"
+#include "xAODJet/JetConstituentVector.h"
+#include "xAODCaloEvent/CaloClusterContainer.h"
+#include "xAODBase/IParticle.h"
+#include "xAODBase/IParticleContainer.h"
+#include "JetRec/JetModifierBase.h"
+
+#include "StoreGate/ReadHandle.h"
+
+HIJetConstituentModifierTool::HIJetConstituentModifierTool(const std::string& myname): JetModifierBase(myname)
+{
+}
+
+StatusCode HIJetConstituentModifierTool::initialize(){
+
+ ATH_MSG_INFO("Initializing HIJetConstituentModifierTool w/ Cluster Key " << m_clusterKey.key() );
+ ATH_CHECK( m_clusterKey.initialize() );
+ return StatusCode::SUCCESS;
+
+}
+
+int HIJetConstituentModifierTool::modifyJet(xAOD::Jet& jet) const {
+
+ float E_min=m_subtractorTool->minEnergyForMoments();
+ const xAOD::JetConstituentVector constituents = jet.getConstituents();
+ std::vector<size_t> cluster_indices;
+ cluster_indices.reserve(constituents.size());
+
+ for (auto citer = constituents.begin(); citer != constituents.end(); ++citer)
+ {
+ cluster_indices.push_back(citer->rawConstituent()->index());
+ }
+
+ /// The accessor for the cluster element links
+ static SG::AuxElement::Accessor< std::vector< ElementLink< xAOD::IParticleContainer > > >
+ constituentAcc( "constituentLinks" );
+ static SG::AuxElement::Accessor< std::vector< float> >
+ constituentWeightAcc( "constituentWeights" );
+
+ if( constituentAcc.isAvailable(jet) ) constituentAcc( jet ).resize(0);
+ if( constituentWeightAcc.isAvailable(jet) ) constituentWeightAcc( jet ).resize(0);
+
+ //save unsubtracted kinematics as moment if they don’t exist already...
+ xAOD::IParticle::FourMom_t unsubtractedP4;
+ unsubtractedP4 = jet.p4();
+ jet.setJetP4(HIJetRec::unsubtractedJetState(),jet.jetP4());
+
+ xAOD::IParticle::FourMom_t subtractedP4;
+ xAOD::IParticle::FourMom_t subtractedOriginCorrP4;
+ xAOD::JetFourMom_t subtractedJet4vec;
+ xAOD::JetFourMom_t subtractedOriginCorrJet4vec;
+ //Cluster container access (from the e-gamma + cluster subtracted deep copy)
+ SG::ReadHandle<xAOD::CaloClusterContainer> readHandleSubtractedClusters ( m_clusterKey );
+
+ const xAOD::CaloClusterContainer* ccl=readHandleSubtractedClusters.get();
+
+ for(auto index : cluster_indices)
+ {
+ const xAOD::CaloCluster* cl=static_cast<const xAOD::CaloCluster*>(ccl->at(index));
+
+ jet.addConstituent(cl);
+ xAOD::IParticle::FourMom_t subtractedClusterP4;
+ subtractedClusterP4=HIJetRec::getClusterP4( cl, HIJetRec::subtractedClusterState());
+ subtractedP4+=subtractedClusterP4;
+ ATH_MSG_DEBUG("Subracted Cluster #: " << cl->index() << " :: E: " << subtractedClusterP4.E() << " :: Eta: " << subtractedClusterP4.Eta() << " :: Phi: " << subtractedClusterP4.Phi());
+
+ if( m_originCorrection ){
+ xAOD::IParticle::FourMom_t subtractedOriginCorrClusterP4;
+ subtractedOriginCorrClusterP4=HIJetRec::getClusterP4( cl, HIJetRec::subtractedOriginCorrectedClusterState());
+ subtractedOriginCorrP4+=subtractedOriginCorrClusterP4;
+ ATH_MSG_DEBUG("Subrtracted OC Cluster #: " << cl->index() << " :: E: " << subtractedOriginCorrClusterP4.E() << " :: Eta: " << subtractedOriginCorrClusterP4.Eta() << " :: Phi: " << subtractedOriginCorrClusterP4.Phi());
+ }
+
+ }
+
+ //Check for subtracted Kinematics
+ if(subtractedP4.E()/std::cosh(subtractedP4.Eta()) < E_min)
+ {
+ subtractedP4=unsubtractedP4;
+ subtractedP4*=1e-7;//ghost scale
+ }
+ //Check for subtracted + Origin Correction Kinematics
+ if(subtractedOriginCorrP4.E()/std::cosh(subtractedOriginCorrP4.Eta()) < E_min && m_originCorrection )
+ {
+ subtractedOriginCorrP4=unsubtractedP4;
+ subtractedOriginCorrP4*=1e-7;//ghost scale
+ }
+
+ subtractedJet4vec.SetCoordinates(subtractedP4.Pt(),subtractedP4.Eta(),
+ subtractedP4.Phi(),subtractedP4.M());
+ jet.setJetP4(HIJetRec::subtractedJetState(), subtractedJet4vec);
+
+ if( m_originCorrection ) {
+ subtractedOriginCorrJet4vec.SetCoordinates(subtractedOriginCorrP4.Pt(),subtractedOriginCorrP4.Eta(),
+ subtractedOriginCorrP4.Phi(),subtractedOriginCorrP4.M());
+ jet.setJetP4(HIJetRec::subtractedOriginCorrectedJetState(), subtractedOriginCorrJet4vec);
+ jet.setJetP4(subtractedOriginCorrJet4vec);
+ jet.setConstituentsSignalState(HIJetRec::subtractedOriginCorrectedConstitState());
+
+ }
+ else {
+ jet.setJetP4(subtractedJet4vec);
+ jet.setConstituentsSignalState(HIJetRec::subtractedConstitState());
+ }
+
+ return 0;
+}
diff --git a/Reconstruction/HeavyIonRec/HIJetRec/Root/HIJetConstituentSubtractionTool.cxx b/Reconstruction/HeavyIonRec/HIJetRec/Root/HIJetConstituentSubtractionTool.cxx
index c4f4202b7606767c7b14e12d8d09e34bcf59928b..f8dc7966ce4030afcede4eb0934907909665d31e 100644
--- a/Reconstruction/HeavyIonRec/HIJetRec/Root/HIJetConstituentSubtractionTool.cxx
+++ b/Reconstruction/HeavyIonRec/HIJetRec/Root/HIJetConstituentSubtractionTool.cxx
@@ -25,7 +25,6 @@ StatusCode HIJetConstituentSubtractionTool::initialize()
{
ATH_MSG_VERBOSE("HIJetConstituentSubtractionTool initialize");
ATH_CHECK( m_eventShapeKey.initialize( !m_eventShapeKey.key().empty()) );
- ATH_CHECK( m_vertexContainer.initialize( !m_vertexContainer.key().empty()) );
return StatusCode::SUCCESS;
}
@@ -39,20 +38,18 @@ StatusCode HIJetConstituentSubtractionTool::modify(xAOD::JetContainer& jets) con
//retrieve UE
//Introduction of a read handle for the HIShapeContainer
- SG::ReadHandle<xAOD::HIEventShapeContainer> read_handle_evtShape ( m_eventShapeKey );
+ SG::ReadHandle<xAOD::HIEventShapeContainer> readHandleEvtShape ( m_eventShapeKey );
const xAOD::HIEventShapeContainer* shape=nullptr;
const HIEventShapeIndex* es_index=nullptr;
if(m_eventShapeKey.key().compare("") != 0)
{
- //if(evtStore()->retrieve(shape,EventShapeKey()).isFailure())
- if (!read_handle_evtShape.isValid())
+ if (!readHandleEvtShape.isValid())
{
ATH_MSG_ERROR("Could not retrieve input HIEventShape " << m_eventShapeKey.key() );
return StatusCode::FAILURE;
}
- shape = read_handle_evtShape.get(); /** TODO check if this is neeed **/
- //HIEventShapeMap is a c++ map that neeeds the string key to identify the paired object
+ shape = readHandleEvtShape.get();
es_index=HIEventShapeMap::getIndex(m_eventShapeKey.key());
if(es_index==nullptr)
{
@@ -69,36 +66,6 @@ StatusCode HIJetConstituentSubtractionTool::modify(xAOD::JetContainer& jets) con
return StatusCode::FAILURE;
}
- const xAOD::Vertex* primVertex=nullptr;
- const xAOD::VertexContainer* vertices=nullptr;
-
- //Introduction of a read handle for the HIShapeContainer
- SG::ReadHandle<xAOD::VertexContainer> read_handle_vertexContainer ( m_vertexContainer );
-
- if(m_originCorrection)
- {
- if(!read_handle_vertexContainer.isValid())
- {
- ATH_MSG_ERROR("Could not retrieve VertexContainer " << m_vertexContainer.key());
- return StatusCode::FAILURE;
- }
- //if(evtStore()->retrieve(vertices,m_vertexContainer).isFailure())
- vertices = read_handle_vertexContainer.get();
- for ( size_t iVertex = 0; iVertex < vertices->size(); ++iVertex )
- {
- if(vertices->at(iVertex)->vertexType() == xAOD::VxType::PriVtx)
- {
- primVertex=vertices->at(iVertex);
- break;
- }
- }
- if(!primVertex)
- {
- ATH_MSG_WARNING("No primary vertices found, using first in container");
- primVertex=vertices->at(0);
- }
- }
-
const xAOD::HIEventShape* eshape = nullptr;
if(m_modulatorTool->getShape(eshape).isFailure())
{
@@ -106,72 +73,22 @@ StatusCode HIJetConstituentSubtractionTool::modify(xAOD::JetContainer& jets) con
return StatusCode::FAILURE;
}
- bool missingMoment=false;
- bool needsUnsubMoment=false;
- if(jets.size() > 0){
- xAOD::JetFourMom_t tmp;
- needsUnsubMoment = !((*jets.begin())->getAttribute<xAOD::JetFourMom_t>(HIJetRec::unsubtractedJetState(),tmp));
- }
-
- //check to see if unsubtracted moment has been stored
for ( xAOD::JetContainer::iterator ijet=jets.begin(); ijet!=jets.end(); ++ijet)
{
xAOD::IParticle::FourMom_t p4_cl;
xAOD::IParticle::FourMom_t p4_subtr;
xAOD::IParticle::FourMom_t p4_unsubtr;
- const xAOD::Vertex* origin=nullptr;
- if(m_originCorrection)
- {
- if( !(*ijet)->getAssociatedObject<xAOD::Vertex>("OriginVertex", origin) )
- {
- origin=primVertex;
- ATH_MSG_DEBUG("Jet has no associated vertex, using PV from container");
- }
- }
-
const xAOD::JetConstituentVector constituents = (*ijet)->getConstituents();
for (xAOD::JetConstituentVector::iterator itr = constituents.begin(); itr != constituents.end(); ++itr)
{
m_subtractorTool->subtract(p4_cl,itr->rawConstituent(),shape,es_index,m_modulatorTool, eshape); //modifies p4_cl to be constituent 4-vector AFTER subtraction
- if(m_originCorrection)
- {
- const xAOD::CaloCluster* cl=static_cast<const xAOD::CaloCluster*>(itr->rawConstituent());
- float mag = 0;
- if(cl->isAvailable<float>("HIMag")) mag=cl->auxdataConst<float>("HIMag");
- else
- {
- double cm_mag=0;
- if(cl->retrieveMoment (xAOD::CaloCluster::CENTER_MAG, cm_mag)) mag=cm_mag;
- }
- if(mag!=0.)
- {
- float eta0=cl->eta0();
- float phi0=cl->phi0();
- float radius=mag/std::cosh(eta0);
- xAOD::IParticle::FourMom_t p4_pos;
- p4_pos.SetX(radius*std::cos(phi0)-origin->x());
- p4_pos.SetY(radius*std::sin(phi0)-origin->y());
- p4_pos.SetZ(radius*std::sinh(eta0)-origin->z());
-
- double deta=p4_pos.Eta()-eta0;
- double dphi=p4_pos.Phi()-phi0;
- //adjust in case eta/phi are flipped in case of neg E clusters
- //this method is agnostic wrt convention
- if(p4_cl.Eta()*eta0 <0.) deta*=-1;
-
- double eta_prime=p4_cl.Eta()+deta;
- double phi_prime=p4_cl.Phi()+dphi;
- double e_subtr=p4_cl.E();
- p4_cl.SetPtEtaPhiE(e_subtr/std::cosh(eta_prime),eta_prime,phi_prime,e_subtr);
- }
- else missingMoment=true;
- }
-
p4_subtr+=p4_cl;
+
if( msgLvl(MSG::DEBUG) )
{
const xAOD::CaloCluster* cl=static_cast<const xAOD::CaloCluster*>(itr->rawConstituent());
+ //here we can still keep cl->p4 because it's taking the unsubtracted state - moreover is debug
p4_unsubtr+=cl->p4(HIJetRec::unsubtractedClusterState());
}
}
@@ -190,9 +107,6 @@ StatusCode HIJetConstituentSubtractionTool::modify(xAOD::JetContainer& jets) con
<< std::setw(10) << std::setprecision(3) << p4_subtr.E()*1e-3
<< std::setw(10) << std::setprecision(3) << p4_subtr.M()*1e-3);
-
-
-
xAOD::JetFourMom_t jet4vec;
//if entire jet has negative E, do no subtraction but set to ghost scale
//prevents cases with large cancellations with small E but pT non-trivial
@@ -203,14 +117,8 @@ StatusCode HIJetConstituentSubtractionTool::modify(xAOD::JetContainer& jets) con
}
jet4vec.SetCoordinates(p4_subtr.Pt(),p4_subtr.Eta(),p4_subtr.Phi(),p4_subtr.M());
-
(*ijet)->setJetP4(momentName(),jet4vec);
- xAOD::JetFourMom_t tmp;
- //if(! (*ijet)->getAttribute<xAOD::JetFourMom_t>(HIJetRec::unsubtractedJetState(),tmp) ){
- if(needsUnsubMoment)
- (*ijet)->setJetP4(HIJetRec::unsubtractedJetState(), (*ijet)->jetP4());
-// }
if(!momentOnly())
{
//hack for now to allow use of pp calib tool skipping pileup subtraction
@@ -218,11 +126,9 @@ StatusCode HIJetConstituentSubtractionTool::modify(xAOD::JetContainer& jets) con
(*ijet)->setJetP4("JetPileupScaleMomentum", jet4vec );
(*ijet)->setJetP4(xAOD::JetEMScaleMomentum, jet4vec);
- (*ijet)->setJetP4(jet4vec);
(*ijet)->setConstituentsSignalState(HIJetRec::subtractedConstitState());
}
}
- if(missingMoment) ATH_MSG_WARNING("No origin correction applied, CENTERMAG missing");
- //Fix from conflict beetween d1493284 (master) and 5af8a733 (21.0)
+
return StatusCode::SUCCESS;
}
diff --git a/Reconstruction/HeavyIonRec/HIJetRec/python/HIJetRecFlags.py b/Reconstruction/HeavyIonRec/HIJetRec/python/HIJetRecFlags.py
index d7a8b43435cc722fc23ed00df820781bb701b629..b3aeaa257c946d3c6460802566a4a211f9b2cce1 100644
--- a/Reconstruction/HeavyIonRec/HIJetRec/python/HIJetRecFlags.py
+++ b/Reconstruction/HeavyIonRec/HIJetRec/python/HIJetRecFlags.py
@@ -100,9 +100,9 @@ class AntiKtRValues(JobProperty):
class DoCellBasedSubtraction(JobProperty):
""" option to use cell based subtraction
"""
- statusOn = True
+ statusOn = False
allowedTypes = ['bool']
- StoredValue = True
+ StoredValue = False
class HarmonicsForSubtraction(JobProperty):
""" List of flow harmonics applied to jet subtraction
diff --git a/Reconstruction/HeavyIonRec/HIJetRec/python/HIJetRecUtils.py b/Reconstruction/HeavyIonRec/HIJetRec/python/HIJetRecUtils.py
index 345e86b517136d9b219598d4b44d2e9bd3f7b4a1..f5b8cc7c53cbc55b18d29b4c6994e7130dafbe9f 100644
--- a/Reconstruction/HeavyIonRec/HIJetRec/python/HIJetRecUtils.py
+++ b/Reconstruction/HeavyIonRec/HIJetRec/python/HIJetRecUtils.py
@@ -145,6 +145,8 @@ def ApplySubtractionToClusters(**kwargs) :
if 'cluster_key' in kwargs.keys() : cluster_key=kwargs['cluster_key']
else : cluster_key=HIJetFlags.HIClusterKey()
+ if 'output_cluster_key' in kwargs.keys() : output_cluster_key=kwargs['output_cluster_key']
+ else : cluster_key=cluster_key+".deepCopy"
if 'modulator' in kwargs.keys() : mod_tool=kwargs['modulator']
else : mod_tool=GetNullModulator()
@@ -152,18 +154,26 @@ def ApplySubtractionToClusters(**kwargs) :
if 'update_only' in kwargs.keys() : update_only = kwargs['update_only']
else : update_only = False
+ if 'apply_origin_correction' in kwargs.keys() : apply_origin_correction=kwargs['apply_origin_correction']
+ else : apply_origin_correction=HIJetFlags.ApplyOriginCorrection()
+
+ do_cluster_moments=False
+ if 'CalculateMoments' in kwargs.keys() : do_cluster_moments=kwargs['CalculateMoments']
+
HIClusterSubtraction=CompFactory.HIClusterSubtraction
toolName='HIClusterSubtraction'
if 'name' in kwargs.keys() : toolName = kwargs['name']
+
theAlg=HIClusterSubtraction(toolName)
theAlg.ClusterKey=cluster_key
+ theAlg.OutClusterKey=output_cluster_key
theAlg.EventShapeKey=event_shape_key
theAlg.Subtractor=GetSubtractorTool(**kwargs)
theAlg.Modulator=mod_tool
theAlg.UpdateOnly=update_only
+ theAlg.SetMoments=do_cluster_moments
+ theAlg.ApplyOriginCorrection=apply_origin_correction
- do_cluster_moments=False
- if 'CalculateMoments' in kwargs.keys() : do_cluster_moments=kwargs['CalculateMoments']
if do_cluster_moments :
CaloClusterMomentsMaker=CompFactory.CaloClusterMomentsMaker
from CaloTools.CaloNoiseToolDefault import CaloNoiseToolDefault
@@ -207,6 +217,25 @@ def ApplySubtractionToClusters(**kwargs) :
jtm.jetrecs += [theAlg]
jtm.HIJetRecs+=[theAlg]
+def GetConstituentsModifierTool(**kwargs) :
+ #For the cluster key, same exact logic as used for ApplySubtractionToClusters
+ if 'cluster_key' in kwargs.keys() : cluster_key=kwargs['cluster_key']
+ else : cluster_key=HIJetFlags.HIClusterKey()
+
+ if 'apply_origin_correction' in kwargs.keys() : apply_origin_correction=kwargs['apply_origin_correction']
+ else : apply_origin_correction=HIJetFlags.ApplyOriginCorrection()
+
+ HIJetConstituentModifierTool=CompFactory.HIJetConstituentModifierTool
+ toolName='HIJetConstituentModifierTool'
+ if 'name' in kwargs.keys() : toolName = kwargs['name']
+
+ cmod=HIJetConstituentModifierTool(toolName)
+ cmod.ClusterKey=cluster_key
+ cmod.Subtractor=GetSubtractorTool(**kwargs)
+ cmod.ApplyOriginCorrection=apply_origin_correction
+
+ jtm.add(cmod)
+ return cmod
def AddIteration(seed_container,shape_name, **kwargs) :
@@ -383,6 +412,7 @@ def GetHIModifierList(coll_name='AntiKt4HIJets',prepend_tools=[],append_tools=[]
if coll_name not in jtm.HICalibMap.keys() :
print ('Calibration for R=%d not available using default R=0.4 calibration')
coll_name='AntiKt4HIJets'
+
mod_list=prepend_tools
mod_list+=[jtm.HICalibMap[coll_name]]
mod_list+=jtm.modifiersMap['HI']
diff --git a/Reconstruction/HeavyIonRec/HIJetRec/share/HIJetRec_jobOptions.py b/Reconstruction/HeavyIonRec/HIJetRec/share/HIJetRec_jobOptions.py
index 6ffa3c145f093350a0b81076446493aaf5645503..1121181189174f2b11c8333d506a3b8ff88844a4 100644
--- a/Reconstruction/HeavyIonRec/HIJetRec/share/HIJetRec_jobOptions.py
+++ b/Reconstruction/HeavyIonRec/HIJetRec/share/HIJetRec_jobOptions.py
@@ -39,8 +39,8 @@ if not jetFlags.useTruth():
jetFlags.useTruth.set_Value_and_Lock(is_mc_or_overlay)
-#Tower level subtraction
-HIJetFlags.DoCellBasedSubtraction.set_Value_and_Lock(False)
+#Tower level subtraction - made it false by default to avoid confusion
+#HIJetFlags.DoCellBasedSubtraction.set_Value_and_Lock(False)
jetFlags.useTracks.set_Value_and_Lock(True)
#HIP mode
@@ -161,11 +161,6 @@ subtr1=MakeSubtractionTool(iter0.OutputEventShapeKey,moment_name="NoIteration",m
#main subtractor
subtr2=MakeSubtractionTool(HIJetFlags.IteratedEventShapeKey(),modulator=modulator1)
-#put subtraction tool at the FRONT of the jet modifiers list
-hi_tools=[subtr1,subtr2]
-hi_tools+=GetFlowMomentTools(iter1.OutputEventShapeKey,iter1.ModulationEventShapeKey)
-
-
#==========#==========#==========#==========#==========#==========
#special addition for egamma
#Downstream egamma jo will call SubtractedCellGetter, it assumes that the container pointed to by
@@ -177,10 +172,18 @@ if not HIJetFlags.DoCellBasedSubtraction():
cell_level_shape_key=iter_egamma.OutputEventShapeKey
#HIJetFlags.IteratedEventShapeKey=iter_egamma.OutputEventShapeKey
-#Subtraction for egamma and to get layers
-ApplySubtractionToClusters(name="HIClusterSubtraction_egamma", event_shape_key=cell_level_shape_key, cluster_key=ClusterKey, modulator=modulator1, CalculateMoments=True, useClusters=False)
+cluster_key_eGamma_deep=ClusterKey+"_eGamma_deep"
+cluster_key_final_deep=cluster_key_eGamma_deep+"_Cluster_deep"
+#Subtraction for egamma and to get layers - here no origin correction yet (done in the next stage)
+ApplySubtractionToClusters(name="HIClusterSubtraction_egamma", event_shape_key=cell_level_shape_key, cluster_key=ClusterKey, output_cluster_key=cluster_key_eGamma_deep, modulator=modulator1, CalculateMoments=True, useClusters=False, apply_origin_correction=False)
#Cluster subtraction for jets
-ApplySubtractionToClusters(event_shape_key=HIJetFlags.IteratedEventShapeKey(), update_only=True, cluster_key=ClusterKey, modulator=modulator1, CalculateMoments=False, useClusters=True)
+ApplySubtractionToClusters(event_shape_key=HIJetFlags.IteratedEventShapeKey(), cluster_key=cluster_key_eGamma_deep, output_cluster_key=cluster_key_final_deep, modulator=modulator1, CalculateMoments=False, useClusters=True, apply_origin_correction=HIJetFlags.ApplyOriginCorrection())
+
+#put subtraction tool at the FRONT of the jet modifiers list
+hi_tools=[subtr1,subtr2]
+hi_tools+=GetFlowMomentTools(iter1.OutputEventShapeKey,iter1.ModulationEventShapeKey)
+hi_tools+=[GetConstituentsModifierTool(name="HIJetConstituentModifierTool", cluster_key=cluster_key_final_deep, apply_origin_correction=HIJetFlags.ApplyOriginCorrection())]
+
###
#subtracted algorithms
#make main jets from unsubtr collections w/ same R, add modifiers for subtraction
diff --git a/Reconstruction/HeavyIonRec/HIJetRec/src/HIClusterSubtraction.cxx b/Reconstruction/HeavyIonRec/HIJetRec/src/HIClusterSubtraction.cxx
index 651947798b5246ee6717ef70ede8f7564d2bba96..ec43691107e963510431a78c3761b3f410ccc212 100644
--- a/Reconstruction/HeavyIonRec/HIJetRec/src/HIClusterSubtraction.cxx
+++ b/Reconstruction/HeavyIonRec/HIJetRec/src/HIClusterSubtraction.cxx
@@ -6,7 +6,12 @@
#include "xAODHIEvent/HIEventShapeContainer.h"
#include "xAODCaloEvent/CaloClusterContainer.h"
#include "HIEventUtils/HIEventShapeMap.h"
+#include "CaloUtils/CaloClusterStoreHelper.h"
#include "HIJetRec/HIJetRecDefs.h"
+#include "xAODBase/IParticleHelpers.h"
+#include "xAODTracking/Vertex.h"
+#include "xAODTracking/VertexContainer.h"
+#include "xAODCaloEvent/CaloClusterAuxContainer.h"
#include "StoreGate/ReadHandle.h"
#include "StoreGate/WriteHandle.h"
@@ -23,7 +28,10 @@ StatusCode HIClusterSubtraction::initialize()
{
//Keys initialization
ATH_CHECK( m_eventShapeKey.initialize() );
- ATH_CHECK( m_clusterKey.initialize() );
+ ATH_CHECK( m_inClusterKey.initialize() );
+ ATH_CHECK( m_outClusterKey.initialize() );
+ //Vertex container needs to be initialized only if origin correction will take place
+ ATH_CHECK( m_vertexContainer.initialize( m_originCorrection ) );
for (auto tool : m_clusterCorrectionTools)
{
@@ -34,74 +42,164 @@ StatusCode HIClusterSubtraction::initialize()
return StatusCode::SUCCESS;
}
+bool HIClusterSubtraction::doOriginCorrection( xAOD::CaloCluster* cl, const xAOD::Vertex* origin, xAOD::IParticle::FourMom_t& p4_cl ) const{
+ //made boolean to return what was "missingMoment" in HIJetConstituentSubtractionTool
+ bool missingMoment = false;
+ float mag = 0;
+ if(cl->isAvailable<float>("HIMag")) mag=cl->auxdataConst<float>("HIMag");
+ else
+ {
+ double cm_mag=0;
+ if(cl->retrieveMoment (xAOD::CaloCluster::CENTER_MAG, cm_mag)) mag=cm_mag;
+ }
+ if(mag!=0.)
+ {
+ float eta0=cl->eta0();
+ float phi0=cl->phi0();
+ float radius=mag/std::cosh(eta0);
+ xAOD::IParticle::FourMom_t p4_pos;
+ p4_pos.SetX(radius*std::cos(phi0)-origin->x());
+ p4_pos.SetY(radius*std::sin(phi0)-origin->y());
+ p4_pos.SetZ(radius*std::sinh(eta0)-origin->z());
+
+ double deta=p4_pos.Eta()-eta0;
+ double dphi=p4_pos.Phi()-phi0;
+ //adjust in case eta/phi are flipped in case of neg E clusters
+ //this method is agnostic wrt convention
+ if(p4_cl.Eta()*eta0 <0.) deta*=-1;
+
+ double eta_prime=p4_cl.Eta()+deta;
+ double phi_prime=p4_cl.Phi()+dphi;
+ double e_subtr=p4_cl.E();
+ p4_cl.SetPtEtaPhiE(e_subtr/std::cosh(eta_prime),eta_prime,phi_prime,e_subtr);
+ }
+ else missingMoment=true;
+
+ return missingMoment;
+}
+
int HIClusterSubtraction::execute() const
{
-
//const jet::cellset_t & badcells = badCellMap.cells() ;
//retrieve UE
- //From here on temporarily commented out code bc decorators needs a dedicated treatment in MT
- //const xAOD::HIEventShapeContainer* shape = 0;
+ const xAOD::HIEventShapeContainer* shape = 0;
SG::ReadHandle<xAOD::HIEventShapeContainer> readHandleEvtShape ( m_eventShapeKey );
- //shape = readHandleEvtShape.get();
- ATH_MSG_WARNING("HIClusterSubtraction not yet migrated to MT. Currently disabled! ");
- //const HIEventShapeIndex* es_index = HIEventShapeMap::getIndex( m_eventShapeKey.key() );
-
+ shape = readHandleEvtShape.cptr();
+ const HIEventShapeIndex* es_index = HIEventShapeMap::getIndex( m_eventShapeKey.key() );
const xAOD::HIEventShape* eshape = nullptr;
- //Hibryd merging between commit c2aeaed0 to master and 5af8a733 to 21.0ss
CHECK(m_modulatorTool->getShape(eshape), 1);
- //This part regards decoration!
- //Needs a more deep implementation to work in MT - to be discussed w/ Bill Balunas
- //Will be fronted in the next step of the migration
- /*
- SG::ReadHandle<xAOD::CaloClusterContainer> read_handle_clusters ( m_clusterKey );
- if (!read_handle_clusters.isValid())
- {
- ATH_MSG_ERROR("Could not retrieve input CaloClusterContainer " << m_clusterKey.key() );
- return 1;
- }
- const xAOD::CaloClusterContainer* ccl=0;
- if(m_updateMode)
+
+ //New implementation: make a deep copy of original HIClusters and apply subtraction to clusters in the new container
+ SG::ReadHandle<xAOD::CaloClusterContainer> readHandleClusters ( m_inClusterKey );
+ // Now a handle to write the deep Copy
+ SG::WriteHandle<xAOD::CaloClusterContainer> writeHandleDeepCopyClusters ( m_outClusterKey );
+ // Preparing keys and container to perfrom the origin correction
+ const xAOD::Vertex* primVertex=nullptr;
+ const xAOD::VertexContainer* vertices=nullptr;
+ // Boolean to flag that at least a vertex is present in the vertex container
+ bool isOriginPossible = true;
+ // Finding the primary vertex in case origin correction has to be performed
+ if(m_originCorrection)
{
- ccl = read_handle_clusters.get();
- //if(evtStore()->retrieve(ccl,m_clusterKey).isFailure())
- std::unique_ptr<std::vector<float> > subtractedE(new std::vector<float>());
- subtractedE->reserve(ccl->size());
- //Decoration TODO: check for migration
- SG::AuxElement::Decorator< float > decorator("HISubtractedE");
-
- for(xAOD::CaloClusterContainer::const_iterator itr=ccl->begin(); itr!=ccl->end(); itr++)
+ // ReadHandle to retrieve the vertex container
+ SG::ReadHandle<xAOD::VertexContainer> readHandleVertexContainer ( m_vertexContainer );
+ vertices = readHandleVertexContainer.get();
+ for ( size_t iVertex = 0; iVertex < vertices->size(); ++iVertex )
+ {
+ if(vertices->at(iVertex)->vertexType() == xAOD::VxType::PriVtx)
+ {
+ primVertex=vertices->at(iVertex);
+ break;
+ }
+ }
+ if(!primVertex && vertices->size() > 0)
{
- const xAOD::CaloCluster* cl=*itr;
- xAOD::IParticle::FourMom_t p4;
- m_subtractorTool->Subtract(p4,cl,shape,es_index,m_modulatorTool,eshape);
- subtractedE->push_back(p4.E());
- decorator(*cl)=p4.E();
+ ATH_MSG_WARNING("No primary vertices found, using first in container");
+ primVertex=vertices->at(0);
}
+ if(!primVertex && vertices->size() == 0)
+ {
+ ATH_MSG_WARNING("No primary vertices found, and vertex container empty. Abortin Origin correction for this event.");
+ isOriginPossible = false;
+ }
}
- else
+ bool missingMoment=false;
+
+ auto originalCluster = readHandleClusters.cptr();
+ // Create the new container and its auxiliary store.
+ xAOD::CaloClusterContainer* copyClusters = new xAOD::CaloClusterContainer();
+ xAOD::AuxContainerBase* copyClustersAux = new xAOD::AuxContainerBase();
+ copyClusters->setStore(copyClustersAux);
+ copyClusters->reserve (originalCluster->size());
+
+ for (const xAOD::CaloCluster* oldCluster : *originalCluster) {
+ xAOD::CaloCluster* newClu=new xAOD::CaloCluster();
+ copyClusters->push_back (newClu);
+ *newClu=*oldCluster;
+ }
+
+ for(xAOD::CaloClusterContainer::iterator itr=copyClusters->begin(); itr!=copyClusters->end(); itr++)
{
- xAOD::CaloClusterContainer* ccl=0;
- ccl = read_handle_clusters.get();
- for(xAOD::CaloClusterContainer::iterator itr=ccl->begin(); itr!=ccl->end(); itr++)
+ xAOD::CaloCluster* cl= *itr;
+ xAOD::IParticle::FourMom_t p4;
+ //Unsubtracted state record done by the subtractor tool functions.
+ if(m_setMoments)
+ {
+ //This flag is set to false for HIJetClustersSubtractorTool and true for HIJetCellSubtractorTool,
+ // but for the second we don't do origin correction. In principle the code is structured to do the same as the
+ //else for m_setMoments=true and HIJetClustersSubtractorTool, therefore we add the code for origin correction also here
+ m_subtractorTool->subtractWithMoments(cl, shape, es_index, m_modulatorTool, eshape);
+ if(isOriginPossible && m_originCorrection)
+ {
+ missingMoment = HIClusterSubtraction::doOriginCorrection( cl, primVertex, p4 );
+ HIJetRec::setClusterP4(p4,cl,HIJetRec::subtractedOriginCorrectedClusterState());
+ }
+ }
+ else
{
- xAOD::CaloCluster* cl=*itr;
- xAOD::IParticle::FourMom_t p4;
- if(m_setMoments) m_subtractorTool->SubtractWithMoments(cl, shape, es_index, m_modulatorTool, eshape);
- else
- {
- m_subtractorTool->Subtract(p4,cl,shape,es_index,m_modulatorTool,eshape);
- HIJetRec::setClusterP4(p4,cl,HIJetRec::subtractedClusterState());
- }
- }
+ m_subtractorTool->subtract(p4,cl,shape,es_index,m_modulatorTool,eshape);
+ HIJetRec::setClusterP4(p4,cl,HIJetRec::subtractedClusterState());
+
+ if(isOriginPossible)
+ {
+ ATH_MSG_DEBUG("Applying origin correction"
+ << std::setw(12) << "Before:"
+ << std::setw(10) << std::setprecision(3) << p4.Pt()*1e-3
+ << std::setw(10) << std::setprecision(3) << p4.Eta()
+ << std::setw(10) << std::setprecision(3) << p4.Phi()
+ << std::setw(10) << std::setprecision(3) << p4.E()*1e-3
+ << std::setw(10) << std::setprecision(3) << p4.M()*1e-3);
+
+ missingMoment = HIClusterSubtraction::doOriginCorrection( cl, primVertex, p4 );
+ HIJetRec::setClusterP4(p4,cl,HIJetRec::subtractedOriginCorrectedClusterState());
+
+ ATH_MSG_DEBUG("Applying origin correction"
+ << std::setw(12) << "After:"
+ << std::setw(10) << std::setprecision(3) << p4.Pt()*1e-3
+ << std::setw(10) << std::setprecision(3) << p4.Eta()
+ << std::setw(10) << std::setprecision(3) << p4.Phi()
+ << std::setw(10) << std::setprecision(3) << p4.E()*1e-3
+ << std::setw(10) << std::setprecision(3) << p4.M()*1e-3);
+ }
+ }
+ }//End of iterator over CaloClusterContainer
+
for(ToolHandleArray<CaloClusterCollectionProcessor>::const_iterator toolIt=m_clusterCorrectionTools.begin();
toolIt != m_clusterCorrectionTools.end(); toolIt++)
{
ATH_MSG_DEBUG(" Applying correction = " << (*toolIt)->name() );
- CHECK((*toolIt)->execute(Gaudi::Hive::currentContext(), ccl), 1);
- //Hibryd merging between commit c2aeaed0 to master and 5af8a733 to 21.0
- //eventually needs to be changed to following for r21:
- //CHECK((*toolIt)->execute(ccl), 1);
+ CHECK((*toolIt)->execute(Gaudi::Hive::currentContext(), copyClusters), 1);
}//End loop over correction tools
- }*/
+
+ if(missingMoment) ATH_MSG_WARNING("No origin correction applied, CENTERMAG missing");
+
+ // Make sure that memory is managed safely
+ std::unique_ptr<xAOD::CaloClusterContainer> outClusters(copyClusters);
+ std::unique_ptr<xAOD::AuxContainerBase> deepAux(copyClustersAux);
+
+ if(writeHandleDeepCopyClusters.record ( std::move(outClusters), std::move(deepAux)).isFailure() ){
+ ATH_MSG_ERROR("Unable to write DeepCopy Copy containers for subtracted clusters with key: " << m_outClusterKey.key());
+ return 1;
+ }
return 0;
}
diff --git a/Reconstruction/HeavyIonRec/HIJetRec/src/HIClusterSubtraction.h b/Reconstruction/HeavyIonRec/HIJetRec/src/HIClusterSubtraction.h
index 44e53ca88ff339e13ff3b80094c5421a288dd44a..7a665337dde5d8b198fee677b3e86381175642fe 100644
--- a/Reconstruction/HeavyIonRec/HIJetRec/src/HIClusterSubtraction.h
+++ b/Reconstruction/HeavyIonRec/HIJetRec/src/HIClusterSubtraction.h
@@ -27,6 +27,7 @@
#include <HIJetRec/IHISubtractorTool.h>
#include <HIJetRec/IHIUEModulatorTool.h>
#include "CaloRec/CaloClusterCollectionProcessor.h"
+#include "xAODTracking/VertexContainer.h"
#include "StoreGate/ReadHandleKey.h"
#include "StoreGate/WriteHandleKey.h"
@@ -45,12 +46,18 @@ public:
virtual int execute() const;
+ bool doOriginCorrection( xAOD::CaloCluster* cl, const xAOD::Vertex* origin, xAOD::IParticle::FourMom_t& p4_cl ) const;
+
private:
/// \brief Name of input cluster container
- SG::ReadHandleKey< xAOD::CaloClusterContainer > m_clusterKey { this, "ClusterKey", "ClusterKey", "Name of the input Cluster Container"};
+ SG::ReadHandleKey< xAOD::CaloClusterContainer > m_inClusterKey { this, "ClusterKey", "ClusterKey", "Name of the input Cluster Container"};
+ /// |brief New writeHandleKey to store the shallow copy used for new CaloClusterTreatment
+ SG::WriteHandleKey< xAOD::CaloClusterContainer > m_outClusterKey { this, "OutClusterKey", "OutClusterKey", "Name of the output Cluster Container (deep Copy)"};
/// \brief Name of HIEventShapeContainer defining background
SG::ReadHandleKey< xAOD::HIEventShapeContainer > m_eventShapeKey { this, "EventShapeKey", "EventShapeKey", "Name of HIEventShapeContainer defining background"};
+ /// |brief Name of Vertex Container for origin correction
+ SG::ReadHandleKey< xAOD::VertexContainer > m_vertexContainer { this, "VertexContainer", "PrimaryVertices", "Vertex container for primary vertices"};
// Tool handles
ToolHandle<IHISubtractorTool> m_subtractorTool { this, "Subtractor", "HIJetSubtractorToolBase", "Handle to IHISubtractorTool which does calculates subtracted kinematics" };
@@ -60,6 +67,8 @@ private:
// Booleans
Gaudi::Property< bool > m_setMoments { this, "SetMoments", true, "Set Moments boolean switch"};
Gaudi::Property< bool > m_updateMode { this, "UpdateOnly", false, "Update Mode boolean switch"};
+ //Moving the origin correction directly here
+ Gaudi::Property< bool > m_originCorrection { this, "ApplyOriginCorrection", false, "Apply Origin Correction boolean switch"};
};
#endif
diff --git a/Reconstruction/HeavyIonRec/HIJetRec/src/HIJetCellSubtractorTool.cxx b/Reconstruction/HeavyIonRec/HIJetRec/src/HIJetCellSubtractorTool.cxx
index 920fc265b04440aaf1fe148c297fd7c1c93cac03..5ca285b3b034790ae97445e6035b3d8445879708 100644
--- a/Reconstruction/HeavyIonRec/HIJetRec/src/HIJetCellSubtractorTool.cxx
+++ b/Reconstruction/HeavyIonRec/HIJetRec/src/HIJetCellSubtractorTool.cxx
@@ -30,9 +30,9 @@ void HIJetCellSubtractorTool::subtract(xAOD::IParticle::FourMom_t& subtr_mom, co
const float phi0=cl->phi0();
float mod=modulator->getModulation(phi0, eshape);
-
//unsigned int eta_phi_index=HICaloCellHelper::FindEtaPhiBin(cl->eta0(),cl->phi0());
xAOD::CaloCluster::const_cell_iterator cellIterEnd = cl->cell_end();
+
for(xAOD::CaloCluster::const_cell_iterator cellIter=cl->cell_begin(); cellIter != cellIterEnd; cellIter++ )
{
CxxUtils::prefetchNext(cellIter, cellIterEnd);
@@ -55,6 +55,7 @@ void HIJetCellSubtractorTool::subtract(xAOD::IParticle::FourMom_t& subtr_mom, co
phi_cl+=cell_E_w*phi;
}
+
if(E_cl!=0.)
{
eta_cl/=E_cl;
@@ -134,6 +135,12 @@ void HIJetCellSubtractorTool::subtractWithMoments(xAOD::CaloCluster* cl, const x
std::vector<float> E_sample(CaloSampling::Unknown,0);
uint32_t samplingPattern=0;
+
+ auto cellLink = cl->getCellLinks();
+ if( cellLink == NULL){
+ ATH_MSG_ERROR("HIJetCellSubtraction: cellLink null - returning");
+ return;
+ }
//unsigned int eta_phi_index=HICaloCellHelper::FindEtaPhiBin(cl->eta0(),cl->phi0());
xAOD::CaloCluster::cell_iterator cellIterEnd = cl->cell_end();
for(xAOD::CaloCluster::cell_iterator cellIter=cl->cell_begin(); cellIter != cellIterEnd; cellIter++ )
@@ -166,8 +173,6 @@ void HIJetCellSubtractorTool::subtractWithMoments(xAOD::CaloCluster* cl, const x
float cell_z=(*cellIter)->z();
etot2+=abs_weight;
er2+=std::sqrt(cell_x*cell_x+cell_y*cell_y+cell_z*cell_z)*abs_weight;
-
-
}
if(E_cl!=0.)
{
diff --git a/Reconstruction/HeavyIonRec/HIJetRec/src/HISubtractedCellMakerTool.cxx b/Reconstruction/HeavyIonRec/HIJetRec/src/HISubtractedCellMakerTool.cxx
index 04b80f0cd74240e888feea240d998353f7be8699..d80316b06030c42a113ebc4173bda48b475a8569 100644
--- a/Reconstruction/HeavyIonRec/HIJetRec/src/HISubtractedCellMakerTool.cxx
+++ b/Reconstruction/HeavyIonRec/HIJetRec/src/HISubtractedCellMakerTool.cxx
@@ -56,13 +56,7 @@ StatusCode HISubtractedCellMakerTool::process (CaloCellContainer* theCells,
return StatusCode::SUCCESS;
}
- // FIXME: m_modulatorTool->retrieveShape() is non-const.
- // It should be made const in order to be able to safely call it from here.
- // However, this method already needs updating to work in MT (and is checked
- // above), so just use a const_cast for now to allow this to compile
- // when ToolHandle restrictions are enabled.
-
- IHIUEModulatorTool* modtool_nc = const_cast<IHIUEModulatorTool*> (m_modulatorTool.get());
+ auto modtool_nc = m_modulatorTool.get();
CHECK(modtool_nc->retrieveShape());
for(auto pCell : *theCells)
@@ -76,8 +70,8 @@ StatusCode HISubtractedCellMakerTool::process (CaloCellContainer* theCells,
{
if( std::abs(eta) - HICaloRange::getRange().getRangeMax(sample) )
{
- double fp_round=(eta > 0.) ? -5e-3 : 5e-3;
- bin=index->getIndex(eta+fp_round,sample);
+ double fp_round=(eta > 0.) ? -5e-3 : 5e-3;
+ bin=index->getIndex(eta+fp_round,sample);
}
}
diff --git a/Reconstruction/HeavyIonRec/HIJetRec/src/components/HIJetRec_entries.cxx b/Reconstruction/HeavyIonRec/HIJetRec/src/components/HIJetRec_entries.cxx
index ac457b88e17794efd6b3fc80df057d6b3ac1e5a0..e34a0313363b73fcf0f14d017830f3b1ec06fc5b 100644
--- a/Reconstruction/HeavyIonRec/HIJetRec/src/components/HIJetRec_entries.cxx
+++ b/Reconstruction/HeavyIonRec/HIJetRec/src/components/HIJetRec_entries.cxx
@@ -8,6 +8,7 @@
#include "HIJetRec/HIEventShapeJetIteration.h"
#include "HIJetRec/HIJetClusterSubtractorTool.h"
#include "HIJetRec/HIJetConstituentSubtractionTool.h"
+#include "HIJetRec/HIJetConstituentModifierTool.h"
#include "HIJetRec/HIJetDRAssociationTool.h"
#include "HIJetRec/HIJetMaxOverMeanTool.h"
#include "HIJetRec/HIJetDiscriminatorTool.h"
@@ -21,7 +22,7 @@ DECLARE_COMPONENT( HIClusterSubtraction )
DECLARE_COMPONENT( HISubtractedCellMakerTool )
#endif
-
+DECLARE_COMPONENT( HIJetConstituentModifierTool )
DECLARE_COMPONENT( HIEventShapeJetIteration )
DECLARE_COMPONENT( HIJetConstituentSubtractionTool )
DECLARE_COMPONENT( HIJetClusterSubtractorTool )