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/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/ServiceRegistry/interface/Service.h"
0006 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0007 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0008 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0009
0010 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0011 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.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 PixelTopology *pixTopo = static_cast<const PixelTopology *>(&(pgdu->specificTopology()));
0109         std::vector<SiPixelCluster::Pixel> pixels(hit.cluster()->pixels());
0110         bool pixelOnEdge = false;
0111         for (std::vector<SiPixelCluster::Pixel>::const_iterator pixel = pixels.begin(); pixel != pixels.end();
0112              ++pixel) {
0113           int pixelX = pixel->x;
0114           int pixelY = pixel->y;
0115           if (pixTopo->isItEdgePixelInX(pixelX) || pixTopo->isItEdgePixelInY(pixelY)) {
0116             pixelOnEdge = true;
0117             break;
0118           }
0119         }
0120         if (pixelOnEdge)
0121           continue;
0122       }
0123
0124       LocalPoint lpos = LocalPoint(hit.localPosition().x(), hit.localPosition().y(), hit.localPosition().z());
0125       GlobalPoint gpos = pgdu->toGlobal(lpos);
0126       VertexHit vh;
0127       vh.z = gpos.z();
0128       vh.r = gpos.perp();
0129       vh.w = hit.cluster()->sizeY();
0130       vhits.push_back(vh);
0131     }
0132
0133     // estimate z-position from cluster lengths
0134     double zest = 0.0;
0135     int nhits = 0, nhits_max = 0;
0136     double chi = 0, chi_max = 1e+9;
0137     for (double z0 = minZ_; z0 <= maxZ_; z0 += zStep_) {
0138       nhits = getContainedHits(vhits, z0, chi);
0139       hClusterVtx->Fill(z0, nhits);
0140       if (nhits == 0)
0141         continue;
0142       if (nhits > nhits_max) {
0143         chi_max = 1e+9;
0144         nhits_max = nhits;
0145       }
0146       if (nhits >= nhits_max && chi < chi_max) {
0147         chi_max = chi;
0148         zest = z0;
0149       }
0150     }
0151
0152     LogTrace("MinBiasTracking") << "  [vertex position] estimated = " << zest
0153                                 << " | pixel barrel hits = " << vhits.size();
0154   }
0155 }
0156
0157 /*****************************************************************************/
0158 int HIPixelClusterVtxAnalyzer::getContainedHits(const std::vector<VertexHit> &hits, double z0, double &chi) {
0159   // Calculate number of hits contained in v-shaped window in cluster y-width vs. z-position.
0160   int n = 0;
0161   chi = 0.;
0162
0163   for (std::vector<VertexHit>::const_iterator hit = hits.begin(); hit != hits.end(); hit++) {
0164     double p = 2 * fabs(hit->z - z0) / hit->r + 0.5;  // FIXME
0165     if (TMath::Abs(p - hit->w) <= 1.) {
0166       chi += fabs(p - hit->w);
0167       n++;
0168     }
0169   }
0170   return n;
0171 }
0172
0173 #include "FWCore/Framework/interface/MakerMacros.h"
0174 DEFINE_FWK_MODULE(HIPixelClusterVtxAnalyzer);