TkAlCaRecoProducers/src/CalibrationTrackSelector.cc

0001 #include "Calibration/TkAlCaRecoProducers/interface/CalibrationTrackSelector.h"
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005
0006 #include "FWCore/Framework/interface/Event.h"
0007
0008 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0009 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
0010 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0011 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0012 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0013
0014 #include "DataFormats/DetId/interface/DetId.h"
0015 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0016 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0017
0018 const int kBPIX = PixelSubdetector::PixelBarrel;
0019 const int kFPIX = PixelSubdetector::PixelEndcap;
0020
0021 // constructor ----------------------------------------------------------------
0022
0023 CalibrationTrackSelector::CalibrationTrackSelector(const edm::ParameterSet &cfg, edm::ConsumesCollector &iC)
0024     : applyBasicCuts_(cfg.getParameter<bool>("applyBasicCuts")),
0025       applyNHighestPt_(cfg.getParameter<bool>("applyNHighestPt")),
0026       applyMultiplicityFilter_(cfg.getParameter<bool>("applyMultiplicityFilter")),
0027       seedOnlyFromAbove_(cfg.getParameter<int>("seedOnlyFrom")),
0028       applyIsolation_(cfg.getParameter<bool>("applyIsolationCut")),
0029       chargeCheck_(cfg.getParameter<bool>("applyChargeCheck")),
0030       nHighestPt_(cfg.getParameter<int>("nHighestPt")),
0031       minMultiplicity_(cfg.getParameter<int>("minMultiplicity")),
0032       maxMultiplicity_(cfg.getParameter<int>("maxMultiplicity")),
0033       multiplicityOnInput_(cfg.getParameter<bool>("multiplicityOnInput")),
0034       ptMin_(cfg.getParameter<double>("ptMin")),
0035       ptMax_(cfg.getParameter<double>("ptMax")),
0036       etaMin_(cfg.getParameter<double>("etaMin")),
0037       etaMax_(cfg.getParameter<double>("etaMax")),
0038       phiMin_(cfg.getParameter<double>("phiMin")),
0039       phiMax_(cfg.getParameter<double>("phiMax")),
0040       nHitMin_(cfg.getParameter<double>("nHitMin")),
0041       nHitMax_(cfg.getParameter<double>("nHitMax")),
0042       chi2nMax_(cfg.getParameter<double>("chi2nMax")),
0043       minHitChargeStrip_(cfg.getParameter<double>("minHitChargeStrip")),
0044       minHitIsolation_(cfg.getParameter<double>("minHitIsolation")),
0045       rphirecHitsTag_(cfg.getParameter<edm::InputTag>("rphirecHits")),
0046       matchedrecHitsTag_(cfg.getParameter<edm::InputTag>("matchedrecHits")),
0047       nHitMin2D_(cfg.getParameter<unsigned int>("nHitMin2D")),
0048       // Ugly to use the same getParameter 6 times, but this allows const cut
0049       // variables...
0050       minHitsinTIB_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inTIB")),
0051       minHitsinTOB_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inTOB")),
0052       minHitsinTID_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inTID")),
0053       minHitsinTEC_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inTEC")),
0054       minHitsinBPIX_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inBPIX")),
0055       minHitsinFPIX_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inFPIX")),
0056       rphirecHitsToken_(iC.consumes<SiStripRecHit2DCollection>(rphirecHitsTag_)),
0057       matchedrecHitsToken_(iC.consumes<SiStripMatchedRecHit2DCollection>(matchedrecHitsTag_)) {
0058   if (applyBasicCuts_)
0059     edm::LogInfo("CalibrationTrackSelector")
0060         << "applying basic track cuts ..."
0061         << "\nptmin,ptmax:     " << ptMin_ << "," << ptMax_ << "\netamin,etamax:   " << etaMin_ << "," << etaMax_
0062         << "\nphimin,phimax:   " << phiMin_ << "," << phiMax_ << "\nnhitmin,nhitmax: " << nHitMin_ << "," << nHitMax_
0063         << "\nnhitmin2D:       " << nHitMin2D_ << "\nchi2nmax:        " << chi2nMax_;
0064
0065   if (applyNHighestPt_)
0066     edm::LogInfo("CalibrationTrackSelector") << "filter N tracks with highest Pt N=" << nHighestPt_;
0067
0068   if (applyMultiplicityFilter_)
0069     edm::LogInfo("CalibrationTrackSelector")
0070         << "apply multiplicity filter N>= " << minMultiplicity_ << "and N<= " << maxMultiplicity_ << " on "
0071         << (multiplicityOnInput_ ? "input" : "output");
0072
0073   if (applyIsolation_)
0074     edm::LogInfo("CalibrationTrackSelector")
0075         << "only retain tracks isolated at least by " << minHitIsolation_ << " cm from other rechits";
0076
0077   if (chargeCheck_)
0078     edm::LogInfo("CalibrationTrackSelector")
0079         << "only retain hits with at least " << minHitChargeStrip_ << " ADC counts of total cluster charge";
0080
0081   edm::LogInfo("CalibrationTrackSelector")
0082       << "Minimum number of hits in TIB/TID/TOB/TEC/BPIX/FPIX = " << minHitsinTIB_ << "/" << minHitsinTID_ << "/"
0083       << minHitsinTOB_ << "/" << minHitsinTEC_ << "/" << minHitsinBPIX_ << "/" << minHitsinFPIX_;
0084 }
0085
0086 // fillDescriptions ----------------------------------------------------------
0087
0088 void CalibrationTrackSelector::fillPSetDescription(edm::ParameterSetDescription &desc) {
0089   desc.add<double>("minHitChargeStrip", 20.0);
0090   desc.add<edm::InputTag>("rphirecHits", edm::InputTag("siStripMatchedRecHits", "rphiRecHit"));
0091   desc.add<bool>("applyMultiplicityFilter", false);
0092   desc.add<edm::InputTag>("matchedrecHits", edm::InputTag("siStripMatchedRecHits", "matchedRecHit"));
0093   desc.add<double>("etaMin", -2.6);
0094   desc.add<double>("etaMax", 2.6);
0095   desc.add<double>("minHitIsolation", 0.01);
0096   desc.add<double>("phiMax", 3.1416);
0097   desc.add<double>("phiMin", -3.1416);
0098   desc.add<double>("ptMin", 10.0);
0099   desc.add<int>("minMultiplicity", 1);
0100   desc.add<double>("nHitMin", 0.0);
0101   desc.add<double>("ptMax", 999.0);
0102   desc.add<double>("nHitMax", 999.0);
0103   desc.add<bool>("applyNHighestPt", false);
0104   desc.add<bool>("applyChargeCheck", false);
0105
0106   // Nested PSet for minHitsPerSubDet
0107   edm::ParameterSetDescription minHitsPerSubDetDesc;
0108   minHitsPerSubDetDesc.add<int>("inTEC", 0);
0109   minHitsPerSubDetDesc.add<int>("inTOB", 0);
0110   minHitsPerSubDetDesc.add<int>("inFPIX", 0);
0111   minHitsPerSubDetDesc.add<int>("inTID", 0);
0112   minHitsPerSubDetDesc.add<int>("inBPIX", 0);
0113   minHitsPerSubDetDesc.add<int>("inTIB", 0);
0114   desc.add<edm::ParameterSetDescription>("minHitsPerSubDet", minHitsPerSubDetDesc);
0115
0116   desc.add<int>("nHighestPt", 2);
0117   desc.add<unsigned int>("nHitMin2D", 0);
0118
0119   desc.add<bool>("applyIsolationCut", false);
0120   desc.add<bool>("multiplicityOnInput", false);
0121   desc.add<int>("maxMultiplicity", 999999);
0122   desc.add<int>("seedOnlyFrom", 0);
0123   desc.add<double>("chi2nMax", 999999.0);
0124
0125   // Settings for base TrackSelector
0126   desc.add<bool>("applyBasicCuts", true);
0127 }
0128
0129 // destructor -----------------------------------------------------------------
0130
0131 CalibrationTrackSelector::~CalibrationTrackSelector() {}
0132
0133 // do selection ---------------------------------------------------------------
0134
0135 CalibrationTrackSelector::Tracks CalibrationTrackSelector::select(const Tracks &tracks, const edm::Event &evt) const {
0136   if (applyMultiplicityFilter_ && multiplicityOnInput_ &&
0137       (tracks.size() < static_cast<unsigned int>(minMultiplicity_) ||
0138        tracks.size() > static_cast<unsigned int>(maxMultiplicity_))) {
0139     return Tracks();  // empty collection
0140   }
0141
0142   Tracks result = tracks;
0143   // apply basic track cuts (if selected)
0144   if (applyBasicCuts_)
0145     result = this->basicCuts(result, evt);
0146
0147   // filter N tracks with highest Pt (if selected)
0148   if (applyNHighestPt_)
0149     result = this->theNHighestPtTracks(result);
0150
0151   // apply minimum multiplicity requirement (if selected)
0152   if (applyMultiplicityFilter_ && !multiplicityOnInput_) {
0153     if (result.size() < static_cast<unsigned int>(minMultiplicity_) ||
0154         result.size() > static_cast<unsigned int>(maxMultiplicity_)) {
0155       result.clear();
0156     }
0157   }
0158
0159   return result;
0160 }
0161
0162 // make basic cuts ------------------------------------------------------------
0163
0164 CalibrationTrackSelector::Tracks CalibrationTrackSelector::basicCuts(const Tracks &tracks,
0165                                                                      const edm::Event &evt) const {
0166   Tracks result;
0167
0168   for (Tracks::const_iterator it = tracks.begin(); it != tracks.end(); ++it) {
0169     const reco::Track *trackp = *it;
0170     float pt = trackp->pt();
0171     float eta = trackp->eta();
0172     float phi = trackp->phi();
0173     int nhit = trackp->numberOfValidHits();
0174     float chi2n = trackp->normalizedChi2();
0175
0176     if (pt > ptMin_ && pt < ptMax_ && eta > etaMin_ && eta < etaMax_ && phi > phiMin_ && phi < phiMax_ &&
0177         nhit >= nHitMin_ && nhit <= nHitMax_ && chi2n < chi2nMax_) {
0178       if (this->detailedHitsCheck(trackp, evt))
0179         result.push_back(trackp);
0180     }
0181   }
0182
0183   return result;
0184 }
0185
0186 //-----------------------------------------------------------------------------
0187
0188 bool CalibrationTrackSelector::detailedHitsCheck(const reco::Track *trackp, const edm::Event &evt) const {
0189   // checking hit requirements beyond simple number of valid hits
0190
0191   if (minHitsinTIB_ || minHitsinTOB_ || minHitsinTID_ || minHitsinTEC_ || minHitsinFPIX_ || minHitsinBPIX_ ||
0192       nHitMin2D_ || chargeCheck_ || applyIsolation_) {  // any detailed hit cut is active, so have to check
0193
0194     int nhitinTIB = 0, nhitinTOB = 0, nhitinTID = 0;
0195     int nhitinTEC = 0, nhitinBPIX = 0, nhitinFPIX = 0;
0196     unsigned int nHit2D = 0;
0197     unsigned int thishit = 0;
0198
0199     for (auto const &hit : trackp->recHits()) {
0200       thishit++;
0201       int type = hit->geographicalId().subdetId();
0202
0203       // *** thishit == 1 means last hit in CTF ***
0204       // (FIXME: assumption might change or not be valid for all tracking
0205       // algorthms)
0206       // ==> for cosmics
0207       // seedOnlyFrom = 1 is TIB-TOB-TEC tracks only
0208       // seedOnlyFrom = 2 is TOB-TEC tracks only
0209       if (seedOnlyFromAbove_ == 1 && thishit == 1 &&
0210           (type == int(StripSubdetector::TOB) || type == int(StripSubdetector::TEC)))
0211         return false;
0212
0213       if (seedOnlyFromAbove_ == 2 && thishit == 1 && type == int(StripSubdetector::TIB))
0214         return false;
0215
0216       if (!hit->isValid())
0217         continue;  // only real hits count as in trackp->numberOfValidHits()
0218       const DetId detId(hit->geographicalId());
0219       if (detId.det() != DetId::Tracker) {
0220         edm::LogError("DetectorMismatch")
0221             << "@SUB=CalibrationTrackSelector::detailedHitsCheck"
0222             << "DetId.det() != DetId::Tracker (=" << DetId::Tracker << "), but " << detId.det() << ".";
0223       }
0224       if (chargeCheck_ && !(this->isOkCharge(hit)))
0225         return false;
0226       if (applyIsolation_ && (!this->isIsolated(hit, evt)))
0227         return false;
0228       if (StripSubdetector::TIB == detId.subdetId())
0229         ++nhitinTIB;
0230       else if (StripSubdetector::TOB == detId.subdetId())
0231         ++nhitinTOB;
0232       else if (StripSubdetector::TID == detId.subdetId())
0233         ++nhitinTID;
0234       else if (StripSubdetector::TEC == detId.subdetId())
0235         ++nhitinTEC;
0236       else if (kBPIX == detId.subdetId())
0237         ++nhitinBPIX;
0238       else if (kFPIX == detId.subdetId())
0239         ++nhitinFPIX;
0240       // Do not call isHit2D(..) if already enough 2D hits for performance
0241       // reason:
0242       if (nHit2D < nHitMin2D_ && this->isHit2D(*hit))
0243         ++nHit2D;
0244     }  // end loop on hits
0245     return (nhitinTIB >= minHitsinTIB_ && nhitinTOB >= minHitsinTOB_ && nhitinTID >= minHitsinTID_ &&
0246             nhitinTEC >= minHitsinTEC_ && nhitinBPIX >= minHitsinBPIX_ && nhitinFPIX >= minHitsinFPIX_ &&
0247             nHit2D >= nHitMin2D_);
0248   } else {  // no cuts set, so we are just fine and can avoid loop on hits
0249     return true;
0250   }
0251 }
0252
0253 //-----------------------------------------------------------------------------
0254
0255 bool CalibrationTrackSelector::isHit2D(const TrackingRecHit &hit) const {
0256   if (hit.dimension() < 2) {
0257     return false;  // some (muon...) stuff really has RecHit1D
0258   } else {
0259     const DetId detId(hit.geographicalId());
0260     if (detId.det() == DetId::Tracker) {
0261       if (detId.subdetId() == kBPIX || detId.subdetId() == kFPIX) {
0262         return true;  // pixel is always 2D
0263       } else {        // should be SiStrip now
0264         if (dynamic_cast<const SiStripRecHit2D *>(&hit))
0265           return false;  // normal hit
0266         else if (dynamic_cast<const SiStripMatchedRecHit2D *>(&hit))
0267           return true;  // matched is 2D
0268         else if (dynamic_cast<const ProjectedSiStripRecHit2D *>(&hit))
0269           return false;  // crazy hit...
0270         else {
0271           edm::LogError("UnknownType") << "@SUB=CalibrationTrackSelector::isHit2D"
0272                                        << "Tracker hit not in pixel and neither SiStripRecHit2D nor "
0273                                        << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
0274           return false;
0275         }
0276       }
0277     } else {  // not tracker??
0278       edm::LogWarning("DetectorMismatch") << "@SUB=CalibrationTrackSelector::isHit2D"
0279                                           << "Hit not in tracker with 'official' dimension >=2.";
0280       return true;  // dimension() >= 2 so accept that...
0281     }
0282   }
0283   // never reached...
0284 }
0285
0286 //-----------------------------------------------------------------------------
0287
0288 bool CalibrationTrackSelector::isOkCharge(const TrackingRecHit *therechit) const {
0289   float charge1 = 0;
0290   float charge2 = 0;
0291   const SiStripMatchedRecHit2D *matchedhit = dynamic_cast<const SiStripMatchedRecHit2D *>(therechit);
0292   const SiStripRecHit2D *hit = dynamic_cast<const SiStripRecHit2D *>(therechit);
0293   const ProjectedSiStripRecHit2D *unmatchedhit = dynamic_cast<const ProjectedSiStripRecHit2D *>(therechit);
0294
0295   if (matchedhit) {
0296     const SiStripCluster &monocluster = matchedhit->monoCluster();
0297     const std::vector<uint16_t> amplitudesmono(monocluster.amplitudes().begin(), monocluster.amplitudes().end());
0298     for (size_t ia = 0; ia < amplitudesmono.size(); ++ia) {
0299       charge1 += amplitudesmono[ia];
0300     }
0301
0302     const SiStripCluster &stereocluster = matchedhit->stereoCluster();
0303     const std::vector<uint16_t> amplitudesstereo(stereocluster.amplitudes().begin(), stereocluster.amplitudes().end());
0304     for (size_t ia = 0; ia < amplitudesstereo.size(); ++ia) {
0305       charge2 += amplitudesstereo[ia];
0306     }
0307     if (charge1 < minHitChargeStrip_ || charge2 < minHitChargeStrip_)
0308       return false;
0309   } else if (hit) {
0310     const SiStripCluster *cluster = &*(hit->cluster());
0311     const std::vector<uint16_t> amplitudes(cluster->amplitudes().begin(), cluster->amplitudes().end());
0312     for (size_t ia = 0; ia < amplitudes.size(); ++ia) {
0313       charge1 += amplitudes[ia];
0314     }
0315     if (charge1 < minHitChargeStrip_)
0316       return false;
0317   } else if (unmatchedhit) {
0318     const SiStripRecHit2D &orighit = unmatchedhit->originalHit();
0319     const SiStripCluster *origcluster = &*(orighit.cluster());
0320     const std::vector<uint16_t> amplitudes(origcluster->amplitudes().begin(), origcluster->amplitudes().end());
0321     for (size_t ia = 0; ia < amplitudes.size(); ++ia) {
0322       charge1 += amplitudes[ia];
0323     }
0324     if (charge1 < minHitChargeStrip_)
0325       return false;
0326   }
0327   return true;
0328 }
0329
0330 //-----------------------------------------------------------------------------
0331
0332 bool CalibrationTrackSelector::isIsolated(const TrackingRecHit *therechit, const edm::Event &evt) const {
0333   const edm::Handle<SiStripRecHit2DCollection> &rphirecHits = evt.getHandle(rphirecHitsToken_);
0334   const edm::Handle<SiStripMatchedRecHit2DCollection> &matchedrecHits = evt.getHandle(matchedrecHitsToken_);
0335
0336   SiStripRecHit2DCollection::DataContainer::const_iterator istripSt;
0337   SiStripMatchedRecHit2DCollection::DataContainer::const_iterator istripStm;
0338   const SiStripRecHit2DCollection &stripcollSt = *rphirecHits;
0339   const SiStripMatchedRecHit2DCollection &stripcollStm = *matchedrecHits;
0340
0341   DetId idet = therechit->geographicalId();
0342
0343   // FIXME: instead of looping the full hit collection, we should explore the
0344   // features of SiStripRecHit2DCollection::rangeRphi = rphirecHits.get(idet)
0345   // and loop only from rangeRphi.first until rangeRphi.second
0346   for (istripSt = stripcollSt.data().begin(); istripSt != stripcollSt.data().end(); ++istripSt) {
0347     const SiStripRecHit2D *aHit = &*(istripSt);
0348     DetId mydet1 = aHit->geographicalId();
0349     if (idet.rawId() != mydet1.rawId())
0350       continue;
0351     float theDistance = (therechit->localPosition() - aHit->localPosition()).mag();
0352     if (theDistance > 0.001 && theDistance < minHitIsolation_)
0353       return false;
0354   }
0355
0356   // FIXME: see above
0357   for (istripStm = stripcollStm.data().begin(); istripStm != stripcollStm.data().end(); ++istripStm) {
0358     const SiStripMatchedRecHit2D *aHit = &*(istripStm);
0359     DetId mydet2 = aHit->geographicalId();
0360     if (idet.rawId() != mydet2.rawId())
0361       continue;
0362     float theDistance = (therechit->localPosition() - aHit->localPosition()).mag();
0363     if (theDistance > 0.001 && theDistance < minHitIsolation_)
0364       return false;
0365   }
0366   return true;
0367 }
0368 //-----------------------------------------------------------------------------
0369
0370 CalibrationTrackSelector::Tracks CalibrationTrackSelector::theNHighestPtTracks(const Tracks &tracks) const {
0371   Tracks sortedTracks = tracks;
0372   Tracks result;
0373
0374   // sort in pt
0375   std::sort(sortedTracks.begin(), sortedTracks.end(), ptComparator);
0376
0377   // copy theTrackMult highest pt tracks to result vector
0378   int n = 0;
0379   for (Tracks::const_iterator it = sortedTracks.begin(); it != sortedTracks.end(); ++it) {
0380     if (n < nHighestPt_) {
0381       result.push_back(*it);
0382       n++;
0383     }
0384   }
0385
0386   return result;
0387 }