backport of #14058: speed up avoiding eta phi computations, and #14151: Slim down PFRecHit

CondTools/Geometry/plugins/calowriters.cc

@@ -75,7 +75,7 @@ CaloGeometryDBEP<HcalGeometry, CaloGeometryDBWriter>::produceAligned( const type
 
     ptr->fillDefaultNamedParameters() ;
 
-    ptr->allocateCorners( hcalTopology->ncells() ) ;
+    ptr->allocateCorners( hcalTopology->ncells() + hcalTopology->getHFSize() ) ;
 
     ptr->allocatePar(    dvec.size() ,
 			 HcalGeometry::k_NumberOfParametersPerShape ) ;

DataFormats/Math/interface/PtEtaPhiMass.h

@@ -34,5 +34,24 @@ class PtEtaPhiMass {
 };
 
 
+class RhoEtaPhi {
+private:
+  float rho_, eta_, phi_;
+public:
+  // default constructor (unitialized)
+  RhoEtaPhi() {}
+
+  //positional constructor (still compatible with Root, c++03)
+  RhoEtaPhi(float irho, float ieta, float iphi):
+    rho_(irho), eta_(ieta), phi_(iphi) {}
+
+  /// transverse momentum
+  float rho() const { return rho_;}
+  /// momentum pseudorapidity
+  float eta() const { return eta_; }
+  /// momentum azimuthal angle
+  float phi() const { return phi_; }
+};
+
 
 #endif

DataFormats/ParticleFlowReco/interface/PFCluster.h

@@ -90,7 +90,7 @@ namespace reco {
     void         setDepth(double depth) {depth_ = depth;}
 
     /// cluster position: rho, eta, phi
-    const REPPoint&       positionREP() const {return posrep_;}
+    const REPPoint& positionREP() const {return posrep_;}
 
     /// computes posrep_ once and for all
     void calculatePositionREP() {

DataFormats/ParticleFlowReco/interface/PFRecHit.h

@@ -7,20 +7,16 @@
 #include <iostream>
 
 #include "DataFormats/Math/interface/Point3D.h"
-#include "Rtypes.h" 
 #include "DataFormats/Math/interface/Vector3D.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"
 
+#include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
+
 
 namespace reco {
 
@@ -35,35 +31,36 @@ namespace reco {
   class PFRecHit {
 
   public:
-
-    // Next typedef uses double in ROOT 6 rather than Double32_t due to a bug in ROOT 5,
-    // which otherwise would make ROOT5 files unreadable in ROOT6.  This does not increase
-    // the size on disk, because due to the bug, double was actually stored on disk in ROOT 5.
-
-    typedef ROOT::Math::PositionVector3D<ROOT::Math::CylindricalEta3D<double> > REPPoint;
-
-    typedef std::vector<REPPoint> REPPointVector;
-
+    using PositionType = GlobalPoint::BasicVectorType;
+    using REPPoint = RhoEtaPhi;
+    using RepCorners = CaloCellGeometry::RepCorners;
+    using REPPointVector = RepCorners;
+    using CornersVec = CaloCellGeometry::CornersVec;
+
+    struct Neighbours {
+      using Pointer = unsigned int const *;
+      Neighbours(){}
+      Neighbours(Pointer ib, unsigned int n) : b(ib), e(ib+n){}
+      Pointer b, e;
+      Pointer begin() const {return b;}
+      Pointer end() const {return e;}
+      unsigned int size() const { return e-b;}
+    };
+
     enum {
       NONE=0
     };
     /// default constructor. Sets energy and position to zero
-    PFRecHit();
+    PFRecHit(){}
 
-    /// constructor from values
-    PFRecHit(unsigned detId,
+    PFRecHit(CaloCellGeometry const * caloCell, unsigned int detId,
              PFLayer::Layer layer,
-             double energy, 
-             const math::XYZPoint& posxyz, 
-             const math::XYZVector& axisxyz, 
-             const std::vector< math::XYZPoint >& cornersxyz);
+             float energy) :
+        caloCell_(caloCell),  detId_(detId),
+        layer_(layer), energy_(energy){}
 
-    PFRecHit(unsigned detId,
-             PFLayer::Layer layer,
-             double energy, 
-             double posx, double posy, double posz, 
-             double axisx, double axisy, double axisz);    
 
+
     /// copy
     PFRecHit(const PFRecHit& other) = default;
     PFRecHit(PFRecHit&& other) = default;
@@ -74,81 +71,71 @@ namespace reco {
     /// destructor
     ~PFRecHit()=default;
 
-    void setEnergy( double energy) { energy_ = energy; }
+    void setEnergy( float energy) { energy_ = energy; }
 
-    void calculatePositionREP();
 
-    void addNeighbour(short x,short y, short z,const PFRecHitRef&);
-    const PFRecHitRef getNeighbour(short x,short y, short z);
+    void addNeighbour(short x,short y, short z, unsigned int);
+    unsigned int getNeighbour(short x,short y, short z) const;
     void setTime( double time) { time_ = time; }
     void setDepth( int depth) { depth_ = depth; }
     void clearNeighbours() {
       neighbours_.clear();
+      neighbourInfos_.clear();
+      neighbours4_ = neighbours8_ = 0;
     }
 
-    const PFRecHitRefVector& neighbours4() const {
-      return neighbours4_;
+    Neighbours neighbours4() const {
+      return buildNeighbours(neighbours4_);
     }
-    const PFRecHitRefVector& neighbours8() const {
-      return neighbours8_;
+    Neighbours neighbours8() const {
+	return buildNeighbours(neighbours8_);
     }
 
-    const PFRecHitRefVector& neighbours() const {
-      return neighbours_;
+    Neighbours neighbours() const {
+      return buildNeighbours(neighbours_.size());
     }
 
     const std::vector<unsigned short>& neighbourInfos() {
       return neighbourInfos_;
     }
 
 
-    void      setNWCorner( double posx, double posy, double posz );
-    void      setSWCorner( double posx, double posy, double posz );
-    void      setSECorner( double posx, double posy, double posz );
-    void      setNECorner( double posx, double posy, double posz );
-
+    /// calo cell
+    CaloCellGeometry const & caloCell() const { return  *caloCell_; }
+    bool hasCaloCell() const { return caloCell_; }
+
     /// rechit detId
     unsigned detId() const {return detId_;}
 
     /// rechit layer
     PFLayer::Layer layer() const { return layer_; }
 
     /// rechit energy
-    double energy() const { return energy_; }
+    float energy() const { return energy_; }
 
 
     /// timing for cleaned hits
-    double time() const { return time_; }
+    float time() const { return time_; }
 
     /// depth for segemntation
     int  depth() const { return depth_; }
 
     /// rechit momentum transverse to the beam, squared.
     double pt2() const { return energy_ * energy_ *
-			   ( position_.X()*position_.X() + 
-			     position_.Y()*position_.Y() ) / 
-			   ( position_.X()*position_.X() +
-			     position_.Y()*position_.Y() + 
-			     position_.Z()*position_.Z()) ; }
+	( position().perp2()/ position().mag2());}
 
 
     /// rechit cell centre x, y, z
-    const math::XYZPoint& position() const { return position_; }
-
-    const REPPoint& positionREP() const { return positionrep_; }
-
-
-    /// rechit cell axis x, y, z
-    const math::XYZVector& getAxisXYZ() const { return axisxyz_; }    
+    PositionType const & position() const { return caloCell().getPosition().basicVector(); }
+
+    RhoEtaPhi const &  positionREP() const { return caloCell().repPos(); }
 
     /// rechit corners
-    const std::vector< math::XYZPoint >& getCornersXYZ() const 
-      { return cornersxyz_; }    
-
-    const std::vector<REPPoint>& getCornersREP() const { return cornersrep_; }
-
-    void size(double& deta, double& dphi) const;
+    CornersVec const & getCornersXYZ() const { return caloCell().getCorners(); }    
 
+    RepCorners const & getCornersREP() const { return caloCell().getCornersREP();}
+
+
     /// comparison >= operator
     bool operator>=(const PFRecHit& rhs) const { return (energy_>=rhs.energy_); }
 
@@ -161,70 +148,40 @@ namespace reco {
     /// comparison < operator
     bool operator< (const PFRecHit& rhs) const { return (energy_< rhs.energy_); }
 
-    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_;             
+    Neighbours buildNeighbours(unsigned int n) const { return  Neighbours(&neighbours_.front(),n);}
+
+    /// cell geometry
+    CaloCellGeometry const * caloCell_=nullptr;
+
+    ///cell detid
+    unsigned  int        detId_=0;             
 
     /// rechit layer
-    PFLayer::Layer      layer_;
+    PFLayer::Layer      layer_=PFLayer::NONE;
 
     /// rechit energy 
-    double              energy_;
+    float              energy_=0;
 
     /// time
-    double              time_;
-
+    float              time_=-1;
 
     /// depth
-    int      depth_;
-
-    /// rechit cell centre: x, y, z
-    math::XYZPoint      position_;
-
-    /// rechit cell centre: rho, eta, phi (transient)
-    REPPoint positionrep_;
+    int      depth_=0;
 
-    /// rechit cell axisxyz
-    math::XYZVector     axisxyz_;
-
-    /// rechit cell corners
-    std::vector< math::XYZPoint > cornersxyz_;
-
-    /// rechit cell corners rho/eta/phi
-    std::vector< REPPoint > cornersrep_;
 
     /// indices to existing neighbours (1 common side)
-    PFRecHitRefVector   neighbours_;
+    std::vector< unsigned int > neighbours_;
     std::vector< unsigned short >   neighbourInfos_;
 
     //Caching the neighbours4/8 per request of Lindsey
-    PFRecHitRefVector   neighbours4_;
-    PFRecHitRefVector   neighbours8_;
-
-
-    /// number of corners
-    static const unsigned    nCorners_;
-
-    /// set position of one of the corners
-    void      setCorner( unsigned i, double posx, double posy, double posz );
+    unsigned int neighbours4_ = 0;
+    unsigned int neighbours8_ = 0;
   };
-  
+
 }
+std::ostream& operator<<(std::ostream& out, const reco::PFRecHit& hit);
 
 #endif

DataFormats/ParticleFlowReco/src/PFCluster.cc

@@ -1,5 +1,47 @@
 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
 
+#include "vdt/vdtMath.h"
+#include "Math/GenVector/etaMax.h"
+
+namespace {
+
+  // an implementation of Eta_FromRhoZ of root libraries using vdt
+  template<typename Scalar>
+  inline Scalar Eta_FromRhoZ_fast(Scalar rho, Scalar z) {    
+    using namespace ROOT::Math;
+    // value to control Taylor expansion of sqrt
+    const Scalar big_z_scaled =
+      std::pow(std::numeric_limits<Scalar>::epsilon(),static_cast<Scalar>(-.25));
+    if (rho > 0) {      
+      Scalar z_scaled = z/rho;
+      if (std::fabs(z_scaled) < big_z_scaled) {
+        return vdt::fast_log(z_scaled+std::sqrt(z_scaled*z_scaled+1.0));
+      } else {
+        // apply correction using first order Taylor expansion of sqrt
+        return  z>0 ? vdt::fast_log(2.0*z_scaled + 0.5/z_scaled) : -vdt::fast_log(-2.0*z_scaled);
+      }
+    }
+    // case vector has rho = 0
+    else if (z==0) {
+      return 0;
+    }
+    else if (z>0) {
+      return z + etaMax<Scalar>();
+    }
+    else {
+      return z - etaMax<Scalar>();
+    }    
+  }
+
+  inline void calculateREP(const math::XYZPoint& pos, double& rho, double& eta, double& phi) {
+    const double z = pos.z();
+    rho = pos.Rho();    
+    eta = Eta_FromRhoZ_fast<double>(rho,z);
+    phi = (pos.x()==0 && pos.y()==0) ? 0 : vdt::fast_atan2(pos.y(), pos.x());
+  }
+
+}
+
 using namespace std;
 using namespace reco;