PixelLowPtUtilities/plugins/ClusterShapeSeedComparitor.cc

0001 #include "RecoTracker/TkSeedingLayers/interface/SeedComparitor.h"
0002 #include "RecoTracker/TkSeedingLayers/interface/SeedComparitorFactory.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "RecoTracker/PixelLowPtUtilities/interface/ClusterShapeHitFilter.h"
0006 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0007 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0008 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
0009 #include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
0010 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0011 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
0012 #include "RecoTracker/TkSeedingLayers/interface/SeedingHitSet.h"
0013 #include "MagneticField/Engine/interface/MagneticField.h"
0014 #include "FWCore/Framework/interface/Event.h"
0015 #include "FWCore/Framework/interface/ConsumesCollector.h"
0016 #include "FWCore/Utilities/interface/EDGetToken.h"
0017 #include "FWCore/Utilities/interface/ESGetToken.h"
0018 #include "FWCore/Utilities/interface/ESInputTag.h"
0019 #include "DataFormats/Common/interface/Handle.h"
0020 #include "DataFormats/SiPixelCluster/interface/SiPixelClusterShapeCache.h"
0021 #include <cstdio>
0022 #include <cassert>
0023
0024 #include <iostream>
0025
0026 class PixelClusterShapeSeedComparitor : public SeedComparitor {
0027 public:
0028   PixelClusterShapeSeedComparitor(const edm::ParameterSet &cfg, edm::ConsumesCollector &iC);
0029   ~PixelClusterShapeSeedComparitor() override;
0030   void init(const edm::Event &ev, const edm::EventSetup &es) override;
0031   bool compatible(const SeedingHitSet &hits) const override { return true; }
0032   bool compatible(const TrajectoryStateOnSurface &, SeedingHitSet::ConstRecHitPointer hit) const override;
0033   bool compatible(const SeedingHitSet &hits,
0034                   const GlobalTrajectoryParameters &helixStateAtVertex,
0035                   const FastHelix &helix) const override;
0036
0037 private:
0038   bool compatibleHit(const TrackingRecHit &hit, const GlobalVector &direction) const;
0039
0040   std::string filterName_;
0041   const edm::ESGetToken<ClusterShapeHitFilter, CkfComponentsRecord> clusterShapeHitFilterESToken_;
0042   const ClusterShapeHitFilter *clusterShapeHitFilter_;
0043   edm::EDGetTokenT<SiPixelClusterShapeCache> pixelClusterShapeCacheToken_;
0044   const SiPixelClusterShapeCache *pixelClusterShapeCache_;
0045   const bool filterAtHelixStage_;
0046   const bool filterPixelHits_, filterStripHits_;
0047 };
0048
0049 PixelClusterShapeSeedComparitor::PixelClusterShapeSeedComparitor(const edm::ParameterSet &cfg,
0050                                                                  edm::ConsumesCollector &iC)
0051     : filterName_(cfg.getParameter<std::string>("ClusterShapeHitFilterName")),
0052       clusterShapeHitFilterESToken_(iC.esConsumes(edm::ESInputTag("", filterName_))),
0053       pixelClusterShapeCache_(nullptr),
0054       filterAtHelixStage_(cfg.getParameter<bool>("FilterAtHelixStage")),
0055       filterPixelHits_(cfg.getParameter<bool>("FilterPixelHits")),
0056       filterStripHits_(cfg.getParameter<bool>("FilterStripHits")) {
0057   if (filterPixelHits_) {
0058     pixelClusterShapeCacheToken_ =
0059         iC.consumes<SiPixelClusterShapeCache>(cfg.getParameter<edm::InputTag>("ClusterShapeCacheSrc"));
0060   }
0061 }
0062
0063 PixelClusterShapeSeedComparitor::~PixelClusterShapeSeedComparitor() {}
0064
0065 void PixelClusterShapeSeedComparitor::init(const edm::Event &ev, const edm::EventSetup &es) {
0066   clusterShapeHitFilter_ = &es.getData(clusterShapeHitFilterESToken_);
0067   if (filterPixelHits_) {
0068     edm::Handle<SiPixelClusterShapeCache> hcache;
0069     ev.getByToken(pixelClusterShapeCacheToken_, hcache);
0070     pixelClusterShapeCache_ = hcache.product();
0071   }
0072 }
0073
0074 bool PixelClusterShapeSeedComparitor::compatible(const TrajectoryStateOnSurface &tsos,
0075                                                  SeedingHitSet::ConstRecHitPointer hit) const {
0076   if (filterAtHelixStage_)
0077     return true;
0078   assert(hit->isValid() && tsos.isValid());
0079   return compatibleHit(*hit, tsos.globalDirection());
0080 }
0081
0082 bool PixelClusterShapeSeedComparitor::compatible(const SeedingHitSet &hits,
0083                                                  const GlobalTrajectoryParameters &helixStateAtVertex,
0084                                                  const FastHelix &helix) const {
0085   if (!filterAtHelixStage_)
0086     return true;
0087
0088   if (!helix.isValid()              //check still if it's a straight line, which are OK
0089       && !helix.circle().isLine())  //complain if it's not even a straight line
0090     edm::LogWarning("InvalidHelix") << "PixelClusterShapeSeedComparitor helix is not valid, result is bad";
0091
0092   float xc = helix.circle().x0(), yc = helix.circle().y0();
0093
0094   GlobalPoint vertex = helixStateAtVertex.position();
0095   GlobalVector momvtx = helixStateAtVertex.momentum();
0096   float x0 = vertex.x(), y0 = vertex.y();
0097   for (unsigned int i = 0, n = hits.size(); i < n; ++i) {
0098     auto const &hit = *hits[i];
0099     GlobalPoint pos = hit.globalPosition();
0100     float x1 = pos.x(), y1 = pos.y(), dx1 = x1 - xc, dy1 = y1 - yc;
0101
0102     // now figure out the proper tangent vector
0103     float perpx = -dy1, perpy = dx1;
0104     if (perpx * (x1 - x0) + perpy * (y1 - y0) < 0) {
0105       perpy = -perpy;
0106       perpx = -perpx;
0107     }
0108
0109     // now normalize (perpx, perpy, 1.0) to momentum (px, py, pz)
0110     float perp2 = perpx * perpx + perpy * perpy;
0111     float pmom2 = momvtx.x() * momvtx.x() + momvtx.y() * momvtx.y(), momz2 = momvtx.z() * momvtx.z(),
0112           mom2 = pmom2 + momz2;
0113     float perpscale = sqrt(pmom2 / mom2 / perp2), zscale = sqrt((1 - pmom2 / mom2));
0114     GlobalVector gdir(perpx * perpscale, perpy * perpscale, (momvtx.z() > 0 ? zscale : -zscale));
0115
0116     if (!compatibleHit(hit, gdir)) {
0117       return false;  // not yet
0118     }
0119   }
0120   return true;
0121 }
0122
0123 bool PixelClusterShapeSeedComparitor::compatibleHit(const TrackingRecHit &hit, const GlobalVector &direction) const {
0124   if (hit.geographicalId().subdetId() <= 2) {
0125     if (!filterPixelHits_)
0126       return true;
0127     const SiPixelRecHit *pixhit = dynamic_cast<const SiPixelRecHit *>(&hit);
0128     if (pixhit == nullptr)
0129       throw cms::Exception("LogicError", "Found a valid hit on the pixel detector which is not a SiPixelRecHit\n");
0130     //printf("Cheching hi hit on detid %10d, local direction is x = %9.6f, y = %9.6f, z = %9.6f\n", hit.geographicalId().rawId(), direction.x(), direction.y(), direction.z());
0131     return clusterShapeHitFilter_->isCompatible(*pixhit, direction, *pixelClusterShapeCache_);
0132   } else {
0133     if (!filterStripHits_)
0134       return true;
0135     const std::type_info &tid = typeid(*&hit);
0136     if (tid == typeid(SiStripMatchedRecHit2D)) {
0137       const SiStripMatchedRecHit2D *matchedHit = dynamic_cast<const SiStripMatchedRecHit2D *>(&hit);
0138       assert(matchedHit != nullptr);
0139       return (
0140           clusterShapeHitFilter_->isCompatible(DetId(matchedHit->monoId()), matchedHit->monoCluster(), direction) &&
0141           clusterShapeHitFilter_->isCompatible(DetId(matchedHit->stereoId()), matchedHit->stereoCluster(), direction));
0142     } else if (tid == typeid(SiStripRecHit2D)) {
0143       const SiStripRecHit2D *recHit = dynamic_cast<const SiStripRecHit2D *>(&hit);
0144       assert(recHit != nullptr);
0145       return clusterShapeHitFilter_->isCompatible(*recHit, direction);
0146     } else if (tid == typeid(ProjectedSiStripRecHit2D)) {
0147       const ProjectedSiStripRecHit2D *precHit = dynamic_cast<const ProjectedSiStripRecHit2D *>(&hit);
0148       assert(precHit != nullptr);
0149       return clusterShapeHitFilter_->isCompatible(precHit->originalHit(), direction);
0150     } else {
0151       //printf("Questo e' un %s, che ci fo?\n", tid.name());
0152       return true;
0153     }
0154   }
0155 }
0156
0157 DEFINE_EDM_PLUGIN(SeedComparitorFactory, PixelClusterShapeSeedComparitor, "PixelClusterShapeSeedComparitor");