ATLAS Offline Software
GetLCClassification.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 //-----------------------------------------------------------------------
6 // File and Version Information:
7 // $Id: GetLCClassification.cxx,v 1.1.1.1 2008-11-04 08:56:11 menke Exp $
8 //
9 // Description: see GetLCClassification.h
10 //
11 // Environment:
12 // Software developed for the ATLAS Detector at CERN LHC
13 //
14 // Author List:
15 // Sven Menke
16 //
17 //-----------------------------------------------------------------------
18 
19 //-----------------------
20 // This Class's Header --
21 //-----------------------
23 
24 //---------------
25 // C++ Headers --
26 //---------------
29 #include "StoreGate/ReadHandle.h"
30 
32 #include "GaudiKernel/ISvcLocator.h"
33 #include "GaudiKernel/StatusCode.h"
34 
35 #include "TFile.h"
36 #include "TH2F.h"
37 #include "TString.h"
38 #include <iterator>
39 #include <cmath>
40 
41 //###############################################################################
43  ISvcLocator* pSvcLocator)
44  : AthAlgorithm(name, pSvcLocator),
45  m_outputFile(nullptr),
46  m_clusterCollName("CaloTopoClusters"),
47  m_NormalizationType("Lin"),
48  m_NormalizationTypeNumber(0),
49  m_ClassificationType("None"),
50  m_ClassificationTypeNumber(0)
51 {
52 
53  std::vector<Gaudi::Histo1DDef> dims(6);
54  dims[0] = Gaudi::Histo1DDef("side",-1.5,1.5,1);
55  dims[1] = Gaudi::Histo1DDef("|eta|",0.,5.,25);
56  dims[2] = Gaudi::Histo1DDef("phi",-M_PI,M_PI,1);
57  dims[3] = Gaudi::Histo1DDef("log10(E_clus (MeV))",log10(200),log10(1e6),14);
58  dims[4] = Gaudi::Histo1DDef("log10(<rho_cell (MeV/mm^3)>)-log10(E_clus (MeV))",-9.0,-4.0,20);
59  dims[5] = Gaudi::Histo1DDef("log10(lambda_clus (mm))",0.0,4.0,20);
60 
61  mapinsert(dims);
62 
63  // Dimensions to use for classification
64  declareProperty("ClassificationDimensions",m_dimensionsmap);
65 
66  // Name of output file to save histograms in
67  declareProperty("OutputFileName",m_outputFileName);
68 
69  // Name of ClusterContainer to use
70  declareProperty("ClusterCollectionName",m_clusterCollName);
71 
72  // Normalization type "Const", "Lin", "Log", "NClus"
73  declareProperty("NormalizationType",m_NormalizationType);
74 
75  // Classification type "None" for single pion MC or
76  // "ParticleID_EM" for ParticleID based em-type clusters
77  // "ParticleID_HAD" for ParticleID based had-type clusters
78  declareProperty("ClassificationType", m_ClassificationType);
79 
80 }
81 
82 //###############################################################################
83 
85 { }
86 
87 //###############################################################################
88 
90 {
91  m_outputFile = std::make_unique<TFile>(m_outputFileName.c_str(),"RECREATE");
92  m_outputFile->cd();
93  m_hclus.resize(0);
94  mapparse();
95 
96  if ( m_NormalizationType == "Lin" ) {
97  ATH_MSG_INFO( "Using weighting proportional to E_calib" );
99  }
100  else if ( m_NormalizationType == "Log" ) {
101  ATH_MSG_INFO( "Using weighting proportional to log(E_calib)" );
103  }
104  else if ( m_NormalizationType == "NClus" ) {
105  ATH_MSG_INFO( "Using weighting proportional to 1/N_Clus_E_calib>0" );
107  }
108  else {
109  ATH_MSG_INFO( "Using constant weighting" );
111  }
112 
113  if ( m_ClassificationType == "None" ) {
114  ATH_MSG_INFO( "Expecting single particle input" );
116  }
117  else if ( m_ClassificationType == "ParticleID_EM" ) {
118  ATH_MSG_INFO( "Expecting ParticleID simulation as input -- use EM type clusters only" );
120  }
121  else if ( m_ClassificationType == "ParticleID_HAD" ) {
122  ATH_MSG_INFO( "Expecting ParticleID simulation as input -- use HAD type clusters only" );
124  }
125  else {
126  ATH_MSG_WARNING( " unknown classification type " << m_ClassificationType << " given! Using None instead" );
128  }
129 
130  int iside(-1);
131  int ieta(-1);
132  int iphi(-1);
133  int ilogE(-1);
134  int ilogrho(-1);
135  int iloglambda(-1);
136 
137  m_isampmap.resize(4,-1);
138  for(unsigned int idim=0;idim<m_dimensions.size();idim++) {
139  if ( m_dimensions[idim].title() == "side" ) {
140  iside = idim;
141  m_isampmap[0] = iside;
142  }
143  else if ( m_dimensions[idim].title() == "|eta|" ) {
144  ieta = idim;
145  m_isampmap[1] = ieta;
146  }
147  else if ( m_dimensions[idim].title() == "phi" ) {
148  iphi = idim;
149  m_isampmap[2] = iphi;
150  }
151  else if ( m_dimensions[idim].title() == "log10(E_clus (MeV))" ) {
152  ilogE = idim;
153  m_isampmap[3] = ilogE;
154  }
155  else if ( m_dimensions[idim].title() == "log10(<rho_cell (MeV/mm^3)>)-log10(E_clus (MeV))" )
156  ilogrho = idim;
157  else if ( m_dimensions[idim].title() == "log10(lambda_clus (mm))" )
158  iloglambda = idim;
159  }
160  if ( ilogE < 0 || ilogrho < 0 || iloglambda < 0 ) {
161  ATH_MSG_FATAL( " Mandatory dimension log10E, log10rho or log10lambda missing ..." );
162  return StatusCode::FAILURE;
163  }
164  int nside = 1;
165  if (iside>=0) nside = m_dimensions[iside].bins();
166  //
167  int neta = 1;
168  if (ieta>=0) neta = m_dimensions[ieta].bins();
169  //
170  int nphi = 1;
171  if (iphi>=0) nphi = m_dimensions[iphi].bins();
172  int nlogE = m_dimensions[ilogE].bins();
173  m_hclus.resize(nside*neta*nphi*nlogE,nullptr);
174  for ( int jside=0;jside<nside;jside++) {
175  for ( int jeta=0;jeta<neta;jeta++) {
176  for ( int jphi=0;jphi<nphi;jphi++) {
177  for(int jlogE=0;jlogE<nlogE;jlogE++) {
178  TString hname("hclus");
179  hname += "_iside_";
180  hname += jside;
181  hname += "_[";
182  hname += (iside>=0?m_dimensions[iside].lowEdge():-1);
183  hname += ",";
184  hname += (iside>=0?m_dimensions[iside].highEdge():-1);
185  hname += ",";
186  hname += nside;
187  hname += "]";
188  hname += "_ieta_";
189  hname += jeta;
190  hname += "_[";
191  hname += (ieta>=0?m_dimensions[ieta].lowEdge():-1);
192  hname += ",";
193  hname += (ieta>=0?m_dimensions[ieta].highEdge():-1);
194  hname += ",";
195  hname += neta;
196  hname += "]";
197  hname += "_iphi_";
198  hname += jphi;
199  hname += "_[";
200  hname += (iphi>=0?m_dimensions[iphi].lowEdge():-1);
201  hname += ",";
202  hname += (iphi>=0?m_dimensions[iphi].highEdge():-1);
203  hname += ",";
204  hname += nphi;
205  hname += "]";
206  hname += "_ilogE_";
207  hname += jlogE;
208  hname += "_[";
209  hname += m_dimensions[ilogE].lowEdge();
210  hname += ",";
211  hname += m_dimensions[ilogE].highEdge();
212  hname += ",";
213  hname += nlogE;
214  hname += "]";
215  int iH = jlogE*nphi*neta*nside+jphi*neta*nside+jeta*nside+jside;
216  m_hclus[iH] = new TH2F(hname,hname,
217  m_dimensions[ilogrho].bins(),
218  m_dimensions[ilogrho].lowEdge(),
219  m_dimensions[ilogrho].highEdge(),
220  m_dimensions[iloglambda].bins(),
221  m_dimensions[iloglambda].lowEdge(),
222  m_dimensions[iloglambda].highEdge());
223  m_hclus[iH]->Sumw2();
224  m_hclus[iH]->SetXTitle("log10(<#rho_{cell}> (MeV/mm^{3})) - log10(E_{clus} (MeV))");
225  m_hclus[iH]->SetYTitle("log10(#lambda_{clus} (mm))");
226  m_hclus[iH]->SetZTitle("Number of Clusters");
227  }
228  }
229  }
230  }
231 
233 
234  return StatusCode::SUCCESS;
235 }
236 
237 //###############################################################################
238 
240 {
241  ATH_MSG_INFO( "Writing out histograms" );
242  m_outputFile->cd();
243  for(unsigned int i=0;i<m_hclus.size();i++) {
244  m_hclus[i]->Write();
245  }
246  m_outputFile->Close();
247 
248  return StatusCode::SUCCESS;
249 }
250 
251 //###############################################################################
252 
254 {
256 
257  // total calib hit energy of all clusters
258  double eCalibTot(0.);
259  double nClusECalibGt0 = 0.0;
260 
261  for (const xAOD::CaloCluster * theCluster : *cc) {
262  double eC=999;
263  if (!theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_TOT,eC)) {
264  ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_TOT" );
265  return StatusCode::FAILURE;
266  }
268  double emFrac = 0;
269  if (!theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_FRAC_EM,emFrac)) {
270  ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_FRAC_EM" );
271  return StatusCode::FAILURE;
272  }
274  eC = 0;
276  eC = 0;
277  }
278  eCalibTot += eC;
279  if ( eC > 0 ) {
280  nClusECalibGt0++;
281  }
282  }
283 
284  if ( eCalibTot > 0 ) {
285  const double inv_eCalibTot = 1. / eCalibTot;
286  const double inv_nClusECalibGt0 = 1. / nClusECalibGt0;
287  for (const xAOD::CaloCluster * pClus : *cc) {
288  double eng = pClus->e();
289  if ( eng > 0 ) {
291  double emFrac = 0;
292  if (pClus->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_FRAC_EM,emFrac)) {
293  ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_FRAC_EM" );
294  return StatusCode::FAILURE;
295  }
297  continue;
299  continue;
300  }
301 
302  double eta = fabs(pClus->eta());
303 
304  int iside = 0;
305  int ieta = 0;
306  int iphi = 0;
307  int ilogE = 0;
308  int nside = 1;
309  int neta = 1;
310  int nphi = 1;
311  int nlogE = 1;
312 
313  if ( m_isampmap[0] >= 0 ) {
314  const Gaudi::Histo1DDef & hd = m_dimensions[m_isampmap[0]];
315  nside = hd.bins();
316  iside = (int)(nside*(((pClus->eta()<0?-1.0:1.0) - hd.lowEdge())
317  /(hd.highEdge()-hd.lowEdge())));
318  if ( iside < 0 || iside > nside-1 ) {
319  ATH_MSG_WARNING( " Side index out of bounds " <<
320  iside << " not in [0," << nside-1 << "]" );
321  iside = -1;
322  }
323  }
324 
325  if ( m_isampmap[1] >= 0 ) {
326  const Gaudi::Histo1DDef & hd = m_dimensions[m_isampmap[1]];
327  neta = hd.bins();
328  ieta = (int)(neta*((eta - hd.lowEdge())
329  /(hd.highEdge()-hd.lowEdge())));
330  if ( ieta < 0 || ieta > neta-1 ) {
331  ATH_MSG_WARNING( " Eta index out of bounds " <<
332  ieta << " not in [0," << neta-1 << "]" );
333  ieta = -1;
334  }
335  }
336  if ( m_isampmap[2] >= 0 ) {
337  const Gaudi::Histo1DDef & hd = m_dimensions[m_isampmap[2]];
338  nphi = hd.bins();
339  iphi = (int)(nphi*((pClus->phi() - hd.lowEdge())
340  /(hd.highEdge()-hd.lowEdge())));
341  if ( iphi < 0 || iphi > nphi-1 ) {
342  ATH_MSG_WARNING( " Phi index out of bounds " <<
343  iphi << " not in [0," << nphi-1 << "]" );
344  iphi = -1;
345  }
346  }
347 
348  const Gaudi::Histo1DDef & hd = m_dimensions[m_isampmap[3]];
349  nlogE = hd.bins();
350  ilogE = (int)(nlogE*((log10(eng) - hd.lowEdge())
351  /(hd.highEdge()-hd.lowEdge())));
352  if ( ilogE >= 0 && ilogE < nlogE ) {
353  double dens=0,lamb=0,ecal=0;
354  if (!pClus->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_TOT,ecal)) {
355  ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_TOT" );
356  return StatusCode::FAILURE;
357  }
358  if (!pClus->retrieveMoment(xAOD::CaloCluster::FIRST_ENG_DENS,dens)) {
359  ATH_MSG_ERROR( "Failed to retrieve cluster moment FIRST_ENG_DENS" );
360  return StatusCode::FAILURE;
361  }
362  if (!pClus->retrieveMoment(xAOD::CaloCluster::CENTER_LAMBDA,lamb)) {
363  ATH_MSG_ERROR( "Failed to retrieve cluster moment CENTER_LAMBDA" );
364  return StatusCode::FAILURE;
365  }
366  if ( dens > 0 &&
367  lamb > 0 &&
368  ecal > 0 )
369  {
370  int iH = ilogE*nphi*neta*nside+iphi*neta*nside+ieta*nside+iside;
371  if ( m_hclus[iH]) {
372  double norm = 0.0;
374  norm = ecal*inv_eCalibTot;
375  }
377  // cluster has to have at least 1% of the calib hit E
378  norm = log10(ecal*inv_eCalibTot)+2.0;
379  }
381  norm = inv_nClusECalibGt0;
382  }
383  else {
384  norm = 1.0;
385  }
386  if ( norm > 0 ) {
387  m_hclus[iH]->Fill(log10(dens)-log10(eng),log10(lamb),norm);
388  }
389  }
390  }
391  }
392  }
393  }
394  }
395 
396  return StatusCode::SUCCESS;
397 }
398 
399 void GetLCClassification::mapinsert(const std::vector<Gaudi::Histo1DDef> & dims) {
400  for (unsigned int i=0;i<dims.size();i++) {
401  m_dimensionsmap[dims[i].title()] = dims[i];
402  }
403 }
404 
406 
407  for (std::pair<const std::string, Gaudi::Histo1DDef>& p : m_dimensionsmap) {
408  m_dimensions.push_back(p.second);
409  ATH_MSG_DEBUG( " New Dimension: "
410  << p.second.title() << ", [" << p.second.lowEdge()
411  << ", " << p.second.highEdge()
412  << ", " << p.second.bins()
413  << "]" );
414  }
415 }
416 
417 //###############################################################################
GetLCDefs::CONST
@ CONST
Definition: GetLCDefs.h:21
GetLCClassification::~GetLCClassification
virtual ~GetLCClassification()
Definition: GetLCClassification.cxx:84
PlotCalibFromCool.norm
norm
Definition: PlotCalibFromCool.py:100
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
GetLCClassification::m_NormalizationType
std::string m_NormalizationType
string to choose different normalization types
Definition: GetLCClassification.h:123
GetLCClassification::GetLCClassification
GetLCClassification(const std::string &name, ISvcLocator *pSvcLocator)
Definition: GetLCClassification.cxx:42
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
dqt_zlumi_pandas.hname
string hname
Definition: dqt_zlumi_pandas.py:279
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:83
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
GetLCDefs.h
M_PI
#define M_PI
Definition: ActiveFraction.h:11
xAOD::CaloCluster_v1::CENTER_LAMBDA
@ CENTER_LAMBDA
Shower depth at Cluster Centroid.
Definition: CaloCluster_v1.h:136
GetLCClassification::mapinsert
void mapinsert(const std::vector< Gaudi::Histo1DDef > &dims)
Definition: GetLCClassification.cxx:399
GetLCClassification::m_hclus
std::vector< TH2F * > m_hclus
Vector of actual histograms.
Definition: GetLCClassification.h:93
GetLCClassification::m_dimensionsmap
std::map< std::string, Gaudi::Histo1DDef > m_dimensionsmap
property to set all dimensions introduced above
Definition: GetLCClassification.h:81
python.TrigEgammaMonitorHelper.TH2F
def TH2F(name, title, nxbins, bins_par2, bins_par3, bins_par4, bins_par5=None, bins_par6=None, path='', **kwargs)
Definition: TrigEgammaMonitorHelper.py:45
GetLCDefs::NONE
@ NONE
Definition: GetLCDefs.h:25
GetLCDefs::NCLUS
@ NCLUS
Definition: GetLCDefs.h:21
GetLCDefs::PARTICLEID_EM
@ PARTICLEID_EM
Definition: GetLCDefs.h:25
GetLCClassification::initialize
virtual StatusCode initialize()
Definition: GetLCClassification.cxx:89
xAOD::CaloCluster_v1
Description of a calorimeter cluster.
Definition: CaloCluster_v1.h:59
GetLCClassification::m_dimensions
std::vector< Gaudi::Histo1DDef > m_dimensions
definition of all dimensions used for classification
Definition: GetLCClassification.h:73
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
GetLCClassification::finalize
virtual StatusCode finalize()
Definition: GetLCClassification.cxx:239
lumiFormat.i
int i
Definition: lumiFormat.py:85
GetLCClassification::m_ClassificationTypeNumber
int m_ClassificationTypeNumber
Definition: GetLCClassification.h:138
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
xAOD::CaloCluster_v1::ENG_CALIB_TOT
@ ENG_CALIB_TOT
Calibration Hit energy inside the cluster.
Definition: CaloCluster_v1.h:195
GetLCClassification::mapparse
void mapparse()
Definition: GetLCClassification.cxx:405
GetLCClassification::execute
virtual StatusCode execute()
Definition: GetLCClassification.cxx:253
covarianceTool.title
title
Definition: covarianceTool.py:542
xAOD::CaloCluster_v1::ENG_CALIB_FRAC_EM
@ ENG_CALIB_FRAC_EM
Calibration Hit energy inside the cluster caused by e/gamma/pi0.
Definition: CaloCluster_v1.h:248
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
plotting.yearwise_luminosity_vs_mu.bins
bins
Definition: yearwise_luminosity_vs_mu.py:30
SG::VarHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:103
xAOD::CaloCluster_v1::FIRST_ENG_DENS
@ FIRST_ENG_DENS
First Moment in E/V.
Definition: CaloCluster_v1.h:143
AthAlgorithm
Definition: AthAlgorithm.h:47
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
errorcheck.h
Helpers for checking error return status codes and reporting errors.
GetLCClassification::m_clusterCollName
SG::ReadHandleKey< xAOD::CaloClusterContainer > m_clusterCollName
Name of the CaloClusterContainer to use.
Definition: GetLCClassification.h:111
GetLCClassification::m_ClassificationType
std::string m_ClassificationType
string to choose different classification types
Definition: GetLCClassification.h:136
GetLCClassification::m_outputFile
std::unique_ptr< TFile > m_outputFile
Output file to save histograms in.
Definition: GetLCClassification.h:106
GetLCDefs::PARTICLEID_HAD
@ PARTICLEID_HAD
Definition: GetLCDefs.h:25
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
GetLCClassification::m_outputFileName
std::string m_outputFileName
Name of the output file to save histograms in.
Definition: GetLCClassification.h:100
GetLCDefs::LOG
@ LOG
Definition: GetLCDefs.h:21
CaloClusterContainer.h
GetLCDefs::LIN
@ LIN
Definition: GetLCDefs.h:21
GetLCClassification::m_NormalizationTypeNumber
int m_NormalizationTypeNumber
Definition: GetLCClassification.h:125
GetLCClassification.h
GetLCClassification::m_isampmap
std::vector< int > m_isampmap
Vector of indices in m_dimensions.
Definition: GetLCClassification.h:87
ReadHandle.h
Handle class for reading from StoreGate.
python.handimod.cc
int cc
Definition: handimod.py:523