Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2025-02-05 23:50:58

0001 #include "Alignment/CommonAlignmentProducer/interface/AlignmentTrackSelector.h"
0002 
0003 #include "FWCore/Framework/interface/ESHandle.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0006 
0007 #include "FWCore/Framework/interface/Event.h"
0008 
0009 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0010 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0011 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit1D.h"
0012 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0013 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
0014 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0015 #include "DataFormats/TrackerRecHit2D/interface/SiTrackerMultiRecHit.h"
0016 #include "DataFormats/SiStripCluster/interface/SiStripCluster.h"
0017 #include "DataFormats/Alignment/interface/AlignmentClusterFlag.h"
0018 #include "DataFormats/Alignment/interface/AliClusterValueMap.h"
0019 #include "DataFormats/DetId/interface/DetId.h"
0020 #include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0021 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0022 
0023 #include <cmath>
0024 
0025 const int kBPIX = PixelSubdetector::PixelBarrel;
0026 const int kFPIX = PixelSubdetector::PixelEndcap;
0027 
0028 // constructor ----------------------------------------------------------------
0029 
0030 AlignmentTrackSelector::AlignmentTrackSelector(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
0031     : tTopoToken_(iC.esConsumes()),
0032       applyBasicCuts_(cfg.getParameter<bool>("applyBasicCuts")),
0033       applyNHighestPt_(cfg.getParameter<bool>("applyNHighestPt")),
0034       applyMultiplicityFilter_(cfg.getParameter<bool>("applyMultiplicityFilter")),
0035       seedOnlyFromAbove_(cfg.getParameter<int>("seedOnlyFrom")),
0036       applyIsolation_(cfg.getParameter<bool>("applyIsolationCut")),
0037       chargeCheck_(cfg.getParameter<bool>("applyChargeCheck")),
0038       nHighestPt_(cfg.getParameter<int>("nHighestPt")),
0039       minMultiplicity_(cfg.getParameter<int>("minMultiplicity")),
0040       maxMultiplicity_(cfg.getParameter<int>("maxMultiplicity")),
0041       multiplicityOnInput_(cfg.getParameter<bool>("multiplicityOnInput")),
0042       ptMin_(cfg.getParameter<double>("ptMin")),
0043       ptMax_(cfg.getParameter<double>("ptMax")),
0044       pMin_(cfg.getParameter<double>("pMin")),
0045       pMax_(cfg.getParameter<double>("pMax")),
0046       etaMin_(cfg.getParameter<double>("etaMin")),
0047       etaMax_(cfg.getParameter<double>("etaMax")),
0048       phiMin_(cfg.getParameter<double>("phiMin")),
0049       phiMax_(cfg.getParameter<double>("phiMax")),
0050       nHitMin_(cfg.getParameter<double>("nHitMin")),
0051       nHitMax_(cfg.getParameter<double>("nHitMax")),
0052       chi2nMax_(cfg.getParameter<double>("chi2nMax")),
0053       d0Min_(cfg.getParameter<double>("d0Min")),
0054       d0Max_(cfg.getParameter<double>("d0Max")),
0055       dzMin_(cfg.getParameter<double>("dzMin")),
0056       dzMax_(cfg.getParameter<double>("dzMax")),
0057       theCharge_(cfg.getParameter<int>("theCharge")),
0058       minHitChargeStrip_(cfg.getParameter<double>("minHitChargeStrip")),
0059       minHitIsolation_(cfg.getParameter<double>("minHitIsolation")),
0060       countStereoHitAs2D_(cfg.getParameter<bool>("countStereoHitAs2D")),
0061       nHitMin2D_(cfg.getParameter<unsigned int>("nHitMin2D")),
0062       // Ugly to use the same getParameter n times, but this allows const cut variables...
0063       minHitsinTIB_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inTIB")),
0064       minHitsinTOB_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inTOB")),
0065       minHitsinTID_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inTID")),
0066       minHitsinTEC_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inTEC")),
0067       minHitsinBPIX_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inBPIX")),
0068       minHitsinFPIX_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inFPIX")),
0069       minHitsinPIX_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inPIXEL")),
0070       minHitsinTIDplus_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inTIDplus")),
0071       minHitsinTIDminus_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inTIDminus")),
0072       minHitsinTECplus_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inTECplus")),
0073       minHitsinTECminus_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inTECminus")),
0074       minHitsinFPIXplus_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inFPIXplus")),
0075       minHitsinFPIXminus_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inFPIXminus")),
0076       minHitsinENDCAP_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inENDCAP")),
0077       minHitsinENDCAPplus_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inENDCAPplus")),
0078       minHitsinENDCAPminus_(cfg.getParameter<edm::ParameterSet>("minHitsPerSubDet").getParameter<int>("inENDCAPminus")),
0079       maxHitDiffEndcaps_(cfg.getParameter<double>("maxHitDiffEndcaps")),
0080       nLostHitMax_(cfg.getParameter<double>("nLostHitMax")),
0081       RorZofFirstHitMin_(cfg.getParameter<std::vector<double>>("RorZofFirstHitMin")),
0082       RorZofFirstHitMax_(cfg.getParameter<std::vector<double>>("RorZofFirstHitMax")),
0083       RorZofLastHitMin_(cfg.getParameter<std::vector<double>>("RorZofLastHitMin")),
0084       RorZofLastHitMax_(cfg.getParameter<std::vector<double>>("RorZofLastHitMax")),
0085       clusterValueMapTag_(cfg.getParameter<edm::InputTag>("hitPrescaleMapTag")),
0086       minPrescaledHits_(cfg.getParameter<int>("minPrescaledHits")),
0087       applyPrescaledHitsFilter_(!clusterValueMapTag_.encode().empty() && minPrescaledHits_ > 0) {
0088   if (applyIsolation_) {
0089     rphirecHitsToken_ = iC.consumes<SiStripRecHit2DCollection>(cfg.getParameter<edm::InputTag>("rphirecHits"));
0090     matchedrecHitsToken_ =
0091         iC.consumes<SiStripMatchedRecHit2DCollection>(cfg.getParameter<edm::InputTag>("matchedrecHits"));
0092   }
0093 
0094   if (applyPrescaledHitsFilter_) {
0095     clusterValueMapToken_ = iC.consumes<AliClusterValueMap>(clusterValueMapTag_);
0096   }
0097 
0098   //convert track quality from string to enum
0099   std::vector<std::string> trkQualityStrings(cfg.getParameter<std::vector<std::string>>("trackQualities"));
0100   std::string qualities;
0101   if (!trkQualityStrings.empty()) {
0102     applyTrkQualityCheck_ = true;
0103     for (unsigned int i = 0; i < trkQualityStrings.size(); ++i) {
0104       (qualities += trkQualityStrings[i]) += ", ";
0105       trkQualities_.push_back(reco::TrackBase::qualityByName(trkQualityStrings[i]));
0106     }
0107   } else
0108     applyTrkQualityCheck_ = false;
0109 
0110   std::vector<std::string> trkIterStrings(cfg.getParameter<std::vector<std::string>>("iterativeTrackingSteps"));
0111   if (!trkIterStrings.empty()) {
0112     applyIterStepCheck_ = true;
0113     std::string tracksteps;
0114     for (unsigned int i = 0; i < trkIterStrings.size(); ++i) {
0115       (tracksteps += trkIterStrings[i]) += ", ";
0116       trkSteps_.push_back(reco::TrackBase::algoByName(trkIterStrings[i]));
0117     }
0118   } else
0119     applyIterStepCheck_ = false;
0120 
0121   if (applyBasicCuts_) {
0122     edm::LogInfo("AlignmentTrackSelector")
0123         << "applying basic track cuts ..."
0124         << "\nptmin,ptmax:     " << ptMin_ << "," << ptMax_ << "\npmin,pmax:       " << pMin_ << "," << pMax_
0125         << "\netamin,etamax:   " << etaMin_ << "," << etaMax_ << "\nphimin,phimax:   " << phiMin_ << "," << phiMax_
0126         << "\nnhitmin,nhitmax: " << nHitMin_ << "," << nHitMax_ << "\nnlosthitmax:     " << nLostHitMax_
0127         << "\nnhitmin2D:       " << nHitMin2D_ << (countStereoHitAs2D_ ? "," : ", not")
0128         << " counting hits on SiStrip stereo modules as 2D"
0129         << "\nchi2nmax:        " << chi2nMax_;
0130 
0131     if (applyIsolation_)
0132       edm::LogInfo("AlignmentTrackSelector")
0133           << "only retain tracks isolated at least by " << minHitIsolation_ << " cm from other rechits";
0134 
0135     if (chargeCheck_)
0136       edm::LogInfo("AlignmentTrackSelector")
0137           << "only retain hits with at least " << minHitChargeStrip_ << " ADC counts of total cluster charge";
0138 
0139     edm::LogInfo("AlignmentTrackSelector")
0140         << "Minimum number of hits in TIB/TID/TOB/TEC/BPIX/FPIX/PIXEL = " << minHitsinTIB_ << "/" << minHitsinTID_
0141         << "/" << minHitsinTOB_ << "/" << minHitsinTEC_ << "/" << minHitsinBPIX_ << "/" << minHitsinFPIX_ << "/"
0142         << minHitsinPIX_;
0143 
0144     edm::LogInfo("AlignmentTrackSelector")
0145         << "Minimum number of hits in TID+/TID-/TEC+/TEC-/FPIX+/FPIX- = " << minHitsinTIDplus_ << "/"
0146         << minHitsinTIDminus_ << "/" << minHitsinTECplus_ << "/" << minHitsinTECminus_ << "/" << minHitsinFPIXplus_
0147         << "/" << minHitsinFPIXminus_;
0148 
0149     edm::LogInfo("AlignmentTrackSelector")
0150         << "Minimum number of hits in EndCap (TID+TEC)/EndCap+/EndCap- = " << minHitsinENDCAP_ << "/"
0151         << minHitsinENDCAPplus_ << "/" << minHitsinENDCAPminus_;
0152 
0153     edm::LogInfo("AlignmentTrackSelector")
0154         << "Max value of |nHitsinENDCAPplus - nHitsinENDCAPminus| = " << maxHitDiffEndcaps_;
0155 
0156     if (!trkQualityStrings.empty()) {
0157       edm::LogInfo("AlignmentTrackSelector") << "Select tracks with these qualities: " << qualities;
0158     }
0159   }
0160 
0161   if (applyNHighestPt_)
0162     edm::LogInfo("AlignmentTrackSelector") << "filter N tracks with highest Pt N=" << nHighestPt_;
0163 
0164   if (applyMultiplicityFilter_)
0165     edm::LogInfo("AlignmentTrackSelector")
0166         << "apply multiplicity filter N>= " << minMultiplicity_ << "and N<= " << maxMultiplicity_ << " on "
0167         << (multiplicityOnInput_ ? "input" : "output");
0168 
0169   if (applyPrescaledHitsFilter_) {
0170     edm::LogInfo("AlignmentTrackSelector") << "apply cut on number of prescaled hits N>= " << minPrescaledHits_
0171                                            << " (prescale info from " << clusterValueMapTag_ << ")";
0172   }
0173 
0174   // Checking whether cuts on positions of first and last track hits are defined properly
0175   if (RorZofFirstHitMin_.size() != 2) {
0176     throw cms::Exception("BadConfig") << "@SUB=AlignmentTrackSelector::AlignmentTrackSelector"
0177                                       << "Wrong configuration of 'RorZofFirstHitMin'."
0178                                       << " Must have exactly 2 values instead of configured "
0179                                       << RorZofFirstHitMin_.size() << ")";
0180   } else {
0181     RorZofFirstHitMin_.at(0) = std::fabs(RorZofFirstHitMin_.at(0));
0182     RorZofFirstHitMin_.at(1) = std::fabs(RorZofFirstHitMin_.at(1));
0183   }
0184   if (RorZofFirstHitMax_.size() != 2) {
0185     throw cms::Exception("BadConfig") << "@SUB=AlignmentTrackSelector::AlignmentTrackSelector"
0186                                       << "Wrong configuration of 'RorZofFirstHitMax'."
0187                                       << " Must have exactly 2 values instead of configured "
0188                                       << RorZofFirstHitMax_.size() << ")";
0189   } else {
0190     RorZofFirstHitMax_.at(0) = std::fabs(RorZofFirstHitMax_.at(0));
0191     RorZofFirstHitMax_.at(1) = std::fabs(RorZofFirstHitMax_.at(1));
0192   }
0193   if (RorZofLastHitMin_.size() != 2) {
0194     throw cms::Exception("BadConfig") << "@SUB=AlignmentTrackSelector::AlignmentTrackSelector"
0195                                       << "Wrong configuration of 'RorZofLastHitMin'."
0196                                       << " Must have exactly 2 values instead of configured "
0197                                       << RorZofLastHitMin_.size() << ")";
0198   } else {
0199     RorZofLastHitMin_.at(0) = std::fabs(RorZofLastHitMin_.at(0));
0200     RorZofLastHitMin_.at(1) = std::fabs(RorZofLastHitMin_.at(1));
0201   }
0202   if (RorZofLastHitMax_.size() != 2) {
0203     throw cms::Exception("BadConfig") << "@SUB=AlignmentTrackSelector::AlignmentTrackSelector"
0204                                       << "Wrong configuration of 'RorZofLastHitMax'."
0205                                       << " Must have exactly 2 values instead of configured "
0206                                       << RorZofLastHitMax_.size() << ")";
0207   } else {
0208     RorZofLastHitMax_.at(0) = std::fabs(RorZofLastHitMax_.at(0));
0209     RorZofLastHitMax_.at(1) = std::fabs(RorZofLastHitMax_.at(1));
0210   }
0211   // If first hit set to be at larger distance then the last hit
0212   if (RorZofFirstHitMin_.at(0) > RorZofLastHitMax_.at(0) && RorZofFirstHitMin_.at(1) > RorZofLastHitMax_.at(1)) {
0213     throw cms::Exception("BadConfig") << "@SUB=AlignmentTrackSelector::AlignmentTrackSelector"
0214                                       << "Position of the first hit is set to larger distance than the last hit:."
0215                                       << " First hit(min): [" << RorZofFirstHitMin_.at(0) << ", "
0216                                       << RorZofFirstHitMin_.at(1) << "]; Last hit(max): [" << RorZofLastHitMax_.at(0)
0217                                       << ", " << RorZofLastHitMax_.at(1) << "];";
0218   }
0219 }
0220 
0221 void AlignmentTrackSelector::fillPSetDescription(edm::ParameterSetDescription& desc) {
0222   // Base TrackSelector settings
0223   desc.add<bool>("applyBasicCuts", true);
0224   desc.add<double>("ptMin", 0.0);
0225   desc.add<double>("ptMax", 999.0);
0226   desc.add<double>("pMin", 0.0);
0227   desc.add<double>("pMax", 9999.0);
0228   desc.add<double>("etaMin", -2.6);
0229   desc.add<double>("etaMax", 2.6);
0230   desc.add<double>("phiMax", 3.1416);
0231   desc.add<double>("phiMin", -3.1416);
0232   desc.add<double>("chi2nMax", 999999.0);
0233   desc.add<int>("theCharge", 0);  // -1: neg charge, +1: pos charge, 0: all charges
0234   desc.add<double>("d0Min", -999999.0);
0235   desc.add<double>("d0Max", 999999.0);
0236   desc.add<double>("dzMin", -999999.0);
0237   desc.add<double>("dzMax", 999999.0);
0238   desc.add<double>("nHitMin", 0.0);
0239   desc.add<double>("nHitMax", 999.0);
0240   desc.add<double>("nLostHitMax", 999.0);
0241   desc.add<unsigned int>("nHitMin2D", 0);
0242   desc.add<std::vector<double>>("RorZofFirstHitMin", {0.0, 0.0});
0243   desc.add<std::vector<double>>("RorZofFirstHitMax", {999.0, 999.0});
0244   desc.add<std::vector<double>>("RorZofLastHitMin", {0.0, 0.0});
0245   desc.add<std::vector<double>>("RorZofLastHitMax", {999.0, 999.0});
0246   desc.add<bool>("countStereoHitAs2D", true);
0247 
0248   // Nested PSet for minHitsPerSubDet
0249   edm::ParameterSetDescription minHitsPerSubDetDesc;
0250   minHitsPerSubDetDesc.add<int>("inTEC", 0);
0251   minHitsPerSubDetDesc.add<int>("inTOB", 0);
0252   minHitsPerSubDetDesc.add<int>("inFPIX", 0);
0253   minHitsPerSubDetDesc.add<int>("inTID", 0);
0254   minHitsPerSubDetDesc.add<int>("inBPIX", 0);
0255   minHitsPerSubDetDesc.add<int>("inTIB", 0);
0256   minHitsPerSubDetDesc.add<int>("inPIXEL", 0);
0257   minHitsPerSubDetDesc.add<int>("inTIDplus", 0);
0258   minHitsPerSubDetDesc.add<int>("inTIDminus", 0);
0259   minHitsPerSubDetDesc.add<int>("inTECplus", 0);
0260   minHitsPerSubDetDesc.add<int>("inTECminus", 0);
0261   minHitsPerSubDetDesc.add<int>("inFPIXplus", 0);
0262   minHitsPerSubDetDesc.add<int>("inFPIXminus", 0);
0263   minHitsPerSubDetDesc.add<int>("inENDCAP", 0);
0264   minHitsPerSubDetDesc.add<int>("inENDCAPplus", 0);
0265   minHitsPerSubDetDesc.add<int>("inENDCAPminus", 0);
0266   desc.add<edm::ParameterSetDescription>("minHitsPerSubDet", minHitsPerSubDetDesc);
0267 
0268   desc.add<double>("maxHitDiffEndcaps", 999.0);
0269   desc.add<int>("seedOnlyFrom", 0);
0270 
0271   // Multiplicity filtering
0272   desc.add<bool>("applyMultiplicityFilter", false);
0273   desc.add<int>("minMultiplicity", 1);
0274   desc.add<int>("maxMultiplicity", 999999);
0275   desc.add<bool>("multiplicityOnInput", false);
0276 
0277   // NHighestPt settings
0278   desc.add<bool>("applyNHighestPt", false);
0279   desc.add<int>("nHighestPt", 2);
0280 
0281   // Hit-related parameters
0282   desc.add<edm::InputTag>("rphirecHits", edm::InputTag("siStripMatchedRecHits", "rphiRecHit"));
0283   desc.add<edm::InputTag>("matchedrecHits", edm::InputTag("siStripMatchedRecHits", "matchedRecHit"));
0284   desc.add<bool>("applyIsolationCut", false);
0285   desc.add<double>("minHitIsolation", 0.01);
0286   desc.add<bool>("applyChargeCheck", false);
0287   desc.add<double>("minHitChargeStrip", 20.0);
0288 
0289   // String-based settings
0290   desc.add<std::vector<std::string>>("trackQualities", {});          // Take all if empty
0291   desc.add<std::vector<std::string>>("iterativeTrackingSteps", {});  // Take all if empty
0292 
0293   // Filtering on hits for Skim&Prescale workflow
0294   desc.add<edm::InputTag>("hitPrescaleMapTag", edm::InputTag(""));  // Ignore prescale map if empty
0295   desc.add<int>("minPrescaledHits", -1);
0296 }
0297 
0298 // destructor -----------------------------------------------------------------
0299 
0300 AlignmentTrackSelector::~AlignmentTrackSelector() = default;
0301 
0302 // do selection ---------------------------------------------------------------
0303 
0304 AlignmentTrackSelector::Tracks AlignmentTrackSelector::select(const Tracks& tracks,
0305                                                               const edm::Event& evt,
0306                                                               const edm::EventSetup& eSetup) const {
0307   if (applyMultiplicityFilter_ && multiplicityOnInput_ &&
0308       (tracks.size() < static_cast<unsigned int>(minMultiplicity_) ||
0309        tracks.size() > static_cast<unsigned int>(maxMultiplicity_))) {
0310     return Tracks();  // empty collection
0311   }
0312 
0313   Tracks result = tracks;
0314   // apply basic track cuts (if selected)
0315   if (applyBasicCuts_)
0316     result = this->basicCuts(result, evt, eSetup);
0317 
0318   // filter N tracks with highest Pt (if selected)
0319   if (applyNHighestPt_)
0320     result = this->theNHighestPtTracks(result);
0321 
0322   // apply minimum multiplicity requirement (if selected)
0323   if (applyMultiplicityFilter_ && !multiplicityOnInput_) {
0324     if (result.size() < static_cast<unsigned int>(minMultiplicity_) ||
0325         result.size() > static_cast<unsigned int>(maxMultiplicity_)) {
0326       result.clear();
0327     }
0328   }
0329 
0330   if (applyPrescaledHitsFilter_) {
0331     result = this->checkPrescaledHits(result, evt);
0332   }
0333 
0334   return result;
0335 }
0336 
0337 ///returns if any of the Filters is used.
0338 bool AlignmentTrackSelector::useThisFilter() {
0339   return applyMultiplicityFilter_ || applyBasicCuts_ || applyNHighestPt_ || applyPrescaledHitsFilter_;
0340 }
0341 
0342 // make basic cuts ------------------------------------------------------------
0343 
0344 AlignmentTrackSelector::Tracks AlignmentTrackSelector::basicCuts(const Tracks& tracks,
0345                                                                  const edm::Event& evt,
0346                                                                  const edm::EventSetup& eSetup) const {
0347   Tracks result;
0348 
0349   for (Tracks::const_iterator it = tracks.begin(); it != tracks.end(); ++it) {
0350     const reco::Track* trackp = *it;
0351     float pt = trackp->pt();
0352     float p = trackp->p();
0353     float eta = trackp->eta();
0354     float phi = trackp->phi();
0355     int nhit = trackp->numberOfValidHits();
0356     int nlosthit = trackp->numberOfLostHits();
0357     float chi2n = trackp->normalizedChi2();
0358 
0359     int q = trackp->charge();
0360     bool isChargeOk = false;
0361     if (theCharge_ == -1 && q < 0)
0362       isChargeOk = true;
0363     else if (theCharge_ == 1 && q > 0)
0364       isChargeOk = true;
0365     else if (theCharge_ == 0)
0366       isChargeOk = true;
0367 
0368     float d0 = trackp->d0();
0369     float dz = trackp->dz();
0370 
0371     // edm::LogDebug("AlignmentTrackSelector") << " pt,eta,phi,nhit: "
0372     //  <<pt<<","<<eta<<","<<phi<<","<<nhit;
0373 
0374     if (pt > ptMin_ && pt < ptMax_ && p > pMin_ && p < pMax_ && eta > etaMin_ && eta < etaMax_ && phi > phiMin_ &&
0375         phi < phiMax_ && nhit >= nHitMin_ && nhit <= nHitMax_ && nlosthit <= nLostHitMax_ && chi2n < chi2nMax_ &&
0376         isChargeOk && d0 >= d0Min_ && d0 <= d0Max_ && dz >= dzMin_ && dz <= dzMax_) {
0377       bool trkQualityOk = false;
0378       if (!applyTrkQualityCheck_ && !applyIterStepCheck_)
0379         trkQualityOk = true;  // nothing required
0380       else
0381         trkQualityOk = this->isOkTrkQuality(trackp);
0382 
0383       bool hitsCheckOk = this->detailedHitsCheck(trackp, evt, eSetup);
0384 
0385       if (trkQualityOk && hitsCheckOk)
0386         result.push_back(trackp);
0387     }
0388   }
0389 
0390   return result;
0391 }
0392 
0393 //-----------------------------------------------------------------------------
0394 
0395 bool AlignmentTrackSelector::detailedHitsCheck(const reco::Track* trackp,
0396                                                const edm::Event& evt,
0397                                                const edm::EventSetup& eSetup) const {
0398   //Retrieve tracker topology from geometry
0399   const TrackerTopology* const tTopo = &eSetup.getData(tTopoToken_);
0400 
0401   // checking hit requirements beyond simple number of valid hits
0402 
0403   if (minHitsinTIB_ || minHitsinTOB_ || minHitsinTID_ || minHitsinTEC_ || minHitsinENDCAP_ || minHitsinENDCAPplus_ ||
0404       minHitsinENDCAPminus_ || (maxHitDiffEndcaps_ < 999) || minHitsinTIDplus_ || minHitsinTIDminus_ ||
0405       minHitsinFPIXplus_ || minHitsinFPIXminus_ || minHitsinTECplus_ || minHitsinTECminus_ || minHitsinFPIX_ ||
0406       minHitsinBPIX_ || minHitsinPIX_ || nHitMin2D_ || chargeCheck_ || applyIsolation_ ||
0407       (seedOnlyFromAbove_ == 1 || seedOnlyFromAbove_ == 2) || !RorZofFirstHitMin_.empty() ||
0408       !RorZofFirstHitMax_.empty() || !RorZofLastHitMin_.empty() || !RorZofLastHitMax_.empty()) {
0409     // any detailed hit cut is active, so have to check
0410 
0411     int nhitinTIB = 0, nhitinTOB = 0, nhitinTID = 0;
0412     int nhitinTEC = 0, nhitinBPIX = 0, nhitinFPIX = 0, nhitinPIXEL = 0;
0413     int nhitinENDCAP = 0, nhitinENDCAPplus = 0, nhitinENDCAPminus = 0;
0414     int nhitinTIDplus = 0, nhitinTIDminus = 0;
0415     int nhitinFPIXplus = 0, nhitinFPIXminus = 0;
0416     int nhitinTECplus = 0, nhitinTECminus = 0;
0417     unsigned int nHit2D = 0;
0418     unsigned int thishit = 0;
0419 
0420     for (auto const& hit : trackp->recHits()) {
0421       thishit++;
0422       const DetId detId(hit->geographicalId());
0423       const int subdetId = detId.subdetId();
0424 
0425       // *** thishit == 1 means last hit in CTF ***
0426       // (FIXME: assumption might change or not be valid for all tracking algorthms)
0427       // ==> for cosmics
0428       // seedOnlyFrom = 1 is TIB-TOB-TEC tracks only
0429       // seedOnlyFrom = 2 is TOB-TEC tracks only
0430       if (seedOnlyFromAbove_ == 1 && thishit == 1 &&
0431           (subdetId == int(SiStripDetId::TOB) || subdetId == int(SiStripDetId::TEC))) {
0432         return false;
0433       }
0434       if (seedOnlyFromAbove_ == 2 && thishit == 1 && subdetId == int(SiStripDetId::TIB)) {
0435         return false;
0436       }
0437 
0438       if (!hit->isValid())
0439         continue;  // only real hits count as in trackp->numberOfValidHits()
0440       if (detId.det() != DetId::Tracker) {
0441         edm::LogError("DetectorMismatch")
0442             << "@SUB=AlignmentTrackSelector::detailedHitsCheck"
0443             << "DetId.det() != DetId::Tracker (=" << DetId::Tracker << "), but " << detId.det() << ".";
0444       }
0445       if (chargeCheck_ && !(this->isOkCharge(hit)))
0446         return false;
0447       if (applyIsolation_ && (!this->isIsolated(hit, evt)))
0448         return false;
0449       if (SiStripDetId::TIB == subdetId)
0450         ++nhitinTIB;
0451       else if (SiStripDetId::TOB == subdetId)
0452         ++nhitinTOB;
0453       else if (SiStripDetId::TID == subdetId) {
0454         ++nhitinTID;
0455         ++nhitinENDCAP;
0456 
0457         if (tTopo->tidIsZMinusSide(detId)) {
0458           ++nhitinTIDminus;
0459           ++nhitinENDCAPminus;
0460         } else if (tTopo->tidIsZPlusSide(detId)) {
0461           ++nhitinTIDplus;
0462           ++nhitinENDCAPplus;
0463         }
0464       } else if (SiStripDetId::TEC == subdetId) {
0465         ++nhitinTEC;
0466         ++nhitinENDCAP;
0467 
0468         if (tTopo->tecIsZMinusSide(detId)) {
0469           ++nhitinTECminus;
0470           ++nhitinENDCAPminus;
0471         } else if (tTopo->tecIsZPlusSide(detId)) {
0472           ++nhitinTECplus;
0473           ++nhitinENDCAPplus;
0474         }
0475       } else if (kBPIX == subdetId) {
0476         ++nhitinBPIX;
0477         ++nhitinPIXEL;
0478       } else if (kFPIX == subdetId) {
0479         ++nhitinFPIX;
0480         ++nhitinPIXEL;
0481 
0482         if (tTopo->pxfSide(detId) == 1)
0483           ++nhitinFPIXminus;
0484         else if (tTopo->pxfSide(detId) == 2)
0485           ++nhitinFPIXplus;
0486       }
0487       // Do not call isHit2D(..) if already enough 2D hits for performance reason:
0488       if (nHit2D < nHitMin2D_ && this->isHit2D(*hit))
0489         ++nHit2D;
0490     }  // end loop on hits
0491 
0492     // Checking whether the track satisfies requirement of the first and last hit positions
0493     bool passedLastHitPositionR = true;
0494     bool passedLastHitPositionZ = true;
0495     bool passedFirstHitPositionR = true;
0496     bool passedFirstHitPositionZ = true;
0497 
0498     if (RorZofFirstHitMin_.at(0) != 0.0 || RorZofFirstHitMin_.at(1) != 0.0 || RorZofFirstHitMax_.at(0) != 999.0 ||
0499         RorZofFirstHitMax_.at(1) != 999.0) {
0500       const reco::TrackBase::Point& firstPoint(trackp->innerPosition());
0501 
0502       if ((std::fabs(firstPoint.R()) < RorZofFirstHitMin_.at(0)))
0503         passedFirstHitPositionR = false;
0504       if ((std::fabs(firstPoint.R()) > RorZofFirstHitMax_.at(0)))
0505         passedFirstHitPositionR = false;
0506       if ((std::fabs(firstPoint.Z()) < RorZofFirstHitMin_.at(1)))
0507         passedFirstHitPositionZ = false;
0508       if ((std::fabs(firstPoint.Z()) > RorZofFirstHitMax_.at(1)))
0509         passedFirstHitPositionZ = false;
0510     }
0511 
0512     if (RorZofLastHitMin_.at(0) != 0.0 || RorZofLastHitMin_.at(1) != 0.0 || RorZofLastHitMax_.at(0) != 999.0 ||
0513         RorZofLastHitMax_.at(1) != 999.0) {
0514       const reco::TrackBase::Point& lastPoint(trackp->outerPosition());
0515 
0516       if ((std::fabs(lastPoint.R()) < RorZofLastHitMin_.at(0)))
0517         passedLastHitPositionR = false;
0518       if ((std::fabs(lastPoint.R()) > RorZofLastHitMax_.at(0)))
0519         passedLastHitPositionR = false;
0520       if ((std::fabs(lastPoint.Z()) < RorZofLastHitMin_.at(1)))
0521         passedLastHitPositionZ = false;
0522       if ((std::fabs(lastPoint.Z()) > RorZofLastHitMax_.at(1)))
0523         passedLastHitPositionZ = false;
0524     }
0525 
0526     bool passedFirstHitPosition = passedFirstHitPositionR || passedFirstHitPositionZ;
0527     bool passedLastHitPosition = passedLastHitPositionR || passedLastHitPositionZ;
0528 
0529     return (nhitinTIB >= minHitsinTIB_ && nhitinTOB >= minHitsinTOB_ && nhitinTID >= minHitsinTID_ &&
0530             nhitinTEC >= minHitsinTEC_ && nhitinENDCAP >= minHitsinENDCAP_ &&
0531             nhitinENDCAPplus >= minHitsinENDCAPplus_ && nhitinENDCAPminus >= minHitsinENDCAPminus_ &&
0532             std::abs(nhitinENDCAPplus - nhitinENDCAPminus) <= maxHitDiffEndcaps_ &&
0533             nhitinTIDplus >= minHitsinTIDplus_ && nhitinTIDminus >= minHitsinTIDminus_ &&
0534             nhitinFPIXplus >= minHitsinFPIXplus_ && nhitinFPIXminus >= minHitsinFPIXminus_ &&
0535             nhitinTECplus >= minHitsinTECplus_ && nhitinTECminus >= minHitsinTECminus_ &&
0536             nhitinBPIX >= minHitsinBPIX_ && nhitinFPIX >= minHitsinFPIX_ && nhitinPIXEL >= minHitsinPIX_ &&
0537             nHit2D >= nHitMin2D_ && passedFirstHitPosition && passedLastHitPosition);
0538   } else {  // no cuts set, so we are just fine and can avoid loop on hits
0539     return true;
0540   }
0541 }
0542 
0543 //-----------------------------------------------------------------------------
0544 
0545 bool AlignmentTrackSelector::isHit2D(const TrackingRecHit& hit) const {
0546   // we count SiStrip stereo modules as 2D if selected via countStereoHitAs2D_
0547   // (since they provide theta information)
0548   if (!hit.isValid() || (hit.dimension() < 2 && !countStereoHitAs2D_ && !dynamic_cast<const SiStripRecHit1D*>(&hit))) {
0549     return false;  // real RecHit1D - but SiStripRecHit1D depends on countStereoHitAs2D_
0550   } else {
0551     const DetId detId(hit.geographicalId());
0552     if (detId.det() == DetId::Tracker) {
0553       if (detId.subdetId() == kBPIX || detId.subdetId() == kFPIX) {
0554         return true;  // pixel is always 2D
0555       } else {        // should be SiStrip now
0556         const SiStripDetId stripId(detId);
0557         if (stripId.stereo())
0558           return countStereoHitAs2D_;  // stereo modules
0559         else if (dynamic_cast<const SiStripRecHit1D*>(&hit) || dynamic_cast<const SiStripRecHit2D*>(&hit))
0560           return false;  // rphi modules hit
0561                          //the following two are not used any more since ages...
0562         else if (dynamic_cast<const SiStripMatchedRecHit2D*>(&hit))
0563           return true;  // matched is 2D
0564         else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(&hit)) {
0565           const ProjectedSiStripRecHit2D* pH = static_cast<const ProjectedSiStripRecHit2D*>(&hit);
0566           return (countStereoHitAs2D_ && this->isHit2D(pH->originalHit()));  // depends on original...
0567         } else {
0568           edm::LogError("UnknownType") << "@SUB=AlignmentTrackSelector::isHit2D"
0569                                        << "Tracker hit not in pixel, neither SiStripRecHit[12]D nor "
0570                                        << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
0571           return false;
0572         }
0573       }
0574     } else {  // not tracker??
0575       edm::LogWarning("DetectorMismatch") << "@SUB=AlignmentTrackSelector::isHit2D"
0576                                           << "Hit not in tracker with 'official' dimension >=2.";
0577       return true;  // dimension() >= 2 so accept that...
0578     }
0579   }
0580   // never reached...
0581 }
0582 
0583 //-----------------------------------------------------------------------------
0584 
0585 bool AlignmentTrackSelector::isOkCharge(const TrackingRecHit* hit) const {
0586   if (!hit || !hit->isValid())
0587     return true;
0588 
0589   // check det and subdet
0590   const DetId id(hit->geographicalId());
0591   if (id.det() != DetId::Tracker) {
0592     edm::LogWarning("DetectorMismatch") << "@SUB=isOkCharge"
0593                                         << "Hit not in tracker!";
0594     return true;
0595   }
0596   if (id.subdetId() == kFPIX || id.subdetId() == kBPIX) {
0597     return true;  // might add some requirement...
0598   }
0599 
0600   // We are in SiStrip now, so test normal hit:
0601   const std::type_info& type = typeid(*hit);
0602 
0603   if (type == typeid(SiStripRecHit2D)) {
0604     const SiStripRecHit2D* stripHit2D = dynamic_cast<const SiStripRecHit2D*>(hit);
0605     if (stripHit2D) {
0606       return this->isOkChargeStripHit(*stripHit2D);
0607     }
0608   } else if (type == typeid(SiStripRecHit1D)) {
0609     const SiStripRecHit1D* stripHit1D = dynamic_cast<const SiStripRecHit1D*>(hit);
0610     if (stripHit1D) {
0611       return this->isOkChargeStripHit(*stripHit1D);
0612     }
0613   }
0614 
0615   else if (type ==
0616            typeid(SiStripMatchedRecHit2D)) {  // or matched (should not occur anymore due to hit splitting since 20X)
0617     const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
0618     if (matchedHit) {
0619       return (this->isOkChargeStripHit(matchedHit->monoHit()) && this->isOkChargeStripHit(matchedHit->stereoHit()));
0620     }
0621   } else if (type == typeid(ProjectedSiStripRecHit2D)) {
0622     // or projected (should not occur anymore due to hit splitting since 20X):
0623     const ProjectedSiStripRecHit2D* projHit = dynamic_cast<const ProjectedSiStripRecHit2D*>(hit);
0624     if (projHit) {
0625       return this->isOkChargeStripHit(projHit->originalHit());
0626     }
0627   } else {
0628     edm::LogError("AlignmentTrackSelector") << "@SUB=isOkCharge"
0629                                             << "Unknown type of a valid tracker hit in Strips "
0630                                             << " SubDet = " << id.subdetId();
0631     return false;
0632   }
0633 
0634   // and now? SiTrackerMultiRecHit? Not here I guess!
0635   // Should we throw instead?
0636   edm::LogError("AlignmentTrackSelector") << "@SUB=isOkCharge"
0637                                           << "Unknown valid tracker hit not in pixel, subdet " << id.subdetId()
0638                                           << ", SiTrackerMultiRecHit " << dynamic_cast<const SiTrackerMultiRecHit*>(hit)
0639                                           << ", BaseTrackerRecHit " << dynamic_cast<const BaseTrackerRecHit*>(hit);
0640 
0641   return true;
0642 }
0643 
0644 //-----------------------------------------------------------------------------
0645 
0646 bool AlignmentTrackSelector::isOkChargeStripHit(const SiStripRecHit2D& siStripRecHit2D) const {
0647   double charge = 0.;
0648 
0649   SiStripRecHit2D::ClusterRef cluster(siStripRecHit2D.cluster());
0650   const auto& amplitudes = cluster->amplitudes();
0651 
0652   for (size_t ia = 0; ia < amplitudes.size(); ++ia) {
0653     charge += amplitudes[ia];
0654   }
0655 
0656   return (charge >= minHitChargeStrip_);
0657 }
0658 
0659 //-----------------------------------------------------------------------------
0660 
0661 bool AlignmentTrackSelector::isOkChargeStripHit(const SiStripRecHit1D& siStripRecHit1D) const {
0662   double charge = 0.;
0663 
0664   SiStripRecHit1D::ClusterRef cluster(siStripRecHit1D.cluster());
0665   const auto& amplitudes = cluster->amplitudes();
0666 
0667   for (size_t ia = 0; ia < amplitudes.size(); ++ia) {
0668     charge += amplitudes[ia];
0669   }
0670 
0671   return (charge >= minHitChargeStrip_);
0672 }
0673 
0674 //-----------------------------------------------------------------------------
0675 
0676 bool AlignmentTrackSelector::isIsolated(const TrackingRecHit* hit, const edm::Event& evt) const {
0677   // FIXME:
0678   // adapt to changes after introduction of SiStripRecHit1D...
0679   //
0680   // edm::ESHandle<TrackerGeometry> tracker;
0681   edm::Handle<SiStripRecHit2DCollection> rphirecHits;
0682   edm::Handle<SiStripMatchedRecHit2DCollection> matchedrecHits;
0683   // es.get<TrackerDigiGeometryRecord>().get(tracker);
0684   evt.getByToken(rphirecHitsToken_, rphirecHits);
0685   evt.getByToken(matchedrecHitsToken_, matchedrecHits);
0686 
0687   SiStripRecHit2DCollection::DataContainer::const_iterator istripSt;
0688   SiStripMatchedRecHit2DCollection::DataContainer::const_iterator istripStm;
0689   const SiStripRecHit2DCollection::DataContainer& stripcollSt = rphirecHits->data();
0690   const SiStripMatchedRecHit2DCollection::DataContainer& stripcollStm = matchedrecHits->data();
0691 
0692   DetId idet = hit->geographicalId();
0693 
0694   // FIXME: instead of looping the full hit collection, we should explore the features of
0695   // SiStripRecHit2DCollection::rangeRphi = rphirecHits.get(idet) and loop
0696   // only from rangeRphi.first until rangeRphi.second
0697   for (istripSt = stripcollSt.begin(); istripSt != stripcollSt.end(); ++istripSt) {
0698     const SiStripRecHit2D* aHit = &*(istripSt);
0699     DetId mydet1 = aHit->geographicalId();
0700     if (idet.rawId() != mydet1.rawId())
0701       continue;
0702     float theDistance = (hit->localPosition() - aHit->localPosition()).mag();
0703     // std::cout << "theDistance1 = " << theDistance << "\n";
0704     if (theDistance > 0.001 && theDistance < minHitIsolation_)
0705       return false;
0706   }
0707 
0708   // FIXME: see above
0709   for (istripStm = stripcollStm.begin(); istripStm != stripcollStm.end(); ++istripStm) {
0710     const SiStripMatchedRecHit2D* aHit = &*(istripStm);
0711     DetId mydet2 = aHit->geographicalId();
0712     if (idet.rawId() != mydet2.rawId())
0713       continue;
0714     float theDistance = (hit->localPosition() - aHit->localPosition()).mag();
0715     // std::cout << "theDistance1 = " << theDistance << "\n";
0716     if (theDistance > 0.001 && theDistance < minHitIsolation_)
0717       return false;
0718   }
0719   return true;
0720 }
0721 
0722 //-----------------------------------------------------------------------------
0723 
0724 AlignmentTrackSelector::Tracks AlignmentTrackSelector::theNHighestPtTracks(const Tracks& tracks) const {
0725   Tracks sortedTracks = tracks;
0726   Tracks result;
0727 
0728   // sort in pt
0729   std::sort(sortedTracks.begin(), sortedTracks.end(), ptComparator);
0730 
0731   // copy theTrackMult highest pt tracks to result vector
0732   int n = 0;
0733   for (Tracks::const_iterator it = sortedTracks.begin(); it != sortedTracks.end(); ++it) {
0734     if (n < nHighestPt_) {
0735       result.push_back(*it);
0736       n++;
0737     }
0738   }
0739 
0740   return result;
0741 }
0742 
0743 //---------
0744 
0745 AlignmentTrackSelector::Tracks AlignmentTrackSelector::checkPrescaledHits(const Tracks& tracks,
0746                                                                           const edm::Event& evt) const {
0747   Tracks result;
0748 
0749   //take Cluster-Flag Assomap
0750   edm::Handle<AliClusterValueMap> fMap;
0751   evt.getByToken(clusterValueMapToken_, fMap);
0752   const AliClusterValueMap& flagMap = *fMap;
0753 
0754   //for each track loop on hits and count the number of taken hits
0755   for (auto const& trackp : tracks) {
0756     int ntakenhits = 0;
0757     //    float pt=trackp->pt();
0758 
0759     for (auto const& hit : trackp->recHits()) {
0760       if (!hit->isValid())
0761         continue;
0762       DetId detid = hit->geographicalId();
0763       int subDet = detid.subdetId();
0764       AlignmentClusterFlag flag;
0765 
0766       bool isPixelHit = (subDet == kFPIX || subDet == kBPIX);
0767 
0768       if (!isPixelHit) {
0769         const std::type_info& type = typeid(*hit);
0770 
0771         if (type == typeid(SiStripRecHit2D)) {
0772           const SiStripRecHit2D* striphit = dynamic_cast<const SiStripRecHit2D*>(hit);
0773           if (striphit != nullptr) {
0774             SiStripRecHit2D::ClusterRef stripclust(striphit->cluster());
0775             flag = flagMap[stripclust];
0776           }
0777         } else if (type == typeid(SiStripRecHit1D)) {
0778           const SiStripRecHit1D* striphit = dynamic_cast<const SiStripRecHit1D*>(hit);
0779           if (striphit != nullptr) {
0780             SiStripRecHit1D::ClusterRef stripclust(striphit->cluster());
0781             flag = flagMap[stripclust];
0782           }
0783         } else {
0784           edm::LogError("AlignmentTrackSelector")
0785               << "ERROR in <AlignmentTrackSelector::checkPrescaledHits>: Dynamic cast of Strip RecHit failed!"
0786               << " Skipping this hit.";
0787           continue;
0788         }
0789 
0790       }  //end if hit in Strips
0791       else {  // test explicitely BPIX/FPIX
0792         const SiPixelRecHit* pixelhit = dynamic_cast<const SiPixelRecHit*>(hit);
0793         if (pixelhit != nullptr) {
0794           SiPixelRecHit::ClusterRef pixclust(pixelhit->cluster());
0795           flag = flagMap[pixclust];
0796         } else {
0797           edm::LogError("AlignmentTrackSelector")
0798               << "ERROR in <AlignmentTrackSelector::checkPrescaledHits>: Dynamic cast of Pixel RecHit failed!  ";
0799         }
0800       }  //end else hit is in Pixel
0801 
0802       if (flag.isTaken())
0803         ntakenhits++;
0804 
0805     }  //end loop on hits
0806     if (ntakenhits >= minPrescaledHits_)
0807       result.push_back(trackp);
0808   }  //end loop on tracks
0809 
0810   return result;
0811 }  //end checkPrescaledHits
0812 
0813 //-----------------------------------------------------------------
0814 
0815 bool AlignmentTrackSelector::isOkTrkQuality(const reco::Track* track) const {
0816   bool qualityOk = false;
0817   bool iterStepOk = false;
0818 
0819   //check iterative step
0820   if (applyIterStepCheck_) {
0821     for (unsigned int i = 0; i < trkSteps_.size(); ++i) {
0822       if (track->algo() == (trkSteps_[i])) {
0823         iterStepOk = true;
0824       }
0825     }
0826   } else
0827     iterStepOk = true;
0828 
0829   //check track quality
0830   if (applyTrkQualityCheck_) {
0831     for (unsigned int i = 0; i < trkQualities_.size(); ++i) {
0832       if (track->quality(trkQualities_[i])) {
0833         qualityOk = true;
0834       }
0835     }
0836   } else
0837     qualityOk = true;
0838 
0839   return qualityOk && iterStepOk;
0840 }  //end check on track quality