Back to home page

Project CMSSW displayed by LXR

 
 

    


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     :  // Thickness correction to dEdx weights
0019       // (100um, 200um, 300um silicon)
0020       // See RecoLocalCalo.HGCalRecProducers.HGCalRecHit_cfi
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;  // convert in GeV
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   // very important - to reset
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       // initial calculation, take only core hits
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       // use a cylinder, include all hits
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     // Select halo hits or not
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   // Make dummy DetId to get abs(z) for layer
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 }