     tempfhd[xyz] = ((double)HostMem.fhd[xyz]/(0xFFFFFFFF*photons_finished));
   }
 
+  // Normalize and write output matrix
+  //for (int xyz = 0; xyz < timegrid_size; xyz++) {
+  //  tgrid[xyz] = ((double)HostMem.temp_xyt[xyz]/(0xFFFFFFFF*photons_finished));
+  //}
+
   printf ("Photons simulated: %llu\n\n", photons_finished);
   FreeMemStructs(&HostMem, &DeviceMem);
   return photons_finished;

__device__ __constant__ IncStruct inclusion_dc[1];
__device__ __constant__ unsigned int ignoreAdetection_dc[1];
__device__ __constant__ unsigned int fhd_activated_dc[1];
__device__ __constant__ unsigned int bulk_method_dc[1];
__device__ __constant__ float xi_dc[1];
__device__ __constant__ float yi_dc[1];


  const int fhd_size = num_x * num_y * num_z; //x + HEIGHT* (y + WIDTH* z)


  // Store in local variables the number of time detectors

  const int num_x_tdet = simulations[0].det.x_temp_numdets;
  const int num_y_tdet = simulations[0].det.y_temp_numdets;
  const long num_tbins = simulations[0].det.temp_bins;

  const long timegrid_size = num_x_tdet * num_y_tdet * num_tbins;

  // FHD simulation
  // Run a simulation

