File indexing completed on 2024-04-06 12:28:28
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
0042 float p = 2 * fabs(hit->z - z0) / hit->r + 0.5;
0043
0044 if (fabs(p - hit->w) <= 1.) {
0045 chi += fabs(p - hit->w);
0046 n++;
0047 }
0048 }
0049
0050 return n;
0051 }
0052 }
0053
0054
0055 PixelVertexProducerClusters::PixelVertexProducerClusters(const edm::ParameterSet& ps)
0056 : geomToken_(esConsumes()), pixelToken_(consumes<SiPixelRecHitCollection>(edm::InputTag("siPixelRecHits"))) {
0057
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
0067 const TrackerGeometry* theTracker = &es.getData(geomToken_);
0068
0069
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
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 }