File indexing completed on 2023-11-21 00:29:05
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
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
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
0078
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
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
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
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 {
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
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
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
0180 const float maxDepth = _param_X0 * (clusterT0 + vdt::fast_log(cl_energy_float));
0181 const float maxToFront = center_cell->getPosition().mag();
0182
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
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 }