 
   // 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;


  // Get intermediate variables used in the computation of the first output
  float oneMinusRatios[nBins];
  for (unsigned int b=0; b<nBins; b++)
    oneMinusRatios[b] = 1 - (pPhoCount[proj_idx * nBins + b] / expectedCounts[b]);

  // Form an intermediate variable used for the gradient of the cost function,
  // (the derivation of the exponential implies that a m_MaterialAttenuations
  // gets out), by equivalent of element-wise product with implicit extension
  float intermForGradient[nEnergies * nMaterials];
  for (unsigned int e=0; e<nEnergies; e++)
    for (unsigned int m=0; m<nMaterials; m++)
      intermForGradient[IDX2D(e,m,nMaterials)] = c_materialAttenuations[IDX2D(e,m,nMaterials)] * attenuationFactors[e];

  // Multiply by the spectrum
  float interm2ForGradient[nBins * nMaterials];
