Execute on an entire collection of clusters.
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;
287 for (const SG::ReadHandleKey<CaloCalibrationHitContainer>& key :
289 {
290 SG::ReadHandle<CaloCalibrationHitContainer> cchc (key, ctx);
291 if ( !cchc.isValid() ) {
293
294
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;
305 for (const SG::ReadHandleKey<CaloCalibrationHitContainer>& key :
307 {
308 SG::ReadHandle<CaloCalibrationHitContainer> cchc (key, ctx);
309 if ( !cchc.isValid() ) {
311
312
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
324
325 if ( !foundAllContainers ) return StatusCode::SUCCESS;
327
328
330
332
333
334
335
336
339 int iClus = 0;
340 for( ;clusIter!=clusIterEnd;++clusIter,++iClus) {
342
343
346 for(; cellIter != cellIterEnd; cellIter++ ){
347 const CaloCell* pCell = (*cellIter);
348 Identifier myId = pCell->
ID();
349
351 CellInfoSet_t::iterator bookmark = cellInfo.lower_bound( myId );
352 if(bookmark == cellInfo.end() || bookmark->first != myId) {
353
354 if (bookmark != cellInfo.begin()) --bookmark;
355 cellInfo.emplace_hint (bookmark, myId, std::move(info));
356 }else{
357
358 bookmark->second.Add( info );
359 }
360 }
361 }
362
363
364
365
366
367 unsigned int nHitsTotal = 0;
368 unsigned int nHitsWithoutParticleUID = 0;
370
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
378 for (
const std::pair<int, double>& p :
pos->second) {
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
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
409
411
412 if (doCalibFrac && !truthParticleContainerReadHandle.isValid()){
414 doCalibFrac = false;
415 }
416
417 std::vector<double> engCalibOut[3];
418
419
420
421
422 std::vector<std::vector <int > > clusListL;
423 std::vector<std::vector <int > > clusListM;
424 std::vector<std::vector <int > > clusListT;
431
432 for(unsigned int ii=0;ii<3;ii++) {
436 engCalibOut[ii].resize(theClusColl->size(),0);
437 iClus = -1;
439 ++iClus;
441
442 if ( clusInfo.engCalibIn.engTot > 0 ) {
443 int iEtaSign = 1;
444 if ( theCluster->
eta() < 0 ) iEtaSign = -1;
450 unsigned int iEtaBin(jeta);
451 if ( jeta < 0 ) iEtaBin = abs(jeta)-1;
452 const std::vector<CalibHitIPhiIEtaRange>&
bins
455 int jp =
range.iPhi+jphi;
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++ ) {
464 }
465 }
466 }
467 }
468 }
469 }
470 }
471
472
473
474
477
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
484 const CaloDetDescrElement* myCDDE =
488 if ( myCDDE ) {
494 for (unsigned int ii=0;ii<3;ii++) {
495 std::vector<std::vector <int > > *pClusList=nullptr;
497 pClusList = &clusListL;
499 pClusList = &clusListM;
501 pClusList = &clusListT;
502 if ( pClusList) {
503
504 double hitClusNorm(0.);
505 std::vector<int> hitClusIndex;
506 std::vector<double > hitClusEffEnergy;
507
511
512
513 auto pos = clusInfo.engCalibParticle.find(uniqueID);
514 if(pos!=clusInfo.engCalibParticle.end()) {
515
516
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 }
535 }
536 }
537 }
538
539
540
541
542
543
544 std::vector<std::vector <int > > *pClusList=nullptr;
546 pClusList = &clusListL;
548 pClusList = &clusListM;
550 pClusList = &clusListT;
551 }
553
555 for (const CaloCalibrationHit* hit : *dmcchc) {
556 Identifier myId = hit->cellID();
558 CaloDmDescrElement* myCDDE(nullptr);
560 if ( myCDDE ) {
563
569
570
571
572
575 int hitClusVecIndex = 0;
576 std::vector<int> hitClusIndex(theClusColl->size());
577 std::vector<double > hitClusEffEnergy(theClusColl->size());
578 double hitClusNorm = 0.0;
579
583
585
586
587 auto pos = clusInfo.engCalibParticle.find(uniqueID);
588 if(pos!=clusInfo.engCalibParticle.end()) {
589 double engClusTruthUniqueIDCalib =
pos->second.engTot;
590
592 double sum_smp_energy = 0.0;
593
598
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){
606 }
else if (phi_diff >
M_PI){
608 }
609 double eta_diff = (myCDDE->
eta()-theCluster->
eta());
610 float distance=sqrt(eta_diff * eta_diff + phi_diff * phi_diff);
611
613
614 hitClusIndex [hitClusVecIndex] = iClus;
615 hitClusEffEnergy [hitClusVecIndex++]= effEner;
616 hitClusNorm += effEner;
617 }
618 }
619 }
620 }
621 hitClusIndex.resize(hitClusVecIndex);
622 hitClusEffEnergy.resize(hitClusVecIndex);
623
624
625
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;
637 }
638 }
639
640 }
641 }
642 }
643 }
644 }
645
646 }
647
648
649 std::vector<double> engCalibFrac;
651
652 std::map<unsigned int,int> truthIDToPdgCodeMap;
653
654
655 for ( const auto *thisTruthParticle : *truthParticleContainerReadHandle){
656
657 if (!thisTruthParticle){
659 continue;
660 }
661
662 truthIDToPdgCodeMap[
HepMC::uniqueID(thisTruthParticle)] = thisTruthParticle->pdgId();
663 }
664
665
666 for( clusIter = theClusColl->begin(),iClus=0; clusIter!=clusIterEnd;++clusIter,++iClus) {
669
670
672 + clusInfo.engCalibIn.engSmp[CaloSampling::PreSamplerB]
673 + clusInfo.engCalibIn.engSmp[CaloSampling::PreSamplerE]
674 + clusInfo.engCalibIn.engSmp[CaloSampling::TileGap3];
675
678 + clusInfo.engCalibIn.engSmp[CaloSampling::PreSamplerB];
679
681
683 + clusInfo.engCalibIn.engSmp[CaloSampling::TileGap3];
684
687 + clusInfo.engCalibIn.engSmp[CaloSampling::PreSamplerE];
688
690
693
695
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
702
703
705 if(clusInfo.engCalibIn.engTot > 0.0) {
706
707 for (const auto& p : clusInfo.engCalibParticle) {
708 int pdg_id = 0;
709 if (
auto it = truthIDToPdgCodeMap.find(
p.first); it != truthIDToPdgCodeMap.end()) {
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) {
717 } else if( pdg_id == 111 || pdg_id == 22 || std::abs(pdg_id)==11) {
719 } else {
721 }
722 }
723 for(
size_t i=0;
i<engCalibFrac.size();
i++) engCalibFrac[i] = engCalibFrac[i]/clusInfo.engCalibIn.engTot;
724 }
725 }
726
729
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
736 switch (vMomentsIter->second) {
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;
788 break;
791 break;
794 break;
795 default:
796
797 break;
798 }
799
800 theCluster->
insertMoment(vMomentsIter->second, myMoments[iMoment]);
801 }
802 }
803 }
804
805 return StatusCode::SUCCESS;
806}
#define ATH_MSG_WARNING(x)
CaloCalibrationHitContainer
static const std::vector< std::string > bins
constexpr int pow(int base, int exp) noexcept
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
std::vector< MyClusInfo > ClusInfo_t
moment_name_set m_validMoments
set of moments which will be calculated.
std::map< Identifier, MyCellInfo > CellInfoSet_t
SG::ReadHandleKey< xAOD::TruthParticleContainer > m_truthParticleContainerKey
ReadHandleKey for truth particle container.
SG::ReadCondHandleKey< CaloDetDescrManager > m_caloDetDescrMgrKey
Conditions Handle Key to access the CaloDetDescrManager.
Identifier ID() const
get ID (from cached data member) non-virtual and inline for fast access
weight_t weight() const
Accessor for weight associated to this cell.
float eta() const
cell eta
float phi() const
cell phi
const CaloDetDescrElement * get_element(const Identifier &cellId) const
get element by its identifier
std::vector< short > m_CaloSampleNeighbours
std::vector< float > m_CaloSampleEtaMin
std::vector< float > m_CaloSampleEtaMax
DataModel_detail::iterator< DataVector > iterator
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.
const_cell_iterator cell_begin() const
Iterator of the underlying CaloClusterCellLink (const version)
float distance(const Amg::Vector3D &p1, const Amg::Vector3D &p2)
calculates the distance between two point in 3D space
constexpr int INVALID_PARTICLE_ID
constexpr int UNDEFINED_ID
CaloCluster_v1 CaloCluster
Define the latest version of the calorimeter cluster class.
1.14.0