plugins/kdtrees/KDTreeLinkerTrackEcal.cc

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 // This class is used to find all links between Tracks and ECAL clusters
0009 // using a KDTree algorithm.
0010 // It is used in PFBlockAlgo.cc in the function links().
0011 class KDTreeLinkerTrackEcal : public KDTreeLinkerBase {
0012 public:
0013   KDTreeLinkerTrackEcal(const edm::ParameterSet &conf);
0014   ~KDTreeLinkerTrackEcal() override;
0015
0016   // With this method, we create the list of track that we want to link.
0017   void insertTargetElt(reco::PFBlockElement *track) override;
0018
0019   // Here, we create the list of ecalCluster that we want to link. From ecalCluster
0020   // and fraction, we will create a second list of rechits that will be used to
0021   // build the KDTree.
0022   void insertFieldClusterElt(reco::PFBlockElement *ecalCluster) override;
0023
0024   // The KDTree building from rechits list.
0025   void buildTree() override;
0026
0027   // Here we will iterate over all tracks. For each track intersection point with ECAL,
0028   // we will search the closest rechits in the KDTree, from rechits we will find the
0029   // ecalClusters and after that we will check the links between the track and
0030   // all closest ecalClusters.
0031   void searchLinks() override;
0032
0033   // Here, we will store all Track/ECAL founded links in the PFBlockElement class
0034   // of each psCluster in the PFmultilinks field.
0035   void updatePFBlockEltWithLinks() override;
0036
0037   // Here we free all allocated structures.
0038   void clear() override;
0039
0040 private:
0041   // Data used by the KDTree algorithm : sets of Tracks and ECAL clusters.
0042   BlockEltSet targetSet_;
0043   BlockEltSet fieldClusterSet_;
0044
0045   // Sets of rechits that compose the ECAL clusters.
0046   RecHitSet rechitsSet_;
0047
0048   // Map of linked Track/ECAL clusters.
0049   BlockElt2BlockEltMap cluster2TargetLinks_;
0050
0051   // Map of the ECAL clusters associated to a rechit.
0052   RecHit2BlockEltMap rechit2ClusterLinks_;
0053
0054   // KD trees
0055   KDTreeLinkerAlgo<reco::PFRecHit const *> tree_;
0056 };
0057
0058 // the text name is different so that we can easily
0059 // construct it when calling the factory
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   // This test is more or less done in PFBlockAlgo.h. In others cases, it should be switch on.
0076   //   if (!((clusterref->layer() == PFLayer::ECAL_ENDCAP) ||
0077   //    (clusterref->layer() == PFLayer::ECAL_BARREL)))
0078   //     return;
0079
0080   const std::vector<reco::PFRecHitFraction> &fraction = clusterref->recHitFractions();
0081
0082   // We create a list of ecalCluster
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     // We save the links rechit to EcalClusters
0094     rechit2ClusterLinks_[&rechit].insert(ecalCluster);
0095
0096     // We create a liste of rechits
0097     rechitsSet_.insert(&rechit);
0098   }
0099 }
0100
0101 void KDTreeLinkerTrackEcal::buildTree() {
0102   // List of pseudo-rechits that will be used to create the KDTree
0103   std::vector<KDTreeNodeInfo<reco::PFRecHit const *, 2>> eltList;
0104
0105   // Filling of this list
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     // Here we solve the problem of phi circular set by duplicating some rechits
0113     // too close to -Pi (or to Pi) and adding (substracting) to them 2 * Pi.
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   // Here we define the upper/lower bounds of the 2D space (eta/phi).
0128   float phimin = -1.0 * M_PI - phiOffset_;
0129   float phimax = M_PI + phiOffset_;
0130
0131   // etamin-etamax, phimin-phimax
0132   KDTreeBox region(-3.0f, 3.0f, phimin, phimax);
0133
0134   // We may now build the KDTree
0135   tree_.build(eltList, region);
0136 }
0137
0138 void KDTreeLinkerTrackEcal::searchLinks() {
0139   // Most of the code has been taken from LinkByRecHit.cc
0140
0141   // We iterate over the tracks.
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     // The track didn't reach ecal
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     // Estimate the maximal envelope in phi/eta that will be used to find rechit candidates.
0161     // Same envelope for cap et barrel rechits.
0162     float range = cristalPhiEtaMaxSize_ * (2.0 + 1.0 / std::min(1., trackPt / 2.));
0163
0164     // We search for all candidate recHits, ie all recHits contained in the maximal size envelope.
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     // Here we check all rechit candidates using the non-approximated method.
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       // Find all clusters associated to given rechit
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) {  // BARREL
0195           // Check if the track is in the barrel
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           // Check if the track and the cluster are linked
0203           if (deta < (_rhsizeeta / 2.) && dphi < (_rhsizephi / 2.))
0204             cluster2TargetLinks_[*clusterIt].insert(*it);
0205
0206         } else {  // ENDCAP
0207
0208           // Check if the track is in the cap
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           // Check if the track and the cluster are linked
0230           if (isinside)
0231             cluster2TargetLinks_[*clusterIt].insert(*it);
0232         }
0233       }
0234     }
0235   }
0236 }
0237
0238 void KDTreeLinkerTrackEcal::updatePFBlockEltWithLinks() {
0239   //TODO YG : Check if cluster positionREP() is valid ?
0240
0241   // Here we save in each ECAL cluster the list of phi/eta values of linked tracks.
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       // We set the multilinks flag of the track (for links to ECAL) to true. It will allow us to
0254       // use it in an optimized way in prefilter
0255       trackElt->setIsValidMultilinks(true, _fieldType);
0256     }
0257
0258     // We set multilinks of the ECAL element (for links to tracks)
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 }