File indexing completed on 2023-03-17 11:17:54
0001 #include "RecoEgamma/EgammaTools/interface/EcalRegressionData.h"
0002
0003 #include "RecoEcal/EgammaCoreTools/interface/EcalTools.h"
0004 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"
0005 #include "DataFormats/Math/interface/deltaR.h"
0006 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0007 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0008 #include "RecoEgamma/EgammaTools/interface/EcalClusterLocal.h"
0009
0010 #include "DataFormats/VertexReco/interface/Vertex.h"
0011 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0012
0013 #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h"
0014 #include "FWCore/Utilities/interface/isFinite.h"
0015
0016 float EcalRegressionData::seedLeftRightAsym() const {
0017 float eLeftRightSum = eLeft() + eRight();
0018 float eLeftRightDiff = eLeft() - eRight();
0019 return eLeftRightSum != 0 ? eLeftRightDiff / eLeftRightSum : 0.;
0020 }
0021
0022 float EcalRegressionData::seedTopBottomAsym() const {
0023 float eTopBottomSum = eTop() + eBottom();
0024 float eTopBottomDiff = eTop() - eBottom();
0025 return eTopBottomSum != 0 ? eTopBottomDiff / eTopBottomSum : 0.;
0026 }
0027
0028 float EcalRegressionData::subClusRawEnergy(size_t clusNr) const {
0029 if (clusNr < subClusRawEnergy_.size())
0030 return subClusRawEnergy_[clusNr];
0031 else
0032 return 0.;
0033 }
0034
0035 float EcalRegressionData::subClusDEta(size_t clusNr) const {
0036 if (clusNr < subClusDEta_.size())
0037 return subClusDEta_[clusNr];
0038 else
0039 return 0.;
0040 }
0041
0042 float EcalRegressionData::subClusDPhi(size_t clusNr) const {
0043 if (clusNr < subClusDPhi_.size())
0044 return subClusDPhi_[clusNr];
0045 else
0046 return 0.;
0047 }
0048
0049 void EcalRegressionData::fill(const reco::SuperCluster& superClus,
0050 const EcalRecHitCollection* ebRecHits,
0051 const EcalRecHitCollection* eeRecHits,
0052 const CaloGeometry* geom,
0053 const CaloTopology* topology,
0054 int nrVertices) {
0055 clear();
0056
0057 const DetId& seedid = superClus.seed()->hitsAndFractions().at(0).first;
0058 isEB_ = (seedid.subdetId() == EcalBarrel);
0059
0060
0061 if (EcalTools::isHGCalDet(seedid.det()))
0062 return;
0063
0064 const EcalRecHitCollection* recHits = isEB_ ? ebRecHits : eeRecHits;
0065
0066 scRawEnergy_ = superClus.rawEnergy();
0067 scCalibEnergy_ = superClus.correctedEnergy();
0068 scPreShowerEnergy_ = superClus.preshowerEnergy();
0069 scEta_ = superClus.eta();
0070 scPhi_ = superClus.phi();
0071 scEtaWidth_ = superClus.etaWidth();
0072 scPhiWidth_ = superClus.phiWidth();
0073 scNrAdditionalClusters_ = superClus.clustersSize() - 1;
0074
0075 seedClusEnergy_ = superClus.seed()->energy();
0076 eMax_ = EcalClusterTools::eMax(*superClus.seed(), recHits);
0077 e2nd_ = EcalClusterTools::e2nd(*superClus.seed(), recHits);
0078 e3x3_ = EcalClusterTools::e3x3(*superClus.seed(), recHits, topology);
0079 eTop_ = EcalClusterTools::eTop(*superClus.seed(), recHits, topology);
0080 eBottom_ = EcalClusterTools::eBottom(*superClus.seed(), recHits, topology);
0081 eLeft_ = EcalClusterTools::eLeft(*superClus.seed(), recHits, topology);
0082 eRight_ = EcalClusterTools::eRight(*superClus.seed(), recHits, topology);
0083 const auto& localCovs = EcalClusterTools::localCovariances(*superClus.seed(), recHits, topology);
0084 sigmaIEtaIEta_ = edm::isNotFinite(localCovs[0]) ? 0. : std::sqrt(localCovs[0]);
0085 sigmaIPhiIPhi_ = edm::isNotFinite(localCovs[2]) ? 0. : std::sqrt(localCovs[2]);
0086
0087 if (sigmaIEtaIEta_ * sigmaIPhiIPhi_ > 0)
0088 sigmaIEtaIPhi_ = localCovs[1] / (sigmaIEtaIEta_ * sigmaIPhiIPhi_);
0089 else if (localCovs[1] > 0)
0090 sigmaIEtaIPhi_ = 1.;
0091 else
0092 sigmaIEtaIPhi_ = -1.;
0093
0094 auto localCoordsFunc = isEB() ? egammaTools::localEcalClusterCoordsEB : egammaTools::localEcalClusterCoordsEE;
0095
0096 float dummy;
0097 localCoordsFunc(
0098 *superClus.seed(), *geom, seedCrysEtaOrX_, seedCrysPhiOrY_, seedCrysIEtaOrIX_, seedCrysIPhiOrIY_, dummy, dummy);
0099
0100 for (auto clus = superClus.clustersBegin() + 1; clus != superClus.clustersEnd(); ++clus) {
0101 const float dEta = (*clus)->eta() - superClus.seed()->eta();
0102 const float dPhi = reco::deltaPhi((*clus)->phi(), superClus.seed()->phi());
0103 const float dR2 = dEta * dEta + dPhi * dPhi;
0104 if (dR2 > maxSubClusDR2_ || maxSubClusDR2_ == 998001.) {
0105 maxSubClusDR2_ = dR2;
0106 maxSubClusDRDEta_ = dEta;
0107 maxSubClusDRDPhi_ = dPhi;
0108 maxSubClusDRRawEnergy_ = (*clus)->energy();
0109 }
0110 subClusRawEnergy_.push_back((*clus)->energy());
0111 subClusDEta_.push_back(dEta);
0112 subClusDPhi_.push_back(dPhi);
0113 }
0114
0115 nrVtx_ = nrVertices;
0116 }
0117
0118 void EcalRegressionData::clear() {
0119 isEB_ = false;
0120 scRawEnergy_ = 0.;
0121 scCalibEnergy_ = 0.;
0122 scPreShowerEnergy_ = 0.;
0123 scEta_ = 0.;
0124 scPhi_ = 0.;
0125 scEtaWidth_ = 0.;
0126 scPhiWidth_ = 0.;
0127 scNrAdditionalClusters_ = 0;
0128
0129 seedClusEnergy_ = 0.;
0130 eMax_ = 0.;
0131 e2nd_ = 0.;
0132 e3x3_ = 0.;
0133 eTop_ = 0.;
0134 eBottom_ = 0.;
0135 eLeft_ = 0.;
0136 eRight_ = 0.;
0137 sigmaIEtaIEta_ = 0.;
0138 sigmaIEtaIPhi_ = 0.;
0139 sigmaIPhiIPhi_ = 0.;
0140
0141 seedCrysPhiOrY_ = 0.;
0142 seedCrysEtaOrX_ = 0.;
0143 seedCrysIEtaOrIX_ = 0;
0144 seedCrysIPhiOrIY_ = 0;
0145
0146 maxSubClusDR2_ = 998001.;
0147 maxSubClusDRDPhi_ = 999.;
0148 maxSubClusDRDEta_ = 999;
0149 maxSubClusDRRawEnergy_ = 0.;
0150
0151 subClusRawEnergy_.clear();
0152 subClusDPhi_.clear();
0153 subClusDEta_.clear();
0154
0155 nrVtx_ = 0;
0156 }
0157
0158 void EcalRegressionData::fillVec(std::vector<float>& inputVec) const {
0159 if (isEB())
0160 fillVecEB_(inputVec);
0161 else
0162 fillVecEE_(inputVec);
0163 }
0164
0165 void EcalRegressionData::fillVecEB_(std::vector<float>& inputVec) const {
0166 inputVec.clear();
0167 inputVec.resize(33);
0168 inputVec[0] = nrVtx();
0169 inputVec[1] = scEta();
0170 inputVec[2] = scPhi();
0171 inputVec[3] = scEtaWidth();
0172 inputVec[4] = scPhiWidth();
0173 inputVec[5] = scSeedR9();
0174 inputVec[6] = seedClusEnergyOverSCRawEnergy();
0175 inputVec[7] = eMaxOverSCRawEnergy();
0176 inputVec[8] = e2ndOverSCRawEnergy();
0177 inputVec[9] = seedLeftRightAsym();
0178 inputVec[10] = seedTopBottomAsym();
0179 inputVec[11] = sigmaIEtaIEta();
0180 inputVec[12] = sigmaIEtaIPhi();
0181 inputVec[13] = sigmaIPhiIPhi();
0182 inputVec[14] = scNrAdditionalClusters();
0183 inputVec[15] = maxSubClusDR();
0184 inputVec[16] = maxSubClusDRDPhi();
0185 inputVec[17] = maxSubClusDRDEta();
0186 inputVec[18] = maxSubClusDRRawEnergyOverSCRawEnergy();
0187 inputVec[19] = subClusRawEnergyOverSCRawEnergy(SubClusNr::C1);
0188 inputVec[20] = subClusRawEnergyOverSCRawEnergy(
0189 SubClusNr::C2);
0190
0191 inputVec[21] = subClusRawEnergyOverSCRawEnergy(SubClusNr::C3);
0192 inputVec[22] = subClusDPhi(SubClusNr::C1);
0193 inputVec[23] = subClusDPhi(SubClusNr::C2);
0194 inputVec[24] = subClusDPhi(SubClusNr::C3);
0195 inputVec[25] = subClusDEta(SubClusNr::C1);
0196 inputVec[26] = subClusDEta(SubClusNr::C2);
0197 inputVec[27] = subClusDEta(SubClusNr::C3);
0198 inputVec[28] = seedCrysEtaOrX();
0199 inputVec[29] = seedCrysPhiOrY();
0200 inputVec[30] = seedCrysIEtaOrIX();
0201 inputVec[31] = seedCrysIPhiOrIY();
0202 inputVec[32] = scCalibEnergy();
0203 }
0204
0205 void EcalRegressionData::fillVecEE_(std::vector<float>& inputVec) const {
0206 inputVec.clear();
0207 inputVec.resize(33);
0208 inputVec[0] = nrVtx();
0209 inputVec[1] = scEta();
0210 inputVec[2] = scPhi();
0211 inputVec[3] = scEtaWidth();
0212 inputVec[4] = scPhiWidth();
0213 inputVec[5] = scSeedR9();
0214 inputVec[6] = seedClusEnergyOverSCRawEnergy();
0215 inputVec[7] = eMaxOverSCRawEnergy();
0216 inputVec[8] = e2ndOverSCRawEnergy();
0217 inputVec[9] = seedLeftRightAsym();
0218 inputVec[10] = seedTopBottomAsym();
0219 inputVec[11] = sigmaIEtaIEta();
0220 inputVec[12] = sigmaIEtaIPhi();
0221 inputVec[13] = sigmaIPhiIPhi();
0222 inputVec[14] = scNrAdditionalClusters();
0223 inputVec[15] = maxSubClusDR();
0224 inputVec[16] = maxSubClusDRDPhi();
0225 inputVec[17] = maxSubClusDRDEta();
0226 inputVec[18] = maxSubClusDRRawEnergyOverSCRawEnergy();
0227 inputVec[19] = subClusRawEnergyOverSCRawEnergy(SubClusNr::C1);
0228 inputVec[20] = subClusRawEnergyOverSCRawEnergy(SubClusNr::C2);
0229 inputVec[21] = subClusRawEnergyOverSCRawEnergy(SubClusNr::C3);
0230 inputVec[22] = subClusDPhi(SubClusNr::C1);
0231 inputVec[23] = subClusDPhi(SubClusNr::C2);
0232 inputVec[24] = subClusDPhi(SubClusNr::C3);
0233 inputVec[25] = subClusDEta(SubClusNr::C1);
0234 inputVec[26] = subClusDEta(SubClusNr::C2);
0235 inputVec[27] = subClusDEta(SubClusNr::C3);
0236 inputVec[28] = scPreShowerEnergyOverSCRawEnergy();
0237 inputVec[29] = scCalibEnergy();
0238 }