/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/RecoEcal/EgammaClusterAlgos/src/PFECALSuperClusterAlgo.cc

Go to the documentation of this file.

#include "RecoEcal/EgammaClusterAlgos/interface/PFECALSuperClusterAlgo.h"
#include "RecoParticleFlow/PFClusterTools/interface/PFClusterWidthAlgo.h"
#include "RecoParticleFlow/PFClusterTools/interface/LinkByRecHit.h"
#include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "RecoEcal/EgammaCoreTools/interface/Mustache.h"
#include "Math/GenVector/VectorUtil.h"
#include "TFile.h"
#include "TH2F.h"
#include "TROOT.h"
#include "TMath.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include <stdexcept>
#include <string>
#include <sstream>
#include <cmath>

using namespace std;
namespace MK = reco::MustacheKernel;


namespace {
  typedef edm::View<reco::PFCluster> PFClusterView;
  typedef edm::Ptr<reco::PFCluster> PFClusterPtr;
  typedef edm::PtrVector<reco::PFCluster> PFClusterPtrVector;
  typedef PFECALSuperClusterAlgo::CalibratedClusterPtr CalibClusterPtr;
  typedef PFECALSuperClusterAlgo::CalibratedClusterPtrVector CalibClusterPtrVector;

  typedef std::binary_function<const CalibClusterPtr&,
                               const CalibClusterPtr&, 
                               bool> ClusBinaryFunction;
  
  typedef std::unary_function<const CalibClusterPtr&, 
                              bool> ClusUnaryFunction;  

  struct SumPSEnergy : public std::binary_function<double,
                                                   const PFClusterPtr&,
                                                   double> {
    PFLayer::Layer _thelayer;
    SumPSEnergy(PFLayer::Layer layer) : _thelayer(layer) {}
    double operator()(double a,
                      const PFClusterPtr& b) {
      return a + (_thelayer == b->layer())*b->energy();
    }
  };

  struct GreaterByE : public ClusBinaryFunction {
    bool operator()(const CalibClusterPtr& x, 
                    const CalibClusterPtr& y) { 
      return x->energy() > y->energy() ; 
    }
  };

  struct GreaterByEt : public ClusBinaryFunction {
    bool operator()(const CalibClusterPtr& x, 
                    const CalibClusterPtr& y) { 
      return x->energy()/std::cosh(x->eta()) > y->energy()/std::cosh(y->eta());
    }
  };

  struct IsASeed : public ClusUnaryFunction {
    double threshold;
    IsASeed(double thresh) : threshold(thresh) {}
    bool operator()(const CalibClusterPtr& x) { 
      return x->energy() > threshold; 
    }
  };

  struct IsLinkedByRecHit : public ClusUnaryFunction {
    const CalibClusterPtr the_seed;
    const double _threshold, _majority;
    const double _maxSatelliteDEta, _maxSatelliteDPhi;
    double x_rechits_tot, x_rechits_match;
    IsLinkedByRecHit(const CalibClusterPtr& s, const double threshold,
                     const double majority, const double maxDEta,
                     const double maxDPhi) : 
      the_seed(s),_threshold(threshold),_majority(majority), 
      _maxSatelliteDEta(maxDEta), _maxSatelliteDPhi(maxDPhi) {}
    bool operator()(const CalibClusterPtr& x) {      
      if( the_seed->energy_nocalib() < _threshold ) return false; 
      const double dEta = std::abs(the_seed->eta()-x->eta());
      const double dPhi = 
        std::abs(TVector2::Phi_mpi_pi(the_seed->phi() - x->phi())); 
      if( _maxSatelliteDEta < dEta || _maxSatelliteDPhi < dPhi) return false;
      // now see if the clusters overlap in rechits
      const auto& seedHitsAndFractions = 
        the_seed->the_ptr()->hitsAndFractions();
      const auto& xHitsAndFractions = 
        x->the_ptr()->hitsAndFractions();      
      x_rechits_tot   = xHitsAndFractions.size();
      x_rechits_match = 0.0;      
      for( const std::pair<DetId, float>& seedHit : seedHitsAndFractions ) {
        for( const std::pair<DetId, float>& xHit : xHitsAndFractions ) {
          if( seedHit.first == xHit.first ) {       
            x_rechits_match += 1.0;
          }
        }       
      }      
      return x_rechits_match/x_rechits_tot > _majority;
    }
  };

  struct IsClustered : public ClusUnaryFunction {
    const CalibClusterPtr the_seed;    
    PFECALSuperClusterAlgo::clustering_type _type;
    bool dynamic_dphi;
    double etawidthSuperCluster_, phiwidthSuperCluster_;
    IsClustered(const CalibClusterPtr s, 
                PFECALSuperClusterAlgo::clustering_type ct,
                const bool dyn_dphi) : 
      the_seed(s), _type(ct), dynamic_dphi(dyn_dphi) {}
    bool operator()(const CalibClusterPtr& x) { 
      const double dphi = 
        std::abs(TVector2::Phi_mpi_pi(the_seed->phi() - x->phi()));  
      const bool isEB = ( PFLayer::ECAL_BARREL == x->the_ptr()->layer() );
      const bool passes_dphi = 
        ( (!dynamic_dphi && dphi < phiwidthSuperCluster_ ) || 
          (dynamic_dphi && MK::inDynamicDPhiWindow(isEB,the_seed->phi(),
                                                   x->energy_nocalib(),
                                                   x->eta(),
                                                   x->phi()) ) );

      switch( _type ) {
      case PFECALSuperClusterAlgo::kBOX:
        return ( std::abs(the_seed->eta()-x->eta())<etawidthSuperCluster_ && 
                 passes_dphi   );
        break;
      case PFECALSuperClusterAlgo::kMustache:
        return ( passes_dphi && 
                 MK::inMustache(the_seed->eta(), 
                                the_seed->phi(),
                                x->energy_nocalib(),
                                x->eta(),
                                x->phi()            ));
        break;
      default: 
        return false;
      }
      return false;
    }
  };

  double testPreshowerDistance(const edm::Ptr<reco::PFCluster>& eeclus,
                               const edm::Ptr<reco::PFCluster>& psclus) {
    if( psclus.isNull() || eeclus.isNull() ) return -1.0;
    /* 
    // commented out since PFCluster::layer() uses a lot of CPU
    // and since 
    if( PFLayer::ECAL_ENDCAP != eeclus->layer() ) return -1.0;
    if( PFLayer::PS1 != psclus->layer() &&
        PFLayer::PS2 != psclus->layer()    ) {
      throw cms::Exception("testPreshowerDistance")
        << "The second argument passed to this function was "
        << "not a preshower cluster!" << std::endl;
    } 
    */
    const reco::PFCluster::REPPoint& pspos = psclus->positionREP();
    const reco::PFCluster::REPPoint& eepos = eeclus->positionREP();
    // lazy continue based on geometry
    if( eeclus->z()*psclus->z() < 0 ) return -1.0;
    const double dphi= std::abs(TVector2::Phi_mpi_pi(eepos.phi() - 
                                                     pspos.phi()));
    if( dphi > 0.6 ) return -1.0;    
    const double deta= std::abs(eepos.eta() - pspos.eta());    
    if( deta > 0.3 ) return -1.0; 
    return LinkByRecHit::testECALAndPSByRecHit(*eeclus,*psclus,false);
  }
}

PFECALSuperClusterAlgo::PFECALSuperClusterAlgo() { }

void PFECALSuperClusterAlgo::
setPFClusterCalibration(const std::shared_ptr<PFEnergyCalibration>& calib) {
  _pfEnergyCalibration = calib;
}

void PFECALSuperClusterAlgo::
loadAndSortPFClusters(const PFClusterView& clusters,
                      const PFClusterView& psclusters) { 
  // reset the system for running
  superClustersEB_.reset(new reco::SuperClusterCollection);
  _clustersEB.clear();
  superClustersEE_.reset(new reco::SuperClusterCollection);  
  _clustersEE.clear();
  _psclustersforee.clear();
  
  auto clusterPtrs = clusters.ptrVector(); 
  //Select PF clusters available for the clustering
  for ( auto& cluster : clusterPtrs ){
    LogDebug("PFClustering") 
      << "Loading PFCluster i="<<cluster.key()
      <<" energy="<<cluster->energy()<<std::endl;
    
    double Ecorr = _pfEnergyCalibration->energyEm(*cluster,
                                                  0.0,0.0,
                                                  applyCrackCorrections_);
    CalibratedClusterPtr calib_cluster(new CalibratedPFCluster(cluster,Ecorr));
    switch( cluster->layer() ) {
    case PFLayer::ECAL_BARREL:
      if( calib_cluster->energy() > threshPFClusterBarrel_ ) {
        _clustersEB.push_back(calib_cluster);   
      }
      break;
    case PFLayer::ECAL_ENDCAP:
      if( calib_cluster->energy() > threshPFClusterEndcap_ ) {
        _clustersEE.push_back(calib_cluster);
        _psclustersforee.emplace(calib_cluster->the_ptr(),
                                 edm::PtrVector<reco::PFCluster>());
      }
      break;
    default:
      break;
    }
  }
  // make the association map of ECAL clusters to preshower clusters  
  edm::PtrVector<reco::PFCluster> clusterPtrsPS = psclusters.ptrVector();
  double dist = -1.0, min_dist = -1.0;
  // match PS clusters to EE clusters, minimum distance to EE is ensured
  // since the inner loop is over the EE clusters
  for( const auto& psclus : clusterPtrsPS ) {   
    if( psclus->energy() < threshPFClusterES_ ) continue;        
    switch( psclus->layer() ) { // just in case this isn't the ES...
    case PFLayer::PS1:
    case PFLayer::PS2:
      break;
    default:
      continue;
    }    
    edm::Ptr<reco::PFCluster> eematch;
    dist = min_dist = -1.0; // reset
    for( const auto& eeclus : _clustersEE ) {
      dist = testPreshowerDistance(eeclus->the_ptr(),psclus);      
      if( dist == -1.0 || (min_dist != -1.0 && dist > min_dist) ) continue;
      if( dist < min_dist || min_dist == -1.0 ) {
        eematch = eeclus->the_ptr();
        min_dist = dist;
      }
    } // loop on EE clusters      
    if( eematch.isNonnull() ) _psclustersforee[eematch].push_back(psclus);
  } // loop on PS clusters

  // sort full cluster collections by their calibrated energy
  // this will put all the seeds first by construction
  GreaterByEt greater;
  std::sort(_clustersEB.begin(), _clustersEB.end(), greater);
  std::sort(_clustersEE.begin(), _clustersEE.end(), greater);  
}

void PFECALSuperClusterAlgo::run() {  
  // clusterize the EB
  buildAllSuperClusters(_clustersEB, threshPFClusterSeedBarrel_);
  // clusterize the EE
  buildAllSuperClusters(_clustersEE, threshPFClusterSeedEndcap_);
}

void PFECALSuperClusterAlgo::
buildAllSuperClusters(CalibClusterPtrVector& clusters,
                      double seedthresh) {
  IsASeed seedable(seedthresh);
  // make sure only seeds appear at the front of the list of clusters
  std::stable_partition(clusters.begin(),clusters.end(),seedable);
  // in each iteration we are working on a list that is already sorted
  // in the cluster energy and remains so through each iteration
  // NB: since clusters is sorted in loadClusters any_of has O(1)
  //     timing until you run out of seeds!  
  while( std::any_of(clusters.cbegin(), clusters.cend(), seedable) ) {    
    buildSuperCluster(clusters.front(),clusters);
  }
}

void PFECALSuperClusterAlgo::
buildSuperCluster(CalibClusterPtr& seed,
                  CalibClusterPtrVector& clusters) {
  IsClustered IsClusteredWithSeed(seed,_clustype,_useDynamicDPhi);
  IsLinkedByRecHit MatchesSeedByRecHit(seed,satelliteThreshold_,
                                       fractionForMajority_,0.1,0.2);
  bool isEE = false;
  SumPSEnergy sumps1(PFLayer::PS1), sumps2(PFLayer::PS2);  
  switch( seed->the_ptr()->layer() ) {
  case PFLayer::ECAL_BARREL:
    IsClusteredWithSeed.phiwidthSuperCluster_ = phiwidthSuperClusterBarrel_;
    IsClusteredWithSeed.etawidthSuperCluster_ = etawidthSuperClusterBarrel_;
    edm::LogInfo("PFClustering") << "Building SC number "  
                                 << superClustersEB_->size() + 1
                                 << " in the ECAL barrel!";
    break;
  case PFLayer::ECAL_ENDCAP:    
    IsClusteredWithSeed.phiwidthSuperCluster_ = phiwidthSuperClusterEndcap_; 
    IsClusteredWithSeed.etawidthSuperCluster_ = etawidthSuperClusterEndcap_;
    edm::LogInfo("PFClustering") << "Building SC number "  
                                 << superClustersEE_->size() + 1
                                 << " in the ECAL endcap!" << std::endl;
    isEE = true;
    break;
  default:
    break;
  }
  
  // this function shuffles the list of clusters into a list
  // where all clustered sub-clusters are at the front 
  // and returns a pointer to the first unclustered cluster.
  // The relative ordering of clusters is preserved 
  // (i.e. both resulting sub-lists are sorted by energy).
  auto not_clustered = std::stable_partition(clusters.begin(),clusters.end(),
                                             IsClusteredWithSeed);
  // satellite cluster merging
  // it was found that large clusters can split!
  if( doSatelliteClusterMerge_ ) {    
    not_clustered = std::stable_partition(not_clustered,clusters.end(),
                                          MatchesSeedByRecHit);
  }

  if(verbose_) {
    edm::LogInfo("PFClustering") << "Dumping cluster detail";
    edm::LogVerbatim("PFClustering")
      << "\tPassed seed: e = " << seed->energy_nocalib() 
      << " eta = " << seed->eta() << " phi = " << seed->phi() 
      << std::endl;  
    for( auto clus = clusters.cbegin(); clus != not_clustered; ++clus ) {
      edm::LogVerbatim("PFClustering") 
        << "\t\tClustered cluster: e = " << (*clus)->energy_nocalib() 
        << " eta = " << (*clus)->eta() << " phi = " << (*clus)->phi() 
        << std::endl;
    }
    for( auto clus = not_clustered; clus != clusters.end(); ++clus ) {
      edm::LogVerbatim("PFClustering") 
        << "\tNon-Clustered cluster: e = " << (*clus)->energy_nocalib() 
        << " eta = " << (*clus)->eta() << " phi = " << (*clus)->phi() 
        << std::endl;
    }    
  }
  // move the clustered clusters out of available cluster list
  // and into a temporary vector for building the SC  
  CalibratedClusterPtrVector clustered(clusters.begin(),not_clustered);
  clusters.erase(clusters.begin(),not_clustered);    
  // need the vector of raw pointers for a PF width class
  std::vector<const reco::PFCluster*> bare_ptrs;
  // calculate necessary parameters and build the SC
  double posX(0), posY(0), posZ(0),
    rawSCEnergy(0), corrSCEnergy(0), clusterCorrEE(0), 
    PS1_clus_sum(0), PS2_clus_sum(0);  
  for( auto& clus : clustered ) {
    bare_ptrs.push_back(clus->the_ptr().get());
      
    const double cluseraw = clus->energy_nocalib();
    const math::XYZPoint& cluspos = clus->the_ptr()->position();
    posX += cluseraw * cluspos.X();
    posY += cluseraw * cluspos.Y();
    posZ += cluseraw * cluspos.Z();
    // update EE calibrated super cluster energies
    if( isEE ) {
      const auto& psclusters = _psclustersforee[clus->the_ptr()];
      PS1_clus_sum = std::accumulate(psclusters.begin(),psclusters.end(),
                                     0.0,sumps1);
      PS2_clus_sum = std::accumulate(psclusters.begin(),psclusters.end(),
                                     0.0,sumps2);
      clusterCorrEE = 
        _pfEnergyCalibration->energyEm(*(clus->the_ptr()),
                                       PS1_clus_sum,PS2_clus_sum,
                                       applyCrackCorrections_);
      clus->resetCalibratedEnergy(clusterCorrEE);
    }

    rawSCEnergy  += cluseraw;
    corrSCEnergy += clus->energy();    
  }
  posX /= rawSCEnergy;
  posY /= rawSCEnergy;
  posZ /= rawSCEnergy;    
  
  // now build the supercluster
  reco::SuperCluster new_sc(corrSCEnergy,math::XYZPoint(posX,posY,posZ)); 
  double ps1_energy(0.0), ps2_energy(0.0), ps_energy(0.0);
  new_sc.setSeed(clustered.front()->the_ptr());
  for( const auto& clus : clustered ) {
    new_sc.addCluster(clus->the_ptr());
    auto& hits_and_fractions = clus->the_ptr()->hitsAndFractions();
    for( auto& hit_and_fraction : hits_and_fractions ) {
      new_sc.addHitAndFraction(hit_and_fraction.first,hit_and_fraction.second);
    }
    const auto& cluspsassociation = _psclustersforee[clus->the_ptr()];     
    // EE rechits should be uniquely matched to sets of pre-shower
    // clusters at this point, so we throw an exception if otherwise
    for( const auto& psclus : cluspsassociation ) {
      auto found_pscluster = std::find(new_sc.preshowerClustersBegin(),
                                       new_sc.preshowerClustersEnd(),
                                       reco::CaloClusterPtr(psclus));
      if( found_pscluster == new_sc.preshowerClustersEnd() ) {
        const double psenergy = psclus->energy();
        new_sc.addPreshowerCluster(psclus);
        ps1_energy += (PFLayer::PS1 == psclus->layer())*psenergy;
        ps2_energy += (PFLayer::PS2 == psclus->layer())*psenergy;
        ps_energy  += psenergy;
      } else {
        throw cms::Exception("PFECALSuperClusterAlgo::buildSuperCluster")
          << "Found a PS cluster matched to more than one EE cluster!" 
          << std::endl << std::hex << psclus.get() << " == " 
          << found_pscluster->get() << std::dec << std::endl;
      }
    }
  }
  new_sc.setPreshowerEnergy(ps_energy); 
  new_sc.setPreshowerEnergyPlane1(ps1_energy);
  new_sc.setPreshowerEnergyPlane2(ps2_energy);
  
  // calculate linearly weighted cluster widths
  PFClusterWidthAlgo pfwidth(bare_ptrs);
  new_sc.setEtaWidth(pfwidth.pflowEtaWidth());
  new_sc.setPhiWidth(pfwidth.pflowPhiWidth());
  
  // cache the value of the raw energy  
  new_sc.rawEnergy();

  // save the super cluster to the appropriate list
  switch( seed->the_ptr()->layer() ) {
  case PFLayer::ECAL_BARREL:
    superClustersEB_->push_back(new_sc);
    break;
  case PFLayer::ECAL_ENDCAP:    
    superClustersEE_->push_back(new_sc);
    break;
  default:
    break;
  }
}