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/PFProducer/interface/KDTreeLinkerBase.h"
0004 #include "CommonTools/RecoAlgos/interface/KDTreeLinkerAlgo.h"
0005
0006 #include "TMath.h"
0007
0008
0009
0010
0011 class KDTreeLinkerTrackEcal : public KDTreeLinkerBase {
0012 public:
0013 KDTreeLinkerTrackEcal(const edm::ParameterSet &conf);
0014 ~KDTreeLinkerTrackEcal() override;
0015
0016
0017 void insertTargetElt(reco::PFBlockElement *track) override;
0018
0019
0020
0021
0022 void insertFieldClusterElt(reco::PFBlockElement *ecalCluster) 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 cluster2TargetLinks_;
0050
0051
0052 RecHit2BlockEltMap rechit2ClusterLinks_;
0053
0054
0055 KDTreeLinkerAlgo<reco::PFRecHit const *> tree_;
0056 };
0057
0058
0059
0060 DEFINE_EDM_PLUGIN(KDTreeLinkerFactory, KDTreeLinkerTrackEcal, "KDTreeTrackAndECALLinker");
0061
0062 KDTreeLinkerTrackEcal::KDTreeLinkerTrackEcal(const edm::ParameterSet &conf) : KDTreeLinkerBase(conf) {}
0063
0064 KDTreeLinkerTrackEcal::~KDTreeLinkerTrackEcal() { clear(); }
0065
0066 void KDTreeLinkerTrackEcal::insertTargetElt(reco::PFBlockElement *track) {
0067 if (track->trackRefPF()->extrapolatedPoint(reco::PFTrajectoryPoint::ECALShowerMax).isValid()) {
0068 targetSet_.insert(track);
0069 }
0070 }
0071
0072 void KDTreeLinkerTrackEcal::insertFieldClusterElt(reco::PFBlockElement *ecalCluster) {
0073 const reco::PFClusterRef &clusterref = ecalCluster->clusterRef();
0074
0075
0076
0077
0078
0079
0080 const std::vector<reco::PFRecHitFraction> &fraction = clusterref->recHitFractions();
0081
0082
0083 fieldClusterSet_.insert(ecalCluster);
0084 for (size_t rhit = 0; rhit < fraction.size(); ++rhit) {
0085 const reco::PFRecHitRef &rh = fraction[rhit].recHitRef();
0086 double fract = fraction[rhit].fraction();
0087
0088 if ((rh.isNull()) || (fract < cutOffFrac))
0089 continue;
0090
0091 const reco::PFRecHit &rechit = *rh;
0092
0093
0094 rechit2ClusterLinks_[&rechit].insert(ecalCluster);
0095
0096
0097 rechitsSet_.insert(&rechit);
0098 }
0099 }
0100
0101 void KDTreeLinkerTrackEcal::buildTree() {
0102
0103 std::vector<KDTreeNodeInfo<reco::PFRecHit const *, 2>> eltList;
0104
0105
0106 for (RecHitSet::const_iterator it = rechitsSet_.begin(); it != rechitsSet_.end(); it++) {
0107 const reco::PFRecHit::REPPoint &posrep = (*it)->positionREP();
0108
0109 KDTreeNodeInfo<reco::PFRecHit const *, 2> rh1(*it, posrep.eta(), posrep.phi());
0110 eltList.push_back(rh1);
0111
0112
0113
0114 if (rh1.dims[1] > (M_PI - phiOffset_)) {
0115 float phi = rh1.dims[1] - 2 * M_PI;
0116 KDTreeNodeInfo<reco::PFRecHit const *, 2> rh2(*it, float(posrep.eta()), phi);
0117 eltList.push_back(rh2);
0118 }
0119
0120 if (rh1.dims[1] < (M_PI * -1.0 + phiOffset_)) {
0121 float phi = rh1.dims[1] + 2 * M_PI;
0122 KDTreeNodeInfo<reco::PFRecHit const *, 2> rh3(*it, float(posrep.eta()), phi);
0123 eltList.push_back(rh3);
0124 }
0125 }
0126
0127
0128 float phimin = -1.0 * M_PI - phiOffset_;
0129 float phimax = M_PI + phiOffset_;
0130
0131
0132 KDTreeBox region(-3.0f, 3.0f, phimin, phimax);
0133
0134
0135 tree_.build(eltList, region);
0136 }
0137
0138 void KDTreeLinkerTrackEcal::searchLinks() {
0139
0140
0141
0142 for (BlockEltSet::iterator it = targetSet_.begin(); it != targetSet_.end(); it++) {
0143 reco::PFRecTrackRef trackref = (*it)->trackRefPF();
0144
0145 const reco::PFTrajectoryPoint &atECAL = trackref->extrapolatedPoint(reco::PFTrajectoryPoint::ECALShowerMax);
0146
0147
0148 if (!atECAL.isValid())
0149 continue;
0150
0151 const reco::PFTrajectoryPoint &atVertex = trackref->extrapolatedPoint(reco::PFTrajectoryPoint::ClosestApproach);
0152
0153 double trackPt = sqrt(atVertex.momentum().Vect().Perp2());
0154 float tracketa = atECAL.positionREP().eta();
0155 float trackphi = atECAL.positionREP().phi();
0156 double trackx = atECAL.position().X();
0157 double tracky = atECAL.position().Y();
0158 double trackz = atECAL.position().Z();
0159
0160
0161
0162 float range = cristalPhiEtaMaxSize_ * (2.0 + 1.0 / std::min(1., trackPt / 2.));
0163
0164
0165 std::vector<reco::PFRecHit const *> recHits;
0166 KDTreeBox trackBox(tracketa - range, tracketa + range, trackphi - range, trackphi + range);
0167 tree_.search(trackBox, recHits);
0168
0169
0170 for (auto const &recHit : recHits) {
0171 const auto &cornersxyz = recHit->getCornersXYZ();
0172 const auto &posxyz = recHit->position();
0173 const auto &rhrep = recHit->positionREP();
0174 const auto &corners = recHit->getCornersREP();
0175
0176 double rhsizeeta = fabs(corners[3].eta() - corners[1].eta());
0177 double rhsizephi = fabs(corners[3].phi() - corners[1].phi());
0178 if (rhsizephi > M_PI)
0179 rhsizephi = 2. * M_PI - rhsizephi;
0180
0181 double deta = fabs(rhrep.eta() - tracketa);
0182 double dphi = fabs(rhrep.phi() - trackphi);
0183 if (dphi > M_PI)
0184 dphi = 2. * M_PI - dphi;
0185
0186
0187 RecHit2BlockEltMap::iterator ret = rechit2ClusterLinks_.find(recHit);
0188
0189 for (BlockEltSet::const_iterator clusterIt = ret->second.begin(); clusterIt != ret->second.end(); clusterIt++) {
0190 reco::PFClusterRef clusterref = (*clusterIt)->clusterRef();
0191 double clusterz = clusterref->position().z();
0192 int fracsNbr = clusterref->recHitFractions().size();
0193
0194 if (clusterref->layer() == PFLayer::ECAL_BARREL) {
0195
0196 if (fabs(trackz) > 300.)
0197 continue;
0198
0199 double _rhsizeeta = rhsizeeta * (2.00 + 1.0 / (fracsNbr * std::min(1., trackPt / 2.)));
0200 double _rhsizephi = rhsizephi * (2.00 + 1.0 / (fracsNbr * std::min(1., trackPt / 2.)));
0201
0202
0203 if (deta < (_rhsizeeta / 2.) && dphi < (_rhsizephi / 2.))
0204 cluster2TargetLinks_[*clusterIt].insert(*it);
0205
0206 } else {
0207
0208
0209 if (fabs(trackz) < 300.)
0210 continue;
0211 if (trackz * clusterz < 0.)
0212 continue;
0213
0214 double x[5];
0215 double y[5];
0216 for (unsigned jc = 0; jc < 4; ++jc) {
0217 auto cornerposxyz = cornersxyz[jc];
0218 x[3 - jc] = cornerposxyz.x() +
0219 (cornerposxyz.x() - posxyz.x()) * (1.00 + 0.50 / fracsNbr / std::min(1., trackPt / 2.));
0220 y[3 - jc] = cornerposxyz.y() +
0221 (cornerposxyz.y() - posxyz.y()) * (1.00 + 0.50 / fracsNbr / std::min(1., trackPt / 2.));
0222 }
0223
0224 x[4] = x[0];
0225 y[4] = y[0];
0226
0227 bool isinside = TMath::IsInside(trackx, tracky, 5, x, y);
0228
0229
0230 if (isinside)
0231 cluster2TargetLinks_[*clusterIt].insert(*it);
0232 }
0233 }
0234 }
0235 }
0236 }
0237
0238 void KDTreeLinkerTrackEcal::updatePFBlockEltWithLinks() {
0239
0240
0241
0242 for (BlockElt2BlockEltMap::iterator it = cluster2TargetLinks_.begin(); it != cluster2TargetLinks_.end(); ++it) {
0243 const auto &ecalElt = it->first;
0244 const auto &trackEltSet = it->second;
0245 reco::PFMultiLinksTC multitracks(true);
0246
0247 for (const auto &trackElt : trackEltSet) {
0248 const reco::PFRecTrackRef &trackref = trackElt->trackRefPF();
0249
0250 reco::PFMultilink multiLink(trackref);
0251 multitracks.linkedPFObjects.push_back(multiLink);
0252
0253
0254
0255 trackElt->setIsValidMultilinks(true, _fieldType);
0256 }
0257
0258
0259 ecalElt->setMultilinks(multitracks, _targetType);
0260 }
0261 }
0262
0263 void KDTreeLinkerTrackEcal::clear() {
0264 targetSet_.clear();
0265 fieldClusterSet_.clear();
0266
0267 rechitsSet_.clear();
0268
0269 rechit2ClusterLinks_.clear();
0270 cluster2TargetLinks_.clear();
0271
0272 tree_.clear();
0273 }