PFTracking/src/PFCheckHitPattern.cc

0001 #include "RecoParticleFlow/PFTracking/interface/PFCheckHitPattern.h"
0002
0003 // To get Tracker Geometry
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 #include "FWCore/Utilities/interface/Exception.h"
0006 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0007 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0008 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0009
0010 // To convert detId to subdet/layer number.
0011 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0012 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0013 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0014
0015 #include <map>
0016
0017 using namespace reco;
0018 using namespace std;
0019
0020 void PFCheckHitPattern::init(const TrackerTopology* tkerTopo, const TrackerGeometry* tkerGeom) {
0021   //
0022   // Note min/max radius (z) of each barrel layer (endcap disk).
0023   //
0024
0025   geomInitDone_ = true;
0026
0027   // Get Tracker geometry
0028   const TrackingGeometry::DetContainer& dets = tkerGeom->dets();
0029
0030   // Loop over all modules in the Tracker.
0031   for (unsigned int i = 0; i < dets.size(); i++) {
0032     // Get subdet and layer of this module
0033     auto detId = dets[i]->geographicalId();
0034     auto detInfo = DetInfo(detId.subdetId(), tkerTopo->layer(detId));
0035     uint32_t subDet = detInfo.first;
0036
0037     // Note r (or z) of module if barrel (or endcap).
0038     double r_or_z;
0039     if (this->barrel(subDet)) {
0040       r_or_z = dets[i]->position().perp();
0041     } else {
0042       r_or_z = fabs(dets[i]->position().z());
0043     }
0044
0045     // Recover min/max r/z value of this layer/disk found so far.
0046     double minRZ = 999.;
0047     double maxRZ = 0.;
0048     if (rangeRorZ_.find(detInfo) != rangeRorZ_.end()) {
0049       minRZ = rangeRorZ_[detInfo].first;
0050       maxRZ = rangeRorZ_[detInfo].second;
0051     }
0052
0053     // Update with those of this module if it exceeds them.
0054     if (minRZ > r_or_z)
0055       minRZ = r_or_z;
0056     if (maxRZ < r_or_z)
0057       maxRZ = r_or_z;
0058     rangeRorZ_[detInfo] = pair<double, double>(minRZ, maxRZ);
0059   }
0060
0061 #if 0
0062   //def DEBUG_CHECKHITPATTERN
0063   RZrangeMap::const_iterator d;
0064   for (d = rangeRorZ_.begin(); d != rangeRorZ_.end(); d++) {
0065     DetInfo detInfo = d->first;
0066     pair<double, double> rangeRZ = d->second;
0067   }
0068 #endif
0069 }
0070
0071 bool PFCheckHitPattern::barrel(uint32_t subDet) {
0072   // Determines if given sub-detector is in the barrel.
0073   return (subDet == StripSubdetector::TIB || subDet == StripSubdetector::TOB ||
0074           subDet == PixelSubdetector::PixelBarrel);
0075 }
0076
0077 pair<PFCheckHitPattern::PFTrackHitInfo, PFCheckHitPattern::PFTrackHitInfo> PFCheckHitPattern::analyze(
0078     const TrackerTopology* tkerTopo,
0079     const TrackerGeometry* tkerGeom,
0080     const TrackBaseRef track,
0081     const TransientVertex& vert) {
0082   // PFCheck if hit pattern of this track is consistent with it being produced
0083   // at given vertex. Pair.first gives number of hits on track in front of vertex.
0084   // Pair.second gives number of missing hits between vertex and innermost hit
0085   // on track.
0086
0087   // Initialise geometry info if not yet done.
0088   if (!geomInitDone_)
0089     this->init(tkerTopo, tkerGeom);
0090
0091   // Get hit patterns of this track
0092   const reco::HitPattern& hp = track.get()->hitPattern();
0093
0094   // Count number of valid hits on track definately in front of the vertex,
0095   // taking into account finite depth of each layer.
0096   unsigned int nHitBefore = 0;
0097   unsigned int nHitAfter = 0;
0098
0099   for (int i = 0; i < hp.numberOfAllHits(HitPattern::TRACK_HITS); i++) {
0100     uint32_t hit = hp.getHitPattern(HitPattern::TRACK_HITS, i);
0101     if (hp.trackerHitFilter(hit) && hp.validHitFilter(hit)) {
0102       uint32_t subDet = hp.getSubStructure(hit);
0103       uint32_t layer = hp.getLayer(hit);
0104       DetInfo detInfo(subDet, layer);
0105       double maxRZ = rangeRorZ_[detInfo].second;
0106
0107       if (this->barrel(subDet)) {
0108         if (vert.position().perp() > maxRZ)
0109           nHitBefore++;
0110         else
0111           nHitAfter++;
0112       } else {
0113         if (fabs(vert.position().z()) > maxRZ)
0114           nHitBefore++;
0115         else
0116           nHitAfter++;
0117       }
0118     }
0119   }
0120
0121   // Count number of missing hits before the innermost hit on the track,
0122   // taking into account finite depth of each layer.
0123   unsigned int nMissHitAfter = 0;
0124   unsigned int nMissHitBefore = 0;
0125
0126   for (int i = 0; i < hp.numberOfAllHits(HitPattern::MISSING_INNER_HITS); i++) {
0127     uint32_t hit = hp.getHitPattern(HitPattern::MISSING_INNER_HITS, i);
0128     if (reco::HitPattern::trackerHitFilter(hit) && reco::HitPattern::missingHitFilter(hit)) {
0129       uint32_t subDet = reco::HitPattern::getSubStructure(hit);
0130       uint32_t layer = reco::HitPattern::getLayer(hit);
0131       DetInfo detInfo(subDet, layer);
0132       double minRZ = rangeRorZ_[detInfo].first;
0133
0134       //      cout << "subDet = " << subDet << " layer = " << layer << " minRZ = " << minRZ << endl;
0135
0136       if (this->barrel(subDet)) {
0137         if (vert.position().perp() < minRZ)
0138           nMissHitAfter++;
0139         else
0140           nMissHitBefore++;
0141       } else {
0142         if (fabs(vert.position().z()) < minRZ)
0143           nMissHitAfter++;
0144         else
0145           nMissHitBefore++;
0146       }
0147     }
0148   }
0149
0150   PFTrackHitInfo trackToVertex(nHitBefore, nMissHitBefore);
0151   PFTrackHitInfo trackFromVertex(nHitAfter, nMissHitAfter);
0152
0153   return pair<PFTrackHitInfo, PFTrackHitInfo>(trackToVertex, trackFromVertex);
0154 }
0155
0156 void PFCheckHitPattern::print(const TrackBaseRef track) const {
0157   // Get hit patterns of this track
0158   const reco::HitPattern& hp = track.get()->hitPattern();
0159
0160   cout << "=== Hits on Track ===" << endl;
0161   this->print(reco::HitPattern::TRACK_HITS, hp);
0162   cout << "=== Hits before track ===" << endl;
0163   this->print(reco::HitPattern::MISSING_INNER_HITS, hp);
0164 }
0165
0166 void PFCheckHitPattern::print(const reco::HitPattern::HitCategory category, const reco::HitPattern& hp) const {
0167   for (int i = 0; i < hp.numberOfAllHits(category); i++) {
0168     uint32_t hit = hp.getHitPattern(category, i);
0169     if (hp.trackerHitFilter(hit)) {
0170       uint32_t subdet = hp.getSubStructure(hit);
0171       uint32_t layer = hp.getLayer(hit);
0172       cout << "hit " << i << " subdet=" << subdet << " layer=" << layer << " type " << hp.getHitType(hit) << endl;
0173     }
0174   }
0175 }