File indexing completed on 2023-03-17 11:17:42
0001 #include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHadTower.h"
0002 #include "CondFormats/HcalObjects/interface/HcalChannelStatus.h"
0003 #include "CondFormats/HcalObjects/interface/HcalChannelQuality.h"
0004 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0005 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0006
0007 #include <iostream>
0008 #include <algorithm>
0009
0010 CaloTowerDetId egamma::towerOf(const reco::CaloCluster& cluster, CaloTowerConstituentsMap const& towerMap) {
0011 DetId detid = cluster.seed();
0012 if (detid.det() != DetId::Ecal) {
0013
0014
0015 detid = cluster.hitsAndFractions()[0].first;
0016 if (detid.det() != DetId::Ecal) {
0017 CaloTowerDetId tower;
0018 return tower;
0019 }
0020 }
0021 CaloTowerDetId id(towerMap.towerOf(detid));
0022 return id;
0023 }
0024
0025 std::vector<CaloTowerDetId> egamma::towersOf(const reco::SuperCluster& sc,
0026 CaloTowerConstituentsMap const& towerMap,
0027 HoeMode mode) {
0028 constexpr unsigned int nMaxClusters = 4;
0029
0030 std::vector<CaloTowerDetId> towers;
0031 std::vector<reco::CaloClusterPtr> orderedClusters;
0032
0033
0034 if (mode == HoeMode::SingleTower) {
0035 towers.push_back(towerOf(*sc.seed(), towerMap));
0036 }
0037
0038
0039 if (mode == HoeMode::TowersBehindCluster) {
0040
0041 for (auto it = sc.clustersBegin(); it != sc.clustersEnd(); ++it) {
0042 orderedClusters.push_back(*it);
0043 }
0044 std::sort(orderedClusters.begin(), orderedClusters.end(), [](auto& c1, auto& c2) { return (*c1 > *c2); });
0045 unsigned nclusters = orderedClusters.size();
0046 for (unsigned iclus = 0; iclus < nclusters && iclus < nMaxClusters; ++iclus) {
0047
0048 CaloTowerDetId id = towerOf(*(orderedClusters[iclus]), towerMap);
0049 #ifdef EDM_ML_DEBUG
0050 std::cout << "CaloTowerId " << id << std::endl;
0051 #endif
0052 if (std::find(towers.begin(), towers.end(), id) == towers.end()) {
0053 towers.push_back(id);
0054 }
0055 }
0056 }
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070 return towers;
0071 }
0072
0073 double egamma::depth1HcalESum(const std::vector<CaloTowerDetId>& towers, CaloTowerCollection const& towerCollection) {
0074 double esum = 0.;
0075 for (auto const& tower : towerCollection) {
0076 if (std::find(towers.begin(), towers.end(), tower.id()) != towers.end()) {
0077 esum += tower.ietaAbs() < 18 || tower.ietaAbs() > 29 ? tower.hadEnergy() : tower.hadEnergyHeInnerLayer();
0078 }
0079 }
0080 return esum;
0081 }
0082
0083 double egamma::depth2HcalESum(const std::vector<CaloTowerDetId>& towers, CaloTowerCollection const& towerCollection) {
0084 double esum = 0.;
0085 for (auto const& tower : towerCollection) {
0086 if (std::find(towers.begin(), towers.end(), tower.id()) != towers.end()) {
0087 esum += tower.hadEnergyHeOuterLayer();
0088 }
0089 }
0090 return esum;
0091 }
0092
0093 bool egamma::hasActiveHcal(const std::vector<CaloTowerDetId>& towers,
0094 CaloTowerConstituentsMap const& towerMap,
0095 const HcalChannelQuality& hcalQuality,
0096 HcalTopology const& hcalTopology) {
0097 bool active = false;
0098 int statusMask = ((1 << HcalChannelStatus::HcalCellOff) | (1 << HcalChannelStatus::HcalCellMask) |
0099 (1 << HcalChannelStatus::HcalCellDead));
0100 #ifdef EDM_ML_DEBUG
0101 std::cout << "DEBUG: hasActiveHcal called with " << towers.size() << " detids. First tower detid ieta "
0102 << towers.front().ieta() << " iphi " << towers.front().iphi() << std::endl;
0103 #endif
0104 for (auto towerid : towers) {
0105 unsigned int ngood = 0, nbad = 0;
0106 for (DetId id : towerMap.constituentsOf(towerid)) {
0107 if (id.det() != DetId::Hcal) {
0108 continue;
0109 }
0110 HcalDetId hid(id);
0111 if (hid.subdet() != HcalBarrel && hid.subdet() != HcalEndcap)
0112 continue;
0113 #ifdef EDM_ML_DEBUG
0114 std::cout << "EgammaHadTower DetId " << std::hex << id.rawId() << " hid.rawId " << hid.rawId() << std::dec
0115 << " sub " << hid.subdet() << " ieta " << hid.ieta() << " iphi " << hid.iphi() << " depth "
0116 << hid.depth() << std::endl;
0117 #endif
0118
0119
0120 int status =
0121 hcalQuality.getValues((DetId)(hcalTopology.idFront(HcalDetId(id))), true)->getValue();
0122
0123 #ifdef EDM_ML_DEBUG
0124 std::cout << "channels status = " << std::hex << status << std::dec << " int value = " << status << std::endl;
0125 #endif
0126
0127 if (status & statusMask) {
0128 #ifdef EDM_ML_DEBUG
0129 std::cout << " BAD!" << std::endl;
0130 #endif
0131 nbad++;
0132 } else {
0133 ngood++;
0134 }
0135 }
0136 #ifdef EDM_ML_DEBUG
0137 std::cout << " overall ngood " << ngood << " nbad " << nbad << "\n";
0138 #endif
0139 if (nbad == 0 || (ngood > 0 && nbad < ngood)) {
0140 active = true;
0141 }
0142 }
0143 return active;
0144 }