diff --git a/RecoEgamma/EgammaHLTProducers/plugins/BuildFile.xml b/RecoEgamma/EgammaHLTProducers/plugins/BuildFile.xml index 7ccea71c41574..62fe9ec9b516c 100644 --- a/RecoEgamma/EgammaHLTProducers/plugins/BuildFile.xml +++ b/RecoEgamma/EgammaHLTProducers/plugins/BuildFile.xml @@ -4,6 +4,7 @@ + diff --git a/RecoEgamma/EgammaHLTProducers/plugins/EgammaHLTHGCalIDVarProducer.cc b/RecoEgamma/EgammaHLTProducers/plugins/EgammaHLTHGCalIDVarProducer.cc new file mode 100644 index 0000000000000..11e5242e1a6e8 --- /dev/null +++ b/RecoEgamma/EgammaHLTProducers/plugins/EgammaHLTHGCalIDVarProducer.cc @@ -0,0 +1,139 @@ + + +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +#include "DataFormats/RecoCandidate/interface/RecoEcalCandidateIsolation.h" +#include "DataFormats/EgammaReco/interface/SuperCluster.h" +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" + +#include "DataFormats/CaloRecHit/interface/CaloCluster.h" +#include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h" + +#include "DataFormats/EgammaReco/interface/SuperCluster.h" +#include "DataFormats/EgammaReco/interface/SuperClusterFwd.h" + +#include "DataFormats/RecoCandidate/interface/RecoEcalCandidate.h" +#include "DataFormats/RecoCandidate/interface/RecoEcalCandidateFwd.h" + +#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h" + +#include "RecoEgamma/EgammaTools/interface/HGCalShowerShapeHelper.h" +#include "RecoEgamma/EgammaTools/interface/HGCalClusterTools.h" + +class EgammaHLTHGCalIDVarProducer : public edm::stream::EDProducer<> { +public: + explicit EgammaHLTHGCalIDVarProducer(const edm::ParameterSet&); + ~EgammaHLTHGCalIDVarProducer() override; + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + void produce(edm::Event&, const edm::EventSetup&) override; + + class PCAAssocMap { + public: + PCAAssocMap(double HGCalShowerShapeHelper::ShowerWidths::*var, const std::string& name) : var_(var), name_(name) {} + + std::unique_ptr& initMap( + const edm::Handle& candHandle) { + assocMap_ = std::make_unique(candHandle); + return assocMap_; + } + + void insert(reco::RecoEcalCandidateRef& ref, const HGCalShowerShapeHelper::ShowerWidths& showerWidths) { + assocMap_->insert(ref, showerWidths.*var_); + } + + std::unique_ptr& operator()() { return assocMap_; } + const std::string& name() const { return name_; } + + private: + double HGCalShowerShapeHelper::ShowerWidths::*var_; + std::string name_; + std::unique_ptr assocMap_; + }; + +private: + // ----------member data --------------------------- + float rCylinder_; + float hOverECone_; + std::vector pcaAssocMaps_; + const edm::EDGetTokenT recoEcalCandidateToken_; + const edm::EDGetTokenT hgcalRecHitToken_; + const edm::EDGetTokenT layerClusterToken_; +}; + +EgammaHLTHGCalIDVarProducer::EgammaHLTHGCalIDVarProducer(const edm::ParameterSet& config) + : rCylinder_(config.getParameter("rCylinder")), + hOverECone_(config.getParameter("hOverECone")), + recoEcalCandidateToken_( + consumes(config.getParameter("recoEcalCandidateProducer"))), + hgcalRecHitToken_(consumes(config.getParameter("hgcalRecHits"))), + layerClusterToken_(consumes(config.getParameter("layerClusters"))) { + pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2xx, "sigma2xx")); + pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2yy, "sigma2yy")); + pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2zz, "sigma2zz")); + pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2xy, "sigma2xy")); + pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2yz, "sigma2yz")); + pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2zx, "sigma2zx")); + pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2uu, "sigma2uu")); + pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2vv, "sigma2vv")); + pcaAssocMaps_.emplace_back(PCAAssocMap(&HGCalShowerShapeHelper::ShowerWidths::sigma2ww, "sigma2ww")); + + produces("rVar"); + produces("hForHOverE"); + for (auto& var : pcaAssocMaps_) { + produces(var.name()); + } +} + +EgammaHLTHGCalIDVarProducer::~EgammaHLTHGCalIDVarProducer() {} + +void EgammaHLTHGCalIDVarProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + desc.add("recoEcalCandidateProducer", edm::InputTag("hltL1SeededRecoEcalCandidate")); + desc.add("hgcalRecHits", edm::InputTag("hgcalRecHits")); + desc.add("layerClusters", edm::InputTag("layerClusters")); + desc.add("rCylinder", 2.8); + desc.add("hOverECone", 0.15); + descriptions.add(("hltEgammaHLTHGCalIDVarProducer"), desc); +} + +void EgammaHLTHGCalIDVarProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) { + auto recoEcalCandHandle = iEvent.getHandle(recoEcalCandidateToken_); + auto hgcalRecHitHandle = iEvent.getHandle(hgcalRecHitToken_); + auto layerClustersHandle = iEvent.getHandle(layerClusterToken_); + + HGCalShowerShapeHelper ssCalc; + ssCalc.initPerEvent(iSetup, *hgcalRecHitHandle); + + auto rVarMap = std::make_unique(recoEcalCandHandle); + auto hForHoverEMap = std::make_unique(recoEcalCandHandle); + for (auto& pcaMap : pcaAssocMaps_) { + pcaMap.initMap(recoEcalCandHandle); + } + + for (size_t candNr = 0; candNr < recoEcalCandHandle->size(); candNr++) { + reco::RecoEcalCandidateRef candRef(recoEcalCandHandle, candNr); + ssCalc.initPerObject(candRef->superCluster()->hitsAndFractions()); + rVarMap->insert(candRef, ssCalc.getRvar(rCylinder_, candRef->superCluster()->energy())); + + float hForHoverE = HGCalClusterTools::hadEnergyInCone( + candRef->superCluster()->eta(), candRef->superCluster()->phi(), *layerClustersHandle, 0., hOverECone_, 0., 0.); + hForHoverEMap->insert(candRef, hForHoverE); + auto pcaWidths = ssCalc.getPCAWidths(rCylinder_); + for (auto& pcaMap : pcaAssocMaps_) { + pcaMap.insert(candRef, pcaWidths); + } + } + iEvent.put(std::move(rVarMap), "rVar"); + iEvent.put(std::move(hForHoverEMap), "hForHOverE"); + for (auto& pcaMap : pcaAssocMaps_) { + iEvent.put(std::move(pcaMap()), pcaMap.name()); + } +} + +DEFINE_FWK_MODULE(EgammaHLTHGCalIDVarProducer); diff --git a/RecoEgamma/EgammaTools/interface/HGCalClusterTools.h b/RecoEgamma/EgammaTools/interface/HGCalClusterTools.h new file mode 100644 index 0000000000000..bfb94fa14d8e2 --- /dev/null +++ b/RecoEgamma/EgammaTools/interface/HGCalClusterTools.h @@ -0,0 +1,44 @@ +#ifndef RecoEgamma_EgammaTools_HGCalClusterTools_h +#define RecoEgamma_EgammaTools_HGCalClusterTools_h + +#include "DataFormats/CaloRecHit/interface/CaloClusterFwd.h" +#include + +class HGCalClusterTools { +public: + enum class EType { ET, ENERGY }; + + static float energyInCone(const float eta, + const float phi, + const std::vector& layerClusters, + const float minDR, + const float maxDR, + const float minEt, + const float minEnergy, + const std::vector& subDets, + const HGCalClusterTools::EType& eType = EType::ENERGY); + + static float hadEnergyInCone(const float eta, + const float phi, + const std::vector& layerClusters, + const float minDR, + const float maxDR, + const float minEt, + const float minEnergy, + const HGCalClusterTools::EType& eType = EType::ENERGY) { + return energyInCone( + eta, phi, layerClusters, minDR, maxDR, minEt, minEnergy, {DetId::HGCalHSi, DetId::HGCalHSc}, eType); + } + static float emEnergyInCone(const float eta, + const float phi, + const std::vector& layerClusters, + const float minDR, + const float maxDR, + const float minEt, + const float minEnergy, + const HGCalClusterTools::EType& eType = EType::ENERGY) { + return energyInCone(eta, phi, layerClusters, minDR, maxDR, minEt, minEnergy, {DetId::HGCalEE}, eType); + } +}; + +#endif diff --git a/RecoEgamma/EgammaTools/interface/HGCalShowerShapeHelper.h b/RecoEgamma/EgammaTools/interface/HGCalShowerShapeHelper.h new file mode 100644 index 0000000000000..a86e084a4b838 --- /dev/null +++ b/RecoEgamma/EgammaTools/interface/HGCalShowerShapeHelper.h @@ -0,0 +1,131 @@ +#ifndef RecoEgamma_EgammaTools_HGCalShowerShapeHelper_h +#define RecoEgamma_EgammaTools_HGCalShowerShapeHelper_h + +// system include files +#include + +// user include files +#include "DataFormats/CaloRecHit/interface/CaloCluster.h" +#include "DataFormats/EgammaCandidates/interface/GsfElectron.h" +#include "DataFormats/FWLite/interface/ESHandle.h" +#include "DataFormats/ForwardDetId/interface/HGCEEDetId.h" +#include "DataFormats/ForwardDetId/interface/HGCalDetId.h" +#include "DataFormats/GsfTrackReco/interface/GsfTrack.h" +#include "DataFormats/HGCRecHit/interface/HGCRecHit.h" +#include "DataFormats/Math/interface/LorentzVector.h" +#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h" +#include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/MakerMacros.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Utilities/interface/StreamID.h" +#include "Geometry/CaloTopology/interface/HGCalTopology.h" +#include "Geometry/Records/interface/IdealGeometryRecord.h" +#include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h" +#include "RecoParticleFlow/PFClusterProducer/interface/InitialClusteringStepBase.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + +class HGCalShowerShapeHelper { +private: + // Good to filter/compute/store this stuff beforehand as they are common to the shower shape variables. + // No point in filtering, computing layer-wise centroids, etc. for each variable again and again. + // Once intitialized, one can the calculate different variables one after another for a given object. + // If a different set of preselections (E, ET, etc.) is required for a given object, then reinitialize using initPerObject(...). + + // In principle should consider the HGCalHSi and HGCalHSc hits (leakage) also. + // Can have subdetector dependent thresholds and layer selection. + // To be implemented. + + double minHitE_; + double minHitET_; + double minHitET2_; + int minLayer_; + int maxLayer_; + int nLayer_; + DetId::Detector subDet_; + + hgcal::RecHitTools recHitTools_; + + std::unordered_map pfRecHitPtrMap_; + std::vector > hitsAndFracs_; + std::vector hitEnergies_; + std::vector hitEnergiesWithFracs_; + + ROOT::Math::XYZVector centroid_; + std::vector layerEnergies_; + std::vector layerCentroids_; + + void setPFRecHitPtrMap(const std::vector &recHits); + + void setFilteredHitsAndFractions(const std::vector > &hitsAndFracs); + +public: + static constexpr double kLDWaferCellSize_ = 0.698; + static constexpr double kHDWaferCellSize_ = 0.465; + + void setLayerWiseStuff(); + + struct ShowerWidths { + double sigma2xx; + double sigma2yy; + double sigma2zz; + + double sigma2xy; + double sigma2yz; + double sigma2zx; + + double sigma2uu; + double sigma2vv; + double sigma2ww; + + ShowerWidths() + : sigma2xx(0.0), + sigma2yy(0.0), + sigma2zz(0.0), + sigma2xy(0.0), + sigma2yz(0.0), + sigma2zx(0.0), + sigma2uu(0.0), + sigma2vv(0.0), + sigma2ww(0.0) {} + }; + + void initPerEvent(const edm::EventSetup &iSetup, const std::vector &recHits); + + void initPerObject(const std::vector > &hitsAndFracs, + double minHitE = 0, + double minHitET = 0, + int minLayer = 1, + int maxLayer = 28, + DetId::Detector subDet = DetId::HGCalEE); + + const double getCellSize(DetId detId); + + // Compute Rvar in a cylinder around the layer centroids + const double getRvar(double cylinderR, double energyNorm, bool useFractions = true, bool useCellSize = true); + + // Compute PCA widths around the layer centroids + const ShowerWidths getPCAWidths(double cylinderR, bool useFractions = false); +}; + +#endif diff --git a/RecoEgamma/EgammaTools/src/HGCalClusterTools.cc b/RecoEgamma/EgammaTools/src/HGCalClusterTools.cc new file mode 100644 index 0000000000000..abc06ae0ad512 --- /dev/null +++ b/RecoEgamma/EgammaTools/src/HGCalClusterTools.cc @@ -0,0 +1,56 @@ +#include "RecoEgamma/EgammaTools/interface/HGCalClusterTools.h" + +#include "DataFormats/CaloRecHit/interface/CaloCluster.h" +#include "DataFormats/Math/interface/deltaR.h" +#include "DataFormats/Math/interface/deltaPhi.h" + +float HGCalClusterTools::energyInCone(const float eta, + const float phi, + const std::vector& layerClusters, + const float minDR, + const float maxDR, + const float minEt, + const float minEnergy, + const std::vector& subDets, + const HGCalClusterTools::EType& eType) { + float hadValue = 0.; + + const float minDR2 = minDR * minDR; + const float maxDR2 = maxDR * maxDR; + + for (auto& clus : layerClusters) { + if (clus.energy() < minEnergy) { + continue; + } + + if (std::find(subDets.begin(), subDets.end(), clus.seed().det()) == subDets.end()) { + continue; + } + + float clusEt = clus.energy() * std::sin(clus.position().theta()); + if (clusEt < minEt) { + continue; + } + + //this is a prefilter on the clusters before we calculuate + //the expensive eta() of the cluster + float dPhi = reco::deltaPhi(phi, clus.phi()); + if (dPhi > maxDR) { + continue; + } + + float dR2 = reco::deltaR2(eta, phi, clus.eta(), clus.phi()); + if (dR2 < minDR2 || dR2 > maxDR2) { + continue; + } + switch (eType) { + case EType::ET: + hadValue += clusEt; + break; + case EType::ENERGY: + hadValue += clus.energy(); + break; + } + } + return hadValue; +} diff --git a/RecoEgamma/EgammaTools/src/HGCalShowerShapeHelper.cc b/RecoEgamma/EgammaTools/src/HGCalShowerShapeHelper.cc new file mode 100644 index 0000000000000..a9937e1a71b4c --- /dev/null +++ b/RecoEgamma/EgammaTools/src/HGCalShowerShapeHelper.cc @@ -0,0 +1,284 @@ +#include "RecoEgamma/EgammaTools/interface/HGCalShowerShapeHelper.h" + +void HGCalShowerShapeHelper::initPerEvent(const edm::EventSetup &iSetup, const std::vector &pfRecHits) { + edm::ESHandle geom; + iSetup.get().get(geom); + recHitTools_.setGeometry(*(geom.product())); + + setPFRecHitPtrMap(pfRecHits); +} + +void HGCalShowerShapeHelper::initPerObject(const std::vector > &hitsAndFracs, + double minHitE, + double minHitET, + int minLayer, + int maxLayer, + DetId::Detector subDet) { + // Safety checks + nLayer_ = maxLayer - minLayer + 1; + assert(nLayer_ > 0); + + minHitE_ = minHitE; + minHitET_ = minHitET; + minHitET2_ = minHitET * minHitET; + minLayer_ = minLayer; + maxLayer_ = maxLayer; + subDet_ = subDet; + + setFilteredHitsAndFractions(hitsAndFracs); + + setLayerWiseStuff(); +} + +void HGCalShowerShapeHelper::setPFRecHitPtrMap(const std::vector &recHits) { + pfRecHitPtrMap_.clear(); + + for (const auto &recHit : recHits) { + pfRecHitPtrMap_[recHit.detId()] = &recHit; + } +} + +void HGCalShowerShapeHelper::setFilteredHitsAndFractions(const std::vector > &hitsAndFracs) { + hitsAndFracs_.clear(); + hitEnergies_.clear(); + hitEnergiesWithFracs_.clear(); + + for (const auto &hnf : hitsAndFracs) { + DetId hitId = hnf.first; + float hitEfrac = hnf.second; + + int hitLayer = recHitTools_.getLayer(hitId); + + if (hitLayer > nLayer_) { + continue; + } + + if (hitId.det() != subDet_) { + continue; + } + + if (pfRecHitPtrMap_.find(hitId.rawId()) == pfRecHitPtrMap_.end()) { + continue; + } + + reco::PFRecHit recHit = *pfRecHitPtrMap_[hitId.rawId()]; + + if (recHit.energy() < minHitE_) { + continue; + } + + if (recHit.pt2() < minHitET2_) { + continue; + } + + // Fill the vectors + hitsAndFracs_.push_back(hnf); + hitEnergies_.push_back(recHit.energy()); + hitEnergiesWithFracs_.push_back(recHit.energy() * hitEfrac); + } +} + +void HGCalShowerShapeHelper::setLayerWiseStuff() { + layerEnergies_.clear(); + layerEnergies_.resize(nLayer_); + + layerCentroids_.clear(); + layerCentroids_.resize(nLayer_); + + centroid_.SetXYZ(0, 0, 0); + + int iHit = -1; + double totalW = 0.0; + + // Compute the centroid per layer + for (const auto &hnf : hitsAndFracs_) { + iHit++; + + DetId hitId = hnf.first; + + double weight = hitEnergies_[iHit]; + totalW += weight; + + const auto &hitPos = recHitTools_.getPosition(hitId); + ROOT::Math::XYZVector hitXYZ(hitPos.x(), hitPos.y(), hitPos.z()); + + centroid_ += weight * hitXYZ; + + int hitLayer = recHitTools_.getLayer(hitId) - 1; + + layerEnergies_[hitLayer] += weight; + layerCentroids_[hitLayer] += weight * hitXYZ; + } + + int iLayer = -1; + + for (auto ¢roid : layerCentroids_) { + iLayer++; + + if (layerEnergies_[iLayer]) { + centroid /= layerEnergies_[iLayer]; + } + } + + if (totalW) { + centroid_ /= totalW; + } +} + +const double HGCalShowerShapeHelper::getCellSize(DetId detId) { + double siThickness = recHitTools_.getSiThickness(detId); + + return siThickness < 150 ? kHDWaferCellSize_ : kLDWaferCellSize_; +} + +const double HGCalShowerShapeHelper::getRvar(double cylinderR, double energyNorm, bool useFractions, bool useCellSize) { + if (hitsAndFracs_.empty()) { + return 0; + } + + if (energyNorm <= 0) { + edm::LogWarning("HGCalShowerShapeHelper") + << "Encountered negative or zero energy for HGCal R-variable denomintor: " << energyNorm << std::endl; + } + + double cylinderR2 = cylinderR * cylinderR; + + double Rvar = 0; + + auto hitEnergyIter = useFractions ? hitEnergiesWithFracs_.begin() : hitEnergies_.begin(); + + hitEnergyIter--; + + for (const auto &hnf : hitsAndFracs_) { + hitEnergyIter++; + + DetId hitId = hnf.first; + + int hitLayer = recHitTools_.getLayer(hitId) - 1; + + const auto &hitPos = recHitTools_.getPosition(hitId); + ROOT::Math::XYZVector hitXYZ(hitPos.x(), hitPos.y(), hitPos.z()); + + auto distXYZ = hitXYZ - layerCentroids_[hitLayer]; + + double r2 = distXYZ.x() * distXYZ.x() + distXYZ.y() * distXYZ.y(); + + // Including the cell size seems to make the variable less sensitive to the HD/LD transition region + if (useCellSize) { + if (std::sqrt(r2) > cylinderR + getCellSize(hitId)) { + continue; + } + } + + else if (r2 > cylinderR2) { + continue; + } + + Rvar += *hitEnergyIter; + } + + Rvar /= energyNorm; + + return Rvar; +} + +const HGCalShowerShapeHelper::ShowerWidths HGCalShowerShapeHelper::getPCAWidths(double cylinderR, bool useFractions) { + if (hitsAndFracs_.empty()) { + return ShowerWidths(); + } + + double cylinderR2 = cylinderR * cylinderR; + + TMatrixD covMat(3, 3); + + double dxdx = 0.0; + double dydy = 0.0; + double dzdz = 0.0; + + double dxdy = 0.0; + double dydz = 0.0; + double dzdx = 0.0; + + double totalW = 0; + + auto hitEnergyIter = useFractions ? hitEnergiesWithFracs_.begin() : hitEnergies_.begin(); + + int nHit = 0; + hitEnergyIter--; + + for (const auto &hnf : hitsAndFracs_) { + hitEnergyIter++; + + DetId hitId = hnf.first; + + const auto &hitPos = recHitTools_.getPosition(hitId); + ROOT::Math::XYZVector hitXYZ(hitPos.x(), hitPos.y(), hitPos.z()); + + int hitLayer = recHitTools_.getLayer(hitId) - 1; + + ROOT::Math::XYZVector radXYZ = hitXYZ - layerCentroids_[hitLayer]; + + double r2 = radXYZ.x() * radXYZ.x() + radXYZ.y() * radXYZ.y(); + + if (r2 > cylinderR2) { + continue; + } + + ROOT::Math::XYZVector dXYZ = hitXYZ - centroid_; + + double weight = *hitEnergyIter; + totalW += weight; + + dxdx += weight * dXYZ.x() * dXYZ.x(); + dydy += weight * dXYZ.y() * dXYZ.y(); + dzdz += weight * dXYZ.z() * dXYZ.z(); + + dxdy += weight * dXYZ.x() * dXYZ.y(); + dydz += weight * dXYZ.y() * dXYZ.z(); + dzdx += weight * dXYZ.z() * dXYZ.x(); + + nHit++; + } + + if (!totalW || nHit < 2) { + return ShowerWidths(); + } + + dxdx /= totalW; + dydy /= totalW; + dzdz /= totalW; + + dxdy /= totalW; + dydz /= totalW; + dzdx /= totalW; + + covMat(0, 0) = dxdx; + covMat(1, 1) = dydy; + covMat(2, 2) = dzdz; + + covMat(0, 1) = covMat(1, 0) = dxdy; + covMat(0, 2) = covMat(2, 0) = dzdx; + covMat(1, 2) = covMat(2, 1) = dydz; + + // Get eigen values and vectors + TVectorD eigVals(3); + TMatrixD eigVecMat(3, 3); + + eigVecMat = covMat.EigenVectors(eigVals); + + ShowerWidths returnWidths; + + returnWidths.sigma2xx = dxdx; + returnWidths.sigma2yy = dydy; + returnWidths.sigma2zz = dzdz; + + returnWidths.sigma2xy = dxdy; + returnWidths.sigma2yz = dydz; + returnWidths.sigma2zx = dzdx; + + returnWidths.sigma2uu = eigVals(1); + returnWidths.sigma2vv = eigVals(2); + returnWidths.sigma2ww = eigVals(0); + + return returnWidths; +}