#include <PFClusterWidthAlgo.h>

Public Member Functions
	PFClusterWidthAlgo (const std::vector< const reco::PFCluster * > &pfclust, const EBRecHitCollection ebRecHits=0, const EERecHitCollection eeRecHits=0)
double	pflowEtaWidth () const
double	pflowPhiWidth () const
double	pflowSigmaEtaEta () const
	~PFClusterWidthAlgo ()
Private Attributes
double	etaWidth_
double	phiWidth_
double	sigmaEtaEta_

Detailed Description

Definition at line 6 of file PFClusterWidthAlgo.h.

Constructor & Destructor Documentation

PFClusterWidthAlgo::PFClusterWidthAlgo	(	const std::vector< const reco::PFCluster * > &	pfclust,
		const EBRecHitCollection *	ebRecHits = `0`,
		const EERecHitCollection *	eeRecHits = `0`
	)

Definition at line 13 of file PFClusterWidthAlgo.cc.

References ExpressReco_HICollisions_FallBack::denominator, dPhi(), ExpressReco_HICollisions_FallBack::e, DetId::Ecal, EcalBarrel, EcalEndcap, edm::SortedCollection< T, SORT >::end(), edm::SortedCollection< T, SORT >::find(), getHLTprescales::index, edm::Ref< C, T, F >::isAvailable(), ExpressReco_HICollisions_FallBack::PFRecHits, Pi, mathSSE::sqrt(), and TwoPi.

                                                                            {


  double numeratorEtaWidth = 0.;
  double numeratorPhiWidth = 0.;
  double sclusterE = 0.;
  double posX = 0.;
  double posY = 0.;
  double posZ = 0.;
  sigmaEtaEta_ = 0.;

  unsigned int nclust= pfclust.size();
  if(nclust == 0 ) {
    etaWidth_ = 0.;
    phiWidth_ = 0.;
    sigmaEtaEta_ = 0.;
  }
  else {

    for(unsigned int icl=0;icl<nclust;++icl) {
      double e = pfclust[icl]->energy();
      sclusterE += e;
      posX += e * pfclust[icl]->position().X();
      posY += e * pfclust[icl]->position().Y();
      posZ += e * pfclust[icl]->position().Z();   
    }
    
    posX /=sclusterE;
    posY /=sclusterE;
    posZ /=sclusterE;
    
    double denominator = sclusterE;
    
    math::XYZPoint pflowSCPos(posX,posY,posZ);
    
    double scEta    = pflowSCPos.eta();
    double scPhi    = pflowSCPos.phi();
    
    double SeedClusEnergy = -1.;
    unsigned int SeedDetID = 0;
    double SeedEta = -1.;
    double SeedPhi = -1.;

    for(unsigned int icl=0; icl<nclust; ++icl) {
      const std::vector< reco::PFRecHitFraction >& PFRecHits =  pfclust[icl]->recHitFractions();
      
      
      for ( std::vector< reco::PFRecHitFraction >::const_iterator it = PFRecHits.begin(); 
            it != PFRecHits.end(); ++it) {
        const PFRecHitRef& RefPFRecHit = it->recHitRef(); 
        double energyHit=0;
        // if the PFRecHit is not available, try to get it from the collections
        //      if(!RefPFRecHit.isAvailable() && ebRecHits && eeRecHits )
        if(ebRecHits && eeRecHits )
          {
            //      std::cout << " Recomputing " << std::endl;
            unsigned index=it-PFRecHits.begin();
            DetId id=pfclust[icl]->hitsAndFractions()[index].first;
            // look for the hit; do not forget to multiply by the fraction
            if(id.det()==DetId::Ecal && id.subdetId()==EcalBarrel)
              {
                EBRecHitCollection::const_iterator itcheck=ebRecHits->find(id);
                if(itcheck!=ebRecHits->end())
                  energyHit= itcheck->energy();
                // The cluster shapes do not take into account the RecHit fraction !            
                // * pfclust[icl]->hitsAndFractions()[index].second;
              }
            if(id.det()==DetId::Ecal && id.subdetId()==EcalEndcap)
              {
                EERecHitCollection::const_iterator itcheck=eeRecHits->find(id);
                if(itcheck!=eeRecHits->end())
                  energyHit= itcheck->energy();
                // The cluster shapes do not take into account the RecHit fraction !            
                //               * pfclust[icl]->hitsAndFractions()[index].second;
              }
          }
        else
          energyHit = RefPFRecHit->energy();

        //only for the first cluster (from GSF) find the seed
        if(icl==0) {
          if (energyHit > SeedClusEnergy) {
            SeedClusEnergy = energyHit;
            SeedEta = RefPFRecHit->position().eta();
            SeedPhi =  RefPFRecHit->position().phi();
            SeedDetID = RefPFRecHit->detId();
          }
        }


        double dPhi = RefPFRecHit->position().phi() - scPhi;
        if (dPhi > + TMath::Pi()) { dPhi = TMath::TwoPi() - dPhi; }
        if (dPhi < - TMath::Pi()) { dPhi = TMath::TwoPi() + dPhi; }
        double dEta = RefPFRecHit->position().eta() - scEta;
        if ( energyHit > 0 ) {
          numeratorEtaWidth += energyHit * dEta * dEta;
          numeratorPhiWidth += energyHit * dPhi * dPhi;
        }
      }
    } // end for ncluster

    //for the first cluster (from GSF) computed sigmaEtaEta
    const std::vector< reco::PFRecHitFraction >& PFRecHits =  pfclust[0]->recHitFractions();
    for ( std::vector< reco::PFRecHitFraction >::const_iterator it = PFRecHits.begin(); 
          it != PFRecHits.end(); ++it) {
      const PFRecHitRef& RefPFRecHit = it->recHitRef(); 
      if(!RefPFRecHit.isAvailable()) 
        return;
      double energyHit = RefPFRecHit->energy();
      if (RefPFRecHit->detId() != SeedDetID) {
        float diffEta =  RefPFRecHit->position().eta() - SeedEta;
        sigmaEtaEta_ += (diffEta*diffEta) * (energyHit/SeedClusEnergy);
      }
    }
    if (sigmaEtaEta_ == 0.) sigmaEtaEta_ = 0.00000001;

    etaWidth_ = sqrt(numeratorEtaWidth / denominator);
    phiWidth_ = sqrt(numeratorPhiWidth / denominator);
    

  } // endif ncluster > 0
}