Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:23

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 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;  // convert in GeV
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   // very important - to reset
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       // initial calculation, take only core hits
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       // use a cylinder, include all hits
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     // Select halo hits or not
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   // Make dummy DetId to get abs(z) for layer
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 }