File indexing completed on 2024-04-06 12:27:28
0001 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0002 #include "RecoParticleFlow/PFProducer/interface/KDTreeLinkerBase.h"
0003 #include "CommonTools/RecoAlgos/interface/KDTreeLinkerAlgo.h"
0004
0005 #include "TMath.h"
0006
0007
0008
0009
0010 class KDTreeLinkerPSEcal : public KDTreeLinkerBase {
0011 public:
0012 KDTreeLinkerPSEcal(const edm::ParameterSet &conf);
0013 ~KDTreeLinkerPSEcal() override;
0014
0015
0016 void insertTargetElt(reco::PFBlockElement *psCluster) override;
0017
0018
0019
0020
0021 void insertFieldClusterElt(reco::PFBlockElement *ecalCluster) override;
0022
0023
0024 void buildTree() override;
0025
0026
0027
0028
0029 void searchLinks() override;
0030
0031
0032
0033 void updatePFBlockEltWithLinks() override;
0034
0035
0036 void clear() override;
0037
0038 private:
0039
0040 void buildTree(const RecHitSet &rechitsSet, KDTreeLinkerAlgo<reco::PFRecHit const *> &tree);
0041
0042 private:
0043
0044 const double resPSpitch_;
0045 const double resPSlength_;
0046 const double ps1ToEcal_;
0047 const double ps2ToEcal_;
0048
0049
0050 BlockEltSet targetSet_;
0051 BlockEltSet fieldClusterSet_;
0052
0053
0054
0055 RecHitSet rechitsNegSet_;
0056 RecHitSet rechitsPosSet_;
0057
0058
0059 BlockElt2BlockEltMap target2ClusterLinks_;
0060
0061
0062 RecHit2BlockEltMap rechit2ClusterLinks_;
0063
0064
0065 KDTreeLinkerAlgo<reco::PFRecHit const *> treeNeg_;
0066 KDTreeLinkerAlgo<reco::PFRecHit const *> treePos_;
0067 };
0068
0069
0070
0071 DEFINE_EDM_PLUGIN(KDTreeLinkerFactory, KDTreeLinkerPSEcal, "KDTreePreshowerAndECALLinker");
0072
0073 KDTreeLinkerPSEcal::KDTreeLinkerPSEcal(const edm::ParameterSet &conf)
0074 : KDTreeLinkerBase(conf), resPSpitch_(0.19), resPSlength_(6.1), ps1ToEcal_(1.072), ps2ToEcal_(1.057) {}
0075
0076 KDTreeLinkerPSEcal::~KDTreeLinkerPSEcal() { clear(); }
0077
0078 void KDTreeLinkerPSEcal::insertTargetElt(reco::PFBlockElement *psCluster) {
0079
0080
0081
0082 targetSet_.insert(psCluster);
0083 }
0084
0085 void KDTreeLinkerPSEcal::insertFieldClusterElt(reco::PFBlockElement *ecalCluster) {
0086 const reco::PFClusterRef &clusterref = ecalCluster->clusterRef();
0087
0088 if (clusterref->layer() != PFLayer::ECAL_ENDCAP)
0089 return;
0090
0091 const std::vector<reco::PFRecHitFraction> &fraction = clusterref->recHitFractions();
0092
0093
0094 fieldClusterSet_.insert(ecalCluster);
0095
0096 double clusterz = clusterref->position().Z();
0097 RecHitSet &rechitsSet = (clusterz < 0) ? rechitsNegSet_ : rechitsPosSet_;
0098
0099 for (size_t rhit = 0; rhit < fraction.size(); ++rhit) {
0100 const reco::PFRecHitRef &rh = fraction[rhit].recHitRef();
0101 double fract = fraction[rhit].fraction();
0102
0103 if ((rh.isNull()) || (fract < cutOffFrac))
0104 continue;
0105
0106 const reco::PFRecHit &rechit = *rh;
0107
0108
0109 rechit2ClusterLinks_[&rechit].insert(ecalCluster);
0110
0111
0112 rechitsSet.insert(&rechit);
0113 }
0114 }
0115
0116 void KDTreeLinkerPSEcal::buildTree() {
0117 buildTree(rechitsNegSet_, treeNeg_);
0118 buildTree(rechitsPosSet_, treePos_);
0119 }
0120
0121 void KDTreeLinkerPSEcal::buildTree(const RecHitSet &rechitsSet, KDTreeLinkerAlgo<reco::PFRecHit const *> &tree) {
0122
0123 std::vector<KDTreeNodeInfo<reco::PFRecHit const *, 2>> eltList;
0124
0125
0126 for (RecHitSet::const_iterator it = rechitsSet.begin(); it != rechitsSet.end(); it++) {
0127 const reco::PFRecHit *rh = *it;
0128 const auto &posxyz = rh->position();
0129
0130 KDTreeNodeInfo<reco::PFRecHit const *, 2> rhinfo{rh, posxyz.x(), posxyz.y()};
0131 eltList.push_back(rhinfo);
0132 }
0133
0134
0135 KDTreeBox region{-150.f, 150.f, -150.f, 150.f};
0136
0137
0138 tree.build(eltList, region);
0139 }
0140
0141 void KDTreeLinkerPSEcal::searchLinks() {
0142
0143
0144
0145 for (BlockEltSet::iterator it = targetSet_.begin(); it != targetSet_.end(); it++) {
0146 reco::PFClusterRef clusterPSRef = (*it)->clusterRef();
0147 const reco::PFCluster &clusterPS = *clusterPSRef;
0148
0149
0150 double zPS = clusterPS.position().Z();
0151 double xPS = clusterPS.position().X();
0152 double yPS = clusterPS.position().Y();
0153
0154 double etaPS = fabs(clusterPS.positionREP().eta());
0155 double deltaX = 0.;
0156 double deltaY = 0.;
0157 float xPSonEcal = xPS;
0158 float yPSonEcal = yPS;
0159
0160 if (clusterPS.layer() == PFLayer::PS1) {
0161
0162
0163 deltaX = resPSpitch_;
0164 deltaY = resPSlength_;
0165 xPSonEcal *= ps1ToEcal_;
0166 yPSonEcal *= ps1ToEcal_;
0167
0168 } else {
0169
0170
0171 deltaY = resPSpitch_;
0172 deltaX = resPSlength_;
0173 xPSonEcal *= ps2ToEcal_;
0174 yPSonEcal *= ps2ToEcal_;
0175 }
0176
0177
0178
0179
0180 double maxEcalRadius = cristalXYMaxSize_ / 2.;
0181
0182
0183 double inflation = 2.4 - (etaPS - 1.6);
0184 float rangeX = maxEcalRadius * (1 + (0.05 + 1.0 / maxEcalRadius * deltaX / 2.)) * inflation;
0185 float rangeY = maxEcalRadius * (1 + (0.05 + 1.0 / maxEcalRadius * deltaY / 2.)) * inflation;
0186
0187
0188 std::vector<reco::PFRecHit const *> recHits;
0189 KDTreeBox trackBox(xPSonEcal - rangeX, xPSonEcal + rangeX, yPSonEcal - rangeY, yPSonEcal + rangeY);
0190
0191 if (zPS < 0)
0192 treeNeg_.search(trackBox, recHits);
0193 else
0194 treePos_.search(trackBox, recHits);
0195
0196 for (auto const &recHit : recHits) {
0197 const auto &corners = recHit->getCornersXYZ();
0198
0199
0200 RecHit2BlockEltMap::iterator ret = rechit2ClusterLinks_.find(recHit);
0201
0202 for (BlockEltSet::const_iterator clusterIt = ret->second.begin(); clusterIt != ret->second.end(); clusterIt++) {
0203 reco::PFClusterRef clusterref = (*clusterIt)->clusterRef();
0204 double clusterz = clusterref->position().z();
0205
0206 const auto &posxyz = recHit->position() * zPS / clusterz;
0207
0208 double x[5];
0209 double y[5];
0210 for (unsigned jc = 0; jc < 4; ++jc) {
0211 auto cornerpos = corners[jc].basicVector() * zPS / clusterz;
0212 x[3 - jc] = cornerpos.x() +
0213 (cornerpos.x() - posxyz.x()) * (0.05 + 1.0 / fabs((cornerpos.x() - posxyz.x())) * deltaX / 2.);
0214 y[3 - jc] = cornerpos.y() +
0215 (cornerpos.y() - posxyz.y()) * (0.05 + 1.0 / fabs((cornerpos.y() - posxyz.y())) * deltaY / 2.);
0216 }
0217
0218 x[4] = x[0];
0219 y[4] = y[0];
0220
0221 bool isinside = TMath::IsInside(xPS, yPS, 5, x, y);
0222
0223
0224 if (isinside)
0225 target2ClusterLinks_[*it].insert(*clusterIt);
0226 }
0227 }
0228 }
0229 }
0230
0231 void KDTreeLinkerPSEcal::updatePFBlockEltWithLinks() {
0232
0233
0234
0235 for (BlockElt2BlockEltMap::iterator it = target2ClusterLinks_.begin(); it != target2ClusterLinks_.end(); ++it) {
0236 const auto &psElt = it->first;
0237 const auto &ecalEltSet = it->second;
0238 reco::PFMultiLinksTC multitracks(true);
0239
0240 for (const auto &ecalElt : ecalEltSet) {
0241 reco::PFMultilink multiLink(ecalElt->clusterRef());
0242 multitracks.linkedPFObjects.push_back(multiLink);
0243
0244
0245
0246 ecalElt->setIsValidMultilinks(true, _targetType);
0247 }
0248
0249
0250 psElt->setMultilinks(multitracks, _fieldType);
0251 }
0252 }
0253
0254 void KDTreeLinkerPSEcal::clear() {
0255 targetSet_.clear();
0256 fieldClusterSet_.clear();
0257
0258 rechitsNegSet_.clear();
0259 rechitsPosSet_.clear();
0260
0261 rechit2ClusterLinks_.clear();
0262 target2ClusterLinks_.clear();
0263
0264 treeNeg_.clear();
0265 treePos_.clear();
0266 }