RecoUtils/src/CheckHitPattern.cc

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