Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:33:15

0001 #include "PhysicsTools/RecoUtils/interface/CheckHitPattern.h"
0002 
0003 // To get Tracker Geometry
0004 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0005 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0006 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0007 
0008 // To convert detId to subdet/layer number.
0009 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0010 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0011 //#include "DataFormats/SiStripDetId/interface/SiStripDetId.h"
0012 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0013 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0014 
0015 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0016 #include "DataFormats/Math/interface/deltaPhi.h"
0017 
0018 #include "FWCore/Utilities/interface/Exception.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 
0021 void CheckHitPattern::init(const edm::EventSetup& iSetup) {
0022   //Retrieve tracker topology from geometry
0023   edm::ESHandle<TrackerTopology> tTopoHandle;
0024   iSetup.get<TrackerTopologyRcd>().get(tTopoHandle);
0025   const TrackerTopology* const tTopo = tTopoHandle.product();
0026 
0027   iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder", trkTool_);  // Needed for vertex fits
0028 
0029   //
0030   // Note min/max radius (z) of each barrel layer (endcap disk).
0031   //
0032 
0033   geomInitDone_ = true;
0034 
0035   // Get Tracker geometry
0036   edm::ESHandle<TrackerGeometry> trackerGeometry;
0037   iSetup.get<TrackerDigiGeometryRecord>().get(trackerGeometry);
0038   const TrackingGeometry::DetContainer& dets = trackerGeometry->dets();
0039 
0040   // Loop over all modules in the Tracker.
0041   for (unsigned int i = 0; i < dets.size(); i++) {
0042     // Get subdet and layer of this module
0043     DetInfo detInfo = this->interpretDetId(dets[i]->geographicalId(), tTopo);
0044     uint32_t subDet = detInfo.first;
0045 
0046     // Note r (or z) of module if barrel (or endcap).
0047     double r_or_z;
0048     if (this->barrel(subDet)) {
0049       r_or_z = dets[i]->position().perp();
0050     } else {
0051       r_or_z = fabs(dets[i]->position().z());
0052     }
0053 
0054     // Recover min/max r/z value of this layer/disk found so far.
0055     double minRZ = 999.;
0056     double maxRZ = 0.;
0057     if (rangeRorZ_.find(detInfo) != rangeRorZ_.end()) {
0058       minRZ = rangeRorZ_[detInfo].first;
0059       maxRZ = rangeRorZ_[detInfo].second;
0060     }
0061 
0062     // Update with those of this module if it exceeds them.
0063     if (minRZ > r_or_z)
0064       minRZ = r_or_z;
0065     if (maxRZ < r_or_z)
0066       maxRZ = r_or_z;
0067     rangeRorZ_[detInfo] = std::pair<double, double>(minRZ, maxRZ);
0068   }
0069 
0070 #ifdef DEBUG_CHECKHITPATTERN
0071   RZrangeMap::const_iterator d;
0072   for (d = rangeRorZ_.begin(); d != rangeRorZ_.end(); d++) {
0073     DetInfo detInfo = d->first;
0074     std::pair<double, double> rangeRZ = d->second;
0075     std::std::cout << "CHECKHITPATTERN: Tracker subdetector type=" << detInfo.first << " layer=" << detInfo.second
0076                    << " has min r (or z) =" << rangeRZ.first << " and max r (or z) = " << rangeRZ.second
0077                    << std::std::endl;
0078   }
0079 #endif
0080 }
0081 
0082 CheckHitPattern::DetInfo CheckHitPattern::interpretDetId(DetId detId, const TrackerTopology* tTopo) {
0083   // Convert detId to a pair<uint32, uint32> consisting of the numbers used by HitPattern
0084   // to identify subdetector and layer number respectively.
0085   return DetInfo(detId.subdetId(), tTopo->layer(detId));
0086 }
0087 
0088 bool CheckHitPattern::barrel(uint32_t subDet) {
0089   // Determines if given sub-detector is in the barrel.
0090   return (subDet == StripSubdetector::TIB || subDet == StripSubdetector::TOB ||
0091           subDet == PixelSubdetector::PixelBarrel);
0092 }
0093 
0094 CheckHitPattern::Result CheckHitPattern::operator()(const reco::Track& track, const VertexState& vert) const {
0095   // Check if hit pattern of this track is consistent with it being produced
0096   // at given vertex.
0097 
0098   // Initialise geometry info if not yet done.
0099   if (!geomInitDone_)
0100     throw cms::Exception("CheckHitPattern::operator() called before CheckHitPattern::init");
0101 
0102   // Optionally set vertex position to zero for debugging.
0103   // VertexState vertDebug( GlobalPoint(0.,0.,0.) , GlobalError(1e-8, 0., 1e-8, 0., 0., 1e-8) );
0104 
0105   // Evaluate track parameters at vertex.
0106   reco::TransientTrack t_trk = trkTool_->build(track);
0107   GlobalVector p3_trk = t_trk.trajectoryStateClosestToPoint(vert.position()).momentum();
0108   bool trkGoesInsideOut =
0109       fabs(reco::deltaPhi<const GlobalVector, const GlobalPoint>(p3_trk, vert.position())) < 0.5 * M_PI;
0110 
0111   LogDebug("CHP") << "TRACK: in-->out ? " << trkGoesInsideOut << " dxy=" << track.dxy() << " sz=" << track.dz()
0112                   << " eta=" << track.eta() << " barrel hits=" << track.hitPattern().numberOfValidPixelHits() << "/"
0113                   << track.hitPattern().numberOfValidStripTIBHits() << "/"
0114                   << track.hitPattern().numberOfValidStripTOBHits();
0115   LogDebug("CHP") << "VERT: r=" << vert.position().perp() << " z=" << vert.position().z();
0116   //  if (vert.position().perp() < 3.5 && fabs(vert.position().z()) < 10. && fabs(track.eta()) < 1 && fabs(track.dxy()) < 2 && fabs(track.dz()) < 2 && track.hitPattern().numberOfValidPixelHits() == 0 && track.hitPattern().numberOfValidStripTIBHits() == 0) LogDebug("CHP")<<"LOOKATTHISTRACK";
0117   // Get hit patterns of this track
0118   const reco::HitPattern& hp = track.hitPattern();
0119 
0120   // Count number of valid hits on track definately in front of the vertex,
0121   // taking into account finite depth of each layer.
0122   unsigned int nHitBefore = 0;
0123   for (int i = 0; i < hp.numberOfAllHits(reco::HitPattern::TRACK_HITS); i++) {
0124     uint32_t hit = hp.getHitPattern(reco::HitPattern::TRACK_HITS, i);
0125     if (reco::HitPattern::trackerHitFilter(hit) && reco::HitPattern::validHitFilter(hit)) {
0126       uint32_t subDet = reco::HitPattern::getSubStructure(hit);
0127       uint32_t layer = reco::HitPattern::getLayer(hit);
0128       DetInfo detInfo(subDet, layer);
0129       auto maxRZ = (*rangeRorZ_.find(detInfo)).second.second;
0130 
0131       if (this->barrel(subDet)) {
0132         // Be careful. If the track starts by going outside-->in, it is allowed to have hits before the vertex !
0133         if (vert.position().perp() > maxRZ && trkGoesInsideOut)
0134           nHitBefore++;
0135       } else {
0136         if (fabs(vert.position().z()) > maxRZ)
0137           nHitBefore++;
0138       }
0139     }
0140   }
0141 
0142   // Count number of missing hits before the innermost hit on the track,
0143   // taking into account finite depth of each layer.
0144   unsigned int nMissHitAfter = 0;
0145   for (int i = 0; i < hp.numberOfAllHits(reco::HitPattern::MISSING_INNER_HITS); i++) {
0146     uint32_t hit = hp.getHitPattern(reco::HitPattern::MISSING_INNER_HITS, i);
0147     //    if (hp.trackerHitFilter(hit)) {
0148     if (reco::HitPattern::trackerHitFilter(hit) && reco::HitPattern::missingHitFilter(hit)) {
0149       uint32_t subDet = reco::HitPattern::getSubStructure(hit);
0150       uint32_t layer = reco::HitPattern::getLayer(hit);
0151       DetInfo detInfo(subDet, layer);
0152       auto minRZ = (*rangeRorZ_.find(detInfo)).second.first;
0153 
0154       if (this->barrel(subDet)) {
0155         // Be careful. If the track starts by going outside-->in, then it misses hits
0156         // in all layers it crosses  before its innermost valid hit.
0157         if (vert.position().perp() < minRZ || !trkGoesInsideOut)
0158           nMissHitAfter++;
0159       } else {
0160         if (fabs(vert.position().z()) < minRZ)
0161           nMissHitAfter++;
0162       }
0163     }
0164   }
0165 
0166   Result result;
0167   result.hitsInFrontOfVert = nHitBefore;
0168   result.missHitsAfterVert = nMissHitAfter;
0169   return result;
0170 }
0171 
0172 void CheckHitPattern::print(const reco::Track& track) {
0173   // Get hit patterns of this track
0174   const reco::HitPattern& hp = track.hitPattern();
0175   std::cout << "=== Hits on Track ===" << std::endl;
0176   print(reco::HitPattern::TRACK_HITS, hp);
0177   std::cout << "=== Hits before track ===" << std::endl;
0178   print(reco::HitPattern::MISSING_INNER_HITS, hp);
0179 }
0180 
0181 void CheckHitPattern::print(const reco::HitPattern::HitCategory category, const reco::HitPattern& hp) {
0182   for (int i = 0; i < hp.numberOfAllHits(category); i++) {
0183     uint32_t hit = hp.getHitPattern(category, i);
0184     if (reco::HitPattern::trackerHitFilter(hit)) {
0185       uint32_t subdet = reco::HitPattern::getSubStructure(hit);
0186       uint32_t layer = reco::HitPattern::getLayer(hit);
0187       std::cout << "hit " << i << " subdet=" << subdet << " layer=" << layer << " type " << hp.getHitType(hit)
0188                 << std::endl;
0189     }
0190   }
0191 }