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       caloGeometryToken_{sumes.esConsumes()} {}
0026
0027 void ClusterTools::getEvent(const edm::Event& ev) {
0028   eerh_ = &ev.get(eetok);
0029   fhrh_ = &ev.get(fhtok);
0030   bhrh_ = &ev.get(bhtok);
0031 }
0032
0033 void ClusterTools::getEventSetup(const edm::EventSetup& es) { rhtools_.setGeometry(es.getData(caloGeometryToken_)); }
0034
0035 float ClusterTools::getClusterHadronFraction(const reco::CaloCluster& clus) const {
0036   float energy = 0.f, energyHad = 0.f;
0037   const auto& hits = clus.hitsAndFractions();
0038   for (const auto& hit : hits) {
0039     const auto& id = hit.first;
0040     const float fraction = hit.second;
0041     if (id.det() == DetId::HGCalEE) {
0042       energy += eerh_->find(id)->energy() * fraction;
0043     } else if (id.det() == DetId::HGCalHSi) {
0044       const float temp = fhrh_->find(id)->energy();
0045       energy += temp * fraction;
0046       energyHad += temp * fraction;
0047     } else if (id.det() == DetId::HGCalHSc) {
0048       const float temp = bhrh_->find(id)->energy();
0049       energy += temp * fraction;
0050       energyHad += temp * fraction;
0051     } else if (id.det() == DetId::Forward) {
0052       switch (id.subdetId()) {
0053         case HGCEE:
0054           energy += eerh_->find(id)->energy() * fraction;
0055           break;
0056         case HGCHEF: {
0057           const float temp = fhrh_->find(id)->energy();
0058           energy += temp * fraction;
0059           energyHad += temp * fraction;
0060         } break;
0061         default:
0062           throw cms::Exception("HGCalClusterTools") << " Cluster contains hits that are not from HGCal! " << std::endl;
0063       }
0064     } else if (id.det() == DetId::Hcal && id.subdetId() == HcalEndcap) {
0065       const float temp = bhrh_->find(id)->energy();
0066       energy += temp * fraction;
0067       energyHad += temp * fraction;
0068     } else {
0069       throw cms::Exception("HGCalClusterTools") << " Cluster contains hits that are not from HGCal! " << std::endl;
0070     }
0071   }
0072   float fraction = -1.f;
0073   if (energy > 0.f) {
0074     fraction = energyHad / energy;
0075   }
0076   return fraction;
0077 }
0078
0079 math::XYZPoint ClusterTools::getMultiClusterPosition(const reco::HGCalMultiCluster& clu) const {
0080   if (clu.clusters().empty())
0081     return math::XYZPoint();
0082
0083   double acc_x = 0.0;
0084   double acc_y = 0.0;
0085   double acc_z = 0.0;
0086   double totweight = 0.;
0087   double mcenergy = getMultiClusterEnergy(clu);
0088   for (const auto& ptr : clu.clusters()) {
0089     if (mcenergy != 0) {
0090       if (ptr->energy() < .01 * mcenergy)
0091         continue;  //cutoff < 1% layer contribution
0092     }
0093     const double weight = ptr->energy();  // weigth each corrdinate only by the total energy of the layer cluster
0094     acc_x += ptr->x() * weight;
0095     acc_y += ptr->y() * weight;
0096     acc_z += ptr->z() * weight;
0097     totweight += weight;
0098   }
0099   if (totweight != 0) {
0100     acc_x /= totweight;
0101     acc_y /= totweight;
0102     acc_z /= totweight;
0103   }
0104   // return x/y/z in absolute coordinates
0105   return math::XYZPoint(acc_x, acc_y, acc_z);
0106 }
0107
0108 int ClusterTools::getLayer(const DetId detid) const { return rhtools_.getLayerWithOffset(detid); }
0109
0110 double ClusterTools::getMultiClusterEnergy(const reco::HGCalMultiCluster& clu) const {
0111   double acc = 0.0;
0112   for (const auto& ptr : clu.clusters()) {
0113     acc += ptr->energy();
0114   }
0115   return acc;
0116 }
0117
0118 bool ClusterTools::getWidths(const reco::CaloCluster& clus,
0119                              double& sigmaetaeta,
0120                              double& sigmaphiphi,
0121                              double& sigmaetaetal,
0122                              double& sigmaphiphil) const {
0123   if (getLayer(clus.hitsAndFractions()[0].first) > (int)rhtools_.lastLayerEE())
0124     return false;
0125   const math::XYZPoint& position(clus.position());
0126   unsigned nhit = clus.hitsAndFractions().size();
0127
0128   sigmaetaeta = 0.;
0129   sigmaphiphi = 0.;
0130   sigmaetaetal = 0.;
0131   sigmaphiphil = 0.;
0132
0133   double sumw = 0.;
0134   double sumlogw = 0.;
0135
0136   for (unsigned int ih = 0; ih < nhit; ++ih) {
0137     const DetId& id = (clus.hitsAndFractions())[ih].first;
0138     if ((clus.hitsAndFractions())[ih].second == 0.)
0139       continue;
0140
0141     if ((id.det() == DetId::HGCalEE) || (id.det() == DetId::Forward && id.subdetId() == HGCEE)) {
0142       const HGCRecHit* theHit = &(*eerh_->find(id));
0143
0144       GlobalPoint cellPos = rhtools_.getPosition(id);
0145       double weight = theHit->energy();
0146       // take w0=2 To be optimized
0147       double logweight = 0;
0148       if (clus.energy() != 0) {
0149         logweight = std::max(0., 2 + log(theHit->energy() / clus.energy()));
0150       }
0151       double deltaetaeta2 = (cellPos.eta() - position.eta()) * (cellPos.eta() - position.eta());
0152       double deltaphiphi2 = (cellPos.phi() - position.phi()) * (cellPos.phi() - position.phi());
0153       sigmaetaeta += deltaetaeta2 * weight;
0154       sigmaphiphi += deltaphiphi2 * weight;
0155       sigmaetaetal += deltaetaeta2 * logweight;
0156       sigmaphiphil += deltaphiphi2 * logweight;
0157       sumw += weight;
0158       sumlogw += logweight;
0159     }
0160   }
0161
0162   if (sumw <= 0.)
0163     return false;
0164
0165   sigmaetaeta /= sumw;
0166   sigmaetaeta = std::sqrt(sigmaetaeta);
0167   sigmaphiphi /= sumw;
0168   sigmaphiphi = std::sqrt(sigmaphiphi);
0169
0170   if (sumlogw != 0) {
0171     sigmaetaetal /= sumlogw;
0172     sigmaetaetal = std::sqrt(sigmaetaetal);
0173     sigmaphiphil /= sumlogw;
0174     sigmaphiphil = std::sqrt(sigmaphiphil);
0175   }
0176
0177   return true;
0178 }