PixelLowPtUtilities/src/LowPtClusterShapeSeedComparitor.cc

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     // Determine circle
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     // Get 2d points
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     // Calculate globalDirs
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 }  // namespace
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 //(const reco::Track* track, const vector<const TrackingRecHit *> & recHits) const
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   // Get global positions
0103   GlobalPoint globalPoss[3];
0104   getGlobalPos(hits, globalPoss);
0105
0106   // Get global directions
0107   GlobalVector globalDirs[3];
0108
0109   bool ok = getGlobalDirs(globalPoss, globalDirs);
0110
0111   // Check whether shape of pixel cluster is compatible
0112   // with local track direction
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 }