From 4868db2bbc91627680a2d0bd8b6d2b6ef211fdfd Mon Sep 17 00:00:00 2001 From: Erica Brondolin Date: Mon, 30 Jan 2023 17:22:35 +0100 Subject: [PATCH] Avoid dividing by phi that could be zero --- include/CLUEAlgo.h | 12 ++++++------ include/CLUEAlgoGPU.h | 10 +++++----- include/Points.h | 2 +- src/CLUEAlgo.cc | 26 +++++++++++++------------- src/CLUEAlgoGPU.cu | 2 +- src/ClueGaudiAlgorithmWrapper.cpp | 14 ++++++++------ src/ClueGaudiAlgorithmWrapper.h | 2 +- 7 files changed, 35 insertions(+), 33 deletions(-) diff --git a/include/CLUEAlgo.h b/include/CLUEAlgo.h index c8cdedf..d750524 100644 --- a/include/CLUEAlgo.h +++ b/include/CLUEAlgo.h @@ -33,7 +33,7 @@ class CLUEAlgo_T { Points points_; - bool setPoints(int n, float* x, float* y, int* layer, float* weight, float* phi = NULL) { + bool setPoints(int n, float* x, float* y, int* layer, float* weight, float* r = NULL) { points_.clear(); // input variables for(int i=0; i x; std::vector y; - std::vector phi; + std::vector r; std::vector layer; std::vector weight; diff --git a/src/CLUEAlgo.cc b/src/CLUEAlgo.cc index e0c09c0..64ef916 100644 --- a/src/CLUEAlgo.cc +++ b/src/CLUEAlgo.cc @@ -46,7 +46,7 @@ template void CLUEAlgo_T::prepareDataStructures( TILES & allLayerTiles ){ for (int i=0; i::calculateLocalDensity( TILES & allLayerTiles ){ // loop over all points for(unsigned i = 0; i < points_.n; i++) { auto lt = allLayerTiles[points_.layer[i]]; - float ri = points_.x[i]/points_.phi[i]; + float ri = points_.r[i]; + float phi_i = points_.x[i]/(1.*ri); // get search box search_box = lt.searchBox(points_.x[i]-dc_, points_.x[i]+dc_, points_.y[i]-dc_, points_.y[i]+dc_); if(!TILES::constants_type_t::endcap){ float dc_phi = dc_/ri; - search_box = lt.searchBoxPhiZ(points_.phi[i]-dc_phi, points_.phi[i]+dc_phi, points_.y[i]-dc_, points_.y[i]+dc_); + search_box = lt.searchBoxPhiZ(phi_i-dc_phi, phi_i+dc_phi, points_.y[i]-dc_, points_.y[i]+dc_); } // loop over bins in the search box @@ -85,10 +86,8 @@ void CLUEAlgo_T::calculateLocalDensity( TILES & allLayerTiles ){ for (unsigned int binIter = 0; binIter < binSize; binIter++) { int j = lt[binId][binIter]; // query N_{dc_}(i) - float dist_ij = distance(i, j); - if(!TILES::constants_type_t::endcap){ - dist_ij = distance(i, j, true, ri); - } + float dist_ij = TILES::constants_type_t::endcap ? + distance(i, j) : distance(i, j, true, ri); if(dist_ij <= dc_) { // sum weights within N_{dc_}(i) points_.rho[i] += (i == j ? 1.f : 0.5f) * points_.weight[j]; @@ -103,7 +102,6 @@ void CLUEAlgo_T::calculateLocalDensity( TILES & allLayerTiles ){ template void CLUEAlgo_T::calculateDistanceToHigher( TILES & allLayerTiles ){ - // loop over all points float dm = outlierDeltaFactor_ * dc_; for(unsigned i = 0; i < points_.n; i++) { @@ -112,7 +110,8 @@ void CLUEAlgo_T::calculateDistanceToHigher( TILES & allLayerTiles ){ int nearestHigher_i = -1; float xi = points_.x[i]; float yi = points_.y[i]; - float ri = points_.x[i]/points_.phi[i]; + float ri = points_.r[i]; + float phi_i = points_.x[i]/(1.*ri); float rho_i = points_.rho[i]; //get search box @@ -120,7 +119,7 @@ void CLUEAlgo_T::calculateDistanceToHigher( TILES & allLayerTiles ){ float dm_phi = dm/ri; std::array search_box = TILES::constants_type_t::endcap ? lt.searchBox(xi-dm, xi+dm, yi-dm, yi+dm): - lt.searchBoxPhiZ(points_.phi[i]-dm_phi, points_.phi[i]+dm_phi, points_.y[i]-dm, points_.y[i]+dm); + lt.searchBoxPhiZ(phi_i-dm_phi, phi_i+dm_phi, points_.y[i]-dm, points_.y[i]+dm); // loop over all bins in the search box for(int xBin = search_box[0]; xBin < search_box[1]+1; ++xBin) { @@ -143,8 +142,7 @@ void CLUEAlgo_T::calculateDistanceToHigher( TILES & allLayerTiles ){ // in the rare case where rho is the same, use detid foundHigher = foundHigher || ((points_.rho[j] == rho_i) && (j>i) ); float dist_ij = TILES::constants_type_t::endcap ? - dist_ij = distance(i, j): - dist_ij = distance(i, j, true, ri);; + distance(i, j) : distance(i, j, true, ri); if(foundHigher && dist_ij <= dm) { // definition of N'_{dm}(i) // find the nearest point within N'_{dm}(i) if (dist_ij < delta_i) { @@ -234,7 +232,9 @@ inline float CLUEAlgo_T::distance(int i, int j, bool isPhi, float r ) con // 2-d distance on the layer if(points_.layer[i] == points_.layer[j] ) { if (isPhi) { - const float drphi = r * reco::deltaPhi(points_.phi[i], points_.phi[j]); + const float phi_i = points_.x[i]/(1.*points_.r[i]); + const float phi_j = points_.x[j]/(1.*points_.r[j]); + const float drphi = r * reco::deltaPhi(phi_i, phi_j); const float dy = points_.y[i] - points_.y[j]; return std::sqrt(dy * dy + drphi * drphi); } else { diff --git a/src/CLUEAlgoGPU.cu b/src/CLUEAlgoGPU.cu index 52e2b72..25e85a0 100644 --- a/src/CLUEAlgoGPU.cu +++ b/src/CLUEAlgoGPU.cu @@ -21,7 +21,7 @@ __global__ void kernel_compute_histogram( TILES &d_hist, int i = blockIdx.x * blockDim.x + threadIdx.x; if(i < numberOfPoints) { // push index of points into tiles - d_hist.fill(d_points.layer[i], d_points.x[i], d_points.y[i], d_points.phi[i], i); + d_hist.fill(d_points.layer[i], d_points.x[i], d_points.y[i], d_points.x[i]/(1.*d_points.r[i]), i); } } //kernel diff --git a/src/ClueGaudiAlgorithmWrapper.cpp b/src/ClueGaudiAlgorithmWrapper.cpp index ae007f3..eadcb47 100644 --- a/src/ClueGaudiAlgorithmWrapper.cpp +++ b/src/ClueGaudiAlgorithmWrapper.cpp @@ -41,12 +41,12 @@ void ClueGaudiAlgorithmWrapper::fillCLUEPoints(std::vector > ClueGaudiAlgorithmWrapper::runAlgo(std::vector< // Run CLUE info() << "Running CLUEAlgo ... " << endmsg; if(isBarrel){ + info() << "... in the barrel" << endmsg; CLDBarrelCLUEAlgo clueAlgo(dc, rhoc, outlierDeltaFactor, true); - if(clueAlgo.setPoints(x.size(), &x[0], &y[0], &layer[0], &weight[0], &phi[0])) + if(clueAlgo.setPoints(x.size(), &x[0], &y[0], &layer[0], &weight[0], &r[0])) throw error() << "Error in setting the clue points for the barrel." << endmsg; clueAlgo.makeClusters(); clueClusters = clueAlgo.getClusters(); cluePoints = clueAlgo.getPoints(); } else { + info() << "... in the endcap" << endmsg; CLDEndcapCLUEAlgo clueAlgo(dc, rhoc, outlierDeltaFactor, true); - if(clueAlgo.setPoints(x.size(), &x[0], &y[0], &layer[0], &weight[0], &phi[0])) + if(clueAlgo.setPoints(x.size(), &x[0], &y[0], &layer[0], &weight[0], &r[0])) throw error() << "Error in setting the clue points for the endcap." << endmsg; clueAlgo.makeClusters(); clueClusters = clueAlgo.getClusters(); @@ -117,7 +119,7 @@ std::map > ClueGaudiAlgorithmWrapper::runAlgo(std::vector< void ClueGaudiAlgorithmWrapper::cleanCLUEPoints(){ x.clear(); y.clear(); - phi.clear(); + r.clear(); layer.clear(); weight.clear(); } diff --git a/src/ClueGaudiAlgorithmWrapper.h b/src/ClueGaudiAlgorithmWrapper.h index 2d7efc2..48d79b2 100644 --- a/src/ClueGaudiAlgorithmWrapper.h +++ b/src/ClueGaudiAlgorithmWrapper.h @@ -43,7 +43,7 @@ class ClueGaudiAlgorithmWrapper : public GaudiAlgorithm { clue::CLUECalorimeterHitCollection clue_hit_coll; std::vector x; std::vector y; - std::vector phi; + std::vector r; std::vector layer; std::vector weight;