Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Speedups to HGCal clustering #16997

Merged
merged 3 commits into from Dec 22, 2016
Merged

Speedups to HGCal clustering #16997

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
6e45d59
speedups to HGCal clustering
lgray Dec 13, 2016
cfeda61
use deltar2 in multiclusters as well since there is only ever a compa…
lgray Dec 14, 2016
f921ac6
code review
lgray Dec 20, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
9 changes: 8 additions & 1 deletion RecoLocalCalo/HGCalRecAlgos/interface/HGCalImagingAlgo.h
View file
Open in desktop
Expand Up @@ -194,7 +194,14 @@ class HGCalImagingAlgo
//@@EM todo: the number of layers should be obtained programmatically - the range is 1-n instead of 0-n-1...

//these functions should be in a helper class.
double distance(const Hexel &pt1, const Hexel &pt2); //2-d distance on the layer (x-y)
inline double distance2(const Hexel &pt1, const Hexel &pt2) { //distance squared
const double dx = pt1.x - pt2.x;
const double dy = pt1.y - pt2.y;
return (dx*dx + dy*dy);
} //distance squaredq
inline double distance(const Hexel &pt1, const Hexel &pt2) { //2-d distance on the layer (x-y)
return std::sqrt(distance2(pt1,pt2));
}
double calculateLocalDensity(std::vector<KDNode> &, KDTree &); //return max density
double calculateDistanceToHigher(std::vector<KDNode> &, KDTree &);
int findAndAssignClusters(std::vector<KDNode> &, KDTree &, double, KDTreeBox &);
Expand Down
12 changes: 8 additions & 4 deletions RecoLocalCalo/HGCalRecAlgos/src/ClusterTools.cc
View file
Open in desktop
Expand Up @@ -11,6 +11,8 @@
#include "FWCore/Framework/interface/EventSetup.h"
#include "DataFormats/Common/interface/Handle.h"

#include "vdt/vdtMath.h"

using namespace hgcal;

ClusterTools::ClusterTools(const edm::ParameterSet& conf,
Expand Down Expand Up @@ -83,14 +85,16 @@ math::XYZPoint ClusterTools::getMultiClusterPosition(const reco::HGCalMultiClust
for( const auto& ptr : clu.clusters() ) {
const double x = ptr->x();
const double y = ptr->y();
const float point_r = std::sqrt(x*x + y*y);
const double point_r2 = (x*x + y*y);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is double-precision well justifiable for the math here?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

everything else is already in double, repeated conversion from double to float costs time.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

and this is in a tight loop

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

the reason for x,y,z and energy to be in double is debatable on the storage side.
there is nothing else inside this loop that requires double.

OK for now

const double point_z = ptr->z()-vz;
const double point_h = std::sqrt(point_r2 + point_z*point_z);
const double weight = ptr->energy() * ptr->size();
assert((y != 0. || x != 0.) && "Cluster position somehow in beampipe.");
assert(point_z != 0. && "Layer-cluster position given as reference point.");
assert(point_z != 0.f && "Layer-cluster position given as reference point.");
const double point_r = std::sqrt(point_r2);
acc_rho += point_r * weight;
acc_phi += std::atan2(y,x) * weight;
acc_eta += -1. * std::log(std::tan(0.5*std::atan2(point_r,point_z))) * weight;
acc_phi += vdt::fast_atan2(y,x) * weight;
acc_eta += -1. * vdt::fast_log(point_r/(point_z + point_h)) * weight;
totweight += weight;
}
const double invweight = 1.0/totweight;
Expand Down
12 changes: 7 additions & 5 deletions RecoLocalCalo/HGCalRecAlgos/src/HGCalDepthPreClusterer.cc
View file
Open in desktop
Expand Up @@ -17,10 +17,10 @@ namespace {
return idx;
}

float dist(const edm::Ptr<reco::BasicCluster> &a,
const edm::Ptr<reco::BasicCluster> &b) {
return reco::deltaR(*a,*b);
}
float dist2(const edm::Ptr<reco::BasicCluster> &a,
const edm::Ptr<reco::BasicCluster> &b) {
return reco::deltaR2(*a,*b);
}
}

std::vector<reco::HGCalMultiCluster> HGCalDepthPreClusterer::makePreClusters(const reco::HGCalMultiCluster::ClusterCollection &thecls) const {
Expand All @@ -29,6 +29,8 @@ std::vector<reco::HGCalMultiCluster> HGCalDepthPreClusterer::makePreClusters(con
std::vector<size_t> es = sorted_indices(thecls);
std::vector<int> vused(es.size(),0);
unsigned int used = 0;
const float radius2 = radius*radius;

for(unsigned int i = 0; i < es.size(); ++i) {
if(vused[i]==0) {
reco::HGCalMultiCluster temp;
Expand All @@ -37,7 +39,7 @@ std::vector<reco::HGCalMultiCluster> HGCalDepthPreClusterer::makePreClusters(con
++used;
for(unsigned int j = i+1; j < es.size(); ++j) {
if(vused[j]==0) {
if( dist(thecls[es[i]],thecls[es[j]])<radius && int(thecls[es[i]]->z()*vused[i])>0 ) {
if( dist2(thecls[es[i]],thecls[es[j]]) < radius2 && int(thecls[es[i]]->z()*vused[i])>0 ) {
temp.push_back(thecls[es[j]]);
vused[j]=vused[i];
++used;
Expand Down
27 changes: 10 additions & 17 deletions RecoLocalCalo/HGCalRecAlgos/src/HGCalImagingAlgo.cc
View file
Open in desktop
Expand Up @@ -171,13 +171,6 @@ math::XYZPoint HGCalImagingAlgo::calculatePosition(std::vector<KDNode> &v){
z/total_weight );
}

double HGCalImagingAlgo::distance(const Hexel &pt1, const Hexel &pt2){
const double dx = pt1.x - pt2.x;
const double dy = pt1.y - pt2.y;
return std::sqrt(dx*dx + dy*dy);
}


double HGCalImagingAlgo::calculateLocalDensity(std::vector<KDNode> &nd, KDTree &lp){
double maxdensity = 0.;
for(unsigned int i = 0; i < nd.size(); ++i){
Expand Down Expand Up @@ -209,37 +202,37 @@ double HGCalImagingAlgo::calculateDistanceToHigher(std::vector<KDNode> &nd, KDTr
maxdensity = nd[rs[0]].data.rho;
else
return maxdensity; // there are no hits
double dist = 50.0;
double dist2 = 2500.0;
//start by setting delta for the highest density hit to
//the most distant hit - this is a convention

for(unsigned int j = 0; j < nd.size(); j++){
double tmp = distance(nd[rs[0]].data, nd[j].data);
dist = tmp > dist ? tmp : dist;
double tmp = distance2(nd[rs[0]].data, nd[j].data);
dist2 = tmp > dist2 ? tmp : dist2;
}
nd[rs[0]].data.delta = dist;
nd[rs[0]].data.delta = std::sqrt(dist2);
nd[rs[0]].data.nearestHigher = nearestHigher;

//now we save the largest distance as a starting point

double max_dist = dist;
const double max_dist2 = dist2;

for(unsigned int oi = 1; oi < nd.size(); ++oi){ // start from second-highest density
dist = max_dist;
dist2 = max_dist2;
unsigned int i = rs[oi];
// we only need to check up to oi since hits
// are ordered by decreasing density
// and all points coming BEFORE oi are guaranteed to have higher rho
// and the ones AFTER to have lower rho
for(unsigned int oj = 0; oj < oi; oj++){
unsigned int j = rs[oj];
double tmp = distance(nd[i].data, nd[j].data);
if(tmp <= dist){ //this "<=" instead of "<" addresses the (rare) case when there are only two hits
dist = tmp;
double tmp = distance2(nd[i].data, nd[j].data);
if(tmp <= dist2){ //this "<=" instead of "<" addresses the (rare) case when there are only two hits
dist2 = tmp;
nearestHigher = j;
}
}
nd[i].data.delta = dist;
nd[i].data.delta = std::sqrt(dist2);
nd[i].data.nearestHigher = nearestHigher; //this uses the original unsorted hitlist
}
return maxdensity;
Expand Down