Skip to content

Commit

Permalink
Merge pull request #15294 from fwyzard/backport_14058+14151_80x
Browse files Browse the repository at this point in the history 
backport of #14058: speed up avoiding eta phi computations, and #14151: Slim down PFRecHit
  • Loading branch information
@davidlange6
davidlange6 committed Aug 26, 2016
2 parents 86b00f3 + 299d487 commit ced6db6
Show file tree
Hide file tree
Showing 49 changed files with 676 additions and 922 deletions.
2 changes: 1 addition & 1 deletion CondTools/Geometry/plugins/calowriters.cc
View file
Expand Up @@ -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 ) ;
Expand Down
19 changes: 19 additions & 0 deletions DataFormats/Math/interface/PtEtaPhiMass.h
View file
Expand Up @@ -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
2 changes: 1 addition & 1 deletion DataFormats/ParticleFlowReco/interface/PFCluster.h
View file
Expand Up @@ -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() {
Expand Down
175 changes: 66 additions & 109 deletions DataFormats/ParticleFlowReco/interface/PFRecHit.h
View file
Expand Up @@ -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 {

Expand All @@ -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;
Expand All @@ -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_); }

Expand All @@ -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
42 changes: 42 additions & 0 deletions DataFormats/ParticleFlowReco/src/PFCluster.cc
View file
@@ -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;

Expand Down

0 comments on commit ced6db6

Please sign in to comment.