File indexing completed on 2024-04-06 12:28:29
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
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
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();
0057 float coshEta = sqrt(1 + sqr(cotTheta));
0058
0059
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
0098 if (recHits.size() <= 2)
0099 return true;
0100
0101
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
0109 vector<GlobalPoint> globalPoss = getGlobalPoss(recHits);
0110
0111
0112 vector<GlobalVector> globalDirs = getGlobalDirs(globalPoss);
0113 if (globalDirs.empty())
0114 return false;
0115
0116 bool ok = true;
0117
0118
0119
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 }