PFClusterProducer/plugins/ECAL2DPositionCalcWithDepthCorr.cc

0001 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0002 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0003 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h"
0004 #include "FWCore/Framework/interface/ConsumesCollector.h"
0005 #include "FWCore/Framework/interface/ESHandle.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "FWCore/Utilities/interface/isFinite.h"
0009 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0010 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0011 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0012 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0013 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0014 #include "RecoParticleFlow/PFClusterProducer/interface/CaloRecHitResolutionProvider.h"
0015 #include "RecoParticleFlow/PFClusterProducer/interface/PFCPositionCalculatorBase.h"
0016
0017 #include "vdt/vdtMath.h"
0018
0019 #include <cmath>
0020 #include <memory>
0021 #include <unordered_map>
0022
0023 /// This is EGM version of the ECAL position + depth correction calculation
0024 class ECAL2DPositionCalcWithDepthCorr : public PFCPositionCalculatorBase {
0025 public:
0026   ECAL2DPositionCalcWithDepthCorr(const edm::ParameterSet& conf, edm::ConsumesCollector& cc)
0027       : PFCPositionCalculatorBase(conf, cc),
0028         _param_T0_EB(conf.getParameter<double>("T0_EB")),
0029         _param_T0_EE(conf.getParameter<double>("T0_EE")),
0030         _param_T0_ES(conf.getParameter<double>("T0_ES")),
0031         _param_W0(conf.getParameter<double>("W0")),
0032         _param_X0(conf.getParameter<double>("X0")),
0033         _minAllowedNorm(conf.getParameter<double>("minAllowedNormalization")),
0034         _ebGeom(nullptr),
0035         _eeGeom(nullptr),
0036         _esGeom(nullptr),
0037         _esPlus(false),
0038         _esMinus(false),
0039         _geomToken(cc.esConsumes<edm::Transition::BeginRun>()) {
0040     _timeResolutionCalc.reset(nullptr);
0041     if (conf.exists("timeResolutionCalc")) {
0042       const edm::ParameterSet& timeResConf = conf.getParameterSet("timeResolutionCalc");
0043       _timeResolutionCalc = std::make_unique<CaloRecHitResolutionProvider>(timeResConf);
0044     }
0045   }
0046   ECAL2DPositionCalcWithDepthCorr(const ECAL2DPositionCalcWithDepthCorr&) = delete;
0047   ECAL2DPositionCalcWithDepthCorr& operator=(const ECAL2DPositionCalcWithDepthCorr&) = delete;
0048
0049   void update(const edm::EventSetup& es) override;
0050
0051   void calculateAndSetPosition(reco::PFCluster&, const HcalPFCuts*) override;
0052   void calculateAndSetPositions(reco::PFClusterCollection&, const HcalPFCuts*) override;
0053
0054 private:
0055   const double _param_T0_EB;
0056   const double _param_T0_EE;
0057   const double _param_T0_ES;
0058   const double _param_W0;
0059   const double _param_X0;
0060   const double _minAllowedNorm;
0061
0062   //const CaloGeometryRecord  _caloGeom;
0063   const CaloSubdetectorGeometry* _ebGeom;
0064   const CaloSubdetectorGeometry* _eeGeom;
0065   const CaloSubdetectorGeometry* _esGeom;
0066   bool _esPlus, _esMinus;
0067
0068   std::unique_ptr<CaloRecHitResolutionProvider> _timeResolutionCalc;
0069
0070   void calculateAndSetPositionActual(reco::PFCluster&) const;
0071
0072   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> _geomToken;
0073 };
0074
0075 DEFINE_EDM_PLUGIN(PFCPositionCalculatorFactory, ECAL2DPositionCalcWithDepthCorr, "ECAL2DPositionCalcWithDepthCorr");
0076
0077 // faithful reimplementation of RecoEcal/EgammaCoreTools PositionCalc
0078 // sorry Stefano
0079
0080 void ECAL2DPositionCalcWithDepthCorr::update(const edm::EventSetup& es) {
0081   edm::ESHandle<CaloGeometry> geohandle = es.getHandle(_geomToken);
0082   _ebGeom = geohandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0083   _eeGeom = geohandle->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0084   _esGeom = geohandle->getSubdetectorGeometry(DetId::Ecal, EcalPreshower);
0085   if (_esGeom) {
0086     // ripped from RecoEcal/EgammaCoreTools
0087     for (uint32_t ic = 0; ic < _esGeom->getValidDetIds().size() && (!_esPlus || !_esMinus); ++ic) {
0088       const double z = _esGeom->getGeometry(_esGeom->getValidDetIds()[ic])->getPosition().z();
0089       _esPlus = _esPlus || (0 < z);
0090       _esMinus = _esMinus || (0 > z);
0091     }
0092   }
0093 }
0094
0095 void ECAL2DPositionCalcWithDepthCorr::calculateAndSetPosition(reco::PFCluster& cluster, const HcalPFCuts*) {
0096   calculateAndSetPositionActual(cluster);
0097 }
0098
0099 void ECAL2DPositionCalcWithDepthCorr::calculateAndSetPositions(reco::PFClusterCollection& clusters, const HcalPFCuts*) {
0100   for (reco::PFCluster& cluster : clusters) {
0101     calculateAndSetPositionActual(cluster);
0102   }
0103 }
0104
0105 void ECAL2DPositionCalcWithDepthCorr::calculateAndSetPositionActual(reco::PFCluster& cluster) const {
0106   constexpr double preshowerStartEta = 1.653;
0107   constexpr double preshowerEndEta = 2.6;
0108   if (!cluster.seed()) {
0109     throw cms::Exception("ClusterWithNoSeed") << " Found a cluster with no seed: " << cluster;
0110   }
0111   double cl_energy = 0;
0112   double cl_energy_float = 0;
0113   double cl_time = 0;
0114   double cl_timeweight = 0.0;
0115   double max_e = 0.0;
0116   double clusterT0 = 0.0;
0117   PFLayer::Layer max_e_layer = PFLayer::NONE;
0118   reco::PFRecHitRef refmax;
0119   // find the seed and max layer
0120   for (const reco::PFRecHitFraction& rhf : cluster.recHitFractions()) {
0121     const reco::PFRecHitRef& refhit = rhf.recHitRef();
0122     const double rh_fraction = rhf.fraction();
0123     const double rh_rawenergy = refhit->energy();
0124     const double rh_energy = rh_rawenergy * rh_fraction;
0125     const double rh_energyf = ((float)rh_rawenergy) * ((float)rh_fraction);
0126     if (!edm::isFinite(rh_energy)) {
0127       throw cms::Exception("PFClusterAlgo") << "rechit " << refhit->detId() << " has non-finite energy... "
0128                                             << "The input of the particle flow clustering seems to be corrupted.";
0129     }
0130     cl_energy += rh_energy;
0131     cl_energy_float += rh_energyf;
0132     // If time resolution is given, calculate weighted average
0133     if (_timeResolutionCalc) {
0134       const double res2 = 1. / _timeResolutionCalc->timeResolution2(rh_rawenergy);
0135       cl_time += rh_fraction * refhit->time() * res2;
0136       cl_timeweight += rh_fraction * res2;
0137     } else {  // assume resolution ~ 1/E**2
0138       const double rh_rawenergy2 = rh_rawenergy * rh_rawenergy;
0139       cl_timeweight += rh_rawenergy2 * rh_fraction;
0140       cl_time += rh_rawenergy2 * rh_fraction * refhit->time();
0141     }
0142     if (rh_energy > max_e) {
0143       max_e = rh_energy;
0144       max_e_layer = rhf.recHitRef()->layer();
0145       refmax = refhit;
0146     }
0147   }
0148   cluster.setEnergy(cl_energy);
0149   cluster.setTime(cl_time / cl_timeweight);
0150   if (_timeResolutionCalc) {
0151     cluster.setTimeError(std::sqrt(1.0f / float(cl_timeweight)));
0152   }
0153   cluster.setLayer(max_e_layer);
0154   const CaloSubdetectorGeometry* ecal_geom = nullptr;
0155   // get seed geometry information
0156   switch (max_e_layer) {
0157     case PFLayer::ECAL_BARREL:
0158       ecal_geom = _ebGeom;
0159       clusterT0 = _param_T0_EB;
0160       break;
0161     case PFLayer::ECAL_ENDCAP:
0162       ecal_geom = _eeGeom;
0163       clusterT0 = _param_T0_EE;
0164       break;
0165     default:
0166       throw cms::Exception("InvalidLayer") << "ECAL Position Calc only accepts ECAL_BARREL or ECAL_ENDCAP";
0167   }
0168
0169   auto center_cell = ecal_geom->getGeometry(refmax->detId());
0170   const double ctreta = center_cell->etaPos();
0171   const double actreta = std::abs(ctreta);
0172   // need to change T0 if in ES
0173   if (actreta > preshowerStartEta && actreta < preshowerEndEta) {
0174     if (ctreta > 0 && _esPlus)
0175       clusterT0 = _param_T0_ES;
0176     if (ctreta < 0 && _esMinus)
0177       clusterT0 = _param_T0_ES;
0178   }
0179   // floats to reproduce exactly the EGM code
0180   const float maxDepth = _param_X0 * (clusterT0 + vdt::fast_log(cl_energy_float));
0181   const float maxToFront = center_cell->getPosition().mag();
0182   // calculate the position
0183   const double logETot_inv = -vdt::fast_log(cl_energy_float);
0184   double position_norm = 0.0;
0185   double x(0.0), y(0.0), z(0.0);
0186   for (const reco::PFRecHitFraction& rhf : cluster.recHitFractions()) {
0187     double weight = 0.0;
0188     const reco::PFRecHitRef& refhit = rhf.recHitRef();
0189     const double rh_energy = ((float)refhit->energy()) * ((float)rhf.fraction());
0190     if (rh_energy > 0.0)
0191       weight = std::max(0.0, (_param_W0 + vdt::fast_log(rh_energy) + logETot_inv));
0192     auto cell = ecal_geom->getGeometry(refhit->detId());
0193     const float depth = maxDepth + maxToFront - cell->getPosition().mag();
0194     const GlobalPoint pos = static_cast<const TruncatedPyramid*>(cell.get())->getPosition(depth);
0195
0196     x += weight * pos.x();
0197     y += weight * pos.y();
0198     z += weight * pos.z();
0199
0200     position_norm += weight;
0201   }
0202
0203   // FALL BACK to LINEAR WEIGHTS
0204   if (position_norm == 0.) {
0205     for (const reco::PFRecHitFraction& rhf : cluster.recHitFractions()) {
0206       double weight = 0.0;
0207       const reco::PFRecHitRef& refhit = rhf.recHitRef();
0208       const double rh_energy = ((float)refhit->energy()) * ((float)rhf.fraction());
0209       if (rh_energy > 0.0)
0210         weight = rh_energy / cluster.energy();
0211
0212       auto cell = ecal_geom->getGeometry(refhit->detId());
0213       const float depth = maxDepth + maxToFront - cell->getPosition().mag();
0214       const GlobalPoint pos = cell->getPosition(depth);
0215
0216       x += weight * pos.x();
0217       y += weight * pos.y();
0218       z += weight * pos.z();
0219
0220       position_norm += weight;
0221     }
0222   }
0223
0224   if (position_norm < _minAllowedNorm) {
0225     edm::LogError("WeirdClusterNormalization") << "PFCluster too far from seeding cell: set position to (0,0,0).";
0226     cluster.setPosition(math::XYZPoint(0, 0, 0));
0227   } else {
0228     const double norm_inverse = 1.0 / position_norm;
0229     x *= norm_inverse;
0230     y *= norm_inverse;
0231     z *= norm_inverse;
0232
0233     cluster.setPosition(math::XYZPoint(x, y, z));
0234     cluster.calculatePositionREP();
0235   }
0236 }