File indexing completed on 2024-04-06 12:27:22
0001
0002 #ifndef RecoParticleFlow_PFClusterProducer_PFRecHitCaloNavigatorWithTime_h
0003 #define RecoParticleFlow_PFClusterProducer_PFRecHitCaloNavigatorWithTime_h
0004
0005 #include <memory>
0006
0007 #include "FWCore/Framework/interface/ConsumesCollector.h"
0008 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0009 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0010 #include "RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h"
0011
0012 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
0013 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0014 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0015 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0016 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0017
0018 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0019 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0020 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
0021 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0022
0023 #include "Geometry/CaloTopology/interface/CaloTowerTopology.h"
0024 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0025
0026 #include "RecoParticleFlow/PFClusterProducer/interface/CaloRecHitResolutionProvider.h"
0027
0028 template <typename D, typename T, bool ownsTopo = true>
0029 class PFRecHitCaloNavigatorWithTime : public PFRecHitNavigatorBase {
0030 public:
0031 PFRecHitCaloNavigatorWithTime(const edm::ParameterSet& iConfig, edm::ConsumesCollector& cc) {
0032 sigmaCut2_ = pow(iConfig.getParameter<double>("sigmaCut"), 2);
0033 const edm::ParameterSet& timeResConf = iConfig.getParameterSet("timeResolutionCalc");
0034 _timeResolutionCalc = std::make_unique<CaloRecHitResolutionProvider>(timeResConf);
0035 }
0036
0037 ~PFRecHitCaloNavigatorWithTime() override {
0038 if (!ownsTopo) {
0039 topology_.release();
0040 }
0041 }
0042
0043 void associateNeighbours(reco::PFRecHit& hit,
0044 std::unique_ptr<reco::PFRecHitCollection>& hits,
0045 edm::RefProd<reco::PFRecHitCollection>& refProd) override {
0046 DetId detid(hit.detId());
0047
0048 CaloNavigator<D> navigator(detid, topology_.get());
0049
0050 DetId N(0);
0051 DetId E(0);
0052 DetId S(0);
0053 DetId W(0);
0054 DetId NW(0);
0055 DetId NE(0);
0056 DetId SW(0);
0057 DetId SE(0);
0058
0059 N = navigator.north();
0060 associateNeighbour(N, hit, hits, refProd, 0, 1);
0061
0062 if (N != DetId(0)) {
0063 NE = navigator.east();
0064 } else {
0065 navigator.home();
0066 E = navigator.east();
0067 NE = navigator.north();
0068 }
0069 associateNeighbour(NE, hit, hits, refProd, 1, 1);
0070 navigator.home();
0071
0072 S = navigator.south();
0073 associateNeighbour(S, hit, hits, refProd, 0, -1);
0074
0075 if (S != DetId(0)) {
0076 SW = navigator.west();
0077 } else {
0078 navigator.home();
0079 W = navigator.west();
0080 SW = navigator.south();
0081 }
0082 associateNeighbour(SW, hit, hits, refProd, -1, -1);
0083 navigator.home();
0084
0085 E = navigator.east();
0086 associateNeighbour(E, hit, hits, refProd, 1, 0);
0087
0088 if (E != DetId(0)) {
0089 SE = navigator.south();
0090 } else {
0091 navigator.home();
0092 S = navigator.south();
0093 SE = navigator.east();
0094 }
0095 associateNeighbour(SE, hit, hits, refProd, 1, -1);
0096 navigator.home();
0097
0098 W = navigator.west();
0099 associateNeighbour(W, hit, hits, refProd, -1, 0);
0100
0101 if (W != DetId(0)) {
0102 NW = navigator.north();
0103 } else {
0104 navigator.home();
0105 N = navigator.north();
0106 NW = navigator.west();
0107 }
0108 associateNeighbour(NW, hit, hits, refProd, -1, 1);
0109 }
0110
0111 protected:
0112 double sigmaCut2_;
0113 std::unique_ptr<const T> topology_;
0114 std::unique_ptr<CaloRecHitResolutionProvider> _timeResolutionCalc;
0115
0116 void associateNeighbour(const DetId& id,
0117 reco::PFRecHit& hit,
0118 std::unique_ptr<reco::PFRecHitCollection>& hits,
0119 edm::RefProd<reco::PFRecHitCollection>& refProd,
0120 short eta,
0121 short phi) {
0122 double sigma2 = 10000.0;
0123
0124 auto found_hit = std::lower_bound(hits->begin(), hits->end(), id, [](const reco::PFRecHit& a, DetId b) {
0125 if (DetId(a.detId()).det() == DetId::Hcal || b.det() == DetId::Hcal)
0126 return (HcalDetId)(a.detId()) < (HcalDetId)(b);
0127
0128 else
0129 return a.detId() < b.rawId();
0130 });
0131
0132 if (found_hit != hits->end() &&
0133 ((found_hit->detId() == id.rawId()) ||
0134 ((id.det() == DetId::Hcal) && ((HcalDetId)(found_hit->detId()) == (HcalDetId)(id))))) {
0135 sigma2 = _timeResolutionCalc->timeResolution2(hit.energy()) +
0136 _timeResolutionCalc->timeResolution2(found_hit->energy());
0137 const double deltaTime = hit.time() - found_hit->time();
0138 if (deltaTime * deltaTime < sigmaCut2_ * sigma2) {
0139 hit.addNeighbour(eta, phi, 0, found_hit - hits->begin());
0140 }
0141 }
0142 }
0143 };
0144
0145 #endif
