HiTracking/plugins/HIPixelClusterVtxProducer.cc

0001 #include "HIPixelClusterVtxProducer.h"
0002
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 #include "FWCore/Framework/interface/Frameworkfwd.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006
0007 #include "DataFormats/VertexReco/interface/Vertex.h"
0008 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0009
0010 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0011 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0012 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0013
0014 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0015 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0016
0017 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0018
0019 #include <fstream>
0020 #include <iostream>
0021 #include <vector>
0022 #include <algorithm>
0023
0024 /*****************************************************************************/
0025 HIPixelClusterVtxProducer::HIPixelClusterVtxProducer(const edm::ParameterSet &ps)
0026     : srcPixelsString_(ps.getParameter<std::string>("pixelRecHits")),
0027       trackerToken_(esConsumes()),
0028       minZ_(ps.getParameter<double>("minZ")),
0029       maxZ_(ps.getParameter<double>("maxZ")),
0030       zStep_(ps.getParameter<double>("zStep"))
0031
0032 {
0033   // Constructor
0034   produces<reco::VertexCollection>();
0035   srcPixels_ = (consumes<SiPixelRecHitCollection>(srcPixelsString_));
0036 }
0037
0038 /*****************************************************************************/
0039 HIPixelClusterVtxProducer::~HIPixelClusterVtxProducer() {
0040   // Destructor
0041 }
0042
0043 /*****************************************************************************/
0044 void HIPixelClusterVtxProducer::produce(edm::Event &ev, const edm::EventSetup &es) {
0045   // new vertex collection
0046   auto vertices = std::make_unique<reco::VertexCollection>();
0047
0048   // get pixel rechits
0049   edm::Handle<SiPixelRecHitCollection> hRecHits;
0050   ev.getByToken(srcPixels_, hRecHits);
0051
0052   // get tracker geometry
0053   if (hRecHits.isValid()) {
0054     const TrackerGeometry *tgeo = &es.getData(trackerToken_);
0055     const SiPixelRecHitCollection *hits = hRecHits.product();
0056
0057     // loop over pixel rechits
0058     std::vector<VertexHit> vhits;
0059     for (SiPixelRecHitCollection::DataContainer::const_iterator hit = hits->data().begin(), end = hits->data().end();
0060          hit != end;
0061          ++hit) {
0062       if (!hit->isValid())
0063         continue;
0064       DetId id(hit->geographicalId());
0065       if (id.subdetId() != int(PixelSubdetector::PixelBarrel))
0066         continue;
0067       const PixelGeomDetUnit *pgdu = static_cast<const PixelGeomDetUnit *>(tgeo->idToDet(id));
0068       if (true) {
0069         const PixelTopology *pixTopo = &(pgdu->specificTopology());
0070         std::vector<SiPixelCluster::Pixel> pixels(hit->cluster()->pixels());
0071         bool pixelOnEdge = false;
0072         for (std::vector<SiPixelCluster::Pixel>::const_iterator pixel = pixels.begin(); pixel != pixels.end();
0073              ++pixel) {
0074           int pixelX = pixel->x;
0075           int pixelY = pixel->y;
0076           if (pixTopo->isItEdgePixelInX(pixelX) || pixTopo->isItEdgePixelInY(pixelY)) {
0077             pixelOnEdge = true;
0078             break;
0079           }
0080         }
0081         if (pixelOnEdge)
0082           continue;
0083       }
0084
0085       LocalPoint lpos = LocalPoint(hit->localPosition().x(), hit->localPosition().y(), hit->localPosition().z());
0086       GlobalPoint gpos = pgdu->toGlobal(lpos);
0087       VertexHit vh;
0088       vh.z = gpos.z();
0089       vh.r = gpos.perp();
0090       vh.w = hit->cluster()->sizeY();
0091       vhits.push_back(vh);
0092     }
0093
0094     // estimate z-position from cluster lengths
0095     double zest = 0.0;
0096     int nhits = 0, nhits_max = 0;
0097     double chi = 0, chi_max = 1e+9;
0098     for (double z0 = minZ_; z0 <= maxZ_; z0 += zStep_) {
0099       nhits = getContainedHits(vhits, z0, chi);
0100       if (nhits == 0)
0101         continue;
0102       if (nhits > nhits_max) {
0103         chi_max = 1e+9;
0104         nhits_max = nhits;
0105       }
0106       if (nhits >= nhits_max && chi < chi_max) {
0107         chi_max = chi;
0108         zest = z0;
0109       }
0110     }
0111
0112     LogTrace("MinBiasTracking") << "  [vertex position] estimated = " << zest
0113                                 << " | pixel barrel hits = " << vhits.size();
0114
0115     // put 1-d vertex and dummy errors into collection
0116     reco::Vertex::Error err;
0117     err(2, 2) = 0.6 * 0.6;
0118     reco::Vertex ver(reco::Vertex::Point(0, 0, zest), err, 0, 1, 1);
0119     vertices->push_back(ver);
0120   }
0121
0122   ev.put(std::move(vertices));
0123 }
0124
0125 /*****************************************************************************/
0126 int HIPixelClusterVtxProducer::getContainedHits(const std::vector<VertexHit> &hits, double z0, double &chi) {
0127   // Calculate number of hits contained in v-shaped window in cluster y-width vs. z-position.
0128   int n = 0;
0129   chi = 0.;
0130
0131   for (std::vector<VertexHit>::const_iterator hit = hits.begin(); hit != hits.end(); hit++) {
0132     double p = 2 * fabs(hit->z - z0) / hit->r + 0.5;  // FIXME
0133     if (TMath::Abs(p - hit->w) <= 1.) {
0134       chi += fabs(p - hit->w);
0135       n++;
0136     }
0137   }
0138   return n;
0139 }