-
Notifications
You must be signed in to change notification settings - Fork 4.2k
/
ElectronEnergyCalibrator.cc
220 lines (184 loc) · 9.91 KB
/
ElectronEnergyCalibrator.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
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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
#include "RecoEgamma/EgammaTools/interface/ElectronEnergyCalibrator.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Utilities/interface/RandomNumberGenerator.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/Utilities/interface/isFinite.h"
#include <CLHEP/Random/RandGaussQ.h>
const EnergyScaleCorrection::ScaleCorrection ElectronEnergyCalibrator::defaultScaleCorr_;
const EnergyScaleCorrection::SmearCorrection ElectronEnergyCalibrator::defaultSmearCorr_;
ElectronEnergyCalibrator::ElectronEnergyCalibrator(const EpCombinationTool &combinator,
const std::string& correctionFile) :
correctionRetriever_(correctionFile),
epCombinationTool_(&combinator),
rng_(nullptr),
minEt_(1.0),
useSmearCorrEcalEnergyErrInComb_(false)
{
}
void ElectronEnergyCalibrator::initPrivateRng(TRandom *rnd)
{
rng_ = rnd;
}
std::array<float,EGEnergySysIndex::kNrSysErrs> ElectronEnergyCalibrator::
calibrate(reco::GsfElectron &ele,
const unsigned int runNumber,
const EcalRecHitCollection *recHits,
edm::StreamID const & id,
const ElectronEnergyCalibrator::EventType eventType) const
{
return calibrate(ele,runNumber,recHits,gauss(id),eventType);
}
std::array<float,EGEnergySysIndex::kNrSysErrs> ElectronEnergyCalibrator::
calibrate(reco::GsfElectron &ele, unsigned int runNumber,
const EcalRecHitCollection *recHits,
const float smearNrSigma,
const ElectronEnergyCalibrator::EventType eventType) const
{
const float scEtaAbs = std::abs(ele.superCluster()->eta());
const float et = ele.ecalEnergy() / cosh(scEtaAbs);
if (et < minEt_ || edm::isNotFinite(et) ) {
std::array<float,EGEnergySysIndex::kNrSysErrs> retVal;
retVal.fill(ele.energy());
retVal[EGEnergySysIndex::kScaleValue] = 1.0;
retVal[EGEnergySysIndex::kSmearValue] = 0.0;
retVal[EGEnergySysIndex::kSmearNrSigma] = smearNrSigma;
retVal[EGEnergySysIndex::kEcalPreCorr] = ele.ecalEnergy();
retVal[EGEnergySysIndex::kEcalErrPreCorr] = ele.ecalEnergyError();
retVal[EGEnergySysIndex::kEcalPostCorr] = ele.ecalEnergy();
retVal[EGEnergySysIndex::kEcalErrPostCorr] = ele.ecalEnergyError();
retVal[EGEnergySysIndex::kEcalTrkPreCorr] = ele.energy();
retVal[EGEnergySysIndex::kEcalTrkErrPreCorr] = ele.corrections().combinedP4Error;
retVal[EGEnergySysIndex::kEcalTrkPostCorr] = ele.energy();
retVal[EGEnergySysIndex::kEcalTrkErrPostCorr] = ele.corrections().combinedP4Error;
return retVal;
}
const DetId seedDetId = ele.superCluster()->seed()->seed();
EcalRecHitCollection::const_iterator seedRecHit = recHits->find(seedDetId);
unsigned int gainSeedSC = 12;
if (seedRecHit != recHits->end()) {
if(seedRecHit->checkFlag(EcalRecHit::kHasSwitchToGain6)) gainSeedSC = 6;
if(seedRecHit->checkFlag(EcalRecHit::kHasSwitchToGain1)) gainSeedSC = 1;
}
const EnergyScaleCorrection::ScaleCorrection* scaleCorr = correctionRetriever_.getScaleCorr(runNumber, et, scEtaAbs, ele.full5x5_r9(), gainSeedSC);
const EnergyScaleCorrection::SmearCorrection* smearCorr = correctionRetriever_.getSmearCorr(runNumber, et, scEtaAbs, ele.full5x5_r9(), gainSeedSC);
if(scaleCorr==nullptr) scaleCorr=&defaultScaleCorr_;
if(smearCorr==nullptr) smearCorr=&defaultSmearCorr_;
std::array<float,EGEnergySysIndex::kNrSysErrs> uncertainties{};
uncertainties[EGEnergySysIndex::kScaleValue] = scaleCorr->scale();
uncertainties[EGEnergySysIndex::kSmearValue] = smearCorr->sigma(et); //even though we use scale = 1.0, we still store the value returned for MC
uncertainties[EGEnergySysIndex::kSmearNrSigma] = smearNrSigma;
//MC central values are not scaled (scale = 1.0), data is not smeared (smearNrSigma = 0)
//the smearing (or resolution extra parameter as it might better be called)
//still has a second order effect on data as it enters the E/p combination as an adjustment
//to the estimate of the resolution contained in caloEnergyError
//MC gets all the scale systematics
if(eventType == EventType::DATA){
setEnergyAndSystVarations(scaleCorr->scale(),0.,et,*scaleCorr,*smearCorr,ele,uncertainties);
}else if(eventType == EventType::MC){
setEnergyAndSystVarations(1.0,smearNrSigma,et,*scaleCorr,*smearCorr,ele,uncertainties);
}
return uncertainties;
}
void ElectronEnergyCalibrator::
setEnergyAndSystVarations(const float scale,const float smearNrSigma,const float et,
const EnergyScaleCorrection::ScaleCorrection& scaleCorr,
const EnergyScaleCorrection::SmearCorrection& smearCorr,
reco::GsfElectron& ele,
std::array<float,EGEnergySysIndex::kNrSysErrs>& energyData)const
{
const float smear = smearCorr.sigma(et);
const float smearRhoUp = smearCorr.sigma(et,1,0);
const float smearRhoDn = smearCorr.sigma(et,-1,0);
const float smearPhiUp = smearCorr.sigma(et,0,1);
const float smearPhiDn = smearCorr.sigma(et,0,-1);
const float smearUp = smearRhoUp;
const float smearDn = smearRhoDn;
const float corr = scale + smear * smearNrSigma;
const float corrRhoUp = scale + smearRhoUp * smearNrSigma;
const float corrRhoDn = scale + smearRhoDn * smearNrSigma;
const float corrPhiUp = scale + smearPhiUp * smearNrSigma;
const float corrPhiDn = scale + smearPhiDn * smearNrSigma;
const float corrUp = corrRhoUp;
const float corrDn = corrRhoDn;
const float corrScaleStatUp = corr+scaleCorr.scaleErrStat();
const float corrScaleStatDn = corr-scaleCorr.scaleErrStat();
const float corrScaleSystUp = corr+scaleCorr.scaleErrSyst();
const float corrScaleSystDn = corr-scaleCorr.scaleErrSyst();
const float corrScaleGainUp = corr+scaleCorr.scaleErrGain();
const float corrScaleGainDn = corr-scaleCorr.scaleErrGain();
const float corrScaleUp = corr+scaleCorr.scaleErr(EnergyScaleCorrection::kErrStatSystGain);
const float corrScaleDn = corr-scaleCorr.scaleErr(EnergyScaleCorrection::kErrStatSystGain);
const math::XYZTLorentzVector oldP4 = ele.p4();
energyData[EGEnergySysIndex::kEcalTrkPreCorr] = ele.energy();
energyData[EGEnergySysIndex::kEcalTrkErrPreCorr] = ele.corrections().combinedP4Error;
energyData[EGEnergySysIndex::kEcalPreCorr] = ele.ecalEnergy();
energyData[EGEnergySysIndex::kEcalErrPreCorr] = ele.ecalEnergyError();
energyData[EGEnergySysIndex::kScaleStatUp] = calCombinedMom(ele,corrScaleStatUp,smear).first;
energyData[EGEnergySysIndex::kScaleStatDown] = calCombinedMom(ele,corrScaleStatDn,smear).first;
energyData[EGEnergySysIndex::kScaleSystUp] = calCombinedMom(ele,corrScaleSystUp,smear).first;
energyData[EGEnergySysIndex::kScaleSystDown] = calCombinedMom(ele,corrScaleSystDn,smear).first;
energyData[EGEnergySysIndex::kScaleGainUp] = calCombinedMom(ele,corrScaleGainUp,smear).first;
energyData[EGEnergySysIndex::kScaleGainDown] = calCombinedMom(ele,corrScaleGainDn,smear).first;
energyData[EGEnergySysIndex::kSmearRhoUp] = calCombinedMom(ele,corrRhoUp,smearRhoUp).first;
energyData[EGEnergySysIndex::kSmearRhoDown] = calCombinedMom(ele,corrRhoDn,smearRhoDn).first;
energyData[EGEnergySysIndex::kSmearPhiUp] = calCombinedMom(ele,corrPhiUp,smearPhiUp).first;
energyData[EGEnergySysIndex::kSmearPhiDown] = calCombinedMom(ele,corrPhiDn,smearPhiDn).first;
energyData[EGEnergySysIndex::kScaleUp] = calCombinedMom(ele,corrScaleUp,smear).first;
energyData[EGEnergySysIndex::kScaleDown] = calCombinedMom(ele,corrScaleDn,smear).first;
energyData[EGEnergySysIndex::kSmearUp] = calCombinedMom(ele,corrUp,smearUp).first;
energyData[EGEnergySysIndex::kSmearDown] = calCombinedMom(ele,corrDn,smearDn).first;
const std::pair<float, float> combinedMomentum = calCombinedMom(ele,corr,smear);
setEcalEnergy(ele,corr,smear);
const float energyCorr = combinedMomentum.first / oldP4.t();
const math::XYZTLorentzVector newP4(oldP4.x() * energyCorr,
oldP4.y() * energyCorr,
oldP4.z() * energyCorr,
combinedMomentum.first);
ele.correctMomentum(newP4, ele.trackMomentumError(), combinedMomentum.second);
energyData[EGEnergySysIndex::kEcalTrkPostCorr] = ele.energy();
energyData[EGEnergySysIndex::kEcalTrkErrPostCorr] = ele.corrections().combinedP4Error;
energyData[EGEnergySysIndex::kEcalPostCorr] = ele.ecalEnergy();
energyData[EGEnergySysIndex::kEcalErrPostCorr] = ele.ecalEnergyError();
}
void ElectronEnergyCalibrator::setEcalEnergy(reco::GsfElectron& ele,
const float scale,
const float smear)const
{
const float oldEcalEnergy = ele.ecalEnergy();
const float oldEcalEnergyErr = ele.ecalEnergyError();
ele.setCorrectedEcalEnergy( oldEcalEnergy*scale );
ele.setCorrectedEcalEnergyError(std::hypot( oldEcalEnergyErr*scale, oldEcalEnergy*smear*scale ) );
}
std::pair<float,float> ElectronEnergyCalibrator::calCombinedMom(reco::GsfElectron& ele,
const float scale,
const float smear)const
{
const float oldEcalEnergy = ele.ecalEnergy();
const float oldEcalEnergyErr = ele.ecalEnergyError();
const auto oldP4 = ele.p4();
const float oldP4Err = ele.p4Error(reco::GsfElectron::P4_COMBINATION);
const float oldTrkMomErr = ele.trackMomentumError();
setEcalEnergy(ele,scale,smear);
float ecalEnergyErrForComb = useSmearCorrEcalEnergyErrInComb_ ? ele.correctedEcalEnergyError() : oldEcalEnergyErr*scale;
const auto& combinedMomentum = epCombinationTool_->combine(ele,ecalEnergyErrForComb);
ele.setCorrectedEcalEnergy(oldEcalEnergy);
ele.setCorrectedEcalEnergyError(oldEcalEnergyErr);
ele.correctMomentum(oldP4,oldTrkMomErr,oldP4Err);
return combinedMomentum;
}
double ElectronEnergyCalibrator::gauss(edm::StreamID const& id) const
{
if (rng_) {
return rng_->Gaus();
} else {
edm::Service<edm::RandomNumberGenerator> rng;
if ( !rng.isAvailable() ) {
throw cms::Exception("Configuration")
<< "XXXXXXX requires the RandomNumberGeneratorService\n"
"which is not present in the configuration file. You must add the service\n"
"in the configuration file or remove the modules that require it.";
}
CLHEP::RandGaussQ gaussDistribution(rng->getEngine(id), 0.0, 1.0);
return gaussDistribution.fire();
}
}