Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:56:17 

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 // destructor -----------------------------------------------------------------
0222 
0223 AlignmentTrackSelector::~AlignmentTrackSelector() = default;
0224 
0225 // do selection ---------------------------------------------------------------
0226 
0227 AlignmentTrackSelector::Tracks AlignmentTrackSelector::select(const Tracks& tracks,
0228                                                               const edm::Event& evt,
0229                                                               const edm::EventSetup& eSetup) const {
0230   if (applyMultiplicityFilter_ && multiplicityOnInput_ &&
0231       (tracks.size() < static_cast<unsigned int>(minMultiplicity_) ||
0232        tracks.size() > static_cast<unsigned int>(maxMultiplicity_))) {
0233     return Tracks();  // empty collection
0234   }
0235 
0236   Tracks result = tracks;
0237   // apply basic track cuts (if selected)
0238   if (applyBasicCuts_)
0239     result = this->basicCuts(result, evt, eSetup);
0240 
0241   // filter N tracks with highest Pt (if selected)
0242   if (applyNHighestPt_)
0243     result = this->theNHighestPtTracks(result);
0244 
0245   // apply minimum multiplicity requirement (if selected)
0246   if (applyMultiplicityFilter_ && !multiplicityOnInput_) {
0247     if (result.size() < static_cast<unsigned int>(minMultiplicity_) ||
0248         result.size() > static_cast<unsigned int>(maxMultiplicity_)) {
0249       result.clear();
0250     }
0251   }
0252 
0253   if (applyPrescaledHitsFilter_) {
0254     result = this->checkPrescaledHits(result, evt);
0255   }
0256 
0257   return result;
0258 }
0259 
0260 ///returns if any of the Filters is used.
0261 bool AlignmentTrackSelector::useThisFilter() {
0262   return applyMultiplicityFilter_ || applyBasicCuts_ || applyNHighestPt_ || applyPrescaledHitsFilter_;
0263 }
0264 
0265 // make basic cuts ------------------------------------------------------------
0266 
0267 AlignmentTrackSelector::Tracks AlignmentTrackSelector::basicCuts(const Tracks& tracks,
0268                                                                  const edm::Event& evt,
0269                                                                  const edm::EventSetup& eSetup) const {
0270   Tracks result;
0271 
0272   for (Tracks::const_iterator it = tracks.begin(); it != tracks.end(); ++it) {
0273     const reco::Track* trackp = *it;
0274     float pt = trackp->pt();
0275     float p = trackp->p();
0276     float eta = trackp->eta();
0277     float phi = trackp->phi();
0278     int nhit = trackp->numberOfValidHits();
0279     int nlosthit = trackp->numberOfLostHits();
0280     float chi2n = trackp->normalizedChi2();
0281 
0282     int q = trackp->charge();
0283     bool isChargeOk = false;
0284     if (theCharge_ == -1 && q < 0)
0285       isChargeOk = true;
0286     else if (theCharge_ == 1 && q > 0)
0287       isChargeOk = true;
0288     else if (theCharge_ == 0)
0289       isChargeOk = true;
0290 
0291     float d0 = trackp->d0();
0292     float dz = trackp->dz();
0293 
0294     // edm::LogDebug("AlignmentTrackSelector") << " pt,eta,phi,nhit: "
0295     //  <<pt<<","<<eta<<","<<phi<<","<<nhit;
0296 
0297     if (pt > ptMin_ && pt < ptMax_ && p > pMin_ && p < pMax_ && eta > etaMin_ && eta < etaMax_ && phi > phiMin_ &&
0298         phi < phiMax_ && nhit >= nHitMin_ && nhit <= nHitMax_ && nlosthit <= nLostHitMax_ && chi2n < chi2nMax_ &&
0299         isChargeOk && d0 >= d0Min_ && d0 <= d0Max_ && dz >= dzMin_ && dz <= dzMax_) {
0300       bool trkQualityOk = false;
0301       if (!applyTrkQualityCheck_ && !applyIterStepCheck_)
0302         trkQualityOk = true;  // nothing required
0303       else
0304         trkQualityOk = this->isOkTrkQuality(trackp);
0305 
0306       bool hitsCheckOk = this->detailedHitsCheck(trackp, evt, eSetup);
0307 
0308       if (trkQualityOk && hitsCheckOk)
0309         result.push_back(trackp);
0310     }
0311   }
0312 
0313   return result;
0314 }
0315 
0316 //-----------------------------------------------------------------------------
0317 
0318 bool AlignmentTrackSelector::detailedHitsCheck(const reco::Track* trackp,
0319                                                const edm::Event& evt,
0320                                                const edm::EventSetup& eSetup) const {
0321   //Retrieve tracker topology from geometry
0322   const TrackerTopology* const tTopo = &eSetup.getData(tTopoToken_);
0323 
0324   // checking hit requirements beyond simple number of valid hits
0325 
0326   if (minHitsinTIB_ || minHitsinTOB_ || minHitsinTID_ || minHitsinTEC_ || minHitsinENDCAP_ || minHitsinENDCAPplus_ ||
0327       minHitsinENDCAPminus_ || (maxHitDiffEndcaps_ < 999) || minHitsinTIDplus_ || minHitsinTIDminus_ ||
0328       minHitsinFPIXplus_ || minHitsinFPIXminus_ || minHitsinTECplus_ || minHitsinTECminus_ || minHitsinFPIX_ ||
0329       minHitsinBPIX_ || minHitsinPIX_ || nHitMin2D_ || chargeCheck_ || applyIsolation_ ||
0330       (seedOnlyFromAbove_ == 1 || seedOnlyFromAbove_ == 2) || !RorZofFirstHitMin_.empty() ||
0331       !RorZofFirstHitMax_.empty() || !RorZofLastHitMin_.empty() || !RorZofLastHitMax_.empty()) {
0332     // any detailed hit cut is active, so have to check
0333 
0334     int nhitinTIB = 0, nhitinTOB = 0, nhitinTID = 0;
0335     int nhitinTEC = 0, nhitinBPIX = 0, nhitinFPIX = 0, nhitinPIXEL = 0;
0336     int nhitinENDCAP = 0, nhitinENDCAPplus = 0, nhitinENDCAPminus = 0;
0337     int nhitinTIDplus = 0, nhitinTIDminus = 0;
0338     int nhitinFPIXplus = 0, nhitinFPIXminus = 0;
0339     int nhitinTECplus = 0, nhitinTECminus = 0;
0340     unsigned int nHit2D = 0;
0341     unsigned int thishit = 0;
0342 
0343     for (auto const& hit : trackp->recHits()) {
0344       thishit++;
0345       const DetId detId(hit->geographicalId());
0346       const int subdetId = detId.subdetId();
0347 
0348       // *** thishit == 1 means last hit in CTF ***
0349       // (FIXME: assumption might change or not be valid for all tracking algorthms)
0350       // ==> for cosmics
0351       // seedOnlyFrom = 1 is TIB-TOB-TEC tracks only
0352       // seedOnlyFrom = 2 is TOB-TEC tracks only
0353       if (seedOnlyFromAbove_ == 1 && thishit == 1 &&
0354           (subdetId == int(SiStripDetId::TOB) || subdetId == int(SiStripDetId::TEC))) {
0355         return false;
0356       }
0357       if (seedOnlyFromAbove_ == 2 && thishit == 1 && subdetId == int(SiStripDetId::TIB)) {
0358         return false;
0359       }
0360 
0361       if (!hit->isValid())
0362         continue;  // only real hits count as in trackp->numberOfValidHits()
0363       if (detId.det() != DetId::Tracker) {
0364         edm::LogError("DetectorMismatch")
0365             << "@SUB=AlignmentTrackSelector::detailedHitsCheck"
0366             << "DetId.det() != DetId::Tracker (=" << DetId::Tracker << "), but " << detId.det() << ".";
0367       }
0368       if (chargeCheck_ && !(this->isOkCharge(hit)))
0369         return false;
0370       if (applyIsolation_ && (!this->isIsolated(hit, evt)))
0371         return false;
0372       if (SiStripDetId::TIB == subdetId)
0373         ++nhitinTIB;
0374       else if (SiStripDetId::TOB == subdetId)
0375         ++nhitinTOB;
0376       else if (SiStripDetId::TID == subdetId) {
0377         ++nhitinTID;
0378         ++nhitinENDCAP;
0379 
0380         if (tTopo->tidIsZMinusSide(detId)) {
0381           ++nhitinTIDminus;
0382           ++nhitinENDCAPminus;
0383         } else if (tTopo->tidIsZPlusSide(detId)) {
0384           ++nhitinTIDplus;
0385           ++nhitinENDCAPplus;
0386         }
0387       } else if (SiStripDetId::TEC == subdetId) {
0388         ++nhitinTEC;
0389         ++nhitinENDCAP;
0390 
0391         if (tTopo->tecIsZMinusSide(detId)) {
0392           ++nhitinTECminus;
0393           ++nhitinENDCAPminus;
0394         } else if (tTopo->tecIsZPlusSide(detId)) {
0395           ++nhitinTECplus;
0396           ++nhitinENDCAPplus;
0397         }
0398       } else if (kBPIX == subdetId) {
0399         ++nhitinBPIX;
0400         ++nhitinPIXEL;
0401       } else if (kFPIX == subdetId) {
0402         ++nhitinFPIX;
0403         ++nhitinPIXEL;
0404 
0405         if (tTopo->pxfSide(detId) == 1)
0406           ++nhitinFPIXminus;
0407         else if (tTopo->pxfSide(detId) == 2)
0408           ++nhitinFPIXplus;
0409       }
0410       // Do not call isHit2D(..) if already enough 2D hits for performance reason:
0411       if (nHit2D < nHitMin2D_ && this->isHit2D(*hit))
0412         ++nHit2D;
0413     }  // end loop on hits
0414 
0415     // Checking whether the track satisfies requirement of the first and last hit positions
0416     bool passedLastHitPositionR = true;
0417     bool passedLastHitPositionZ = true;
0418     bool passedFirstHitPositionR = true;
0419     bool passedFirstHitPositionZ = true;
0420 
0421     if (RorZofFirstHitMin_.at(0) != 0.0 || RorZofFirstHitMin_.at(1) != 0.0 || RorZofFirstHitMax_.at(0) != 999.0 ||
0422         RorZofFirstHitMax_.at(1) != 999.0) {
0423       const reco::TrackBase::Point& firstPoint(trackp->innerPosition());
0424 
0425       if ((std::fabs(firstPoint.R()) < RorZofFirstHitMin_.at(0)))
0426         passedFirstHitPositionR = false;
0427       if ((std::fabs(firstPoint.R()) > RorZofFirstHitMax_.at(0)))
0428         passedFirstHitPositionR = false;
0429       if ((std::fabs(firstPoint.Z()) < RorZofFirstHitMin_.at(1)))
0430         passedFirstHitPositionZ = false;
0431       if ((std::fabs(firstPoint.Z()) > RorZofFirstHitMax_.at(1)))
0432         passedFirstHitPositionZ = false;
0433     }
0434 
0435     if (RorZofLastHitMin_.at(0) != 0.0 || RorZofLastHitMin_.at(1) != 0.0 || RorZofLastHitMax_.at(0) != 999.0 ||
0436         RorZofLastHitMax_.at(1) != 999.0) {
0437       const reco::TrackBase::Point& lastPoint(trackp->outerPosition());
0438 
0439       if ((std::fabs(lastPoint.R()) < RorZofLastHitMin_.at(0)))
0440         passedLastHitPositionR = false;
0441       if ((std::fabs(lastPoint.R()) > RorZofLastHitMax_.at(0)))
0442         passedLastHitPositionR = false;
0443       if ((std::fabs(lastPoint.Z()) < RorZofLastHitMin_.at(1)))
0444         passedLastHitPositionZ = false;
0445       if ((std::fabs(lastPoint.Z()) > RorZofLastHitMax_.at(1)))
0446         passedLastHitPositionZ = false;
0447     }
0448 
0449     bool passedFirstHitPosition = passedFirstHitPositionR || passedFirstHitPositionZ;
0450     bool passedLastHitPosition = passedLastHitPositionR || passedLastHitPositionZ;
0451 
0452     return (nhitinTIB >= minHitsinTIB_ && nhitinTOB >= minHitsinTOB_ && nhitinTID >= minHitsinTID_ &&
0453             nhitinTEC >= minHitsinTEC_ && nhitinENDCAP >= minHitsinENDCAP_ &&
0454             nhitinENDCAPplus >= minHitsinENDCAPplus_ && nhitinENDCAPminus >= minHitsinENDCAPminus_ &&
0455             std::abs(nhitinENDCAPplus - nhitinENDCAPminus) <= maxHitDiffEndcaps_ &&
0456             nhitinTIDplus >= minHitsinTIDplus_ && nhitinTIDminus >= minHitsinTIDminus_ &&
0457             nhitinFPIXplus >= minHitsinFPIXplus_ && nhitinFPIXminus >= minHitsinFPIXminus_ &&
0458             nhitinTECplus >= minHitsinTECplus_ && nhitinTECminus >= minHitsinTECminus_ &&
0459             nhitinBPIX >= minHitsinBPIX_ && nhitinFPIX >= minHitsinFPIX_ && nhitinPIXEL >= minHitsinPIX_ &&
0460             nHit2D >= nHitMin2D_ && passedFirstHitPosition && passedLastHitPosition);
0461   } else {  // no cuts set, so we are just fine and can avoid loop on hits
0462     return true;
0463   }
0464 }
0465 
0466 //-----------------------------------------------------------------------------
0467 
0468 bool AlignmentTrackSelector::isHit2D(const TrackingRecHit& hit) const {
0469   // we count SiStrip stereo modules as 2D if selected via countStereoHitAs2D_
0470   // (since they provide theta information)
0471   if (!hit.isValid() || (hit.dimension() < 2 && !countStereoHitAs2D_ && !dynamic_cast<const SiStripRecHit1D*>(&hit))) {
0472     return false;  // real RecHit1D - but SiStripRecHit1D depends on countStereoHitAs2D_
0473   } else {
0474     const DetId detId(hit.geographicalId());
0475     if (detId.det() == DetId::Tracker) {
0476       if (detId.subdetId() == kBPIX || detId.subdetId() == kFPIX) {
0477         return true;  // pixel is always 2D
0478       } else {        // should be SiStrip now
0479         const SiStripDetId stripId(detId);
0480         if (stripId.stereo())
0481           return countStereoHitAs2D_;  // stereo modules
0482         else if (dynamic_cast<const SiStripRecHit1D*>(&hit) || dynamic_cast<const SiStripRecHit2D*>(&hit))
0483           return false;  // rphi modules hit
0484                          //the following two are not used any more since ages...
0485         else if (dynamic_cast<const SiStripMatchedRecHit2D*>(&hit))
0486           return true;  // matched is 2D
0487         else if (dynamic_cast<const ProjectedSiStripRecHit2D*>(&hit)) {
0488           const ProjectedSiStripRecHit2D* pH = static_cast<const ProjectedSiStripRecHit2D*>(&hit);
0489           return (countStereoHitAs2D_ && this->isHit2D(pH->originalHit()));  // depends on original...
0490         } else {
0491           edm::LogError("UnknownType") << "@SUB=AlignmentTrackSelector::isHit2D"
0492                                        << "Tracker hit not in pixel, neither SiStripRecHit[12]D nor "
0493                                        << "SiStripMatchedRecHit2D nor ProjectedSiStripRecHit2D.";
0494           return false;
0495         }
0496       }
0497     } else {  // not tracker??
0498       edm::LogWarning("DetectorMismatch") << "@SUB=AlignmentTrackSelector::isHit2D"
0499                                           << "Hit not in tracker with 'official' dimension >=2.";
0500       return true;  // dimension() >= 2 so accept that...
0501     }
0502   }
0503   // never reached...
0504 }
0505 
0506 //-----------------------------------------------------------------------------
0507 
0508 bool AlignmentTrackSelector::isOkCharge(const TrackingRecHit* hit) const {
0509   if (!hit || !hit->isValid())
0510     return true;
0511 
0512   // check det and subdet
0513   const DetId id(hit->geographicalId());
0514   if (id.det() != DetId::Tracker) {
0515     edm::LogWarning("DetectorMismatch") << "@SUB=isOkCharge"
0516                                         << "Hit not in tracker!";
0517     return true;
0518   }
0519   if (id.subdetId() == kFPIX || id.subdetId() == kBPIX) {
0520     return true;  // might add some requirement...
0521   }
0522 
0523   // We are in SiStrip now, so test normal hit:
0524   const std::type_info& type = typeid(*hit);
0525 
0526   if (type == typeid(SiStripRecHit2D)) {
0527     const SiStripRecHit2D* stripHit2D = dynamic_cast<const SiStripRecHit2D*>(hit);
0528     if (stripHit2D) {
0529       return this->isOkChargeStripHit(*stripHit2D);
0530     }
0531   } else if (type == typeid(SiStripRecHit1D)) {
0532     const SiStripRecHit1D* stripHit1D = dynamic_cast<const SiStripRecHit1D*>(hit);
0533     if (stripHit1D) {
0534       return this->isOkChargeStripHit(*stripHit1D);
0535     }
0536   }
0537 
0538   else if (type ==
0539            typeid(SiStripMatchedRecHit2D)) {  // or matched (should not occur anymore due to hit splitting since 20X)
0540     const SiStripMatchedRecHit2D* matchedHit = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
0541     if (matchedHit) {
0542       return (this->isOkChargeStripHit(matchedHit->monoHit()) && this->isOkChargeStripHit(matchedHit->stereoHit()));
0543     }
0544   } else if (type == typeid(ProjectedSiStripRecHit2D)) {
0545     // or projected (should not occur anymore due to hit splitting since 20X):
0546     const ProjectedSiStripRecHit2D* projHit = dynamic_cast<const ProjectedSiStripRecHit2D*>(hit);
0547     if (projHit) {
0548       return this->isOkChargeStripHit(projHit->originalHit());
0549     }
0550   } else {
0551     edm::LogError("AlignmentTrackSelector") << "@SUB=isOkCharge"
0552                                             << "Unknown type of a valid tracker hit in Strips "
0553                                             << " SubDet = " << id.subdetId();
0554     return false;
0555   }
0556 
0557   // and now? SiTrackerMultiRecHit? Not here I guess!
0558   // Should we throw instead?
0559   edm::LogError("AlignmentTrackSelector") << "@SUB=isOkCharge"
0560                                           << "Unknown valid tracker hit not in pixel, subdet " << id.subdetId()
0561                                           << ", SiTrackerMultiRecHit " << dynamic_cast<const SiTrackerMultiRecHit*>(hit)
0562                                           << ", BaseTrackerRecHit " << dynamic_cast<const BaseTrackerRecHit*>(hit);
0563 
0564   return true;
0565 }
0566 
0567 //-----------------------------------------------------------------------------
0568 
0569 bool AlignmentTrackSelector::isOkChargeStripHit(const SiStripRecHit2D& siStripRecHit2D) const {
0570   double charge = 0.;
0571 
0572   SiStripRecHit2D::ClusterRef cluster(siStripRecHit2D.cluster());
0573   const auto& amplitudes = cluster->amplitudes();
0574 
0575   for (size_t ia = 0; ia < amplitudes.size(); ++ia) {
0576     charge += amplitudes[ia];
0577   }
0578 
0579   return (charge >= minHitChargeStrip_);
0580 }
0581 
0582 //-----------------------------------------------------------------------------
0583 
0584 bool AlignmentTrackSelector::isOkChargeStripHit(const SiStripRecHit1D& siStripRecHit1D) const {
0585   double charge = 0.;
0586 
0587   SiStripRecHit1D::ClusterRef cluster(siStripRecHit1D.cluster());
0588   const auto& amplitudes = cluster->amplitudes();
0589 
0590   for (size_t ia = 0; ia < amplitudes.size(); ++ia) {
0591     charge += amplitudes[ia];
0592   }
0593 
0594   return (charge >= minHitChargeStrip_);
0595 }
0596 
0597 //-----------------------------------------------------------------------------
0598 
0599 bool AlignmentTrackSelector::isIsolated(const TrackingRecHit* hit, const edm::Event& evt) const {
0600   // FIXME:
0601   // adapt to changes after introduction of SiStripRecHit1D...
0602   //
0603   // edm::ESHandle<TrackerGeometry> tracker;
0604   edm::Handle<SiStripRecHit2DCollection> rphirecHits;
0605   edm::Handle<SiStripMatchedRecHit2DCollection> matchedrecHits;
0606   // es.get<TrackerDigiGeometryRecord>().get(tracker);
0607   evt.getByToken(rphirecHitsToken_, rphirecHits);
0608   evt.getByToken(matchedrecHitsToken_, matchedrecHits);
0609 
0610   SiStripRecHit2DCollection::DataContainer::const_iterator istripSt;
0611   SiStripMatchedRecHit2DCollection::DataContainer::const_iterator istripStm;
0612   const SiStripRecHit2DCollection::DataContainer& stripcollSt = rphirecHits->data();
0613   const SiStripMatchedRecHit2DCollection::DataContainer& stripcollStm = matchedrecHits->data();
0614 
0615   DetId idet = hit->geographicalId();
0616 
0617   // FIXME: instead of looping the full hit collection, we should explore the features of
0618   // SiStripRecHit2DCollection::rangeRphi = rphirecHits.get(idet) and loop
0619   // only from rangeRphi.first until rangeRphi.second
0620   for (istripSt = stripcollSt.begin(); istripSt != stripcollSt.end(); ++istripSt) {
0621     const SiStripRecHit2D* aHit = &*(istripSt);
0622     DetId mydet1 = aHit->geographicalId();
0623     if (idet.rawId() != mydet1.rawId())
0624       continue;
0625     float theDistance = (hit->localPosition() - aHit->localPosition()).mag();
0626     // std::cout << "theDistance1 = " << theDistance << "\n";
0627     if (theDistance > 0.001 && theDistance < minHitIsolation_)
0628       return false;
0629   }
0630 
0631   // FIXME: see above
0632   for (istripStm = stripcollStm.begin(); istripStm != stripcollStm.end(); ++istripStm) {
0633     const SiStripMatchedRecHit2D* aHit = &*(istripStm);
0634     DetId mydet2 = aHit->geographicalId();
0635     if (idet.rawId() != mydet2.rawId())
0636       continue;
0637     float theDistance = (hit->localPosition() - aHit->localPosition()).mag();
0638     // std::cout << "theDistance1 = " << theDistance << "\n";
0639     if (theDistance > 0.001 && theDistance < minHitIsolation_)
0640       return false;
0641   }
0642   return true;
0643 }
0644 
0645 //-----------------------------------------------------------------------------
0646 
0647 AlignmentTrackSelector::Tracks AlignmentTrackSelector::theNHighestPtTracks(const Tracks& tracks) const {
0648   Tracks sortedTracks = tracks;
0649   Tracks result;
0650 
0651   // sort in pt
0652   std::sort(sortedTracks.begin(), sortedTracks.end(), ptComparator);
0653 
0654   // copy theTrackMult highest pt tracks to result vector
0655   int n = 0;
0656   for (Tracks::const_iterator it = sortedTracks.begin(); it != sortedTracks.end(); ++it) {
0657     if (n < nHighestPt_) {
0658       result.push_back(*it);
0659       n++;
0660     }
0661   }
0662 
0663   return result;
0664 }
0665 
0666 //---------
0667 
0668 AlignmentTrackSelector::Tracks AlignmentTrackSelector::checkPrescaledHits(const Tracks& tracks,
0669                                                                           const edm::Event& evt) const {
0670   Tracks result;
0671 
0672   //take Cluster-Flag Assomap
0673   edm::Handle<AliClusterValueMap> fMap;
0674   evt.getByToken(clusterValueMapToken_, fMap);
0675   const AliClusterValueMap& flagMap = *fMap;
0676 
0677   //for each track loop on hits and count the number of taken hits
0678   for (auto const& trackp : tracks) {
0679     int ntakenhits = 0;
0680     //    float pt=trackp->pt();
0681 
0682     for (auto const& hit : trackp->recHits()) {
0683       if (!hit->isValid())
0684         continue;
0685       DetId detid = hit->geographicalId();
0686       int subDet = detid.subdetId();
0687       AlignmentClusterFlag flag;
0688 
0689       bool isPixelHit = (subDet == kFPIX || subDet == kBPIX);
0690 
0691       if (!isPixelHit) {
0692         const std::type_info& type = typeid(*hit);
0693 
0694         if (type == typeid(SiStripRecHit2D)) {
0695           const SiStripRecHit2D* striphit = dynamic_cast<const SiStripRecHit2D*>(hit);
0696           if (striphit != nullptr) {
0697             SiStripRecHit2D::ClusterRef stripclust(striphit->cluster());
0698             flag = flagMap[stripclust];
0699           }
0700         } else if (type == typeid(SiStripRecHit1D)) {
0701           const SiStripRecHit1D* striphit = dynamic_cast<const SiStripRecHit1D*>(hit);
0702           if (striphit != nullptr) {
0703             SiStripRecHit1D::ClusterRef stripclust(striphit->cluster());
0704             flag = flagMap[stripclust];
0705           }
0706         } else {
0707           edm::LogError("AlignmentTrackSelector")
0708               << "ERROR in <AlignmentTrackSelector::checkPrescaledHits>: Dynamic cast of Strip RecHit failed!"
0709               << " Skipping this hit.";
0710           continue;
0711         }
0712 
0713       }       //end if hit in Strips
0714       else {  // test explicitely BPIX/FPIX
0715         const SiPixelRecHit* pixelhit = dynamic_cast<const SiPixelRecHit*>(hit);
0716         if (pixelhit != nullptr) {
0717           SiPixelRecHit::ClusterRef pixclust(pixelhit->cluster());
0718           flag = flagMap[pixclust];
0719         } else {
0720           edm::LogError("AlignmentTrackSelector")
0721               << "ERROR in <AlignmentTrackSelector::checkPrescaledHits>: Dynamic cast of Pixel RecHit failed!  ";
0722         }
0723       }  //end else hit is in Pixel
0724 
0725       if (flag.isTaken())
0726         ntakenhits++;
0727 
0728     }  //end loop on hits
0729     if (ntakenhits >= minPrescaledHits_)
0730       result.push_back(trackp);
0731   }  //end loop on tracks
0732 
0733   return result;
0734 }  //end checkPrescaledHits
0735 
0736 //-----------------------------------------------------------------
0737 
0738 bool AlignmentTrackSelector::isOkTrkQuality(const reco::Track* track) const {
0739   bool qualityOk = false;
0740   bool iterStepOk = false;
0741 
0742   //check iterative step
0743   if (applyIterStepCheck_) {
0744     for (unsigned int i = 0; i < trkSteps_.size(); ++i) {
0745       if (track->algo() == (trkSteps_[i])) {
0746         iterStepOk = true;
0747       }
0748     }
0749   } else
0750     iterStepOk = true;
0751 
0752   //check track quality
0753   if (applyTrkQualityCheck_) {
0754     for (unsigned int i = 0; i < trkQualities_.size(); ++i) {
0755       if (track->quality(trkQualities_[i])) {
0756         qualityOk = true;
0757       }
0758     }
0759   } else
0760     qualityOk = true;
0761 
0762   return qualityOk && iterStepOk;
0763 }  //end check on track quality
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!