|
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
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
#ifndef BTauReco_BJetTagTrackProbability_h
#define BTauReco_BJetTagTrackProbability_h
#include "DataFormats/BTauReco/interface/RefMacros.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/BTauReco/interface/JTATagInfo.h"
#include "DataFormats/JetReco/interface/JetTracksAssociation.h"
namespace reco {
class TrackProbabilityTagInfo : public JTATagInfo {
public:
TrackProbabilityTagInfo(const std::vector<double>& probability2d,
const std::vector<double>& probability3d,
const std::vector<int>& trackOrder2d,
const std::vector<int>& trackOrder3d,
const JetTracksAssociationRef& jtaRef)
: JTATagInfo(jtaRef),
m_probability2d(probability2d),
m_probability3d(probability3d),
m_trackOrder2d(trackOrder2d),
m_trackOrder3d(trackOrder3d) {}
TrackProbabilityTagInfo() {}
~TrackProbabilityTagInfo() override {}
int factorial(int n) const {
if (n < 2)
return 1;
else
return n * factorial(n - 1);
}
virtual float probability(size_t n, int ip) const {
if (ip == 0) {
if (n < m_probability3d.size())
return m_probability3d[n];
} else {
if (n < m_probability2d.size())
return m_probability2d[n];
}
return -10.;
}
virtual float jetProbability(int ip, float minTrackProb) const {
const std::vector<double>* vp;
if (ip == 0)
vp = &m_probability3d;
else
vp = &m_probability2d;
const std::vector<double>& v = *vp;
int ngoodtracks = v.size();
double SumJet = 0.;
for (std::vector<double>::const_iterator q = v.begin(); q != v.end(); q++) {
SumJet += (*q > minTrackProb) ? log(*q) : log(minTrackProb);
}
double ProbJet;
double Loginvlog = 0;
if (SumJet < 0.) {
if (ngoodtracks >= 2) {
Loginvlog = log(-SumJet);
}
double Prob = 1.;
for (int l = 1; l != ngoodtracks; l++) {
Prob += exp(l * Loginvlog - log(1. * factorial(l)));
}
double LogProb = log(Prob);
ProbJet = std::min(exp(std::max(LogProb + SumJet, -30.)), 1.);
} else {
ProbJet = 1.;
}
if (ProbJet > 1)
std::cout << "ProbJet too high: " << ProbJet << std::endl;
//double LogProbJet=-log(ProbJet);
//return 1.-ProbJet;
return -log10(ProbJet) / 4.;
}
/**
Recompute discriminator
ipType = 0 means 3d impact parameter
ipType = 1 means transverse impact parameter
minProb is the minimum probability allowed for a single track. Tracks with lower probability
are considered with a probability = minProb.
*/
virtual float discriminator(int ipType, float minProb) const { return jetProbability(ipType, minProb); }
virtual int selectedTracks(int ipType) const {
if (ipType == 0)
return m_probability3d.size();
else
return m_probability2d.size();
}
virtual int trackIndex(size_t n, int ip) const {
if (ip == 0) {
if (n < m_probability3d.size())
return m_trackOrder3d[n];
} else {
if (n < m_probability2d.size())
return m_trackOrder2d[n];
}
return 0;
}
virtual const Track& track(size_t n, int ipType) const { return *tracks()[trackIndex(n, ipType)]; }
TrackProbabilityTagInfo* clone() const override { return new TrackProbabilityTagInfo(*this); }
private:
std::vector<double> m_probability2d; //
std::vector<double> m_probability3d; // create a smarter container instead of
std::vector<int> m_trackOrder2d; // this pair of vectors.
std::vector<int> m_trackOrder3d; //
};
//typedef edm::ExtCollection< TrackProbabilityTagInfo,JetTagCollection> TrackProbabilityExtCollection;
//typedef edm::OneToOneAssociation<JetTagCollection, TrackProbabilityTagInfo> TrackProbabilityExtCollection;
DECLARE_EDM_REFS(TrackProbabilityTagInfo)
} // namespace reco
#endif
|