 
   // Compute the efficient spectrum at the current pixel
   float efficientSpectrum[nBins * nEnergies];
-  for (unsigned int r=0; r<nBins; r++)
-    for (unsigned int c=0; c<nEnergies; c++)
-      efficientSpectrum[c + nEnergies * r] = pSpectrum[c + nEnergies * first_proj_idx] * c_binnedDetectorResponse[c + nEnergies * r];
+  for (unsigned int b=0; b<nBins; b++)
+    for (unsigned int e=0; e<nEnergies; e++)
+      efficientSpectrum[IDX2D(b,e,nEnergies)] = pSpectrum[e + nEnergies * first_proj_idx] * c_binnedDetectorResponse[IDX2D(b,e,nEnergies)];
 
   // Get attenuation factors at each energy from material projections
   float attenuationFactors[nEnergies];

#include <cublas_v2.h>
#include <cuda_runtime.h>

#define IDX2D(r,c,cols) = ((r)*(cols)+(c))

// CONSTANTS //////////////////////////////////////////////////////////////
__constant__ int3 c_projSize;

                         nEnergies);

  // Sum on the bins
  for (unsigned int b=0; b<nBins; b++)
    for (unsigned int c=0; c<nMaterials * nMaterials; c++)
      pOut2[proj_idx * nMaterials * nMaterials + c] += interm2ForHessian[IDX2D(b,c,nMaterials * nMaterials)];

  // Multiply by the projection of ones
  for (unsigned int c=0; c<nMaterials * nMaterials; c++)
