EgammaTools/src/HGCalShowerShapeHelper.cc

0001 #include "RecoEgamma/EgammaTools/interface/HGCalShowerShapeHelper.h"
0002
0003 #include <Math/Vector3D.h>
0004 #include <TMatrixD.h>
0005 #include <TVectorD.h>
0006
0007 #include <algorithm>
0008
0009 const double HGCalShowerShapeHelper::kLDWaferCellSize_ = 0.698;
0010 const double HGCalShowerShapeHelper::kHDWaferCellSize_ = 0.465;
0011
0012 HGCalShowerShapeHelper::ShowerShapeCalc::ShowerShapeCalc(
0013     std::shared_ptr<const hgcal::RecHitTools> recHitTools,
0014     std::shared_ptr<const std::unordered_map<uint32_t, const reco::PFRecHit *>> pfRecHitPtrMap,
0015     const std::vector<std::pair<DetId, float>> &hitsAndFracs,
0016     const double rawEnergy,
0017     const double minHitE,
0018     const double minHitET,
0019     const int minLayer,
0020     const int maxLayer,
0021     const DetId::Detector subDet)
0022     : recHitTools_(recHitTools),
0023       pfRecHitPtrMap_(pfRecHitPtrMap),
0024       rawEnergy_(rawEnergy),
0025       minHitE_(minHitE),
0026       minHitET_(minHitET),
0027       minHitET2_(minHitET * minHitET),
0028       minLayer_(minLayer),
0029       maxLayer_(maxLayer <= 0 ? recHitTools_->lastLayerEE() : maxLayer),
0030       nLayer_(maxLayer_ - minLayer_ + 1),
0031       subDet_(subDet) {
0032   assert(nLayer_ > 0);
0033   setFilteredHitsAndFractions(hitsAndFracs);
0034   setLayerWiseInfo();
0035 }
0036
0037 double HGCalShowerShapeHelper::ShowerShapeCalc::getCellSize(DetId detId) const {
0038   return recHitTools_->getSiThickIndex(detId) == 0 ? kHDWaferCellSize_ : kLDWaferCellSize_;
0039 }
0040
0041 double HGCalShowerShapeHelper::ShowerShapeCalc::getRvar(double cylinderR, bool useFractions, bool useCellSize) const {
0042   if (hitsAndFracs_.empty()) {
0043     return 0.0;
0044   }
0045
0046   if (rawEnergy_ <= 0.0) {
0047     edm::LogWarning("HGCalShowerShapeHelper")
0048         << "Encountered negative or zero energy for HGCal R-variable denominator: " << rawEnergy_ << std::endl;
0049   }
0050
0051   double cylinderR2 = cylinderR * cylinderR;
0052
0053   double rVar = 0.0;
0054
0055   auto hitEnergyIter = useFractions ? hitEnergiesWithFracs_.begin() : hitEnergies_.begin();
0056
0057   hitEnergyIter--;
0058
0059   for (const auto &hnf : hitsAndFracs_) {
0060     hitEnergyIter++;
0061
0062     DetId hitId = hnf.first;
0063
0064     int hitLayer = recHitTools_->getLayer(hitId) - 1;
0065
0066     const auto &hitPos = recHitTools_->getPosition(hitId);
0067     ROOT::Math::XYZVector hitXYZ(hitPos.x(), hitPos.y(), hitPos.z());
0068
0069     auto distXYZ = hitXYZ - layerCentroids_[hitLayer];
0070
0071     double r2 = distXYZ.x() * distXYZ.x() + distXYZ.y() * distXYZ.y();
0072
0073     // Including the cell size seems to make the variable less sensitive to the HD/LD transition region
0074     if (useCellSize) {
0075       if (std::sqrt(r2) > cylinderR + getCellSize(hitId)) {
0076         continue;
0077       }
0078     }
0079
0080     else if (r2 > cylinderR2) {
0081       continue;
0082     }
0083
0084     rVar += *hitEnergyIter;
0085   }
0086
0087   rVar /= rawEnergy_;
0088
0089   return rVar;
0090 }
0091
0092 HGCalShowerShapeHelper::ShowerWidths HGCalShowerShapeHelper::ShowerShapeCalc::getPCAWidths(double cylinderR,
0093                                                                                            bool useFractions) const {
0094   if (hitsAndFracs_.empty()) {
0095     return ShowerWidths();
0096   }
0097
0098   double cylinderR2 = cylinderR * cylinderR;
0099
0100   TMatrixD covMat(3, 3);
0101
0102   double dxdx = 0.0;
0103   double dydy = 0.0;
0104   double dzdz = 0.0;
0105
0106   double dxdy = 0.0;
0107   double dydz = 0.0;
0108   double dzdx = 0.0;
0109
0110   double totalW = 0.0;
0111
0112   auto hitEnergyIter = useFractions ? hitEnergiesWithFracs_.begin() : hitEnergies_.begin();
0113
0114   int nHit = 0;
0115   hitEnergyIter--;
0116
0117   for (const auto &hnf : hitsAndFracs_) {
0118     hitEnergyIter++;
0119
0120     DetId hitId = hnf.first;
0121
0122     const auto &hitPos = recHitTools_->getPosition(hitId);
0123     ROOT::Math::XYZVector hitXYZ(hitPos.x(), hitPos.y(), hitPos.z());
0124
0125     int hitLayer = recHitTools_->getLayer(hitId) - 1;
0126
0127     ROOT::Math::XYZVector radXYZ = hitXYZ - layerCentroids_[hitLayer];
0128
0129     double r2 = radXYZ.x() * radXYZ.x() + radXYZ.y() * radXYZ.y();
0130
0131     if (r2 > cylinderR2) {
0132       continue;
0133     }
0134
0135     ROOT::Math::XYZVector dXYZ = hitXYZ - centroid_;
0136
0137     double weight = *hitEnergyIter;
0138     totalW += weight;
0139
0140     dxdx += weight * dXYZ.x() * dXYZ.x();
0141     dydy += weight * dXYZ.y() * dXYZ.y();
0142     dzdz += weight * dXYZ.z() * dXYZ.z();
0143
0144     dxdy += weight * dXYZ.x() * dXYZ.y();
0145     dydz += weight * dXYZ.y() * dXYZ.z();
0146     dzdx += weight * dXYZ.z() * dXYZ.x();
0147
0148     nHit++;
0149   }
0150
0151   if (!totalW || nHit < 2) {
0152     return ShowerWidths();
0153   }
0154
0155   dxdx /= totalW;
0156   dydy /= totalW;
0157   dzdz /= totalW;
0158
0159   dxdy /= totalW;
0160   dydz /= totalW;
0161   dzdx /= totalW;
0162
0163   covMat(0, 0) = dxdx;
0164   covMat(1, 1) = dydy;
0165   covMat(2, 2) = dzdz;
0166
0167   covMat(0, 1) = covMat(1, 0) = dxdy;
0168   covMat(0, 2) = covMat(2, 0) = dzdx;
0169   covMat(1, 2) = covMat(2, 1) = dydz;
0170
0171   if (!covMat.Sum()) {
0172     return ShowerWidths();
0173   }
0174
0175   // Get eigen values and vectors
0176   TVectorD eigVals(3);
0177   TMatrixD eigVecMat(3, 3);
0178
0179   eigVecMat = covMat.EigenVectors(eigVals);
0180
0181   ShowerWidths returnWidths;
0182
0183   returnWidths.sigma2xx = dxdx;
0184   returnWidths.sigma2yy = dydy;
0185   returnWidths.sigma2zz = dzdz;
0186
0187   returnWidths.sigma2xy = dxdy;
0188   returnWidths.sigma2yz = dydz;
0189   returnWidths.sigma2zx = dzdx;
0190
0191   returnWidths.sigma2uu = eigVals(1);
0192   returnWidths.sigma2vv = eigVals(2);
0193   returnWidths.sigma2ww = eigVals(0);
0194
0195   return returnWidths;
0196 }
0197
0198 std::vector<double> HGCalShowerShapeHelper::ShowerShapeCalc::getEnergyHighestHits(unsigned int nrHits,
0199                                                                                   bool useFractions) const {
0200   std::vector<double> sortedEnergies(nrHits, 0.);
0201   const auto &hits = useFractions ? hitEnergiesWithFracs_ : hitEnergies_;
0202   std::partial_sort_copy(
0203       hits.begin(), hits.end(), sortedEnergies.begin(), sortedEnergies.end(), std::greater<double>());
0204   return sortedEnergies;
0205 }
0206
0207 void HGCalShowerShapeHelper::ShowerShapeCalc::setFilteredHitsAndFractions(
0208     const std::vector<std::pair<DetId, float>> &hitsAndFracs) {
0209   hitsAndFracs_.clear();
0210   hitEnergies_.clear();
0211   hitEnergiesWithFracs_.clear();
0212
0213   for (const auto &hnf : hitsAndFracs) {
0214     DetId hitId = hnf.first;
0215     float hitEfrac = hnf.second;
0216
0217     int hitLayer = recHitTools_->getLayer(hitId);
0218
0219     if (hitLayer > nLayer_) {
0220       continue;
0221     }
0222
0223     if (hitId.det() != subDet_) {
0224       continue;
0225     }
0226     auto hitIt = pfRecHitPtrMap_->find(hitId.rawId());
0227     if (hitIt == pfRecHitPtrMap_->end()) {
0228       continue;
0229     }
0230
0231     const reco::PFRecHit &recHit = *hitIt->second;
0232
0233     if (recHit.energy() < minHitE_) {
0234       continue;
0235     }
0236
0237     if (recHit.pt2() < minHitET2_) {
0238       continue;
0239     }
0240
0241     // Fill the vectors
0242     hitsAndFracs_.push_back(hnf);
0243     hitEnergies_.push_back(recHit.energy());
0244     hitEnergiesWithFracs_.push_back(recHit.energy() * hitEfrac);
0245   }
0246 }
0247
0248 void HGCalShowerShapeHelper::ShowerShapeCalc::setLayerWiseInfo() {
0249   layerEnergies_.clear();
0250   layerEnergies_.resize(nLayer_);
0251
0252   layerCentroids_.clear();
0253   layerCentroids_.resize(nLayer_);
0254
0255   centroid_.SetXYZ(0, 0, 0);
0256
0257   int iHit = -1;
0258   double totalW = 0.0;
0259
0260   // Compute the centroid per layer
0261   for (const auto &hnf : hitsAndFracs_) {
0262     iHit++;
0263
0264     DetId hitId = hnf.first;
0265
0266     double weight = hitEnergies_[iHit];
0267     totalW += weight;
0268
0269     const auto &hitPos = recHitTools_->getPosition(hitId);
0270     ROOT::Math::XYZVector hitXYZ(hitPos.x(), hitPos.y(), hitPos.z());
0271
0272     centroid_ += weight * hitXYZ;
0273
0274     int hitLayer = recHitTools_->getLayer(hitId) - 1;
0275
0276     layerEnergies_[hitLayer] += weight;
0277     layerCentroids_[hitLayer] += weight * hitXYZ;
0278   }
0279
0280   int iLayer = -1;
0281
0282   for (auto &centroid : layerCentroids_) {
0283     iLayer++;
0284
0285     if (layerEnergies_[iLayer]) {
0286       centroid /= layerEnergies_[iLayer];
0287     }
0288   }
0289
0290   if (totalW) {
0291     centroid_ /= totalW;
0292   }
0293 }
0294
0295 HGCalShowerShapeHelper::HGCalShowerShapeHelper()
0296     : recHitTools_(std::make_shared<hgcal::RecHitTools>()),
0297       pfRecHitPtrMap_(std::make_shared<std::unordered_map<uint32_t, const reco::PFRecHit *>>()) {}
0298
0299 HGCalShowerShapeHelper::HGCalShowerShapeHelper(edm::ConsumesCollector &&sumes)
0300     : recHitTools_(std::make_shared<hgcal::RecHitTools>()),
0301       pfRecHitPtrMap_(std::make_shared<std::unordered_map<uint32_t, const reco::PFRecHit *>>()) {
0302   setTokens(sumes);
0303 }
0304
0305 void HGCalShowerShapeHelper::initPerSetup(const edm::EventSetup &iSetup) {
0306   recHitTools_->setGeometry(iSetup.getData(caloGeometryToken_));
0307 }
0308
0309 void HGCalShowerShapeHelper::initPerEvent(const std::vector<reco::PFRecHit> &pfRecHits) {
0310   setPFRecHitPtrMap(pfRecHits);
0311 }
0312
0313 void HGCalShowerShapeHelper::initPerEvent(const edm::EventSetup &iSetup, const std::vector<reco::PFRecHit> &pfRecHits) {
0314   initPerSetup(iSetup);
0315   initPerEvent(pfRecHits);
0316 }
0317
0318 HGCalShowerShapeHelper::ShowerShapeCalc HGCalShowerShapeHelper::createCalc(
0319     const std::vector<std::pair<DetId, float>> &hitsAndFracs,
0320     double rawEnergy,
0321     double minHitE,
0322     double minHitET,
0323     int minLayer,
0324     int maxLayer,
0325     DetId::Detector subDet) const {
0326   return ShowerShapeCalc(
0327       recHitTools_, pfRecHitPtrMap_, hitsAndFracs, rawEnergy, minHitE, minHitET, minLayer, maxLayer, subDet);
0328 }
0329
0330 void HGCalShowerShapeHelper::setPFRecHitPtrMap(const std::vector<reco::PFRecHit> &recHits) {
0331   pfRecHitPtrMap_->clear();
0332
0333   for (const auto &recHit : recHits) {
0334     (*pfRecHitPtrMap_)[recHit.detId()] = &recHit;
0335   }
0336 }