HGCalRecAlgos/src/ClusterTools.cc

0001 #include "RecoLocalCalo/HGCalRecAlgos/interface/ClusterTools.h"
0002
0003 #include "DataFormats/DetId/interface/DetId.h"
0004 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0005 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0006 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0007
0008 #include "FWCore/Framework/interface/ESHandle.h"
0009
0010 #include "FWCore/Framework/interface/Event.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "DataFormats/Common/interface/Handle.h"
0013
0014 #include "vdt/vdtMath.h"
0015
0016 #include <iostream>
0017
0018 using namespace hgcal;
0019 ClusterTools::ClusterTools() {}
0020
0021 ClusterTools::ClusterTools(const edm::ParameterSet& conf, edm::ConsumesCollector& sumes)
0022     : eetok(sumes.consumes<HGCRecHitCollection>(conf.getParameter<edm::InputTag>("HGCEEInput"))),
0023       fhtok(sumes.consumes<HGCRecHitCollection>(conf.getParameter<edm::InputTag>("HGCFHInput"))),
0024       bhtok(sumes.consumes<HGCRecHitCollection>(conf.getParameter<edm::InputTag>("HGCBHInput"))),
0025       hitMapToken_(sumes.consumes<std::unordered_map<DetId, const unsigned int>>(
0026           conf.getParameter<edm::InputTag>("hgcalHitMap"))),
0027       caloGeometryToken_{sumes.esConsumes()} {}
0028
0029 void ClusterTools::getEvent(const edm::Event& ev) {
0030   eerh_ = &ev.get(eetok);
0031   fhrh_ = &ev.get(fhtok);
0032   bhrh_ = &ev.get(bhtok);
0033   hitMap_ = &ev.get(hitMapToken_);
0034   rechitManager_ = std::make_unique<MultiVectorManager<HGCRecHit>>();
0035   rechitManager_->addVector(*eerh_);
0036   rechitManager_->addVector(*fhrh_);
0037   rechitManager_->addVector(*bhrh_);
0038 }
0039
0040 void ClusterTools::getEventSetup(const edm::EventSetup& es) { rhtools_.setGeometry(es.getData(caloGeometryToken_)); }
0041
0042 float ClusterTools::getClusterHadronFraction(const reco::CaloCluster& clus) const {
0043   float energy = 0.f, energyHad = 0.f;
0044   const auto& hits = clus.hitsAndFractions();
0045   const auto& rhmanager = *rechitManager_;
0046   for (const auto& hit : hits) {
0047     const auto& id = hit.first;
0048     const float fraction = hit.second;
0049     auto hitIter = hitMap_->find(id.rawId());
0050     if (hitIter == hitMap_->end()) {
0051       continue;
0052     }
0053     unsigned int rechitIndex = hitIter->second;
0054     float hitEnergy = rhmanager[rechitIndex].energy() * fraction;
0055     energy += hitEnergy;
0056     if (id.det() == DetId::HGCalHSi || id.det() == DetId::HGCalHSc ||
0057         (id.det() == DetId::Forward && id.subdetId() == HGCHEF) ||
0058         (id.det() == DetId::Hcal && id.subdetId() == HcalEndcap)) {
0059       energyHad += hitEnergy;
0060     }
0061   }
0062   float hadronicFraction = -1.f;
0063   if (energy > 0.f) {
0064     hadronicFraction = energyHad / energy;
0065   }
0066   return hadronicFraction;
0067 }
0068
0069 math::XYZPoint ClusterTools::getMultiClusterPosition(const reco::HGCalMultiCluster& clu) const {
0070   if (clu.clusters().empty())
0071     return math::XYZPoint();
0072
0073   double acc_x = 0.0;
0074   double acc_y = 0.0;
0075   double acc_z = 0.0;
0076   double totweight = 0.;
0077   double mcenergy = getMultiClusterEnergy(clu);
0078   for (const auto& ptr : clu.clusters()) {
0079     if (mcenergy != 0) {
0080       if (ptr->energy() < .01 * mcenergy)
0081         continue;  //cutoff < 1% layer contribution
0082     }
0083     const double weight = ptr->energy();  // weigth each corrdinate only by the total energy of the layer cluster
0084     acc_x += ptr->x() * weight;
0085     acc_y += ptr->y() * weight;
0086     acc_z += ptr->z() * weight;
0087     totweight += weight;
0088   }
0089   if (totweight != 0) {
0090     acc_x /= totweight;
0091     acc_y /= totweight;
0092     acc_z /= totweight;
0093   }
0094   // return x/y/z in absolute coordinates
0095   return math::XYZPoint(acc_x, acc_y, acc_z);
0096 }
0097
0098 int ClusterTools::getLayer(const DetId detid) const { return rhtools_.getLayerWithOffset(detid); }
0099
0100 double ClusterTools::getMultiClusterEnergy(const reco::HGCalMultiCluster& clu) const {
0101   double acc = 0.0;
0102   for (const auto& ptr : clu.clusters()) {
0103     acc += ptr->energy();
0104   }
0105   return acc;
0106 }
0107
0108 bool ClusterTools::getWidths(const reco::CaloCluster& clus,
0109                              double& sigmaetaeta,
0110                              double& sigmaphiphi,
0111                              double& sigmaetaetal,
0112                              double& sigmaphiphil) const {
0113   if (getLayer(clus.hitsAndFractions()[0].first) > (int)rhtools_.lastLayerEE())
0114     return false;
0115   const math::XYZPoint& position(clus.position());
0116   unsigned nhit = clus.hitsAndFractions().size();
0117
0118   sigmaetaeta = 0.;
0119   sigmaphiphi = 0.;
0120   sigmaetaetal = 0.;
0121   sigmaphiphil = 0.;
0122
0123   double sumw = 0.;
0124   double sumlogw = 0.;
0125   const auto& rhmanager = *rechitManager_;
0126
0127   for (unsigned int ih = 0; ih < nhit; ++ih) {
0128     const DetId& id = (clus.hitsAndFractions())[ih].first;
0129     if ((clus.hitsAndFractions())[ih].second == 0.)
0130       continue;
0131
0132     if ((id.det() == DetId::HGCalEE) || (id.det() == DetId::Forward && id.subdetId() == HGCEE)) {
0133       auto hitIter = hitMap_->find(id.rawId());
0134       if (hitIter == hitMap_->end()) {
0135         continue;
0136       }
0137       unsigned int rechitIndex = hitIter->second;
0138       float hitEnergy = rhmanager[rechitIndex].energy();
0139
0140       GlobalPoint cellPos = rhtools_.getPosition(id);
0141       double weight = hitEnergy;
0142       // take w0=2 To be optimized
0143       double logweight = 0;
0144       if (clus.energy() != 0) {
0145         logweight = std::max(0., 2 + std::log(hitEnergy / clus.energy()));
0146       }
0147       double deltaetaeta2 = (cellPos.eta() - position.eta()) * (cellPos.eta() - position.eta());
0148       double deltaphiphi2 = (cellPos.phi() - position.phi()) * (cellPos.phi() - position.phi());
0149       sigmaetaeta += deltaetaeta2 * weight;
0150       sigmaphiphi += deltaphiphi2 * weight;
0151       sigmaetaetal += deltaetaeta2 * logweight;
0152       sigmaphiphil += deltaphiphi2 * logweight;
0153       sumw += weight;
0154       sumlogw += logweight;
0155     }
0156   }
0157
0158   if (sumw <= 0.)
0159     return false;
0160
0161   sigmaetaeta /= sumw;
0162   sigmaetaeta = std::sqrt(sigmaetaeta);
0163   sigmaphiphi /= sumw;
0164   sigmaphiphi = std::sqrt(sigmaphiphi);
0165
0166   if (sumlogw != 0) {
0167     sigmaetaetal /= sumlogw;
0168     sigmaetaetal = std::sqrt(sigmaetaetal);
0169     sigmaphiphil /= sumlogw;
0170     sigmaphiphil = std::sqrt(sigmaphiphil);
0171   }
0172
0173   return true;
0174 }