plugins/kdtrees/KDTreeLinkerTrackHcal.cc

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 // This class is used to find all links between Tracks and HCAL clusters
0009 // using a KDTree algorithm.
0010 // It is used in PFBlockAlgo.cc in the function links().
0011 class KDTreeLinkerTrackHcal : public KDTreeLinkerBase {
0012 public:
0013   KDTreeLinkerTrackHcal(const edm::ParameterSet& conf);
0014   ~KDTreeLinkerTrackHcal() 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 hcalCluster that we want to link. From hcalCluster
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* hcalCluster) 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 HCAL,
0028   // we will search the closest rechits in the KDTree, from rechits we will find the
0029   // hcalClusters and after that we will check the links between the track and
0030   // all closest hcalClusters.
0031   void searchLinks() override;
0032
0033   // Here, we will store all Track/HCAL 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 HCAL clusters.
0042   BlockEltSet targetSet_;
0043   BlockEltSet fieldClusterSet_;
0044
0045   // Sets of rechits that compose the HCAL clusters.
0046   RecHitSet rechitsSet_;
0047
0048   // Map of linked Track/HCAL clusters.
0049   BlockElt2BlockEltMap target2ClusterLinks_;
0050
0051   // Map of the HCAL clusters associated to a rechit.
0052   RecHit2BlockEltMap rechit2ClusterLinks_;
0053
0054   // KD trees
0055   KDTreeLinkerAlgo<reco::PFRecHit const*> tree_;
0056
0057   // TrajectoryPoints
0058   std::string trajectoryLayerEntranceString_;
0059   std::string trajectoryLayerExitString_;
0060   reco::PFTrajectoryPoint::LayerType trajectoryLayerEntrance_;
0061   reco::PFTrajectoryPoint::LayerType trajectoryLayerExit_;
0062   bool checkExit_;
0063
0064   // Hcal-track links
0065   int nMaxHcalLinksPerTrack_;
0066 };
0067
0068 // the text name is different so that we can easily
0069 // construct it when calling the factory
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   // Initialization
0078   cristalPhiEtaMaxSize_ = 0.2;
0079   phiOffset_ = 0.32;
0080   // convert TrajectoryLayers info from string to enum
0081   trajectoryLayerEntrance_ = reco::PFTrajectoryPoint::layerTypeByName(trajectoryLayerEntranceString_);
0082   trajectoryLayerExit_ = reco::PFTrajectoryPoint::layerTypeByName(trajectoryLayerExitString_);
0083   // make sure the requested setting is supported
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   // flag if exit layer should be checked or not
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   // This test is more or less done in PFBlockAlgo.h. In others cases, it should be switch on.
0106   //   if (!((clusterref->layer() == PFLayer::HCAL_ENDCAP) ||
0107   //    (clusterref->layer() == PFLayer::HCAL_BARREL1)))
0108   //     return;
0109
0110   const std::vector<reco::PFRecHitFraction>& fraction = clusterref->recHitFractions();
0111
0112   // We create a list of hcalCluster
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     // We save the links rechit to HcalClusters
0124     rechit2ClusterLinks_[&rechit].insert(hcalCluster);
0125
0126     // We create a liste of rechits
0127     rechitsSet_.insert(&rechit);
0128   }
0129 }
0130
0131 void KDTreeLinkerTrackHcal::buildTree() {
0132   // List of pseudo-rechits that will be used to create the KDTree
0133   std::vector<KDTreeNodeInfo<reco::PFRecHit const*, 2>> eltList;
0134
0135   // Filling of this list
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     // Here we solve the problem of phi circular set by duplicating some rechits
0143     // too close to -Pi (or to Pi) and adding (substracting) to them 2 * Pi.
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   // Here we define the upper/lower bounds of the 2D space (eta/phi).
0158   float phimin = -1.0 * M_PI - phiOffset_;
0159   float phimax = M_PI + phiOffset_;
0160
0161   // etamin-etamax, phimin-phimax
0162   KDTreeBox region(-3.0f, 3.0f, phimin, phimax);
0163
0164   // We may now build the KDTree
0165   tree_.build(eltList, region);
0166 }
0167
0168 void KDTreeLinkerTrackHcal::searchLinks() {
0169   // Most of the code has been taken from LinkByRecHit.cc
0170
0171   // We iterate over the tracks.
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     // The track didn't reach hcal
0178     if (!atHCAL.isValid())
0179       continue;
0180
0181     // In case the exit point check is requested, check eta and phi differences between entrance and exit
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     }  // checkExit_
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     // Estimate the maximal envelope in phi/eta that will be used to find rechit candidates.
0203     // Same envelope for cap et barrel rechits.
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     // We search for all candidate recHits, ie all recHits contained in the maximal size envelope.
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     // Here we check all rechit candidates using the non-approximated method.
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       // Find all clusters associated to given rechit
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         // Check if the track and the cluster are linked
0239         if (deta < (_rhsizeeta / 2.) && dphi < (_rhsizephi / 2.))
0240           target2ClusterLinks_[*it].insert(*clusterIt);
0241       }
0242     }
0243   }
0244 }
0245
0246 void KDTreeLinkerTrackHcal::updatePFBlockEltWithLinks() {
0247   //TODO YG : Check if cluster positionREP() is valid ?
0248
0249   // Here we save in each track the list of phi/eta values of linked clusters.
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     // No restriction on the number of HCAL links per track or isLinkedToDisplacedVertex
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         // We set the multilinks flag of the track (for links to ECAL) to true. It will allow us to
0263         // use it in an optimized way in prefilter
0264         hcalElt->setIsValidMultilinks(true, _targetType);
0265       }
0266
0267     }
0268     //
0269     // Store only the N closest HCAL links per track.
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       // Check exit point
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       // Fill the vector of distances between HCAL clusters and the track
0289       for (const auto& hcalElt : hcalEltSet) {
0290         double clusterphi = hcalElt->clusterRef()->positionREP().phi();
0291         double clustereta = hcalElt->clusterRef()->positionREP().eta();
0292
0293         // when checkExit_ is false
0294         if (!checkExit_) {
0295           dist = LinkByRecHit::computeDist(clustereta, clusterphi, tkreppos.Eta(), tkreppos.Phi());
0296         }
0297         // when checkExit_ is true
0298         else {
0299           //special case ! A looper  can exit the barrel inwards and hit the endcap
0300           //In this case calculate the distance based on the first crossing since
0301           //the looper will probably never make it to the endcap
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         }  // checkExit_
0311
0312         vDist.push_back(dist);
0313       }  // loop over hcalEltSet
0314
0315       // Fill multitracks
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         // We set the multilinks flag of the track (for links to ECAL) to true. It will allow us to
0321         // use it in an optimized way in prefilter
0322         (*hcalEltIt)->setIsValidMultilinks(true, _targetType);
0323
0324         if (multitracks.linkedPFObjects.size() >= (unsigned)nMaxHcalLinksPerTrack_)
0325           break;
0326       }
0327     }
0328
0329     // Store multitracks
0330     trackElt->setMultilinks(multitracks, _fieldType);
0331   }  // loop over target2ClusterLinks_
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 }