File indexing completed on 2023-10-25 09:59:02
0001 #ifndef RecoEcal_EgammaCoreTools_EcalClusterTools_h
0002 #define RecoEcal_EgammaCoreTools_EcalClusterTools_h
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0031 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0032 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0033 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0034
0035 #include "DataFormats/Math/interface/Vector3D.h"
0036
0037 #include <vector>
0038 #include <cmath>
0039 #include <TMath.h>
0040 #include <TMatrixT.h>
0041 #include <TMatrixD.h>
0042 #include <TVectorT.h>
0043 #include <TVectorD.h>
0044
0045 #include "DataFormats/DetId/interface/DetId.h"
0046 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0047 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0048 #include "RecoCaloTools/Navigation/interface/CaloRectangle.h"
0049
0050 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0051
0052 #include "CLHEP/Geometry/Transform3D.h"
0053
0054 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0055 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0056 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0057 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0058 #include "CondFormats/EcalObjects/interface/EcalPFRecHitThresholds.h"
0059 #include "CondFormats/DataRecord/interface/EcalPFRecHitThresholdsRcd.h"
0060 #include "RecoEcal/EgammaCoreTools/interface/EgammaLocalCovParamDefaults.h"
0061
0062 class DetId;
0063 class CaloTopology;
0064 class CaloGeometry;
0065
0066 struct Cluster2ndMoments {
0067
0068 float sMaj;
0069 float sMin;
0070
0071 float alpha;
0072 Cluster2ndMoments() : sMaj(0.), sMin(0.), alpha(0.) {}
0073 };
0074
0075 template <bool noZS>
0076 class EcalClusterToolsT {
0077 public:
0078
0079
0080
0081
0082 static float e1x3(const reco::BasicCluster &cluster,
0083 const EcalRecHitCollection *recHits,
0084 const CaloTopology *topology);
0085
0086 static float e3x1(const reco::BasicCluster &cluster,
0087 const EcalRecHitCollection *recHits,
0088 const CaloTopology *topology);
0089
0090 static float e1x5(const reco::BasicCluster &cluster,
0091 const EcalRecHitCollection *recHits,
0092 const CaloTopology *topology);
0093
0094 static float e5x1(const reco::BasicCluster &cluster,
0095 const EcalRecHitCollection *recHits,
0096 const CaloTopology *topology);
0097
0098 static float e2x2(const reco::BasicCluster &cluster,
0099 const EcalRecHitCollection *recHits,
0100 const CaloTopology *topology);
0101
0102 static float e3x2(const reco::BasicCluster &cluster,
0103 const EcalRecHitCollection *recHits,
0104 const CaloTopology *topology);
0105
0106 static float e3x3(const reco::BasicCluster &cluster,
0107 const EcalRecHitCollection *recHits,
0108 const CaloTopology *topology);
0109
0110 static float e4x4(const reco::BasicCluster &cluster,
0111 const EcalRecHitCollection *recHits,
0112 const CaloTopology *topology);
0113
0114 static float e5x5(const reco::BasicCluster &cluster,
0115 const EcalRecHitCollection *recHits,
0116 const CaloTopology *topology);
0117 static int n5x5(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits, const CaloTopology *topology);
0118
0119
0120
0121 static float e2x5Right(const reco::BasicCluster &cluster,
0122 const EcalRecHitCollection *recHits,
0123 const CaloTopology *topology);
0124
0125
0126 static float e2x5Left(const reco::BasicCluster &cluster,
0127 const EcalRecHitCollection *recHits,
0128 const CaloTopology *topology);
0129
0130
0131
0132 static float e2x5Top(const reco::BasicCluster &cluster,
0133 const EcalRecHitCollection *recHits,
0134 const CaloTopology *topology);
0135
0136
0137 static float e2x5Bottom(const reco::BasicCluster &cluster,
0138 const EcalRecHitCollection *recHits,
0139 const CaloTopology *topology);
0140
0141
0142
0143 static float e2x5Max(const reco::BasicCluster &cluster,
0144 const EcalRecHitCollection *recHits,
0145 const CaloTopology *topology);
0146
0147
0148 static float eLeft(const reco::BasicCluster &cluster,
0149 const EcalRecHitCollection *recHits,
0150 const CaloTopology *topology);
0151
0152 static float eRight(const reco::BasicCluster &cluster,
0153 const EcalRecHitCollection *recHits,
0154 const CaloTopology *topology);
0155
0156 static float eTop(const reco::BasicCluster &cluster,
0157 const EcalRecHitCollection *recHits,
0158 const CaloTopology *topology);
0159
0160 static float eBottom(const reco::BasicCluster &cluster,
0161 const EcalRecHitCollection *recHits,
0162 const CaloTopology *topology);
0163
0164
0165 static float eMax(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits);
0166
0167
0168 static float e2nd(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits);
0169
0170
0171 static std::pair<DetId, float> getMaximum(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits);
0172
0173 static std::vector<float> energyBasketFractionEta(const reco::BasicCluster &cluster,
0174 const EcalRecHitCollection *recHits);
0175
0176 static std::vector<float> energyBasketFractionPhi(const reco::BasicCluster &cluster,
0177 const EcalRecHitCollection *recHits);
0178
0179
0180 static std::vector<float> lat(const reco::BasicCluster &cluster,
0181 const EcalRecHitCollection *recHits,
0182 const CaloGeometry *geometry,
0183 bool logW = true,
0184 float w0 = 4.7);
0185
0186
0187 static std::array<float, 3> covariances(const reco::BasicCluster &cluster,
0188 const EcalRecHitCollection *recHits,
0189 const CaloTopology *topology,
0190 const CaloGeometry *geometry,
0191 float w0 = 4.7);
0192
0193
0194
0195
0196
0197
0198
0199
0200
0201 static std::array<float, 3> localCovariances(const reco::BasicCluster &cluster,
0202 const EcalRecHitCollection *recHits,
0203 const CaloTopology *topology,
0204 float w0 = EgammaLocalCovParamDefaults::kRelEnCut,
0205 const EcalPFRecHitThresholds *thresholds = nullptr,
0206 float multEB = 0.0,
0207 float multEE = 0.0);
0208
0209 static std::array<float, 3> scLocalCovariances(const reco::SuperCluster &cluster,
0210 const EcalRecHitCollection *recHits,
0211 const CaloTopology *topology,
0212 float w0 = 4.7);
0213
0214
0215 static Cluster2ndMoments cluster2ndMoments(const reco::BasicCluster &basicCluster,
0216 const EcalRecHitCollection &recHits,
0217 double phiCorrectionFactor = 0.8,
0218 double w0 = 4.7,
0219 bool useLogWeights = true);
0220
0221 static Cluster2ndMoments cluster2ndMoments(const reco::SuperCluster &superCluster,
0222 const EcalRecHitCollection &recHits,
0223 double phiCorrectionFactor = 0.8,
0224 double w0 = 4.7,
0225 bool useLogWeights = true);
0226 static Cluster2ndMoments cluster2ndMoments(const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs,
0227 double phiCorrectionFactor = 0.8,
0228 double w0 = 4.7,
0229 bool useLogWeights = true);
0230
0231 static double zernike20(const reco::BasicCluster &cluster,
0232 const EcalRecHitCollection *recHits,
0233 const CaloGeometry *geometry,
0234 double R0 = 6.6,
0235 bool logW = true,
0236 float w0 = 4.7);
0237 static double zernike42(const reco::BasicCluster &cluster,
0238 const EcalRecHitCollection *recHits,
0239 const CaloGeometry *geometry,
0240 double R0 = 6.6,
0241 bool logW = true,
0242 float w0 = 4.7);
0243
0244
0245
0246 static std::vector<DetId> matrixDetId(const CaloTopology *topology, DetId id, CaloRectangle rectangle);
0247 static std::vector<DetId> matrixDetId(const CaloTopology *topology, DetId id, int size) {
0248 return matrixDetId(topology, id, CaloRectangle{-size, size, -size, size});
0249 }
0250
0251
0252
0253 static float matrixEnergy(const reco::BasicCluster &cluster,
0254 const EcalRecHitCollection *recHits,
0255 const CaloTopology *topology,
0256 DetId id,
0257 CaloRectangle rectangle);
0258 static int matrixSize(const reco::BasicCluster &cluster,
0259 const EcalRecHitCollection *recHits,
0260 const CaloTopology *topology,
0261 DetId id,
0262 CaloRectangle rectangle);
0263
0264 static float getFraction(const std::vector<std::pair<DetId, float>> &v_id, DetId id);
0265
0266 static std::pair<DetId, float> getMaximum(const std::vector<std::pair<DetId, float>> &v_id,
0267 const EcalRecHitCollection *recHits);
0268
0269
0270 static float recHitEnergy(DetId id, const EcalRecHitCollection *recHits);
0271
0272
0273 static std::vector<float> roundnessBarrelSuperClusters(const reco::SuperCluster &superCluster,
0274 const EcalRecHitCollection &recHits,
0275 int weightedPositionMethod = 0,
0276 float energyThreshold = 0.0);
0277 static std::vector<float> roundnessBarrelSuperClustersUserExtended(const reco::SuperCluster &superCluster,
0278 const EcalRecHitCollection &recHits,
0279 int ieta_delta = 0,
0280 int iphi_delta = 0,
0281 float energyRHThresh = 0.00000,
0282 int weightedPositionMethod = 0);
0283 static std::vector<float> roundnessSelectedBarrelRecHits(
0284 const std::vector<std::pair<const EcalRecHit *, float>> &rhVector, int weightedPositionMethod = 0);
0285
0286
0287 static int nrSaturatedCrysIn5x5(const DetId &id, const EcalRecHitCollection *recHits, const CaloTopology *topology);
0288
0289 private:
0290 struct EcalClusterEnergyDeposition {
0291 double deposited_energy;
0292 double r;
0293 double phi;
0294 };
0295
0296 static std::vector<EcalClusterEnergyDeposition> getEnergyDepTopology(const reco::BasicCluster &cluster,
0297 const EcalRecHitCollection *recHits,
0298 const CaloGeometry *geometry,
0299 bool logW,
0300 float w0);
0301
0302 static math::XYZVector meanClusterPosition(const reco::BasicCluster &cluster,
0303 const EcalRecHitCollection *recHits,
0304 const CaloTopology *topology,
0305 const CaloGeometry *geometry);
0306
0307
0308
0309 static std::pair<float, float> mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster,
0310 const EcalRecHitCollection *recHits,
0311 const CaloTopology *topology);
0312
0313 static std::pair<float, float> mean5x5PositionInXY(const reco::BasicCluster &cluster,
0314 const EcalRecHitCollection *recHits,
0315 const CaloTopology *topology);
0316
0317 static double f00(double r) { return 1; }
0318 static double f11(double r) { return r; }
0319 static double f20(double r) { return 2.0 * r * r - 1.0; }
0320 static double f22(double r) { return r * r; }
0321 static double f31(double r) { return 3.0 * r * r * r - 2.0 * r; }
0322 static double f33(double r) { return r * r * r; }
0323 static double f40(double r) { return 6.0 * r * r * r * r - 6.0 * r * r + 1.0; }
0324 static double f42(double r) { return 4.0 * r * r * r * r - 3.0 * r * r; }
0325 static double f44(double r) { return r * r * r * r; }
0326 static double f51(double r) { return 10.0 * pow(r, 5) - 12.0 * pow(r, 3) + 3.0 * r; }
0327 static double f53(double r) { return 5.0 * pow(r, 5) - 4.0 * pow(r, 3); }
0328 static double f55(double r) { return pow(r, 5); }
0329
0330 static double absZernikeMoment(const reco::BasicCluster &cluster,
0331 const EcalRecHitCollection *recHits,
0332 const CaloGeometry *geometry,
0333 int n,
0334 int m,
0335 double R0,
0336 bool logW,
0337 float w0);
0338 static double fast_AbsZernikeMoment(const reco::BasicCluster &cluster,
0339 const EcalRecHitCollection *recHits,
0340 const CaloGeometry *geometry,
0341 int n,
0342 int m,
0343 double R0,
0344 bool logW,
0345 float w0);
0346 static double calc_AbsZernikeMoment(const reco::BasicCluster &cluster,
0347 const EcalRecHitCollection *recHits,
0348 const CaloGeometry *geometry,
0349 int n,
0350 int m,
0351 double R0,
0352 bool logW,
0353 float w0);
0354
0355 static double factorial(int n) {
0356 double res = 1.;
0357 for (int i = 2; i <= n; ++i)
0358 res *= i;
0359 return res;
0360 }
0361
0362
0363 static float getIEta(const DetId &id);
0364 static float getIPhi(const DetId &id);
0365 static float getNormedIX(const DetId &id);
0366 static float getNormedIY(const DetId &id);
0367 static float getDPhiEndcap(const DetId &crysId, float meanX, float meanY);
0368 static float getNrCrysDiffInEta(const DetId &crysId, const DetId &orginId);
0369 static float getNrCrysDiffInPhi(const DetId &crysId, const DetId &orginId);
0370
0371
0372 static int deltaIEta(int seed_ieta, int rh_ieta);
0373 static int deltaIPhi(int seed_iphi, int rh_iphi);
0374 static std::vector<int> getSeedPosition(const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs);
0375 static float getSumEnergy(const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs);
0376 static float computeWeight(float eRH, float energyTotal, int weightedPositionMethod);
0377 };
0378
0379
0380 template <bool noZS>
0381 float EcalClusterToolsT<noZS>::getFraction(const std::vector<std::pair<DetId, float>> &v_id, DetId id) {
0382 if (noZS)
0383 return 1.0;
0384 float frac = 0.0;
0385 for (size_t i = 0; i < v_id.size(); ++i) {
0386 if (v_id[i].first.rawId() == id.rawId()) {
0387 frac = v_id[i].second;
0388 break;
0389 }
0390 }
0391 return frac;
0392 }
0393
0394 template <bool noZS>
0395 std::pair<DetId, float> EcalClusterToolsT<noZS>::getMaximum(const std::vector<std::pair<DetId, float>> &v_id,
0396 const EcalRecHitCollection *recHits) {
0397 float max = 0;
0398 DetId id(0);
0399 for (size_t i = 0; i < v_id.size(); ++i) {
0400 float energy = recHitEnergy(v_id[i].first, recHits) * (noZS ? 1.0 : v_id[i].second);
0401 if (energy > max) {
0402 max = energy;
0403 id = v_id[i].first;
0404 }
0405 }
0406 return std::pair<DetId, float>(id, max);
0407 }
0408
0409 template <bool noZS>
0410 std::pair<DetId, float> EcalClusterToolsT<noZS>::getMaximum(const reco::BasicCluster &cluster,
0411 const EcalRecHitCollection *recHits) {
0412 return getMaximum(cluster.hitsAndFractions(), recHits);
0413 }
0414
0415 template <bool noZS>
0416 float EcalClusterToolsT<noZS>::recHitEnergy(DetId id, const EcalRecHitCollection *recHits) {
0417 if (id == DetId(0)) {
0418 return 0;
0419 } else {
0420 EcalRecHitCollection::const_iterator it = recHits->find(id);
0421 if (it != recHits->end()) {
0422 if (noZS && (it->checkFlag(EcalRecHit::kTowerRecovered) || it->checkFlag(EcalRecHit::kWeird) ||
0423 (it->detid().subdetId() == EcalBarrel && it->checkFlag(EcalRecHit::kDiWeird))))
0424 return 0.0;
0425 else
0426 return it->energy();
0427 } else {
0428
0429
0430 return 0;
0431 }
0432 }
0433 return 0;
0434 }
0435
0436
0437
0438
0439
0440
0441
0442
0443
0444
0445
0446
0447
0448 template <bool noZS>
0449 float EcalClusterToolsT<noZS>::matrixEnergy(const reco::BasicCluster &cluster,
0450 const EcalRecHitCollection *recHits,
0451 const CaloTopology *topology,
0452 DetId id,
0453 CaloRectangle rectangle) {
0454 float energy = 0;
0455 auto const &vId = cluster.hitsAndFractions();
0456
0457 for (auto const &detId : rectangle(id, *topology)) {
0458 energy += recHitEnergy(detId, recHits) * getFraction(vId, detId);
0459 }
0460
0461 return energy;
0462 }
0463
0464 template <bool noZS>
0465 int EcalClusterToolsT<noZS>::matrixSize(const reco::BasicCluster &cluster,
0466 const EcalRecHitCollection *recHits,
0467 const CaloTopology *topology,
0468 DetId id,
0469 CaloRectangle rectangle) {
0470 int result = 0;
0471
0472 for (auto const &detId : rectangle(id, *topology)) {
0473 const float energy = recHitEnergy(detId, recHits);
0474 const float frac = getFraction(cluster.hitsAndFractions(), detId);
0475 if (energy * frac > 0)
0476 result++;
0477 }
0478 return result;
0479 }
0480
0481 template <bool noZS>
0482 std::vector<DetId> EcalClusterToolsT<noZS>::matrixDetId(const CaloTopology *topology,
0483 DetId id,
0484 CaloRectangle rectangle) {
0485 std::vector<DetId> v;
0486 for (auto const &detId : rectangle(id, *topology)) {
0487 if (detId != DetId(0))
0488 v.push_back(detId);
0489 }
0490 return v;
0491 }
0492
0493 template <bool noZS>
0494 float EcalClusterToolsT<noZS>::e2x2(const reco::BasicCluster &cluster,
0495 const EcalRecHitCollection *recHits,
0496 const CaloTopology *topology) {
0497 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0498 std::list<float> energies;
0499 float max_E = matrixEnergy(cluster, recHits, topology, id, {-1, 0, -1, 0});
0500 max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-1, 0, 0, 1}));
0501 max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {0, 1, 0, 1}));
0502 max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {0, 1, -1, 0}));
0503 return max_E;
0504 }
0505
0506 template <bool noZS>
0507 float EcalClusterToolsT<noZS>::e3x2(const reco::BasicCluster &cluster,
0508 const EcalRecHitCollection *recHits,
0509 const CaloTopology *topology) {
0510 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0511 float max_E = matrixEnergy(cluster, recHits, topology, id, {-1, 1, -1, 0});
0512 max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {0, 1, -1, 1}));
0513 max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-1, 1, 0, 1}));
0514 max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-1, 0, -1, 1}));
0515 return max_E;
0516 }
0517
0518 template <bool noZS>
0519 float EcalClusterToolsT<noZS>::e3x3(const reco::BasicCluster &cluster,
0520 const EcalRecHitCollection *recHits,
0521 const CaloTopology *topology) {
0522 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0523 return matrixEnergy(cluster, recHits, topology, id, {-1, 1, -1, 1});
0524 }
0525
0526 template <bool noZS>
0527 float EcalClusterToolsT<noZS>::e4x4(const reco::BasicCluster &cluster,
0528 const EcalRecHitCollection *recHits,
0529 const CaloTopology *topology) {
0530 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0531 float max_E = matrixEnergy(cluster, recHits, topology, id, {-1, 2, -2, 1});
0532 max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-2, 1, -2, 1}));
0533 max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-2, 1, -1, 2}));
0534 max_E = std::max(max_E, matrixEnergy(cluster, recHits, topology, id, {-1, 2, -1, 2}));
0535 return max_E;
0536 }
0537
0538 template <bool noZS>
0539 float EcalClusterToolsT<noZS>::e5x5(const reco::BasicCluster &cluster,
0540 const EcalRecHitCollection *recHits,
0541 const CaloTopology *topology) {
0542 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0543 return matrixEnergy(cluster, recHits, topology, id, {-2, 2, -2, 2});
0544 }
0545
0546 template <bool noZS>
0547 int EcalClusterToolsT<noZS>::n5x5(const reco::BasicCluster &cluster,
0548 const EcalRecHitCollection *recHits,
0549 const CaloTopology *topology) {
0550 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0551 return matrixSize(cluster, recHits, topology, id, {-2, 2, -2, 2});
0552 }
0553
0554 template <bool noZS>
0555 float EcalClusterToolsT<noZS>::eMax(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits) {
0556 return getMaximum(cluster.hitsAndFractions(), recHits).second;
0557 }
0558
0559 template <bool noZS>
0560 float EcalClusterToolsT<noZS>::e2nd(const reco::BasicCluster &cluster, const EcalRecHitCollection *recHits) {
0561 std::vector<float> energies;
0562 const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
0563 energies.reserve(v_id.size());
0564 if (v_id.size() < 2)
0565 return 0;
0566 for (size_t i = 0; i < v_id.size(); ++i) {
0567 energies.push_back(recHitEnergy(v_id[i].first, recHits) * (noZS ? 1.0 : v_id[i].second));
0568 }
0569 std::partial_sort(energies.begin(), energies.begin() + 2, energies.end(), std::greater<float>());
0570 return energies[1];
0571 }
0572
0573 template <bool noZS>
0574 float EcalClusterToolsT<noZS>::e2x5Right(const reco::BasicCluster &cluster,
0575 const EcalRecHitCollection *recHits,
0576 const CaloTopology *topology) {
0577 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0578 return matrixEnergy(cluster, recHits, topology, id, {1, 2, -2, 2});
0579 }
0580
0581 template <bool noZS>
0582 float EcalClusterToolsT<noZS>::e2x5Left(const reco::BasicCluster &cluster,
0583 const EcalRecHitCollection *recHits,
0584 const CaloTopology *topology) {
0585 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0586 return matrixEnergy(cluster, recHits, topology, id, {-2, -1, -2, 2});
0587 }
0588
0589 template <bool noZS>
0590 float EcalClusterToolsT<noZS>::e2x5Top(const reco::BasicCluster &cluster,
0591 const EcalRecHitCollection *recHits,
0592 const CaloTopology *topology) {
0593 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0594 return matrixEnergy(cluster, recHits, topology, id, {-2, 2, 1, 2});
0595 }
0596
0597 template <bool noZS>
0598 float EcalClusterToolsT<noZS>::e2x5Bottom(const reco::BasicCluster &cluster,
0599 const EcalRecHitCollection *recHits,
0600 const CaloTopology *topology) {
0601 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0602 return matrixEnergy(cluster, recHits, topology, id, {-2, 2, -2, -1});
0603 }
0604
0605
0606
0607 template <bool noZS>
0608 float EcalClusterToolsT<noZS>::e2x5Max(const reco::BasicCluster &cluster,
0609 const EcalRecHitCollection *recHits,
0610 const CaloTopology *topology) {
0611 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0612
0613
0614 float left = matrixEnergy(cluster, recHits, topology, id, {-1, -1, -2, 2});
0615
0616 float right = matrixEnergy(cluster, recHits, topology, id, {1, 1, -2, 2});
0617
0618 float centre = matrixEnergy(cluster, recHits, topology, id, {0, 0, -2, 2});
0619
0620
0621 return left > right ? left + centre : right + centre;
0622 }
0623
0624 template <bool noZS>
0625 float EcalClusterToolsT<noZS>::e1x5(const reco::BasicCluster &cluster,
0626 const EcalRecHitCollection *recHits,
0627 const CaloTopology *topology) {
0628 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0629 return matrixEnergy(cluster, recHits, topology, id, {0, 0, -2, 2});
0630 }
0631
0632 template <bool noZS>
0633 float EcalClusterToolsT<noZS>::e5x1(const reco::BasicCluster &cluster,
0634 const EcalRecHitCollection *recHits,
0635 const CaloTopology *topology) {
0636 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0637 return matrixEnergy(cluster, recHits, topology, id, {-2, 2, 0, 0});
0638 }
0639
0640 template <bool noZS>
0641 float EcalClusterToolsT<noZS>::e1x3(const reco::BasicCluster &cluster,
0642 const EcalRecHitCollection *recHits,
0643 const CaloTopology *topology) {
0644 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0645 return matrixEnergy(cluster, recHits, topology, id, {0, 0, -1, 1});
0646 }
0647
0648 template <bool noZS>
0649 float EcalClusterToolsT<noZS>::e3x1(const reco::BasicCluster &cluster,
0650 const EcalRecHitCollection *recHits,
0651 const CaloTopology *topology) {
0652 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0653 return matrixEnergy(cluster, recHits, topology, id, {-1, 1, 0, 0});
0654 }
0655
0656 template <bool noZS>
0657 float EcalClusterToolsT<noZS>::eLeft(const reco::BasicCluster &cluster,
0658 const EcalRecHitCollection *recHits,
0659 const CaloTopology *topology) {
0660 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0661 return matrixEnergy(cluster, recHits, topology, id, {-1, -1, 0, 0});
0662 }
0663
0664 template <bool noZS>
0665 float EcalClusterToolsT<noZS>::eRight(const reco::BasicCluster &cluster,
0666 const EcalRecHitCollection *recHits,
0667 const CaloTopology *topology) {
0668 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0669 return matrixEnergy(cluster, recHits, topology, id, {1, 1, 0, 0});
0670 }
0671
0672 template <bool noZS>
0673 float EcalClusterToolsT<noZS>::eTop(const reco::BasicCluster &cluster,
0674 const EcalRecHitCollection *recHits,
0675 const CaloTopology *topology) {
0676 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0677 return matrixEnergy(cluster, recHits, topology, id, {0, 0, 1, 1});
0678 }
0679
0680 template <bool noZS>
0681 float EcalClusterToolsT<noZS>::eBottom(const reco::BasicCluster &cluster,
0682 const EcalRecHitCollection *recHits,
0683 const CaloTopology *topology) {
0684 DetId id = getMaximum(cluster.hitsAndFractions(), recHits).first;
0685 return matrixEnergy(cluster, recHits, topology, id, {0, 0, -1, -1});
0686 }
0687
0688 template <bool noZS>
0689 std::vector<float> EcalClusterToolsT<noZS>::energyBasketFractionEta(const reco::BasicCluster &cluster,
0690 const EcalRecHitCollection *recHits) {
0691 std::vector<float> basketFraction(2 * EBDetId::kModulesPerSM);
0692 float clusterEnergy = cluster.energy();
0693 const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
0694 if (v_id[0].first.subdetId() != EcalBarrel) {
0695 edm::LogWarning("EcalClusterToolsT<noZS>::energyBasketFractionEta")
0696 << "Trying to get basket fraction for endcap basic-clusters. Basket fractions can be obtained ONLY for barrel "
0697 "basic-clusters. Returning empty vector.";
0698 return basketFraction;
0699 }
0700 for (size_t i = 0; i < v_id.size(); ++i) {
0701 basketFraction[EBDetId(v_id[i].first).im() - 1 + EBDetId(v_id[i].first).positiveZ() * EBDetId::kModulesPerSM] +=
0702 recHitEnergy(v_id[i].first, recHits) * v_id[i].second / clusterEnergy;
0703 }
0704 std::sort(basketFraction.rbegin(), basketFraction.rend());
0705 return basketFraction;
0706 }
0707
0708 template <bool noZS>
0709 std::vector<float> EcalClusterToolsT<noZS>::energyBasketFractionPhi(const reco::BasicCluster &cluster,
0710 const EcalRecHitCollection *recHits) {
0711 std::vector<float> basketFraction(2 * (EBDetId::MAX_IPHI / EBDetId::kCrystalsInPhi));
0712 float clusterEnergy = cluster.energy();
0713 const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
0714 if (v_id[0].first.subdetId() != EcalBarrel) {
0715 edm::LogWarning("EcalClusterToolsT<noZS>::energyBasketFractionPhi")
0716 << "Trying to get basket fraction for endcap basic-clusters. Basket fractions can be obtained ONLY for barrel "
0717 "basic-clusters. Returning empty vector.";
0718 return basketFraction;
0719 }
0720 for (size_t i = 0; i < v_id.size(); ++i) {
0721 basketFraction[(EBDetId(v_id[i].first).iphi() - 1) / EBDetId::kCrystalsInPhi +
0722 EBDetId(v_id[i].first).positiveZ() * EBDetId::kTowersInPhi] +=
0723 recHitEnergy(v_id[i].first, recHits) * (noZS ? 1.0 : v_id[i].second) / clusterEnergy;
0724 }
0725 std::sort(basketFraction.rbegin(), basketFraction.rend());
0726 return basketFraction;
0727 }
0728
0729 template <bool noZS>
0730 std::vector<typename EcalClusterToolsT<noZS>::EcalClusterEnergyDeposition>
0731 EcalClusterToolsT<noZS>::getEnergyDepTopology(const reco::BasicCluster &cluster,
0732 const EcalRecHitCollection *recHits,
0733 const CaloGeometry *geometry,
0734 bool logW,
0735 float w0) {
0736 std::vector<typename EcalClusterToolsT<noZS>::EcalClusterEnergyDeposition> energyDistribution;
0737
0738
0739 CLHEP::Hep3Vector clVect(cluster.position().x(), cluster.position().y(), cluster.position().z());
0740 CLHEP::Hep3Vector clDir(clVect);
0741 clDir *= 1.0 / clDir.mag();
0742
0743 CLHEP::Hep3Vector theta_axis(clDir.y(), -clDir.x(), 0.0);
0744 theta_axis *= 1.0 / theta_axis.mag();
0745 CLHEP::Hep3Vector phi_axis = theta_axis.cross(clDir);
0746
0747 const std::vector<std::pair<DetId, float>> &clusterDetIds = cluster.hitsAndFractions();
0748
0749 EcalClusterEnergyDeposition clEdep;
0750 EcalRecHit testEcalRecHit;
0751 std::vector<std::pair<DetId, float>>::const_iterator posCurrent;
0752
0753 for (posCurrent = clusterDetIds.begin(); posCurrent != clusterDetIds.end(); ++posCurrent) {
0754 EcalRecHitCollection::const_iterator itt = recHits->find((*posCurrent).first);
0755 testEcalRecHit = *itt;
0756
0757 if (((*posCurrent).first != DetId(0)) && (recHits->find((*posCurrent).first) != recHits->end())) {
0758 clEdep.deposited_energy = testEcalRecHit.energy() * (noZS ? 1.0 : (*posCurrent).second);
0759
0760 if (logW) {
0761
0762
0763 double weight = std::max(0.0, w0 + log(std::abs(clEdep.deposited_energy) / cluster.energy()));
0764 if (weight == 0) {
0765 LogDebug("ClusterShapeAlgo") << "Crystal has insufficient energy: E = " << clEdep.deposited_energy
0766 << " GeV; skipping... ";
0767 continue;
0768 } else
0769 LogDebug("ClusterShapeAlgo") << "===> got crystal. Energy = " << clEdep.deposited_energy << " GeV. ";
0770 }
0771 DetId id_ = (*posCurrent).first;
0772 const CaloSubdetectorGeometry *geo = geometry->getSubdetectorGeometry(id_);
0773 auto this_cell = geo->getGeometry(id_);
0774 const GlobalPoint &cellPos = this_cell->getPosition();
0775 CLHEP::Hep3Vector gblPos(cellPos.x(), cellPos.y(), cellPos.z());
0776
0777 CLHEP::Hep3Vector diff = gblPos - clVect;
0778
0779
0780
0781 CLHEP::Hep3Vector DigiVect = diff - diff.dot(clDir) * clDir;
0782 clEdep.r = DigiVect.mag();
0783 LogDebug("ClusterShapeAlgo") << "E = " << clEdep.deposited_energy << "\tdiff = " << diff.mag()
0784 << "\tr = " << clEdep.r;
0785 clEdep.phi = DigiVect.angle(theta_axis);
0786 if (DigiVect.dot(phi_axis) < 0)
0787 clEdep.phi = 2 * M_PI - clEdep.phi;
0788 energyDistribution.push_back(clEdep);
0789 }
0790 }
0791 return energyDistribution;
0792 }
0793
0794 template <bool noZS>
0795 std::vector<float> EcalClusterToolsT<noZS>::lat(const reco::BasicCluster &cluster,
0796 const EcalRecHitCollection *recHits,
0797 const CaloGeometry *geometry,
0798 bool logW,
0799 float w0) {
0800 std::vector<EcalClusterToolsT::EcalClusterEnergyDeposition> energyDistribution =
0801 getEnergyDepTopology(cluster, recHits, geometry, logW, w0);
0802
0803 std::vector<float> lat;
0804 double r, redmoment = 0;
0805 double phiRedmoment = 0;
0806 double etaRedmoment = 0;
0807 int n, n1, n2, tmp;
0808 int clusterSize = energyDistribution.size();
0809 float etaLat_, phiLat_, lat_;
0810 if (clusterSize < 3) {
0811 etaLat_ = 0.0;
0812 lat_ = 0.0;
0813 lat.push_back(0.);
0814 lat.push_back(0.);
0815 lat.push_back(0.);
0816 return lat;
0817 }
0818
0819 n1 = 0;
0820 n2 = 1;
0821 if (energyDistribution[1].deposited_energy > energyDistribution[0].deposited_energy) {
0822 tmp = n2;
0823 n2 = n1;
0824 n1 = tmp;
0825 }
0826 for (int i = 2; i < clusterSize; i++) {
0827 n = i;
0828 if (energyDistribution[i].deposited_energy > energyDistribution[n1].deposited_energy) {
0829 tmp = n2;
0830 n2 = n1;
0831 n1 = i;
0832 n = tmp;
0833 } else {
0834 if (energyDistribution[i].deposited_energy > energyDistribution[n2].deposited_energy) {
0835 tmp = n2;
0836 n2 = i;
0837 n = tmp;
0838 }
0839 }
0840
0841 r = energyDistribution[n].r;
0842 redmoment += r * r * energyDistribution[n].deposited_energy;
0843 double rphi = r * cos(energyDistribution[n].phi);
0844 phiRedmoment += rphi * rphi * energyDistribution[n].deposited_energy;
0845 double reta = r * sin(energyDistribution[n].phi);
0846 etaRedmoment += reta * reta * energyDistribution[n].deposited_energy;
0847 }
0848 double e1 = energyDistribution[n1].deposited_energy;
0849 double e2 = energyDistribution[n2].deposited_energy;
0850
0851 lat_ = redmoment / (redmoment + 2.19 * 2.19 * (e1 + e2));
0852 phiLat_ = phiRedmoment / (phiRedmoment + 2.19 * 2.19 * (e1 + e2));
0853 etaLat_ = etaRedmoment / (etaRedmoment + 2.19 * 2.19 * (e1 + e2));
0854
0855 lat.push_back(etaLat_);
0856 lat.push_back(phiLat_);
0857 lat.push_back(lat_);
0858 return lat;
0859 }
0860
0861 template <bool noZS>
0862 math::XYZVector EcalClusterToolsT<noZS>::meanClusterPosition(const reco::BasicCluster &cluster,
0863 const EcalRecHitCollection *recHits,
0864 const CaloTopology *topology,
0865 const CaloGeometry *geometry) {
0866
0867 math::XYZVector meanPosition(0.0, 0.0, 0.0);
0868 const std::vector<std::pair<DetId, float>> &hsAndFs = cluster.hitsAndFractions();
0869 std::vector<DetId> v_id = matrixDetId(topology, getMaximum(cluster, recHits).first, 2);
0870 for (const std::pair<DetId, float> &hitAndFrac : hsAndFs) {
0871 for (std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it) {
0872 if (hitAndFrac.first != *it && !noZS)
0873 continue;
0874 const CaloSubdetectorGeometry *geo = geometry->getSubdetectorGeometry(*it);
0875 GlobalPoint positionGP = geo->getGeometry(*it)->getPosition();
0876 math::XYZVector position(positionGP.x(), positionGP.y(), positionGP.z());
0877 meanPosition = meanPosition + recHitEnergy(*it, recHits) * position * hitAndFrac.second;
0878 }
0879 if (noZS)
0880 break;
0881 }
0882 return meanPosition / e5x5(cluster, recHits, topology);
0883 }
0884
0885
0886
0887
0888
0889
0890
0891
0892 template <bool noZS>
0893 std::pair<float, float> EcalClusterToolsT<noZS>::mean5x5PositionInLocalCrysCoord(const reco::BasicCluster &cluster,
0894 const EcalRecHitCollection *recHits,
0895 const CaloTopology *topology) {
0896 DetId seedId = getMaximum(cluster, recHits).first;
0897 float meanDEta = 0.;
0898 float meanDPhi = 0.;
0899 float energySum = 0.;
0900
0901 const std::vector<std::pair<DetId, float>> &hsAndFs = cluster.hitsAndFractions();
0902 std::vector<DetId> v_id = matrixDetId(topology, seedId, 2);
0903 for (const std::pair<DetId, float> &hAndF : hsAndFs) {
0904 for (std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it) {
0905 if (hAndF.first != *it && !noZS)
0906 continue;
0907 float energy = recHitEnergy(*it, recHits) * hAndF.second;
0908 if (energy < 0.)
0909 continue;
0910 meanDEta += energy * getNrCrysDiffInEta(*it, seedId);
0911 meanDPhi += energy * getNrCrysDiffInPhi(*it, seedId);
0912 energySum += energy;
0913 }
0914 if (noZS)
0915 break;
0916 }
0917 meanDEta /= energySum;
0918 meanDPhi /= energySum;
0919 return std::pair<float, float>(meanDEta, meanDPhi);
0920 }
0921
0922
0923
0924
0925
0926
0927
0928 template <bool noZS>
0929 std::pair<float, float> EcalClusterToolsT<noZS>::mean5x5PositionInXY(const reco::BasicCluster &cluster,
0930 const EcalRecHitCollection *recHits,
0931 const CaloTopology *topology) {
0932 DetId seedId = getMaximum(cluster, recHits).first;
0933
0934 std::pair<float, float> meanXY(0., 0.);
0935 if (seedId.subdetId() == EcalBarrel)
0936 return meanXY;
0937
0938 float energySum = 0.;
0939
0940 const std::vector<std::pair<DetId, float>> &hsAndFs = cluster.hitsAndFractions();
0941 std::vector<DetId> v_id = matrixDetId(topology, seedId, 2);
0942 for (const std::pair<DetId, float> &hAndF : hsAndFs) {
0943 for (std::vector<DetId>::const_iterator it = v_id.begin(); it != v_id.end(); ++it) {
0944 if (hAndF.first != *it && !noZS)
0945 continue;
0946 float energy = recHitEnergy(*it, recHits) * hAndF.second;
0947 if (energy < 0.)
0948 continue;
0949 meanXY.first += energy * getNormedIX(*it);
0950 meanXY.second += energy * getNormedIY(*it);
0951 energySum += energy;
0952 }
0953 if (noZS)
0954 break;
0955 }
0956 meanXY.first /= energySum;
0957 meanXY.second /= energySum;
0958 return meanXY;
0959 }
0960
0961 template <bool noZS>
0962 std::array<float, 3> EcalClusterToolsT<noZS>::covariances(const reco::BasicCluster &cluster,
0963 const EcalRecHitCollection *recHits,
0964 const CaloTopology *topology,
0965 const CaloGeometry *geometry,
0966 float w0) {
0967 float e_5x5 = e5x5(cluster, recHits, topology);
0968 float covEtaEta, covEtaPhi, covPhiPhi;
0969 if (e_5x5 >= 0.) {
0970
0971 const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
0972 math::XYZVector meanPosition = meanClusterPosition(cluster, recHits, topology, geometry);
0973
0974
0975 double numeratorEtaEta = 0;
0976 double numeratorEtaPhi = 0;
0977 double numeratorPhiPhi = 0;
0978 double denominator = 0;
0979
0980 DetId id = getMaximum(v_id, recHits).first;
0981 CaloRectangle rectangle{-2, 2, -2, 2};
0982 for (auto const &detId : rectangle(id, *topology)) {
0983 float frac = getFraction(v_id, detId);
0984 float energy = recHitEnergy(detId, recHits) * frac;
0985
0986 if (energy <= 0)
0987 continue;
0988
0989 const CaloSubdetectorGeometry *geo = geometry->getSubdetectorGeometry(detId);
0990 GlobalPoint position = geo->getGeometry(detId)->getPosition();
0991
0992 double dPhi = position.phi() - meanPosition.phi();
0993 if (dPhi > +Geom::pi()) {
0994 dPhi = Geom::twoPi() - dPhi;
0995 }
0996 if (dPhi < -Geom::pi()) {
0997 dPhi = Geom::twoPi() + dPhi;
0998 }
0999
1000 double dEta = position.eta() - meanPosition.eta();
1001 double w = 0.;
1002 w = std::max(0.0f, w0 + std::log(energy / e_5x5));
1003
1004 denominator += w;
1005 numeratorEtaEta += w * dEta * dEta;
1006 numeratorEtaPhi += w * dEta * dPhi;
1007 numeratorPhiPhi += w * dPhi * dPhi;
1008 }
1009
1010 if (denominator != 0.0) {
1011 covEtaEta = numeratorEtaEta / denominator;
1012 covEtaPhi = numeratorEtaPhi / denominator;
1013 covPhiPhi = numeratorPhiPhi / denominator;
1014 } else {
1015 covEtaEta = 999.9;
1016 covEtaPhi = 999.9;
1017 covPhiPhi = 999.9;
1018 }
1019
1020 } else {
1021
1022
1023 covEtaEta = 0;
1024 covEtaPhi = 0;
1025 covPhiPhi = 0;
1026 }
1027 std::array<float, 3> v{{covEtaEta, covEtaPhi, covPhiPhi}};
1028 return v;
1029 }
1030
1031
1032
1033
1034
1035 template <bool noZS>
1036 std::array<float, 3> EcalClusterToolsT<noZS>::localCovariances(const reco::BasicCluster &cluster,
1037 const EcalRecHitCollection *recHits,
1038 const CaloTopology *topology,
1039 float w0,
1040 const EcalPFRecHitThresholds *thresholds,
1041 float multEB,
1042 float multEE) {
1043 float e_5x5 = e5x5(cluster, recHits, topology);
1044 float covEtaEta, covEtaPhi, covPhiPhi;
1045
1046 if (e_5x5 >= 0.) {
1047
1048 const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
1049 std::pair<float, float> mean5x5PosInNrCrysFromSeed = mean5x5PositionInLocalCrysCoord(cluster, recHits, topology);
1050 std::pair<float, float> mean5x5XYPos = mean5x5PositionInXY(cluster, recHits, topology);
1051
1052
1053 double numeratorEtaEta = 0;
1054 double numeratorEtaPhi = 0;
1055 double numeratorPhiPhi = 0;
1056 double denominator = 0;
1057
1058
1059
1060 const double barrelCrysSize = 0.01745;
1061 const double endcapCrysSize = 0.0447;
1062
1063 DetId seedId = getMaximum(v_id, recHits).first;
1064
1065 bool isBarrel = seedId.subdetId() == EcalBarrel;
1066 const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1067 float mult = isBarrel ? multEB : multEE;
1068
1069
1070
1071
1072 CaloRectangle rectangle{-2, 2, -2, 2};
1073 for (auto const &detId : rectangle(seedId, *topology)) {
1074 float frac = getFraction(v_id, detId);
1075 float energy = recHitEnergy(detId, recHits) * frac;
1076
1077 if ((thresholds == nullptr) && (mult != 0.0)) {
1078 throw cms::Exception("EmptyPFRechHitThresColl")
1079 << "In EcalClusterTools::localCovariances, if EcalPFRecHitThresholds==nulptr, then multEB and multEE "
1080 "should be 0 as well.";
1081 }
1082 float rhThres = 0.0;
1083 if (thresholds != nullptr) {
1084 rhThres = (*thresholds)[detId];
1085 }
1086
1087 if (energy <= (rhThres * mult))
1088 continue;
1089
1090 float dEta = getNrCrysDiffInEta(detId, seedId) - mean5x5PosInNrCrysFromSeed.first;
1091 float dPhi = 0;
1092
1093 if (isBarrel)
1094 dPhi = getNrCrysDiffInPhi(detId, seedId) - mean5x5PosInNrCrysFromSeed.second;
1095 else
1096 dPhi = getDPhiEndcap(detId, mean5x5XYPos.first, mean5x5XYPos.second);
1097
1098 double w = std::max(0.0f, w0 + std::log(energy / e_5x5));
1099
1100 denominator += w;
1101 numeratorEtaEta += w * dEta * dEta;
1102 numeratorEtaPhi += w * dEta * dPhi;
1103 numeratorPhiPhi += w * dPhi * dPhi;
1104 }
1105
1106
1107 if (denominator != 0.0) {
1108 covEtaEta = crysSize * crysSize * numeratorEtaEta / denominator;
1109 covEtaPhi = crysSize * crysSize * numeratorEtaPhi / denominator;
1110 covPhiPhi = crysSize * crysSize * numeratorPhiPhi / denominator;
1111 } else {
1112 covEtaEta = 999.9;
1113 covEtaPhi = 999.9;
1114 covPhiPhi = 999.9;
1115 }
1116
1117 } else {
1118
1119
1120 covEtaEta = 0;
1121 covEtaPhi = 0;
1122 covPhiPhi = 0;
1123 }
1124 std::array<float, 3> v{{covEtaEta, covEtaPhi, covPhiPhi}};
1125 return v;
1126 }
1127
1128 template <bool noZS>
1129 double EcalClusterToolsT<noZS>::zernike20(const reco::BasicCluster &cluster,
1130 const EcalRecHitCollection *recHits,
1131 const CaloGeometry *geometry,
1132 double R0,
1133 bool logW,
1134 float w0) {
1135 return absZernikeMoment(cluster, recHits, geometry, 2, 0, R0, logW, w0);
1136 }
1137
1138 template <bool noZS>
1139 double EcalClusterToolsT<noZS>::zernike42(const reco::BasicCluster &cluster,
1140 const EcalRecHitCollection *recHits,
1141 const CaloGeometry *geometry,
1142 double R0,
1143 bool logW,
1144 float w0) {
1145 return absZernikeMoment(cluster, recHits, geometry, 4, 2, R0, logW, w0);
1146 }
1147
1148 template <bool noZS>
1149 double EcalClusterToolsT<noZS>::absZernikeMoment(const reco::BasicCluster &cluster,
1150 const EcalRecHitCollection *recHits,
1151 const CaloGeometry *geometry,
1152 int n,
1153 int m,
1154 double R0,
1155 bool logW,
1156 float w0) {
1157
1158 if ((m > n) || ((n - m) % 2 != 0) || (n < 0) || (m < 0))
1159 return -1;
1160
1161
1162
1163 if ((n > 20) || (R0 <= 2.19))
1164 return -1;
1165 if (n <= 5)
1166 return fast_AbsZernikeMoment(cluster, recHits, geometry, n, m, R0, logW, w0);
1167 else
1168 return calc_AbsZernikeMoment(cluster, recHits, geometry, n, m, R0, logW, w0);
1169 }
1170
1171 template <bool noZS>
1172 double EcalClusterToolsT<noZS>::fast_AbsZernikeMoment(const reco::BasicCluster &cluster,
1173 const EcalRecHitCollection *recHits,
1174 const CaloGeometry *geometry,
1175 int n,
1176 int m,
1177 double R0,
1178 bool logW,
1179 float w0) {
1180 double r, ph, e, Re = 0, Im = 0;
1181 double TotalEnergy = cluster.energy();
1182 int index = (n / 2) * (n / 2) + (n / 2) + m;
1183 std::vector<EcalClusterEnergyDeposition> energyDistribution =
1184 getEnergyDepTopology(cluster, recHits, geometry, logW, w0);
1185 int clusterSize = energyDistribution.size();
1186 if (clusterSize < 3)
1187 return 0.0;
1188
1189 for (int i = 0; i < clusterSize; i++) {
1190 r = energyDistribution[i].r / R0;
1191 if (r < 1) {
1192 std::vector<double> pol;
1193 pol.push_back(f00(r));
1194 pol.push_back(f11(r));
1195 pol.push_back(f20(r));
1196 pol.push_back(f22(r));
1197 pol.push_back(f31(r));
1198 pol.push_back(f33(r));
1199 pol.push_back(f40(r));
1200 pol.push_back(f42(r));
1201 pol.push_back(f44(r));
1202 pol.push_back(f51(r));
1203 pol.push_back(f53(r));
1204 pol.push_back(f55(r));
1205 ph = (energyDistribution[i]).phi;
1206 e = energyDistribution[i].deposited_energy;
1207 Re = Re + e / TotalEnergy * pol[index] * cos((double)m * ph);
1208 Im = Im - e / TotalEnergy * pol[index] * sin((double)m * ph);
1209 }
1210 }
1211 return sqrt(Re * Re + Im * Im);
1212 }
1213
1214 template <bool noZS>
1215 double EcalClusterToolsT<noZS>::calc_AbsZernikeMoment(const reco::BasicCluster &cluster,
1216 const EcalRecHitCollection *recHits,
1217 const CaloGeometry *geometry,
1218 int n,
1219 int m,
1220 double R0,
1221 bool logW,
1222 float w0) {
1223 double r, ph, e, Re = 0, Im = 0, f_nm;
1224 double TotalEnergy = cluster.energy();
1225 std::vector<EcalClusterEnergyDeposition> energyDistribution =
1226 getEnergyDepTopology(cluster, recHits, geometry, logW, w0);
1227 int clusterSize = energyDistribution.size();
1228 if (clusterSize < 3)
1229 return 0.0;
1230
1231 for (int i = 0; i < clusterSize; ++i) {
1232 r = energyDistribution[i].r / R0;
1233 if (r < 1) {
1234 ph = energyDistribution[i].phi;
1235 e = energyDistribution[i].deposited_energy;
1236 f_nm = 0;
1237 for (int s = 0; s <= (n - m) / 2; s++) {
1238 if (s % 2 == 0) {
1239 f_nm = f_nm + factorial(n - s) * pow(r, (double)(n - 2 * s)) /
1240 (factorial(s) * factorial((n + m) / 2 - s) * factorial((n - m) / 2 - s));
1241 } else {
1242 f_nm = f_nm - factorial(n - s) * pow(r, (double)(n - 2 * s)) /
1243 (factorial(s) * factorial((n + m) / 2 - s) * factorial((n - m) / 2 - s));
1244 }
1245 }
1246 Re = Re + e / TotalEnergy * f_nm * cos((double)m * ph);
1247 Im = Im - e / TotalEnergy * f_nm * sin((double)m * ph);
1248 }
1249 }
1250 return sqrt(Re * Re + Im * Im);
1251 }
1252
1253
1254
1255
1256
1257 template <bool noZS>
1258 float EcalClusterToolsT<noZS>::getIEta(const DetId &id) {
1259 if (id.det() == DetId::Ecal) {
1260 if (id.subdetId() == EcalBarrel) {
1261 EBDetId ebId(id);
1262 return ebId.ieta();
1263 } else if (id.subdetId() == EcalEndcap) {
1264 float iXNorm = getNormedIX(id);
1265 float iYNorm = getNormedIY(id);
1266
1267 return std::sqrt(iXNorm * iXNorm + iYNorm * iYNorm);
1268 }
1269 }
1270 return 0.;
1271 }
1272
1273
1274
1275
1276
1277 template <bool noZS>
1278 float EcalClusterToolsT<noZS>::getIPhi(const DetId &id) {
1279 if (id.det() == DetId::Ecal) {
1280 if (id.subdetId() == EcalBarrel) {
1281 EBDetId ebId(id);
1282 return ebId.iphi();
1283 }
1284 }
1285 return 0.;
1286 }
1287
1288
1289 template <bool noZS>
1290 float EcalClusterToolsT<noZS>::getNormedIX(const DetId &id) {
1291 if (id.det() == DetId::Ecal && id.subdetId() == EcalEndcap) {
1292 EEDetId eeId(id);
1293 int iXNorm = eeId.ix() - 50;
1294 if (iXNorm <= 0)
1295 iXNorm--;
1296 return iXNorm;
1297 }
1298 return 0;
1299 }
1300
1301
1302 template <bool noZS>
1303 float EcalClusterToolsT<noZS>::getNormedIY(const DetId &id) {
1304 if (id.det() == DetId::Ecal && id.subdetId() == EcalEndcap) {
1305 EEDetId eeId(id);
1306 int iYNorm = eeId.iy() - 50;
1307 if (iYNorm <= 0)
1308 iYNorm--;
1309 return iYNorm;
1310 }
1311 return 0;
1312 }
1313
1314
1315 template <bool noZS>
1316 float EcalClusterToolsT<noZS>::getNrCrysDiffInEta(const DetId &crysId, const DetId &orginId) {
1317 float crysIEta = getIEta(crysId);
1318 float orginIEta = getIEta(orginId);
1319 bool isBarrel = orginId.subdetId() == EcalBarrel;
1320
1321 float nrCrysDiff = crysIEta - orginIEta;
1322
1323
1324
1325 if (isBarrel) {
1326 if (crysIEta * orginIEta < 0) {
1327 if (crysIEta > 0)
1328 nrCrysDiff--;
1329 else
1330 nrCrysDiff++;
1331 }
1332 }
1333 return nrCrysDiff;
1334 }
1335
1336
1337 template <bool noZS>
1338 float EcalClusterToolsT<noZS>::getNrCrysDiffInPhi(const DetId &crysId, const DetId &orginId) {
1339 float crysIPhi = getIPhi(crysId);
1340 float orginIPhi = getIPhi(orginId);
1341 bool isBarrel = orginId.subdetId() == EcalBarrel;
1342
1343 float nrCrysDiff = crysIPhi - orginIPhi;
1344
1345 if (isBarrel) {
1346 if (nrCrysDiff > +180) {
1347 nrCrysDiff = nrCrysDiff - 360;
1348 }
1349 if (nrCrysDiff < -180) {
1350 nrCrysDiff = nrCrysDiff + 360;
1351 }
1352 }
1353 return nrCrysDiff;
1354 }
1355
1356
1357 template <bool noZS>
1358 float EcalClusterToolsT<noZS>::getDPhiEndcap(const DetId &crysId, float meanX, float meanY) {
1359 float iXNorm = getNormedIX(crysId);
1360 float iYNorm = getNormedIY(crysId);
1361
1362 float hitLocalR2 = (iXNorm - meanX) * (iXNorm - meanX) + (iYNorm - meanY) * (iYNorm - meanY);
1363 float hitR2 = iXNorm * iXNorm + iYNorm * iYNorm;
1364 float meanR2 = meanX * meanX + meanY * meanY;
1365 float hitR = sqrt(hitR2);
1366 float meanR = sqrt(meanR2);
1367
1368 float tmp = (hitR2 + meanR2 - hitLocalR2) / (2 * hitR * meanR);
1369 if (tmp < -1)
1370 tmp = -1;
1371 if (tmp > 1)
1372 tmp = 1;
1373 float phi = acos(tmp);
1374 float dPhi = hitR * phi;
1375
1376 return dPhi;
1377 }
1378
1379 template <bool noZS>
1380 std::array<float, 3> EcalClusterToolsT<noZS>::scLocalCovariances(const reco::SuperCluster &cluster,
1381 const EcalRecHitCollection *recHits,
1382 const CaloTopology *topology,
1383 float w0) {
1384 const reco::BasicCluster bcluster = *(cluster.seed());
1385
1386 float e_5x5 = e5x5(bcluster, recHits, topology);
1387 float covEtaEta, covEtaPhi, covPhiPhi;
1388
1389 if (e_5x5 >= 0.) {
1390 const std::vector<std::pair<DetId, float>> &v_id = cluster.hitsAndFractions();
1391 std::pair<float, float> mean5x5PosInNrCrysFromSeed = mean5x5PositionInLocalCrysCoord(bcluster, recHits, topology);
1392 std::pair<float, float> mean5x5XYPos = mean5x5PositionInXY(cluster, recHits, topology);
1393
1394 double numeratorEtaEta = 0;
1395 double numeratorEtaPhi = 0;
1396 double numeratorPhiPhi = 0;
1397 double denominator = 0;
1398
1399 const double barrelCrysSize = 0.01745;
1400 const double endcapCrysSize = 0.0447;
1401
1402 DetId seedId = getMaximum(v_id, recHits).first;
1403 bool isBarrel = seedId.subdetId() == EcalBarrel;
1404
1405 const double crysSize = isBarrel ? barrelCrysSize : endcapCrysSize;
1406
1407 for (size_t i = 0; i < v_id.size(); ++i) {
1408 float frac = getFraction(v_id, v_id[i].first);
1409 float energy = recHitEnergy(v_id[i].first, recHits) * frac;
1410
1411 if (energy <= 0)
1412 continue;
1413
1414 float dEta = getNrCrysDiffInEta(v_id[i].first, seedId) - mean5x5PosInNrCrysFromSeed.first;
1415 float dPhi = 0;
1416 if (isBarrel)
1417 dPhi = getNrCrysDiffInPhi(v_id[i].first, seedId) - mean5x5PosInNrCrysFromSeed.second;
1418 else
1419 dPhi = getDPhiEndcap(v_id[i].first, mean5x5XYPos.first, mean5x5XYPos.second);
1420
1421 double w = 0.;
1422 w = std::max(0.0f, w0 + std::log(energy / e_5x5));
1423
1424 denominator += w;
1425 numeratorEtaEta += w * dEta * dEta;
1426 numeratorEtaPhi += w * dEta * dPhi;
1427 numeratorPhiPhi += w * dPhi * dPhi;
1428 }
1429
1430
1431 if (denominator != 0.0) {
1432 covEtaEta = crysSize * crysSize * numeratorEtaEta / denominator;
1433 covEtaPhi = crysSize * crysSize * numeratorEtaPhi / denominator;
1434 covPhiPhi = crysSize * crysSize * numeratorPhiPhi / denominator;
1435 } else {
1436 covEtaEta = 999.9;
1437 covEtaPhi = 999.9;
1438 covPhiPhi = 999.9;
1439 }
1440
1441 } else {
1442
1443
1444 covEtaEta = 0;
1445 covEtaPhi = 0;
1446 covPhiPhi = 0;
1447 }
1448
1449 std::array<float, 3> v{{covEtaEta, covEtaPhi, covPhiPhi}};
1450 return v;
1451 }
1452
1453
1454
1455
1456
1457
1458 template <bool noZS>
1459 Cluster2ndMoments EcalClusterToolsT<noZS>::cluster2ndMoments(const reco::BasicCluster &basicCluster,
1460 const EcalRecHitCollection &recHits,
1461 double phiCorrectionFactor,
1462 double w0,
1463 bool useLogWeights) {
1464 if (recHits.empty())
1465 return Cluster2ndMoments();
1466
1467 std::vector<std::pair<const EcalRecHit *, float>> RH_ptrs_fracs;
1468
1469 const std::vector<std::pair<DetId, float>> &myHitsPair = basicCluster.hitsAndFractions();
1470
1471 for (unsigned int i = 0; i < myHitsPair.size(); i++) {
1472
1473 EcalRecHitCollection::const_iterator myRH = recHits.find(myHitsPair[i].first);
1474 if (myRH != recHits.end()) {
1475 RH_ptrs_fracs.push_back(std::make_pair(&(*myRH), myHitsPair[i].second));
1476 }
1477 }
1478
1479 return EcalClusterToolsT<noZS>::cluster2ndMoments(RH_ptrs_fracs, phiCorrectionFactor, w0, useLogWeights);
1480 }
1481
1482 template <bool noZS>
1483 Cluster2ndMoments EcalClusterToolsT<noZS>::cluster2ndMoments(const reco::SuperCluster &superCluster,
1484 const EcalRecHitCollection &recHits,
1485 double phiCorrectionFactor,
1486 double w0,
1487 bool useLogWeights) {
1488 return EcalClusterToolsT<noZS>::cluster2ndMoments(
1489 *(superCluster.seed()), recHits, phiCorrectionFactor, w0, useLogWeights);
1490 }
1491
1492 template <bool noZS>
1493 Cluster2ndMoments EcalClusterToolsT<noZS>::cluster2ndMoments(
1494 const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs,
1495 double phiCorrectionFactor,
1496 double w0,
1497 bool useLogWeights) {
1498 if (RH_ptrs_fracs.empty())
1499 return Cluster2ndMoments();
1500
1501 double mid_eta(0), mid_phi(0), mid_x(0), mid_y(0);
1502
1503 double Etot = EcalClusterToolsT<noZS>::getSumEnergy(RH_ptrs_fracs);
1504
1505 double max_phi = -10.;
1506 double min_phi = 100.;
1507
1508 std::vector<double> etaDetId;
1509 std::vector<double> phiDetId;
1510 std::vector<double> xDetId;
1511 std::vector<double> yDetId;
1512 std::vector<double> wiDetId;
1513
1514 unsigned int nCry = 0;
1515 double denominator = 0.;
1516 bool isBarrel(true);
1517
1518
1519 for (std::vector<std::pair<const EcalRecHit *, float>>::const_iterator rhf_ptr = RH_ptrs_fracs.begin();
1520 rhf_ptr != RH_ptrs_fracs.end();
1521 rhf_ptr++) {
1522 const EcalRecHit *rh_ptr = rhf_ptr->first;
1523
1524
1525 double temp_eta(0), temp_phi(0), temp_x(0), temp_y(0);
1526 isBarrel = rh_ptr->detid().subdetId() == EcalBarrel;
1527
1528 if (isBarrel) {
1529 temp_eta = (getIEta(rh_ptr->detid()) > 0. ? getIEta(rh_ptr->detid()) + 84.5 : getIEta(rh_ptr->detid()) + 85.5);
1530 temp_phi = getIPhi(rh_ptr->detid()) - 0.5;
1531 } else {
1532 temp_eta = getIEta(rh_ptr->detid());
1533 temp_x = getNormedIX(rh_ptr->detid());
1534 temp_y = getNormedIY(rh_ptr->detid());
1535 }
1536
1537 double temp_ene = rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
1538
1539 double temp_wi = ((useLogWeights) ? std::max(0.0, w0 + std::log(std::abs(temp_ene) / Etot)) : temp_ene);
1540
1541 if (temp_phi > max_phi)
1542 max_phi = temp_phi;
1543 if (temp_phi < min_phi)
1544 min_phi = temp_phi;
1545 etaDetId.push_back(temp_eta);
1546 phiDetId.push_back(temp_phi);
1547 xDetId.push_back(temp_x);
1548 yDetId.push_back(temp_y);
1549 wiDetId.push_back(temp_wi);
1550 denominator += temp_wi;
1551 nCry++;
1552 }
1553
1554 if (isBarrel) {
1555
1556 if (max_phi == 359.5 && min_phi == 0.5) {
1557 for (unsigned int i = 0; i < nCry; i++) {
1558 if (phiDetId[i] - 179. > 0.)
1559 phiDetId[i] -= 360.;
1560 mid_phi += phiDetId[i] * wiDetId[i];
1561 mid_eta += etaDetId[i] * wiDetId[i];
1562 }
1563 } else {
1564 for (unsigned int i = 0; i < nCry; i++) {
1565 mid_phi += phiDetId[i] * wiDetId[i];
1566 mid_eta += etaDetId[i] * wiDetId[i];
1567 }
1568 }
1569 } else {
1570 for (unsigned int i = 0; i < nCry; i++) {
1571 mid_eta += etaDetId[i] * wiDetId[i];
1572 mid_x += xDetId[i] * wiDetId[i];
1573 mid_y += yDetId[i] * wiDetId[i];
1574 }
1575 }
1576
1577 mid_eta /= denominator;
1578 mid_phi /= denominator;
1579 mid_x /= denominator;
1580 mid_y /= denominator;
1581
1582
1583
1584
1585 double See = 0.;
1586 double Spp = 0.;
1587 double Sep = 0.;
1588 double deta(0), dphi(0);
1589
1590 for (unsigned int i = 0; i < nCry; i++) {
1591 if (isBarrel) {
1592 deta = etaDetId[i] - mid_eta;
1593 dphi = phiDetId[i] - mid_phi;
1594 } else {
1595 deta = etaDetId[i] - mid_eta;
1596 float hitLocalR2 = (xDetId[i] - mid_x) * (xDetId[i] - mid_x) + (yDetId[i] - mid_y) * (yDetId[i] - mid_y);
1597 float hitR2 = xDetId[i] * xDetId[i] + yDetId[i] * yDetId[i];
1598 float meanR2 = mid_x * mid_x + mid_y * mid_y;
1599 float hitR = sqrt(hitR2);
1600 float meanR = sqrt(meanR2);
1601 float phi = acos((hitR2 + meanR2 - hitLocalR2) / (2 * hitR * meanR));
1602 dphi = hitR * phi;
1603 }
1604 See += (wiDetId[i] * deta * deta) / denominator;
1605 Spp += phiCorrectionFactor * (wiDetId[i] * dphi * dphi) / denominator;
1606 Sep += sqrt(phiCorrectionFactor) * (wiDetId[i] * deta * dphi) / denominator;
1607 }
1608
1609 Cluster2ndMoments returnMoments;
1610
1611
1612 returnMoments.sMaj = ((See + Spp) + sqrt((See - Spp) * (See - Spp) + 4. * Sep * Sep)) / 2.;
1613 returnMoments.sMin = ((See + Spp) - sqrt((See - Spp) * (See - Spp) + 4. * Sep * Sep)) / 2.;
1614
1615 returnMoments.alpha = atan((See - Spp + sqrt((Spp - See) * (Spp - See) + 4. * Sep * Sep)) / (2. * Sep));
1616
1617 return returnMoments;
1618 }
1619
1620
1621
1622
1623
1624
1625
1626
1627 template <bool noZS>
1628 std::vector<float> EcalClusterToolsT<noZS>::roundnessBarrelSuperClusters(
1629 const reco::SuperCluster &superCluster,
1630 const EcalRecHitCollection &recHits,
1631 int weightedPositionMethod,
1632 float energyThreshold) {
1633 std::vector<std::pair<const EcalRecHit *, float>> RH_ptrs_fracs;
1634 const std::vector<std::pair<DetId, float>> &myHitsPair = superCluster.hitsAndFractions();
1635 for (unsigned int i = 0; i < myHitsPair.size(); ++i) {
1636
1637 EcalRecHitCollection::const_iterator myRH = recHits.find(myHitsPair[i].first);
1638 if (myRH != recHits.end() && myRH->energy() * (noZS ? 1.0 : myHitsPair[i].second) > energyThreshold) {
1639
1640 RH_ptrs_fracs.push_back(std::make_pair(&(*myRH), myHitsPair[i].second));
1641 }
1642 }
1643 std::vector<float> temp =
1644 EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits(RH_ptrs_fracs, weightedPositionMethod);
1645 return temp;
1646 }
1647
1648
1649
1650
1651
1652 template <bool noZS>
1653 std::vector<float> EcalClusterToolsT<noZS>::roundnessBarrelSuperClustersUserExtended(
1654 const reco::SuperCluster &superCluster,
1655 const EcalRecHitCollection &recHits,
1656 int ieta_delta,
1657 int iphi_delta,
1658 float energyRHThresh,
1659 int weightedPositionMethod) {
1660 std::vector<std::pair<const EcalRecHit *, float>> RH_ptrs_fracs;
1661 const std::vector<std::pair<DetId, float>> &myHitsPair = superCluster.hitsAndFractions();
1662 for (unsigned int i = 0; i < myHitsPair.size(); ++i) {
1663
1664 EcalRecHitCollection::const_iterator myRH = recHits.find(myHitsPair[i].first);
1665 if (myRH != recHits.end() && myRH->energy() * (noZS ? 1.0 : myHitsPair[i].second) > energyRHThresh)
1666 RH_ptrs_fracs.push_back(std::make_pair(&(*myRH), myHitsPair[i].second));
1667 }
1668
1669 std::vector<int> seedPosition = EcalClusterToolsT<noZS>::getSeedPosition(RH_ptrs_fracs);
1670
1671 for (EcalRecHitCollection::const_iterator rh = recHits.begin(); rh != recHits.end(); rh++) {
1672 EBDetId EBdetIdi(rh->detid());
1673 float the_fraction = 0;
1674
1675 bool inEtaWindow = (abs(deltaIEta(seedPosition[0], EBdetIdi.ieta())) <= ieta_delta);
1676 bool inPhiWindow = (abs(deltaIPhi(seedPosition[1], EBdetIdi.iphi())) <= iphi_delta);
1677 bool passEThresh = (rh->energy() > energyRHThresh);
1678 bool alreadyCounted = false;
1679
1680
1681 bool is_SCrh_inside_recHits = false;
1682 for (unsigned int i = 0; i < myHitsPair.size(); i++) {
1683 EcalRecHitCollection::const_iterator SCrh = recHits.find(myHitsPair[i].first);
1684 if (SCrh != recHits.end()) {
1685 the_fraction = myHitsPair[i].second;
1686 is_SCrh_inside_recHits = true;
1687 if (rh->detid() == SCrh->detid())
1688 alreadyCounted = true;
1689 }
1690 }
1691
1692 if (is_SCrh_inside_recHits && !alreadyCounted && passEThresh && inEtaWindow && inPhiWindow) {
1693 RH_ptrs_fracs.push_back(std::make_pair(&(*rh), the_fraction));
1694 }
1695
1696 }
1697 return EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits(RH_ptrs_fracs, weightedPositionMethod);
1698 }
1699
1700
1701
1702
1703 template <bool noZS>
1704 std::vector<float> EcalClusterToolsT<noZS>::roundnessSelectedBarrelRecHits(
1705 const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs,
1706 int weightedPositionMethod) {
1707
1708
1709 std::vector<float> shapes;
1710
1711
1712 if (RH_ptrs_fracs.size() < 2) {
1713 shapes.push_back(-3);
1714 shapes.push_back(-3);
1715 return shapes;
1716 }
1717
1718
1719 std::vector<int> seedPosition = EcalClusterToolsT<noZS>::getSeedPosition(RH_ptrs_fracs);
1720 int tempInt = seedPosition[0];
1721 if (tempInt < 0)
1722 tempInt++;
1723 float energyTotal = EcalClusterToolsT<noZS>::getSumEnergy(RH_ptrs_fracs);
1724
1725
1726 float centerIEta = 0.;
1727 float centerIPhi = 0.;
1728 float denominator = 0.;
1729
1730 for (std::vector<std::pair<const EcalRecHit *, float>>::const_iterator rhf_ptr = RH_ptrs_fracs.begin();
1731 rhf_ptr != RH_ptrs_fracs.end();
1732 rhf_ptr++) {
1733 const EcalRecHit *rh_ptr = rhf_ptr->first;
1734
1735 EBDetId EBdetIdi(rh_ptr->detid());
1736 if (fabs(energyTotal) < 0.0001) {
1737
1738 shapes.push_back(-2);
1739 shapes.push_back(-2);
1740 return shapes;
1741 }
1742 float rh_energy = rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
1743 float weight = 0;
1744 if (std::abs(weightedPositionMethod) < 0.0001) {
1745 weight = rh_energy / energyTotal;
1746 } else {
1747 weight = std::max(0.0, 4.2 + log(rh_energy / energyTotal));
1748 }
1749 denominator += weight;
1750 centerIEta += weight * deltaIEta(seedPosition[0], EBdetIdi.ieta());
1751 centerIPhi += weight * deltaIPhi(seedPosition[1], EBdetIdi.iphi());
1752 }
1753 if (fabs(denominator) < 0.0001) {
1754
1755 shapes.push_back(-2);
1756 shapes.push_back(-2);
1757 return shapes;
1758 }
1759 centerIEta = centerIEta / denominator;
1760 centerIPhi = centerIPhi / denominator;
1761
1762
1763 TMatrixDSym inertia(2);
1764 double inertia00 = 0.;
1765 double inertia01 = 0.;
1766 double inertia11 = 0.;
1767 for (std::vector<std::pair<const EcalRecHit *, float>>::const_iterator rhf_ptr = RH_ptrs_fracs.begin();
1768 rhf_ptr != RH_ptrs_fracs.end();
1769 rhf_ptr++) {
1770 const EcalRecHit *rh_ptr = rhf_ptr->first;
1771
1772 EBDetId EBdetIdi(rh_ptr->detid());
1773
1774 if (fabs(energyTotal) < 0.0001) {
1775
1776 shapes.push_back(-2);
1777 shapes.push_back(-2);
1778 return shapes;
1779 }
1780 float rh_energy = rh_ptr->energy() * (noZS ? 1.0 : rhf_ptr->second);
1781 float weight = 0;
1782 if (std::abs(weightedPositionMethod) < 0.0001) {
1783 weight = rh_energy / energyTotal;
1784 } else {
1785 weight = std::max(0.0, 4.2 + log(rh_energy / energyTotal));
1786 }
1787
1788 float ieta_rh_to_center = deltaIEta(seedPosition[0], EBdetIdi.ieta()) - centerIEta;
1789 float iphi_rh_to_center = deltaIPhi(seedPosition[1], EBdetIdi.iphi()) - centerIPhi;
1790
1791 inertia00 += weight * iphi_rh_to_center * iphi_rh_to_center;
1792 inertia01 -= weight * iphi_rh_to_center * ieta_rh_to_center;
1793 inertia11 += weight * ieta_rh_to_center * ieta_rh_to_center;
1794 }
1795
1796 inertia[0][0] = inertia00;
1797 inertia[0][1] = inertia01;
1798 inertia[1][0] = inertia01;
1799 inertia[1][1] = inertia11;
1800
1801
1802 TMatrixD eVectors(2, 2);
1803 TVectorD eValues(2);
1804
1805 eVectors = inertia.EigenVectors(eValues);
1806
1807
1808
1809
1810 TVectorD smallerAxis(2);
1811 smallerAxis[0] = eVectors[0][1];
1812 smallerAxis[1] = eVectors[1][1];
1813
1814
1815 Double_t temp = fabs(smallerAxis[0]);
1816 if (fabs(eValues[0]) < 0.0001) {
1817
1818 shapes.push_back(-2);
1819 shapes.push_back(-2);
1820 return shapes;
1821 }
1822
1823 float Roundness = eValues[1] / eValues[0];
1824 float Angle = acos(temp);
1825
1826 if (-0.00001 < Roundness && Roundness < 0)
1827 Roundness = 0.;
1828 if (-0.00001 < Angle && Angle < 0)
1829 Angle = 0.;
1830
1831 shapes.push_back(Roundness);
1832 shapes.push_back(Angle);
1833 return shapes;
1834 }
1835
1836 template <bool noZS>
1837 int EcalClusterToolsT<noZS>::nrSaturatedCrysIn5x5(const DetId &id,
1838 const EcalRecHitCollection *recHits,
1839 const CaloTopology *topology) {
1840 int nrSat = 0;
1841 CaloRectangle rectangle{-2, 2, -2, 2};
1842 for (auto const &detId : rectangle(id, *topology)) {
1843 auto recHitIt = recHits->find(detId);
1844 if (recHitIt != recHits->end() && recHitIt->checkFlag(EcalRecHit::kSaturated)) {
1845 nrSat++;
1846 }
1847 }
1848 return nrSat;
1849 }
1850
1851
1852
1853
1854
1855 template <bool noZS>
1856 int EcalClusterToolsT<noZS>::deltaIPhi(int seed_iphi, int rh_iphi) {
1857 int rel_iphi = rh_iphi - seed_iphi;
1858
1859 if (rel_iphi > 180)
1860 rel_iphi = rel_iphi - 360;
1861 if (rel_iphi < -180)
1862 rel_iphi = rel_iphi + 360;
1863 return rel_iphi;
1864 }
1865
1866
1867
1868
1869 template <bool noZS>
1870 int EcalClusterToolsT<noZS>::deltaIEta(int seed_ieta, int rh_ieta) {
1871
1872 if (seed_ieta < 0)
1873 seed_ieta++;
1874 if (rh_ieta < 0)
1875 rh_ieta++;
1876 int rel_ieta = rh_ieta - seed_ieta;
1877 return rel_ieta;
1878 }
1879
1880
1881 template <bool noZS>
1882 std::vector<int> EcalClusterToolsT<noZS>::getSeedPosition(
1883 const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs) {
1884 std::vector<int> seedPosition;
1885 float eSeedRH = 0;
1886 int iEtaSeedRH = 0;
1887 int iPhiSeedRH = 0;
1888
1889 for (auto const &rhf : RH_ptrs_fracs) {
1890 const EcalRecHit *rh_ptr = rhf.first;
1891
1892 EBDetId EBdetIdi(rh_ptr->detid());
1893 float rh_energy = rh_ptr->energy() * (noZS ? 1.0 : rhf.second);
1894
1895 if (eSeedRH < rh_energy) {
1896 eSeedRH = rh_energy;
1897 iEtaSeedRH = EBdetIdi.ieta();
1898 iPhiSeedRH = EBdetIdi.iphi();
1899 }
1900
1901 }
1902
1903 seedPosition.push_back(iEtaSeedRH);
1904 seedPosition.push_back(iPhiSeedRH);
1905 return seedPosition;
1906 }
1907
1908
1909 template <bool noZS>
1910 float EcalClusterToolsT<noZS>::getSumEnergy(const std::vector<std::pair<const EcalRecHit *, float>> &RH_ptrs_fracs) {
1911 float sumE = 0.;
1912 for (const auto &hAndF : RH_ptrs_fracs) {
1913 sumE += hAndF.first->energy() * (noZS ? 1.0 : hAndF.second);
1914 }
1915 return sumE;
1916 }
1917
1918 typedef EcalClusterToolsT<false> EcalClusterTools;
1919
1920 namespace noZS {
1921 typedef EcalClusterToolsT<true> EcalClusterTools;
1922 }
1923
1924 #endif