ATLAS Offline Software
Loading...
Searching...
No Matches
GetLCWeights.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 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//###############################################################################
45GetLCWeights::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
158
159//###############################################################################
160
161StatusCode 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
338StatusCode 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
355StatusCode GetLCWeights::execute(const EventContext& ctx)
356{
357 SG::ReadHandle<xAOD::CaloClusterContainer> cc (m_clusterCollName, ctx);
358
359 std::vector<const CaloCalibrationHitContainer *> v_cchc;
360 for (const SG::ReadHandleKey<CaloCalibrationHitContainer>& k : m_CalibrationHitContainerNames) {
361 SG::ReadHandle<CaloCalibrationHitContainer> cchc (k, ctx);
362 v_cchc.push_back(cchc.cptr());
363 }
364
365 std::vector<ClusWeight *> cellVector[CaloCell_ID::NSUBCALO];
366
367 const CaloCell_ID* calo_id = nullptr;
368 ATH_CHECK(detStore()->retrieve(calo_id,"CaloCell_ID"));
369
370 for(int ic=0;ic<CaloCell_ID::NSUBCALO; ic++) {
371 unsigned int maxHashSize(0);
372 IdentifierHash myHashMin,myHashMax;
373 calo_id->calo_cell_hash_range (ic,myHashMin,myHashMax);
374 maxHashSize = myHashMax-myHashMin;
375 cellVector[ic].resize(maxHashSize,nullptr);
376 }
377
378 // loop over all cell members of all clusters and fill cell vector
379 // for used cells
380
381 // total calib hit energy of all clusters
382 double eCalibTot(0.);
383
384 unsigned int iClus = 0;
385 double nClusECalibGt0 = 0.0;
386 for (const xAOD::CaloCluster* theCluster : *cc) {
387 double eC=999;
388 if (!theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_TOT,eC)) {
389 ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_TOT");
390 return StatusCode::FAILURE;
391 }
392 if ( m_ClassificationTypeNumber != GetLCDefs::NONE ) {
393 double emFrac=-999;
394 if (!theCluster->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_FRAC_EM,emFrac)){
395 ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_FAC_EM");
396 return StatusCode::FAILURE;
397 }
398 if (m_ClassificationTypeNumber == GetLCDefs::PARTICLEID_EM && emFrac < 0.5 )
399 eC = 0;
400 if (m_ClassificationTypeNumber == GetLCDefs::PARTICLEID_HAD && emFrac > 0.5 )
401 eC = 0;
402 }
403 eCalibTot += eC;
404 if ( eC > 0 ) {
405 nClusECalibGt0++;
406 }
407 xAOD::CaloCluster::const_cell_iterator cellIter = theCluster->cell_begin();
408 xAOD::CaloCluster::const_cell_iterator cellIterEnd = theCluster->cell_end();
409 for(; cellIter != cellIterEnd; cellIter++ ){
410 const CaloCell* pCell = (*cellIter);
411 Identifier myId = pCell->ID();
412 IdentifierHash myHashId;
413 int otherSubDet(0);
414 myHashId = calo_id->subcalo_cell_hash(myId,otherSubDet);
415 ClusWeight * myClus = new ClusWeight();
416 myClus->iClus = iClus;
417 myClus->weight = cellIter.weight();
418 myClus->next = nullptr;
419 myClus->eCalibTot = 0;
420 ClusWeight * theList = cellVector[otherSubDet][(unsigned int)myHashId];
421 if ( theList ) {
422 while ( theList->next )
423 theList = theList->next;
424 theList->next = myClus;
425 }
426 else {
427 cellVector[otherSubDet][(unsigned int)myHashId] = myClus;
428 }
429 }
430 }
431
432 // calculate total calib energy of each cell in each cluster
433
434 for (const CaloCalibrationHitContainer* cchc : v_cchc) {
435 for (const CaloCalibrationHit* hit : *cchc) {
436 Identifier myId = hit->cellID();
437 int otherSubDet;
438 IdentifierHash myHashId = calo_id->subcalo_cell_hash(myId,otherSubDet);
439 if ( myHashId != CaloCell_ID::NOT_VALID ) {
440 ClusWeight * theList = cellVector[otherSubDet][(unsigned int)myHashId];
441 while ( theList ) {
442 theList->eCalibTot += hit->energyTotal();
443 theList = theList->next;
444 }
445 }
446 }
447 }
448
449 // fill weight histos for all cells in all clusters
450 // since pions can be split in several clusters the weights are
451 // calculated by importance of the cluster for the current pion - i.e.
452 // the weights get a weight itself proportinal to the calib hit energy
453 // sum of the cluster over the total calib hit energy
454
455 if ( eCalibTot > 0 && nClusECalibGt0 > 0 ) {
456 const double inv_eCalibTot = 1. / eCalibTot;
457 const double inv_nClusECalibGt0 = 1. / nClusECalibGt0;
458 for (unsigned int j=0;j<cc->size();j++) {
459 const xAOD::CaloCluster * pClus = cc->at(j);
460 double eng = pClus->e();
461 double eCalib=-999;
462 if (!pClus->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_TOT,eCalib)) {
463 ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_TOT");
464 return StatusCode::FAILURE;
465 }
466 if ( eng > 0 && eCalib > 0 ) {
467 if ( m_ClassificationTypeNumber != GetLCDefs::NONE ) {
468 double emFrac=-999;
469 if (!pClus->retrieveMoment(xAOD::CaloCluster::ENG_CALIB_FRAC_EM,emFrac)) {
470 ATH_MSG_ERROR( "Failed to retrieve cluster moment ENG_CALIB_FAC_EM");
471 return StatusCode::FAILURE;
472 }
473 if (m_ClassificationTypeNumber == GetLCDefs::PARTICLEID_EM && emFrac < 0.5 )
474 continue;
475 if (m_ClassificationTypeNumber == GetLCDefs::PARTICLEID_HAD && emFrac > 0.5 )
476 continue;
477 }
478 xAOD::CaloCluster::const_cell_iterator cellIter = pClus->cell_begin();
479 xAOD::CaloCluster::const_cell_iterator cellIterEnd = pClus->cell_end();
480 for(; cellIter != cellIterEnd; cellIter++ ){
481 const CaloCell* pCell = (*cellIter);
482 const Identifier myId = pCell->ID();
483 const CaloDetDescrElement* myCDDE=pCell->caloDDE();
484 const int caloSample = myCDDE->getSampling();//m_calo_id->calo_sample(myId);
485 if ( !m_isampmap[caloSample].empty() &&
486 !m_weight[caloSample].empty() ) {
487 int isideCell = 0;
488 int ietaCell = 0;
489 int iphiCell = 0;
490 int nside = 1;
491 int neta = 1;
492 int nphi = 1;
493 if ( m_isampmap[caloSample][1] >= 0 ) {
494 const Gaudi::Histo1DDef & hd = m_dimensions[m_isampmap[caloSample][0]][m_isampmap[caloSample][1]];
495 nside = hd.bins();
496 isideCell = (int)(nside*(((pCell->eta()<0?-1.0:1.0) - hd.lowEdge())
497 /(hd.highEdge()-hd.lowEdge())));
498 if ( isideCell < 0 || isideCell > nside-1 ) {
499 ATH_MSG_WARNING( " Side index out of bounds " <<
500 isideCell << " not in [0," << nside-1 << "] for "
501 << "Sampl=" << caloSample );
502 isideCell = -1;
503 }
504 }
505 if ( m_isampmap[caloSample][2] >= 0 ) {
506 const Gaudi::Histo1DDef & hd = m_dimensions[m_isampmap[caloSample][0]][m_isampmap[caloSample][2]];
507 neta = hd.bins();
508 ietaCell = (int)(neta*((fabs(pCell->eta()) - hd.lowEdge())
509 /(hd.highEdge()-hd.lowEdge())));
510 if ( ietaCell < 0 || ietaCell > neta-1 ) {
511 ATH_MSG_WARNING( " Eta index out of bounds " <<
512 ietaCell << " not in [0," << neta-1 << "] for "
513 << "Sampl=" << caloSample );
514 ietaCell = -1;
515 }
516 }
517 if ( m_isampmap[caloSample][3] >= 0 ) {
518 const Gaudi::Histo1DDef & hd = m_dimensions[m_isampmap[caloSample][0]][m_isampmap[caloSample][3]];
519 nphi = hd.bins();
520 iphiCell = (int)(nphi*((pCell->phi() - hd.lowEdge())
521 /(hd.highEdge()-hd.lowEdge())));
522 if ( iphiCell < 0 || iphiCell > nphi-1 ) {
523 ATH_MSG_WARNING( " Phi index out of bounds " <<
524 iphiCell << " not in [0," << nphi-1 << "] for "
525 << "Sampl=" << caloSample );
526 iphiCell = -1;
527 }
528 }
529 if ( isideCell >=0 && ietaCell >=0 && iphiCell >= 0 ) {
530 if ( myCDDE->volume() > 0 ) {
531 IdentifierHash myHashId;
532 int otherSubDet(0);
533 myHashId = calo_id->subcalo_cell_hash(myId,otherSubDet);
534 unsigned int iW = iphiCell*neta*nside+ietaCell*nside+isideCell;
535 if ( iW >= m_weight[caloSample].size() ) {
536 ATH_MSG_WARNING( " Index out of bounds " <<
537 iW << " > " << m_weight[caloSample].size()-1 << " for "
538 << "Sampl=" << caloSample
539 << ", iphi=" << iphiCell
540 << ", ieta=" << ietaCell
541 << ", iside=" << isideCell );
542 }
543 else {
544 ClusWeight * theList = cellVector[otherSubDet][(unsigned int)myHashId];
545 while ( theList && theList->iClus != j )
546 theList = theList->next;
547
548 if (m_weight[caloSample][iW] && theList && eCalibTot > 0) {
549 double norm = 0.0;
550 if ( m_NormalizationTypeNumber == GetLCDefs::LIN ) {
551 norm = eCalib * inv_eCalibTot;
552 }
553 else if ( m_NormalizationTypeNumber == GetLCDefs::LOG ) {
554 if ( eCalib > 0 ) {
555 // cluster has to have at least 1% of the calib hit E
556 norm = log10(eCalib * inv_eCalibTot)+2.0;
557 }
558 }
559 else if ( m_NormalizationTypeNumber == GetLCDefs::NCLUS ) {
560 if ( eCalib > 0 ) {
561 norm = inv_nClusECalibGt0;
562 }
563 }
564 else {
565 norm = 1.0;
566 }
567 if ( norm > 0 ) {
568 if ( !m_useInversionMethod && pCell->e()>0 ) {
569 m_weight[caloSample][iW]->Fill(log10(eng),
570 log10(pCell->e()/myCDDE->volume()),
571 theList->eCalibTot/pCell->e(),norm);
572
573 }
574 else if (m_useInversionMethod && theList->eCalibTot>0 ) {
575 m_weight[caloSample][iW]->Fill(log10(eng),
576 log10(theList->eCalibTot/myCDDE->volume()),
577 pCell->e()/theList->eCalibTot,
578 norm);
579 }
580 }
581 }
582 }
583 }
584 }
585 }
586 }
587 }
588 }
589 }
590
591 for(unsigned int ic=0;ic<CaloCell_ID::NSUBCALO; ic++) {
592 for (unsigned int ii = 0;ii<cellVector[ic].size();ii++ ) {
593 ClusWeight * theList = cellVector[ic][ii];
594 ClusWeight * prev = nullptr;
595 while ( theList) {
596 while ( theList->next ) {
597 prev = theList;
598 theList = theList->next;
599 }
600 delete theList;
601 if ( prev )
602 prev->next = nullptr;
603 else
604 cellVector[ic][ii] = nullptr;
605 theList = cellVector[ic][ii];
606 prev = nullptr;
607 }
608 }
609 }
610
611 return StatusCode::SUCCESS;
612}
613
614void GetLCWeights::mapinsert(const std::vector<Gaudi::Histo1DDef> & dims) {
615 for (unsigned int i=0;i<dims.size();i++) {
616 m_dimensionsmap[dims[0].title()+":"+dims[i].title()] = dims[i];
617 }
618}
619
620void GetLCWeights::mapparse() {
621
622 std::vector<int> theUsedSamplings(CaloSampling::Unknown,-1);
623
624 int nsamp(-1);
625
626 for (const std::pair<const std::string, Gaudi::Histo1DDef>& p : m_dimensionsmap) {
627 std::string_view dimname = std::string_view(p.first).substr(0,p.first.find(':'));
628 int theSampling(CaloSampling::Unknown);
629 for (unsigned int jsamp = 0;jsamp< CaloSampling::Unknown; jsamp++) {
630 if ( dimname == CaloSamplingHelper::getSamplingName((CaloSampling::CaloSample)jsamp)) {
631 theSampling = jsamp;
632 break;
633 }
634 }
635 if ( theSampling == CaloSampling::Unknown ) {
636 msg(MSG::ERROR) << "Calorimeter sampling " << dimname
637 << " is not a valid Calorimeter sampling name and will be ignored! "
638 << "Valid names are: ";
639 for (unsigned int jsamp = 0;jsamp< CaloSampling::Unknown; jsamp++) {
640 msg() << CaloSamplingHelper::getSamplingName((CaloSampling::CaloSample)jsamp);
641 if ( jsamp < CaloSampling::Unknown-1)
642 msg() << ", ";
643 else
644 msg() << ".";
645 }
646 msg() << endmsg;
647 }
648 else {
649 if ( theUsedSamplings[theSampling] == -1 ) {
650 nsamp++;
651 theUsedSamplings[theSampling] = nsamp;
652 m_dimensions.resize(nsamp+1);
653 m_dimensions[nsamp].resize(0);
654 }
655 m_dimensions[theUsedSamplings[theSampling]].push_back(p.second);
656 ATH_MSG_DEBUG(" New Dimension for " << dimname << ": "
657 << p.second.title() << ", [" << p.second.lowEdge()
658 << ", " << p.second.highEdge()
659 << ", " << p.second.bins()
660 << "]");
661 }
662 }
663}
664
665//###############################################################################
M_PI
#define M_PI
Definition ActiveFraction.h:11
endmsg
#define endmsg
Definition AnalysisConfig_Ntuple.cxx:63
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition AthMsgStreamMacros.h:33
ATH_MSG_FATAL
#define ATH_MSG_FATAL(x)
Definition AthMsgStreamMacros.h:34
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition AthMsgStreamMacros.h:31
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
CaloCell_ID.h
errorcheck.h
Helpers for checking error return status codes and reporting errors.
bins
static const std::vector< std::string > bins
Definition MuonTriggerSFFilesTest.cxx:22
ReadHandle.h
Handle class for reading from StoreGate.
size
size_t size() const
Number of registered mappings.
empty
static const Attributes_t empty
Definition XmlStreamer.cxx:16
AthAlgorithm::AthAlgorithm
AthAlgorithm(const std::string &name, ISvcLocator *pSvcLocator)
Constructor.
Definition AthAlgorithm.cxx:10
AthCommonAlgorithm< Gaudi::Algorithm >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
Definition AthCommonDataStore.h:145
AthCommonAlgorithm< Gaudi::Algorithm >::detStore
const ServiceHandle< StoreGateSvc > & detStore() const
Definition AthCommonDataStore.h:95
CaloCalibrationHit
Class to store calorimeter calibration hit.
Definition CaloCalibrationHit.h:24
CaloCell_Base_ID::NOT_VALID
@ NOT_VALID
Definition CaloCell_Base_ID.h:45
CaloCell_Base_ID::NSUBCALO
@ NSUBCALO
Definition CaloCell_Base_ID.h:45
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
CaloCell_Base_ID::subcalo_cell_hash
IdentifierHash subcalo_cell_hash(const Identifier cellId, int &subCalo) const
create hash id from 'global' cell id
CaloCell_ID
Helper class for offline cell identifiers.
Definition CaloCell_ID.h:34
CaloCell
Data object for each calorimeter readout cell.
Definition CaloCell.h:57
CaloCell::e
virtual double e() const override final
get energy (data member) (synonym to method energy()
Definition CaloCell.h:333
CaloCell::phi
virtual double phi() const override final
get phi (through CaloDetDescrElement)
Definition CaloCell.h:375
CaloCell::caloDDE
const CaloDetDescrElement * caloDDE() const
get pointer to CaloDetDescrElement (data member)
Definition CaloCell.h:321
CaloCell::eta
virtual double eta() const override final
get eta (through CaloDetDescrElement)
Definition CaloCell.h:382
CaloCell::ID
Identifier ID() const
get ID (from cached data member) non-virtual and inline for fast access
Definition CaloCell.h:295
CaloClusterCellLink::const_iterator::weight
weight_t weight() const
Accessor for weight associated to this cell.
Definition CaloClusterCellLink.h:82
CaloDetDescrElement
This class groups all DetDescr information related to a CaloCell.
Definition Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:66
CaloDetDescrElement::getSampling
CaloCell_ID::CaloSample getSampling() const
cell sampling
Definition Calorimeter/CaloDetDescr/CaloDetDescr/CaloDetDescrElement.h:395
CaloSamplingHelper::getSamplingName
static const std::string & getSamplingName(const CaloSampling::CaloSample theSample)
Returns a string (name) for each CaloSampling.
Definition CaloUtils/src/CaloSamplingHelper.cxx:44
ClusWeight::iClus
unsigned int iClus
Definition GetLCWeights.h:35
ClusWeight::eCalibTot
double eCalibTot
Definition GetLCWeights.h:37
ClusWeight::weight
double weight
Definition GetLCWeights.h:36
ClusWeight::next
ClusWeight * next
Definition GetLCWeights.h:38
GetLCWeights::m_NormalizationType
std::string m_NormalizationType
string to choose different normalization types
Definition GetLCWeights.h:162
GetLCWeights::m_dimensions
std::vector< std::vector< Gaudi::Histo1DDef > > m_dimensions
definition of all dimensions used for each sampling
Definition GetLCWeights.h:92
GetLCWeights::m_useInversionMethod
bool m_useInversionMethod
flag to switch on/off the use of the inversion method
Definition GetLCWeights.h:150
GetLCWeights::execute
virtual StatusCode execute(const EventContext &ctx)
Execute method.
Definition GetLCWeights.cxx:355
GetLCWeights::mapinsert
void mapinsert(const std::vector< Gaudi::Histo1DDef > &dims)
Definition GetLCWeights.cxx:614
GetLCWeights::m_outputFile
std::unique_ptr< TFile > m_outputFile
Output file to save histograms in.
Definition GetLCWeights.h:126
GetLCWeights::m_dimensionsmap
std::map< std::string, Gaudi::Histo1DDef > m_dimensionsmap
property to set all dimensions introduced above
Definition GetLCWeights.h:100
GetLCWeights::m_clusterCollName
SG::ReadHandleKey< xAOD::CaloClusterContainer > m_clusterCollName
Name of the CaloClusterContainer to use.
Definition GetLCWeights.h:131
GetLCWeights::GetLCWeights
GetLCWeights(const std::string &name, ISvcLocator *pSvcLocator)
Definition GetLCWeights.cxx:45
GetLCWeights::initialize
virtual StatusCode initialize()
Definition GetLCWeights.cxx:161
GetLCWeights::finalize
virtual StatusCode finalize()
Definition GetLCWeights.cxx:338
GetLCWeights::m_NormalizationTypeNumber
int m_NormalizationTypeNumber
Definition GetLCWeights.h:164
GetLCWeights::~GetLCWeights
virtual ~GetLCWeights()
Definition GetLCWeights.cxx:156
GetLCWeights::m_ClassificationTypeNumber
int m_ClassificationTypeNumber
Definition GetLCWeights.h:176
GetLCWeights::m_isampmap
std::vector< std::vector< int > > m_isampmap
Vector of indices in m_dimensions for each sampling.
Definition GetLCWeights.h:106
GetLCWeights::m_ClassificationType
std::string m_ClassificationType
string to choose different classification types
Definition GetLCWeights.h:175
GetLCWeights::m_outputFileName
std::string m_outputFileName
Name of the output file to save histograms in.
Definition GetLCWeights.h:120
GetLCWeights::m_weight
std::vector< std::vector< TProfile2D * > > m_weight
Vector of vector of actual histograms.
Definition GetLCWeights.h:113
GetLCWeights::m_CalibrationHitContainerNames
SG::ReadHandleKeyArray< CaloCalibrationHitContainer > m_CalibrationHitContainerNames
vector of calibration hit container names to use.
Definition GetLCWeights.h:138
IdentifierHash
This is a "hash" representation of an Identifier.
Definition IdentifierHash.h:25
SG::ReadHandleKey
Property holding a SG store/key/clid from which a ReadHandle is made.
Definition StoreGate/StoreGate/ReadHandleKey.h:40
SG::ReadHandle::cptr
const_pointer_type cptr()
Dereference the pointer.
xAOD::CaloCluster_v1::retrieveMoment
bool retrieveMoment(MomentType type, double &value) const
Retrieve individual moment.
Definition CaloCluster_v1.cxx:662
xAOD::CaloCluster_v1::e
virtual double e() const
The total energy of the particle.
Definition CaloCluster_v1.cxx:265
xAOD::CaloCluster_v1::const_cell_iterator
CaloClusterCellLink::const_iterator const_cell_iterator
Iterator of the underlying CaloClusterCellLink (explicitly const version).
Definition CaloCluster_v1.h:796
xAOD::CaloCluster_v1::cell_end
const_cell_iterator cell_end() const
Definition CaloCluster_v1.h:816
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:251
xAOD::CaloCluster_v1::ENG_CALIB_TOT
@ ENG_CALIB_TOT
Calibration Hit energy inside the cluster.
Definition CaloCluster_v1.h:198
xAOD::CaloCluster_v1::cell_begin
const_cell_iterator cell_begin() const
Iterator of the underlying CaloClusterCellLink (const version).
Definition CaloCluster_v1.h:815
GetLCDefs::PARTICLEID_EM
@ PARTICLEID_EM
Definition GetLCDefs.h:25
GetLCDefs::PARTICLEID_HAD
@ PARTICLEID_HAD
Definition GetLCDefs.h:25
xAOD::CaloCluster
CaloCluster_v1 CaloCluster
Define the latest version of the calorimeter cluster class.
Definition Event/xAOD/xAODCaloEvent/xAODCaloEvent/CaloCluster.h:19
Generated on for ATLAS Offline Software by doxygen 1.16.1