ATLAS Offline Software
GetLCWeights.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: GetLCWeights.cxx,v 1.3 2009-03-06 14:43:23 pospelov Exp $
8 //
9 // Description: see GetLCWeights.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/GetLCWeights.h"
23 
24 //---------------
25 // C++ Headers --
26 //---------------
27 #include "CaloLocalHadCalib/GetLCDefs.h"
28 #include "xAODCaloEvent/CaloClusterContainer.h"
29 #include "CaloEvent/CaloCell.h"
30 #include "CaloUtils/CaloSamplingHelper.h"
31 #include "CaloSimEvent/CaloCalibrationHit.h"
32 #include "CaloSimEvent/CaloCalibrationHitContainer.h"
33 #include "CaloIdentifier/CaloCell_ID.h"
34 #include "StoreGate/ReadHandle.h"
35 
36 #include "AthenaKernel/errorcheck.h"
37 
38 #include "TFile.h"
39 #include "TProfile2D.h"
40 #include "TString.h"
41 #include <iterator>
42 #include <cmath>
43 
44 //###############################################################################
45 GetLCWeights::GetLCWeights(const std::string& name,
46  ISvcLocator* pSvcLocator)
47  : AthAlgorithm(name, pSvcLocator),
48  m_outputFile(nullptr),
49  m_clusterCollName("CaloTopoClusters"),
50  m_useInversionMethod(true),
51  m_NormalizationType("Lin"),
52  m_NormalizationTypeNumber(0),
53  m_ClassificationType("None"),
54  m_ClassificationTypeNumber(0)
55 {
56 
57  std::vector<Gaudi::Histo1DDef> dims(7);
58  dims[1] = Gaudi::Histo1DDef("side",-1.5,1.5,1);
59  dims[2] = Gaudi::Histo1DDef("|eta|",0.,1.6,16);
60  dims[3] = Gaudi::Histo1DDef("phi",-M_PI,M_PI,1);
61  dims[4] = Gaudi::Histo1DDef("log10(E_clus (MeV))",log10(200),log10(1e6),14);
62  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-4.6,1.1,14);
63  dims[6] = Gaudi::Histo1DDef("weight",-2,3,1);
64 
65  dims[0] = Gaudi::Histo1DDef("EMB1",0.5,1.5,1);
66  mapinsert(dims);
67  dims[0] = Gaudi::Histo1DDef("EMB2",1.5,2.5,1);
68  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-4.8,0.5,14);
69  mapinsert(dims);
70  dims[0] = Gaudi::Histo1DDef("EMB3",2.5,3.5,1);
71  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-4.9,0.6,14);
72  mapinsert(dims);
73  dims[0] = Gaudi::Histo1DDef("EME1",4.5,5.5,1);
74  dims[2] = Gaudi::Histo1DDef("|eta|",1.2,3.2,20);
75  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-4.5,1.7,14);
76  mapinsert(dims);
77  dims[0] = Gaudi::Histo1DDef("EME2",5.5,6.5,1);
78  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-4.8,1.2,14);
79  mapinsert(dims);
80  dims[0] = Gaudi::Histo1DDef("EME3",6.5,7.5,1);
81  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-4.5,1.0,14);
82  mapinsert(dims);
83  dims[0] = Gaudi::Histo1DDef("HEC0",7.5,8.5,1);
84  dims[2] = Gaudi::Histo1DDef("|eta|",1.4,3.4,20);
85  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-4.8,0.6,14);
86  mapinsert(dims);
87  dims[0] = Gaudi::Histo1DDef("HEC1",8.5,9.5,1);
88  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-4.9,-0.2,14);
89  mapinsert(dims);
90  dims[0] = Gaudi::Histo1DDef("HEC2",9.5,10.5,1);
91  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-4.8,-0.2,14);
92  mapinsert(dims);
93  dims[0] = Gaudi::Histo1DDef("HEC3",10.5,11.5,1);
94  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-4.7,-0.2,14);
95  mapinsert(dims);
96  dims[0] = Gaudi::Histo1DDef("TileBar0",11.5,12.5,1);
97  dims[2] = Gaudi::Histo1DDef("|eta|",0.,1.2,12);
98  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-6,-1,14);
99  mapinsert(dims);
100  dims[0] = Gaudi::Histo1DDef("TileBar1",12.5,13.5,1);
101  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-6.7,-1.2,14);
102  mapinsert(dims);
103  dims[0] = Gaudi::Histo1DDef("TileBar2",13.5,14.5,1);
104  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-7,-1.9,14);
105  mapinsert(dims);
106  dims[0] = Gaudi::Histo1DDef("TileGap1",14.5,15.5,1);
107  dims[2] = Gaudi::Histo1DDef("|eta|",0.8,1.8,10);
108  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-6,-1,14);
109  mapinsert(dims);
110  dims[0] = Gaudi::Histo1DDef("TileGap2",15.5,16.5,1);
111  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-6.5,-1.5,14);
112  mapinsert(dims);
113  dims[0] = Gaudi::Histo1DDef("TileExt0",17.5,18.5,1);
114  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-6.5,-1,14);
115  mapinsert(dims);
116  dims[0] = Gaudi::Histo1DDef("TileExt1",18.5,19.5,1);
117  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-6.6,-1.5,14);
118  mapinsert(dims);
119  dims[0] = Gaudi::Histo1DDef("TileExt2",19.5,20.5,1);
120  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-7.2,-2.4,14);
121  mapinsert(dims);
122  dims[0] = Gaudi::Histo1DDef("FCal1",20.5,21.5,1);
123  dims[2] = Gaudi::Histo1DDef("|eta|",2.8,5.0,22);
124  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-3.5,2,14);
125  mapinsert(dims);
126  dims[0] = Gaudi::Histo1DDef("FCal2",21.5,22.5,1);
127  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-3.5,1.5,14);
128  mapinsert(dims);
129  dims[0] = Gaudi::Histo1DDef("FCal3",22.5,23.5,1);
130  dims[5] = Gaudi::Histo1DDef("log10(rho_cell (MeV/mm^3))",-3.8,1,14);
131  mapinsert(dims);
132 
133  // Dimensions to use for all samplings
134  declareProperty("SamplingDimensions",m_dimensionsmap);
135 
136  // Name of output file to save histograms in
137  declareProperty("OutputFileName",m_outputFileName);
138 
139  // Name of ClusterContainer to use
140  declareProperty("ClusterCollectionName",m_clusterCollName);
141 
142  // Names of CalibrationHitContainers to use
143  declareProperty("CalibrationHitContainerNames",m_CalibrationHitContainerNames);
144  // Use Inversion Method
145  declareProperty("UseInversionMethod",m_useInversionMethod);
146  // Normalization type "Const", "Lin", "Log", "NClus"
147  declareProperty("NormalizationType",m_NormalizationType);
148  // Classification type "None" for single pion MC or
149  // "ParticleID_EM" for ParticleID based em-type clusters
150  // "ParticleID_HAD" for ParticleID based had-type clusters
151  declareProperty("ClassificationType", m_ClassificationType);
152 }
153 
154 //###############################################################################
155 
156 GetLCWeights::~GetLCWeights()
157 { }
158 
159 //###############################################################################
160 
161 StatusCode GetLCWeights::initialize()
162 {
163 
164  m_outputFile = std::make_unique<TFile>(m_outputFileName.c_str(),"RECREATE");
165  m_outputFile->cd();
166 
167  m_weight.resize(CaloSampling::Unknown);
168  m_isampmap.resize(CaloSampling::Unknown);
169  mapparse();
170 
171 
172  if ( m_NormalizationType == "Lin" ) {
173  ATH_MSG_INFO( "Using weighting proportional to E_calib" );
174  m_NormalizationTypeNumber = GetLCDefs::LIN;
175  }
176  else if ( m_NormalizationType == "Log" ) {
177  ATH_MSG_INFO( "Using weighting proportional to log(E_calib)" );
178  m_NormalizationTypeNumber = GetLCDefs::LOG;
179  }
180  else if ( m_NormalizationType == "NClus" ) {
181  ATH_MSG_INFO( "Using weighting proportional to 1/N_Clus_E_calib>0" );
182  m_NormalizationTypeNumber = GetLCDefs::NCLUS;
183  }
184  else {
185  ATH_MSG_INFO( "Using constant weighting" );
186  m_NormalizationTypeNumber = GetLCDefs::CONST;
187  }
188 
189  if ( m_ClassificationType == "None" ) {
190  ATH_MSG_INFO( "Expecting single particle input" );
191  m_ClassificationTypeNumber = GetLCDefs::NONE;
192  }
193  else if ( m_ClassificationType == "ParticleID_EM" ) {
194  ATH_MSG_INFO( "Expecting ParticleID simulation as input -- use EM type clusters only" );
195  m_ClassificationTypeNumber = GetLCDefs::PARTICLEID_EM;
196  }
197  else if ( m_ClassificationType == "ParticleID_HAD" ) {
198  ATH_MSG_INFO( "Expecting ParticleID simulation as input -- use HAD type clusters only" );
199  m_ClassificationTypeNumber = GetLCDefs::PARTICLEID_HAD;
200  }
201  else {
202  ATH_MSG_WARNING( " unknown classification type " << m_ClassificationType << " given! Using None instead" );
203  m_ClassificationTypeNumber = GetLCDefs::NONE;
204  }
205 
206  for(unsigned int isamp=0;isamp<m_dimensions.size();isamp++) {
207  int theSampling(CaloSampling::Unknown);
208  for (unsigned int jsamp = 0;jsamp< CaloSampling::Unknown; jsamp++) {
209  if ( m_dimensions[isamp][0].title() == CaloSamplingHelper::getSamplingName((CaloSampling::CaloSample)jsamp)) {
210  theSampling = jsamp;
211  break;
212  }
213  }
214  if ( theSampling == CaloSampling::Unknown ) {
215  ATH_MSG_ERROR( "Calorimeter sampling "
216  << m_dimensions[isamp][0].title()
217  << " is not a valid Calorimeter sampling name and will be ignored! "
218  << "Valid names are: ";
219  for (unsigned int jsamp = 0;jsamp< CaloSampling::Unknown; jsamp++) {
220  msg() << CaloSamplingHelper::getSamplingName((CaloSampling::CaloSample)jsamp);
221  if ( jsamp < CaloSampling::Unknown-1)
222  msg() << ", ";
223  else
224  msg() << ".";
225  }
226  msg() );
227  }
228  else {
229  m_isampmap[theSampling].resize(4,-1);
230  m_isampmap[theSampling][0] = isamp;
231  int iside(-1);
232  int ieta(-1);
233  int iphi(-1);
234  int ilogE(-1);
235  int ilogrho(-1);
236  int iweight(-1);
237 
238  for(unsigned int idim=1;idim<m_dimensions[isamp].size();idim++) {
239  if ( m_dimensions[isamp][idim].title() == "side" ) {
240  iside = idim;
241  m_isampmap[theSampling][1] = iside;
242  }
243  else if ( m_dimensions[isamp][idim].title() == "|eta|" ) {
244  ieta = idim;
245  m_isampmap[theSampling][2] = ieta;
246  }
247  else if ( m_dimensions[isamp][idim].title() == "phi" ) {
248  iphi = idim;
249  m_isampmap[theSampling][3] = iphi;
250  }
251  else if ( m_dimensions[isamp][idim].title() == "log10(E_clus (MeV))" )
252  ilogE = idim;
253  else if ( m_dimensions[isamp][idim].title() == "log10(rho_cell (MeV/mm^3))" )
254  ilogrho = idim;
255  else if ( m_dimensions[isamp][idim].title() == "weight" )
256  iweight = idim;
257  }
258  if ( ilogE < 0 || ilogrho < 0 || iweight < 0 ) {
259  ATH_MSG_FATAL( " Mandatory dimension log10E, log10rho or weight missing ..." );
260  return StatusCode::FAILURE;
261  }
262  int nside = (iside>=0?m_dimensions[isamp][iside].bins():1);
263  int neta = (ieta>=0?m_dimensions[isamp][ieta].bins():1);
264  int nphi = (iphi>=0?m_dimensions[isamp][iphi].bins():1);
265  m_weight[theSampling].resize(nside*neta*nphi,nullptr);
266  for ( int jside=0;jside<nside;jside++) {
267  for ( int jeta=0;jeta<neta;jeta++) {
268  for ( int jphi=0;jphi<nphi;jphi++) {
269  TString wname("");
270  if ( m_useInversionMethod )
271  wname += "inv_weight";
272  else
273  wname += "weight";
274  wname += "_isamp_";
275  wname += theSampling;
276  wname += "_iside_";
277  wname += jside;
278  wname += "_[";
279  wname += (iside>=0?m_dimensions[isamp][iside].lowEdge():-1);
280  wname += ",";
281  wname += (iside>=0?m_dimensions[isamp][iside].highEdge():-1);
282  wname += ",";
283  wname += nside;
284  wname += "]";
285  wname += "_ieta_";
286  wname += jeta;
287  wname += "_[";
288  wname += (ieta>=0?m_dimensions[isamp][ieta].lowEdge():-1);
289  wname += ",";
290  wname += (ieta>=0?m_dimensions[isamp][ieta].highEdge():-1);
291  wname += ",";
292  wname += neta;
293  wname += "]";
294  wname += "_iphi_";
295  wname += jphi;
296  wname += "_[";
297  wname += (iphi>=0?m_dimensions[isamp][iphi].lowEdge():-1);
298  wname += ",";
299  wname += (iphi>=0?m_dimensions[isamp][iphi].highEdge():-1);
300  wname += ",";
301  wname += nphi;
302  wname += "]";
303  int iW = jphi*neta*nside+jeta*nside+jside;
304  m_weight[theSampling][iW]=
305  new TProfile2D(wname,wname,
306  m_dimensions[isamp][ilogE].bins(),
307  m_dimensions[isamp][ilogE].lowEdge(),
308  m_dimensions[isamp][ilogE].highEdge(),
309  m_dimensions[isamp][ilogrho].bins(),
310  m_dimensions[isamp][ilogrho].lowEdge(),
311  m_dimensions[isamp][ilogrho].highEdge(),
312  m_dimensions[isamp][iweight].lowEdge(),
313  m_dimensions[isamp][iweight].highEdge(),
314  "spread");
315  if ( m_useInversionMethod ) {
316  m_weight[theSampling][iW]->SetYTitle("log10(#rho_{cell}^{true} (MeV/mm^{3}))");
317  m_weight[theSampling][iW]->SetZTitle("E_{reco}/E_{tot}");
318  }
319  else {
320  m_weight[theSampling][iW]->SetYTitle("log10(#rho_{cell} (MeV/mm^{3}))");
321  m_weight[theSampling][iW]->SetZTitle("E_{tot}/E_{reco}");
322  }
323  m_weight[theSampling][iW]->SetXTitle("log10(E_{clus} (MeV))");
324  }
325  }
326  }
327  }
328  }
329 
330  ATH_CHECK( m_clusterCollName.initialize() );
331  ATH_CHECK( m_CalibrationHitContainerNames.initialize() );
332 
333  return StatusCode::SUCCESS;
334 }
335 
336 //###############################################################################
337 
338 StatusCode GetLCWeights::finalize()
339 {
340  ATH_MSG_INFO( "Writing out histograms" );
341  m_outputFile->cd();
342  for(unsigned int i=0;i<m_weight.size();i++) {
343  for(unsigned int j=0;j<m_weight[i].size();j++) {
344  if ( m_weight[i][j] )
345  m_weight[i][j]->Write();
346  }
347  }
348  m_outputFile->Close();
349 
350  return StatusCode::SUCCESS;
351 }
352 
353 //###############################################################################
354 
355 StatusCode GetLCWeights::execute()
356 {
357  const EventContext& ctx = getContext();
358  SG::ReadHandle<xAOD::CaloClusterContainer> cc (m_clusterCollName, ctx);
359 
360  std::vector<const CaloCalibrationHitContainer *> v_cchc;
361  for (const SG::ReadHandleKey<CaloCalibrationHitContainer>& k : m_CalibrationHitContainerNames) {
362  SG::ReadHandle<CaloCalibrationHitContainer> cchc (k, ctx);
363  v_cchc.push_back(cchc.cptr());
364  }
365 
366  std::vector<ClusWeight *> cellVector[CaloCell_ID::NSUBCALO];
367 
368  const CaloCell_ID* calo_id = nullptr;
369  ATH_CHECK(detStore()->retrieve(calo_id,"CaloCell_ID"));
370 
371  for(int ic=0;ic<CaloCell_ID::NSUBCALO; ic++) {
372  unsigned int maxHashSize(0);
373  IdentifierHash myHashMin,myHashMax;
374  calo_id->calo_cell_hash_range (ic,myHashMin,myHashMax);
375  maxHashSize = myHashMax-myHashMin;
376  cellVector[ic].resize(maxHashSize,nullptr);
377  }
378 
379  // loop over all cell members of all clusters and fill cell vector
380  // for used cells
381 
382  // total calib hit energy of all clusters
383  double eCalibTot(0.);
384 
385  unsigned int iClus = 0;
386  double nClusECalibGt0 = 0.0;
387  for (const xAOD::CaloCluster* theCluster : *cc) {
388  double eC=999;
389  if (!theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_TOT,eC)) {
390  ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_TOT");
391  return StatusCode::FAILURE;
392  }
393  if ( m_ClassificationTypeNumber != GetLCDefs::NONE ) {
394  double emFrac=-999;
395  if (!theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_FRAC_EM,emFrac)){
396  ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_FAC_EM");
397  return StatusCode::FAILURE;
398  }
399  if (m_ClassificationTypeNumber == GetLCDefs::PARTICLEID_EM && emFrac < 0.5 )
400  eC = 0;
401  if (m_ClassificationTypeNumber == GetLCDefs::PARTICLEID_HAD && emFrac > 0.5 )
402  eC = 0;
403  }
404  eCalibTot += eC;
405  if ( eC > 0 ) {
406  nClusECalibGt0++;
407  }
408  xAOD::CaloCluster::const_cell_iterator cellIter = theCluster->cell_begin();
409  xAOD::CaloCluster::const_cell_iterator cellIterEnd = theCluster->cell_end();
410  for(; cellIter != cellIterEnd; cellIter++ ){
411  const CaloCell* pCell = (*cellIter);
412  Identifier myId = pCell->ID();
413  IdentifierHash myHashId;
414  int otherSubDet(0);
415  myHashId = calo_id->subcalo_cell_hash(myId,otherSubDet);
416  ClusWeight * myClus = new ClusWeight();
417  myClus->iClus = iClus;
418  myClus->weight = cellIter.weight();
419  myClus->next = nullptr;
420  myClus->eCalibTot = 0;
421  ClusWeight * theList = cellVector[otherSubDet][(unsigned int)myHashId];
422  if ( theList ) {
423  while ( theList->next )
424  theList = theList->next;
425  theList->next = myClus;
426  }
427  else {
428  cellVector[otherSubDet][(unsigned int)myHashId] = myClus;
429  }
430  }
431  }
432 
433  // calculate total calib energy of each cell in each cluster
434 
435  for (const CaloCalibrationHitContainer* cchc : v_cchc) {
436  for (const CaloCalibrationHit* hit : *cchc) {
437  Identifier myId = hit->cellID();
438  int otherSubDet;
439  IdentifierHash myHashId = calo_id->subcalo_cell_hash(myId,otherSubDet);
440  if ( myHashId != CaloCell_ID::NOT_VALID ) {
441  ClusWeight * theList = cellVector[otherSubDet][(unsigned int)myHashId];
442  while ( theList ) {
443  theList->eCalibTot += hit->energyTotal();
444  theList = theList->next;
445  }
446  }
447  }
448  }
449 
450  // fill weight histos for all cells in all clusters
451  // since pions can be split in several clusters the weights are
452  // calculated by importance of the cluster for the current pion - i.e.
453  // the weights get a weight itself proportinal to the calib hit energy
454  // sum of the cluster over the total calib hit energy
455 
456  if ( eCalibTot > 0 ) {
457  const double inv_eCalibTot = 1. / eCalibTot;
458  const double inv_nClusECalibGt0 = 1. / nClusECalibGt0;
459  for (unsigned int j=0;j<cc->size();j++) {
460  const xAOD::CaloCluster * pClus = cc->at(j);
461  double eng = pClus->e();
462  double eCalib=-999;
463  if (!pClus->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_TOT,eCalib)) {
464  ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_TOT");
465  return StatusCode::FAILURE;
466  }
467  if ( eng > 0 && eCalib > 0 ) {
468  if ( m_ClassificationTypeNumber != GetLCDefs::NONE ) {
469  double emFrac=-999;
470  if (!pClus->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_FRAC_EM,emFrac)) {
471  ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_FAC_EM");
472  return StatusCode::FAILURE;
473  }
474  if (m_ClassificationTypeNumber == GetLCDefs::PARTICLEID_EM && emFrac < 0.5 )
475  continue;
476  if (m_ClassificationTypeNumber == GetLCDefs::PARTICLEID_HAD && emFrac > 0.5 )
477  continue;
478  }
479  xAOD::CaloCluster::const_cell_iterator cellIter = pClus->cell_begin();
480  xAOD::CaloCluster::const_cell_iterator cellIterEnd = pClus->cell_end();
481  for(; cellIter != cellIterEnd; cellIter++ ){
482  const CaloCell* pCell = (*cellIter);
483  const Identifier myId = pCell->ID();
484  const CaloDetDescrElement* myCDDE=pCell->caloDDE();
485  const int caloSample = myCDDE->getSampling();//m_calo_id->calo_sample(myId);
486  if ( !m_isampmap[caloSample].empty() &&
487  !m_weight[caloSample].empty() ) {
488  int isideCell = 0;
489  int ietaCell = 0;
490  int iphiCell = 0;
491  int nside = 1;
492  int neta = 1;
493  int nphi = 1;
494  if ( m_isampmap[caloSample][1] >= 0 ) {
495  const Gaudi::Histo1DDef & hd = m_dimensions[m_isampmap[caloSample][0]][m_isampmap[caloSample][1]];
496  nside = hd.bins();
497  isideCell = (int)(nside*(((pCell->eta()<0?-1.0:1.0) - hd.lowEdge())
498  /(hd.highEdge()-hd.lowEdge())));
499  if ( isideCell < 0 || isideCell > nside-1 ) {
500  ATH_MSG_WARNING( " Side index out of bounds " <<
501  isideCell << " not in [0," << nside-1 << "] for "
502  << "Sampl=" << caloSample );
503  isideCell = -1;
504  }
505  }
506  if ( m_isampmap[caloSample][2] >= 0 ) {
507  const Gaudi::Histo1DDef & hd = m_dimensions[m_isampmap[caloSample][0]][m_isampmap[caloSample][2]];
508  neta = hd.bins();
509  ietaCell = (int)(neta*((fabs(pCell->eta()) - hd.lowEdge())
510  /(hd.highEdge()-hd.lowEdge())));
511  if ( ietaCell < 0 || ietaCell > neta-1 ) {
512  ATH_MSG_WARNING( " Eta index out of bounds " <<
513  ietaCell << " not in [0," << neta-1 << "] for "
514  << "Sampl=" << caloSample );
515  ietaCell = -1;
516  }
517  }
518  if ( m_isampmap[caloSample][3] >= 0 ) {
519  const Gaudi::Histo1DDef & hd = m_dimensions[m_isampmap[caloSample][0]][m_isampmap[caloSample][3]];
520  nphi = hd.bins();
521  iphiCell = (int)(nphi*((pCell->phi() - hd.lowEdge())
522  /(hd.highEdge()-hd.lowEdge())));
523  if ( iphiCell < 0 || iphiCell > nphi-1 ) {
524  ATH_MSG_WARNING( " Phi index out of bounds " <<
525  iphiCell << " not in [0," << nphi-1 << "] for "
526  << "Sampl=" << caloSample );
527  iphiCell = -1;
528  }
529  }
530  if ( isideCell >=0 && ietaCell >=0 && iphiCell >= 0 ) {
531  if ( myCDDE->volume() > 0 ) {
532  IdentifierHash myHashId;
533  int otherSubDet(0);
534  myHashId = calo_id->subcalo_cell_hash(myId,otherSubDet);
535  unsigned int iW = iphiCell*neta*nside+ietaCell*nside+isideCell;
536  if ( iW >= m_weight[caloSample].size() ) {
537  ATH_MSG_WARNING( " Index out of bounds " <<
538  iW << " > " << m_weight[caloSample].size()-1 << " for "
539  << "Sampl=" << caloSample
540  << ", iphi=" << iphiCell
541  << ", ieta=" << ietaCell
542  << ", iside=" << isideCell );
543  }
544  else {
545  ClusWeight * theList = cellVector[otherSubDet][(unsigned int)myHashId];
546  while ( theList && theList->iClus != j )
547  theList = theList->next;
548 
549  if (m_weight[caloSample][iW] && theList && eCalibTot > 0) {
550  double norm = 0.0;
551  if ( m_NormalizationTypeNumber == GetLCDefs::LIN ) {
552  norm = eCalib * inv_eCalibTot;
553  }
554  else if ( m_NormalizationTypeNumber == GetLCDefs::LOG ) {
555  if ( eCalib > 0 ) {
556  // cluster has to have at least 1% of the calib hit E
557  norm = log10(eCalib * inv_eCalibTot)+2.0;
558  }
559  }
560  else if ( m_NormalizationTypeNumber == GetLCDefs::NCLUS ) {
561  if ( eCalib > 0 ) {
562  norm = inv_nClusECalibGt0;
563  }
564  }
565  else {
566  norm = 1.0;
567  }
568  if ( norm > 0 ) {
569  if ( !m_useInversionMethod && pCell->e()>0 ) {
570  m_weight[caloSample][iW]->Fill(log10(eng),
571  log10(pCell->e()/myCDDE->volume()),
572  theList->eCalibTot/pCell->e(),norm);
573 
574  }
575  else if (m_useInversionMethod && theList->eCalibTot>0 ) {
576  m_weight[caloSample][iW]->Fill(log10(eng),
577  log10(theList->eCalibTot/myCDDE->volume()),
578  pCell->e()/theList->eCalibTot,
579  norm);
580  }
581  }
582  }
583  }
584  }
585  }
586  }
587  }
588  }
589  }
590  }
591 
592  for(unsigned int ic=0;ic<CaloCell_ID::NSUBCALO; ic++) {
593  for (unsigned int ii = 0;ii<cellVector[ic].size();ii++ ) {
594  ClusWeight * theList = cellVector[ic][ii];
595  ClusWeight * prev = nullptr;
596  while ( theList) {
597  while ( theList->next ) {
598  prev = theList;
599  theList = theList->next;
600  }
601  delete theList;
602  if ( prev )
603  prev->next = nullptr;
604  else
605  cellVector[ic][ii] = nullptr;
606  theList = cellVector[ic][ii];
607  prev = nullptr;
608  }
609  }
610  }
611 
612  return StatusCode::SUCCESS;
613 }
614 
615 void GetLCWeights::mapinsert(const std::vector<Gaudi::Histo1DDef> & dims) {
616  for (unsigned int i=0;i<dims.size();i++) {
617  m_dimensionsmap[dims[0].title()+":"+dims[i].title()] = dims[i];
618  }
619 }
620 
621 void GetLCWeights::mapparse() {
622 
623  std::vector<int> theUsedSamplings(CaloSampling::Unknown,-1);
624 
625  int nsamp(-1);
626 
627  for (const std::pair<const std::string, Gaudi::Histo1DDef>& p : m_dimensionsmap) {
628  std::string_view dimname = std::string_view(p.first).substr(0,p.first.find(':'));
629  int theSampling(CaloSampling::Unknown);
630  for (unsigned int jsamp = 0;jsamp< CaloSampling::Unknown; jsamp++) {
631  if ( dimname == CaloSamplingHelper::getSamplingName((CaloSampling::CaloSample)jsamp)) {
632  theSampling = jsamp;
633  break;
634  }
635  }
636  if ( theSampling == CaloSampling::Unknown ) {
637  msg(MSG::ERROR) << "Calorimeter sampling " << dimname
638  << " is not a valid Calorimeter sampling name and will be ignored! "
639  << "Valid names are: ";
640  for (unsigned int jsamp = 0;jsamp< CaloSampling::Unknown; jsamp++) {
641  msg() << CaloSamplingHelper::getSamplingName((CaloSampling::CaloSample)jsamp);
642  if ( jsamp < CaloSampling::Unknown-1)
643  msg() << ", ";
644  else
645  msg() << ".";
646  }
647  msg() << endmsg;
648  }
649  else {
650  if ( theUsedSamplings[theSampling] == -1 ) {
651  nsamp++;
652  theUsedSamplings[theSampling] = nsamp;
653  m_dimensions.resize(nsamp+1);
654  m_dimensions[nsamp].resize(0);
655  }
656  m_dimensions[theUsedSamplings[theSampling]].push_back(p.second);
657  ATH_MSG_DEBUG(" New Dimension for " << dimname << ": "
658  << p.second.title() << ", [" << p.second.lowEdge()
659  << ", " << p.second.highEdge()
660  << ", " << p.second.bins()
661  << "]");
662  }
663  }
664 }
665 
666 //###############################################################################
python.PyKernel.retrieve
def retrieve(aClass, aKey=None)
Definition: PyKernel.py:110
GetLCDefs::CONST
@ CONST
Definition: GetLCDefs.h:21
GetLCWeights::m_dimensions
std::vector< std::vector< Gaudi::Histo1DDef > > m_dimensions
definition of all dimensions used for each sampling
Definition: GetLCWeights.h:92
CaloCalibrationHitContainer
Definition: CaloCalibrationHitContainer.h:25
GetLCDefs::Unknown
@ Unknown
Definition: GetLCDefs.h:21
PlotCalibFromCool.norm
norm
Definition: PlotCalibFromCool.py:100
xAOD::CaloCluster_v1::cell_begin
const_cell_iterator cell_begin() const
Iterator of the underlying CaloClusterCellLink (const version)
Definition: CaloCluster_v1.h:812
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition: AthMsgStreamMacros.h:34
CaloCalibrationHit.h
CaloCell::phi
virtual double phi() const override final
get phi (through CaloDetDescrElement)
Definition: CaloCell.h:359
python.PerfMonSerializer.p
def p
Definition: PerfMonSerializer.py:743
GetLCWeights::execute
virtual StatusCode execute()
Definition: GetLCWeights.cxx:355
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
SG::ReadHandle::cptr
const_pointer_type cptr()
Dereference the pointer.
CaloCellPos2Ntuple.int
int
Definition: CaloCellPos2Ntuple.py:24
python.App.bins
bins
Definition: App.py:410
ClusWeight::weight
double weight
Definition: GetLCWeights.h:36
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
CaloDetDescrElement
This class groups all DetDescr information related to a CaloCell. Provides a generic interface for al...
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:66
GetLCDefs.h
GetLCWeights::m_useInversionMethod
bool m_useInversionMethod
flag to switch on/off the use of the inversion method
Definition: GetLCWeights.h:150
ClusWeight::eCalibTot
double eCalibTot
Definition: GetLCWeights.h:37
GetLCWeights::mapinsert
void mapinsert(const std::vector< Gaudi::Histo1DDef > &dims)
Definition: GetLCWeights.cxx:615
TProfile2D
Definition: rootspy.cxx:531
M_PI
#define M_PI
Definition: ActiveFraction.h:11
GetLCWeights::GetLCWeights
GetLCWeights(const std::string &name, ISvcLocator *pSvcLocator)
Definition: GetLCWeights.cxx:45
CaloCell.h
CaloCell::e
virtual double e() const override final
get energy (data member) (synonym to method energy()
Definition: CaloCell.h:317
GetLCWeights::m_outputFileName
std::string m_outputFileName
Name of the output file to save histograms in.
Definition: GetLCWeights.h:120
SG::ReadHandleKey< CaloCalibrationHitContainer >
ClusWeight::iClus
unsigned int iClus
Definition: GetLCWeights.h:35
empty
bool empty(TH1 *h)
Definition: computils.cxx:294
GetLCWeights::m_NormalizationType
std::string m_NormalizationType
string to choose different normalization types
Definition: GetLCWeights.h:162
CaloCalibrationHitContainer.h
CaloClusterCellLink::const_iterator::weight
weight_t weight() const
Accessor for weight associated to this cell.
Definition: CaloClusterCellLink.h:82
CaloCell_ID.h
AthCommonDataStore< AthCommonMsg< Algorithm > >::detStore
const ServiceHandle< StoreGateSvc > & detStore() const
The standard StoreGateSvc/DetectorStore Returns (kind of) a pointer to the StoreGateSvc.
Definition: AthCommonDataStore.h:95
GetLCDefs::NONE
@ NONE
Definition: GetLCDefs.h:25
GetLCDefs::NCLUS
@ NCLUS
Definition: GetLCDefs.h:21
GetLCDefs::PARTICLEID_EM
@ PARTICLEID_EM
Definition: GetLCDefs.h:25
python.setupRTTAlg.size
int size
Definition: setupRTTAlg.py:39
GetLCWeights::finalize
virtual StatusCode finalize()
Definition: GetLCWeights.cxx:338
xAOD::CaloCluster_v1
Description of a calorimeter cluster.
Definition: CaloCluster_v1.h:59
GetLCWeights::m_NormalizationTypeNumber
int m_NormalizationTypeNumber
Definition: GetLCWeights.h:164
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
lumiFormat.i
int i
Definition: lumiFormat.py:92
CaloSampling::CaloSample
CaloSample
Definition: Calorimeter/CaloGeoHelpers/CaloGeoHelpers/CaloSampling.h:22
CaloCell_Base_ID::calo_cell_hash_range
void calo_cell_hash_range(const Identifier id, IdentifierHash &caloCellMin, IdentifierHash &caloCellMax) const
to loop on 'global' cell hashes of one sub-calorimeter alone
Identifier
Definition: DetectorDescription/Identifier/Identifier/Identifier.h:32
endmsg
#define endmsg
Definition: AnalysisConfig_Ntuple.cxx:63
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
CaloCell::caloDDE
const CaloDetDescrElement * caloDDE() const
get pointer to CaloDetDescrElement (data member)
Definition: CaloCell.h:305
xAOD::CaloCluster_v1::ENG_CALIB_TOT
@ ENG_CALIB_TOT
Calibration Hit energy inside the cluster.
Definition: CaloCluster_v1.h:195
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
CaloSamplingHelper.h
CaloCell_Base_ID::subcalo_cell_hash
IdentifierHash subcalo_cell_hash(const Identifier cellId, int &subCalo) const
create hash id from 'global' cell id
xAOD::CaloCluster_v1::retrieveMoment
bool retrieveMoment(MomentType type, double &value) const
Retrieve individual moment.
Definition: CaloCluster_v1.cxx:738
CaloCell_ID
Helper class for offline cell identifiers.
Definition: CaloCell_ID.h:34
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
GetLCWeights::m_outputFile
std::unique_ptr< TFile > m_outputFile
Output file to save histograms in.
Definition: GetLCWeights.h:126
GetLCWeights::m_CalibrationHitContainerNames
SG::ReadHandleKeyArray< CaloCalibrationHitContainer > m_CalibrationHitContainerNames
vector of calibration hit container names to use.
Definition: GetLCWeights.h:138
AthAlgorithm
Definition: AthAlgorithm.h:47
grepfile.ic
int ic
Definition: grepfile.py:33
CaloCalibrationHit
Class to store calorimeter calibration hit.
Definition: CaloCalibrationHit.h:21
GetLCWeights::m_ClassificationTypeNumber
int m_ClassificationTypeNumber
Definition: GetLCWeights.h:176
CaloDetDescrElement::volume
float volume() const
cell volume
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:381
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:195
errorcheck.h
Helpers for checking error return status codes and reporting errors.
GetLCWeights::mapparse
void mapparse()
Definition: GetLCWeights.cxx:621
GetLCWeights::~GetLCWeights
virtual ~GetLCWeights()
Definition: GetLCWeights.cxx:156
CaloCell::ID
Identifier ID() const
get ID (from cached data member) non-virtual and inline for fast access
Definition: CaloCell.h:279
GetLCDefs::PARTICLEID_HAD
@ PARTICLEID_HAD
Definition: GetLCDefs.h:25
GetLCWeights::initialize
virtual StatusCode initialize()
Definition: GetLCWeights.cxx:161
GetLCWeights::m_weight
std::vector< std::vector< TProfile2D * > > m_weight
Vector of vector of actual histograms.
Definition: GetLCWeights.h:113
CaloSamplingHelper::getSamplingName
static const std::string & getSamplingName(const CaloSampling::CaloSample theSample)
Returns a string (name) for each CaloSampling.
Definition: CaloUtils/src/CaloSamplingHelper.cxx:42
CaloCell
Data object for each calorimeter readout cell.
Definition: CaloCell.h:57
CaloDetDescrElement::getSampling
CaloCell_ID::CaloSample getSampling() const
cell sampling
Definition: Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:395
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
xAOD::CaloCluster_v1::cell_end
const_cell_iterator cell_end() const
Definition: CaloCluster_v1.h:813
AthCommonMsg< Algorithm >::msg
MsgStream & msg() const
Definition: AthCommonMsg.h:24
GetLCDefs::LOG
@ LOG
Definition: GetLCDefs.h:21
ClusWeight
Definition: GetLCWeights.h:33
CaloClusterContainer.h
GetLCDefs::LIN
@ LIN
Definition: GetLCDefs.h:21
CaloCell_Base_ID::NSUBCALO
@ NSUBCALO
Definition: CaloCell_Base_ID.h:46
ClusWeight::next
ClusWeight * next
Definition: GetLCWeights.h:38
GetLCWeights.h
GetLCWeights::m_dimensionsmap
std::map< std::string, Gaudi::Histo1DDef > m_dimensionsmap
property to set all dimensions introduced above
Definition: GetLCWeights.h:100
ReadHandle.h
Handle class for reading from StoreGate.
IdentifierHash
Definition: IdentifierHash.h:38
CaloCell_Base_ID::NOT_VALID
@ NOT_VALID
Definition: CaloCell_Base_ID.h:46
GetLCWeights::m_ClassificationType
std::string m_ClassificationType
string to choose different classification types
Definition: GetLCWeights.h:175
CaloClusterCellLink::const_iterator
const_iterator to loop over cells belonging to a cluster
Definition: CaloClusterCellLink.h:57
GetLCWeights::m_isampmap
std::vector< std::vector< int > > m_isampmap
Vector of indices in m_dimensions for each sampling.
Definition: GetLCWeights.h:106
xAOD::CaloCluster_v1::e
virtual double e() const
The total energy of the particle.
Definition: CaloCluster_v1.cxx:265
GetLCWeights::m_clusterCollName
SG::ReadHandleKey< xAOD::CaloClusterContainer > m_clusterCollName
Name of the CaloClusterContainer to use.
Definition: GetLCWeights.h:131
CaloCell::eta
virtual double eta() const override final
get eta (through CaloDetDescrElement)
Definition: CaloCell.h:366
fitman.k
k
Definition: fitman.py:528
python.handimod.cc
int cc
Definition: handimod.py:523
Generated on Fri Jul 5 2024 21:16:28 for ATLAS Offline Software by   doxygen 1.8.18