ParticleFlow/src/PFClusterWidthAlgo.cc

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     //first loop, compute supercluster position at ecal face, and energy sum from rechit loop
0028     //in order to be consistent with variance calculation
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         //compute rechit energy taking into account fractions
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     }  // end for ncluster
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     //second loop, compute variances
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         //compute rechit energy taking into account fractions
0066         double energyHit = RefPFRecHit->energy() * it->fraction();
0067
0068         //only for the first cluster (from GSF) find the seed
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     }  // end for ncluster
0083
0084     //for the first cluster (from GSF) computed sigmaEtaEta
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   }  // endif ncluster > 0
0103 }
0104 PFClusterWidthAlgo::~PFClusterWidthAlgo() {}