HiTracking/plugins/HIPixelMedianVtxProducer.cc

0001 #include "HIPixelMedianVtxProducer.h"
0002
0003 #include "FWCore/Framework/interface/MakerMacros.h"
0004 #include "FWCore/Framework/interface/Frameworkfwd.h"
0005 #include "FWCore/Framework/interface/Event.h"
0006
0007 #include "DataFormats/VertexReco/interface/Vertex.h"
0008 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0009
0010 #include <fstream>
0011 #include <iostream>
0012 #include <vector>
0013 #include <algorithm>
0014
0015 #include "TROOT.h"
0016 #include "TH1F.h"
0017 #include "TF1.h"
0018 #include "TMinuitMinimizer.h"
0019
0020 /*****************************************************************************/
0021 HIPixelMedianVtxProducer::HIPixelMedianVtxProducer(const edm::ParameterSet& ps)
0022     : theTrackCollection(consumes<reco::TrackCollection>(ps.getParameter<edm::InputTag>("TrackCollection"))),
0023       thePtMin(ps.getParameter<double>("PtMin")),
0024       thePeakFindThresh(ps.getParameter<unsigned int>("PeakFindThreshold")),
0025       thePeakFindMaxZ(ps.getParameter<double>("PeakFindMaxZ")),
0026       thePeakFindBinning(ps.getParameter<int>("PeakFindBinsPerCm")),
0027       theFitThreshold(ps.getParameter<int>("FitThreshold")),
0028       theFitMaxZ(ps.getParameter<double>("FitMaxZ")),
0029       theFitBinning(ps.getParameter<int>("FitBinsPerCm")) {
0030   produces<reco::VertexCollection>();
0031
0032   //In order to make fitting ROOT histograms thread safe
0033   // one must call this undocumented function
0034   TMinuitMinimizer::UseStaticMinuit(false);
0035 }
0036
0037 /*****************************************************************************/
0038 void HIPixelMedianVtxProducer::produce(edm::Event& ev, const edm::EventSetup& es) {
0039   // Get pixel tracks
0040   edm::Handle<reco::TrackCollection> trackCollection;
0041   ev.getByToken(theTrackCollection, trackCollection);
0042   const reco::TrackCollection tracks_ = *(trackCollection.product());
0043
0044   // Select tracks above minimum pt
0045   std::vector<const reco::Track*> tracks;
0046   for (unsigned int i = 0; i < tracks_.size(); i++) {
0047     if (tracks_[i].pt() < thePtMin && std::fabs(tracks_[i].vz()) < 100000.)
0048       continue;
0049     reco::TrackRef recTrack(trackCollection, i);
0050     tracks.push_back(&(*recTrack));
0051   }
0052
0053   LogTrace("MinBiasTracking") << " [VertexProducer] selected tracks: " << tracks.size() << " (out of " << tracks_.size()
0054                               << ") above pt = " << thePtMin;
0055
0056   // Output vertex collection
0057   auto vertices = std::make_unique<reco::VertexCollection>();
0058
0059   // No tracks -> return empty collection
0060   if (tracks.empty()) {
0061     ev.put(std::move(vertices));
0062     return;
0063   }
0064
0065   // Sort tracks according to vertex z position
0066   std::sort(tracks.begin(), tracks.end(), ComparePairs());
0067
0068   // Calculate median vz
0069   float med;
0070   if (tracks.size() % 2 == 0)
0071     med = (tracks[tracks.size() / 2 - 1]->vz() + tracks[tracks.size() / 2]->vz()) / 2;
0072   else
0073     med = tracks[tracks.size() / 2]->vz();
0074
0075   LogTrace("MinBiasTracking") << "  [vertex position] median    = " << med << " cm";
0076
0077   // In high multiplicity events, fit around most probable position
0078   if (tracks.size() > thePeakFindThresh) {
0079     // Find maximum bin
0080     TH1F hmax("hmax", "hmax", thePeakFindBinning * 2.0 * thePeakFindMaxZ, -1.0 * thePeakFindMaxZ, thePeakFindMaxZ);
0081
0082     for (std::vector<const reco::Track*>::const_iterator track = tracks.begin(); track != tracks.end(); track++)
0083       if (fabs((*track)->vz()) < thePeakFindMaxZ)
0084         hmax.Fill((*track)->vz());
0085
0086     int maxBin = hmax.GetMaximumBin();
0087
0088     LogTrace("MinBiasTracking") << "  [vertex position] most prob = " << hmax.GetBinCenter(maxBin) << " cm";
0089
0090     // Find 3-bin weighted average
0091     float num = 0.0, denom = 0.0;
0092     for (int i = -1; i <= 1; i++) {
0093       num += hmax.GetBinContent(maxBin + i) * hmax.GetBinCenter(maxBin + i);
0094       denom += hmax.GetBinContent(maxBin + i);
0095     }
0096
0097     if (denom == 0) {
0098       reco::Vertex::Error err;
0099       err(2, 2) = 0.1 * 0.1;
0100       reco::Vertex ver(reco::Vertex::Point(0, 0, 99999.), err, 0, 0, 1);
0101       vertices->push_back(ver);
0102       ev.put(std::move(vertices));
0103       return;
0104     }
0105
0106     float nBinAvg = num / denom;
0107
0108     // Center fit at 3-bin weighted average around max bin
0109     med = nBinAvg;
0110
0111     LogTrace("MinBiasTracking") << "  [vertex position] 3-bin weighted average = " << nBinAvg << " cm";
0112   }
0113
0114   // Bin vz-values around most probable value (or median) for fitting
0115   TH1F histo("histo", "histo", theFitBinning * 2.0 * theFitMaxZ, -1.0 * theFitMaxZ, theFitMaxZ);
0116   histo.Sumw2();
0117
0118   for (std::vector<const reco::Track*>::const_iterator track = tracks.begin(); track != tracks.end(); track++)
0119     if (fabs((*track)->vz() - med) < theFitMaxZ)
0120       histo.Fill((*track)->vz() - med);
0121
0122   LogTrace("MinBiasTracking") << "  [vertex position] most prob for fit = "
0123                               << med + histo.GetBinCenter(histo.GetMaximumBin()) << " cm";
0124
0125   // If there are very few entries, don't do the fit
0126   if (histo.GetEntries() <= theFitThreshold) {
0127     LogTrace("MinBiasTracking") << "  [vertex position] Fewer than" << theFitThreshold
0128                                 << " entries in fit histogram. Returning median.";
0129
0130     reco::Vertex::Error err;
0131     err(2, 2) = 0.1 * 0.1;
0132     reco::Vertex ver(reco::Vertex::Point(0, 0, med), err, 0, 1, 1);
0133     vertices->push_back(ver);
0134     ev.put(std::move(vertices));
0135     return;
0136   }
0137
0138   // Otherwise, there are enough entries to refine the estimate with a fit
0139   TF1 f1("f1", "[0]*exp(-0.5 * ((x-[1])/[2])^2) + [3]");
0140   f1.SetParameters(10., 0., 0.02, 0.002 * tracks.size());
0141   f1.SetParLimits(1, -0.1, 0.1);
0142   f1.SetParLimits(2, 0.001, 0.05);
0143   f1.SetParLimits(3, 0.0, 0.005 * tracks.size());
0144
0145   histo.Fit(&f1, "QN SERIAL");
0146
0147   LogTrace("MinBiasTracking") << "  [vertex position] fitted    = " << med + f1.GetParameter(1) << " +- "
0148                               << f1.GetParError(1) << " cm";
0149
0150   reco::Vertex::Error err;
0151   err(2, 2) = f1.GetParError(1) * f1.GetParError(1);
0152   reco::Vertex ver(reco::Vertex::Point(0, 0, med + f1.GetParameter(1)), err, 0, 1, 1);
0153   vertices->push_back(ver);
0154
0155   ev.put(std::move(vertices));
0156   return;
0157 }