EgammaCoreTools/src/ClusterShapeAlgo.cc

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     // first find energy-weighted mean position - doing it when filling the energy map might save time
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     // now we can calculate the covariances
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     // Warn the user if there was no energy in the cells and return zeroes.
0399     //       std::cout << "\ClusterShapeAlgo::Calculate_Covariances:  no energy in supplied cells.\n";
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     //std::cout << "No basket correction if no hits or for endacap!" << std::endl;
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   // Calculate only the moments which go into the default cluster shape
0520   // (moments with m>=2 are the only sensitive to azimuthal shape)
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   // 1. Check if n,m are correctly
0527   if ((m > n) || ((n - m) % 2 != 0) || (n < 0) || (m < 0))
0528     return -1;
0529
0530   // 2. Check if n,R0 are within validity Range :
0531   // n>20 or R0<2.19cm  just makes no sense !
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   // resets the energy distribution
0624   energyDistribution_.clear();
0625
0626   // init a map of the energy deposition centered on the
0627   // cluster centroid. This is for momenta calculation only.
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   // in the transverse plane, axis perpendicular to clusterDir
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   // loop over crystals
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       // if logarithmic weight is requested, apply cut on minimum energy of the recHit
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());  //surface position?
0668       // Evaluate the distance from the cluster centroid
0669       CLHEP::Hep3Vector diff = gblPos - clVect;
0670       // Important: for the moment calculation, only the "lateral distance" is important
0671       // "lateral distance" r_i = distance of the digi position from the axis Origin-Cluster Center
0672       // ---> subtract the projection on clDir
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 }