PixelLowPtUtilities/plugins/PixelVertexProducerMedian.cc

0001 #include "PixelVertexProducerMedian.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/TrackReco/interface/Track.h"
0009 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0010
0011 #include "DataFormats/VertexReco/interface/Vertex.h"
0012 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0013
0014 #include <fstream>
0015 #include <iostream>
0016 #include <vector>
0017 #include <algorithm>
0018
0019 #include "TROOT.h"
0020 #include "TH1F.h"
0021 #include "TF1.h"
0022
0023 /*****************************************************************************/
0024 struct ComparePairs {
0025   bool operator()(const reco::Track* t1, const reco::Track* t2) { return (t1->vz() < t2->vz()); };
0026 };
0027
0028 /*****************************************************************************/
0029 PixelVertexProducerMedian::PixelVertexProducerMedian(const edm::ParameterSet& ps) : theConfig(ps) {
0030   thePtMin = theConfig.getParameter<double>("PtMin");
0031   produces<reco::VertexCollection>();
0032 }
0033
0034 /*****************************************************************************/
0035 PixelVertexProducerMedian::~PixelVertexProducerMedian() {}
0036
0037 /*****************************************************************************/
0038 void PixelVertexProducerMedian::produce(edm::StreamID, edm::Event& ev, const edm::EventSetup& es) const {
0039   // Get pixel tracks
0040   edm::Handle<reco::TrackCollection> trackCollection;
0041   std::string trackCollectionName = theConfig.getParameter<std::string>("TrackCollection");
0042   ev.getByLabel(trackCollectionName, trackCollection);
0043   const reco::TrackCollection tracks_ = *(trackCollection.product());
0044
0045   // Select tracks
0046   std::vector<const reco::Track*> tracks;
0047   for (unsigned int i = 0; i < tracks_.size(); i++) {
0048     if (tracks_[i].pt() > thePtMin) {
0049       reco::TrackRef recTrack(trackCollection, i);
0050       tracks.push_back(&(*recTrack));
0051     }
0052   }
0053
0054   LogTrace("MinBiasTracking") << " [VertexProducer] selected tracks: " << tracks.size() << " (out of " << tracks_.size()
0055                               << ")";
0056
0057   auto vertices = std::make_unique<reco::VertexCollection>();
0058
0059   if (!tracks.empty()) {
0060     // Sort along vertex z position
0061     std::sort(tracks.begin(), tracks.end(), ComparePairs());
0062
0063     // Median
0064     float med;
0065     if (tracks.size() % 2 == 0)
0066       med = (tracks[tracks.size() / 2 - 1]->vz() + tracks[tracks.size() / 2]->vz()) / 2;
0067     else
0068       med = tracks[tracks.size() / 2]->vz();
0069
0070     LogTrace("MinBiasTracking") << "  [vertex position] median    = " << med << " cm";
0071
0072     if (tracks.size() > 10) {
0073       // Binning around med, halfWidth
0074       int nBin = 100;
0075       float halfWidth = 0.1;  // cm
0076
0077       // Most probable
0078       TH1F histo("histo", "histo", nBin, -halfWidth, halfWidth);
0079
0080       for (std::vector<const reco::Track*>::const_iterator track = tracks.begin(); track != tracks.end(); track++)
0081         if (fabs((*track)->vz() - med) < halfWidth)
0082           histo.Fill((*track)->vz() - med);
0083
0084       LogTrace("MinBiasTracking") << "  [vertex position] most prob = "
0085                                   << med + histo.GetBinCenter(histo.GetMaximumBin()) << " cm";
0086
0087       // Fit above max/2
0088       histo.Sumw2();
0089
0090       TF1 f1("f1", "[0]*exp(-0.5 * ((x-[1])/[2])^2) + [3]");
0091       f1.SetParameters(10., 0., 0.01, 1.);
0092
0093       histo.Fit(&f1, "QN");
0094
0095       LogTrace("MinBiasTracking") << "  [vertex position] fitted    = " << med + f1.GetParameter(1) << " +- "
0096                                   << f1.GetParError(1) << " cm";
0097
0098       // Store
0099       reco::Vertex::Error err;
0100       err(2, 2) = f1.GetParError(1) * f1.GetParError(1);
0101       reco::Vertex ver(reco::Vertex::Point(0, 0, med + f1.GetParameter(1)), err, 0, 1, 1);
0102       vertices->push_back(ver);
0103     } else {
0104       // Store
0105       reco::Vertex::Error err;
0106       err(2, 2) = 0.1 * 0.1;
0107       reco::Vertex ver(reco::Vertex::Point(0, 0, med), err, 0, 1, 1);
0108       vertices->push_back(ver);
0109     }
0110   }
0111   ev.put(std::move(vertices));
0112 }