                                          , d_radii2
                                          , n_radii
                                          , d_in_radius);
-      // compute pops / radius
+      // compute pops per radius
       for (unsigned int r=0; r < n_radii; ++r) {
-        //TODO define range / offset
-        reduce_sum<BSIZE> <<< rng
+        // pops stored col-wise -> just set an offset ...
+        offset = r*n_rows;
+        reduce_sum<BSIZE> <<< n_rows_ext
                             , BSIZE
                             , 0
                             , streams[i_stream] >>> (offset

    float* d_radii2;
    cudaMalloc((void**) &d_radii2
             , sizeof(float) * n_radii);
    for (float& r: radii) {
      r *= r;
    }
    cudaMemcpy(d_radii2
             , radii.data()
             , sizeof(float) * n_radii
             , cudaMemcpyHostToDevice);
    // tmp buffer for in/out info & reference coords (per stream)

                              , int i_gpu) {
    unsigned int n_rows_ext = Tools::min_multiplicator(n_rows, BSIZE) * BSIZE;
    unsigned int n_radii = radii.size();
    // make sure radii are in descending order
    std::sort(radii.begin()
            , radii.end()
            , [](float lhs, float rhs) -> bool {
                return lhs > rhs;
              });
    // setup device & streams
    cudaSetDevice(i_gpu);
    cudaStream_t streams[N_STREAMS];
