Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 10:43:52 

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