-
Notifications
You must be signed in to change notification settings - Fork 4.3k
/
PVClusterComparer.cc
102 lines (90 loc) · 3.57 KB
/
PVClusterComparer.cc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
#include "RecoPixelVertexing/PixelVertexFinding/interface/PVClusterComparer.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "TMath.h"
PVClusterComparer::PVClusterComparer()
: track_pT_min_(2.5), track_pT_max_(10.), track_chi2_max_(9999999.), track_prob_min_(-1.) {
setChisquareQuantile();
}
PVClusterComparer::PVClusterComparer(double track_pt_min,
double track_pt_max,
double track_chi2_max,
double track_prob_min)
: track_pT_min_(track_pt_min),
track_pT_max_(track_pt_max),
track_chi2_max_(track_chi2_max),
track_prob_min_(track_prob_min) {
setChisquareQuantile();
}
double PVClusterComparer::pTSquaredSum(const PVCluster &v) {
double sum = 0;
for (size_t i = 0; i < v.tracks().size(); ++i) {
double pt = v.tracks()[i]->pt();
if (pt < track_pT_min_)
continue; // Don't count tracks below track_pT_min_ (2.5 GeV)
// RM : exclude badly reconstructed tracks from the sum
// if (track_prob_min_ >= 0. && track_prob_min_ <= 1.)
// if (TMath::Prob(v.tracks()[i]->chi2(),v.tracks()[i]->ndof()) < track_prob_min_) continue ;
if (track_prob_min_ >= 0. && track_prob_min_ <= 1.) {
size_t ndof = v.tracks()[i]->ndof();
if (ndof >= maxChi2_.size())
updateChisquareQuantile(ndof);
// cut on chi2 which corresponds to the configured probability
if (v.tracks()[i]->chi2() > maxChi2_[ndof])
continue;
}
if (v.tracks()[i]->normalizedChi2() > track_chi2_max_)
continue;
if (pt > track_pT_max_)
pt = track_pT_max_;
sum += pt * pt;
}
return sum;
}
double PVClusterComparer::pTSquaredSum(const reco::Vertex &v) {
double sum = 0;
for (reco::Vertex::trackRef_iterator i = v.tracks_begin(), ie = v.tracks_end(); i != ie; ++i) {
double pt = (*i)->pt();
if (pt < track_pT_min_)
continue; // Don't count tracks below track_pT_min_ (2.5 GeV)
// RM : exclude badly reconstructed tracks from the sum
// if (track_prob_min_ >= 0. && track_prob_min_ <= 1.)
// if (TMath::Prob((*i)->chi2(),(*i)->ndof()) < track_prob_min_) continue ;
if (track_prob_min_ >= 0. && track_prob_min_ <= 1.) {
unsigned int ndof = (*i)->ndof();
if (ndof >= maxChi2_.size())
updateChisquareQuantile(ndof);
// cut on chi2 which corresponds to the configured probability
if ((*i)->chi2() > maxChi2_[ndof])
continue;
}
if ((*i)->normalizedChi2() > track_chi2_max_)
continue;
if (pt > track_pT_max_)
pt = track_pT_max_;
sum += pt * pt;
}
return sum;
}
void PVClusterComparer::setChisquareQuantile() {
std::vector<double> maxChi2(20, 0.);
if (track_prob_min_ >= 0. && track_prob_min_ <= 1.)
for (size_t ndof = 0; ndof < maxChi2.size(); ++ndof)
// http://root.cern.ch/root/html/TMath.html#TMath:ChisquareQuantile
maxChi2[ndof] = TMath::ChisquareQuantile(1 - track_prob_min_, ndof);
maxChi2_ = maxChi2;
}
void PVClusterComparer::updateChisquareQuantile(size_t ndof) {
size_t oldsize = maxChi2_.size();
// maxChi2_.resize(ndof+1);
for (size_t i = oldsize; i <= ndof; ++i) {
double chi2 = TMath::ChisquareQuantile(1 - track_prob_min_, i);
maxChi2_.push_back(chi2);
}
}
bool PVClusterComparer::operator()(const PVCluster &v1, const PVCluster &v2) {
return (pTSquaredSum(v1) > pTSquaredSum(v2));
}
bool PVClusterComparer::operator()(const reco::Vertex &v1, const reco::Vertex &v2) {
return (pTSquaredSum(v1) > pTSquaredSum(v2));
}