File indexing completed on 2025-02-05 23:51:00
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 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
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
0126 desc.add<bool>("applyBasicCuts", true);
0127 }
0128
0129
0130
0131 CalibrationTrackSelector::~CalibrationTrackSelector() {}
0132
0133
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();
0140 }
0141
0142 Tracks result = tracks;
0143
0144 if (applyBasicCuts_)
0145 result = this->basicCuts(result, evt);
0146
0147
0148 if (applyNHighestPt_)
0149 result = this->theNHighestPtTracks(result);
0150
0151
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
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
0190
0191 if (minHitsinTIB_ || minHitsinTOB_ || minHitsinTID_ || minHitsinTEC_ || minHitsinFPIX_ || minHitsinBPIX_ ||
0192 nHitMin2D_ || chargeCheck_ || applyIsolation_) {
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
0204
0205
0206
0207
0208
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;
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
0241
0242 if (nHit2D < nHitMin2D_ && this->isHit2D(*hit))
0243 ++nHit2D;
0244 }
0245 return (nhitinTIB >= minHitsinTIB_ && nhitinTOB >= minHitsinTOB_ && nhitinTID >= minHitsinTID_ &&
0246 nhitinTEC >= minHitsinTEC_ && nhitinBPIX >= minHitsinBPIX_ && nhitinFPIX >= minHitsinFPIX_ &&
0247 nHit2D >= nHitMin2D_);
0248 } else {
0249 return true;
0250 }
0251 }
0252
0253
0254
0255 bool CalibrationTrackSelector::isHit2D(const TrackingRecHit &hit) const {
0256 if (hit.dimension() < 2) {
0257 return false;
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;
0263 } else {
0264 if (dynamic_cast<const SiStripRecHit2D *>(&hit))
0265 return false;
0266 else if (dynamic_cast<const SiStripMatchedRecHit2D *>(&hit))
0267 return true;
0268 else if (dynamic_cast<const ProjectedSiStripRecHit2D *>(&hit))
0269 return false;
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 {
0278 edm::LogWarning("DetectorMismatch") << "@SUB=CalibrationTrackSelector::isHit2D"
0279 << "Hit not in tracker with 'official' dimension >=2.";
0280 return true;
0281 }
0282 }
0283
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
0344
0345
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
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
0375 std::sort(sortedTracks.begin(), sortedTracks.end(), ptComparator);
0376
0377
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 }