plugins/kdtrees/KDTreeLinkerPSEcal.cc

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 // This class is used to find all links between PreShower clusters and ECAL clusters
0008 // using a KDTree algorithm.
0009 // It is used in PFBlockAlgo.cc in the function links().
0010 class KDTreeLinkerPSEcal : public KDTreeLinkerBase {
0011 public:
0012   KDTreeLinkerPSEcal(const edm::ParameterSet &conf);
0013   ~KDTreeLinkerPSEcal() override;
0014
0015   // With this method, we create the list of psCluster that we want to link.
0016   void insertTargetElt(reco::PFBlockElement *psCluster) override;
0017
0018   // Here, we create the list of ecalCluster that we want to link. From ecalCluster
0019   // and fraction, we will create a second list of rechits that will be used to
0020   // build the KDTree.
0021   void insertFieldClusterElt(reco::PFBlockElement *ecalCluster) override;
0022
0023   // The KDTree building from rechits list.
0024   void buildTree() override;
0025
0026   // Here we will iterate over all psCluster. For each one, we will search the closest
0027   // rechits in the KDTree, from rechits we will find the ecalClusters and after that
0028   // we will check the links between the psCluster and all closest ecalClusters.
0029   void searchLinks() override;
0030
0031   // Here, we will store all PS/ECAL founded links in the PFBlockElement class
0032   // of each psCluster in the PFmultilinks field.
0033   void updatePFBlockEltWithLinks() override;
0034
0035   // Here we free all allocated structures.
0036   void clear() override;
0037
0038 private:
0039   // This method allows us to build the "tree" from the "rechitsSet".
0040   void buildTree(const RecHitSet &rechitsSet, KDTreeLinkerAlgo<reco::PFRecHit const *> &tree);
0041
0042 private:
0043   // Some const values.
0044   const double resPSpitch_;
0045   const double resPSlength_;
0046   const double ps1ToEcal_;  // ratio : zEcal / zPS1
0047   const double ps2ToEcal_;  // ration : zEcal / zPS2
0048
0049   // Data used by the KDTree algorithm : sets of PS and ECAL clusters.
0050   BlockEltSet targetSet_;
0051   BlockEltSet fieldClusterSet_;
0052
0053   // Sets of rechits that compose the ECAL clusters. We differenctiate
0054   // the rechits by their Z value.
0055   RecHitSet rechitsNegSet_;
0056   RecHitSet rechitsPosSet_;
0057
0058   // Map of linked PS/ECAL clusters.
0059   BlockElt2BlockEltMap target2ClusterLinks_;
0060
0061   // Map of the ECAL clusters associated to a rechit.
0062   RecHit2BlockEltMap rechit2ClusterLinks_;
0063
0064   // KD trees
0065   KDTreeLinkerAlgo<reco::PFRecHit const *> treeNeg_;
0066   KDTreeLinkerAlgo<reco::PFRecHit const *> treePos_;
0067 };
0068
0069 // the text name is different so that we can easily
0070 // construct it when calling the factory
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   // This test is more or less done in PFBlockAlgo.h. In others cases, it should be switch on.
0080   //if (!((psCluster->clusterRef()->layer() == PFLayer::PS1) || (psCluster->clusterRef()->layer() == PFLayer::PS2)))
0081   //  return;
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   // We create a list of cluster
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     // We save the links rechit to Clusters
0109     rechit2ClusterLinks_[&rechit].insert(ecalCluster);
0110
0111     // We create a liste of rechits
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   // List of pseudo-rechits that will be used to create the KDTree
0123   std::vector<KDTreeNodeInfo<reco::PFRecHit const *, 2>> eltList;
0124
0125   // Filling of this eltList
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   // xmin-xmax, ymain-ymax
0135   KDTreeBox region{-150.f, 150.f, -150.f, 150.f};
0136
0137   // We may now build the KDTree
0138   tree.build(eltList, region);
0139 }
0140
0141 void KDTreeLinkerPSEcal::searchLinks() {
0142   // Most of the code has been taken from LinkByRecHit.cc
0143
0144   // We iterate over the PS clusters.
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     // PS cluster position, extrapolated to ECAL
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) {  // PS1
0161
0162       // vertical strips, measure x with pitch precision
0163       deltaX = resPSpitch_;
0164       deltaY = resPSlength_;
0165       xPSonEcal *= ps1ToEcal_;
0166       yPSonEcal *= ps1ToEcal_;
0167
0168     } else {  // PS2
0169
0170       // horizontal strips, measure y with pitch precision
0171       deltaY = resPSpitch_;
0172       deltaX = resPSlength_;
0173       xPSonEcal *= ps2ToEcal_;
0174       yPSonEcal *= ps2ToEcal_;
0175     }
0176
0177     // Estimate the maximal envelope in phi/eta that will be used to find rechit candidates.
0178     // Same envelope for cap et barrel rechits.
0179
0180     double maxEcalRadius = cristalXYMaxSize_ / 2.;
0181
0182     // The inflation factor includes the approximate projection from Preshower to ECAL
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     // We search for all candidate recHits, ie all recHits contained in the maximal size envelope.
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       // Find all clusters associated to given rechit
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         // Check if the track and the cluster are linked
0224         if (isinside)
0225           target2ClusterLinks_[*it].insert(*clusterIt);
0226       }
0227     }
0228   }
0229 }
0230
0231 void KDTreeLinkerPSEcal::updatePFBlockEltWithLinks() {
0232   //TODO YG : Check if cluster positionREP() is valid ?
0233
0234   // Here we save in each PS the list of phi/eta values of linked ECAL clusters.
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       // We set the multilinks flag of the ECAL element (for links to PS) to true. It will allow us to
0245       // use it in an optimized way in prefilter
0246       ecalElt->setIsValidMultilinks(true, _targetType);
0247     }
0248
0249     // We set multilinks of the PS element (for links to ECAL)
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 }