ATLAS Offline Software
GetLCClassification.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2021 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 //-----------------------
22 #include "CaloLocalHadCalib/GetLCClassification.h"
23 
24 //---------------
25 // C++ Headers --
26 //---------------
27 #include "CaloLocalHadCalib/GetLCDefs.h"
28 #include "xAODCaloEvent/CaloClusterContainer.h"
29 #include "StoreGate/ReadHandle.h"
30 
31 #include "AthenaKernel/errorcheck.h"
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 //###############################################################################
42 GetLCClassification::GetLCClassification(const std::string& name,
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 
84 GetLCClassification::~GetLCClassification()
85 { }
86 
87 //###############################################################################
88 
89 StatusCode GetLCClassification::initialize()
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" );
98  m_NormalizationTypeNumber = GetLCDefs::LIN;
99  }
100  else if ( m_NormalizationType == "Log" ) {
101  ATH_MSG_INFO( "Using weighting proportional to log(E_calib)" );
102  m_NormalizationTypeNumber = GetLCDefs::LOG;
103  }
104  else if ( m_NormalizationType == "NClus" ) {
105  ATH_MSG_INFO( "Using weighting proportional to 1/N_Clus_E_calib>0" );
106  m_NormalizationTypeNumber = GetLCDefs::NCLUS;
107  }
108  else {
109  ATH_MSG_INFO( "Using constant weighting" );
110  m_NormalizationTypeNumber = GetLCDefs::CONST;
111  }
112 
113  if ( m_ClassificationType == "None" ) {
114  ATH_MSG_INFO( "Expecting single particle input" );
115  m_ClassificationTypeNumber = GetLCDefs::NONE;
116  }
117  else if ( m_ClassificationType == "ParticleID_EM" ) {
118  ATH_MSG_INFO( "Expecting ParticleID simulation as input -- use EM type clusters only" );
119  m_ClassificationTypeNumber = GetLCDefs::PARTICLEID_EM;
120  }
121  else if ( m_ClassificationType == "ParticleID_HAD" ) {
122  ATH_MSG_INFO( "Expecting ParticleID simulation as input -- use HAD type clusters only" );
123  m_ClassificationTypeNumber = GetLCDefs::PARTICLEID_HAD;
124  }
125  else {
126  ATH_MSG_WARNING( " unknown classification type " << m_ClassificationType << " given! Using None instead" );
127  m_ClassificationTypeNumber = GetLCDefs::NONE;
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 = (iside>=0?m_dimensions[iside].bins():1);
165  int neta = (ieta>=0?m_dimensions[ieta].bins():1);
166  int nphi = (iphi>=0?m_dimensions[iphi].bins():1);
167  int nlogE = m_dimensions[ilogE].bins();
168  m_hclus.resize(nside*neta*nphi*nlogE,nullptr);
169  for ( int jside=0;jside<nside;jside++) {
170  for ( int jeta=0;jeta<neta;jeta++) {
171  for ( int jphi=0;jphi<nphi;jphi++) {
172  for(int jlogE=0;jlogE<nlogE;jlogE++) {
173  TString hname("hclus");
174  hname += "_iside_";
175  hname += jside;
176  hname += "_[";
177  hname += (iside>=0?m_dimensions[iside].lowEdge():-1);
178  hname += ",";
179  hname += (iside>=0?m_dimensions[iside].highEdge():-1);
180  hname += ",";
181  hname += nside;
182  hname += "]";
183  hname += "_ieta_";
184  hname += jeta;
185  hname += "_[";
186  hname += (ieta>=0?m_dimensions[ieta].lowEdge():-1);
187  hname += ",";
188  hname += (ieta>=0?m_dimensions[ieta].highEdge():-1);
189  hname += ",";
190  hname += neta;
191  hname += "]";
192  hname += "_iphi_";
193  hname += jphi;
194  hname += "_[";
195  hname += (iphi>=0?m_dimensions[iphi].lowEdge():-1);
196  hname += ",";
197  hname += (iphi>=0?m_dimensions[iphi].highEdge():-1);
198  hname += ",";
199  hname += nphi;
200  hname += "]";
201  hname += "_ilogE_";
202  hname += jlogE;
203  hname += "_[";
204  hname += m_dimensions[ilogE].lowEdge();
205  hname += ",";
206  hname += m_dimensions[ilogE].highEdge();
207  hname += ",";
208  hname += nlogE;
209  hname += "]";
210  int iH = jlogE*nphi*neta*nside+jphi*neta*nside+jeta*nside+jside;
211  m_hclus[iH] = new TH2F(hname,hname,
212  m_dimensions[ilogrho].bins(),
213  m_dimensions[ilogrho].lowEdge(),
214  m_dimensions[ilogrho].highEdge(),
215  m_dimensions[iloglambda].bins(),
216  m_dimensions[iloglambda].lowEdge(),
217  m_dimensions[iloglambda].highEdge());
218  m_hclus[iH]->Sumw2();
219  m_hclus[iH]->SetXTitle("log10(<#rho_{cell}> (MeV/mm^{3})) - log10(E_{clus} (MeV))");
220  m_hclus[iH]->SetYTitle("log10(#lambda_{clus} (mm))");
221  m_hclus[iH]->SetZTitle("Number of Clusters");
222  }
223  }
224  }
225  }
226 
227  ATH_CHECK( m_clusterCollName.initialize() );
228 
229  return StatusCode::SUCCESS;
230 }
231 
232 //###############################################################################
233 
234 StatusCode GetLCClassification::finalize()
235 {
236  ATH_MSG_INFO( "Writing out histograms" );
237  m_outputFile->cd();
238  for(unsigned int i=0;i<m_hclus.size();i++) {
239  m_hclus[i]->Write();
240  }
241  m_outputFile->Close();
242 
243  return StatusCode::SUCCESS;
244 }
245 
246 //###############################################################################
247 
248 StatusCode GetLCClassification::execute()
249 {
250  SG::ReadHandle<xAOD::CaloClusterContainer> cc (m_clusterCollName);
251 
252  // total calib hit energy of all clusters
253  double eCalibTot(0.);
254  double nClusECalibGt0 = 0.0;
255 
256  for (const xAOD::CaloCluster * theCluster : *cc) {
257  double eC=999;
258  if (!theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_TOT,eC)) {
259  ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_TOT" );
260  return StatusCode::FAILURE;
261  }
262  if ( m_ClassificationTypeNumber != GetLCDefs::NONE ) {
263  double emFrac = 0;
264  if (!theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_FRAC_EM,emFrac)) {
265  ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_FRAC_EM" );
266  return StatusCode::FAILURE;
267  }
268  if (m_ClassificationTypeNumber == GetLCDefs::PARTICLEID_EM && emFrac < 0.5 )
269  eC = 0;
270  if (m_ClassificationTypeNumber == GetLCDefs::PARTICLEID_HAD && emFrac > 0.5 )
271  eC = 0;
272  }
273  eCalibTot += eC;
274  if ( eC > 0 ) {
275  nClusECalibGt0++;
276  }
277  }
278 
279  if ( eCalibTot > 0 ) {
280  const double inv_eCalibTot = 1. / eCalibTot;
281  const double inv_nClusECalibGt0 = 1. / nClusECalibGt0;
282  for (const xAOD::CaloCluster * pClus : *cc) {
283  double eng = pClus->e();
284  if ( eng > 0 ) {
285  if ( m_ClassificationTypeNumber != GetLCDefs::NONE ) {
286  double emFrac = 0;
287  if (pClus->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_FRAC_EM,emFrac)) {
288  ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_FRAC_EM" );
289  return StatusCode::FAILURE;
290  }
291  if (m_ClassificationTypeNumber == GetLCDefs::PARTICLEID_EM && emFrac < 0.5 )
292  continue;
293  if (m_ClassificationTypeNumber == GetLCDefs::PARTICLEID_HAD && emFrac > 0.5 )
294  continue;
295  }
296 
297  double eta = fabs(pClus->eta());
298 
299  int iside = 0;
300  int ieta = 0;
301  int iphi = 0;
302  int ilogE = 0;
303  int nside = 1;
304  int neta = 1;
305  int nphi = 1;
306  int nlogE = 1;
307 
308  if ( m_isampmap[0] >= 0 ) {
309  const Gaudi::Histo1DDef & hd = m_dimensions[m_isampmap[0]];
310  nside = hd.bins();
311  iside = (int)(nside*(((pClus->eta()<0?-1.0:1.0) - hd.lowEdge())
312  /(hd.highEdge()-hd.lowEdge())));
313  if ( iside < 0 || iside > nside-1 ) {
314  ATH_MSG_WARNING( " Side index out of bounds " <<
315  iside << " not in [0," << nside-1 << "]" );
316  iside = -1;
317  }
318  }
319 
320  if ( m_isampmap[1] >= 0 ) {
321  const Gaudi::Histo1DDef & hd = m_dimensions[m_isampmap[1]];
322  neta = hd.bins();
323  ieta = (int)(neta*((eta - hd.lowEdge())
324  /(hd.highEdge()-hd.lowEdge())));
325  if ( ieta < 0 || ieta > neta-1 ) {
326  ATH_MSG_WARNING( " Eta index out of bounds " <<
327  ieta << " not in [0," << neta-1 << "]" );
328  ieta = -1;
329  }
330  }
331  if ( m_isampmap[2] >= 0 ) {
332  const Gaudi::Histo1DDef & hd = m_dimensions[m_isampmap[2]];
333  nphi = hd.bins();
334  iphi = (int)(nphi*((pClus->phi() - hd.lowEdge())
335  /(hd.highEdge()-hd.lowEdge())));
336  if ( iphi < 0 || iphi > nphi-1 ) {
337  ATH_MSG_WARNING( " Phi index out of bounds " <<
338  iphi << " not in [0," << nphi-1 << "]" );
339  iphi = -1;
340  }
341  }
342 
343  const Gaudi::Histo1DDef & hd = m_dimensions[m_isampmap[3]];
344  nlogE = hd.bins();
345  ilogE = (int)(nlogE*((log10(eng) - hd.lowEdge())
346  /(hd.highEdge()-hd.lowEdge())));
347  if ( ilogE >= 0 && ilogE < nlogE ) {
348  double dens=0,lamb=0,ecal=0;
349  if (!pClus->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_TOT,ecal)) {
350  ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_TOT" );
351  return StatusCode::FAILURE;
352  }
353  if (!pClus->retrieveMoment(xAOD::CaloCluster::FIRST_ENG_DENS,dens)) {
354  ATH_MSG_ERROR( "Failed to retrieve cluster moment FIRST_ENG_DENS" );
355  return StatusCode::FAILURE;
356  }
357  if (!pClus->retrieveMoment(xAOD::CaloCluster::CENTER_LAMBDA,lamb)) {
358  ATH_MSG_ERROR( "Failed to retrieve cluster moment CENTER_LAMBDA" );
359  return StatusCode::FAILURE;
360  }
361  if ( dens > 0 &&
362  lamb > 0 &&
363  ecal > 0 )
364  {
365  int iH = ilogE*nphi*neta*nside+iphi*neta*nside+ieta*nside+iside;
366  if ( m_hclus[iH]) {
367  double norm = 0.0;
368  if ( m_NormalizationTypeNumber == GetLCDefs::LIN ) {
369  norm = ecal*inv_eCalibTot;
370  }
371  else if ( m_NormalizationTypeNumber == GetLCDefs::LOG ) {
372  // cluster has to have at least 1% of the calib hit E
373  norm = log10(ecal*inv_eCalibTot)+2.0;
374  }
375  else if ( m_NormalizationTypeNumber == GetLCDefs::NCLUS ) {
376  norm = inv_nClusECalibGt0;
377  }
378  else {
379  norm = 1.0;
380  }
381  if ( norm > 0 ) {
382  m_hclus[iH]->Fill(log10(dens)-log10(eng),log10(lamb),norm);
383  }
384  }
385  }
386  }
387  }
388  }
389  }
390 
391  return StatusCode::SUCCESS;
392 }
393 
394 void GetLCClassification::mapinsert(const std::vector<Gaudi::Histo1DDef> & dims) {
395  for (unsigned int i=0;i<dims.size();i++) {
396  m_dimensionsmap[dims[i].title()] = dims[i];
397  }
398 }
399 
400 void GetLCClassification::mapparse() {
401 
402  for (std::pair<const std::string, Gaudi::Histo1DDef>& p : m_dimensionsmap) {
403  m_dimensions.push_back(p.second);
404  ATH_MSG_DEBUG( " New Dimension: "
405  << p.second.title() << ", [" << p.second.lowEdge()
406  << ", " << p.second.highEdge()
407  << ", " << p.second.bins()
408  << "]" );
409  }
410 }
411 
412 //###############################################################################
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
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
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
python.App.bins
bins
Definition: App.py:410
dqt_zlumi_pandas.hname
string hname
Definition: dqt_zlumi_pandas.py:272
eta
Scalar eta() const
pseudorapidity method
Definition: AmgMatrixBasePlugin.h:79
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:394
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
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
GetLCClassification::finalize
virtual StatusCode finalize()
Definition: GetLCClassification.cxx:234
lumiFormat.i
int i
Definition: lumiFormat.py:92
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:400
GetLCClassification::execute
virtual StatusCode execute()
Definition: GetLCClassification.cxx:248
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
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:195
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
Generated on Fri Jul 5 2024 21:16:27 for ATLAS Offline Software by   doxygen 1.8.18