File indexing completed on 2023-03-17 11:17:55
0001
0002 #include "RecoEgamma/EgammaTools/interface/EgammaPCAHelper.h"
0003
0004 #include "DataFormats/HGCRecHit/interface/HGCRecHit.h"
0005 #include "DataFormats/Math/interface/deltaPhi.h"
0006 #include "DataFormats/ForwardDetId/interface/HGCalDetId.h"
0007 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0008
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010
0011 #include <algorithm>
0012 #include <iostream>
0013 #include <memory>
0014
0015 using namespace hgcal;
0016
0017 EGammaPCAHelper::EGammaPCAHelper()
0018 :
0019
0020
0021 invThicknessCorrection_({1. / 1.132, 1. / 1.092, 1. / 1.084}),
0022 pca_(new TPrincipal(3, "D")) {
0023 hitMap_ = nullptr;
0024 debug_ = false;
0025 }
0026
0027 void EGammaPCAHelper::setHitMap(const std::unordered_map<DetId, const HGCRecHit*>* hitMap) {
0028 hitMap_ = hitMap;
0029 pcaIteration_ = 0;
0030 }
0031
0032 void EGammaPCAHelper::setRecHitTools(const hgcal::RecHitTools* recHitTools) {
0033 recHitTools_ = recHitTools;
0034 maxlayer_ = recHitTools_->lastLayerBH();
0035 }
0036
0037 void EGammaPCAHelper::storeRecHits(const reco::HGCalMultiCluster& cluster) {
0038 theCluster_ = &cluster;
0039 std::vector<std::pair<DetId, float>> result;
0040 for (reco::HGCalMultiCluster::component_iterator it = cluster.begin(); it != cluster.end(); it++) {
0041 const std::vector<std::pair<DetId, float>>& hf = (*it)->hitsAndFractions();
0042 result.insert(result.end(), hf.begin(), hf.end());
0043 }
0044 storeRecHits(result);
0045 }
0046
0047 void EGammaPCAHelper::storeRecHits(const reco::CaloCluster& cluster) {
0048 theCluster_ = &cluster;
0049 storeRecHits(cluster.hitsAndFractions());
0050 }
0051
0052 void EGammaPCAHelper::storeRecHits(const std::vector<std::pair<DetId, float>>& hf) {
0053 std::vector<double> pcavars;
0054 pcavars.resize(3, 0.);
0055 theSpots_.clear();
0056 pcaIteration_ = 0;
0057
0058 sigu_ = 0.;
0059 sigv_ = 0.;
0060 sigp_ = 0.;
0061 sige_ = 0.;
0062
0063 unsigned hfsize = hf.size();
0064 if (debug_)
0065 std::cout << "The seed cluster constains " << hfsize << " hits " << std::endl;
0066
0067 if (hfsize == 0)
0068 return;
0069
0070 for (unsigned int j = 0; j < hfsize; j++) {
0071 unsigned int layer = recHitTools_->getLayerWithOffset(hf[j].first);
0072
0073 const DetId rh_detid = hf[j].first;
0074 std::unordered_map<DetId, const HGCRecHit*>::const_iterator itcheck = hitMap_->find(rh_detid);
0075 if (itcheck == hitMap_->end()) {
0076 edm::LogWarning("EgammaPCAHelper") << " Big problem, unable to find a hit " << rh_detid.rawId() << " "
0077 << rh_detid.det() << " " << HGCalDetId(rh_detid) << std::endl;
0078 continue;
0079 }
0080 if (debug_) {
0081 std::cout << "DetId " << rh_detid.rawId() << " " << layer << " " << itcheck->second->energy() << std::endl;
0082 std::cout << " Hit " << itcheck->second << " " << itcheck->second->energy() << std::endl;
0083 }
0084 float fraction = hf[j].second;
0085
0086 int thickIndex = recHitTools_->getSiThickIndex(rh_detid);
0087 double mip = dEdXWeights_[layer] * 0.001;
0088 if (thickIndex > -1 and thickIndex < 3)
0089 mip *= invThicknessCorrection_[thickIndex];
0090
0091 pcavars[0] = recHitTools_->getPosition(rh_detid).x();
0092 pcavars[1] = recHitTools_->getPosition(rh_detid).y();
0093 pcavars[2] = recHitTools_->getPosition(rh_detid).z();
0094 if (pcavars[2] == 0.)
0095 edm::LogWarning("EgammaPCAHelper") << " Problem, hit with z =0 ";
0096 else {
0097 Spot mySpot(rh_detid, itcheck->second->energy(), pcavars, layer, fraction, mip);
0098 theSpots_.push_back(mySpot);
0099 }
0100 }
0101 if (debug_) {
0102 std::cout << " Stored " << theSpots_.size() << " hits " << std::endl;
0103 }
0104 }
0105
0106 void EGammaPCAHelper::computePCA(float radius, bool withHalo) {
0107
0108 pca_ = std::make_unique<TPrincipal>(3, "D");
0109 bool initialCalculation = radius < 0;
0110 if (debug_)
0111 std::cout << " Initial calculation " << initialCalculation << std::endl;
0112 if (initialCalculation && withHalo) {
0113 edm::LogWarning("EGammaPCAHelper") << "Warning - in the first iteration, the halo hits are excluded " << std::endl;
0114 withHalo = false;
0115 }
0116
0117 float radius2 = radius * radius;
0118 if (!initialCalculation) {
0119 math::XYZVector mainAxis(axis_);
0120 mainAxis.unit();
0121 math::XYZVector phiAxis(barycenter_.x(), barycenter_.y(), 0);
0122 math::XYZVector udir(mainAxis.Cross(phiAxis));
0123 udir = udir.unit();
0124 trans_ = Transform3D(Point(barycenter_),
0125 Point(barycenter_ + axis_),
0126 Point(barycenter_ + udir),
0127 Point(0, 0, 0),
0128 Point(0., 0., 1.),
0129 Point(1., 0., 0.));
0130 }
0131
0132 std::set<int> layers;
0133 for (const auto& spot : theSpots_) {
0134 if (spot.layer() > recHitTools_->lastLayerEE())
0135 continue;
0136 if (!withHalo && (!spot.isCore()))
0137 continue;
0138 if (initialCalculation) {
0139
0140 if (!spot.isCore())
0141 continue;
0142 layers.insert(spot.layer());
0143 for (int i = 0; i < spot.multiplicity(); ++i)
0144 pca_->AddRow(spot.row());
0145 } else {
0146
0147 math::XYZPoint local = trans_(Point(spot.row()[0], spot.row()[1], spot.row()[2]));
0148 if (local.Perp2() > radius2)
0149 continue;
0150 layers.insert(spot.layer());
0151 for (int i = 0; i < spot.multiplicity(); ++i)
0152 pca_->AddRow(spot.row());
0153 }
0154 }
0155 if (debug_)
0156 std::cout << " Nlayers " << layers.size() << std::endl;
0157 if (layers.size() < 3) {
0158 pcaIteration_ = -1;
0159 return;
0160 }
0161 pca_->MakePrincipals();
0162 ++pcaIteration_;
0163 const TVectorD& means = *(pca_->GetMeanValues());
0164 const TMatrixD& eigens = *(pca_->GetEigenVectors());
0165
0166 barycenter_ = math::XYZPoint(means[0], means[1], means[2]);
0167 axis_ = math::XYZVector(eigens(0, 0), eigens(1, 0), eigens(2, 0));
0168 if (axis_.z() * barycenter_.z() < 0.0) {
0169 axis_ = -1. * axis_;
0170 }
0171 }
0172
0173 void EGammaPCAHelper::computeShowerWidth(float radius, bool withHalo) {
0174 sigu_ = 0.;
0175 sigv_ = 0.;
0176 sigp_ = 0.;
0177 sige_ = 0.;
0178 double cyl_ene = 0.;
0179
0180 float radius2 = radius * radius;
0181 for (const auto& spot : theSpots_) {
0182 Point globalPoint(spot.row()[0], spot.row()[1], spot.row()[2]);
0183 math::XYZPoint local = trans_(globalPoint);
0184 if (local.Perp2() > radius2)
0185 continue;
0186
0187
0188 if (withHalo && spot.fraction() < 0)
0189 continue;
0190 if (!withHalo && !(spot.isCore()))
0191 continue;
0192
0193 sige_ += (globalPoint.eta() - theCluster_->eta()) * (globalPoint.eta() - theCluster_->eta()) * spot.energy();
0194 sigp_ += deltaPhi(globalPoint.phi(), theCluster_->phi()) * deltaPhi(globalPoint.phi(), theCluster_->phi()) *
0195 spot.energy();
0196
0197 sigu_ += local.x() * local.x() * spot.energy();
0198 sigv_ += local.y() * local.y() * spot.energy();
0199 cyl_ene += spot.energy();
0200 }
0201
0202 if (cyl_ene > 0.) {
0203 const double inv_cyl_ene = 1. / cyl_ene;
0204 sigu_ = sigu_ * inv_cyl_ene;
0205 sigv_ = sigv_ * inv_cyl_ene;
0206 sigp_ = sigp_ * inv_cyl_ene;
0207 sige_ = sige_ * inv_cyl_ene;
0208 }
0209 sigu_ = std::sqrt(sigu_);
0210 sigv_ = std::sqrt(sigv_);
0211 sigp_ = std::sqrt(sigp_);
0212 sige_ = std::sqrt(sige_);
0213 }
0214
0215 bool EGammaPCAHelper::checkIteration() const {
0216 if (pcaIteration_ == 0) {
0217 if (debug_)
0218 std::cout << " The PCA has not been run yet " << std::endl;
0219 return false;
0220 } else if (pcaIteration_ == 1) {
0221 if (debug_)
0222 std::cout << " The PCA has been run only once - careful " << std::endl;
0223 return false;
0224 } else if (pcaIteration_ == -1) {
0225 if (debug_)
0226 std::cout << " Not enough layers to perform PCA " << std::endl;
0227 return false;
0228 }
0229 return true;
0230 }
0231
0232 void EGammaPCAHelper::clear() {
0233 theSpots_.clear();
0234 pcaIteration_ = 0;
0235 sigu_ = 0.;
0236 sigv_ = 0.;
0237 sigp_ = 0.;
0238 sige_ = 0.;
0239 }
0240
0241 LongDeps EGammaPCAHelper::energyPerLayer(float radius, bool withHalo) {
0242 if (debug_)
0243 checkIteration();
0244 std::set<int> layers;
0245 float radius2 = radius * radius;
0246 std::vector<float> energyPerLayer(maxlayer_ + 1, 0.f);
0247 math::XYZVector mainAxis(axis_);
0248 mainAxis.unit();
0249 math::XYZVector phiAxis(barycenter_.x(), barycenter_.y(), 0);
0250 math::XYZVector udir(mainAxis.Cross(phiAxis));
0251 udir = udir.unit();
0252 trans_ = Transform3D(Point(barycenter_),
0253 Point(barycenter_ + axis_),
0254 Point(barycenter_ + udir),
0255 Point(0, 0, 0),
0256 Point(0., 0., 1.),
0257 Point(1., 0., 0.));
0258 float energyEE = 0.;
0259 float energyFH = 0.;
0260 float energyBH = 0.;
0261
0262 for (const auto& spot : theSpots_) {
0263 if (!withHalo && !spot.isCore())
0264 continue;
0265 math::XYZPoint local = trans_(Point(spot.row()[0], spot.row()[1], spot.row()[2]));
0266 if (local.Perp2() > radius2)
0267 continue;
0268 energyPerLayer[spot.layer()] += spot.energy();
0269 layers.insert(spot.layer());
0270 if (spot.detId().det() == DetId::HGCalEE or spot.subdet() == HGCEE) {
0271 energyEE += spot.energy();
0272 } else if (spot.detId().det() == DetId::HGCalHSi or spot.subdet() == HGCHEF) {
0273 energyFH += spot.energy();
0274 } else if (spot.detId().det() == DetId::HGCalHSc or spot.subdet() == HGCHEB) {
0275 energyBH += spot.energy();
0276 }
0277 }
0278 return LongDeps(radius, energyPerLayer, energyEE, energyFH, energyBH, layers);
0279 }
0280
0281 void EGammaPCAHelper::printHits(float radius) const {
0282 unsigned nSpots = theSpots_.size();
0283 float radius2 = radius * radius;
0284 for (unsigned i = 0; i < nSpots; ++i) {
0285 Spot spot(theSpots_[i]);
0286 math::XYZPoint local = trans_(Point(spot.row()[0], spot.row()[1], spot.row()[2]));
0287 if (local.Perp2() < radius2) {
0288 std::cout << i << " " << theSpots_[i].detId().rawId() << " " << theSpots_[i].layer() << " "
0289 << theSpots_[i].energy() << " " << theSpots_[i].isCore();
0290 std::cout << " " << std::sqrt(local.Perp2()) << std::endl;
0291 }
0292 }
0293 }
0294
0295 float EGammaPCAHelper::findZFirstLayer(const LongDeps& ld) const {
0296 unsigned int firstLayer = 0;
0297 for (unsigned il = 1; il <= maxlayer_; ++il) {
0298 if (ld.energyPerLayer()[il] > 0.) {
0299 firstLayer = il;
0300 break;
0301 }
0302 }
0303
0304 return recHitTools_->getPositionLayer(firstLayer).z();
0305 }
0306
0307 float EGammaPCAHelper::clusterDepthCompatibility(const LongDeps& ld,
0308 float& measuredDepth,
0309 float& expectedDepth,
0310 float& expectedSigma) {
0311 expectedDepth = -999.;
0312 expectedSigma = -999.;
0313 measuredDepth = -999.;
0314 if (!checkIteration())
0315 return -999.;
0316
0317 float z = findZFirstLayer(ld);
0318 math::XYZVector dir = axis_.unit();
0319 measuredDepth = std::abs((z - std::abs(barycenter_.z())) / dir.z());
0320 return showerDepth_.getClusterDepthCompatibility(measuredDepth, ld.energyEE(), expectedDepth, expectedSigma);
0321 }
