File indexing completed on 2024-04-06 12:25:20
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
0033
0034 TMinuitMinimizer::UseStaticMinuit(false);
0035 }
0036
0037
0038 void HIPixelMedianVtxProducer::produce(edm::Event& ev, const edm::EventSetup& es) {
0039
0040 edm::Handle<reco::TrackCollection> trackCollection;
0041 ev.getByToken(theTrackCollection, trackCollection);
0042 const reco::TrackCollection tracks_ = *(trackCollection.product());
0043
0044
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
0057 auto vertices = std::make_unique<reco::VertexCollection>();
0058
0059
0060 if (tracks.empty()) {
0061 ev.put(std::move(vertices));
0062 return;
0063 }
0064
0065
0066 std::sort(tracks.begin(), tracks.end(), ComparePairs());
0067
0068
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
0078 if (tracks.size() > thePeakFindThresh) {
0079
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
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
0109 med = nBinAvg;
0110
0111 LogTrace("MinBiasTracking") << " [vertex position] 3-bin weighted average = " << nBinAvg << " cm";
0112 }
0113
0114
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
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
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 }