Rewrite and Improvements for PFRecHit and PFCluster Producers

DataFormats/ParticleFlowReco/interface/PFCluster.h

@@ -9,10 +9,12 @@
 #include "Rtypes.h" 
 
 #include "DataFormats/ParticleFlowReco/interface/PFRecHitFraction.h"
+#include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
 #include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
 
 #include <iostream>
 #include <vector>
+#include <algorithm>
 #if !defined(__CINT__) && !defined(__MAKECINT__) && !defined(__REFLEX__)
 #include <atomic>
 #endif
@@ -48,30 +50,38 @@ namespace reco {
     typedef std::vector<std::pair<CaloClusterPtr::key_type,edm::Ptr<PFCluster> > > EEtoPSAssociation;
     typedef ROOT::Math::PositionVector3D<ROOT::Math::CylindricalEta3D<Double32_t> > REPPoint;
 
-    PFCluster() : CaloCluster(CaloCluster::particleFlow), color_(1) {}
+    PFCluster() : CaloCluster(CaloCluster::particleFlow), layer_(PFLayer::NONE), color_(1) {}
 
     /// constructor
     PFCluster(PFLayer::Layer layer, double energy,
 	      double x, double y, double z );
 
     /// resets clusters parameters
     void reset();
+
+    /// reset only hits and fractions
+    void resetHitsAndFractions();
 
     /// add a given fraction of the rechit
     void addRecHitFraction( const reco::PFRecHitFraction& frac);
 
     /// vector of rechit fractions
     const std::vector< reco::PFRecHitFraction >& recHitFractions() const 
-      { return rechits_; }
+      { return rechits_; }    
 
     /// set layer
     void setLayer( PFLayer::Layer layer);
 
     /// cluster layer, see PFLayer.h in this directory
-    PFLayer::Layer  layer() const;     
+    PFLayer::Layer  layer() const;
 
     /// cluster energy
     double        energy() const {return energy_;}
+
+    /// cluster time
+    double        time() const {return time_;}
+
+    void         setTime(double time) {time_ = time;}
 
     /// cluster position: rho, eta, phi
     const REPPoint&       positionREP() const {return posrep_;}
@@ -136,13 +146,35 @@ namespace reco {
     double vy() const { return vertex().y(); }
     double vz() const { return vertex().z(); }    
 
+#if !defined(__CINT__) && !defined(__MAKECINT__) && !defined(__REFLEX__)
+    template<typename pruner>
+      void pruneUsing(pruner prune) {
+      hitsAndFractions_.clear();
+      std::vector<reco::PFRecHitFraction>::iterator iter = 
+	std::stable_partition(rechits_.begin(),rechits_.end(),prune);
+      rechits_.erase(iter,rechits_.end());
+      hitsAndFractions_.reserve(rechits_.size());
+      for( const auto& hitfrac : rechits_ ) {
+	hitsAndFractions_.emplace_back(hitfrac.recHitRef()->detId(),
+				       hitfrac.fraction());
+      }
+    }
+#endif
+
   private:
 
     /// vector of rechit fractions (transient)
     std::vector< reco::PFRecHitFraction >  rechits_;
 
     /// cluster position: rho, eta, phi (transient)
     REPPoint            posrep_;
+
+    ///Michalis :Add timing information
+    double time_;
+
+    /// transient layer
+    PFLayer::Layer layer_; 
+
 #if !defined(__CINT__) && !defined(__MAKECINT__) && !defined(__REFLEX__)
      /// \todo move to PFClusterTools
     static std::atomic<int>    depthCorMode_;
@@ -159,7 +191,6 @@ namespace reco {
     /// \todo move to PFClusterTools
     static std::atomic<double> depthCorBp_;
 #else
-
     /// \todo move to PFClusterTools
     static int    depthCorMode_;
 
@@ -180,8 +211,6 @@ namespace reco {
 
     /// color (transient)
     int                 color_;
-
-    friend class ::PFClusterAlgo;
   };
 }
 

DataFormats/ParticleFlowReco/interface/PFRecHit.h

@@ -9,16 +9,18 @@
 #include "DataFormats/Math/interface/Point3D.h"
 #include "Rtypes.h" 
 #include "DataFormats/Math/interface/Vector3D.h"
-// #include "DataFormats/DetId/interface/DetId.h"
 #include "Math/GenVector/PositionVector3D.h"
-
 #include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
+#include "DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h"
 
 //C decide what is the default rechit index. 
 //C maybe 0 ? -> compression 
 //C then the position is index-1. 
 //C provide a helper class to access the rechit. 
 
+#include "DataFormats/CaloRecHit/interface/CaloRecHit.h"
+#include "DataFormats/Common/interface/RefToBase.h"
+
 
 namespace reco {
 
@@ -27,19 +29,20 @@ namespace reco {
 
      \author Colin Bernet
      \date   July 2006
+
+     Feb 2014 [Michalis: 8 years later!Modifying the class to be able to generalize the neighbours for 3D calorimeters ]
   */
   class PFRecHit {
 
   public:
-
-    enum {
-      NONE=0
-    };
 
     typedef ROOT::Math::PositionVector3D<ROOT::Math::CylindricalEta3D<Double32_t> > REPPoint;
 
     typedef std::vector<REPPoint> REPPointVector;
-
+
+    enum {
+      NONE=0
+    };
     /// default constructor. Sets energy and position to zero
     PFRecHit();
 
@@ -63,38 +66,31 @@ namespace reco {
     /// destructor
     virtual ~PFRecHit();
 
-
-    /// NICHOLAS JIN
-    /// Zeros the energy
-
     void setEnergy( double energy) { energy_ = energy; }
 
-    /// calculates rho eta phi position once and for all
     void calculatePositionREP();
 
-    //C neighbours must be initialized correctly !!
-    /*     void setNeighbours( const std::vector< unsigned >& neighbours ); */
-    void add4Neighbour( unsigned index );
-    void add8Neighbour( unsigned index );
-    void setEnergyUp( double eUp) { energyUp_ = eUp; }
-    void setRescale( double factor) { rescale_ = factor; }
-
+    void addNeighbour(short x,short y, short z,const PFRecHitRef&);
+    void setTime( double time) { time_ = time; }
     void clearNeighbours() {
-      neighbours4_.clear();
-      neighbours8_.clear();
-      neighbours4_.reserve(4);
-      neighbours8_.reserve(8);    
+      neighbours_.clear();
     }
 
-
+    const PFRecHitRefVector& neighbours4() const {
+      return neighbours4_;
+    }
+    const PFRecHitRefVector& neighbours8() const {
+      return neighbours8_;
+    }
+
+    const PFRecHitRefVector& neighbours() const {
+      return neighbours_;
+    }
+
+    const std::vector<unsigned short>& neighbourInfos() {
+      return neighbourInfos_;
+    }
 
-    /// \brief search for pointers to neighbours, using neighbours' DetId.
-    /// 
-    /// pointers to neighbours are not persistent, in contrary to the DetId's 
-    /// of the neighbours. This function searches a map of rechits 
-    /// for the DetId's stored in neighboursIds4_ and  neighboursIds8_. 
-    /// The corresponding pointers are stored in neighbours4_ and neighbours8_.
-    // void      findPtrsToNeighbours( const std::map<unsigned,  reco::PFRecHit* >& allhits );
 
     void      setNWCorner( double posx, double posy, double posz );
     void      setSWCorner( double posx, double posy, double posz );
@@ -110,14 +106,9 @@ namespace reco {
     /// rechit energy
     double energy() const { return energy_; }
 
-    /// rescaling factor (obsolete)
-    double rescale() const { return rescale_; }
 
     /// timing for cleaned hits
-    double time() const { return rescale_; }
-
-    /// For HF hits: rechit energy (and neighbour's) in the other HF layer
-    double energyUp() const { return energyUp_; }
+    double time() const { return time_; }
 
     /// rechit momentum transverse to the beam, squared.
     double pt2() const { return energy_ * energy_ *
@@ -127,21 +118,12 @@ namespace reco {
 			     position_.Y()*position_.Y() + 
 			     position_.Z()*position_.Z()) ; }
 
-    //C remove cause I want to be able to run on const rechits
-    /// \return seed state (-1:unknown, 0:no, 1 yes)
-    // int  seedState() const { return seedState_; }
-
-    /// is seed ? 
-    // bool isSeed() const { return (seedState_>0) ? true : false; }
-
-    /// set seed status
-    // void youAreSeed(int seedstate=1) {seedState_ = seedstate;} 
 
     /// rechit cell centre x, y, z
     const math::XYZPoint& position() const { return position_; }
 
-    /// rechit cell centre rho, eta, phi. call calculatePositionREP before !
-    const REPPoint& positionREP() const;
+    const REPPoint& positionREP() const { return positionrep_; }
+
 
     /// rechit cell axis x, y, z
     const math::XYZVector& getAxisXYZ() const { return axisxyz_; }    
@@ -150,36 +132,7 @@ namespace reco {
     const std::vector< math::XYZPoint >& getCornersXYZ() const 
       { return cornersxyz_; }    
 
-    /// rechit corners
-    const REPPointVector& getCornersREP() const 
-      { return cornersrep_; }    
-
-    const std::vector< unsigned >& neighbours4() const 
-      {return neighbours4_;}  
-
-    const std::vector< unsigned >& neighbours8() const 
-      {return neighbours8_;}  
-
-    const std::vector< unsigned >& neighboursIds4() const 
-      {return neighboursIds4_;}  
-
-    const std::vector< unsigned >& neighboursIds8() const 
-      {return neighboursIds8_;}  
-
-    /*     const std::vector< unsigned >& getNeighboursIds4() const  */
-    /*       {return neighboursIds4_;}   */
-
-    /*     const std::vector< unsigned >& getNeighboursIds8() const  */
-    /*       {return neighboursIds8_;}   */
-
-    /// is rechit 'id' a direct neighbour of this ? 
-    /// id is the rechit index ! not the detId
-    bool  isNeighbour4(unsigned id) const;
-
-    /// is rechit 'id' a neighbour of this ? 
-    /// id is the rechit index ! not the detId
-    bool  isNeighbour8(unsigned id) const;
-
+    const std::vector<REPPoint>& getCornersREP() const { return cornersrep_; }
 
     void size(double& deta, double& dphi) const;
 
@@ -198,54 +151,56 @@ namespace reco {
     friend std::ostream& operator<<(std::ostream& out, 
                                     const reco::PFRecHit& hit);
 
+    const edm::RefToBase<CaloRecHit>& originalRecHit() const {
+      return originalRecHit_;
+    }
+
+    template<typename T> 
+    void setOriginalRecHit(const T& rh) {
+      originalRecHit_ = edm::RefToBase<CaloRecHit>(rh);
+    }
+
   private:
 
+    // original rechit
+    edm::RefToBase<CaloRecHit> originalRecHit_;
+
     ///C cell detid - should be detid or index in collection ?
     unsigned            detId_;             
 
     /// rechit layer
-    PFLayer::Layer                 layer_;
+    PFLayer::Layer      layer_;
 
     /// rechit energy 
     double              energy_;
 
-    /// Internal rescaling factor of the energy (1. = default, 0 = killed channels, x = rescaled)
-    double              rescale_;
-
-    /// For HF hits : hit energy in the other layer (EM for HAD, and HAD for EM)
-    double              energyUp_;
+    /// time
+    double              time_;
 
-    /// is this a seed ? (-1:unknown, 0:no, 1 yes) (transient)
-    // int                 seedState_;
-
     /// rechit cell centre: x, y, z
     math::XYZPoint      position_;
 
     /// rechit cell centre: rho, eta, phi (transient)
-    REPPoint            posrep_;
+    REPPoint positionrep_;
 
     /// rechit cell axisxyz
     math::XYZVector     axisxyz_;
 
     /// rechit cell corners
     std::vector< math::XYZPoint > cornersxyz_;
-    REPPointVector cornersrep_;
+
+    /// rechit cell corners rho/eta/phi
+    std::vector< REPPoint > cornersrep_;
 
     /// indices to existing neighbours (1 common side)
-    std::vector< unsigned >   neighbours4_;
-
-    /// indices to existing neighbours (1 common side or diagonal) 
-    std::vector< unsigned >   neighbours8_;
+    PFRecHitRefVector   neighbours_;
+    std::vector< unsigned short >   neighbourInfos_;
 
-    /// detids of existing neighbours (1 common side)
-    std::vector< unsigned >   neighboursIds4_;
+    //Caching the neighbours4/8 per request of Lindsey
+    PFRecHitRefVector   neighbours4_;
+    PFRecHitRefVector   neighbours8_;
 
-    /// detids of existing neighbours (1 common side or diagonal) 
-    std::vector< unsigned >   neighboursIds8_;
 
-    /// number of neighbours
-    static const unsigned    nNeighbours_;
-
     /// number of corners
     static const unsigned    nCorners_;
 

DataFormats/ParticleFlowReco/interface/PFRecHitFwd.h

@@ -22,6 +22,9 @@ namespace reco {
 
   /// iterator over a vector of references to PFRecHit objects
   typedef PFRecHitRefVector::iterator basicRecHit_iterator;
+
+  /// ref to base vector for dealing with views
+  typedef edm::RefToBaseVector<reco::PFRecHit> PFRecHitBaseRefVector;
 }
 
 #endif