File indexing completed on 2023-03-31 23:02:11
0001 #include <cmath>
0002
0003 #include "DataFormats/GeometryVector/interface/Basic2DVector.h"
0004 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0005 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0006 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0007 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0008 #include "FWCore/Framework/interface/ConsumesCollector.h"
0009 #include "FWCore/Framework/interface/Event.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "FWCore/Utilities/interface/ESInputTag.h"
0014 #include "RecoTracker/PixelLowPtUtilities/interface/HitInfo.h"
0015 #include "RecoTracker/PixelLowPtUtilities/interface/LowPtClusterShapeSeedComparitor.h"
0016 #include "RecoTracker/PixelTrackFitting/interface/CircleFromThreePoints.h"
0017 #include "RecoTracker/TkSeedingLayers/interface/SeedingHitSet.h"
0018
0019 namespace {
0020 typedef Basic2DVector<float> Vector2D;
0021
0022 inline float sqr(float x) { return x * x; }
0023
0024
0025 inline float areaParallelogram(const Vector2D& a, const Vector2D& b) { return a.x() * b.y() - a.y() * b.x(); }
0026
0027
0028
0029 inline bool getGlobalDirs(GlobalPoint const* g, GlobalVector* globalDirs) {
0030
0031 CircleFromThreePoints circle(g[0], g[1], g[2]);
0032
0033 float curvature = circle.curvature();
0034 if (0.f == curvature) {
0035 LogDebug("LowPtClusterShapeSeedComparitor")
0036 << "the curvature is null:"
0037 << "\n point1: " << g[0] << "\n point2: " << g[1] << "\n point3: " << g[2];
0038 return false;
0039 }
0040
0041
0042 Vector2D p[3];
0043 Vector2D c = circle.center();
0044 for (int i = 0; i != 3; i++)
0045 p[i] = g[i].basicVector().xy() - c;
0046
0047 float area = std::abs(areaParallelogram(p[1] - p[0], p[1]));
0048
0049 float a12 = std::asin(std::min(area * curvature * curvature, 1.f));
0050
0051 float slope = (g[1].z() - g[0].z()) / a12;
0052
0053
0054
0055 float cotTheta = slope * curvature;
0056 float sinTheta = 1.f / std::sqrt(1.f + sqr(cotTheta));
0057 float cosTheta = cotTheta * sinTheta;
0058
0059 if (areaParallelogram(p[0], p[1]) < 0)
0060 sinTheta = -sinTheta;
0061
0062 for (int i = 0; i != 3; i++) {
0063 Vector2D vl = p[i] * (curvature * sinTheta);
0064 globalDirs[i] = GlobalVector(-vl.y(), vl.x(), cosTheta);
0065 }
0066 return true;
0067 }
0068
0069
0070
0071 inline void getGlobalPos(const SeedingHitSet& hits, GlobalPoint* globalPoss) {
0072 for (unsigned int i = 0; i != hits.size(); ++i)
0073 globalPoss[i] = hits[i]->globalPosition();
0074 }
0075
0076 }
0077
0078
0079 LowPtClusterShapeSeedComparitor::LowPtClusterShapeSeedComparitor(const edm::ParameterSet& ps,
0080 edm::ConsumesCollector& iC)
0081 : thePixelClusterShapeCacheToken(
0082 iC.consumes<SiPixelClusterShapeCache>(ps.getParameter<edm::InputTag>("clusterShapeCacheSrc"))),
0083 theShapeFilterLabel_(ps.getParameter<std::string>("clusterShapeHitFilter")),
0084 clusterShapeHitFilterESToken_(iC.esConsumes(edm::ESInputTag("", theShapeFilterLabel_))),
0085 trackerTopologyESToken_(iC.esConsumes()) {}
0086
0087
0088 void LowPtClusterShapeSeedComparitor::init(const edm::Event& e, const edm::EventSetup& es) {
0089 clusterShapeHitFilter_ = &es.getData(clusterShapeHitFilterESToken_);
0090 trackerTopology_ = &es.getData(trackerTopologyESToken_);
0091 e.getByToken(thePixelClusterShapeCacheToken, thePixelClusterShapeCache);
0092 }
0093
0094 bool LowPtClusterShapeSeedComparitor::compatible(const SeedingHitSet& hits) const
0095
0096 {
0097 assert(hits.size() == 3);
0098
0099 if (clusterShapeHitFilter_ == nullptr)
0100 throw cms::Exception("LogicError") << "LowPtClusterShapeSeedComparitor: init(EventSetup) method was not called";
0101
0102
0103 GlobalPoint globalPoss[3];
0104 getGlobalPos(hits, globalPoss);
0105
0106
0107 GlobalVector globalDirs[3];
0108
0109 bool ok = getGlobalDirs(globalPoss, globalDirs);
0110
0111
0112
0113
0114 if (!ok) {
0115 LogDebug("LowPtClusterShapeSeedComparitor") << "curvarture 0:"
0116 << "\nnHits: " << hits.size() << " will say the seed is good anyway.";
0117 return true;
0118 }
0119
0120 for (int i = 0; i < 3; i++) {
0121 const SiPixelRecHit* pixelRecHit = dynamic_cast<const SiPixelRecHit*>(hits[i]->hit());
0122
0123 if (!pixelRecHit) {
0124 edm::LogError("LowPtClusterShapeSeedComparitor") << "this is not a pixel cluster";
0125 ok = false;
0126 break;
0127 }
0128
0129 if (!pixelRecHit->isValid()) {
0130 ok = false;
0131 break;
0132 }
0133
0134 LogDebug("LowPtClusterShapeSeedComparitor") << "about to compute compatibility."
0135 << "hit ptr: " << pixelRecHit << "global direction:" << globalDirs[i];
0136
0137 if (!clusterShapeHitFilter_->isCompatible(*pixelRecHit, globalDirs[i], *thePixelClusterShapeCache)) {
0138 LogTrace("LowPtClusterShapeSeedComparitor")
0139 << " clusShape is not compatible" << HitInfo::getInfo(*hits[i]->hit(), trackerTopology_);
0140
0141 ok = false;
0142 break;
0143 }
0144 }
0145
0146 return ok;
0147 }