CMS 3D CMS Logo

RegressionHelper.cc
Go to the documentation of this file.
5 #include "TVector2.h"
6 #include "TFile.h"
7 
9  bool useEcalReg,
10  bool useCombinationReg,
12  : caloTopology{cc.esConsumes()}, caloGeometry{cc.esConsumes()} {
13  if (useEcalReg && cfg.ecalWeightsFromDB) {
14  ecalRegBarrel = cc.esConsumes<GBRForest, GBRWrapperRcd>(edm::ESInputTag("", cfg.ecalRegressionWeightLabels[0]));
15  ecalRegEndcap = cc.esConsumes<GBRForest, GBRWrapperRcd>(edm::ESInputTag("", cfg.ecalRegressionWeightLabels[1]));
16  ecalRegErrorBarrel =
17  cc.esConsumes<GBRForest, GBRWrapperRcd>(edm::ESInputTag("", cfg.ecalRegressionWeightLabels[2]));
18  ecalRegErrorEndcap =
19  cc.esConsumes<GBRForest, GBRWrapperRcd>(edm::ESInputTag("", cfg.ecalRegressionWeightLabels[3]));
20  }
21  if (useCombinationReg && cfg.combinationWeightsFromDB) {
22  combinationReg =
23  cc.esConsumes<GBRForest, GBRWrapperRcd>(edm::ESInputTag("", cfg.combinationRegressionWeightLabels[0]));
24  }
25 }
26 
28  bool useEcalReg,
29  bool useCombinationReg,
31  : cfg_(config), esGetTokens_{cfg_, useEcalReg, useCombinationReg, cc} {}
32 
35  const EcalRecHitCollection& rechitsEB,
36  const EcalRecHitCollection& rechitsEE) const {
37  double cor, err;
38  getEcalRegression(*ele.superCluster(), vertices, rechitsEB, rechitsEE, cor, err);
39  ele.setCorrectedEcalEnergy(cor * ele.superCluster()->correctedEnergy());
40  ele.setCorrectedEcalEnergyError(err * ele.superCluster()->correctedEnergy());
41 }
42 
46 
47  // Ecal regression
48 
49  // if at least one of the set of weights come from the DB
50  if (cfg_.ecalWeightsFromDB) {
51  ecalRegBarrel_ = &es.getData(esGetTokens_.ecalRegBarrel); // ECAL barrel
52  ecalRegEndcap_ = &es.getData(esGetTokens_.ecalRegEndcap); // ECAL endcaps
53  ecalRegErrorBarrel_ = &es.getData(esGetTokens_.ecalRegErrorBarrel); // ECAL barrel error
54  ecalRegErrorEndcap_ = &es.getData(esGetTokens_.ecalRegErrorEndcap); // ECAL endcap error
56  }
58  // Combination
61  }
62 
63  // read weights from file - for debugging. Even if it is one single files, 4 files should b set in the vector
65  TFile file0(edm::FileInPath(cfg_.ecalRegressionWeightFiles[0].c_str()).fullPath().c_str());
66  ecalRegBarrel_ = (const GBRForest*)file0.Get(cfg_.ecalRegressionWeightLabels[0].c_str());
67  file0.Close();
68  TFile file1(edm::FileInPath(cfg_.ecalRegressionWeightFiles[1].c_str()).fullPath().c_str());
70  file1.Close();
71  TFile file2(edm::FileInPath(cfg_.ecalRegressionWeightFiles[2].c_str()).fullPath().c_str());
73  file2.Close();
74  TFile file3(edm::FileInPath(cfg_.ecalRegressionWeightFiles[3].c_str()).fullPath().c_str());
75  ecalRegErrorEndcap_ = (const GBRForest*)file3.Get(cfg_.ecalRegressionWeightLabels[3].c_str());
77  file3.Close();
78  }
79 
82  TFile file0(edm::FileInPath(cfg_.combinationRegressionWeightFiles[0].c_str()).fullPath().c_str());
83  combinationReg_ = (const GBRForest*)file0.Get(cfg_.combinationRegressionWeightLabels[0].c_str());
85  file0.Close();
86  }
87 }
88 
91  const EcalRecHitCollection& rechitsEB,
92  const EcalRecHitCollection& rechitsEE,
93  double& energyFactor,
94  double& errorFactor) const {
95  energyFactor = -999.;
96  errorFactor = -999.;
97 
98  std::vector<float> rInputs;
99  EcalRegressionData regData;
100  regData.fill(sc, &rechitsEB, &rechitsEE, caloGeometry_, caloTopology_, &vertices);
101  regData.fillVec(rInputs);
102  if (sc.seed()->hitsAndFractions()[0].first.subdetId() == EcalBarrel) {
103  energyFactor = ecalRegBarrel_->GetResponse(&rInputs[0]);
104  errorFactor = ecalRegErrorBarrel_->GetResponse(&rInputs[0]);
105  } else if (sc.seed()->hitsAndFractions()[0].first.subdetId() == EcalEndcap) {
106  energyFactor = ecalRegEndcap_->GetResponse(&rInputs[0]);
107  errorFactor = ecalRegErrorEndcap_->GetResponse(&rInputs[0]);
108  } else {
109  throw cms::Exception("RegressionHelper::calculateRegressedEnergy")
110  << "Supercluster seed is either EB nor EE!" << std::endl;
111  }
112 }
113 
115  float energy = ele.correctedEcalEnergy();
116  float energyError = ele.correctedEcalEnergyError();
117  float momentum = ele.trackMomentumAtVtx().R();
118  float momentumError = ele.trackMomentumError();
119  int elClass = -1;
120 
121  switch (ele.classification()) {
123  elClass = 0;
124  break;
126  elClass = 1;
127  break;
129  elClass = 2;
130  break;
132  elClass = 3;
133  break;
135  elClass = 4;
136  break;
137  default:
138  elClass = -1;
139  }
140 
141  bool isEcalDriven = ele.ecalDriven();
142  bool isTrackerDriven = ele.trackerDrivenSeed();
143  bool isEB = ele.isEB();
144 
145  // compute relative errors and ratio of errors
146  float energyRelError = energyError / energy;
147  float momentumRelError = momentumError / momentum;
148  float errorRatio = energyRelError / momentumRelError;
149 
150  // calculate E/p and corresponding error
151  float eOverP = energy / momentum;
152  float eOverPerror = eOverP * std::hypot(energyRelError, momentumRelError);
153 
154  // fill input variables
155  std::vector<float> regressionInputs;
156  regressionInputs.resize(11, 0.);
157 
158  regressionInputs[0] = energy;
159  regressionInputs[1] = energyRelError;
160  regressionInputs[2] = momentum;
161  regressionInputs[3] = momentumRelError;
162  regressionInputs[4] = errorRatio;
163  regressionInputs[5] = eOverP;
164  regressionInputs[6] = eOverPerror;
165  regressionInputs[7] = static_cast<float>(isEcalDriven);
166  regressionInputs[8] = static_cast<float>(isTrackerDriven);
167  regressionInputs[9] = static_cast<float>(elClass);
168  regressionInputs[10] = static_cast<float>(isEB);
169 
170  // retrieve combination weight
171  float weight = 0.;
172  if (eOverP > 0.025 &&
173  fabs(momentum - energy) < 15. * sqrt(momentumError * momentumError +
174  energyError * energyError)) // protection against crazy track measurement
175  {
176  weight = combinationReg_->GetResponse(regressionInputs.data());
177  if (weight > 1.)
178  weight = 1.;
179  else if (weight < 0.)
180  weight = 0.;
181  }
182 
183  float combinedMomentum = weight * momentum + (1. - weight) * energy;
184  float combinedMomentumError =
185  sqrt(weight * weight * momentumError * momentumError + (1. - weight) * (1. - weight) * energyError * energyError);
186 
187  // FIXME : pure tracker electrons have track momentum error of 999.
188  // If the combination try to combine such electrons then the original combined momentum is kept
189  if (momentumError != 999. || weight == 0.) {
190  math::XYZTLorentzVector oldMomentum = ele.p4();
191  math::XYZTLorentzVector newMomentum(oldMomentum.x() * combinedMomentum / oldMomentum.t(),
192  oldMomentum.y() * combinedMomentum / oldMomentum.t(),
193  oldMomentum.z() * combinedMomentum / oldMomentum.t(),
194  combinedMomentum);
195 
196  ele.setP4(reco::GsfElectron::P4_COMBINATION, newMomentum, combinedMomentumError, true);
197  ele.setMass(0.0);
198  }
199 }
std::vector< std::string > ecalRegressionWeightFiles
void applyCombinationRegression(reco::GsfElectron &ele) const
void setP4(P4Kind kind, const LorentzVector &p4, float p4Error, bool setCandidate)
Definition: GsfElectron.cc:194
std::vector< std::string > combinationRegressionWeightLabels
edm::ESGetToken< GBRForest, GBRWrapperRcd > ecalRegErrorBarrel
RegressionHelper(Configuration const &, bool useEcalReg, bool useCombinationReg, edm::ConsumesCollector &cc)
MPlex< T, D1, D2, N > hypot(const MPlex< T, D1, D2, N > &a, const MPlex< T, D1, D2, N > &b)
Definition: Matriplex.h:616
const GBRForest * ecalRegEndcap_
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
edm::ESGetToken< CaloGeometry, CaloGeometryRecord > caloGeometry
void setMass(double m) final
set particle mass
edm::ESGetToken< GBRForest, GBRWrapperRcd > ecalRegEndcap
std::vector< std::string > ecalRegressionWeightLabels
bool trackerDrivenSeed() const
Definition: GsfElectron.h:159
void fill(const reco::SuperCluster &superClus, const EcalRecHitCollection *ebRecHits, const EcalRecHitCollection *eeRecHits, const CaloGeometry *geom, const CaloTopology *topology, const reco::VertexCollection *vertices)
uint32_t cc[maxCellsPerHit]
Definition: gpuFishbone.h:49
float trackMomentumError() const
Definition: GsfElectron.h:884
std::vector< Vertex > VertexCollection
collection of Vertex objects
Definition: VertexFwd.h:9
Definition: weight.py:1
edm::ESGetToken< CaloTopology, CaloTopologyRecord > caloTopology
Definition: config.py:1
const Configuration cfg_
bool ecalDriven() const
Definition: GsfElectron.cc:168
Classification classification() const
Definition: GsfElectron.h:805
bool isEB() const
Definition: GsfElectron.h:328
const CaloTopology * caloTopology_
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
const GBRForest * ecalRegBarrel_
const GBRForest * ecalRegErrorBarrel_
void setCorrectedEcalEnergyError(float newEnergyError)
Definition: GsfElectron.cc:170
const ESGetTokens esGetTokens_
float correctedEcalEnergy() const
Definition: GsfElectron.h:882
edm::ESGetToken< GBRForest, GBRWrapperRcd > ecalRegBarrel
edm::ESGetToken< GBRForest, GBRWrapperRcd > combinationReg
T sqrt(T t)
Definition: SSEVec.h:23
ESGetTokens(Configuration const &cfg, bool useEcalReg, bool useCombinationReg, edm::ConsumesCollector &cc)
math::XYZVectorF trackMomentumAtVtx() const
Definition: GsfElectron.h:268
const LorentzVector & p4(P4Kind kind) const
Definition: GsfElectron.cc:217
std::vector< std::string > combinationRegressionWeightFiles
void fillVec(std::vector< float > &inputVec) const
void applyEcalRegression(reco::GsfElectron &electron, const reco::VertexCollection &vertices, const EcalRecHitCollection &rechitsEB, const EcalRecHitCollection &rechitsEE) const
const GBRForest * ecalRegErrorEndcap_
edm::ESGetToken< GBRForest, GBRWrapperRcd > ecalRegErrorEndcap
double GetResponse(const float *vector) const
Definition: GBRForest.h:48
void setCorrectedEcalEnergy(float newEnergy)
Definition: GsfElectron.cc:174
const CaloClusterPtr & seed() const
seed BasicCluster
Definition: SuperCluster.h:79
float correctedEcalEnergyError() const
Definition: GsfElectron.h:883
bool combinationRegressionInitialized_
void getEcalRegression(const reco::SuperCluster &sc, const reco::VertexCollection &vertices, const EcalRecHitCollection &rechitsEB, const EcalRecHitCollection &rechitsEE, double &energyFactor, double &errorFactor) const
void checkSetup(const edm::EventSetup &)
SuperClusterRef superCluster() const override
reference to a SuperCluster
Definition: GsfElectron.h:155
const CaloGeometry * caloGeometry_
const GBRForest * combinationReg_