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