File indexing completed on 2024-05-10 02:21:06
0001 #include <iostream>
0002
0003 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0004 #include "RecoEcal/EgammaCoreTools/interface/ClusterShapeAlgo.h"
0005 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
0006 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0007 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0008 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0009 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
0010
0011 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0012 #include "DataFormats/Math/interface/Point3D.h"
0013 #include "DataFormats/Math/interface/Vector3D.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015
0016 #include <CLHEP/Geometry/Transform3D.h>
0017 #include <CLHEP/Units/SystemOfUnits.h>
0018
0019 ClusterShapeAlgo::ClusterShapeAlgo(const edm::ParameterSet& par) : parameterSet_(par) {}
0020
0021 reco::ClusterShape ClusterShapeAlgo::Calculate(const reco::BasicCluster& passedCluster,
0022 const EcalRecHitCollection* hits,
0023 const CaloSubdetectorGeometry* geometry,
0024 const CaloSubdetectorTopology* topology) {
0025 Calculate_TopEnergy(passedCluster, hits);
0026 Calculate_2ndEnergy(passedCluster, hits);
0027 Create_Map(hits, topology);
0028 Calculate_e2x2();
0029 Calculate_e3x2();
0030 Calculate_e3x3();
0031 Calculate_e4x4();
0032 Calculate_e5x5();
0033 Calculate_e2x5Right();
0034 Calculate_e2x5Left();
0035 Calculate_e2x5Top();
0036 Calculate_e2x5Bottom();
0037 Calculate_Covariances(passedCluster, hits, geometry);
0038 Calculate_BarrelBasketEnergyFraction(passedCluster, hits, Eta, geometry);
0039 Calculate_BarrelBasketEnergyFraction(passedCluster, hits, Phi, geometry);
0040 Calculate_EnergyDepTopology(passedCluster, hits, geometry, true);
0041 Calculate_lat(passedCluster);
0042 Calculate_ComplexZernikeMoments(passedCluster);
0043
0044 return reco::ClusterShape(covEtaEta_,
0045 covEtaPhi_,
0046 covPhiPhi_,
0047 eMax_,
0048 eMaxId_,
0049 e2nd_,
0050 e2ndId_,
0051 e2x2_,
0052 e3x2_,
0053 e3x3_,
0054 e4x4_,
0055 e5x5_,
0056 e2x5Right_,
0057 e2x5Left_,
0058 e2x5Top_,
0059 e2x5Bottom_,
0060 e3x2Ratio_,
0061 lat_,
0062 etaLat_,
0063 phiLat_,
0064 A20_,
0065 A42_,
0066 energyBasketFractionEta_,
0067 energyBasketFractionPhi_);
0068 }
0069
0070 void ClusterShapeAlgo::Calculate_TopEnergy(const reco::BasicCluster& passedCluster, const EcalRecHitCollection* hits) {
0071 double eMax = 0;
0072 DetId eMaxId(0);
0073
0074 const std::vector<std::pair<DetId, float> >& clusterDetIds = passedCluster.hitsAndFractions();
0075
0076 EcalRecHit testEcalRecHit;
0077
0078 for (auto const& posCurrent : clusterDetIds) {
0079 if ((posCurrent.first != DetId(0)) && (hits->find(posCurrent.first) != hits->end())) {
0080 EcalRecHitCollection::const_iterator itt = hits->find(posCurrent.first);
0081 testEcalRecHit = *itt;
0082
0083 if (testEcalRecHit.energy() * posCurrent.second > eMax) {
0084 eMax = testEcalRecHit.energy() * posCurrent.second;
0085 eMaxId = testEcalRecHit.id();
0086 }
0087 }
0088 }
0089
0090 eMax_ = eMax;
0091 eMaxId_ = eMaxId;
0092 }
0093
0094 void ClusterShapeAlgo::Calculate_2ndEnergy(const reco::BasicCluster& passedCluster, const EcalRecHitCollection* hits) {
0095 double e2nd = 0;
0096 DetId e2ndId(0);
0097
0098 const std::vector<std::pair<DetId, float> >& clusterDetIds = passedCluster.hitsAndFractions();
0099
0100 EcalRecHit testEcalRecHit;
0101
0102 for (auto const& posCurrent : clusterDetIds) {
0103 if ((posCurrent.first != DetId(0)) && (hits->find(posCurrent.first) != hits->end())) {
0104 EcalRecHitCollection::const_iterator itt = hits->find(posCurrent.first);
0105 testEcalRecHit = *itt;
0106
0107 if (testEcalRecHit.energy() * posCurrent.second > e2nd && testEcalRecHit.id() != eMaxId_) {
0108 e2nd = testEcalRecHit.energy() * posCurrent.second;
0109 e2ndId = testEcalRecHit.id();
0110 }
0111 }
0112 }
0113
0114 e2nd_ = e2nd;
0115 e2ndId_ = e2ndId;
0116 }
0117
0118 void ClusterShapeAlgo::Create_Map(const EcalRecHitCollection* hits, const CaloSubdetectorTopology* topology) {
0119 EcalRecHit tempEcalRecHit;
0120 CaloNavigator<DetId> posCurrent = CaloNavigator<DetId>(eMaxId_, topology);
0121
0122 for (int x = 0; x < 5; x++)
0123 for (int y = 0; y < 5; y++) {
0124 posCurrent.home();
0125 posCurrent.offsetBy(-2 + x, -2 + y);
0126
0127 if ((*posCurrent != DetId(0)) && (hits->find(*posCurrent) != hits->end())) {
0128 EcalRecHitCollection::const_iterator itt = hits->find(*posCurrent);
0129 tempEcalRecHit = *itt;
0130 energyMap_[y][x] = std::make_pair(tempEcalRecHit.id(), tempEcalRecHit.energy());
0131 } else
0132 energyMap_[y][x] = std::make_pair(DetId(0), 0);
0133 }
0134 }
0135
0136 void ClusterShapeAlgo::Calculate_e2x2() {
0137 double e2x2Max = 0;
0138 double e2x2Test = 0;
0139
0140 int deltaX = 0, deltaY = 0;
0141
0142 for (int corner = 0; corner < 4; corner++) {
0143 switch (corner) {
0144 case 0:
0145 deltaX = -1;
0146 deltaY = -1;
0147 break;
0148 case 1:
0149 deltaX = -1;
0150 deltaY = 1;
0151 break;
0152 case 2:
0153 deltaX = 1;
0154 deltaY = -1;
0155 break;
0156 case 3:
0157 deltaX = 1;
0158 deltaY = 1;
0159 break;
0160 }
0161
0162 e2x2Test = energyMap_[2][2].second;
0163 e2x2Test += energyMap_[2 + deltaY][2].second;
0164 e2x2Test += energyMap_[2][2 + deltaX].second;
0165 e2x2Test += energyMap_[2 + deltaY][2 + deltaX].second;
0166
0167 if (e2x2Test > e2x2Max) {
0168 e2x2Max = e2x2Test;
0169 e2x2_Diagonal_X_ = 2 + deltaX;
0170 e2x2_Diagonal_Y_ = 2 + deltaY;
0171 }
0172 }
0173
0174 e2x2_ = e2x2Max;
0175 }
0176
0177 void ClusterShapeAlgo::Calculate_e3x2() {
0178 double e3x2 = 0.0;
0179 double e3x2Ratio = 0.0, e3x2RatioNumerator = 0.0, e3x2RatioDenominator = 0.0;
0180
0181 int e2ndX = 2, e2ndY = 2;
0182 int deltaY = 0, deltaX = 0;
0183
0184 double nextEnergy = -999;
0185 int nextEneryDirection = -1;
0186
0187 for (int cardinalDirection = 0; cardinalDirection < 4; cardinalDirection++) {
0188 switch (cardinalDirection) {
0189 case 0:
0190 deltaX = -1;
0191 deltaY = 0;
0192 break;
0193 case 1:
0194 deltaX = 1;
0195 deltaY = 0;
0196 break;
0197 case 2:
0198 deltaX = 0;
0199 deltaY = -1;
0200 break;
0201 case 3:
0202 deltaX = 0;
0203 deltaY = 1;
0204 break;
0205 }
0206
0207 if (energyMap_[2 + deltaY][2 + deltaX].second >= nextEnergy) {
0208 nextEnergy = energyMap_[2 + deltaY][2 + deltaX].second;
0209 nextEneryDirection = cardinalDirection;
0210
0211 e2ndX = 2 + deltaX;
0212 e2ndY = 2 + deltaY;
0213 }
0214 }
0215
0216 switch (nextEneryDirection) {
0217 case 0:;
0218 case 1:
0219 deltaX = 0;
0220 deltaY = 1;
0221 break;
0222 case 2:;
0223 case 3:
0224 deltaX = 1;
0225 deltaY = 0;
0226 break;
0227 }
0228
0229 for (int sign = -1; sign <= 1; sign++)
0230 e3x2 += (energyMap_[2 + deltaY * sign][2 + deltaX * sign].second +
0231 energyMap_[e2ndY + deltaY * sign][e2ndX + deltaX * sign].second);
0232
0233 e3x2RatioNumerator =
0234 energyMap_[e2ndY + deltaY][e2ndX + deltaX].second + energyMap_[e2ndY - deltaY][e2ndX - deltaX].second;
0235 e3x2RatioDenominator = 0.5 + energyMap_[2 + deltaY][2 + deltaX].second + energyMap_[2 - deltaY][2 - deltaX].second;
0236 e3x2Ratio = e3x2RatioNumerator / e3x2RatioDenominator;
0237
0238 e3x2_ = e3x2;
0239 e3x2Ratio_ = e3x2Ratio;
0240 }
0241
0242 void ClusterShapeAlgo::Calculate_e3x3() {
0243 double e3x3 = 0;
0244
0245 for (int i = 1; i <= 3; i++)
0246 for (int j = 1; j <= 3; j++)
0247 e3x3 += energyMap_[j][i].second;
0248
0249 e3x3_ = e3x3;
0250 }
0251
0252 void ClusterShapeAlgo::Calculate_e4x4() {
0253 double e4x4 = 0;
0254
0255 int startX = -1, startY = -1;
0256
0257 switch (e2x2_Diagonal_X_) {
0258 case 1:
0259 startX = 0;
0260 break;
0261 case 3:
0262 startX = 1;
0263 break;
0264 }
0265
0266 switch (e2x2_Diagonal_Y_) {
0267 case 1:
0268 startY = 0;
0269 break;
0270 case 3:
0271 startY = 1;
0272 break;
0273 }
0274
0275 for (int i = startX; i <= startX + 3; i++)
0276 for (int j = startY; j <= startY + 3; j++)
0277 e4x4 += energyMap_[j][i].second;
0278
0279 e4x4_ = e4x4;
0280 }
0281
0282 void ClusterShapeAlgo::Calculate_e5x5() {
0283 double e5x5 = 0;
0284
0285 for (int i = 0; i <= 4; i++)
0286 for (int j = 0; j <= 4; j++)
0287 e5x5 += energyMap_[j][i].second;
0288
0289 e5x5_ = e5x5;
0290 }
0291
0292 void ClusterShapeAlgo::Calculate_e2x5Right() {
0293 double e2x5R = 0.0;
0294 for (int i = 0; i <= 4; i++) {
0295 for (int j = 0; j <= 4; j++) {
0296 if (j > 2) {
0297 e2x5R += energyMap_[i][j].second;
0298 }
0299 }
0300 }
0301 e2x5Right_ = e2x5R;
0302 }
0303
0304 void ClusterShapeAlgo::Calculate_e2x5Left() {
0305 double e2x5L = 0.0;
0306 for (int i = 0; i <= 4; i++) {
0307 for (int j = 0; j <= 4; j++) {
0308 if (j < 2) {
0309 e2x5L += energyMap_[i][j].second;
0310 }
0311 }
0312 }
0313 e2x5Left_ = e2x5L;
0314 }
0315
0316 void ClusterShapeAlgo::Calculate_e2x5Bottom() {
0317 double e2x5B = 0.0;
0318 for (int i = 0; i <= 4; i++) {
0319 for (int j = 0; j <= 4; j++) {
0320 if (i > 2) {
0321 e2x5B += energyMap_[i][j].second;
0322 }
0323 }
0324 }
0325 e2x5Bottom_ = e2x5B;
0326 }
0327
0328 void ClusterShapeAlgo::Calculate_e2x5Top() {
0329 double e2x5T = 0.0;
0330 for (int i = 0; i <= 4; i++) {
0331 for (int j = 0; j <= 4; j++) {
0332 if (i < 2) {
0333 e2x5T += energyMap_[i][j].second;
0334 }
0335 }
0336 }
0337 e2x5Top_ = e2x5T;
0338 }
0339
0340 void ClusterShapeAlgo::Calculate_Covariances(const reco::BasicCluster& passedCluster,
0341 const EcalRecHitCollection* hits,
0342 const CaloSubdetectorGeometry* geometry) {
0343 if (e5x5_ > 0.) {
0344 double w0_ = parameterSet_.getParameter<double>("W0");
0345
0346
0347 math::XYZVector meanPosition(0.0, 0.0, 0.0);
0348 for (int i = 0; i < 5; ++i) {
0349 for (int j = 0; j < 5; ++j) {
0350 DetId id = energyMap_[i][j].first;
0351 if (id != DetId(0)) {
0352 GlobalPoint positionGP = geometry->getGeometry(id)->getPosition();
0353 math::XYZVector position(positionGP.x(), positionGP.y(), positionGP.z());
0354 meanPosition = meanPosition + energyMap_[i][j].second * position;
0355 }
0356 }
0357 }
0358
0359 meanPosition /= e5x5_;
0360
0361
0362 double numeratorEtaEta = 0;
0363 double numeratorEtaPhi = 0;
0364 double numeratorPhiPhi = 0;
0365 double denominator = 0;
0366
0367 for (int i = 0; i < 5; ++i) {
0368 for (int j = 0; j < 5; ++j) {
0369 DetId id = energyMap_[i][j].first;
0370 if (id != DetId(0)) {
0371 GlobalPoint position = geometry->getGeometry(id)->getPosition();
0372
0373 double dPhi = position.phi() - meanPosition.phi();
0374 if (dPhi > +Geom::pi()) {
0375 dPhi = Geom::twoPi() - dPhi;
0376 }
0377 if (dPhi < -Geom::pi()) {
0378 dPhi = Geom::twoPi() + dPhi;
0379 }
0380
0381 double dEta = position.eta() - meanPosition.eta();
0382 double w = 0.;
0383 if (energyMap_[i][j].second > 0.)
0384 w = std::max(0.0, w0_ + log(energyMap_[i][j].second / e5x5_));
0385
0386 denominator += w;
0387 numeratorEtaEta += w * dEta * dEta;
0388 numeratorEtaPhi += w * dEta * dPhi;
0389 numeratorPhiPhi += w * dPhi * dPhi;
0390 }
0391 }
0392 }
0393
0394 covEtaEta_ = numeratorEtaEta / denominator;
0395 covEtaPhi_ = numeratorEtaPhi / denominator;
0396 covPhiPhi_ = numeratorPhiPhi / denominator;
0397 } else {
0398
0399
0400 covEtaEta_ = 0;
0401 covEtaPhi_ = 0;
0402 covPhiPhi_ = 0;
0403 }
0404 }
0405
0406 void ClusterShapeAlgo::Calculate_BarrelBasketEnergyFraction(const reco::BasicCluster& passedCluster,
0407 const EcalRecHitCollection* hits,
0408 const int EtaPhi,
0409 const CaloSubdetectorGeometry* geometry) {
0410 if ((hits == nullptr) || (((*hits)[0]).id().subdetId() != EcalBarrel)) {
0411
0412 return;
0413 }
0414
0415 std::map<int, double> indexedBasketEnergy;
0416 const std::vector<std::pair<DetId, float> >& clusterDetIds = passedCluster.hitsAndFractions();
0417 const EcalBarrelGeometry* subDetGeometry = dynamic_cast<const EcalBarrelGeometry*>(geometry);
0418
0419 for (auto const& posCurrent : clusterDetIds) {
0420 int basketIndex = 999;
0421
0422 if (EtaPhi == Eta) {
0423 int unsignedIEta = abs(EBDetId(posCurrent.first).ieta());
0424 std::vector<int> etaBasketSize = subDetGeometry->getEtaBaskets();
0425
0426 for (unsigned int i = 0; i < etaBasketSize.size(); i++) {
0427 unsignedIEta -= etaBasketSize[i];
0428 if (unsignedIEta - 1 < 0) {
0429 basketIndex = i;
0430 break;
0431 }
0432 }
0433 basketIndex = (basketIndex + 1) * (EBDetId(posCurrent.first).ieta() > 0 ? 1 : -1);
0434
0435 } else if (EtaPhi == Phi) {
0436 int halfNumBasketInPhi = (EBDetId::MAX_IPHI - EBDetId::MIN_IPHI + 1) / subDetGeometry->getBasketSizeInPhi() / 2;
0437
0438 basketIndex = (EBDetId(posCurrent.first).iphi() - 1) / subDetGeometry->getBasketSizeInPhi() -
0439 (EBDetId((clusterDetIds[0]).first).iphi() - 1) / subDetGeometry->getBasketSizeInPhi();
0440
0441 if (basketIndex >= halfNumBasketInPhi)
0442 basketIndex -= 2 * halfNumBasketInPhi;
0443 else if (basketIndex < -1 * halfNumBasketInPhi)
0444 basketIndex += 2 * halfNumBasketInPhi;
0445
0446 } else
0447 throw(std::runtime_error("\n\nOh No! Calculate_BarrelBasketEnergyFraction called on invalid index.\n\n"));
0448
0449 indexedBasketEnergy[basketIndex] += (hits->find(posCurrent.first))->energy();
0450 }
0451
0452 std::vector<double> energyFraction;
0453 for (std::map<int, double>::iterator posCurrent = indexedBasketEnergy.begin();
0454 posCurrent != indexedBasketEnergy.end();
0455 posCurrent++) {
0456 energyFraction.push_back(indexedBasketEnergy[posCurrent->first] / passedCluster.energy());
0457 }
0458
0459 switch (EtaPhi) {
0460 case Eta:
0461 energyBasketFractionEta_ = energyFraction;
0462 break;
0463 case Phi:
0464 energyBasketFractionPhi_ = energyFraction;
0465 break;
0466 }
0467 }
0468
0469 void ClusterShapeAlgo::Calculate_lat(const reco::BasicCluster& passedCluster) {
0470 double r, redmoment = 0;
0471 double phiRedmoment = 0;
0472 double etaRedmoment = 0;
0473 int n, n1, n2, tmp;
0474 int clusterSize = energyDistribution_.size();
0475 if (clusterSize < 3) {
0476 etaLat_ = 0.0;
0477 lat_ = 0.0;
0478 return;
0479 }
0480
0481 n1 = 0;
0482 n2 = 1;
0483 if (energyDistribution_[1].deposited_energy > energyDistribution_[0].deposited_energy) {
0484 tmp = n2;
0485 n2 = n1;
0486 n1 = tmp;
0487 }
0488 for (int i = 2; i < clusterSize; i++) {
0489 n = i;
0490 if (energyDistribution_[i].deposited_energy > energyDistribution_[n1].deposited_energy) {
0491 tmp = n2;
0492 n2 = n1;
0493 n1 = i;
0494 n = tmp;
0495 } else {
0496 if (energyDistribution_[i].deposited_energy > energyDistribution_[n2].deposited_energy) {
0497 tmp = n2;
0498 n2 = i;
0499 n = tmp;
0500 }
0501 }
0502
0503 r = energyDistribution_[n].r;
0504 redmoment += r * r * energyDistribution_[n].deposited_energy;
0505 double rphi = r * cos(energyDistribution_[n].phi);
0506 phiRedmoment += rphi * rphi * energyDistribution_[n].deposited_energy;
0507 double reta = r * sin(energyDistribution_[n].phi);
0508 etaRedmoment += reta * reta * energyDistribution_[n].deposited_energy;
0509 }
0510 double e1 = energyDistribution_[n1].deposited_energy;
0511 double e2 = energyDistribution_[n2].deposited_energy;
0512
0513 lat_ = redmoment / (redmoment + 2.19 * 2.19 * (e1 + e2));
0514 phiLat_ = phiRedmoment / (phiRedmoment + 2.19 * 2.19 * (e1 + e2));
0515 etaLat_ = etaRedmoment / (etaRedmoment + 2.19 * 2.19 * (e1 + e2));
0516 }
0517
0518 void ClusterShapeAlgo::Calculate_ComplexZernikeMoments(const reco::BasicCluster& passedCluster) {
0519
0520
0521 A20_ = absZernikeMoment(passedCluster, 2, 0);
0522 A42_ = absZernikeMoment(passedCluster, 4, 2);
0523 }
0524
0525 double ClusterShapeAlgo::absZernikeMoment(const reco::BasicCluster& passedCluster, int n, int m, double R0) {
0526
0527 if ((m > n) || ((n - m) % 2 != 0) || (n < 0) || (m < 0))
0528 return -1;
0529
0530
0531
0532 if ((n > 20) || (R0 <= 2.19))
0533 return -1;
0534 if (n <= 5)
0535 return fast_AbsZernikeMoment(passedCluster, n, m, R0);
0536 else
0537 return calc_AbsZernikeMoment(passedCluster, n, m, R0);
0538 }
0539
0540 double ClusterShapeAlgo::f00(double r) { return 1; }
0541
0542 double ClusterShapeAlgo::f11(double r) { return r; }
0543
0544 double ClusterShapeAlgo::f20(double r) { return 2.0 * r * r - 1.0; }
0545
0546 double ClusterShapeAlgo::f22(double r) { return r * r; }
0547
0548 double ClusterShapeAlgo::f31(double r) { return 3.0 * r * r * r - 2.0 * r; }
0549
0550 double ClusterShapeAlgo::f33(double r) { return r * r * r; }
0551
0552 double ClusterShapeAlgo::f40(double r) { return 6.0 * r * r * r * r - 6.0 * r * r + 1.0; }
0553
0554 double ClusterShapeAlgo::f42(double r) { return 4.0 * r * r * r * r - 3.0 * r * r; }
0555
0556 double ClusterShapeAlgo::f44(double r) { return r * r * r * r; }
0557
0558 double ClusterShapeAlgo::f51(double r) { return 10.0 * pow(r, 5) - 12.0 * pow(r, 3) + 3.0 * r; }
0559
0560 double ClusterShapeAlgo::f53(double r) { return 5.0 * pow(r, 5) - 4.0 * pow(r, 3); }
0561
0562 double ClusterShapeAlgo::f55(double r) { return pow(r, 5); }
0563
0564 double ClusterShapeAlgo::fast_AbsZernikeMoment(const reco::BasicCluster& passedCluster, int n, int m, double R0) {
0565 double r, ph, e, Re = 0, Im = 0, result;
0566 double TotalEnergy = passedCluster.energy();
0567 int index = (n / 2) * (n / 2) + (n / 2) + m;
0568 int clusterSize = energyDistribution_.size();
0569 if (clusterSize < 3)
0570 return 0.0;
0571
0572 for (int i = 0; i < clusterSize; i++) {
0573 r = energyDistribution_[i].r / R0;
0574 if (r < 1) {
0575 fcn_.clear();
0576 Calculate_Polynomials(r);
0577 ph = (energyDistribution_[i]).phi;
0578 e = energyDistribution_[i].deposited_energy;
0579 Re = Re + e / TotalEnergy * fcn_[index] * cos((double)m * ph);
0580 Im = Im - e / TotalEnergy * fcn_[index] * sin((double)m * ph);
0581 }
0582 }
0583 result = sqrt(Re * Re + Im * Im);
0584
0585 return result;
0586 }
0587
0588 double ClusterShapeAlgo::calc_AbsZernikeMoment(const reco::BasicCluster& passedCluster, int n, int m, double R0) {
0589 double r, ph, e, Re = 0, Im = 0, f_nm, result;
0590 double TotalEnergy = passedCluster.energy();
0591 int clusterSize = energyDistribution_.size();
0592 if (clusterSize < 3)
0593 return 0.0;
0594
0595 for (int i = 0; i < clusterSize; i++) {
0596 r = energyDistribution_[i].r / R0;
0597 if (r < 1) {
0598 ph = (energyDistribution_[i]).phi;
0599 e = energyDistribution_[i].deposited_energy;
0600 f_nm = 0;
0601 for (int s = 0; s <= (n - m) / 2; s++) {
0602 if (s % 2 == 0) {
0603 f_nm = f_nm + factorial(n - s) * pow(r, (double)(n - 2 * s)) /
0604 (factorial(s) * factorial((n + m) / 2 - s) * factorial((n - m) / 2 - s));
0605 } else {
0606 f_nm = f_nm - factorial(n - s) * pow(r, (double)(n - 2 * s)) /
0607 (factorial(s) * factorial((n + m) / 2 - s) * factorial((n - m) / 2 - s));
0608 }
0609 }
0610 Re = Re + e / TotalEnergy * f_nm * cos((double)m * ph);
0611 Im = Im - e / TotalEnergy * f_nm * sin((double)m * ph);
0612 }
0613 }
0614 result = sqrt(Re * Re + Im * Im);
0615
0616 return result;
0617 }
0618
0619 void ClusterShapeAlgo::Calculate_EnergyDepTopology(const reco::BasicCluster& passedCluster,
0620 const EcalRecHitCollection* hits,
0621 const CaloSubdetectorGeometry* geometry,
0622 bool logW) {
0623
0624 energyDistribution_.clear();
0625
0626
0627
0628 CLHEP::Hep3Vector clVect(passedCluster.position().x(), passedCluster.position().y(), passedCluster.position().z());
0629 CLHEP::Hep3Vector clDir(clVect);
0630 clDir *= 1.0 / clDir.mag();
0631
0632 CLHEP::Hep3Vector theta_axis(clDir.y(), -clDir.x(), 0.0);
0633 theta_axis *= 1.0 / theta_axis.mag();
0634 CLHEP::Hep3Vector phi_axis = theta_axis.cross(clDir);
0635
0636 const std::vector<std::pair<DetId, float> >& clusterDetIds = passedCluster.hitsAndFractions();
0637
0638 EcalClusterEnergyDeposition clEdep;
0639 EcalRecHit testEcalRecHit;
0640
0641 for (auto const& posCurrent : clusterDetIds) {
0642 EcalRecHitCollection::const_iterator itt = hits->find(posCurrent.first);
0643 testEcalRecHit = *itt;
0644
0645 if ((posCurrent.first != DetId(0)) && (hits->find(posCurrent.first) != hits->end())) {
0646 clEdep.deposited_energy = testEcalRecHit.energy();
0647
0648
0649 if (logW) {
0650 double w0_ = parameterSet_.getParameter<double>("W0");
0651
0652 if (clEdep.deposited_energy == 0) {
0653 LogDebug("ClusterShapeAlgo") << "Crystal has zero energy; skipping... ";
0654 continue;
0655 }
0656 double weight = std::max(0.0, w0_ + log(fabs(clEdep.deposited_energy) / passedCluster.energy()));
0657 if (weight == 0) {
0658 LogDebug("ClusterShapeAlgo") << "Crystal has insufficient energy: E = " << clEdep.deposited_energy
0659 << " GeV; skipping... ";
0660 continue;
0661 } else
0662 LogDebug("ClusterShapeAlgo") << "===> got crystal. Energy = " << clEdep.deposited_energy << " GeV. ";
0663 }
0664 DetId id_ = posCurrent.first;
0665 auto this_cell = geometry->getGeometry(id_);
0666 const GlobalPoint& cellPos = this_cell->getPosition();
0667 CLHEP::Hep3Vector gblPos(cellPos.x(), cellPos.y(), cellPos.z());
0668
0669 CLHEP::Hep3Vector diff = gblPos - clVect;
0670
0671
0672
0673 CLHEP::Hep3Vector DigiVect = diff - diff.dot(clDir) * clDir;
0674 clEdep.r = DigiVect.mag();
0675 LogDebug("ClusterShapeAlgo") << "E = " << clEdep.deposited_energy << "\tdiff = " << diff.mag()
0676 << "\tr = " << clEdep.r;
0677 clEdep.phi = DigiVect.angle(theta_axis);
0678 if (DigiVect.dot(phi_axis) < 0)
0679 clEdep.phi = 2 * M_PI - clEdep.phi;
0680 energyDistribution_.push_back(clEdep);
0681 }
0682 }
0683 }
0684
0685 void ClusterShapeAlgo::Calculate_Polynomials(double rho) {
0686 fcn_.push_back(f00(rho));
0687 fcn_.push_back(f11(rho));
0688 fcn_.push_back(f20(rho));
0689 fcn_.push_back(f31(rho));
0690 fcn_.push_back(f22(rho));
0691 fcn_.push_back(f33(rho));
0692 fcn_.push_back(f40(rho));
0693 fcn_.push_back(f51(rho));
0694 fcn_.push_back(f42(rho));
0695 fcn_.push_back(f53(rho));
0696 fcn_.push_back(f44(rho));
0697 fcn_.push_back(f55(rho));
0698 }
0699
0700 double ClusterShapeAlgo::factorial(int n) const {
0701 double res = 1.0;
0702 for (int i = 2; i <= n; i++)
0703 res *= (double)i;
0704 return res;
0705 }