|
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
|
#include "EgammaAnalysis/ElectronTools/interface/ElectronEnergyCalibratorRun2.h"
#include <CLHEP/Random/RandGaussQ.h>
#include "FWCore/AbstractServices/interface/RandomNumberGenerator.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Utilities/interface/Exception.h"
ElectronEnergyCalibratorRun2::ElectronEnergyCalibratorRun2(EpCombinationTool &combinator,
bool isMC,
bool synchronization,
std::string correctionFile)
: epCombinationTool_(&combinator),
isMC_(isMC),
synchronization_(synchronization),
rng_(nullptr),
_correctionRetriever(correctionFile) // here is opening the files and reading the corrections
{
if (isMC_) {
_correctionRetriever.doScale = false;
_correctionRetriever.doSmearings = true;
} else {
_correctionRetriever.doScale = true;
_correctionRetriever.doSmearings = false;
}
}
ElectronEnergyCalibratorRun2::~ElectronEnergyCalibratorRun2() {}
void ElectronEnergyCalibratorRun2::initPrivateRng(TRandom *rnd) { rng_ = rnd; }
void ElectronEnergyCalibratorRun2::calibrate(reco::GsfElectron &electron,
unsigned int runNumber,
edm::StreamID const &id) const {
SimpleElectron simple(electron, runNumber, isMC_);
calibrate(simple, id);
simple.writeTo(electron);
}
void ElectronEnergyCalibratorRun2::calibrate(SimpleElectron &electron, edm::StreamID const &id) const {
assert(isMC_ == electron.isMC());
float smear = 0.0, scale = 1.0;
float aeta = std::abs(electron.getEta()); //, r9 = electron.getR9();
float et = electron.getNewEnergy() / cosh(aeta);
scale = _correctionRetriever.ScaleCorrection(electron.getRunNumber(), electron.isEB(), electron.getR9(), aeta, et);
smear = _correctionRetriever.getSmearingSigma(
electron.getRunNumber(), electron.isEB(), electron.getR9(), aeta, et, 0., 0.);
double newEcalEnergy, newEcalEnergyError;
if (isMC_) {
double corr = 1.0 + smear * gauss(id);
newEcalEnergy = electron.getNewEnergy() * corr;
newEcalEnergyError = std::hypot(electron.getNewEnergyError() * corr, smear * newEcalEnergy);
} else {
newEcalEnergy = electron.getNewEnergy() * scale;
newEcalEnergyError = std::hypot(electron.getNewEnergyError() * scale, smear * newEcalEnergy);
}
electron.setNewEnergy(newEcalEnergy);
electron.setNewEnergyError(newEcalEnergyError);
epCombinationTool_->combine(electron);
}
double ElectronEnergyCalibratorRun2::gauss(edm::StreamID const &id) const {
if (synchronization_)
return 1.0;
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();
}
}
|