EgammaIsolationAlgos/src/EgammaHadTower.cc

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     // Basic clusters of hybrid super-cluster do not have the seed set; take the first DetId instead
0014     // Should be checked . The single Tower Mode should be favoured until fixed
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   // in this mode, check only the tower behind the seed
0034   if (mode == HoeMode::SingleTower) {
0035     towers.push_back(towerOf(*sc.seed(), towerMap));
0036   }
0037
0038   // in this mode check the towers behind each basic cluster
0039   if (mode == HoeMode::TowersBehindCluster) {
0040     // Loop on the basic clusters
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       // Get the tower
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   //  if(towers.size() > 4) {
0058   //    std::cout << " NTOWERS " << towers.size() << " ";
0059   //    for(unsigned i=0; i<towers.size() ; ++i) {
0060   //      std::cout << towers[i] << " ";
0061   //    }
0062   //    std::cout <<  std::endl;
0063   //    for ( unsigned iclus=0 ; iclus < orderedClusters.size(); ++iclus) {
0064   //      // Get the tower
0065   //      CaloTowerDetId id = towerOf(*(orderedClusters[iclus]));
0066   //      std::cout << " Pos " << orderedClusters[iclus]->position() << " " << orderedClusters[iclus]->energy() << " " << id ;
0067   //    }
0068   //    std::cout << std::endl;
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       // Sunanda's fix for 2017 Plan1
0119       // and removed protection
0120       int status =
0121           hcalQuality.getValues((DetId)(hcalTopology.idFront(HcalDetId(id))), /*throwOnFail=*/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 }