special/plugins/HLTPixelClusterShapeFilter.cc

0001 #include "HLTrigger/HLTcore/interface/HLTFilter.h"
0002 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0003 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include "DataFormats/Common/interface/Handle.h"
0009 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0010 #include "DataFormats/HLTReco/interface/TriggerTypeDefs.h"
0011 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0012 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0013 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0014 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0015 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0016 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0017 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0018 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0019 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0020
0021 //
0022 // class declaration
0023 //
0024
0025 class HLTPixelClusterShapeFilter : public HLTFilter {
0026 public:
0027   explicit HLTPixelClusterShapeFilter(const edm::ParameterSet &);
0028   ~HLTPixelClusterShapeFilter() override;
0029   static void fillDescriptions(edm::ConfigurationDescriptions &descriptions);
0030
0031 private:
0032   edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> const trackerDigiGeometryRecordToken_;
0033   edm::EDGetTokenT<SiPixelRecHitCollection> inputToken_;
0034   edm::InputTag inputTag_;  // input tag identifying product containing pixel clusters
0035   double minZ_;             // beginning z-vertex position
0036   double maxZ_;             // end z-vertex position
0037   double zStep_;            // size of steps in z-vertex test
0038
0039   std::vector<double> clusterPars_;  //pixel cluster polynomial pars for vertex compatibility cut
0040   int nhitsTrunc_;                   //maximum pixel clusters to apply compatibility check
0041   double clusterTrunc_;              //maximum vertex compatibility value for event rejection
0042
0043   struct VertexHit {
0044     float z;
0045     float r;
0046     float w;
0047   };
0048
0049   bool hltFilter(edm::Event &,
0050                  const edm::EventSetup &,
0051                  trigger::TriggerFilterObjectWithRefs &filterproduct) const override;
0052   int getContainedHits(const std::vector<VertexHit> &hits, double z0, double &chi) const;
0053 };
0054
0055 //
0056 // constructors and destructor
0057 //
0058
0059 HLTPixelClusterShapeFilter::HLTPixelClusterShapeFilter(const edm::ParameterSet &config)
0060     : HLTFilter(config),
0061       trackerDigiGeometryRecordToken_(esConsumes()),
0062       inputTag_(config.getParameter<edm::InputTag>("inputTag")),
0063       minZ_(config.getParameter<double>("minZ")),
0064       maxZ_(config.getParameter<double>("maxZ")),
0065       zStep_(config.getParameter<double>("zStep")),
0066       clusterPars_(config.getParameter<std::vector<double> >("clusterPars")),
0067       nhitsTrunc_(config.getParameter<int>("nhitsTrunc")),
0068       clusterTrunc_(config.getParameter<double>("clusterTrunc")) {
0069   inputToken_ = consumes<SiPixelRecHitCollection>(inputTag_);
0070   LogDebug("") << "Using the " << inputTag_ << " input collection";
0071 }
0072
0073 HLTPixelClusterShapeFilter::~HLTPixelClusterShapeFilter() = default;
0074
0075 void HLTPixelClusterShapeFilter::fillDescriptions(edm::ConfigurationDescriptions &descriptions) {
0076   edm::ParameterSetDescription desc;
0077   makeHLTFilterDescription(desc);
0078   desc.add<edm::InputTag>("inputTag", edm::InputTag("hltSiPixelRecHits"));
0079   desc.add<double>("minZ", -20.0);
0080   desc.add<double>("maxZ", 20.05);
0081   desc.add<double>("zStep", 0.2);
0082   std::vector<double> temp;
0083   temp.push_back(0.0);
0084   temp.push_back(0.0045);
0085   desc.add<std::vector<double> >("clusterPars", temp);
0086   desc.add<int>("nhitsTrunc", 150.);
0087   desc.add<double>("clusterTrunc", 2.0);
0088   descriptions.add("hltPixelClusterShapeFilter", desc);
0089 }
0090
0091 //
0092 // member functions
0093 //
0094
0095 // ------------ method called to produce the data  ------------
0096 bool HLTPixelClusterShapeFilter::hltFilter(edm::Event &event,
0097                                            const edm::EventSetup &iSetup,
0098                                            trigger::TriggerFilterObjectWithRefs &filterproduct) const {
0099   // All HLT filters must create and fill an HLT filter object,
0100   // recording any reconstructed physics objects satisfying (or not)
0101   // this HLT filter, and place it in the Event.
0102
0103   // The filter object
0104   if (saveTags())
0105     filterproduct.addCollectionTag(inputTag_);
0106   bool accept = true;
0107
0108   // get hold of products from Event
0109   edm::Handle<SiPixelRecHitCollection> hRecHits;
0110   event.getByToken(inputToken_, hRecHits);
0111
0112   // get tracker geometry
0113   if (hRecHits.isValid()) {
0114     auto const &trackerHandle = iSetup.getHandle(trackerDigiGeometryRecordToken_);
0115
0116     const TrackerGeometry *tgeo = trackerHandle.product();
0117     const SiPixelRecHitCollection *hits = hRecHits.product();
0118
0119     // loop over pixel rechits
0120     int nPxlHits = 0;
0121     std::vector<VertexHit> vhits;
0122     for (auto const &hit : hits->data()) {
0123       if (!hit.isValid())
0124         continue;
0125       ++nPxlHits;
0126       DetId id(hit.geographicalId());
0127       if (id.subdetId() != int(PixelSubdetector::PixelBarrel))
0128         continue;
0129       const PixelGeomDetUnit *pgdu = static_cast<const PixelGeomDetUnit *>(tgeo->idToDet(id));
0130       if (true) {
0131         const PixelTopology *pixTopo = &(pgdu->specificTopology());
0132         std::vector<SiPixelCluster::Pixel> pixels(hit.cluster()->pixels());
0133         bool pixelOnEdge = false;
0134         for (std::vector<SiPixelCluster::Pixel>::const_iterator pixel = pixels.begin(); pixel != pixels.end();
0135              ++pixel) {
0136           int pixelX = pixel->x;
0137           int pixelY = pixel->y;
0138           if (pixTopo->isItEdgePixelInX(pixelX) || pixTopo->isItEdgePixelInY(pixelY)) {
0139             pixelOnEdge = true;
0140             break;
0141           }
0142         }
0143         if (pixelOnEdge)
0144           continue;
0145       }
0146
0147       LocalPoint lpos = LocalPoint(hit.localPosition().x(), hit.localPosition().y(), hit.localPosition().z());
0148       GlobalPoint gpos = pgdu->toGlobal(lpos);
0149       VertexHit vh;
0150       vh.z = gpos.z();
0151       vh.r = gpos.perp();
0152       vh.w = hit.cluster()->sizeY();
0153       vhits.push_back(vh);
0154     }
0155
0156     // estimate z-position from cluster lengths
0157     double zest = 0.0;
0158     int nhits = 0, nhits_max = 0;
0159     double chi = 0, chi_max = 1e+9;
0160     for (double z0 = minZ_; z0 <= maxZ_; z0 += zStep_) {
0161       nhits = getContainedHits(vhits, z0, chi);
0162       if (nhits == 0)
0163         continue;
0164       if (nhits > nhits_max) {
0165         chi_max = 1e+9;
0166         nhits_max = nhits;
0167       }
0168       if (nhits >= nhits_max && chi < chi_max) {
0169         chi_max = chi;
0170         zest = z0;
0171       }
0172     }
0173
0174     chi = 0;
0175     int nbest = 0, nminus = 0, nplus = 0;
0176     nbest = getContainedHits(vhits, zest, chi);
0177     nminus = getContainedHits(vhits, zest - 10., chi);
0178     nplus = getContainedHits(vhits, zest + 10., chi);
0179
0180     double clusVtxQual = 0.0;
0181     if ((nminus + nplus) > 0)
0182       clusVtxQual = (2.0 * nbest) / (nminus + nplus);  // A/B
0183     else if (nbest > 0)
0184       clusVtxQual = 1000.0;  // A/0 (set to arbitrarily large number)
0185     else
0186       clusVtxQual = 0;  // 0/0 (already the default)
0187
0188     // construct polynomial cut on cluster vertex quality vs. npixelhits
0189     double polyCut = 0;
0190     for (unsigned int i = 0; i < clusterPars_.size(); i++) {
0191       polyCut += clusterPars_[i] * std::pow((double)nPxlHits, (int)i);
0192     }
0193     if (nPxlHits < nhitsTrunc_)
0194       polyCut = 0;  // don't use cut below nhitsTrunc_ pixel hits
0195     if (polyCut > clusterTrunc_ && clusterTrunc_ > 0)
0196       polyCut = clusterTrunc_;  // no cut above clusterTrunc_
0197
0198     if (clusVtxQual < polyCut)
0199       accept = false;
0200   }
0201
0202   // return with final filter decision
0203   return accept;
0204 }
0205
0206 int HLTPixelClusterShapeFilter::getContainedHits(const std::vector<VertexHit> &hits, double z0, double &chi) const {
0207   // Calculate number of hits contained in v-shaped window in cluster y-width vs. z-position.
0208   int n = 0;
0209   chi = 0.;
0210
0211   for (auto hit : hits) {
0212     double p = 2 * fabs(hit.z - z0) / hit.r + 0.5;  // FIXME
0213     if (fabs(p - hit.w) <= 1.) {
0214       chi += fabs(p - hit.w);
0215       n++;
0216     }
0217   }
0218   return n;
0219 }
0220
0221 // define as a framework module
0222 #include "FWCore/Framework/interface/MakerMacros.h"
0223 DEFINE_FWK_MODULE(HLTPixelClusterShapeFilter);