File indexing completed on 2024-04-06 12:27:28
0001 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0002 #include "DataFormats/ParticleFlowReco/interface/PFBlockElement.h"
0003 #include "RecoParticleFlow/PFClusterTools/interface/LinkByRecHit.h"
0004 #include "RecoParticleFlow/PFProducer/interface/KDTreeLinkerBase.h"
0005 #include "CommonTools/RecoAlgos/interface/KDTreeLinkerAlgo.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008
0009
0010
0011 class KDTreeLinkerTrackHcal : public KDTreeLinkerBase {
0012 public:
0013 KDTreeLinkerTrackHcal(const edm::ParameterSet& conf);
0014 ~KDTreeLinkerTrackHcal() override;
0015
0016
0017 void insertTargetElt(reco::PFBlockElement* track) override;
0018
0019
0020
0021
0022 void insertFieldClusterElt(reco::PFBlockElement* hcalCluster) override;
0023
0024
0025 void buildTree() override;
0026
0027
0028
0029
0030
0031 void searchLinks() override;
0032
0033
0034
0035 void updatePFBlockEltWithLinks() override;
0036
0037
0038 void clear() override;
0039
0040 private:
0041
0042 BlockEltSet targetSet_;
0043 BlockEltSet fieldClusterSet_;
0044
0045
0046 RecHitSet rechitsSet_;
0047
0048
0049 BlockElt2BlockEltMap target2ClusterLinks_;
0050
0051
0052 RecHit2BlockEltMap rechit2ClusterLinks_;
0053
0054
0055 KDTreeLinkerAlgo<reco::PFRecHit const*> tree_;
0056
0057
0058 std::string trajectoryLayerEntranceString_;
0059 std::string trajectoryLayerExitString_;
0060 reco::PFTrajectoryPoint::LayerType trajectoryLayerEntrance_;
0061 reco::PFTrajectoryPoint::LayerType trajectoryLayerExit_;
0062 bool checkExit_;
0063
0064
0065 int nMaxHcalLinksPerTrack_;
0066 };
0067
0068
0069
0070 DEFINE_EDM_PLUGIN(KDTreeLinkerFactory, KDTreeLinkerTrackHcal, "KDTreeTrackAndHCALLinker");
0071
0072 KDTreeLinkerTrackHcal::KDTreeLinkerTrackHcal(const edm::ParameterSet& conf)
0073 : KDTreeLinkerBase(conf),
0074 trajectoryLayerEntranceString_(conf.getParameter<std::string>("trajectoryLayerEntrance")),
0075 trajectoryLayerExitString_(conf.getParameter<std::string>("trajectoryLayerExit")),
0076 nMaxHcalLinksPerTrack_(conf.getParameter<int>("nMaxHcalLinksPerTrack")) {
0077
0078 cristalPhiEtaMaxSize_ = 0.2;
0079 phiOffset_ = 0.32;
0080
0081 trajectoryLayerEntrance_ = reco::PFTrajectoryPoint::layerTypeByName(trajectoryLayerEntranceString_);
0082 trajectoryLayerExit_ = reco::PFTrajectoryPoint::layerTypeByName(trajectoryLayerExitString_);
0083
0084 assert((trajectoryLayerEntrance_ == reco::PFTrajectoryPoint::HCALEntrance &&
0085 trajectoryLayerExit_ == reco::PFTrajectoryPoint::HCALExit) ||
0086 (trajectoryLayerEntrance_ == reco::PFTrajectoryPoint::HCALEntrance &&
0087 trajectoryLayerExit_ == reco::PFTrajectoryPoint::Unknown) ||
0088 (trajectoryLayerEntrance_ == reco::PFTrajectoryPoint::VFcalEntrance &&
0089 trajectoryLayerExit_ == reco::PFTrajectoryPoint::Unknown));
0090
0091 checkExit_ = trajectoryLayerExit_ != reco::PFTrajectoryPoint::Unknown;
0092 }
0093
0094 KDTreeLinkerTrackHcal::~KDTreeLinkerTrackHcal() { clear(); }
0095
0096 void KDTreeLinkerTrackHcal::insertTargetElt(reco::PFBlockElement* track) {
0097 if (track->trackRefPF()->extrapolatedPoint(trajectoryLayerEntrance_).isValid()) {
0098 targetSet_.insert(track);
0099 }
0100 }
0101
0102 void KDTreeLinkerTrackHcal::insertFieldClusterElt(reco::PFBlockElement* hcalCluster) {
0103 const reco::PFClusterRef& clusterref = hcalCluster->clusterRef();
0104
0105
0106
0107
0108
0109
0110 const std::vector<reco::PFRecHitFraction>& fraction = clusterref->recHitFractions();
0111
0112
0113 fieldClusterSet_.insert(hcalCluster);
0114 for (size_t rhit = 0; rhit < fraction.size(); ++rhit) {
0115 const reco::PFRecHitRef& rh = fraction[rhit].recHitRef();
0116 double fract = fraction[rhit].fraction();
0117
0118 if ((rh.isNull()) || (fract < cutOffFrac))
0119 continue;
0120
0121 const reco::PFRecHit& rechit = *rh;
0122
0123
0124 rechit2ClusterLinks_[&rechit].insert(hcalCluster);
0125
0126
0127 rechitsSet_.insert(&rechit);
0128 }
0129 }
0130
0131 void KDTreeLinkerTrackHcal::buildTree() {
0132
0133 std::vector<KDTreeNodeInfo<reco::PFRecHit const*, 2>> eltList;
0134
0135
0136 for (RecHitSet::const_iterator it = rechitsSet_.begin(); it != rechitsSet_.end(); it++) {
0137 const reco::PFRecHit::REPPoint& posrep = (*it)->positionREP();
0138
0139 KDTreeNodeInfo<reco::PFRecHit const*, 2> rh1(*it, posrep.eta(), posrep.phi());
0140 eltList.push_back(rh1);
0141
0142
0143
0144 if (rh1.dims[1] > (M_PI - phiOffset_)) {
0145 float phi = rh1.dims[1] - 2 * M_PI;
0146 KDTreeNodeInfo<reco::PFRecHit const*, 2> rh2(*it, float(posrep.eta()), phi);
0147 eltList.push_back(rh2);
0148 }
0149
0150 if (rh1.dims[1] < (M_PI * -1.0 + phiOffset_)) {
0151 float phi = rh1.dims[1] + 2 * M_PI;
0152 KDTreeNodeInfo<reco::PFRecHit const*, 2> rh3(*it, float(posrep.eta()), phi);
0153 eltList.push_back(rh3);
0154 }
0155 }
0156
0157
0158 float phimin = -1.0 * M_PI - phiOffset_;
0159 float phimax = M_PI + phiOffset_;
0160
0161
0162 KDTreeBox region(-3.0f, 3.0f, phimin, phimax);
0163
0164
0165 tree_.build(eltList, region);
0166 }
0167
0168 void KDTreeLinkerTrackHcal::searchLinks() {
0169
0170
0171
0172 for (BlockEltSet::iterator it = targetSet_.begin(); it != targetSet_.end(); it++) {
0173 reco::PFRecTrackRef trackref = (*it)->trackRefPF();
0174
0175 const reco::PFTrajectoryPoint& atHCAL = trackref->extrapolatedPoint(trajectoryLayerEntrance_);
0176
0177
0178 if (!atHCAL.isValid())
0179 continue;
0180
0181
0182 double dHeta = 0.0;
0183 float dHphi = 0.0;
0184 if (checkExit_) {
0185 const reco::PFTrajectoryPoint& atHCALExit = trackref->extrapolatedPoint(trajectoryLayerExit_);
0186 dHeta = atHCALExit.positionREP().eta() - atHCAL.positionREP().eta();
0187 dHphi = atHCALExit.positionREP().phi() - atHCAL.positionREP().phi();
0188 if (dHphi > M_PI)
0189 dHphi = dHphi - 2. * M_PI;
0190 else if (dHphi < -M_PI)
0191 dHphi = dHphi + 2. * M_PI;
0192 }
0193
0194 float tracketa = atHCAL.positionREP().eta() + 0.1 * dHeta;
0195 float trackphi = atHCAL.positionREP().phi() + 0.1 * dHphi;
0196
0197 if (trackphi > M_PI)
0198 trackphi -= 2 * M_PI;
0199 else if (trackphi < -M_PI)
0200 trackphi += 2 * M_PI;
0201
0202
0203
0204 double inflation = 1.;
0205 float rangeeta = (cristalPhiEtaMaxSize_ * (1.5 + 0.5) + 0.2 * fabs(dHeta)) * inflation;
0206 float rangephi = (cristalPhiEtaMaxSize_ * (1.5 + 0.5) + 0.2 * fabs(dHphi)) * inflation;
0207
0208
0209 std::vector<reco::PFRecHit const*> recHits;
0210 KDTreeBox trackBox(tracketa - rangeeta, tracketa + rangeeta, trackphi - rangephi, trackphi + rangephi);
0211 tree_.search(trackBox, recHits);
0212
0213
0214 for (auto const& recHit : recHits) {
0215 const auto& rhrep = recHit->positionREP();
0216 const auto& corners = recHit->getCornersREP();
0217
0218 double rhsizeeta = fabs(corners[3].eta() - corners[1].eta());
0219 double rhsizephi = fabs(corners[3].phi() - corners[1].phi());
0220 if (rhsizephi > M_PI)
0221 rhsizephi = 2. * M_PI - rhsizephi;
0222
0223 double deta = fabs(rhrep.eta() - tracketa);
0224 double dphi = fabs(rhrep.phi() - trackphi);
0225 if (dphi > M_PI)
0226 dphi = 2. * M_PI - dphi;
0227
0228
0229 RecHit2BlockEltMap::iterator ret = rechit2ClusterLinks_.find(recHit);
0230
0231 for (BlockEltSet::iterator clusterIt = ret->second.begin(); clusterIt != ret->second.end(); clusterIt++) {
0232 const reco::PFClusterRef clusterref = (*clusterIt)->clusterRef();
0233 int fracsNbr = clusterref->recHitFractions().size();
0234
0235 double _rhsizeeta = rhsizeeta * (1.5 + 0.5 / fracsNbr) + 0.2 * fabs(dHeta);
0236 double _rhsizephi = rhsizephi * (1.5 + 0.5 / fracsNbr) + 0.2 * fabs(dHphi);
0237
0238
0239 if (deta < (_rhsizeeta / 2.) && dphi < (_rhsizephi / 2.))
0240 target2ClusterLinks_[*it].insert(*clusterIt);
0241 }
0242 }
0243 }
0244 }
0245
0246 void KDTreeLinkerTrackHcal::updatePFBlockEltWithLinks() {
0247
0248
0249
0250 for (BlockElt2BlockEltMap::iterator it = target2ClusterLinks_.begin(); it != target2ClusterLinks_.end(); ++it) {
0251 const auto& trackElt = it->first;
0252 const auto& hcalEltSet = it->second;
0253 reco::PFMultiLinksTC multitracks(true);
0254
0255
0256
0257 if (nMaxHcalLinksPerTrack_ < 0. || trackElt->isLinkedToDisplacedVertex()) {
0258 for (const auto& hcalElt : hcalEltSet) {
0259 reco::PFMultilink multiLink(hcalElt->clusterRef());
0260 multitracks.linkedPFObjects.push_back(multiLink);
0261
0262
0263
0264 hcalElt->setIsValidMultilinks(true, _targetType);
0265 }
0266
0267 }
0268
0269
0270 else {
0271 const reco::PFRecTrackRef& trackref = trackElt->trackRefPF();
0272 const reco::PFTrajectoryPoint& tkAtHCALEnt = trackref->extrapolatedPoint(trajectoryLayerEntrance_);
0273 const reco::PFCluster::REPPoint& tkreppos = tkAtHCALEnt.positionREP();
0274
0275 double dHEta = 0.;
0276 double dHPhi = 0.;
0277 double dRHCALEx = 0.;
0278 if (checkExit_) {
0279 const reco::PFTrajectoryPoint& tkAtHCALEx = trackref->extrapolatedPoint(trajectoryLayerExit_);
0280 dHEta = (tkAtHCALEx.positionREP().Eta() - tkAtHCALEnt.positionREP().Eta());
0281 dHPhi = reco::deltaPhi(tkAtHCALEx.positionREP().Phi(), tkAtHCALEnt.positionREP().Phi());
0282 dRHCALEx = tkAtHCALEx.position().R();
0283 }
0284
0285 std::vector<double> vDist;
0286 double dist(-1.0);
0287
0288
0289 for (const auto& hcalElt : hcalEltSet) {
0290 double clusterphi = hcalElt->clusterRef()->positionREP().phi();
0291 double clustereta = hcalElt->clusterRef()->positionREP().eta();
0292
0293
0294 if (!checkExit_) {
0295 dist = LinkByRecHit::computeDist(clustereta, clusterphi, tkreppos.Eta(), tkreppos.Phi());
0296 }
0297
0298 else {
0299
0300
0301
0302 if (dRHCALEx < tkAtHCALEnt.position().R()) {
0303 dist = LinkByRecHit::computeDist(clustereta, clusterphi, tkreppos.Eta(), tkreppos.Phi());
0304 edm::LogWarning("KDTreeLinkerTrackHcal ")
0305 << "Special case of linking with track hitting HCAL and looping back in the tracker ";
0306 } else {
0307 dist = LinkByRecHit::computeDist(
0308 clustereta, clusterphi, tkreppos.Eta() + 0.1 * dHEta, tkreppos.Phi() + 0.1 * dHPhi);
0309 }
0310 }
0311
0312 vDist.push_back(dist);
0313 }
0314
0315
0316 for (auto i : sort_indexes(vDist)) {
0317 const BlockEltSet::iterator hcalEltIt = std::next(hcalEltSet.begin(), i);
0318 reco::PFMultilink multiLink((*hcalEltIt)->clusterRef());
0319 multitracks.linkedPFObjects.push_back(multiLink);
0320
0321
0322 (*hcalEltIt)->setIsValidMultilinks(true, _targetType);
0323
0324 if (multitracks.linkedPFObjects.size() >= (unsigned)nMaxHcalLinksPerTrack_)
0325 break;
0326 }
0327 }
0328
0329
0330 trackElt->setMultilinks(multitracks, _fieldType);
0331 }
0332 }
0333
0334 void KDTreeLinkerTrackHcal::clear() {
0335 targetSet_.clear();
0336 fieldClusterSet_.clear();
0337
0338 rechitsSet_.clear();
0339
0340 rechit2ClusterLinks_.clear();
0341 target2ClusterLinks_.clear();
0342
0343 tree_.clear();
0344 }