File indexing completed on 2023-03-17 10:43:52
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
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
0049
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
0087
0088 CalibrationTrackSelector::~CalibrationTrackSelector() {}
0089
0090
0091
0092 CalibrationTrackSelector::Tracks CalibrationTrackSelector::select(const Tracks &tracks, const edm::Event &evt) const {
0093 if (applyMultiplicityFilter_ && multiplicityOnInput_ &&
0094 (tracks.size() < static_cast<unsigned int>(minMultiplicity_) ||
0095 tracks.size() > static_cast<unsigned int>(maxMultiplicity_))) {
0096 return Tracks();
0097 }
0098
0099 Tracks result = tracks;
0100
0101 if (applyBasicCuts_)
0102 result = this->basicCuts(result, evt);
0103
0104
0105 if (applyNHighestPt_)
0106 result = this->theNHighestPtTracks(result);
0107
0108
0109 if (applyMultiplicityFilter_ && !multiplicityOnInput_) {
0110 if (result.size() < static_cast<unsigned int>(minMultiplicity_) ||
0111 result.size() > static_cast<unsigned int>(maxMultiplicity_)) {
0112 result.clear();
0113 }
0114 }
0115
0116 return result;
0117 }
0118
0119
0120
0121 CalibrationTrackSelector::Tracks CalibrationTrackSelector::basicCuts(const Tracks &tracks,
0122 const edm::Event &evt) const {
0123 Tracks result;
0124
0125 for (Tracks::const_iterator it = tracks.begin(); it != tracks.end(); ++it) {
0126 const reco::Track *trackp = *it;
0127 float pt = trackp->pt();
0128 float eta = trackp->eta();
0129 float phi = trackp->phi();
0130 int nhit = trackp->numberOfValidHits();
0131 float chi2n = trackp->normalizedChi2();
0132
0133 if (pt > ptMin_ && pt < ptMax_ && eta > etaMin_ && eta < etaMax_ && phi > phiMin_ && phi < phiMax_ &&
0134 nhit >= nHitMin_ && nhit <= nHitMax_ && chi2n < chi2nMax_) {
0135 if (this->detailedHitsCheck(trackp, evt))
0136 result.push_back(trackp);
0137 }
0138 }
0139
0140 return result;
0141 }
0142
0143
0144
0145 bool CalibrationTrackSelector::detailedHitsCheck(const reco::Track *trackp, const edm::Event &evt) const {
0146
0147
0148 if (minHitsinTIB_ || minHitsinTOB_ || minHitsinTID_ || minHitsinTEC_ || minHitsinFPIX_ || minHitsinBPIX_ ||
0149 nHitMin2D_ || chargeCheck_ || applyIsolation_) {
0150
0151 int nhitinTIB = 0, nhitinTOB = 0, nhitinTID = 0;
0152 int nhitinTEC = 0, nhitinBPIX = 0, nhitinFPIX = 0;
0153 unsigned int nHit2D = 0;
0154 unsigned int thishit = 0;
0155
0156 for (auto const &hit : trackp->recHits()) {
0157 thishit++;
0158 int type = hit->geographicalId().subdetId();
0159
0160
0161
0162
0163
0164
0165
0166 if (seedOnlyFromAbove_ == 1 && thishit == 1 &&
0167 (type == int(StripSubdetector::TOB) || type == int(StripSubdetector::TEC)))
0168 return false;
0169
0170 if (seedOnlyFromAbove_ == 2 && thishit == 1 && type == int(StripSubdetector::TIB))
0171 return false;
0172
0173 if (!hit->isValid())
0174 continue;
0175 const DetId detId(hit->geographicalId());
0176 if (detId.det() != DetId::Tracker) {
0177 edm::LogError("DetectorMismatch")
0178 << "@SUB=CalibrationTrackSelector::detailedHitsCheck"
0179 << "DetId.det() != DetId::Tracker (=" << DetId::Tracker << "), but " << detId.det() << ".";
0180 }
0181 if (chargeCheck_ && !(this->isOkCharge(hit)))
0182 return false;
0183 if (applyIsolation_ && (!this->isIsolated(hit, evt)))
0184 return false;
0185 if (StripSubdetector::TIB == detId.subdetId())
0186 ++nhitinTIB;
0187 else if (StripSubdetector::TOB == detId.subdetId())
0188 ++nhitinTOB;
0189 else if (StripSubdetector::TID == detId.subdetId())
0190 ++nhitinTID;
0191 else if (StripSubdetector::TEC == detId.subdetId())
0192 ++nhitinTEC;
0193 else if (kBPIX == detId.subdetId())
0194 ++nhitinBPIX;
0195 else if (kFPIX == detId.subdetId())
0196 ++nhitinFPIX;
0197
0198
0199 if (nHit2D < nHitMin2D_ && this->isHit2D(*hit))
0200 ++nHit2D;
0201 }
0202 return (nhitinTIB >= minHitsinTIB_ && nhitinTOB >= minHitsinTOB_ && nhitinTID >= minHitsinTID_ &&
0203 nhitinTEC >= minHitsinTEC_ && nhitinBPIX >= minHitsinBPIX_ && nhitinFPIX >= minHitsinFPIX_ &&
0204 nHit2D >= nHitMin2D_);
0205 } else {
0206 return true;
0207 }
0208 }
0209
0210
0211
0212 bool CalibrationTrackSelector::isHit2D(const TrackingRecHit &hit) const {
0213 if (hit.dimension() < 2) {
0214 return false;
0215 } else {
0216 const DetId detId(hit.geographicalId());
0217 if (detId.det() == DetId::Tracker) {
0218 if (detId.subdetId() == kBPIX || detId.subdetId() == kFPIX) {
0219 return true;
0220 } else {
0221 if (dynamic_cast<const SiStripRecHit2D *>(&hit))
0222 return false;
0223 else if (dynamic_cast<const SiStripMatchedRecHit2D *>(&hit))
0224 return true;
0225 else if (dynamic_cast<const ProjectedSiStripRecHit2D *>(&hit))
0226 return false;
0227 else {
0228 edm::LogError("UnknownType") << "@SUB=CalibrationTrackSelector::isHit2D"
0229 << "Tracker hit not in pixel and neither SiStripRecHit2D nor "
0230 << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
0231 return false;
0232 }
0233 }
0234 } else {
0235 edm::LogWarning("DetectorMismatch") << "@SUB=CalibrationTrackSelector::isHit2D"
0236 << "Hit not in tracker with 'official' dimension >=2.";
0237 return true;
0238 }
0239 }
0240
0241 }
0242
0243
0244
0245 bool CalibrationTrackSelector::isOkCharge(const TrackingRecHit *therechit) const {
0246 float charge1 = 0;
0247 float charge2 = 0;
0248 const SiStripMatchedRecHit2D *matchedhit = dynamic_cast<const SiStripMatchedRecHit2D *>(therechit);
0249 const SiStripRecHit2D *hit = dynamic_cast<const SiStripRecHit2D *>(therechit);
0250 const ProjectedSiStripRecHit2D *unmatchedhit = dynamic_cast<const ProjectedSiStripRecHit2D *>(therechit);
0251
0252 if (matchedhit) {
0253 const SiStripCluster &monocluster = matchedhit->monoCluster();
0254 const std::vector<uint16_t> amplitudesmono(monocluster.amplitudes().begin(), monocluster.amplitudes().end());
0255 for (size_t ia = 0; ia < amplitudesmono.size(); ++ia) {
0256 charge1 += amplitudesmono[ia];
0257 }
0258
0259 const SiStripCluster &stereocluster = matchedhit->stereoCluster();
0260 const std::vector<uint16_t> amplitudesstereo(stereocluster.amplitudes().begin(), stereocluster.amplitudes().end());
0261 for (size_t ia = 0; ia < amplitudesstereo.size(); ++ia) {
0262 charge2 += amplitudesstereo[ia];
0263 }
0264 if (charge1 < minHitChargeStrip_ || charge2 < minHitChargeStrip_)
0265 return false;
0266 } else if (hit) {
0267 const SiStripCluster *cluster = &*(hit->cluster());
0268 const std::vector<uint16_t> amplitudes(cluster->amplitudes().begin(), cluster->amplitudes().end());
0269 for (size_t ia = 0; ia < amplitudes.size(); ++ia) {
0270 charge1 += amplitudes[ia];
0271 }
0272 if (charge1 < minHitChargeStrip_)
0273 return false;
0274 } else if (unmatchedhit) {
0275 const SiStripRecHit2D &orighit = unmatchedhit->originalHit();
0276 const SiStripCluster *origcluster = &*(orighit.cluster());
0277 const std::vector<uint16_t> amplitudes(origcluster->amplitudes().begin(), origcluster->amplitudes().end());
0278 for (size_t ia = 0; ia < amplitudes.size(); ++ia) {
0279 charge1 += amplitudes[ia];
0280 }
0281 if (charge1 < minHitChargeStrip_)
0282 return false;
0283 }
0284 return true;
0285 }
0286
0287
0288
0289 bool CalibrationTrackSelector::isIsolated(const TrackingRecHit *therechit, const edm::Event &evt) const {
0290 const edm::Handle<SiStripRecHit2DCollection> &rphirecHits = evt.getHandle(rphirecHitsToken_);
0291 const edm::Handle<SiStripMatchedRecHit2DCollection> &matchedrecHits = evt.getHandle(matchedrecHitsToken_);
0292
0293 SiStripRecHit2DCollection::DataContainer::const_iterator istripSt;
0294 SiStripMatchedRecHit2DCollection::DataContainer::const_iterator istripStm;
0295 const SiStripRecHit2DCollection &stripcollSt = *rphirecHits;
0296 const SiStripMatchedRecHit2DCollection &stripcollStm = *matchedrecHits;
0297
0298 DetId idet = therechit->geographicalId();
0299
0300
0301
0302
0303 for (istripSt = stripcollSt.data().begin(); istripSt != stripcollSt.data().end(); ++istripSt) {
0304 const SiStripRecHit2D *aHit = &*(istripSt);
0305 DetId mydet1 = aHit->geographicalId();
0306 if (idet.rawId() != mydet1.rawId())
0307 continue;
0308 float theDistance = (therechit->localPosition() - aHit->localPosition()).mag();
0309 if (theDistance > 0.001 && theDistance < minHitIsolation_)
0310 return false;
0311 }
0312
0313
0314 for (istripStm = stripcollStm.data().begin(); istripStm != stripcollStm.data().end(); ++istripStm) {
0315 const SiStripMatchedRecHit2D *aHit = &*(istripStm);
0316 DetId mydet2 = aHit->geographicalId();
0317 if (idet.rawId() != mydet2.rawId())
0318 continue;
0319 float theDistance = (therechit->localPosition() - aHit->localPosition()).mag();
0320 if (theDistance > 0.001 && theDistance < minHitIsolation_)
0321 return false;
0322 }
0323 return true;
0324 }
0325
0326
0327 CalibrationTrackSelector::Tracks CalibrationTrackSelector::theNHighestPtTracks(const Tracks &tracks) const {
0328 Tracks sortedTracks = tracks;
0329 Tracks result;
0330
0331
0332 std::sort(sortedTracks.begin(), sortedTracks.end(), ptComparator);
0333
0334
0335 int n = 0;
0336 for (Tracks::const_iterator it = sortedTracks.begin(); it != sortedTracks.end(); ++it) {
0337 if (n < nHighestPt_) {
0338 result.push_back(*it);
0339 n++;
0340 }
0341 }
0342
0343 return result;
0344 }