PFClusterAlgo.h

Go to the documentation of this file.

#ifndef RecoParticleFlow_PFClusterProducer_PFClusterAlgo_h
#define RecoParticleFlow_PFClusterProducer_PFClusterAlgo_h


#include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
#include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h"
#include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
#include "DataFormats/Common/interface/OrphanHandle.h"

// for E/gamma position calc
#include "FWCore/Framework/interface/ESHandle.h"
#include "RecoEcal/EgammaCoreTools/interface/PositionCalc.h"

#include <string>
#include <vector>
#include <map>
#include <set>

#include <memory>

class TFile;
class TH2F;


namespace edm {
  template<> 
  struct StrictWeakOrdering<reco::PFRecHit> {   
    typedef unsigned key_type;
    bool operator()(unsigned a, reco::PFRecHit const& b) const { 
      return a < b.detId(); 
    }
    bool operator()(reco::PFRecHit const& a, unsigned b) const { 
      return a.detId() < b; 
    }
    bool operator()(reco::PFRecHit const& a, reco::PFRecHit const& b) const { 
      return ( a.detId() < b.detId()  ); 
    }
  };
}

class PFClusterAlgo {

 public:

  typedef edm::StrictWeakOrdering<reco::PFRecHit> PFStrictWeakOrdering;
  typedef edm::SortedCollection<reco::PFRecHit> SortedPFRecHitCollection;
  typedef edm::Handle<edm::View<reco::PFCluster> > PFClusterHandle;

  enum PositionCalcType { EGPositionCalc,
              EGPositionFormula,              
              PFPositionCalc,
              kNotDefined };

  PFClusterAlgo();

  virtual ~PFClusterAlgo() {;}

  // void init(const std::map<unsigned, reco::PFRecHit* >& rechits );

  void enableDebugging(bool debug) { debug_ = debug;}
 

  typedef edm::Handle< reco::PFRecHitCollection > PFRecHitHandle;
  
  void doClustering( const reco::PFRecHitCollection& rechits );
  void doClustering( const reco::PFRecHitCollection& rechits, 
             const std::vector<bool> & mask );

  void doClustering( const PFRecHitHandle& rechitsHandle );  
  void doClustering( const PFRecHitHandle& rechitsHandle, 
             const std::vector<bool> & mask );
    
  
  void setThreshBarrel(double thresh) {threshBarrel_ = thresh;}
  void setThreshPtBarrel(double thresh) {threshPtBarrel_ = thresh;}

  void setThreshSeedBarrel(double thresh) {threshSeedBarrel_ = thresh;}
  void setThreshPtSeedBarrel(double thresh) {threshPtSeedBarrel_ = thresh;}

  void setThreshCleanBarrel(double thresh) {threshCleanBarrel_ = thresh;}
  void setS4S1CleanBarrel(const std::vector<double>& coeffs) {minS4S1Barrel_ = coeffs;}

  void setThreshDoubleSpikeBarrel( double thresh ) { threshDoubleSpikeBarrel_ = thresh;}
  void setS6S2DoubleSpikeBarrel( double cut ) { minS6S2DoubleSpikeBarrel_ = cut;}

  void setThreshEndcap(double thresh) {threshEndcap_ = thresh;}
  void setThreshPtEndcap(double thresh) {threshPtEndcap_ = thresh;}

  void setThreshSeedEndcap(double thresh) {threshSeedEndcap_ = thresh;}
  void setThreshPtSeedEndcap(double thresh) {threshPtSeedEndcap_ = thresh;}

  void setThreshCleanEndcap(double thresh) {threshCleanEndcap_ = thresh;}
  void setS4S1CleanEndcap(const std::vector<double>& coeffs) {minS4S1Endcap_ = coeffs;}

  void setThreshDoubleSpikeEndcap( double thresh ) { threshDoubleSpikeEndcap_ = thresh;}
  void setS6S2DoubleSpikeEndcap( double cut ) { minS6S2DoubleSpikeEndcap_ = cut;}

  void setHistos(TFile* file, TH2F* hB, TH2F* hE) {file_=file; hBNeighbour = hB; hENeighbour = hE;}

  void setNNeighbours(int n) { nNeighbours_ = n;}

  // set position calculation variant
  void setPositionCalcType( const PositionCalcType& t ) {which_pos_calc_ = t;}

  // setup E/gamma position calc
  void setEGammaPosCalc( const edm::ParameterSet& conf ) { 
    eg_pos_calc.reset(new PositionCalc(conf));
  }
  void setEBGeom(const CaloSubdetectorGeometry* esh) {
    eb_geom = esh;
  }
  void setEEGeom(const CaloSubdetectorGeometry* esh) {
    ee_geom = esh;
  }
  void setPreshowerGeom(const CaloSubdetectorGeometry* esh) {
    preshower_geom = esh;
  }

  // set W0 parameter for EG-type position formula
  void setPosCalcW0( double w0 ) { param_W0_ = w0; }

  void setPosCalcP1( double p1 ) { posCalcP1_ = p1; }

  void setPosCalcNCrystal(int n) { posCalcNCrystal_ = n;}

  void setShowerSigma( double sigma ) { showerSigma_ = sigma;}

  void setUseCornerCells( bool usecornercells ) { useCornerCells_ = usecornercells;}
  
  void setCleanRBXandHPDs( bool cleanRBXandHPDs) { cleanRBXandHPDs_ = cleanRBXandHPDs; }

 
  double threshBarrel() const {return threshBarrel_;}

  double threshSeedBarrel() const  {return threshSeedBarrel_;}


  double threshEndcap() const {return threshEndcap_;}

  double threshSeedEndcap() const {return threshSeedEndcap_;}


  int nNeighbours() const { return nNeighbours_;}

  double posCalcP1() const { return posCalcP1_; }

  int posCalcNCrystal() const {return  posCalcNCrystal_;}

  double showerSigma() const { return showerSigma_ ;}
  
  void write();



  const reco::PFRecHit& rechit(unsigned i,
                   const reco::PFRecHitCollection& rechits );

  bool masked(unsigned rhi) const;

  enum Color { 
    NONE=0,
    SEED, 
    SPECIAL 
  };

  unsigned color(unsigned rhi) const;

  bool isSeed(unsigned rhi) const;

  std::auto_ptr< std::vector< reco::PFCluster > >& clusters()  
    {return pfClusters_;}

  std::auto_ptr< std::vector< reco::PFRecHit > >& rechitsCleaned()  
    {return pfRecHitsCleaned_;}
  

  enum Parameter { THRESH, 
           SEED_THRESH,
           PT_THRESH,
           SEED_PT_THRESH,
                   CLEAN_THRESH,
                   CLEAN_S4S1,
                   DOUBLESPIKE_THRESH,
                   DOUBLESPIKE_S6S2
  };
  
    double parameter( Parameter paramtype, PFLayer::Layer layer, unsigned iCoeff = 0, int iring0=0) const; 

  
  enum SeedState {
    UNKNOWN=-1,
    NO=0,
    YES=1,
    CLEAN=2
  };


  friend std::ostream& operator<<(std::ostream& out,const PFClusterAlgo& algo);

  typedef std::map<unsigned, unsigned >::const_iterator IDH;
  typedef std::multimap<double, unsigned >::iterator EH;
 

 private:
  void doClusteringWorker( const reco::PFRecHitCollection& rechits );  

  void cleanRBXAndHPD( const reco::PFRecHitCollection& rechits );

  void findSeeds( const reco::PFRecHitCollection& rechits );

  void buildTopoClusters( const reco::PFRecHitCollection& rechits ); 

  void buildTopoCluster( std::vector< unsigned >& cluster, unsigned rhi, 
             const reco::PFRecHitCollection& rechits ); 
  
  void buildPFClusters( const std::vector< unsigned >& cluster, 
            const reco::PFRecHitCollection& rechits ); 

  void calculateClusterPosition( reco::PFCluster& cluster, 
                                 reco::PFCluster& clusterwodepthcor,
                 bool depcor = true,
                 int posCalcNCrystal=0);
  
  reco::PFRecHitRef  createRecHitRef( const reco::PFRecHitCollection& rechits, 
                      unsigned rhi );

  void paint( unsigned rhi, unsigned color=1 );

  std::pair<double,double> dCrack(double phi, double eta);  

  PFRecHitHandle           rechitsHandle_;   
  // for E\gamma Position calc
  std::auto_ptr<SortedPFRecHitCollection> sortedRecHits_; 

  std::set<unsigned>       idUsedRecHits_;

  std::multimap<double, unsigned, std::greater<double> >  eRecHits_;

  std::vector< bool >      mask_;

  std::vector< unsigned >  color_;

  std::vector< SeedState > seedStates_;

  std::vector< bool >      usedInTopo_;

  std::vector< unsigned >  seeds_; 

  std::vector< std::vector< unsigned > > topoClusters_;

  // std::vector< reco::PFCluster >  allClusters_;

  std::auto_ptr< std::vector<reco::PFCluster> > pfClusters_;
  
  std::auto_ptr< std::vector<reco::PFRecHit> > pfRecHitsCleaned_;  

  double threshBarrel_;
  double threshPtBarrel_;

  double threshSeedBarrel_;
  double threshPtSeedBarrel_;

  double threshEndcap_;
  double threshPtEndcap_;

  double threshSeedEndcap_;
  double threshPtSeedEndcap_;

  double threshCleanBarrel_;
  std::vector<double> minS4S1Barrel_;

  double threshDoubleSpikeBarrel_;
  double minS6S2DoubleSpikeBarrel_;

  double threshCleanEndcap_;
  std::vector<double> minS4S1Endcap_;

  double threshDoubleSpikeEndcap_;
  double minS6S2DoubleSpikeEndcap_;

  int    nNeighbours_;

  int    posCalcNCrystal_;

  PositionCalcType which_pos_calc_;
  double posCalcP1_;
  double param_W0_; // for EG position calc variant
  std::auto_ptr<PositionCalc> eg_pos_calc; //directly use E/gamma pos calc
  const CaloSubdetectorGeometry *eb_geom, *ee_geom, *preshower_geom;

  double showerSigma_;

  bool useCornerCells_;

  bool cleanRBXandHPDs_;

  bool   debug_;


  static  unsigned  prodNum_;

  // Histograms
  TH2F* hBNeighbour;
  TH2F* hENeighbour;
  TFile*     file_; 

};

#endif