Avoid dividing by phi that could be zero

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<n; ++i)
@@ -42,14 +42,14 @@ class CLUEAlgo_T {
 	points_.y.push_back(y[i]);
 	points_.layer.push_back(layer[i]);
 	points_.weight.push_back(weight[i]);
-	if(phi != NULL){
-          points_.phi.push_back(phi[i]);
+	if(r != NULL){
+          points_.r.push_back(r[i]);
         } else {
-          // If the layer tile is declared as endcap, the phi info is not used
+          // If the layer tile is declared as endcap, the r info is not used
           if(TILES::constants_type_t::endcap){
-            points_.phi.push_back(0.0);
+            points_.r.push_back(0.0);
           } else {
-            std::cerr << "ERROR: phi info is not present but you are using a barrel LayerTile! " << std::endl;
+            std::cerr << "ERROR: r info is not present but you are using a barrel LayerTile! " << std::endl;
           }
         }
       }

include/CLUEAlgoGPU.h

@@ -10,12 +10,12 @@ static const int localStackSizePerSeed = 20;
 
 struct PointsPtr {
   float *x; 
-  float *y ;
-  float *phi ;
-  int *layer ;
-  float *weight ;
+  float *y;
+  float *r;
+  int *layer;
+  float *weight;
 
-  float *rho ; 
+  float *rho; 
   float *delta; 
   int *nearestHigher;
   int *clusterIndex; 

include/Points.h

@@ -5,7 +5,7 @@ struct Points {
 
   std::vector<float> x;
   std::vector<float> y;
-  std::vector<float> phi;
+  std::vector<float> r;
   std::vector<int> layer;
   std::vector<float> weight;
 

src/CLUEAlgo.cc

@@ -46,7 +46,7 @@ template <typename TILES>
 void CLUEAlgo_T<TILES>::prepareDataStructures( TILES & allLayerTiles ){
   for (int i=0; i<points_.n; i++){
     // push index of points into tiles
-    allLayerTiles.fill( points_.layer[i], points_.x[i], points_.y[i], points_.phi[i], i );
+    allLayerTiles.fill( points_.layer[i], points_.x[i], points_.y[i], points_.x[i]/(1.*points_.r[i]), i );
   }
 }
 
@@ -58,14 +58,15 @@ void CLUEAlgo_T<TILES>::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<TILES>::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<TILES>::calculateLocalDensity( TILES & allLayerTiles ){
 
 template <typename TILES>
 void CLUEAlgo_T<TILES>::calculateDistanceToHigher( TILES & allLayerTiles ){
-
   // loop over all points
   float dm = outlierDeltaFactor_ * dc_;
   for(unsigned i = 0; i < points_.n; i++) {
@@ -112,15 +110,16 @@ void CLUEAlgo_T<TILES>::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
     auto lt = allLayerTiles[points_.layer[i]];
     float dm_phi = dm/ri;
     std::array<int,4> 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<TILES>::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<TILES>::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 {

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
 

src/ClueGaudiAlgorithmWrapper.cpp

@@ -41,12 +41,12 @@ void ClueGaudiAlgorithmWrapper::fillCLUEPoints(std::vector<clue::CLUECalorimeter
     if(ch.inBarrel()){
       x.push_back(ch.getPhi()*ch.getR());
       y.push_back(ch.getPosition().z);
-      phi.push_back(ch.getPhi());
+      r.push_back(ch.getR());
     } else {
       x.push_back(ch.getPosition().x);
       y.push_back(ch.getPosition().y);
-      // For the endcap the phi info is not mandatory because it is not used
-      phi.push_back(ch.getPhi());
+      // For the endcap the r info is not mandatory because it is not used
+      r.push_back(ch.getR());
     }
     layer.push_back(ch.getLayer());
     weight.push_back(ch.getEnergy());
@@ -67,15 +67,17 @@ std::map<int, std::vector<int> > 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<int, std::vector<int> > ClueGaudiAlgorithmWrapper::runAlgo(std::vector<
 void ClueGaudiAlgorithmWrapper::cleanCLUEPoints(){
   x.clear();
   y.clear();
-  phi.clear();
+  r.clear();
   layer.clear();
   weight.clear();
 }

src/ClueGaudiAlgorithmWrapper.h

@@ -43,7 +43,7 @@ class ClueGaudiAlgorithmWrapper : public GaudiAlgorithm {
   clue::CLUECalorimeterHitCollection clue_hit_coll;
   std::vector<float> x;
   std::vector<float> y;
-  std::vector<float> phi;
+  std::vector<float> r;
   std::vector<int> layer;
   std::vector<float> weight;