File indexing completed on 2023-03-17 11:18:16
0001 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/ServiceRegistry/interface/Service.h"
0005 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0006 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0007 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0008
0009 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0010 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0011 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0012 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0013 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0014
0015 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0016 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0017
0018 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
0019
0020 #include "DataFormats/Provenance/interface/RunLumiEventNumber.h"
0021
0022 #include <fstream>
0023 #include <iostream>
0024 #include <vector>
0025 #include <algorithm>
0026
0027
0028 #include <TH1.h>
0029
0030 class HIPixelClusterVtxAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0031 public:
0032 explicit HIPixelClusterVtxAnalyzer(const edm::ParameterSet &ps);
0033
0034 private:
0035 struct VertexHit {
0036 float z;
0037 float r;
0038 float w;
0039 };
0040
0041 virtual void analyze(const edm::Event &ev, const edm::EventSetup &es);
0042 int getContainedHits(const std::vector<VertexHit> &hits, double z0, double &chi);
0043
0044 const edm::EDGetTokenT<SiPixelRecHitCollection> srcPixels_;
0045 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerToken_;
0046
0047 const double minZ_;
0048 const double maxZ_;
0049 const double zStep_;
0050 const int maxHists_;
0051
0052 edm::Service<TFileService> fs;
0053 int counter;
0054 };
0055
0056
0057 HIPixelClusterVtxAnalyzer::HIPixelClusterVtxAnalyzer(const edm::ParameterSet &ps)
0058 : srcPixels_(consumes<SiPixelRecHitCollection>(ps.getParameter<edm::InputTag>("pixelRecHits"))),
0059 trackerToken_(esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()),
0060 minZ_(ps.getParameter<double>("minZ")),
0061 maxZ_(ps.getParameter<double>("maxZ")),
0062 zStep_(ps.getParameter<double>("zStep")),
0063 maxHists_(ps.getParameter<int>("maxHists")),
0064 counter(0)
0065
0066 {
0067
0068 usesResource(TFileService::kSharedResource);
0069 }
0070
0071
0072 void HIPixelClusterVtxAnalyzer::analyze(const edm::Event &ev, const edm::EventSetup &es) {
0073 if (counter > maxHists_)
0074 return;
0075 std::cout << "counter = " << counter << std::endl;
0076 counter++;
0077
0078 edm::EventNumber_t evtnum = ev.id().event();
0079 TH1D *hClusterVtx = fs->make<TH1D>(Form("hClusterVtx_%llu", evtnum),
0080 "compatibility of pixel cluster length with vertex hypothesis; z [cm]",
0081 (int)((maxZ_ - minZ_) / zStep_),
0082 minZ_,
0083 maxZ_);
0084
0085
0086 edm::Handle<SiPixelRecHitCollection> hRecHits;
0087 try {
0088 ev.getByToken(srcPixels_, hRecHits);
0089 } catch (...) {
0090 }
0091
0092
0093 if (hRecHits.isValid()) {
0094 const auto &trackerHandle = es.getHandle(trackerToken_);
0095 const TrackerGeometry *tgeo = trackerHandle.product();
0096 const SiPixelRecHitCollection *hits = hRecHits.product();
0097
0098
0099 std::vector<VertexHit> vhits;
0100 for (auto const &hit : hits->data()) {
0101 if (!hit.isValid())
0102 continue;
0103 DetId id(hit.geographicalId());
0104 if (id.subdetId() != int(PixelSubdetector::PixelBarrel))
0105 continue;
0106 const PixelGeomDetUnit *pgdu = static_cast<const PixelGeomDetUnit *>(tgeo->idToDet(id));
0107 if (1) {
0108 const RectangularPixelTopology *pixTopo =
0109 static_cast<const RectangularPixelTopology *>(&(pgdu->specificTopology()));
0110 std::vector<SiPixelCluster::Pixel> pixels(hit.cluster()->pixels());
0111 bool pixelOnEdge = false;
0112 for (std::vector<SiPixelCluster::Pixel>::const_iterator pixel = pixels.begin(); pixel != pixels.end();
0113 ++pixel) {
0114 int pixelX = pixel->x;
0115 int pixelY = pixel->y;
0116 if (pixTopo->isItEdgePixelInX(pixelX) || pixTopo->isItEdgePixelInY(pixelY)) {
0117 pixelOnEdge = true;
0118 break;
0119 }
0120 }
0121 if (pixelOnEdge)
0122 continue;
0123 }
0124
0125 LocalPoint lpos = LocalPoint(hit.localPosition().x(), hit.localPosition().y(), hit.localPosition().z());
0126 GlobalPoint gpos = pgdu->toGlobal(lpos);
0127 VertexHit vh;
0128 vh.z = gpos.z();
0129 vh.r = gpos.perp();
0130 vh.w = hit.cluster()->sizeY();
0131 vhits.push_back(vh);
0132 }
0133
0134
0135 double zest = 0.0;
0136 int nhits = 0, nhits_max = 0;
0137 double chi = 0, chi_max = 1e+9;
0138 for (double z0 = minZ_; z0 <= maxZ_; z0 += zStep_) {
0139 nhits = getContainedHits(vhits, z0, chi);
0140 hClusterVtx->Fill(z0, nhits);
0141 if (nhits == 0)
0142 continue;
0143 if (nhits > nhits_max) {
0144 chi_max = 1e+9;
0145 nhits_max = nhits;
0146 }
0147 if (nhits >= nhits_max && chi < chi_max) {
0148 chi_max = chi;
0149 zest = z0;
0150 }
0151 }
0152
0153 LogTrace("MinBiasTracking") << " [vertex position] estimated = " << zest
0154 << " | pixel barrel hits = " << vhits.size();
0155 }
0156 }
0157
0158
0159 int HIPixelClusterVtxAnalyzer::getContainedHits(const std::vector<VertexHit> &hits, double z0, double &chi) {
0160
0161 int n = 0;
0162 chi = 0.;
0163
0164 for (std::vector<VertexHit>::const_iterator hit = hits.begin(); hit != hits.end(); hit++) {
0165 double p = 2 * fabs(hit->z - z0) / hit->r + 0.5;
0166 if (TMath::Abs(p - hit->w) <= 1.) {
0167 chi += fabs(p - hit->w);
0168 n++;
0169 }
0170 }
0171 return n;
0172 }
0173
0174 #include "FWCore/Framework/interface/MakerMacros.h"
0175 DEFINE_FWK_MODULE(HIPixelClusterVtxAnalyzer);