HiTracking/src/HICaloCompatibleTrackSelector.cc

0001 /*
0002
0003  Based on analytical track selector
0004
0005  - This track selector assigns a quality bit to tracks deemed compatible with matching calo info
0006  - The default mode is to use the matching provided by particle flow,
0007    but a delta R matching to calorimeter towers is also supported
0008  - No selection is done other then selecting calo-compatible tracks.
0009  - The code always keeps all tracks in the input collection (should make configurable)
0010  - Note that matching by PF candidate only works on the same track collection used as input to PF
0011  - Tower code not up to data
0012
0013    Authors:  Matthew Nguyen, Andre Yoon, Frank Ma (November 4th, 2011)
0014
0015  */
0016
0017 // Basic inclusion
0018 #include "RecoHI/HiTracking/interface/HICaloCompatibleTrackSelector.h"
0019
0020 #include "FWCore/Framework/interface/EventSetup.h"
0021
0022 #include "DataFormats/Math/interface/deltaR.h"
0023 #include <Math/DistFunc.h>
0024 #include "TMath.h"
0025
0026 using reco::modules::HICaloCompatibleTrackSelector;
0027
0028 HICaloCompatibleTrackSelector::HICaloCompatibleTrackSelector(const edm::ParameterSet& cfg)
0029     : srcTracks_(consumes<TrackCollection>(cfg.getParameter<edm::InputTag>("srcTracks"))),
0030       srcPFCands_(consumes<PFCandidateCollection>(cfg.getParameter<edm::InputTag>("srcPFCands"))),
0031       srcTower_(consumes<CaloTowerCollection>(cfg.getParameter<edm::InputTag>("srcTower"))),
0032       usePFCandMatching_(cfg.getUntrackedParameter<bool>("usePFCandMatching", true)),
0033       trkMatchPtMin_(cfg.getUntrackedParameter<double>("trkMatchPtMin", 10.0)),
0034       trkCompPtMin_(cfg.getUntrackedParameter<double>("trkCompPtMin", 35.0)),
0035       trkEtaMax_(cfg.getUntrackedParameter<double>("trkEtaMax", 2.4)),
0036       towerPtMin_(cfg.getUntrackedParameter<double>("towerPtMin", 5.0)),
0037       matchConeRadius_(cfg.getUntrackedParameter<double>("matchConeRadius", 0.087)),
0038       keepAllTracks_(cfg.getUntrackedParameter<bool>("keepAllTracks", true)),
0039       copyExtras_(cfg.getUntrackedParameter<bool>("copyExtras", true)),
0040       copyTrajectories_(cfg.getUntrackedParameter<bool>("copyTrajectories", true)),
0041       qualityToSet_(cfg.getParameter<std::string>("qualityToSet")),
0042       qualityToSkip_(cfg.getParameter<std::string>("qualityToSkip")),
0043       qualityToMatch_(cfg.getParameter<std::string>("qualityToMatch")),
0044       minimumQuality_(cfg.getParameter<std::string>("minimumQuality")),
0045       resetQuality_(cfg.getUntrackedParameter<bool>("resetQuality", true)),
0046       passMuons_(cfg.getUntrackedParameter<bool>("passMuons", true)),
0047       passElectrons_(cfg.getUntrackedParameter<bool>("passElectrons", false)),
0048       funcDeltaRTowerMatch_(cfg.getParameter<std::string>("funcDeltaRTowerMatch")),
0049       funcCaloComp_(cfg.getParameter<std::string>("funcCaloComp")) {
0050   std::string alias(cfg.getParameter<std::string>("@module_label"));
0051   produces<reco::TrackCollection>().setBranchAlias(alias + "Tracks");
0052   if (copyExtras_) {
0053     produces<reco::TrackExtraCollection>().setBranchAlias(alias + "TrackExtras");
0054     produces<TrackingRecHitCollection>().setBranchAlias(alias + "RecHits");
0055   }
0056   if (copyTrajectories_) {
0057     produces<std::vector<Trajectory>>().setBranchAlias(alias + "Trajectories");
0058     produces<TrajTrackAssociationCollection>().setBranchAlias(alias + "TrajectoryTrackAssociations");
0059     srcTrackTrajs_ = (consumes<std::vector<Trajectory>>(cfg.getParameter<edm::InputTag>("srcTracks")));
0060     srcTrackTrajAssoc_ = (consumes<TrajTrackAssociationCollection>(cfg.getParameter<edm::InputTag>("srcTracks")));
0061   }
0062
0063   // pt dependence of delta R matching requirement
0064   fDeltaRTowerMatch = new TF1("fDeltaRTowerMatch", funcDeltaRTowerMatch_.c_str(), 0, 200);
0065   // pt dependance of calo compatibility, i.e., minimum sum Calo Et vs. track pT
0066   fCaloComp = new TF1("fCaloComp", funcCaloComp_.c_str(), 0, 200);  // a parameterization of pt dependent cut
0067 }
0068
0069 HICaloCompatibleTrackSelector::~HICaloCompatibleTrackSelector() {}
0070
0071 void HICaloCompatibleTrackSelector::produce(edm::Event& evt, const edm::EventSetup& es) {
0072   using namespace std;
0073   using namespace edm;
0074   using namespace reco;
0075
0076   LogDebug("HICaloCompatibleTrackSelector") << "min pt for selection = " << trkMatchPtMin_ << endl;
0077
0078   Handle<TrackCollection> hSrcTrack;
0079   Handle<vector<Trajectory>> hTraj;
0080   Handle<vector<Trajectory>> hTrajP;
0081   Handle<TrajTrackAssociationCollection> hTTAss;
0082
0083   evt.getByToken(srcTracks_, hSrcTrack);
0084
0085   selTracks_ = std::make_unique<TrackCollection>();
0086   rTracks_ = evt.getRefBeforePut<TrackCollection>();
0087   if (copyExtras_) {
0088     selTrackExtras_ = std::make_unique<TrackExtraCollection>();
0089     selHits_ = std::make_unique<TrackingRecHitCollection>();
0090     rHits_ = evt.getRefBeforePut<TrackingRecHitCollection>();
0091     rTrackExtras_ = evt.getRefBeforePut<TrackExtraCollection>();
0092   }
0093
0094   if (copyTrajectories_)
0095     trackRefs_.resize(hSrcTrack->size());
0096
0097   Handle<PFCandidateCollection> pfCandidates;
0098   Handle<CaloTowerCollection> towers;
0099
0100   bool isPFThere = false;
0101   bool isTowerThere = false;
0102
0103   if (usePFCandMatching_)
0104     isPFThere = evt.getByToken(srcPFCands_, pfCandidates);
0105   else
0106     isTowerThere = evt.getByToken(srcTower_, towers);
0107
0108   size_t current = 0;
0109   for (TI ti = hSrcTrack->begin(), ed = hSrcTrack->end(); ti != ed; ++ti, ++current) {
0110     const reco::Track& trk = *ti;
0111
0112     bool isSelected;
0113     if (usePFCandMatching_)
0114       isSelected = selectByPFCands(ti, hSrcTrack, pfCandidates, isPFThere);
0115     else
0116       isSelected = selectByTowers(ti, hSrcTrack, towers, isTowerThere);
0117
0118     if (!keepAllTracks_ && !isSelected)
0119       continue;
0120
0121     // Add all tracks to output collection, the rest of the code only sets the quality bit
0122     selTracks_->push_back(reco::Track(trk));  // clone and store
0123
0124     if (isSelected)
0125       selTracks_->back().setQuality(reco::TrackBase::qualityByName(qualityToSet_));
0126
0127     if (copyExtras_) {
0128       // TrackExtras
0129       selTrackExtras_->push_back(TrackExtra(trk.outerPosition(),
0130                                             trk.outerMomentum(),
0131                                             trk.outerOk(),
0132                                             trk.innerPosition(),
0133                                             trk.innerMomentum(),
0134                                             trk.innerOk(),
0135                                             trk.outerStateCovariance(),
0136                                             trk.outerDetId(),
0137                                             trk.innerStateCovariance(),
0138                                             trk.innerDetId(),
0139                                             trk.seedDirection(),
0140                                             trk.seedRef()));
0141       selTracks_->back().setExtra(TrackExtraRef(rTrackExtras_, selTrackExtras_->size() - 1));
0142       TrackExtra& tx = selTrackExtras_->back();
0143       tx.setResiduals(trk.residuals());
0144       // TrackingRecHits
0145       auto const firstHitIndex = selHits_->size();
0146       for (trackingRecHit_iterator hit = trk.recHitsBegin(); hit != trk.recHitsEnd(); ++hit) {
0147         selHits_->push_back((*hit)->clone());
0148       }
0149       tx.setHits(rHits_, firstHitIndex, selHits_->size() - firstHitIndex);
0150     }
0151     if (copyTrajectories_) {
0152       trackRefs_[current] = TrackRef(rTracks_, selTracks_->size() - 1);
0153     }
0154
0155   }  // close track loop
0156
0157   if (copyTrajectories_) {
0158     Handle<vector<Trajectory>> hTraj;
0159     Handle<TrajTrackAssociationCollection> hTTAss;
0160     evt.getByToken(srcTrackTrajs_, hTTAss);
0161     evt.getByToken(srcTrackTrajAssoc_, hTraj);
0162     selTrajs_ = std::make_unique<std::vector<Trajectory>>();
0163     rTrajectories_ = evt.getRefBeforePut<vector<Trajectory>>();
0164     selTTAss_ = std::make_unique<TrajTrackAssociationCollection>();
0165     for (size_t i = 0, n = hTraj->size(); i < n; ++i) {
0166       Ref<vector<Trajectory>> trajRef(hTraj, i);
0167       TrajTrackAssociationCollection::const_iterator match = hTTAss->find(trajRef);
0168       if (match != hTTAss->end()) {
0169         const Ref<TrackCollection>& trkRef = match->val;
0170         short oldKey = static_cast<short>(trkRef.key());
0171         if (trackRefs_[oldKey].isNonnull()) {
0172           selTrajs_->push_back(Trajectory(*trajRef));
0173           selTTAss_->insert(Ref<vector<Trajectory>>(rTrajectories_, selTrajs_->size() - 1), trackRefs_[oldKey]);
0174         }
0175       }
0176     }
0177   }
0178
0179   static const string emptyString;
0180   evt.put(std::move(selTracks_));
0181   if (copyExtras_) {
0182     evt.put(std::move(selTrackExtras_));
0183     evt.put(std::move(selHits_));
0184   }
0185   if (copyTrajectories_) {
0186     evt.put(std::move(selTrajs_));
0187     evt.put(std::move(selTTAss_));
0188   }
0189 }
0190
0191 bool HICaloCompatibleTrackSelector::selectByPFCands(TI ti,
0192                                                     edm::Handle<TrackCollection> hSrcTrack,
0193                                                     edm::Handle<PFCandidateCollection> pfCandidates,
0194                                                     bool isPFThere) {
0195   const reco::Track& trk = *ti;
0196
0197   // If it passes this quality threshold or is under the minimum match pT, automatically save it
0198   if (trk.quality(reco::TrackBase::qualityByName(qualityToSkip_))) {
0199     return true;
0200   } else if (!trk.quality(reco::TrackBase::qualityByName(minimumQuality_))) {
0201     return false;
0202   } else {
0203     double trkPt = trk.pt();
0204     //if(trkPt < trkMatchPtMin_ ) return false;
0205
0206     double caloEt = 0.0;
0207     if (usePFCandMatching_) {
0208       if (isPFThere) {
0209         unsigned int trackKey = ti - hSrcTrack->begin();
0210         caloEt = matchPFCandToTrack(pfCandidates, trackKey, trkPt);
0211       }
0212     }
0213
0214     // Set quality bit based on calo matching
0215     if (!(caloEt > 0.))
0216       return false;
0217
0218     if (trkPt <= trkCompPtMin_) {
0219       if (trk.quality(reco::TrackBase::qualityByName(qualityToMatch_)))
0220         return true;
0221       else
0222         return false;
0223     } else {
0224       // loose cuts are implied in selectors, make configurable?
0225       float compPt =
0226           (fCaloComp->Eval(trkPt) != fCaloComp->Eval(trkPt)) ? 0 : fCaloComp->Eval(trkPt);  // protect agains NaN
0227       if (caloEt > compPt)
0228         return true;
0229       else
0230         return false;
0231     }
0232   }  // else above trkMatchPtMin_
0233
0234   throw cms::Exception("Undefined case in HICaloCompatibleTrackSelector")
0235       << "Undefined case in HICaloCompatibleTrackSelector";
0236 }
0237
0238 bool HICaloCompatibleTrackSelector::selectByTowers(TI ti,
0239                                                    edm::Handle<TrackCollection> hSrcTrack,
0240                                                    edm::Handle<CaloTowerCollection> towers,
0241                                                    bool isTowerThere) {
0242   // Not up to date! use PF towers instead
0243
0244   const reco::Track& trk = *ti;
0245
0246   // If it passes the high purity cuts, then consider it confirmed
0247   if (trk.quality(reco::TrackBase::qualityByName(qualityToSkip_))) {
0248     return true;
0249   } else {
0250     if (trk.pt() < trkMatchPtMin_ || !trk.quality(reco::TrackBase::qualityByName(qualityToMatch_)))
0251       return false;
0252
0253     double caloEt = 0.0;
0254     if (isTowerThere) {
0255       double matchDr;
0256       matchByDrAllowReuse(trk, towers, matchDr, caloEt);
0257       float matchConeRadius_pt = (fDeltaRTowerMatch->Eval(trk.pt()) != fDeltaRTowerMatch->Eval(trk.pt()))
0258                                      ? 0
0259                                      : fDeltaRTowerMatch->Eval(trk.pt());  // protect agains NaN
0260       if (caloEt > 0 && matchDr > matchConeRadius_pt)
0261         caloEt = 0.;
0262     }
0263
0264     if (trk.pt() < trkCompPtMin_ || caloEt > 0.75 * (trk.pt() - trkCompPtMin_))
0265       return true;
0266     else
0267       return false;
0268   }
0269 }
0270
0271 double HICaloCompatibleTrackSelector::matchPFCandToTrack(const edm::Handle<reco::PFCandidateCollection>& pfCandidates,
0272                                                          unsigned it,
0273                                                          double trkPt) {
0274   // This function returns the sum of the calo energy in the block containing the track
0275   // There is another piece of information which could be useful which is the pT assigned to the charged hadron by PF
0276
0277   double sum_ecal = 0.0, sum_hcal = 0.0;
0278
0279   int candType = 0;
0280   reco::PFCandidate matchedCand;
0281
0282   // loop over the PFCandidates until you find the one whose trackRef points to the track
0283
0284   for (CI ci = pfCandidates->begin(); ci != pfCandidates->end(); ++ci) {
0285     const reco::PFCandidate& cand = *ci;
0286
0287     int type = cand.particleId();
0288
0289     // only charged hadrons and leptons can be asscociated with a track
0290     if (!(type == PFCandidate::h ||  //type1
0291           type == PFCandidate::e ||  //type2
0292           type == PFCandidate::mu    //type3
0293           ))
0294       continue;
0295
0296     unsigned candTrackRefKey = cand.trackRef().key();
0297
0298     if (it == candTrackRefKey) {
0299       matchedCand = cand;
0300       candType = type;
0301       break;
0302     }
0303   }
0304   // take all muons as compatible, extend to electrons when validataed
0305   if (passMuons_ && candType == 3)
0306     return 9999.;
0307   if (passElectrons_ && candType == 2)
0308     return 9999.;
0309
0310   if (trkPt < trkMatchPtMin_)
0311     return 0.;
0312
0313   if (candType > 0) {
0314     // Now that we found the matched PF candidate, loop over the elements in the block summing the calo Et
0315     for (unsigned ib = 0; ib < matchedCand.elementsInBlocks().size(); ib++) {
0316       PFBlockRef blockRef = matchedCand.elementsInBlocks()[ib].first;
0317
0318       unsigned indexInBlock = matchedCand.elementsInBlocks()[ib].second;
0319       const edm::OwnVector<reco::PFBlockElement>& elements = (*blockRef).elements();
0320
0321       //This tells you what type of element it is:
0322       //cout<<" block type"<<elements[indexInBlock].type()<<endl;
0323
0324       switch (elements[indexInBlock].type()) {
0325         case PFBlockElement::ECAL: {
0326           reco::PFClusterRef clusterRef = elements[indexInBlock].clusterRef();
0327           sum_ecal += clusterRef->energy() / cosh(clusterRef->eta());
0328           break;
0329         }
0330
0331         case PFBlockElement::HCAL: {
0332           reco::PFClusterRef clusterRef = elements[indexInBlock].clusterRef();
0333           sum_hcal += clusterRef->energy() / cosh(clusterRef->eta());
0334           break;
0335         }
0336         case PFBlockElement::TRACK: {
0337           //Do nothing since track are not normally linked to other tracks
0338           break;
0339         }
0340         default:
0341           break;
0342       }
0343
0344       // We could also add in the PS, HO, ..
0345
0346     }  // end of elementsInBlocks()
0347   }  // end of isCandFound
0348
0349   return sum_ecal + sum_hcal;
0350 }
0351
0352 void HICaloCompatibleTrackSelector::matchByDrAllowReuse(const reco::Track& trk,
0353                                                         const edm::Handle<CaloTowerCollection>& towers,
0354                                                         double& bestdr,
0355                                                         double& bestpt) {
0356   // loop over towers
0357   bestdr = 1e10;
0358   bestpt = 0.;
0359   for (unsigned int i = 0; i < towers->size(); ++i) {
0360     const CaloTower& tower = (*towers)[i];
0361     double pt = tower.pt();
0362     if (pt < towerPtMin_)
0363       continue;
0364     if (fabs(tower.eta()) > trkEtaMax_)
0365       continue;
0366     double dr = reco::deltaR(tower, trk);
0367     if (dr < bestdr) {
0368       bestdr = dr;
0369       bestpt = pt;
0370     }
0371   }
0372 }