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