|
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
|
#include "EgammaAnalysis/ElectronTools/interface/ElectronEPcombinator.h"
#include <iostream>
//Accessor to the combination results
void ElectronEPcombinator::combine(SimpleElectron& electron) {
electron_ = electron;
computeEPcombination();
electron.setCombinedMomentum(combinedMomentum_);
electron.setCombinedMomentumError(combinedMomentumError_);
}
//Core code to compute the EP combination
void ElectronEPcombinator::computeEPcombination() {
if (mode_ == 1) {
scEnergy_ = electron_.getNewEnergy();
scEnergyError_ = electron_.getNewEnergyError();
}
if (mode_ == 2) {
scEnergy_ = electron_.getRegEnergy();
scEnergyError_ = electron_.getRegEnergyError();
}
trackerMomentum_ = electron_.getTrackerMomentum();
trackerMomentumError_ = electron_.getTrackerMomentumError();
elClass_ = electron_.getElClass();
combinedMomentum_ = scEnergy_; // initial
combinedMomentumError_ = 999.;
// first check for large errors
if (trackerMomentumError_ / trackerMomentum_ > 0.5 && scEnergyError_ / scEnergy_ <= 0.5) {
combinedMomentum_ = scEnergy_;
combinedMomentumError_ = scEnergyError_;
} else if (trackerMomentumError_ / trackerMomentum_ <= 0.5 && scEnergyError_ / scEnergy_ > 0.5) {
combinedMomentum_ = trackerMomentum_;
combinedMomentumError_ = trackerMomentumError_;
} else if (trackerMomentumError_ / trackerMomentum_ > 0.5 && scEnergyError_ / scEnergy_ > 0.5) {
if (trackerMomentumError_ / trackerMomentum_ < scEnergyError_ / scEnergy_) {
combinedMomentum_ = trackerMomentum_;
combinedMomentumError_ = trackerMomentumError_;
} else {
combinedMomentum_ = scEnergy_;
combinedMomentumError_ = scEnergyError_;
}
}
// then apply the combination algorithm
else {
// calculate E/p and corresponding error
float eOverP = scEnergy_ / trackerMomentum_;
float errorEOverP = sqrt((scEnergyError_ / trackerMomentum_) * (scEnergyError_ / trackerMomentum_) +
(scEnergy_ * trackerMomentumError_ / trackerMomentum_ / trackerMomentum_) *
(scEnergy_ * trackerMomentumError_ / trackerMomentum_ / trackerMomentum_));
bool eleIsNotInCombination = false;
if ((eOverP > 1 + 2.5 * errorEOverP) || (eOverP < 1 - 2.5 * errorEOverP) || (eOverP < 0.8) || (eOverP > 1.3)) {
eleIsNotInCombination = true;
}
if (eleIsNotInCombination) {
if (eOverP > 1) {
combinedMomentum_ = scEnergy_;
combinedMomentumError_ = scEnergyError_;
} else {
if (elClass_ == 0) // == reco::GsfElectron::GOLDEN)
{
combinedMomentum_ = scEnergy_;
combinedMomentumError_ = scEnergyError_;
}
if (elClass_ == 1) //reco::GsfElectron::BIGBREM)
{
if (scEnergy_ < 36) {
combinedMomentum_ = trackerMomentum_;
combinedMomentumError_ = trackerMomentumError_;
} else {
combinedMomentum_ = scEnergy_;
combinedMomentumError_ = scEnergyError_;
}
}
if (elClass_ == 2) // == reco::GsfElectron::BADTRACK)
{
combinedMomentum_ = scEnergy_;
combinedMomentumError_ = scEnergyError_;
}
if (elClass_ == 3) //reco::GsfElectron::SHOWERING)
{
if (scEnergy_ < 30) {
combinedMomentum_ = trackerMomentum_;
combinedMomentumError_ = trackerMomentumError_;
} else {
combinedMomentum_ = scEnergy_;
combinedMomentumError_ = scEnergyError_;
}
}
if (elClass_ == 4) //reco::GsfElectron::GAP)
{
if (scEnergy_ < 60) {
combinedMomentum_ = trackerMomentum_;
combinedMomentumError_ = trackerMomentumError_;
} else {
combinedMomentum_ = scEnergy_;
combinedMomentumError_ = scEnergyError_;
}
}
}
} else {
// combination
combinedMomentum_ = (scEnergy_ / scEnergyError_ / scEnergyError_ +
trackerMomentum_ / trackerMomentumError_ / trackerMomentumError_) /
(1 / scEnergyError_ / scEnergyError_ + 1 / trackerMomentumError_ / trackerMomentumError_);
float combinedMomentum_Variance =
1 / (1 / scEnergyError_ / scEnergyError_ + 1 / trackerMomentumError_ / trackerMomentumError_);
combinedMomentumError_ = sqrt(combinedMomentum_Variance);
}
}
}
|