CMS 3D CMS Logo

ElectronEnergyCalibrator.cc
Go to the documentation of this file.
2 
7 #include <CLHEP/Random/RandGaussQ.h>
8 
11 
13  const std::string& correctionFile) :
14  correctionRetriever_(correctionFile),
15  epCombinationTool_(&combinator),
16  rng_(nullptr),
17  minEt_(1.0)
18 {
19 
20 }
21 
23 {
24  rng_ = rnd;
25 }
26 
27 std::array<float,EGEnergySysIndex::kNrSysErrs> ElectronEnergyCalibrator::
29  const unsigned int runNumber,
30  const EcalRecHitCollection *recHits,
31  edm::StreamID const & id,
32  const ElectronEnergyCalibrator::EventType eventType) const
33 {
34  return calibrate(ele,runNumber,recHits,gauss(id),eventType);
35 }
36 
37 std::array<float,EGEnergySysIndex::kNrSysErrs> ElectronEnergyCalibrator::
39  const EcalRecHitCollection *recHits,
40  const float smearNrSigma,
41  const ElectronEnergyCalibrator::EventType eventType) const
42 {
43  const float scEtaAbs = std::abs(ele.superCluster()->eta());
44  const float et = ele.ecalEnergy() / cosh(scEtaAbs);
45 
46  if (et < minEt_ || edm::isNotFinite(et) ) {
47  std::array<float,EGEnergySysIndex::kNrSysErrs> retVal;
48  retVal.fill(ele.energy());
49  retVal[EGEnergySysIndex::kScaleValue] = 1.0;
50  retVal[EGEnergySysIndex::kSmearValue] = 0.0;
51  retVal[EGEnergySysIndex::kSmearNrSigma] = smearNrSigma;
60  return retVal;
61  }
62 
63  const DetId seedDetId = ele.superCluster()->seed()->seed();
64  EcalRecHitCollection::const_iterator seedRecHit = recHits->find(seedDetId);
65  unsigned int gainSeedSC = 12;
66  if (seedRecHit != recHits->end()) {
67  if(seedRecHit->checkFlag(EcalRecHit::kHasSwitchToGain6)) gainSeedSC = 6;
68  if(seedRecHit->checkFlag(EcalRecHit::kHasSwitchToGain1)) gainSeedSC = 1;
69  }
70 
71  const EnergyScaleCorrection::ScaleCorrection* scaleCorr = correctionRetriever_.getScaleCorr(runNumber, et, scEtaAbs, ele.full5x5_r9(), gainSeedSC);
72  const EnergyScaleCorrection::SmearCorrection* smearCorr = correctionRetriever_.getSmearCorr(runNumber, et, scEtaAbs, ele.full5x5_r9(), gainSeedSC);
73  if(scaleCorr==nullptr) scaleCorr=&defaultScaleCorr_;
74  if(smearCorr==nullptr) smearCorr=&defaultSmearCorr_;
75 
76  std::array<float,EGEnergySysIndex::kNrSysErrs> uncertainties{};
77 
78  uncertainties[EGEnergySysIndex::kScaleValue] = scaleCorr->scale();
79  uncertainties[EGEnergySysIndex::kSmearValue] = smearCorr->sigma(et); //even though we use scale = 1.0, we still store the value returned for MC
80  uncertainties[EGEnergySysIndex::kSmearNrSigma] = smearNrSigma;
81  //MC central values are not scaled (scale = 1.0), data is not smeared (smearNrSigma = 0)
82  //the smearing (or resolution extra parameter as it might better be called)
83  //still has a second order effect on data as it enters the E/p combination as an adjustment
84  //to the estimate of the resolution contained in caloEnergyError
85  //MC gets all the scale systematics
86  if(eventType == EventType::DATA){
87  setEnergyAndSystVarations(scaleCorr->scale(),0.,et,*scaleCorr,*smearCorr,ele,uncertainties);
88  }else if(eventType == EventType::MC){
89  setEnergyAndSystVarations(1.0,smearNrSigma,et,*scaleCorr,*smearCorr,ele,uncertainties);
90  }
91 
92  return uncertainties;
93 
94 }
95 
97 setEnergyAndSystVarations(const float scale,const float smearNrSigma,const float et,
100  reco::GsfElectron& ele,
101  std::array<float,EGEnergySysIndex::kNrSysErrs>& energyData)const
102 {
103 
104  const float smear = smearCorr.sigma(et);
105  const float smearRhoUp = smearCorr.sigma(et,1,0);
106  const float smearRhoDn = smearCorr.sigma(et,-1,0);
107  const float smearPhiUp = smearCorr.sigma(et,0,1);
108  const float smearPhiDn = smearCorr.sigma(et,0,-1);
109  const float smearUp = smearRhoUp;
110  const float smearDn = smearRhoDn;
111 
112  const float corr = scale + smear * smearNrSigma;
113  const float corrRhoUp = scale + smearRhoUp * smearNrSigma;
114  const float corrRhoDn = scale + smearRhoDn * smearNrSigma;
115  const float corrPhiUp = scale + smearPhiUp * smearNrSigma;
116  const float corrPhiDn = scale + smearPhiDn * smearNrSigma;
117  const float corrUp = corrRhoUp;
118  const float corrDn = corrRhoDn;
119 
120  const float corrScaleStatUp = corr+scaleCorr.scaleErrStat();
121  const float corrScaleStatDn = corr-scaleCorr.scaleErrStat();
122  const float corrScaleSystUp = corr+scaleCorr.scaleErrSyst();
123  const float corrScaleSystDn = corr-scaleCorr.scaleErrSyst();
124  const float corrScaleGainUp = corr+scaleCorr.scaleErrGain();
125  const float corrScaleGainDn = corr-scaleCorr.scaleErrGain();
126  const float corrScaleUp = corr+scaleCorr.scaleErr(EnergyScaleCorrection::kErrStatSystGain);
127  const float corrScaleDn = corr-scaleCorr.scaleErr(EnergyScaleCorrection::kErrStatSystGain);
128 
129  const math::XYZTLorentzVector oldP4 = ele.p4();
130  energyData[EGEnergySysIndex::kEcalTrkPreCorr] = ele.energy();
132  energyData[EGEnergySysIndex::kEcalPreCorr] = ele.ecalEnergy();
134 
135  energyData[EGEnergySysIndex::kScaleStatUp] = calCombinedMom(ele,corrScaleStatUp,smear).first;
136  energyData[EGEnergySysIndex::kScaleStatDown] = calCombinedMom(ele,corrScaleStatDn,smear).first;
137  energyData[EGEnergySysIndex::kScaleSystUp] = calCombinedMom(ele,corrScaleSystUp,smear).first;
138  energyData[EGEnergySysIndex::kScaleSystDown] = calCombinedMom(ele,corrScaleSystDn,smear).first;
139  energyData[EGEnergySysIndex::kScaleGainUp] = calCombinedMom(ele,corrScaleGainUp,smear).first;
140  energyData[EGEnergySysIndex::kScaleGainDown] = calCombinedMom(ele,corrScaleGainDn,smear).first;
141 
142  energyData[EGEnergySysIndex::kSmearRhoUp] = calCombinedMom(ele,corrRhoUp,smearRhoUp).first;
143  energyData[EGEnergySysIndex::kSmearRhoDown] = calCombinedMom(ele,corrRhoDn,smearRhoDn).first;
144  energyData[EGEnergySysIndex::kSmearPhiUp] = calCombinedMom(ele,corrPhiUp,smearPhiUp).first;
145  energyData[EGEnergySysIndex::kSmearPhiDown] = calCombinedMom(ele,corrPhiDn,smearPhiDn).first;
146 
147  energyData[EGEnergySysIndex::kScaleUp] = calCombinedMom(ele,corrScaleUp,smear).first;
148  energyData[EGEnergySysIndex::kScaleDown] = calCombinedMom(ele,corrScaleDn,smear).first;
149  energyData[EGEnergySysIndex::kSmearUp] = calCombinedMom(ele,corrUp,smearUp).first;
150  energyData[EGEnergySysIndex::kSmearDown] = calCombinedMom(ele,corrDn,smearDn).first;
151 
152  const std::pair<float, float> combinedMomentum = calCombinedMom(ele,corr,smear);
153  setEcalEnergy(ele,corr,smear);
154  const float energyCorr = combinedMomentum.first / oldP4.t();
155 
156  const math::XYZTLorentzVector newP4(oldP4.x() * energyCorr,
157  oldP4.y() * energyCorr,
158  oldP4.z() * energyCorr,
159  combinedMomentum.first);
160 
161  ele.correctMomentum(newP4, ele.trackMomentumError(), combinedMomentum.second);
162  energyData[EGEnergySysIndex::kEcalTrkPostCorr] = ele.energy();
164 
165  energyData[EGEnergySysIndex::kEcalPostCorr] = ele.ecalEnergy();
167 
168 }
169 
170 
172  const float scale,
173  const float smear)const
174 {
175  const float oldEcalEnergy = ele.ecalEnergy();
176  const float oldEcalEnergyErr = ele.ecalEnergyError();
177  ele.setCorrectedEcalEnergy( oldEcalEnergy*scale );
178  ele.setCorrectedEcalEnergyError(std::hypot( oldEcalEnergyErr*scale, oldEcalEnergy*smear*scale ) );
179 }
180 
182  const float scale,
183  const float smear)const
184 {
185  const float oldEcalEnergy = ele.ecalEnergy();
186  const float oldEcalEnergyErr = ele.ecalEnergyError();
187 
188  const auto oldP4 = ele.p4();
189  const float oldP4Err = ele.p4Error(reco::GsfElectron::P4_COMBINATION);
190  const float oldTrkMomErr = ele.trackMomentumError();
191 
192  setEcalEnergy(ele,scale,smear);
193  const auto& combinedMomentum = epCombinationTool_->combine(ele,oldEcalEnergyErr*scale);
194 
195  ele.setCorrectedEcalEnergy(oldEcalEnergy);
196  ele.setCorrectedEcalEnergyError(oldEcalEnergyErr);
197  ele.correctMomentum(oldP4,oldTrkMomErr,oldP4Err);
198 
199  return combinedMomentum;
200 }
201 
202 
204 {
205  if (rng_) {
206  return rng_->Gaus();
207  } else {
209  if ( !rng.isAvailable() ) {
210  throw cms::Exception("Configuration")
211  << "XXXXXXX requires the RandomNumberGeneratorService\n"
212  "which is not present in the configuration file. You must add the service\n"
213  "in the configuration file or remove the modules that require it.";
214  }
215  CLHEP::RandGaussQ gaussDistribution(rng->getEngine(id), 0.0, 1.0);
216  return gaussDistribution.fire();
217  }
218 }
float trackMomentumError() const
Definition: GsfElectron.h:848
const Corrections & corrections() const
Definition: GsfElectron.h:853
const LorentzVector & p4(P4Kind kind) const
Definition: GsfElectron.cc:228
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
void setEcalEnergy(reco::GsfElectron &ele, const float scale, const float smear) const
EnergyScaleCorrection correctionRetriever_
void correctMomentum(const LorentzVector &p4, float trackMomentumError, float p4Error)
Definition: GsfElectron.h:870
#define nullptr
std::vector< EcalRecHit >::const_iterator const_iterator
void combine(SimpleElectron &mySimpleElectron) const
float p4Error(P4Kind kind) const
Definition: GsfElectron.cc:240
void calibrate(SimpleElectron &electron, edm::StreamID const &)
virtual CLHEP::HepRandomEngine & getEngine(StreamID const &)=0
Use this engine in event methods.
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
void setCorrectedEcalEnergyError(float newEnergyError)
Definition: GsfElectron.cc:179
const ScaleCorrection * getScaleCorr(unsigned int runnr, double et, double eta, double r9, unsigned int gainSeed) const
double energy() const final
energy
bool isAvailable() const
Definition: Service.h:40
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
const SmearCorrection * getSmearCorr(unsigned int runnr, double et, double eta, double r9, unsigned int gainSeed) const
float ecalEnergyError() const
Definition: GsfElectron.h:861
JetCorrectorParameters corr
Definition: classes.h:5
float scaleErr(const std::bitset< kErrNrBits > &uncBitMask) const
float sigma(const float et, const float nrSigmaRho=0., const float nrSigmaPhi=0.) const
const_iterator end() const
double gauss(edm::StreamID const &id) const
Definition: DetId.h:18
et
define resolution functions of each parameter
float ecalEnergy() const
Definition: GsfElectron.h:860
static const EnergyScaleCorrection::SmearCorrection defaultSmearCorr_
const EpCombinationTool * epCombinationTool_
float full5x5_r9() const
Definition: GsfElectron.h:467
void 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
void setCorrectedEcalEnergy(float newEnergy)
Definition: GsfElectron.cc:182
static const EnergyScaleCorrection::ScaleCorrection defaultScaleCorr_
iterator find(key_type k)
SuperClusterRef superCluster() const override
reference to a SuperCluster
Definition: GsfElectron.h:185
std::pair< float, float > calCombinedMom(reco::GsfElectron &ele, const float scale, const float smear) const