CMS 3D CMS Logo

PixelVertexProducerMedian.cc
Go to the documentation of this file.
2 
7 
10 
13 
14 #include <fstream>
15 #include <iostream>
16 #include <vector>
17 #include <algorithm>
18 
19 #include "TROOT.h"
20 #include "TH1F.h"
21 #include "TF1.h"
22 
23 /*****************************************************************************/
24 struct ComparePairs {
25  bool operator()(const reco::Track* t1, const reco::Track* t2) { return (t1->vz() < t2->vz()); };
26 };
27 
28 /*****************************************************************************/
30  produces<reco::VertexCollection>();
31 }
32 
33 /*****************************************************************************/
35 
36 /*****************************************************************************/
38  // Get pixel tracks
40  std::string trackCollectionName = theConfig.getParameter<std::string>("TrackCollection");
41  ev.getByLabel(trackCollectionName, trackCollection);
42  const reco::TrackCollection tracks_ = *(trackCollection.product());
43 
44  thePtMin = theConfig.getParameter<double>("PtMin");
45 
46  // Select tracks
47  std::vector<const reco::Track*> tracks;
48  for (unsigned int i = 0; i < tracks_.size(); i++) {
49  if (tracks_[i].pt() > thePtMin) {
50  reco::TrackRef recTrack(trackCollection, i);
51  tracks.push_back(&(*recTrack));
52  }
53  }
54 
55  LogTrace("MinBiasTracking") << " [VertexProducer] selected tracks: " << tracks.size() << " (out of " << tracks_.size()
56  << ")";
57 
58  auto vertices = std::make_unique<reco::VertexCollection>();
59 
60  if (!tracks.empty()) {
61  // Sort along vertex z position
62  std::sort(tracks.begin(), tracks.end(), ComparePairs());
63 
64  // Median
65  float med;
66  if (tracks.size() % 2 == 0)
67  med = (tracks[tracks.size() / 2 - 1]->vz() + tracks[tracks.size() / 2]->vz()) / 2;
68  else
69  med = tracks[tracks.size() / 2]->vz();
70 
71  LogTrace("MinBiasTracking") << " [vertex position] median = " << med << " cm";
72 
73  if (tracks.size() > 10) {
74  // Binning around med, halfWidth
75  int nBin = 100;
76  float halfWidth = 0.1; // cm
77 
78  // Most probable
79  TH1F histo("histo", "histo", nBin, -halfWidth, halfWidth);
80 
81  for (std::vector<const reco::Track*>::const_iterator track = tracks.begin(); track != tracks.end(); track++)
82  if (fabs((*track)->vz() - med) < halfWidth)
83  histo.Fill((*track)->vz() - med);
84 
85  LogTrace("MinBiasTracking") << " [vertex position] most prob = "
86  << med + histo.GetBinCenter(histo.GetMaximumBin()) << " cm";
87 
88  // Fit above max/2
89  histo.Sumw2();
90 
91  TF1 f1("f1", "[0]*exp(-0.5 * ((x-[1])/[2])^2) + [3]");
92  f1.SetParameters(10., 0., 0.01, 1.);
93 
94  histo.Fit(&f1, "QN");
95 
96  LogTrace("MinBiasTracking") << " [vertex position] fitted = " << med + f1.GetParameter(1) << " +- "
97  << f1.GetParError(1) << " cm";
98 
99  // Store
101  err(2, 2) = f1.GetParError(1) * f1.GetParError(1);
102  reco::Vertex ver(reco::Vertex::Point(0, 0, med + f1.GetParameter(1)), err, 0, 1, 1);
103  vertices->push_back(ver);
104  } else {
105  // Store
107  err(2, 2) = 0.1 * 0.1;
108  reco::Vertex ver(reco::Vertex::Point(0, 0, med), err, 0, 1, 1);
109  vertices->push_back(ver);
110  }
111  }
112  ev.put(std::move(vertices));
113 }
T getParameter(std::string const &) const
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:131
std::vector< Track > TrackCollection
collection of Tracks
Definition: TrackFwd.h:14
math::Error< dimension >::type Error
covariance error matrix (3x3)
Definition: Vertex.h:43
bool ev
PixelVertexProducerMedian(const edm::ParameterSet &ps)
math::XYZPoint Point
point in the space
Definition: Vertex.h:39
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:488
#define LogTrace(id)
double vz() const
z coordinate of the reference point on track
Definition: TrackBase.h:626
T const * product() const
Definition: Handle.h:69
void produce(edm::Event &ev, const edm::EventSetup &es) override
bool operator()(const reco::Track *t1, const reco::Track *t2)
def move(src, dest)
Definition: eostools.py:511