PixelLowPtUtilities/src/ClusterShapeTrackFilter.cc

0001 #include "DataFormats/SiPixelCluster/interface/SiPixelClusterShapeCache.h"
0002 #include "DataFormats/TrackReco/interface/Track.h"
0003 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0004 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0005 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0008 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0009 #include "RecoTracker/PixelLowPtUtilities/interface/ClusterShapeHitFilter.h"
0010 #include "RecoTracker/PixelLowPtUtilities/interface/ClusterShapeTrackFilter.h"
0011 #include "RecoTracker/PixelLowPtUtilities/interface/HitInfo.h"
0012 #include "RecoTracker/PixelTrackFitting/interface/CircleFromThreePoints.h"
0013
0014 inline float sqr(float x) { return x * x; }
0015
0016 using namespace std;
0017
0018 /*****************************************************************************/
0019 ClusterShapeTrackFilter::ClusterShapeTrackFilter(const SiPixelClusterShapeCache* cache,
0020                                                  double ptmin,
0021                                                  double ptmax,
0022                                                  const TrackerGeometry* tracker,
0023                                                  const ClusterShapeHitFilter* shape,
0024                                                  const TrackerTopology* ttopo)
0025     : theTracker(tracker), theFilter(shape), theClusterShapeCache(cache), tTopo(ttopo), ptMin(ptmin), ptMax(ptmax) {}
0026
0027 /*****************************************************************************/
0028 ClusterShapeTrackFilter::~ClusterShapeTrackFilter() {}
0029
0030 /*****************************************************************************/
0031 float ClusterShapeTrackFilter::areaParallelogram(const Global2DVector& a, const Global2DVector& b) const {
0032   return a.x() * b.y() - a.y() * b.x();
0033 }
0034
0035 /*****************************************************************************/
0036 vector<GlobalVector> ClusterShapeTrackFilter::getGlobalDirs(const vector<GlobalPoint>& g) const {
0037   // Get 2d points
0038   vector<Global2DVector> p;
0039   for (vector<GlobalPoint>::const_iterator ig = g.begin(); ig != g.end(); ig++)
0040     p.push_back(Global2DVector(ig->x(), ig->y()));
0041
0042   //
0043   vector<GlobalVector> globalDirs;
0044
0045   // Determine circle
0046   CircleFromThreePoints circle(g[0], g[1], g[2]);
0047
0048   if (circle.curvature() != 0.) {
0049     Global2DVector c(circle.center().x(), circle.center().y());
0050
0051     float rad2 = (p[0] - c).mag2();
0052     float a12 = asin(fabsf(areaParallelogram(p[0] - c, p[1] - c)) / rad2);
0053
0054     float slope = (g[1].z() - g[0].z()) / a12;
0055
0056     float cotTheta = slope * circle.curvature();  // == sinhEta
0057     float coshEta = sqrt(1 + sqr(cotTheta));      // == 1/sinTheta
0058
0059     // Calculate globalDirs
0060     float sinTheta = 1. / coshEta;
0061     float cosTheta = cotTheta * sinTheta;
0062
0063     int dir;
0064     if (areaParallelogram(p[0] - c, p[1] - c) > 0)
0065       dir = 1;
0066     else
0067       dir = -1;
0068
0069     float curvature = circle.curvature();
0070
0071     for (vector<Global2DVector>::const_iterator ip = p.begin(); ip != p.end(); ip++) {
0072       Global2DVector v = (*ip - c) * curvature * dir;
0073       globalDirs.push_back(GlobalVector(-v.y() * sinTheta, v.x() * sinTheta, cosTheta));
0074     }
0075   }
0076
0077   return globalDirs;
0078 }
0079
0080 /*****************************************************************************/
0081 vector<GlobalPoint> ClusterShapeTrackFilter::getGlobalPoss(const vector<const TrackingRecHit*>& recHits) const {
0082   vector<GlobalPoint> globalPoss;
0083
0084   for (vector<const TrackingRecHit*>::const_iterator recHit = recHits.begin(); recHit != recHits.end(); recHit++) {
0085     DetId detId = (*recHit)->geographicalId();
0086
0087     GlobalPoint gpos = theTracker->idToDet(detId)->toGlobal((*recHit)->localPosition());
0088
0089     globalPoss.push_back(gpos);
0090   }
0091
0092   return globalPoss;
0093 }
0094
0095 /*****************************************************************************/
0096 bool ClusterShapeTrackFilter::operator()(const reco::Track* track, const vector<const TrackingRecHit*>& recHits) const {
0097   // Do not even look at pairs
0098   if (recHits.size() <= 2)
0099     return true;
0100
0101   // Check pt
0102   if (track->pt() < ptMin || track->pt() > ptMax) {
0103     LogTrace("ClusterShapeTrackFilter") << "  [ClusterShapeTrackFilter] pt not in range: " << ptMin << " "
0104                                         << track->pt() << " " << ptMax;
0105     return false;
0106   }
0107
0108   // Get global positions
0109   vector<GlobalPoint> globalPoss = getGlobalPoss(recHits);
0110
0111   // Get global directions
0112   vector<GlobalVector> globalDirs = getGlobalDirs(globalPoss);
0113   if (globalDirs.empty())
0114     return false;
0115
0116   bool ok = true;
0117
0118   // Check whether shape of pixel cluster is compatible
0119   // with local track direction
0120   for (unsigned int i = 0; i < recHits.size(); i++) {
0121     const SiPixelRecHit* pixelRecHit = dynamic_cast<const SiPixelRecHit*>(recHits[i]);
0122
0123     if (!pixelRecHit->isValid()) {
0124       ok = false;
0125       break;
0126     }
0127
0128     if (!theFilter->isCompatible(*pixelRecHit, globalDirs[i], *theClusterShapeCache)) {
0129       LogTrace("ClusterShapeTrackFilter")
0130           << "  [ClusterShapeTrackFilter] clusShape problem" << HitInfo::getInfo(*recHits[i], tTopo);
0131
0132       ok = false;
0133       break;
0134     }
0135   }
0136
0137   return ok;
0138 }