ATLAS Offline Software
HIJetClusterSubtractorTool.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
13 #include <TH3F.h>
14 #include <TFile.h>
15 #include <iomanip>
16 
18  m_init(false),
19  m_h3W(nullptr),
20  m_h3Eta(nullptr),
21  m_h3Phi(nullptr)
22 {
23  setUseCells(false);
24 }
25 
26 void HIJetClusterSubtractorTool::subtract(xAOD::IParticle::FourMom_t& subtr_mom, const xAOD::IParticle* cl_in, const xAOD::HIEventShapeContainer* shape, const HIEventShapeIndex* index, const ToolHandle<IHIUEModulatorTool>& modulator, const xAOD::HIEventShape* eshape) const
27 {
28 
29  const xAOD::CaloCluster* cl=static_cast<const xAOD::CaloCluster*>(cl_in);
30  float E_unsubtr=cl->e(HIJetRec::unsubtractedClusterState());
31  float eta_unsubtr=cl->eta(HIJetRec::unsubtractedClusterState());
32  float phi_unsubtr=cl->phi(HIJetRec::unsubtractedClusterState());
33 
34  float E_subtr=0;
35  float E_eta_subtr=0;
36  float E_phi_subtr=0;
37  const float eta0=cl->eta0();
38  const float phi0=cl->phi0();
39  float mod=1;
40  if(modulator) mod=modulator->getModulation(phi0, eshape);
41 
42  float energy=E_unsubtr;
43  float eta=eta0;
44  float phi=phi0;
45 
46  static const SG::ConstAccessor<float> HIEtaPhiWeightAcc("HIEtaPhiWeight");
47  if(HIEtaPhiWeightAcc.isAvailable(*cl))
48  {
49  float DF_weight=HIEtaPhiWeightAcc(*cl);
50  if(DF_weight!=0.) mod/=DF_weight;
51  }
52 
53  if(shape==nullptr)
54  {
55  ATH_MSG_ERROR("No HIEventShape supplied, cannot do subtraction");
56  return;
57  }
58 
59  if(m_useSamplings)
60  {
61  if(index==nullptr)
62  {
63  ATH_MSG_ERROR("No HIEventShapeIndex supplied, cannot do subtraction");
64  return;
65  }
66 
67  for(unsigned int sample=0; sample<24; sample++)
68  {
70  if(!cl->hasSampling(s)) continue;
71 
72  float nCells=index->getShape(eta0,sample,shape)->nCells();
73  float rho=0;
74  if(nCells!=0.) rho=index->getShape(eta0,sample,shape)->rho()/nCells;
75  rho*=mod*std::cosh(eta0);
76 
77  E_subtr+=rho*getWeight(eta0,phi0,sample);
78  E_eta_subtr+=rho*getWeightEta(eta0,phi0,sample);
79  E_phi_subtr+=rho*getWeightPhi(eta0,phi0,sample);
80  }
81 
82  energy=E_unsubtr-E_subtr;
83  if(energy!=0.)
84  {
85  eta=eta_unsubtr+(eta_unsubtr*E_subtr-E_eta_subtr)/energy;
86  phi=phi_unsubtr+(phi_unsubtr*E_subtr-E_phi_subtr)/energy;
87  }
89  }
90  else
91  {
92  unsigned int es_bin=HI::TowerBins::findBinEta(eta0);
93 
94  float area=shape->at(es_bin)->area();
95  float rho=(area==0.) ? 0. : shape->at(es_bin)->et()/area;
96 
97  rho*=mod*std::cosh(eta0);
99  energy=E_unsubtr-E_subtr;
100  eta=eta_unsubtr;
101  phi=phi_unsubtr;
102 
103  ATH_MSG_VERBOSE("Subtracting"
104  << std::setw(8) << eta0
105  << std::setw(8) << phi0
106  << std::setw(10) << std::setprecision(3) << modulator->getModulation(phi0, eshape)
107  << std::setw(10) << std::setprecision(3) << mod
108  << std::setw(10) << std::setprecision(3) << modulator->getModulation(phi0, eshape)/mod
109  << std::setw(10) << std::setprecision(3) << area
110  << std::setw(10) << std::setprecision(3) << rho*area
111  << std::setw(10) << std::setprecision(3) << E_unsubtr*1e-3
112  << std::setw(10) << std::setprecision(3) << E_subtr*1e-3
113  << std::setw(10) << std::setprecision(3) << E_unsubtr*1e-3/std::cosh(eta0)
114  << std::setw(10) << std::setprecision(3) << E_subtr*1e-3/std::cosh(eta0));
115 
116 
117 
118 
119  }
120  float ET=energy/std::cosh(eta);
121  subtr_mom.SetPxPyPzE(ET*std::cos(phi),ET*std::sin(phi),ET*std::sinh(eta),energy);
122 }
123 
125 {
126 
127  float eta0=cl->eta0();
128  float phi0=cl->phi0();
129 
130  if(m_useSamplings)
131  {
132  for(unsigned int sample=0; sample<24; sample++)
133  {
135  if(!cl->hasSampling(s)) continue;
136  float esamp=cl->eSample(s);
137  float ET=esamp/std::cosh(eta0);
138  xAOD::HIEventShape* slice=shape->at(index->getIndex(eta0,sample));
139  float area_cluster=getWeight(eta0,phi0,sample);
140  //update members
141  updateSlice(slice,ET,phi0,area_cluster);
142  }
143  }
144  else
145  {
146  unsigned int es_bin=HI::TowerBins::findBinEta(eta0);
147  xAOD::HIEventShape* slice=shape->at(es_bin);
148  constexpr float area_cluster=HI::TowerBins::getBinArea();
149  float ET=cl->e(HIJetRec::unsubtractedClusterState())/std::cosh(eta0);
150  static const SG::ConstAccessor<float> HIEtaPhiWeightAcc("HIEtaPhiWeight");
151  if(HIEtaPhiWeightAcc.isAvailable(*cl))
152  {
153  float HI_weight=HIEtaPhiWeightAcc(*cl);
154  ET*=HI_weight;
155  }
156  //update members
157  updateSlice(slice,ET,phi0,area_cluster);
158  }
159 }
160 
161 
162 void HIJetClusterSubtractorTool::subtractWithMoments(xAOD::CaloCluster* cl, const xAOD::HIEventShapeContainer* shape, const HIEventShapeIndex* index, const ToolHandle<IHIUEModulatorTool>& modulator, const xAOD::HIEventShape* eshape) const
163 {
164  xAOD::IParticle::FourMom_t subtr_mom;
165  subtract(subtr_mom, cl, shape, index, modulator, eshape);
167 }
168 
170 {
171  if(m_useSamplings)
172  {
173  std::string local_path=static_cast<std::string>(m_configDir)+m_inputFile;
174  std::string full_path=PathResolverFindCalibFile(local_path);
175  ATH_MSG_INFO("Reading input file "<< m_inputFile << " from " << full_path);
176  TFile* f=TFile::Open(full_path.c_str());
177  if(f==nullptr)
178  {
179  ATH_MSG_FATAL("Cannot read config file " << full_path );
180  return StatusCode::FAILURE;
181  }
182 
183  m_h3W=(TH3F*)f->GetObjectChecked("h3_w","TH3F");
184  if(m_h3W==nullptr)
185  {
186  ATH_MSG_FATAL("Cannot find TH3F m_h3W in config file " << full_path );
187  return StatusCode::FAILURE;
188  }
189 
190  m_h3Eta=(TH3F*)f->GetObjectChecked("h3_eta","TH3F");
191  if(m_h3Eta==nullptr)
192  {
193  ATH_MSG_FATAL("Cannot find TH3F m_h3Eta in config file " << full_path );
194  return StatusCode::FAILURE;
195  }
196 
197  m_h3Phi=(TH3F*)f->GetObjectChecked("h3_phi","TH3F");
198  if(m_h3Phi==nullptr)
199  {
200  ATH_MSG_FATAL("Cannot find TH3F m_h3Phi in config file " << full_path );
201  return StatusCode::FAILURE;
202  }
203  m_h3W->SetDirectory(0);
204  m_h3Eta->SetDirectory(0);
205  m_h3Phi->SetDirectory(0);
206 
207  f->Close();
208  }
209  m_init=true;
210  return StatusCode::SUCCESS;
211 
212 }
213 
214 
216 {
217  return initializeTool();
218 
219 }
220 
221 
222 float HIJetClusterSubtractorTool::getWeight(float eta, float phi, int sample) const
223 {
224  return m_h3W->GetBinContent(m_h3W->FindFixBin(eta,phi,sample));
225 }
227 {
228  return m_h3Eta->GetBinContent(m_h3Eta->FindFixBin(eta,phi,sample));
229 }
231 {
232  return m_h3Phi->GetBinContent(m_h3Phi->FindFixBin(eta,phi,sample));
233 }
234 
235 void HIJetClusterSubtractorTool::updateSlice(xAOD::HIEventShape* slice, float ET, float phi0, float area_cluster) const
236 {
237  if(area_cluster==0.)
238  {
239  area_cluster=1.;
240 
241  if(ET==0.)
242  {
243  ATH_MSG_INFO("Provided cluster area was 0, check if the input file " << m_inputFile << " from " << PathResolverFindCalibFile(static_cast<std::string>(m_configDir)+m_inputFile) << " was generated with the currently used geometry. See ATLHI-472");
244  }
245  else
246  {
247  ATH_MSG_ERROR("Provided cluster area was 0, the HIEventShape will be incorrect! Check if the input file " << m_inputFile << " from " << PathResolverFindCalibFile(static_cast<std::string>(m_configDir)+m_inputFile) << " was generated with the currently used geometry. See ATLHI-472");
248  }
249  }
250 
251  float area_slice=slice->area();
252  float area_sf=1.;
253 
254  if(area_slice!=0.) area_sf-=area_cluster/area_slice;
255  slice->setNCells( std::floor(area_sf*slice->nCells()) );
256  slice->setEt(slice->et()-ET);
257  slice->setArea(area_slice-area_cluster);
258 
259  slice->setRho(slice->rho() - ET/area_cluster);
260 
261 
262  for(unsigned int ih=0; ih<slice->etCos().size(); ih++)
263  {
264  float ih_f=ih+1;
265  float tmp_cos = slice->etCos().at(ih);
266  slice->etCos()[ih] = tmp_cos + cos(ih_f*phi0)*ET;
267 
268  float tmp_sin = slice->etSin().at(ih);
269  slice->etSin()[ih] = tmp_sin + sin(ih_f*phi0)*ET;
270  }
271 }
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
HIJetClusterSubtractorTool::m_h3W
TH3F * m_h3W
Definition: HIJetClusterSubtractorTool.h:64
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
HIEventShapeIndex.h
HIJetClusterSubtractorTool::getWeightEta
float getWeightEta(float eta, float phi, int sample) const
Definition: HIJetClusterSubtractorTool.cxx:226
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:67
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:83
index
Definition: index.py:1
HI::TowerBins::findBinEta
unsigned int findBinEta(float eta)
Definition: HIEventDefs.h:46
InDetAccessor::phi0
@ phi0
Definition: InDetAccessor.h:33
HIJetSubtractorToolBase::setUseCells
void setUseCells(bool v)
Definition: HIJetSubtractorToolBase.h:64
HIJetClusterSubtractorTool::HIJetClusterSubtractorTool
HIJetClusterSubtractorTool(const std::string &myname)
Implements method defined in base First argument is reference to four vector that is updated to refle...
Definition: HIJetClusterSubtractorTool.cxx:17
SG::ConstAccessor< float >
HIJetClusterSubtractorTool::updateSlice
void updateSlice(xAOD::HIEventShape *slice, float ET, float phi0, float area_cluster) const
Definition: HIJetClusterSubtractorTool.cxx:235
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
xAOD::IParticle
Class providing the definition of the 4-vector interface.
Definition: Event/xAOD/xAODBase/xAODBase/IParticle.h:41
HI::TowerBins::getBinArea
constexpr float getBinArea()
Definition: HIEventDefs.h:34
HIJetClusterSubtractorTool::m_h3Eta
TH3F * m_h3Eta
Definition: HIJetClusterSubtractorTool.h:65
xAOD::IParticle::FourMom_t
TLorentzVector FourMom_t
Definition of the 4-momentum type.
Definition: Event/xAOD/xAODBase/xAODBase/IParticle.h:69
HIJetClusterSubtractorTool::m_init
bool m_init
Definition: HIJetClusterSubtractorTool.h:53
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
xAOD::HIEventShape_v2
Interface class for the HI reconstruction EDM.
Definition: HIEventShape_v2.h:32
xAOD::CaloCluster_v1
Description of a calorimeter cluster.
Definition: CaloCluster_v1.h:62
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
HIJetRec::subtractedClusterState
constexpr xAOD::CaloCluster::State subtractedClusterState()
Definition: HIJetRecDefs.h:26
maskDeadModules.mod
mod
Definition: maskDeadModules.py:36
ParticleGun_FastCalo_ChargeFlip_Config.energy
energy
Definition: ParticleGun_FastCalo_ChargeFlip_Config.py:78
FullCPAlgorithmsTest_eljob.sample
sample
Definition: FullCPAlgorithmsTest_eljob.py:116
perfmonmt-refit.slice
slice
Definition: perfmonmt-refit.py:52
xAOD::nCells
setRawEt setRawPhi nCells
Definition: TrigCaloCluster_v1.cxx:33
CaloSampling::CaloSample
CaloSample
Definition: Calorimeter/CaloGeoHelpers/CaloGeoHelpers/CaloSampling.h:22
CaloCluster.h
HIJetClusterSubtractorTool::updateUsingCluster
virtual void updateUsingCluster(xAOD::HIEventShapeContainer *shape, const HIEventShapeIndex *index, const xAOD::CaloCluster *cl) const override
Method to update the shape based on a given cluster two sets of indices are passed by reference and u...
Definition: HIJetClusterSubtractorTool.cxx:124
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
HIJetClusterSubtractorTool::m_useSamplings
Gaudi::Property< bool > m_useSamplings
Definition: HIJetClusterSubtractorTool.h:68
HIJetClusterSubtractorTool::subtractWithMoments
virtual void subtractWithMoments(xAOD::CaloCluster *, const xAOD::HIEventShapeContainer *, const HIEventShapeIndex *index, const ToolHandle< IHIUEModulatorTool > &, const xAOD::HIEventShape *eshape) const override
Definition: HIJetClusterSubtractorTool.cxx:162
HIJetSubtractorToolBase
Abstract base for JetModifiers applying constituent based subtraction.
Definition: HIJetSubtractorToolBase.h:24
hist_file_dump.f
f
Definition: hist_file_dump.py:140
HIJetClusterSubtractorTool::getWeightPhi
float getWeightPhi(float eta, float phi, int sample) const
Definition: HIJetClusterSubtractorTool.cxx:230
DataVector
Derived DataVector<T>.
Definition: DataVector.h:794
PathResolver.h
HIEventShapeIndex
Definition: HIEventShapeIndex.h:16
HIJetClusterSubtractorTool::m_inputFile
Gaudi::Property< std::string > m_inputFile
Definition: HIJetClusterSubtractorTool.h:60
EventInfo.h
PathResolverFindCalibFile
std::string PathResolverFindCalibFile(const std::string &logical_file_name)
Definition: PathResolver.cxx:283
HIJetClusterSubtractorTool::m_h3Phi
TH3F * m_h3Phi
Definition: HIJetClusterSubtractorTool.h:66
HIJetClusterSubtractorTool::initializeTool
virtual StatusCode initializeTool()
Definition: HIJetClusterSubtractorTool.cxx:169
IHIUEModulatorTool.h
HIJetRec::unsubtractedClusterState
constexpr xAOD::CaloCluster::State unsubtractedClusterState()
Definition: HIJetRecDefs.h:25
HIJetClusterSubtractorTool::getWeight
float getWeight(float eta, float phi, int sample) const
Definition: HIJetClusterSubtractorTool.cxx:222
HIJetClusterSubtractorTool::m_configDir
Gaudi::Property< std::string > m_configDir
Definition: HIJetClusterSubtractorTool.h:62
python.SystemOfUnits.s
float s
Definition: SystemOfUnits.py:147
HIJetSubtractorToolBase::setSubtractedEtaPhi
void setSubtractedEtaPhi(float E, float &eta, float &phi, float eta0, float phi0, float sig) const
Definition: HIJetSubtractorToolBase.cxx:23
SG::ConstAccessor::isAvailable
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
HIJetClusterSubtractorTool::initialize
virtual StatusCode initialize() override
Dummy implementation of the initialisation function.
Definition: HIJetClusterSubtractorTool.cxx:215
HIJetRecDefs.h
HIJetClusterSubtractorTool::subtract
virtual void subtract(xAOD::IParticle::FourMom_t &, const xAOD::IParticle *, const xAOD::HIEventShapeContainer *, const HIEventShapeIndex *, const ToolHandle< IHIUEModulatorTool > &, const xAOD::HIEventShape *eshape) const override
Abstract method where particle itself is not modified IParticle::FourMom_t containing kinematics afte...
Definition: HIJetClusterSubtractorTool.cxx:26
DataVector::at
const T * at(size_type n) const
Access an element, as an rvalue.
ConstAccessor.h
Helper class to provide constant type-safe access to aux data.
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
area
double area(double R)
Definition: ConvertStaveServices.cxx:42
dq_make_web_display.cl
cl
print [x.__class__ for x in toList(dqregion.getSubRegions()) ]
Definition: dq_make_web_display.py:25
HIJetRec::setClusterP4
void setClusterP4(const xAOD::CaloCluster::FourMom_t &p, xAOD::CaloCluster *cl, xAOD::CaloCluster::State s)
Definition: HIJetRecDefs.h:36
HIJetClusterSubtractorTool.h
fitman.rho
rho
Definition: fitman.py:532