File indexing completed on 2023-03-17 10:45:23
0001 #include "CommonTools/ParticleFlow/interface/PFClusterWidthAlgo.h"
0002 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
0003 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h"
0004 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0005 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0006 #include "DataFormats/Math/interface/deltaPhi.h"
0007 #include "TMath.h"
0008 #include <algorithm>
0009 using namespace std;
0010 using namespace reco;
0011
0012 PFClusterWidthAlgo::PFClusterWidthAlgo(const std::vector<const reco::PFCluster*>& pfclust) {
0013 double numeratorEtaWidth = 0.;
0014 double numeratorPhiWidth = 0.;
0015 double sclusterE = 0.;
0016 double posX = 0.;
0017 double posY = 0.;
0018 double posZ = 0.;
0019 sigmaEtaEta_ = 0.;
0020
0021 unsigned int nclust = pfclust.size();
0022 if (nclust == 0) {
0023 etaWidth_ = 0.;
0024 phiWidth_ = 0.;
0025 sigmaEtaEta_ = 0.;
0026 } else {
0027
0028
0029 for (unsigned int icl = 0; icl < nclust; ++icl) {
0030 const std::vector<reco::PFRecHitFraction>& PFRecHits = pfclust[icl]->recHitFractions();
0031
0032 for (std::vector<reco::PFRecHitFraction>::const_iterator it = PFRecHits.begin(); it != PFRecHits.end(); ++it) {
0033 const PFRecHitRef& RefPFRecHit = it->recHitRef();
0034
0035 double energyHit = RefPFRecHit->energy() * it->fraction();
0036
0037 sclusterE += energyHit;
0038 posX += energyHit * RefPFRecHit->position().x();
0039 posY += energyHit * RefPFRecHit->position().y();
0040 posZ += energyHit * RefPFRecHit->position().z();
0041 }
0042 }
0043
0044 double denominator = sclusterE;
0045
0046 posX /= sclusterE;
0047 posY /= sclusterE;
0048 posZ /= sclusterE;
0049
0050 math::XYZPoint pflowSCPos(posX, posY, posZ);
0051
0052 double scEta = pflowSCPos.eta();
0053 double scPhi = pflowSCPos.phi();
0054
0055 double SeedClusEnergy = -1.;
0056 unsigned int SeedDetID = 0;
0057 double SeedEta = -1.;
0058
0059
0060 for (unsigned int icl = 0; icl < nclust; ++icl) {
0061 const auto& PFRecHits = pfclust[icl]->recHitFractions();
0062
0063 for (auto it = PFRecHits.begin(); it != PFRecHits.end(); ++it) {
0064 const PFRecHitRef& RefPFRecHit = it->recHitRef();
0065
0066 double energyHit = RefPFRecHit->energy() * it->fraction();
0067
0068
0069 if (icl == 0) {
0070 if (energyHit > SeedClusEnergy) {
0071 SeedClusEnergy = energyHit;
0072 SeedEta = RefPFRecHit->position().eta();
0073 SeedDetID = RefPFRecHit->detId();
0074 }
0075 }
0076
0077 double dPhi = reco::deltaPhi(RefPFRecHit->positionREP().phi(), scPhi);
0078 double dEta = RefPFRecHit->positionREP().eta() - scEta;
0079 numeratorEtaWidth += energyHit * dEta * dEta;
0080 numeratorPhiWidth += energyHit * dPhi * dPhi;
0081 }
0082 }
0083
0084
0085 const auto& PFRecHits = pfclust[0]->recHitFractions();
0086 for (auto it = PFRecHits.begin(); it != PFRecHits.end(); ++it) {
0087 const auto& RefPFRecHit = it->recHitRef();
0088 if (!RefPFRecHit.isAvailable())
0089 return;
0090 double energyHit = RefPFRecHit->energy();
0091 if (RefPFRecHit->detId() != SeedDetID) {
0092 float diffEta = RefPFRecHit->positionREP().eta() - SeedEta;
0093 sigmaEtaEta_ += (diffEta * diffEta) * (energyHit / SeedClusEnergy);
0094 }
0095 }
0096 if (sigmaEtaEta_ == 0.)
0097 sigmaEtaEta_ = 0.00000001;
0098
0099 etaWidth_ = std::sqrt(numeratorEtaWidth / denominator);
0100 phiWidth_ = std::sqrt(numeratorPhiWidth / denominator);
0101
0102 }
0103 }
0104 PFClusterWidthAlgo::~PFClusterWidthAlgo() {}