ATLAS Offline Software
Loading...
Searching...
No Matches
CaloCalibClusterMomentsMaker2.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: CaloCalibClusterMomentsMaker2.cxx,v 1.16 2009-05-18 16:16:49 pospelov Exp $
8//
9// Description: see CaloCalibClusterMomentsMaker2.h
10//
11// Environment:
12// Software developed for the ATLAS Detector at CERN LHC
13//
14// Author List:
15// Sven Menke
16//
17//-----------------------------------------------------------------------
18
19//-----------------------
20// This Class's Header --
21//-----------------------
23
24//---------------
25// C++ Headers --
26//---------------
27#include <iterator>
28#include <sstream>
29#include <set>
30
31#include "CaloEvent/CaloCell.h"
37
39
41
42#include "CLHEP/Units/SystemOfUnits.h"
43
44#include <CLHEP/Vector/LorentzVector.h>
45#include <cmath>
46
47
48using CLHEP::HepLorentzVector;
49using CLHEP::MeV;
50using CLHEP::cm;
51
52
53//###############################################################################
54
56 const std::string& name,
57 const IInterface* parent)
58 : AthAlgTool(type, name, parent),
59 m_calo_id(nullptr),
60 m_caloDM_ID(nullptr),
61 m_caloDmDescrManager(nullptr),
62 m_useParticleID(true),
63 m_energyMin(200*MeV),
64 m_energyMinCalib(20*MeV),
65 m_apars_alpha(0.5),
66 m_apars_r0(0.2),
68{
69 declareInterface<CaloClusterCollectionProcessor> (this);
70 // Name(s) of Moments to calculate
71 declareProperty("MomentsNames",m_momentsNames);
72 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_TOT"),xAOD::CaloCluster::ENG_CALIB_TOT));
73 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_OUT_L"),xAOD::CaloCluster::ENG_CALIB_OUT_L));
74 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_OUT_M"),xAOD::CaloCluster::ENG_CALIB_OUT_M));
75 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_OUT_T"),xAOD::CaloCluster::ENG_CALIB_OUT_T));
76 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_EMB0"),xAOD::CaloCluster::ENG_CALIB_EMB0));
77 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_EME0"),xAOD::CaloCluster::ENG_CALIB_EME0));
78 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_TILEG3"),xAOD::CaloCluster::ENG_CALIB_TILEG3));
79 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_TOT"),xAOD::CaloCluster::ENG_CALIB_DEAD_TOT));
80 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_EMB0"),xAOD::CaloCluster::ENG_CALIB_DEAD_EMB0));
81 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_TILE0"),xAOD::CaloCluster::ENG_CALIB_DEAD_TILE0));
82 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_TILEG3"),xAOD::CaloCluster::ENG_CALIB_DEAD_TILEG3));
83 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_EME0"),xAOD::CaloCluster::ENG_CALIB_DEAD_EME0));
84 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_HEC0"),xAOD::CaloCluster::ENG_CALIB_DEAD_HEC0));
85 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_FCAL"),xAOD::CaloCluster::ENG_CALIB_DEAD_FCAL));
86 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_LEAKAGE"),xAOD::CaloCluster::ENG_CALIB_DEAD_LEAKAGE));
87 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_DEAD_UNCLASS"),xAOD::CaloCluster::ENG_CALIB_DEAD_UNCLASS));
88 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_FRAC_EM"),xAOD::CaloCluster::ENG_CALIB_FRAC_EM));
89 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_FRAC_HAD"),xAOD::CaloCluster::ENG_CALIB_FRAC_HAD));
90 m_validNames.push_back(moment_name_pair(std::string("ENG_CALIB_FRAC_REST"),xAOD::CaloCluster::ENG_CALIB_FRAC_REST));
91
92 // Name(s) of Moments which can be stored on the AOD - all others go to ESD
93 m_momentsNamesAOD.emplace_back("ENG_CALIB_TOT");
94 m_momentsNamesAOD.emplace_back("ENG_CALIB_OUT_L");
95 m_momentsNamesAOD.emplace_back("ENG_CALIB_OUT_M");
96 m_momentsNamesAOD.emplace_back("ENG_CALIB_OUT_T");
97 m_momentsNamesAOD.emplace_back("ENG_CALIB_EMB0");
98 m_momentsNamesAOD.emplace_back("ENG_CALIB_EME0");
99 m_momentsNamesAOD.emplace_back("ENG_CALIB_TILEG3");
100 m_momentsNamesAOD.emplace_back("ENG_CALIB_DEAD_TOT");
101 m_momentsNamesAOD.emplace_back("ENG_CALIB_DEAD_EMB0");
102 m_momentsNamesAOD.emplace_back("ENG_CALIB_DEAD_TILE0");
103 m_momentsNamesAOD.emplace_back("ENG_CALIB_DEAD_TILEG3");
104 m_momentsNamesAOD.emplace_back("ENG_CALIB_DEAD_EME0");
105 m_momentsNamesAOD.emplace_back("ENG_CALIB_DEAD_HEC0");
106 m_momentsNamesAOD.emplace_back("ENG_CALIB_DEAD_FCAL");
107 m_momentsNamesAOD.emplace_back("ENG_CALIB_DEAD_LEAKAGE");
108 m_momentsNamesAOD.emplace_back("ENG_CALIB_DEAD_UNCLASS");
109 m_momentsNamesAOD.emplace_back("ENG_CALIB_FRAC_EM");
110 m_momentsNamesAOD.emplace_back("ENG_CALIB_FRAC_HAD");
111 m_momentsNamesAOD.emplace_back("ENG_CALIB_FRAC_REST");
112
113 declareProperty("AODMomentsNames",m_momentsNamesAOD);
114 declareProperty("CalibrationHitContainerNames",m_CalibrationHitContainerNames);
115 declareProperty("DMCalibrationHitContainerNames",m_DMCalibrationHitContainerNames);
116 m_n_phi_out = 127; // not more than 127 since we store indices (-127,...-1,0,...,126) and have 8 bits only
117 m_n_eta_out = 127;
119 m_out_eta_max = 6;
120
121 m_rmaxOut[0] = 1.0;
122 m_rmaxOut[1] = 0.5;
123 m_rmaxOut[2] = 0.3;
124
125 for( int im=0;im<3;im++) {
126 m_i_phi_eta[im].resize(m_n_eta_out);
127 }
128 m_doDeadEnergySharing = false;
129 m_foundAllContainers = false;
130 m_doOutOfClusterL = false;
131 m_doOutOfClusterM = false;
132 m_doOutOfClusterT = false;
133 m_doDeadL = false;
134 m_doDeadM = false;
135 m_doDeadT = false;
136 m_doCalibFrac = false;
137
138 declareProperty("MatchDmType",m_MatchDmType);
139
140 declareProperty( "UseParticleID",m_useParticleID);
141}
142
143
144
145//###############################################################################
146
148{
149 ATH_MSG_INFO( "Initializing " << name() );
150
151 for (const std::string& name : m_momentsNames) {
152 bool isValid(false);
153 for (const moment_name_pair& vname : m_validNames) {
154 if ( name == vname.first ) {
155 m_validMoments.insert(vname);
156 isValid = true;
157 ATH_MSG_DEBUG( "Inserting " << name );
158 break;
159 }
160 }
161 if ( !isValid) {
162 msg() << MSG::ERROR << "Moment " << name
163 << " is not a valid Moment name and will be ignored! "
164 << "Valid names are:";
165 for (unsigned int i=0;i<m_validNames.size();i++)
166 msg() << (i==0?" ":", ") << m_validNames[i].first;
167 msg() << endmsg;
168 }
169 }
170
171 // to switch on clever DeadMaterial assignment procedure
172 // and check if tight, medium and/or loose versions for out-of-cluster and
173 // simple dead-material assignment are wanted
174 for (const moment_name_pair& vname : m_validNames) {
175 switch (vname.second) {
178 break;
180 m_doOutOfClusterL = true;
181 break;
183 m_doOutOfClusterM = true;
184 break;
186 m_doOutOfClusterT = true;
187 break;
189 m_doCalibFrac = true;
190 break;
192 m_doCalibFrac = true;
193 break;
195 m_doCalibFrac = true;
196 break;
197 default:
198 break;
199 }
200 }
201
203 ATH_MSG_INFO( "Usage of ParticleID was switched off (UseParticleID==False), no ENG_CALIB_FRAC_* moments will be available" );
204 m_doCalibFrac = false;
205 }
206
207 for (const std::string& name : m_momentsNamesAOD) {
208 for (const moment_name_pair& vname : m_validNames) {
209 if ( vname.first == name ) {
210 m_momentsAOD.insert(vname.second);
211 break;
212 }
213 }
214 }
217 }
218
219 // dead material identifier description manager
221
222 ATH_CHECK( detStore()->retrieve(m_calo_id, "CaloCell_ID") );
223 ATH_CHECK( detStore()->retrieve(m_caloDM_ID) );
224
225 // initialize distance tables
226 for(int jeta = 0;jeta<m_n_eta_out;jeta++) {
227 double eta0 = (jeta+0.5) * (m_out_eta_max)/m_n_eta_out;
228 HepLorentzVector middle(1,0,0,1);
229 middle.setREtaPhi(1./cosh(eta0),eta0,0);
230 double x_rmaxOut[3];
231 for (int im=0;im<3;im++) {
232 x_rmaxOut[im] = m_rmaxOut[im]*angle_mollier_factor(eta0);
233 }
234 for (int jp=-m_n_phi_out;jp<m_n_phi_out;jp++) {
235 double phi = (jp+0.5) * m_out_phi_max/m_n_phi_out;
236 int ietaMin[3] = {m_n_eta_out,m_n_eta_out,m_n_eta_out};
237 int ietaMax[3] = {-m_n_eta_out,-m_n_eta_out,-m_n_eta_out};
238 for(int je = -m_n_eta_out;je<m_n_eta_out;je++) {
239 double eta = (je+0.5) * m_out_eta_max/m_n_eta_out;
240 HepLorentzVector cpoint(1,0,0,1);
241 cpoint.setREtaPhi(1./cosh(eta),eta,phi);
242 double r = middle.angle(cpoint.vect());
243 for (int im=0;im<3;im++) {
244 if ( r < x_rmaxOut[im] ) {
245 if ( je < ietaMin[im] )
246 ietaMin[im] = je;
247 if ( je > ietaMax[im] )
248 ietaMax[im] = je;
249 }
250 }
251 }
252 for (int im=0;im<3;im++) {
253 if ( ietaMin[im] <= ietaMax[im] ) {
254 CalibHitIPhiIEtaRange theRange{};
255 theRange.iPhi = (char)jp;
256 theRange.iEtaMin = (char)ietaMin[im];
257 theRange.iEtaMax = (char)ietaMax[im];
258 m_i_phi_eta[im][jeta].push_back(theRange);
259 }
260 }
261 }
262 }
263
266
268
269 ATH_CHECK(m_caloDetDescrMgrKey.initialize());
270
271 return StatusCode::SUCCESS;
272}
273
274//###############################################################################
275
276StatusCode
278 xAOD::CaloClusterContainer *theClusColl) const
279{
280
281 ATH_MSG_DEBUG("Starting CaloCalibClusterMomentsMaker2::execute");
283 const CaloDetDescrManager* calo_dd_man = *caloMgrHandle;
284
285 bool foundAllContainers (true);
286 std::vector<const CaloCalibrationHitContainer *> v_cchc;
289 {
291 if ( !cchc.isValid() ) {
293 // print ERROR message only if there was at least one event with
294 // all containers
295 msg(MSG::ERROR) << "SG does not contain calibration hit container " << key.key() << endmsg;
296 }
297 foundAllContainers = false;
298 }
299 else {
300 v_cchc.push_back(cchc.cptr());
301 }
302 }
303
304 std::vector<const CaloCalibrationHitContainer *> v_dmcchc;
307 {
309 if ( !cchc.isValid() ) {
311 // print ERROR message only if there was at least one event with
312 // all containers
313 ATH_MSG_ERROR("SG does not contain DM calibration hit container " << key.key());
314 }
315 foundAllContainers = false;
316 }
317 else {
318 v_dmcchc.push_back(cchc.cptr());
319 }
320 }
321
322 if ( !m_foundAllContainers && foundAllContainers )
324
325 if ( !foundAllContainers ) return StatusCode::SUCCESS;
326 ATH_MSG_DEBUG("SG has all containers ");
327
328 // will contain detailed info about cluster calibration eneries
329 ClusInfo_t clusInfoVec (theClusColl->size());
330
331 CellInfoSet_t cellInfo;
332
333 /* ********************************************
334 filling the map with info which cell belongs to which cluster and what
335 is the cell weight in the cluster
336 ******************************************** */
337 xAOD::CaloClusterContainer::iterator clusIter = theClusColl->begin();
338 xAOD::CaloClusterContainer::iterator clusIterEnd = theClusColl->end();
339 int iClus = 0;
340 for( ;clusIter!=clusIterEnd;++clusIter,++iClus) {
341 const xAOD::CaloCluster * theCluster = (*clusIter);
342
343 // loop over all cell members and fill cell vector for used cells
344 xAOD::CaloCluster::const_cell_iterator cellIter = theCluster->cell_begin();
345 xAOD::CaloCluster::const_cell_iterator cellIterEnd = theCluster->cell_end();
346 for(; cellIter != cellIterEnd; cellIter++ ){
347 const CaloCell* pCell = (*cellIter);
348 Identifier myId = pCell->ID();
349
350 MyCellInfo info (iClus, cellIter.weight() );
351 CellInfoSet_t::iterator bookmark = cellInfo.lower_bound( myId );
352 if(bookmark == cellInfo.end() || bookmark->first != myId) {
353 // We haven't had a infohit in this cell before. Add it to our set.
354 if (bookmark != cellInfo.begin()) --bookmark;
355 cellInfo.emplace_hint (bookmark, myId, std::move(info));
356 }else{
357 // Update the existing hit.
358 bookmark->second.Add( info );
359 }
360 } // cellIter
361 } // iClus
362
363
364 /* ********************************************
365 calculate total calib energy inside clusters
366 ******************************************** */
367 unsigned int nHitsTotal = 0;
368 unsigned int nHitsWithoutParticleUID = 0;
369 for (const CaloCalibrationHitContainer* cchc : v_cchc) {
370 //loop over cells in calibration container
371 for (const CaloCalibrationHit* hit : *cchc) {
372 Identifier myId = hit->cellID();
373
374 CellInfoSet_t::iterator pos = cellInfo.find( myId );
375 if(pos != cellInfo.end() ) {
376 // i.e. given hit id belongs to one or more clusters
378 for ( const std::pair<int, double>& p : pos->second) {
379 int iClus = p.first;
380 double weight = p.second;
381 if(m_useParticleID){
382 const int uniqueID = HepMC::uniqueID(hit);
383 if (uniqueID == HepMC::INVALID_PARTICLE_ID) {
384 ATH_MSG_ERROR("Invalid uniqueID detected - this sample cannot be properly analysed.");
385 break;
386 }
387 clusInfoVec[iClus].Add(weight * hit->energyTotal(), nsmp, (unsigned int)(HepMC::uniqueID(hit)));
388 }else{
389 clusInfoVec[iClus].Add(weight * hit->energyTotal(), nsmp);
390 }
391 }
392 }
394 nHitsWithoutParticleUID++;
395 }
396 nHitsTotal++;
397 }
398 }
399
400 // if all calibration hits have ParticleUID(i.e GenParticle::id())==0 when simulation was done without ParticleUID
401 bool doCalibFrac = m_doCalibFrac;
402 bool useParticleID = m_useParticleID;
403 if(m_useParticleID && (nHitsTotal == nHitsWithoutParticleUID) ) {
404 ATH_MSG_INFO("Calibration hits do not have ParticleUID, ids of particle-caused hits are always 0. Continuing without ParticleID machinery.");
405 useParticleID = false;
406 }
407
408 // reading particle information for later calcution of calibration enegry fraction caused
409 // by particles of different types
411
412 if (doCalibFrac && !truthParticleContainerReadHandle.isValid()){
413 ATH_MSG_WARNING("Invalid read handle to TruthParticleContainer with key: " << m_truthParticleContainerKey.key());
414 doCalibFrac = false;
415 }
416
417 std::vector<double> engCalibOut[3];
418
419 /* ****************************************************************
420 lists of clusters populating given area [eta*phi][iClus list]
421 **************************************************************** */
422 std::vector<std::vector <int > > clusListL;
423 std::vector<std::vector <int > > clusListM;
424 std::vector<std::vector <int > > clusListT;
426 clusListL.resize((2*m_n_phi_out+1)*(2*m_n_eta_out+1));
428 clusListM.resize((2*m_n_phi_out+1)*(2*m_n_eta_out+1));
430 clusListT.resize((2*m_n_phi_out+1)*(2*m_n_eta_out+1));
431
432 for(unsigned int ii=0;ii<3;ii++) {
436 engCalibOut[ii].resize(theClusColl->size(),0);
437 iClus = -1;
438 for (const xAOD::CaloCluster * theCluster : *theClusColl) {
439 ++iClus;
440 MyClusInfo& clusInfo = clusInfoVec[iClus];
441
442 if ( clusInfo.engCalibIn.engTot > 0 ) {
443 int iEtaSign = 1;
444 if ( theCluster->eta() < 0 ) iEtaSign = -1;
445 int jeta = (int)(floor(m_n_eta_out*(theCluster->eta()/(m_out_eta_max))));
446 int jphi = (int)(floor(m_n_phi_out*(theCluster->phi())/(m_out_phi_max)));
447 if ( jeta >= -m_n_eta_out && jeta < m_n_eta_out ) {
448 if ( jphi < -m_n_phi_out ) jphi += 2*m_n_phi_out;
449 if ( jphi >= m_n_phi_out ) jphi -= 2*m_n_phi_out;
450 unsigned int iEtaBin(jeta);
451 if ( jeta < 0 ) iEtaBin = abs(jeta)-1;
452 const std::vector<CalibHitIPhiIEtaRange>& bins
453 = m_i_phi_eta[ii][iEtaBin];
454 for (const CalibHitIPhiIEtaRange& range : bins) {
455 int jp = range.iPhi+jphi;
456 if ( jp < -m_n_phi_out ) jp += 2*m_n_phi_out;
457 if ( jp >= m_n_phi_out ) jp -= 2*m_n_phi_out;
458 int jEtaMin = iEtaSign<0?-range.iEtaMax-1:range.iEtaMin;
459 int jEtaMax = iEtaSign<0?-range.iEtaMin-1:range.iEtaMax;
460 for( int je = jEtaMin;je<=jEtaMax;je++ ) {
461 if(ii == 0 ) clusListL[(jp+m_n_phi_out)*(2*m_n_eta_out+1)+je+m_n_eta_out].push_back(iClus);
462 else if(ii == 1 ) clusListM[(jp+m_n_phi_out)*(2*m_n_eta_out+1)+je+m_n_eta_out].push_back(iClus);
463 else if(ii == 2 ) clusListT[(jp+m_n_phi_out)*(2*m_n_eta_out+1)+je+m_n_eta_out].push_back(iClus);
464 }
465 }
466 }
467 }
468 }
469 }
470 }
471
472 /* ****************************************************************
473 calculate out-of-cluster energy of clusters
474 **************************************************************** */
476 for (const CaloCalibrationHitContainer* cchc : v_cchc) {
477 //loop over cells in calibration container
478 for (const CaloCalibrationHit* hit : *cchc) {
479 Identifier myId = hit->cellID();
480
481 CellInfoSet_t::iterator pos = cellInfo.find( myId );
482 if(pos == cellInfo.end() ) {
483 // hit is not inside any cluster
484 const CaloDetDescrElement* myCDDE =
485 calo_dd_man->get_element(myId);
486 int uniqueID(HepMC::UNDEFINED_ID);
487 if(useParticleID) uniqueID = HepMC::uniqueID(hit);
488 if ( myCDDE ) {
489 int jeO = (int)floor(m_n_eta_out*(myCDDE->eta()/m_out_eta_max));
490 if ( jeO >= -m_n_eta_out && jeO < m_n_eta_out ) {
491 int jpO = (int)floor(m_n_phi_out*(myCDDE->phi()/m_out_phi_max));
492 if ( jpO < -m_n_phi_out ) jpO += 2*m_n_phi_out;
493 if ( jpO >= m_n_phi_out ) jpO -= 2*m_n_phi_out;
494 for (unsigned int ii=0;ii<3;ii++) {
495 std::vector<std::vector <int > > *pClusList=nullptr;
496 if ( ii == 0 && m_doOutOfClusterL )
497 pClusList = &clusListL;
498 else if ( ii == 1 && m_doOutOfClusterM )
499 pClusList = &clusListM;
500 else if ( ii == 2 && m_doOutOfClusterT )
501 pClusList = &clusListT;
502 if ( pClusList) {
503 // loop over list of potential neighboring clusters
504 double hitClusNorm(0.);
505 std::vector<int> hitClusIndex;
506 std::vector<double > hitClusEffEnergy;
507 // loop over clusters which match given OOC hit
508 for(unsigned int i_cls=0; i_cls<(*pClusList)[(jpO+m_n_phi_out)*(2*m_n_eta_out+1)+jeO+m_n_eta_out].size(); i_cls++){
509 int iClus = (*pClusList)[(jpO+m_n_phi_out)*(2*m_n_eta_out+1)+jeO+m_n_eta_out][i_cls];
510 MyClusInfo& clusInfo = clusInfoVec[iClus];
511 // getting access to calibration energy inside cluster caused by same particleUID (uniqueID)
512 // as given OOC hit
513 auto pos = clusInfo.engCalibParticle.find(uniqueID);
514 if(pos!=clusInfo.engCalibParticle.end()) {
515 // given cluster have some energy inside caused by same particle as given OOC hitClusEffEnergy
516 // so the hit will be assigned to this cluster with some weight
517 hitClusNorm += pos->second.engTot;
518 hitClusEffEnergy.push_back(pos->second.engTot);
519 hitClusIndex.push_back(iClus);
520 }
521 }
522 if ( hitClusNorm > 0 ) {
523 const double inv_hitClusNorm = 1. / hitClusNorm;
524 for(unsigned int i_cls=0; i_cls<hitClusIndex.size(); i_cls++){
525 int iClus = hitClusIndex[i_cls];
526 double w = hitClusEffEnergy[i_cls] * inv_hitClusNorm;
527 engCalibOut[ii][iClus] += w*hit->energyTotal();
528 }
529 }
530 }
531 }
532 }
533 }
534 } // hits not in cluster
535 }
536 }
537 }
538
539 // ------------------------------------------------------------------------
540 // calculate dead-material energy of clusters (way2)
541 // + energy is shared among clusters within certain area
542 // + distance to clusters and energy in specific samplings are used as sharing criteria
543 // + calculations are done separately for different dead material areas
544 std::vector<std::vector <int > > *pClusList=nullptr;
545 if ( m_MatchDmType == kMatchDmLoose ) {
546 pClusList = &clusListL;
547 } else if ( m_MatchDmType == kMatchDmMedium ) {
548 pClusList = &clusListM;
549 } else if ( m_MatchDmType == kMatchDmTight ) {
550 pClusList = &clusListT;
551 }
552 if( m_doDeadEnergySharing && pClusList) {
553
554 for (const CaloCalibrationHitContainer* dmcchc : v_dmcchc) {
555 for (const CaloCalibrationHit* hit : *dmcchc) {
556 Identifier myId = hit->cellID();
557 if (m_calo_id->is_lar_dm(myId) || m_calo_id->is_tile_dm(myId)) {
558 CaloDmDescrElement* myCDDE(nullptr);
559 myCDDE = m_caloDmDescrManager->get_element(myId);
560 if ( myCDDE ) {
561 int uniqueID(HepMC::UNDEFINED_ID);
562 if(useParticleID) uniqueID = HepMC::uniqueID(hit);
563
564 int jeO = (int)floor(m_n_eta_out*(myCDDE->eta()/m_out_eta_max));
565 if ( jeO >= -m_n_eta_out && jeO < m_n_eta_out ) {
566 int jpO = (int)floor(m_n_phi_out*(myCDDE->phi()/m_out_phi_max));
567 if ( jpO < -m_n_phi_out ) jpO += 2*m_n_phi_out;
568 if ( jpO >= m_n_phi_out ) jpO -= 2*m_n_phi_out;
569
570 /* ****************************************
571 share energy of DM hit among different clusters
572 ***************************************** */
573 int nDmArea = m_caloDmDescrManager->get_dm_area(myId);
574 const CaloDmRegion *dmRegion = m_caloDmDescrManager->get_dm_region(myId);
575 int hitClusVecIndex = 0;
576 std::vector<int> hitClusIndex(theClusColl->size());
577 std::vector<double > hitClusEffEnergy(theClusColl->size());
578 double hitClusNorm = 0.0;
579 // loop over clusters which match given DM hit
580 for(unsigned int i_cls=0; i_cls< (*pClusList)[(jpO+m_n_phi_out)*(2*m_n_eta_out+1)+jeO+m_n_eta_out].size(); i_cls++){
581 int iClus = (*pClusList)[(jpO+m_n_phi_out)*(2*m_n_eta_out+1)+jeO+m_n_eta_out][i_cls];
582 xAOD::CaloCluster * theCluster = theClusColl->at(iClus);
583
584 MyClusInfo& clusInfo = clusInfoVec[iClus];
585 // getting access to calibration energy inside cluster caused by same particleUID (uniqueID)
586 // as given OOC hit
587 auto pos = clusInfo.engCalibParticle.find(uniqueID);
588 if(pos!=clusInfo.engCalibParticle.end()) {
589 double engClusTruthUniqueIDCalib = pos->second.engTot;
590
591 if(engClusTruthUniqueIDCalib > m_energyMinCalib && theCluster->e()>m_energyMin) {
592 double sum_smp_energy = 0.0;
593 // loop over calo sampling numbers registered for given DM area
594 for(unsigned int i_smp=0; i_smp<dmRegion->m_CaloSampleNeighbours.size(); i_smp++) {
596 if( (dmRegion->m_CaloSampleEtaMin[i_smp]-0.5) <= theCluster->eta() &&
597 theCluster->eta() <= (dmRegion->m_CaloSampleEtaMax[i_smp]+0.5) ){
598 //sum_smp_energy += theCluster->eSample(nsmp);
599 sum_smp_energy += pos->second.engSmp[nsmp];
600 }
601 }
602 if(sum_smp_energy > 0.0) {
603 double phi_diff=myCDDE->phi()-theCluster->phi();
604 if(phi_diff <= -M_PI){
605 phi_diff += 2.*M_PI;
606 } else if (phi_diff > M_PI){
607 phi_diff -= 2.*M_PI;
608 }
609 double eta_diff = (myCDDE->eta()-theCluster->eta());
610 float distance=sqrt(eta_diff * eta_diff + phi_diff * phi_diff);
611
612 double effEner = pow(sum_smp_energy,m_apars_alpha)*exp(-distance/m_apars_r0);
613
614 hitClusIndex [hitClusVecIndex] = iClus;
615 hitClusEffEnergy [hitClusVecIndex++]= effEner;
616 hitClusNorm += effEner;
617 }
618 } // good energetic cluster
619 } // we ve found a uniqueID
620 } // loop over clusters in the list
621 hitClusIndex.resize(hitClusVecIndex);
622 hitClusEffEnergy.resize(hitClusVecIndex);
623
624
625 // now we have to calculate weight for assignment hit energy to cluster
626 if(hitClusNorm >0.0) {
627 const double inv_hitClusNorm = 1. / hitClusNorm;
628 for(unsigned int i_cls=0; i_cls<hitClusIndex.size(); i_cls++){
629 int iClus = hitClusIndex[i_cls];
630 double dm_weight = hitClusEffEnergy[i_cls] * inv_hitClusNorm;
631 if(dm_weight > 1.0 || dm_weight < 0.0 ){
632 std::cout << "CaloCalibClusterMomentsMaker2::execute() ->Error! Strange weight " << dm_weight<< std::endl;
633 std::cout << hitClusEffEnergy[i_cls] << " " << hitClusNorm << std::endl;
634 }
635 clusInfoVec[iClus].engCalibDeadInArea[nDmArea] += hit->energyTotal()*dm_weight;
636 clusInfoVec[iClus].engCalibDeadInArea[CaloDmDescrArea::DMA_ALL] += hit->energyTotal()*dm_weight;
637 }
638 } // hit clus norm
639
640 } // loop over united ieta*iphi
641 } // myCDDE
642 } // tile and lar
643 } // eol over DM calibration hits
644 } //eol over containers
645
646 } // doDeadEnergySharing
647
648 // fraction of calibration energies caused by different particles
649 std::vector<double> engCalibFrac;
650 engCalibFrac.resize(kCalibFracMax, 0.0);
651
652 std::map<unsigned int,int> truthIDToPdgCodeMap;
653
654 //loop on truth particle container is slow, so put needed information in a map for faster key lookup in later loops
655 for ( const auto *thisTruthParticle : *truthParticleContainerReadHandle){
656
657 if (!thisTruthParticle){
658 ATH_MSG_WARNING("Got invalid pointer to TruthParticle");
659 continue;
660 }
661
662 truthIDToPdgCodeMap[HepMC::uniqueID(thisTruthParticle)] = thisTruthParticle->pdgId();
663 }//truth particle loop
664
665 // assign moments
666 for( clusIter = theClusColl->begin(),iClus=0; clusIter!=clusIterEnd;++clusIter,++iClus) {
667 xAOD::CaloCluster * theCluster = *clusIter;
668 MyClusInfo& clusInfo = clusInfoVec[iClus];
669
670 // total DM energy assigned to cluster
671 double eng_calib_dead_tot = clusInfo.engCalibDeadInArea[CaloDmDescrArea::DMA_ALL]
672 + clusInfo.engCalibIn.engSmp[CaloSampling::PreSamplerB]
673 + clusInfo.engCalibIn.engSmp[CaloSampling::PreSamplerE]
674 + clusInfo.engCalibIn.engSmp[CaloSampling::TileGap3];
675 // DM energy before barrel presampler, inside it, and between presampler and strips
676 double eng_calib_dead_emb0 = clusInfo.engCalibDeadInArea[CaloDmDescrArea::DMA_EMB0]
678 + clusInfo.engCalibIn.engSmp[CaloSampling::PreSamplerB];
679 // DM energy between barrel and tile
680 double eng_calib_dead_tile0 = clusInfo.engCalibDeadInArea[CaloDmDescrArea::DMA_EMB3_TILE0];
681 // DM energy before scintillator and inside scintillator
682 double eng_calib_dead_tileg3 = clusInfo.engCalibDeadInArea[CaloDmDescrArea::DMA_SCN]
683 + clusInfo.engCalibIn.engSmp[CaloSampling::TileGap3];
684 // DM energy beforee endcap presampler, inside it and between presampler and strips
685 double eng_calib_dead_eme0 = clusInfo.engCalibDeadInArea[CaloDmDescrArea::DMA_EME0]
687 + clusInfo.engCalibIn.engSmp[CaloSampling::PreSamplerE];
688 // DM energy between emec and hec
689 double eng_calib_dead_hec0 = clusInfo.engCalibDeadInArea[CaloDmDescrArea::DMA_EME3_HEC0];
690 // DM energy before FCAL and between HEC and FCAL
691 double eng_calib_dead_fcal = clusInfo.engCalibDeadInArea[CaloDmDescrArea::DMA_FCAL0]
693 // DM leakage behind the calorimeter
694 double eng_calib_dead_leakage = clusInfo.engCalibDeadInArea[CaloDmDescrArea::DMA_LEAK];
695 // the rest of DM energy which remains unclassified
696 double eng_calib_dead_unclass = eng_calib_dead_tot - eng_calib_dead_emb0 - eng_calib_dead_tile0
697 - eng_calib_dead_tileg3 - eng_calib_dead_eme0 - eng_calib_dead_hec0 - eng_calib_dead_fcal
698 - eng_calib_dead_leakage;
699
700 if(doCalibFrac){
701/*****************************************************************************
702Calculation of energy fraction caused by particles of different types
703*****************************************************************************/
704 engCalibFrac.assign(kCalibFracMax, 0.0);
705 if(clusInfo.engCalibIn.engTot > 0.0) {
706 // each MyClusInfo has a map of particle's uniqueID (GenParticle::id()) and particle calibration deposits in given cluster
707 for (const auto& p : clusInfo.engCalibParticle) {
708 int pdg_id = 0;
709 if ( auto it = truthIDToPdgCodeMap.find(p.first); it != truthIDToPdgCodeMap.end()) {
710 pdg_id = it->second;
711 } else {
712 ATH_MSG_WARNING("truthIDToPdgCodeMap cannot find an entry with uniqueID " << p.first);
713 continue;
714 }
715 if( std::abs(pdg_id) == 211) {
716 engCalibFrac[kCalibFracHAD] += p.second.engTot;
717 } else if( pdg_id == 111 || pdg_id == 22 || std::abs(pdg_id)==11) {
718 engCalibFrac[kCalibFracEM] += p.second.engTot;
719 } else {
720 engCalibFrac[kCalibFracREST] += p.second.engTot;
721 }
722 }
723 for(size_t i=0; i<engCalibFrac.size(); i++) engCalibFrac[i] = engCalibFrac[i]/clusInfo.engCalibIn.engTot;
724 }
725 }
726
727 if ( !m_momentsNames.empty() ) {
728 std::vector<double> myMoments(m_validMoments.size(),0);
729 // assign moments
730 moment_name_set::const_iterator vMomentsIter = m_validMoments.begin();
731 moment_name_set::const_iterator vMomentsIterEnd = m_validMoments.end();
732
733 int iMoment=0;
734 for(; vMomentsIter!=vMomentsIterEnd; ++vMomentsIter,++iMoment) {
735 // now calculate the actual moments
736 switch (vMomentsIter->second) {
738 ATH_MSG_DEBUG("Inserting ENG_CALIB_TOT");
739 myMoments[iMoment] = clusInfo.engCalibIn.engTot;
740 break;
742 myMoments[iMoment] = engCalibOut[0][iClus];
743 break;
745 myMoments[iMoment] = engCalibOut[1][iClus];
746 break;
748 myMoments[iMoment] = engCalibOut[2][iClus];
749 break;
751 myMoments[iMoment] = clusInfo.engCalibIn.engSmp[CaloSampling::PreSamplerB];
752 break;
754 myMoments[iMoment] = clusInfo.engCalibIn.engSmp[CaloSampling::PreSamplerE];
755 break;
757 myMoments[iMoment] = clusInfo.engCalibIn.engSmp[CaloSampling::TileGap3];
758 break;
760 myMoments[iMoment] = eng_calib_dead_tot;
761 break;
763 myMoments[iMoment] = eng_calib_dead_emb0;
764 break;
766 myMoments[iMoment] = eng_calib_dead_tile0;
767 break;
769 myMoments[iMoment] = eng_calib_dead_tileg3;
770 break;
772 myMoments[iMoment] = eng_calib_dead_eme0;
773 break;
775 myMoments[iMoment] = eng_calib_dead_hec0;
776 break;
778 myMoments[iMoment] = eng_calib_dead_fcal;
779 break;
781 myMoments[iMoment] = eng_calib_dead_leakage;
782 break;
784 myMoments[iMoment] = eng_calib_dead_unclass;
785 break;
787 myMoments[iMoment] = engCalibFrac[kCalibFracEM];
788 break;
790 myMoments[iMoment] = engCalibFrac[kCalibFracHAD];
791 break;
793 myMoments[iMoment] = engCalibFrac[kCalibFracREST];
794 break;
795 default:
796 // nothing to be done for other moments
797 break;
798 }
799
800 theCluster->insertMoment(vMomentsIter->second, myMoments[iMoment]);
801 }
802 }
803 }
804
805 return StatusCode::SUCCESS;
806}
807
808
809
810/* ****************************************************************************
811
812**************************************************************************** */
814{
815 double eta = fabs(x);
816 double ff;
817 if(eta<1.6){
818 ff = atan(5.0*1.7/(200.0*cosh(eta)));
819 }else if(eta<3.2){
820 ff = atan(5.0*1.6/(420./tanh(eta)));
821 }else{
822 ff = atan(5.0*0.95/(505./tanh(eta)));
823 }
824 return ff*(1./atan(5.0*1.7/200.0));
825}
#define M_PI
Scalar eta() const
pseudorapidity method
Scalar phi() const
phi method
#define endmsg
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_ERROR(x)
#define ATH_MSG_INFO(x)
#define ATH_MSG_WARNING(x)
#define ATH_MSG_DEBUG(x)
bool isValid(const T &p)
Av: we implement here an ATLAS-sepcific convention: all particles which are 99xxxxx are fine.
Definition AtlasPID.h:878
Definition of CaloDetDescrManager.
static const std::vector< std::string > bins
Handle class for reading from StoreGate.
#define x
constexpr int pow(int base, int exp) noexcept
AthAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor with parameters:
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
const ServiceHandle< StoreGateSvc > & detStore() const
MsgStream & msg() const
Class to define range of valid bins in eta x phi plane.
Class to store cluster number and weight for calorimeter cells.
Class to store cluster's calibration energies.
std::array< double, CaloDmDescrArea::DMA_MAX > engCalibDeadInArea
SG::ReadHandleKeyArray< CaloCalibrationHitContainer > m_DMCalibrationHitContainerNames
vector of dead material calibration hit container names to use.
moment_name_vector m_validNames
vector holding the names of valid moments which can be calculated.
std::set< xAOD::CaloCluster::MomentType > m_momentsAOD
set holding the list of moment enums which go in the first store i.e.
virtual StatusCode execute(const EventContext &ctx, xAOD::CaloClusterContainer *theClusColl) const override
Execute on an entire collection of clusters.
moment_name_set m_validMoments
set of moments which will be calculated.
SG::ReadHandleKeyArray< CaloCalibrationHitContainer > m_CalibrationHitContainerNames
vector of calibration hit container names to use.
std::map< Identifier, MyCellInfo > CellInfoSet_t
SG::ReadHandleKey< xAOD::TruthParticleContainer > m_truthParticleContainerKey
ReadHandleKey for truth particle container.
std::vector< std::string > m_momentsNames
vector holding the input list of names of moments to calculate.
SG::ReadCondHandleKey< CaloDetDescrManager > m_caloDetDescrMgrKey
Conditions Handle Key to access the CaloDetDescrManager.
std::array< std::vector< std::vector< CalibHitIPhiIEtaRange > >, 3 > m_i_phi_eta
CaloCalibClusterMomentsMaker2(const std::string &type, const std::string &name, const IInterface *parent)
const CaloDmDescrManager * m_caloDmDescrManager
std::vector< std::string > m_momentsNamesAOD
vector holding the list of moment names which go in the first store i.e.
std::pair< std::string, xAOD::CaloCluster::MomentType > moment_name_pair
typedef for a pair to index the enums defined in CaloClusterMoment with a string.
Class to store calorimeter calibration hit.
Data object for each calorimeter readout cell.
Definition CaloCell.h:57
Identifier ID() const
get ID (from cached data member) non-virtual and inline for fast access
Definition CaloCell.h:295
This class groups all DetDescr information related to a CaloCell.
const CaloDetDescrElement * get_element(const Identifier &cellId) const
get element by its identifier
This class provides the client interface for accessing the detector description information common to...
static const CaloDmDescrManager * instance()
std::vector< short > m_CaloSampleNeighbours
std::vector< float > m_CaloSampleEtaMin
std::vector< float > m_CaloSampleEtaMax
const T * at(size_type n) const
Access an element, as an rvalue.
DataModel_detail::iterator< DataVector > iterator
Definition DataVector.h:842
const_iterator end() const noexcept
Return a const_iterator pointing past the end of the collection.
const_iterator begin() const noexcept
Return a const_iterator pointing at the beginning of the collection.
size_type size() const noexcept
Returns the number of elements in the collection.
Property holding a SG store/key/clid from which a ReadHandle is made.
virtual bool isValid() override final
Can the handle be successfully dereferenced?
const_pointer_type cptr()
Dereference the pointer.
virtual double eta() const
The pseudorapidity ( ) of the particle.
virtual double e() const
The total energy of the particle.
void insertMoment(MomentType type, double value)
CaloClusterCellLink::const_iterator const_cell_iterator
Iterator of the underlying CaloClusterCellLink (explicitly const version)
const_cell_iterator cell_end() const
virtual double phi() const
The azimuthal angle ( ) of the particle.
@ ENG_CALIB_OUT_M
Attached Calibration Hit energy outside clusters but inside the calorimeter with medium matching (Ang...
@ ENG_CALIB_OUT_L
Attached Calibration Hit energy outside clusters but inside the calorimeter with loose matching (Angl...
@ ENG_CALIB_DEAD_UNCLASS
Attached Calibration Hit energy in dead material in unclassified areas of the detector.
@ ENG_CALIB_DEAD_HEC0
Attached Calibration Hit energy in dead material between EME3 and HEC0.
@ ENG_CALIB_DEAD_TILEG3
Attached Calibration Hit energy in dead material before scintillator.
@ ENG_CALIB_FRAC_REST
Calibration Hit energy inside the cluster caused by other particles.
@ ENG_CALIB_DEAD_EME0
Attached Calibration Hit energy in dead material before EME0, between EME0 and EME1.
@ ENG_CALIB_DEAD_TILE0
Attached Calibration Hit energy in dead material between EMB3 and TILE0.
@ ENG_CALIB_FRAC_EM
Calibration Hit energy inside the cluster caused by e/gamma/pi0.
@ ENG_CALIB_DEAD_FCAL
Attached Calibration Hit energy in dead material before FCAL, between FCAL and HEC.
@ ENG_CALIB_TOT
Calibration Hit energy inside the cluster.
@ ENG_CALIB_OUT_T
Attached Calibration Hit energy outside clusters but inside the calorimeter with tight matching (Angl...
@ ENG_CALIB_DEAD_LEAKAGE
Attached Calibration Hit energy in dead material behind calorimeters.
@ ENG_CALIB_FRAC_HAD
Calibration Hit energy inside the cluster caused by charged pi+ and pi-.
@ ENG_CALIB_EMB0
Calibration Hit energy inside the cluster barrel presampler.
@ ENG_CALIB_EME0
Calibration Hit energy inside the cluster endcap presampler.
@ ENG_CALIB_DEAD_EMB0
Attached Calibration Hit energy in dead material before EMB0, between EMB0 and EMB1.
@ ENG_CALIB_DEAD_TOT
Attached Calibration Hit energy in dead material.
@ ENG_CALIB_TILEG3
Calibration Hit energy inside the cluster scintillator.
const_cell_iterator cell_begin() const
Iterator of the underlying CaloClusterCellLink (const version)
int r
Definition globals.cxx:22
constexpr int INVALID_PARTICLE_ID
int uniqueID(const T &p)
constexpr int UNDEFINED_ID
CaloCluster_v1 CaloCluster
Define the latest version of the calorimeter cluster class.
CaloClusterContainer_v1 CaloClusterContainer
Define the latest version of the calorimeter cluster container.