|
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
|
#include "EgammaAnalysis/ElectronTools/interface/EpCombinationTool.h"
#include "CondFormats/GBRForest/interface/GBRForest.h"
#include <TFile.h>
#include <TSystem.h>
#include <cmath>
#include <vector>
#include <iostream>
using namespace std;
/*****************************************************************/
EpCombinationTool::EpCombinationTool()
: m_forest(nullptr),
m_ownForest(false)
/*****************************************************************/
{}
/*****************************************************************/
EpCombinationTool::~EpCombinationTool()
/*****************************************************************/
{
if (m_ownForest)
delete m_forest;
}
/*****************************************************************/
bool EpCombinationTool::init(const std::string& regressionFileName, const std::string& bdtName)
/*****************************************************************/
{
TFile* regressionFile = TFile::Open(regressionFileName.c_str());
if (!regressionFile) {
cout << "ERROR: Cannot open regression file " << regressionFileName << "\n";
return false;
}
if (m_ownForest)
delete m_forest;
m_forest = (GBRForest*)regressionFile->Get(bdtName.c_str());
m_ownForest = true;
//regressionFile->GetObject(bdtName.c_str(), m_forest);
if (!m_forest) {
cout << "ERROR: Cannot find forest " << bdtName << " in " << regressionFileName << "\n";
regressionFile->Close();
return false;
}
regressionFile->Close();
return true;
}
bool EpCombinationTool::init(const GBRForest* forest) {
if (m_ownForest)
delete m_forest;
m_forest = forest;
m_ownForest = false;
return true;
}
/*****************************************************************/
void EpCombinationTool::combine(SimpleElectron& mySimpleElectron) const
/*****************************************************************/
{
if (!m_forest) {
cout << "ERROR: The combination tool is not initialized\n";
return;
}
float energy = mySimpleElectron.getNewEnergy();
float energyError = mySimpleElectron.getNewEnergyError();
float momentum = mySimpleElectron.getTrackerMomentum();
float momentumError = mySimpleElectron.getTrackerMomentumError();
int electronClass = mySimpleElectron.getElClass();
bool isEcalDriven = mySimpleElectron.isEcalDriven();
bool isTrackerDriven = mySimpleElectron.isTrackerDriven();
bool isEB = mySimpleElectron.isEB();
// compute relative errors and ratio of errors
float energyRelError = energyError / energy;
float momentumRelError = momentumError / momentum;
float errorRatio = energyRelError / momentumRelError;
// calculate E/p and corresponding error
float eOverP = energy / momentum;
float eOverPerror =
sqrt((energyError / momentum) * (energyError / momentum) +
(energy * momentumError / momentum / momentum) * (energy * momentumError / momentum / momentum));
// fill input variables
float regressionInputs[11];
regressionInputs[0] = energy;
regressionInputs[1] = energyRelError;
regressionInputs[2] = momentum;
regressionInputs[3] = momentumRelError;
regressionInputs[4] = errorRatio;
regressionInputs[5] = eOverP;
regressionInputs[6] = eOverPerror;
regressionInputs[7] = static_cast<float>(isEcalDriven);
regressionInputs[8] = static_cast<float>(isTrackerDriven);
regressionInputs[9] = static_cast<float>(electronClass);
regressionInputs[10] = static_cast<float>(isEB);
// retrieve combination weight
float weight = 0.;
if (eOverP > 0.025 &&
fabs(momentum - energy) < 15. * sqrt(momentumError * momentumError + energyError * energyError) &&
((momentumError < 10. * momentum) || (energy < 200.))) // protection against crazy track measurement
{
weight = m_forest->GetResponse(regressionInputs);
if (weight > 1.)
weight = 1.;
else if (weight < 0.)
weight = 0.;
}
float combinedMomentum = weight * momentum + (1. - weight) * energy;
float combinedMomentumError =
sqrt(weight * weight * momentumError * momentumError + (1. - weight) * (1. - weight) * energyError * energyError);
// FIXME : pure tracker electrons have track momentum error of 999.
// If the combination try to combine such electrons then the original combined momentum is kept
if (momentumError != 999. || weight == 0.) {
mySimpleElectron.setCombinedMomentum(combinedMomentum);
mySimpleElectron.setCombinedMomentumError(combinedMomentumError);
}
}
|