PixelLowPtUtilities/plugins/PixelVertexProducerClusters.cc

0001 #include "PixelVertexProducerClusters.h"
0002
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007
0008 #include "DataFormats/VertexReco/interface/Vertex.h"
0009 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0010
0011 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0012 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0013 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0014
0015 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0016
0017 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0018 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0019
0020 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0021
0022 #include <vector>
0023 #include <algorithm>
0024
0025 using namespace std;
0026
0027 namespace {
0028   class VertexHit {
0029   public:
0030     float z;
0031     float r;
0032     int w;
0033   };
0034
0035   /*****************************************************************************/
0036   int getContainedHits(const vector<VertexHit>& hits, float z0, float& chi) {
0037     int n = 0;
0038     chi = 0.;
0039
0040     for (vector<VertexHit>::const_iterator hit = hits.begin(); hit != hits.end(); hit++) {
0041       // Predicted cluster width in y direction
0042       float p = 2 * fabs(hit->z - z0) / hit->r + 0.5;  // FIXME
0043
0044       if (fabs(p - hit->w) <= 1.) {
0045         chi += fabs(p - hit->w);
0046         n++;
0047       }
0048     }
0049
0050     return n;
0051   }
0052 }  // namespace
0053
0054 /*****************************************************************************/
0055 PixelVertexProducerClusters::PixelVertexProducerClusters(const edm::ParameterSet& ps)
0056     : geomToken_(esConsumes()), pixelToken_(consumes<SiPixelRecHitCollection>(edm::InputTag("siPixelRecHits"))) {
0057   // Product
0058   produces<reco::VertexCollection>();
0059 }
0060
0061 /*****************************************************************************/
0062 PixelVertexProducerClusters::~PixelVertexProducerClusters() {}
0063
0064 /*****************************************************************************/
0065 void PixelVertexProducerClusters::produce(edm::StreamID, edm::Event& ev, const edm::EventSetup& es) const {
0066   // Get tracker geometry
0067   const TrackerGeometry* theTracker = &es.getData(geomToken_);
0068
0069   // Get pixel hit collections
0070   edm::Handle<SiPixelRecHitCollection> pixelColl;
0071   ev.getByToken(pixelToken_, pixelColl);
0072
0073   const SiPixelRecHitCollection* thePixelHits = pixelColl.product();
0074
0075   auto vertices = std::make_unique<reco::VertexCollection>();
0076
0077   if (!thePixelHits->empty()) {
0078     vector<VertexHit> hits;
0079
0080     for (SiPixelRecHitCollection::DataContainer::const_iterator recHit = thePixelHits->data().begin(),
0081                                                                 recHitEnd = thePixelHits->data().end();
0082          recHit != recHitEnd;
0083          ++recHit) {
0084       if (recHit->isValid()) {
0085         //      if(!(recHit->isOnEdge() || recHit->hasBadPixels()))
0086         DetId id = recHit->geographicalId();
0087         const PixelGeomDetUnit* pgdu = dynamic_cast<const PixelGeomDetUnit*>(theTracker->idToDetUnit(id));
0088         const PixelTopology* theTopol = (&(pgdu->specificTopology()));
0089         vector<SiPixelCluster::Pixel> pixels = recHit->cluster()->pixels();
0090
0091         bool pixelOnEdge = false;
0092         for (vector<SiPixelCluster::Pixel>::const_iterator pixel = pixels.begin(); pixel != pixels.end(); pixel++) {
0093           int pos_x = (int)pixel->x;
0094           int pos_y = (int)pixel->y;
0095
0096           if (theTopol->isItEdgePixelInX(pos_x) || theTopol->isItEdgePixelInY(pos_y)) {
0097             pixelOnEdge = true;
0098             break;
0099           }
0100         }
0101
0102         if (!pixelOnEdge)
0103           if (id.subdetId() == int(PixelSubdetector::PixelBarrel)) {
0104             LocalPoint lpos =
0105                 LocalPoint(recHit->localPosition().x(), recHit->localPosition().y(), recHit->localPosition().z());
0106
0107             GlobalPoint gpos = theTracker->idToDet(id)->toGlobal(lpos);
0108
0109             VertexHit hit;
0110             hit.z = gpos.z();
0111             hit.r = gpos.perp();
0112             hit.w = recHit->cluster()->sizeY();
0113
0114             hits.push_back(hit);
0115           }
0116       }
0117     }
0118
0119     int nhits;
0120     int nhits_max = 0;
0121     float chi;
0122     float chi_max = 1e+9;
0123
0124     float zest = 0;
0125
0126     for (float z0 = -15.9; z0 <= 15.95; z0 += 0.1) {
0127       nhits = getContainedHits(hits, z0, chi);
0128
0129       if (nhits > 0) {
0130         if (nhits > nhits_max) {
0131           chi_max = 1e+9;
0132           nhits_max = nhits;
0133         }
0134
0135         if (nhits >= nhits_max)
0136           if (chi < chi_max) {
0137             chi_max = chi;
0138             zest = z0;
0139           }
0140       }
0141     }
0142
0143     LogTrace("MinBiasTracking") << "  [vertex position] estimated = " << zest
0144                                 << " | pixel barrel hits = " << thePixelHits->size();
0145
0146     reco::Vertex::Error err;
0147     err(2, 2) = 0.6 * 0.6;
0148     reco::Vertex ver(reco::Vertex::Point(0, 0, zest), err, 0, 1, 1);
0149     vertices->push_back(ver);
0150   }
0151
0152   ev.put(std::move(vertices));
0153 }