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()
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 && nClusECalibGt0 > 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
615void 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
621void 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//###############################################################################
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.
empty
static const Attributes_t empty
Definition XmlStreamer.cxx:16
AthAlgorithm::AthAlgorithm
AthAlgorithm(const std::string &name, ISvcLocator *pSvcLocator)
Constructor with parameters:
Definition AthAlgorithm.cxx:24
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
Definition AthCommonDataStore.h:145
AthCommonDataStore< AthCommonMsg< Algorithm > >::detStore
const ServiceHandle< StoreGateSvc > & detStore() const
Definition AthCommonDataStore.h:95
AthCommonMsg< Algorithm >::msg
MsgStream & msg() const
Definition AthCommonMsg.h:24
CaloCalibrationHit
Class to store calorimeter calibration hit.
Definition CaloCalibrationHit.h:23
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::mapinsert
void mapinsert(const std::vector< Gaudi::Histo1DDef > &dims)
Definition GetLCWeights.cxx:615
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::execute
virtual StatusCode execute()
Definition GetLCWeights.cxx:355
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.14.0