HiTracking/test/HIPixelClusterVtxAnalyzer.cc

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 // ROOT includes
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_;  //pixel rec hits
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   // Constructor
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   // get pixel rechits
0086   edm::Handle<SiPixelRecHitCollection> hRecHits;
0087   try {
0088     ev.getByToken(srcPixels_, hRecHits);
0089   } catch (...) {
0090   }
0091
0092   // get tracker geometry
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     // loop over pixel rechits
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     // estimate z-position from cluster lengths
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   // Calculate number of hits contained in v-shaped window in cluster y-width vs. z-position.
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;  // FIXME
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);