<commit_before>/*	This file is part of CUDAMCFL.

    CUDAMCFL is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    CUDAMCFL is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with CUDAMCFL.  If not, see <http://www.gnu.org/licenses/>.
*/

// forward declaration of the device code
__global__ void MCd(MemStruct);
__global__ void MCd3D(MemStruct);
__device__ float rand_MWC_oc(unsigned long long*,unsigned int*);
__device__ float rand_MWC_co(unsigned long long*,unsigned int*);
__device__ void LaunchPhoton(PhotonStruct*,unsigned long long*, unsigned int*, MemStruct);
__global__ void LaunchPhoton_Global(MemStruct);
__device__ void Spin(PhotonStruct*, float,unsigned long long*,unsigned int*);
__device__ unsigned int Reflect(PhotonStruct*, int, unsigned long long*, unsigned int*, int);
__device__ unsigned int PhotonSurvive(PhotonStruct*, unsigned long long*, unsigned int*);
__device__ unsigned int MoveToFirstBoundary(PhotonStruct*, unsigned short, short*, float);

__global__ void MCd(MemStruct DeviceMem)
{
  // Block index
  int bx=blockIdx.x;

  // Thread index
  int tx=threadIdx.x;

  // First element processed by the block
  int begin=NUM_THREADS_PER_BLOCK*bx;

  // Bulk thickness
	const float esp=layers_dc[(*n_layers_dc)].z_max;

  // 3D bulk matrix width and height
  const unsigned int num_x = __float2uint_rn(4*esp*__int2float_rn(*grid_size_dc));
  const unsigned int num_y = __float2uint_rn(4*esp*__int2float_rn(*grid_size_dc));

  // Size of output images
  const float size_x = __fdividef(det_dc[0].dx*__int2float_rn(det_dc[0].nx),2.);
  const float size_y = __fdividef(det_dc[0].dy*__int2float_rn(det_dc[0].ny),2.);

  // Size of time array
  const int num_x_tdet = det_dc[0].x_temp_numdets;
  const int num_y_tdet = det_dc[0].y_temp_numdets;
  const long num_tbins = det_dc[0].temp_bins;

  unsigned long long int x=DeviceMem.x[begin+tx];//coherent
	unsigned int a=DeviceMem.a[begin+tx];//coherent

	float s;	// step length
	unsigned int index; // temporal variable to store indexes to arrays
	unsigned int w; // photon weight

	PhotonStruct p = DeviceMem.p[begin+tx];

	int new_layer;

	//First, make sure the thread (photon) is active
	unsigned int ii = 0;
	if(!DeviceMem.thread_active[begin+tx]) ii = NUMSTEPS_GPU;

	for(;ii<NUMSTEPS_GPU;ii++) {
    //this is the main while loop

    // Check if inside spherical inclusion
    if((((p.x-inclusion_dc[0].x)*(p.x-inclusion_dc[0].x)) +
       ((p.y-inclusion_dc[0].y)*(p.y-inclusion_dc[0].y)) +
       ((p.z-inclusion_dc[0].z)*(p.z-inclusion_dc[0].z))) <=
       (inclusion_dc[0].r*inclusion_dc[0].r)) {
      // Inside inclusion
			if(inclusion_dc[0].mutr!=FLT_MAX)
				//s = -__logf(rand_MWC_oc(&x,&a))*inclusion_dc[0].mutr;//sample step length [cm] //HERE AN OPEN_OPEN FUNCTION WOULD BE APPRECIATED
        s = inclusion_dc[0].mutr;
      else
				s = 2.0f;//temporary, say the step in glass is 100 cm.
		}
		else 	{
      // Outside inclusion
			if(layers_dc[p.layer].mutr!=FLT_MAX)
				//s = -__logf(rand_MWC_oc(&x,&a))*layers_dc[p.layer].mutr;//sample step length [cm] //HERE AN OPEN_OPEN FUNCTION WOULD BE APPRECIATED
        s = layers_dc[p.layer].mutr;
      else
				s = 2.0f;//temporary, say the step in glass is 100 cm.
		}

		//Check for layer transitions and in case, calculate s
		new_layer = p.layer;

    //Check for upwards reflection/transmission & calculate new s
		if(p.z+s*p.dz<layers_dc[p.layer].z_min) {
      new_layer--;
      s = __fdividef(layers_dc[p.layer].z_min-p.z,p.dz);
      }

    //Check for downward reflection/transmission
		if(p.z+s*p.dz>=layers_dc[p.layer].z_max){
      new_layer++;
      s = __fdividef(layers_dc[p.layer].z_max-p.z,p.dz);
    }

    //Move photon
    p.x += p.dx*s;
		p.y += p.dy*s;
		p.z += p.dz*s;

    //Update time of flight
    p.tof += (unsigned long)(s/C_CMFS);

    if(p.z>layers_dc[p.layer].z_max)p.z=layers_dc[p.layer].z_max;//needed? TODO
		if(p.z<layers_dc[p.layer].z_min)p.z=layers_dc[p.layer].z_min;//needed? TODO

    // Accumulate photon hitting density if needed
    if (*fhd_activated_dc == 1){
      if(fabsf(p.x)<2*esp && fabsf(p.y)<2*esp && p.z<esp){ //Inside space of fhd
        //Use round to zero so there are no over sampled voxels (for ex: (max_x,0,0) and (0,1,0) should not map to the same voxel)
        index = __float2uint_rz((p.x+2*esp)*__int2float_rn(*grid_size_dc))
                  + num_x * (__float2uint_rz((p.y+2*esp)*__int2float_rn(*grid_size_dc))
                  + num_y * __float2uint_rz((p.z)*__int2float_rn(*grid_size_dc))); //x + HEIGHT* (y + WIDTH* z) Fx[ix + num_x * (iy + iz * num_y)]

        if (DeviceMem.fhd[index] + p.weight < LLONG_MAX) atomicAdd(&DeviceMem.fhd[index], p.weight); // Check for overflow and add atomically //TODO why LLONG_MAX?
      }
    }

    // New layer?
		if(new_layer!=p.layer) {
      // set the remaining step length to 0
			s = 0.0f;

      if(Reflect(&p,new_layer,&x,&a,1)==0u) {
        // Photon is transmitted

        if(new_layer == 0){
          // Diffuse reflectance
          if(fabsf(p.x-det_dc[0].x0)<size_x && fabsf(p.y-det_dc[0].y0)<size_y) {
            // Photon is detectable
            // Use round to zero so there are no over sampled pixels (for ex: (max_x,0) and (0,1) should not map to the same pixel)
            index=__float2uint_rz(__fdividef(p.y-det_dc[0].y0+size_y,det_dc[0].dy)) * det_dc[0].nx +
                  __float2uint_rz(__fdividef(p.x-det_dc[0].x0+size_x,det_dc[0].dx));
    				if ((DeviceMem.Rd_xy[index] + p.weight) < LLONG_MAX) atomicAdd(&DeviceMem.Rd_xy[index], p.weight); // Check for overflow and add atomicall
    			}
          if (*do_temp_sim_dc==1u && det_dc[0].temp_rort==0u && p.tof<det_dc[0].max_temp){
            // Save time value in apropiate bin
            for (int xpos = 0; xpos < num_x_tdet; xpos++){
              for (int ypos = 0; ypos < num_y_tdet; ypos++){
                if ((fabsf(p.x - DeviceMem.tdet_pos_x[xpos]) + fabsf(p.y - DeviceMem.tdet_pos_y[ypos])) < (det_dc[0].temp_det_r) * (det_dc[0].temp_det_r)){
                  // Inside time detector
                  index = xpos + num_x_tdet*(ypos + num_y_tdet*__float2uint_rz(__fdividef(p.tof, det_dc[0].max_temp)*num_tbins));
                  if (DeviceMem.time_xyt[index] + p.weight < LLONG_MAX) atomicAdd(&DeviceMem.time_xyt[index], p.weight); // Check for overflow and add atomically //TODO why LLONG_MAX?
                }
              }
            }
          }
          p.weight = 0; // Set the remaining weight to 0, effectively killing the photon
        }

        if(new_layer > *n_layers_dc) {
          // Diffuse transmitance
          if(fabsf(p.x-det_dc[0].x0)<size_x && fabsf(p.y-det_dc[0].y0)<size_y) {
            // Photon transmitted
            // Use round to zero so there are no over sampled pixels (for ex: (max_x,0) and (0,1) should not map to the same pixel)
            index=__float2uint_rz(__fdividef(p.y-det_dc[0].y0+size_y,det_dc[0].dy)) * det_dc[0].nx +
                  __float2uint_rz(__fdividef(p.x-det_dc[0].x0+size_x,det_dc[0].dx));
            if ((DeviceMem.Tt_xy[index] + p.weight) < LLONG_MAX) atomicAdd(&DeviceMem.Tt_xy[index], p.weight); // Check for overflow and add atomically
          }
          if (*do_temp_sim_dc==1u && det_dc[0].temp_rort==1u && p.tof<det_dc[0].max_temp){
            // Save time value in apropiate bin
            for (int xpos = 0; xpos < num_x_tdet; xpos++){
              for (int ypos = 0; ypos < num_y_tdet; ypos++){
                if ((fabsf(p.x - DeviceMem.tdet_pos_x[xpos]) + fabsf(p.y - DeviceMem.tdet_pos_y[ypos])) < (det_dc[0].temp_det_r) * (det_dc[0].temp_det_r)){
                  // Inside time detector
                  index = xpos + num_x_tdet*(ypos + num_y_tdet*__float2uint_rz(__fdividef(p.tof, det_dc[0].max_temp)*num_tbins));
                  if (DeviceMem.time_xyt[index] + p.weight < LLONG_MAX) atomicAdd(&DeviceMem.time_xyt[index], p.weight); // Check for overflow and add atomically //TODO why LLONG_MAX?
                }
              }
            }
          }


					p.weight = 0; // Set the remaining weight to 0, killing the photon
        }
			}
		}

		w=0;

		if(s > 0.0f) {
			// Drop weight (apparently only when the photon is scattered)
			if((((p.x-inclusion_dc[0].x)*(p.x-inclusion_dc[0].x)) +
          ((p.y-inclusion_dc[0].y)*(p.y-inclusion_dc[0].y)) +
          ((p.z-inclusion_dc[0].z)*(p.z-inclusion_dc[0].z))) <=
          (inclusion_dc[0].r*inclusion_dc[0].r)) {
            // Inside inclusion
				    w = __float2uint_rn(inclusion_dc[0].mua*inclusion_dc[0].mutr*__uint2float_rn(p.weight));
            if (p.weight - w >= 0) // Check for underflow
              p.weight -= w;
            else
              p.weight = 0;
				    Spin(&p,inclusion_dc[0].g,&x,&a);
      }
			else {
        // Outside inclusion
				w = __float2uint_rn(layers_dc[p.layer].mua*layers_dc[p.layer].mutr*__uint2float_rn(p.weight));
        if (p.weight - w >= 0) // Check for underflow
          p.weight -= w;
        else
          p.weight = 0;
				Spin(&p,layers_dc[p.layer].g,&x,&a);
			}
		}

		if(!PhotonSurvive(&p,&x,&a)) // Check if photons survives or not
		{
			if(atomicAdd(DeviceMem.num_terminated_photons,1ULL) < (*num_photons_dc-NUM_THREADS)) {
        // Ok to launch another photon
				LaunchPhoton(&p,&x,&a, DeviceMem);//Launch a new photon
      }
			else {
        // No more photons should be launched.
				DeviceMem.thread_active[begin+tx] = 0u; // Set thread to inactive
				ii = NUMSTEPS_GPU;				// Exit main loop
			}

		}
	}//end main for loop!
	__syncthreads();//necessary?

	//save the state of the MC simulation in global memory before exiting
	DeviceMem.p[begin+tx] = p;	//This one is incoherent!!!
	DeviceMem.x[begin+tx] = x; //this one also seems to be coherent

}//end MCd


__global__ void MCd3D(MemStruct DeviceMem)
{

  // Block index
  int bx=blockIdx.x;

  // Thread index
  int tx=threadIdx.x;

  // First element processed by the block
  int begin=NUM_THREADS_PER_BLOCK*bx;

  // 3D bulk matrix width and height
  const unsigned int num_x = __float2uint_rn(4*(*esp_dc)*__int2float_rn(*grid_size_dc));
  const unsigned int num_y = __float2uint_rn(4*(*esp_dc)*__int2float_rn(*grid_size_dc));

  // Size of output images
  const float size_x = __fdividef(det_dc[0].dx*__int2float_rn(det_dc[0].nx),2.);
  const float size_y = __fdividef(det_dc[0].dy*__int2float_rn(det_dc[0].ny),2.);

  // Size of time array
  const int num_x_tdet = det_dc[0].x_temp_numdets;
  const int num_y_tdet = det_dc[0].y_temp_numdets;
  const long num_tbins = det_dc[0].temp_bins;

  // Last Bulk
  const unsigned short last_bulk = *n_bulks_dc+1;

  unsigned long long int x=DeviceMem.x[begin+tx];//coherent
	unsigned int a=DeviceMem.a[begin+tx];//coherent

	float s;	// step length
	unsigned int index; // temporal variable to store indexes to arrays
	unsigned int w; // photon weight
  //bool in_glass = FALSE;

	PhotonStruct p = DeviceMem.p[begin+tx];

	unsigned int new_bulk; // int storing bulk descriptor of current position

	// First, make sure the thread (photon) is active
	unsigned int ii = 0;
	if(!DeviceMem.thread_active[begin+tx]) ii = NUMSTEPS_GPU;

	for(;ii<NUMSTEPS_GPU;ii++) {
    // Main while loop

  	if(bulks_dc[p.bulkpos].mutr!=FLT_MAX) {
			s = -__logf(rand_MWC_oc(&x,&a))*bulks_dc[p.bulkpos].mutr;//sample step length [cm] //HERE AN OPEN_OPEN FUNCTION WOULD BE APPRECIATED
      //s = bulks_dc[p.bulkpos].mutr;

      new_bulk = p.bulkpos;

      //int side_scape = 0;

      //Check for upwards reflection/transmission and move to surface
  		if(p.z+s*p.dz<0.) {
        new_bulk = 0;
        //side_scape = 0;
        s = __fdividef(-p.z,p.dz);
      }

      //Check for downward reflection/transmission and move to surface
  		if(p.z+s*p.dz>(*esp_dc)){
        new_bulk = last_bulk;
        //side_scape = 0;
        s = __fdividef((*esp_dc)-p.z,p.dz);
      }

      // Move photon TODO:here?
      p.x += p.dx*s;
      p.y += p.dy*s;
      p.z += p.dz*s;

      //Update time of flight
      p.tof += (unsigned long)(s/C_CMFS);

      // Retrieve bulk position
      if(new_bulk!=0 && new_bulk!=last_bulk) {
        if (fabsf(p.x)<2*(*esp_dc) && fabsf(p.y)<2*(*esp_dc) && p.z<(*esp_dc)){
          // Inside space of 3D matrix. Calculate index of the current voxel
          // Use round to zero so there are no over sampled voxels (for ex: (max_x,0,0) and (0,1,0) should not map to the same voxel)
          index = __float2uint_rz((p.x+2*(*esp_dc))*__int2float_rn(*grid_size_dc))
                + num_x * (__float2uint_rz((p.y+2*(*esp_dc))*__int2float_rn(*grid_size_dc))
                + num_y * __float2uint_rz((p.z)*__int2float_rn(*grid_size_dc)));
          new_bulk = DeviceMem.bulk_info[index];
        }
        else {
          // Photon scaped to the sides
          // Outside space of 3D matrix, assume inside homogeneous medium (should we assume it is outside the bulk (bulkpos=0)? TODO)
          new_bulk = 1;
        }
      }

    }
    else {
			//s = 100.0f; //temporary, say the step in glass is 100 cm.
      //in_glass=TRUE;
      s = 100.0f;
      new_bulk = MoveToFirstBoundary(&p, p.bulkpos, DeviceMem.bulk_info, s);
      //printf("chau\n");

    }

		//Check for layer transitions and in case, calculate s


    //if(p.z>(*esp_dc)) p.z=(*esp_dc);//needed? TODO
		//if(p.z<0.) p.z=0.;//needed? TODO

    unsigned int reflected;

		if(new_bulk != p.bulkpos) {
      // If changing voxel do Reflect
      reflected = Reflect(&p,new_bulk,&x,&a,2);
    }

    // Accumulate photon hitting density if needed
    if(fabsf(p.x)<2*(*esp_dc) && fabsf(p.y)<2*(*esp_dc) && p.z<(*esp_dc) && *fhd_activated_dc == 1) {
      // Inside space of 3D matrix. Calculate index of the current voxel
      // Use round to zero so there are no over sampled voxels (for ex: (max_x,0,0) and (0,1,0) should not map to the same voxel)
      index = __float2uint_rz((p.x+2*(*esp_dc))*__int2float_rn(*grid_size_dc))
            + num_x * (__float2uint_rz((p.y+2*(*esp_dc))*__int2float_rn(*grid_size_dc))
            + num_y * __float2uint_rz((p.z)*__int2float_rn(*grid_size_dc)));
      if (DeviceMem.fhd[index] + p.weight < LLONG_MAX) atomicAdd(&DeviceMem.fhd[index], p.weight); // Check for overflow and add atomically
      //TODO why LLONG_MAX?
    }


    if(p.bulkpos == 0){
      // Photon is outside bulk and reflected
      if(fabsf(p.x-det_dc[0].x0)<size_x && fabsf(p.y-det_dc[0].y0)<size_y) {
        // Photon is detectable
        // Use round to zero so there are no over sampled pixels (for ex: (max_x,0) and (0,1) should not map to the same pixel)
        index=__float2uint_rz(__fdividef(p.y-det_dc[0].y0+size_y,det_dc[0].dy)) * det_dc[0].nx +
              __float2uint_rz(__fdividef(p.x-det_dc[0].x0+size_x,det_dc[0].dx));
				if ((DeviceMem.Rd_xy[index] + p.weight) < LLONG_MAX) atomicAdd(&DeviceMem.Rd_xy[index], p.weight); // Check for overflow and add atomicall
				}

      if (*do_temp_sim_dc==1u && det_dc[0].temp_rort==0u && p.tof<det_dc[0].max_temp){
        // Save time value in apropiate bin
        for (int xpos = 0; xpos < num_x_tdet; xpos++){
          for (int ypos = 0; ypos < num_y_tdet; ypos++){
            if ((fabsf(p.x - DeviceMem.tdet_pos_x[xpos]) + fabsf(p.y - DeviceMem.tdet_pos_y[ypos])) < (det_dc[0].temp_det_r) * (det_dc[0].temp_det_r)){
              // Inside time detector
              index = xpos + num_x_tdet*(ypos + num_y_tdet*__float2uint_rz(__fdividef(p.tof, det_dc[0].max_temp)*num_tbins));
              if (DeviceMem.time_xyt[index] + p.weight < LLONG_MAX) atomicAdd(&DeviceMem.time_xyt[index], p.weight); // Check for overflow and add atomically //TODO why LLONG_MAX?
            }
          }
        }
      }

      p.weight = 0u; // Set the remaining weight to 0, effectively killing the photon
      s = 0.0f;
    }

    if(p.bulkpos == last_bulk){
      // Photon is outside bulk and transmitted
      if(fabsf(p.x-det_dc[0].x0)<size_x && fabsf(p.y-det_dc[0].y0)<size_y) {
        // Photon transmitted
        // Use round to zero so there are no over sampled pixels (for ex: (max_x,0) and (0,1) should not map to the same pixel)
        index=__float2uint_rz(__fdividef(p.y-det_dc[0].y0+size_y,det_dc[0].dy)) * det_dc[0].nx +
              __float2uint_rz(__fdividef(p.x-det_dc[0].x0+size_x,det_dc[0].dx));
        if ((DeviceMem.Tt_xy[index] + p.weight) < LLONG_MAX) atomicAdd(&DeviceMem.Tt_xy[index], p.weight); // Check for overflow and add atomically
        }

      if (*do_temp_sim_dc==1u && det_dc[0].temp_rort==1u && p.tof<det_dc[0].max_temp){
        // Save time value in apropiate bin
        for (int xpos = 0; xpos < num_x_tdet; xpos++){
          for (int ypos = 0; ypos < num_y_tdet; ypos++){
            if ((fabsf(p.x - DeviceMem.tdet_pos_x[xpos]) + fabsf(p.y - DeviceMem.tdet_pos_y[ypos])) < (det_dc[0].temp_det_r) * (det_dc[0].temp_det_r)){
              // Inside time detector
              index = xpos + num_x_tdet*(ypos + num_y_tdet*__float2uint_rz(__fdividef(p.tof, det_dc[0].max_temp)*num_tbins));
              if (DeviceMem.time_xyt[index] + p.weight < LLONG_MAX) atomicAdd(&DeviceMem.time_xyt[index], p.weight); // Check for overflow and add atomically //TODO why LLONG_MAX?
            }
          }
        }
      }

      p.weight = 0u; // Set the remaining weight to 0, effectively killing the photon
      s = 0.0f;
      }

		w=0;

		if(s > 0.0f) {
    	// Drop weight (apparently only when the photon is scattered)
			w = __float2uint_ru(bulks_dc[p.bulkpos].mua*bulks_dc[p.bulkpos].mutr*__uint2float_rn(p.weight));
			if (p.weight - w >= 0 && w > 0) // Check for underflow
        p.weight -= w;
      else
        p.weight = 0u;
			Spin(&p,bulks_dc[p.bulkpos].g,&x,&a);
		}

		if(!PhotonSurvive(&p,&x,&a)) {
      // Check if photons survives or not
			if(atomicAdd(DeviceMem.num_terminated_photons,1ULL) < (*num_photons_dc-NUM_THREADS)) {
        // Ok to launch another photon
				LaunchPhoton(&p,&x,&a, DeviceMem);//Launch a new photon
      }
			else {<commit_msg>Fix bug in the calculation of the distance between the photon and the detector<commit_after>/*	This file is part of CUDAMCFL.

    CUDAMCFL is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    CUDAMCFL is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with CUDAMCFL.  If not, see <http://www.gnu.org/licenses/>.
*/

// forward declaration of the device code
__global__ void MCd(MemStruct);
__global__ void MCd3D(MemStruct);
__device__ float rand_MWC_oc(unsigned long long*,unsigned int*);
__device__ float rand_MWC_co(unsigned long long*,unsigned int*);
__device__ void LaunchPhoton(PhotonStruct*,unsigned long long*, unsigned int*, MemStruct);
__global__ void LaunchPhoton_Global(MemStruct);
__device__ void Spin(PhotonStruct*, float,unsigned long long*,unsigned int*);
__device__ unsigned int Reflect(PhotonStruct*, int, unsigned long long*, unsigned int*, int);
__device__ unsigned int PhotonSurvive(PhotonStruct*, unsigned long long*, unsigned int*);
__device__ unsigned int MoveToFirstBoundary(PhotonStruct*, unsigned short, short*, float);

__global__ void MCd(MemStruct DeviceMem)
{
  // Block index
  int bx=blockIdx.x;

  // Thread index
  int tx=threadIdx.x;

  // First element processed by the block
  int begin=NUM_THREADS_PER_BLOCK*bx;

  // Bulk thickness
	const float esp=layers_dc[(*n_layers_dc)].z_max;

  // 3D bulk matrix width and height
  const unsigned int num_x = __float2uint_rn(4*esp*__int2float_rn(*grid_size_dc));
  const unsigned int num_y = __float2uint_rn(4*esp*__int2float_rn(*grid_size_dc));

  // Size of output images
  const float size_x = __fdividef(det_dc[0].dx*__int2float_rn(det_dc[0].nx),2.);
  const float size_y = __fdividef(det_dc[0].dy*__int2float_rn(det_dc[0].ny),2.);

  // Size of time array
  const int num_x_tdet = det_dc[0].x_temp_numdets;
  const int num_y_tdet = det_dc[0].y_temp_numdets;
  const long num_tbins = det_dc[0].temp_bins;

  unsigned long long int x=DeviceMem.x[begin+tx];//coherent
	unsigned int a=DeviceMem.a[begin+tx];//coherent

	float s;	// step length
	unsigned int index; // temporal variable to store indexes to arrays
	unsigned int w; // photon weight

	PhotonStruct p = DeviceMem.p[begin+tx];

	int new_layer;

	//First, make sure the thread (photon) is active
	unsigned int ii = 0;
	if(!DeviceMem.thread_active[begin+tx]) ii = NUMSTEPS_GPU;

	for(;ii<NUMSTEPS_GPU;ii++) {
    //this is the main while loop

    // Check if inside spherical inclusion
    if((((p.x-inclusion_dc[0].x)*(p.x-inclusion_dc[0].x)) +
       ((p.y-inclusion_dc[0].y)*(p.y-inclusion_dc[0].y)) +
       ((p.z-inclusion_dc[0].z)*(p.z-inclusion_dc[0].z))) <=
       (inclusion_dc[0].r*inclusion_dc[0].r)) {
      // Inside inclusion
			if(inclusion_dc[0].mutr!=FLT_MAX)
				//s = -__logf(rand_MWC_oc(&x,&a))*inclusion_dc[0].mutr;//sample step length [cm] //HERE AN OPEN_OPEN FUNCTION WOULD BE APPRECIATED
        s = inclusion_dc[0].mutr;
      else
				s = 2.0f;//temporary, say the step in glass is 100 cm.
		}
		else 	{
      // Outside inclusion
			if(layers_dc[p.layer].mutr!=FLT_MAX)
				//s = -__logf(rand_MWC_oc(&x,&a))*layers_dc[p.layer].mutr;//sample step length [cm] //HERE AN OPEN_OPEN FUNCTION WOULD BE APPRECIATED
        s = layers_dc[p.layer].mutr;
      else
				s = 2.0f;//temporary, say the step in glass is 100 cm.
		}

		//Check for layer transitions and in case, calculate s
		new_layer = p.layer;

    //Check for upwards reflection/transmission & calculate new s
		if(p.z+s*p.dz<layers_dc[p.layer].z_min) {
      new_layer--;
      s = __fdividef(layers_dc[p.layer].z_min-p.z,p.dz);
      }

    //Check for downward reflection/transmission
		if(p.z+s*p.dz>=layers_dc[p.layer].z_max){
      new_layer++;
      s = __fdividef(layers_dc[p.layer].z_max-p.z,p.dz);
    }

    //Move photon
    p.x += p.dx*s;
		p.y += p.dy*s;
		p.z += p.dz*s;

    //Update time of flight
    p.tof += (unsigned long)(s/C_CMFS);

    if(p.z>layers_dc[p.layer].z_max)p.z=layers_dc[p.layer].z_max;//needed? TODO
		if(p.z<layers_dc[p.layer].z_min)p.z=layers_dc[p.layer].z_min;//needed? TODO

    // Accumulate photon hitting density if needed
    if (*fhd_activated_dc == 1){
      if(fabsf(p.x)<2*esp && fabsf(p.y)<2*esp && p.z<esp){ //Inside space of fhd
        //Use round to zero so there are no over sampled voxels (for ex: (max_x,0,0) and (0,1,0) should not map to the same voxel)
        index = __float2uint_rz((p.x+2*esp)*__int2float_rn(*grid_size_dc))
                  + num_x * (__float2uint_rz((p.y+2*esp)*__int2float_rn(*grid_size_dc))
                  + num_y * __float2uint_rz((p.z)*__int2float_rn(*grid_size_dc))); //x + HEIGHT* (y + WIDTH* z) Fx[ix + num_x * (iy + iz * num_y)]

        if (DeviceMem.fhd[index] + p.weight < LLONG_MAX) atomicAdd(&DeviceMem.fhd[index], p.weight); // Check for overflow and add atomically //TODO why LLONG_MAX?
      }
    }

    // New layer?
		if(new_layer!=p.layer) {
      // set the remaining step length to 0
			s = 0.0f;

      if(Reflect(&p,new_layer,&x,&a,1)==0u) {
        // Photon is transmitted

        if(new_layer == 0){
          // Diffuse reflectance
          if(fabsf(p.x-det_dc[0].x0)<size_x && fabsf(p.y-det_dc[0].y0)<size_y) {
            // Photon is detectable
            // Use round to zero so there are no over sampled pixels (for ex: (max_x,0) and (0,1) should not map to the same pixel)
            index=__float2uint_rz(__fdividef(p.y-det_dc[0].y0+size_y,det_dc[0].dy)) * det_dc[0].nx +
                  __float2uint_rz(__fdividef(p.x-det_dc[0].x0+size_x,det_dc[0].dx));
    				if ((DeviceMem.Rd_xy[index] + p.weight) < LLONG_MAX) atomicAdd(&DeviceMem.Rd_xy[index], p.weight); // Check for overflow and add atomicall
    			}
          if (*do_temp_sim_dc==1u && det_dc[0].temp_rort==0u && p.tof<det_dc[0].max_temp){
            // Save time value in apropiate bin
            for (int xpos = 0; xpos < num_x_tdet; xpos++){
              for (int ypos = 0; ypos < num_y_tdet; ypos++){
                if ((powf((p.x - DeviceMem.tdet_pos_x[xpos]),2) + powf((p.y - DeviceMem.tdet_pos_y[ypos]),2)) < powf((det_dc[0].temp_det_r),2)){
                  // Inside time detector ix + xtnum * (iy + it * ytnum);
                  index = xpos + num_x_tdet*(ypos + num_y_tdet*__float2uint_rz(__fdividef(p.tof, det_dc[0].max_temp)*num_tbins));
                  if (DeviceMem.time_xyt[index] + p.weight < LLONG_MAX) atomicAdd(&DeviceMem.time_xyt[index], p.weight); // Check for overflow and add atomically //TODO why LLONG_MAX?
                }
              }
            }
          }
          p.weight = 0; // Set the remaining weight to 0, effectively killing the photon
        }

        if(new_layer > *n_layers_dc) {
          // Diffuse transmitance
          if(fabsf(p.x-det_dc[0].x0)<size_x && fabsf(p.y-det_dc[0].y0)<size_y) {
            // Photon transmitted
            // Use round to zero so there are no over sampled pixels (for ex: (max_x,0) and (0,1) should not map to the same pixel)
            index=__float2uint_rz(__fdividef(p.y-det_dc[0].y0+size_y,det_dc[0].dy)) * det_dc[0].nx +
                  __float2uint_rz(__fdividef(p.x-det_dc[0].x0+size_x,det_dc[0].dx));
            if ((DeviceMem.Tt_xy[index] + p.weight) < LLONG_MAX) atomicAdd(&DeviceMem.Tt_xy[index], p.weight); // Check for overflow and add atomically
          }
          if (*do_temp_sim_dc==1u && det_dc[0].temp_rort==1u && p.tof<det_dc[0].max_temp){
            // Save time value in apropiate bin
            for (int xpos = 0; xpos < num_x_tdet; xpos++){
              for (int ypos = 0; ypos < num_y_tdet; ypos++){
                if ((powf((p.x - DeviceMem.tdet_pos_x[xpos]),2) + powf((p.y - DeviceMem.tdet_pos_y[ypos]),2)) < powf((det_dc[0].temp_det_r),2)){
                  // Inside time detector
                  index = xpos + num_x_tdet*(ypos + num_y_tdet*__float2uint_rz(__fdividef(p.tof, det_dc[0].max_temp)*num_tbins));
                  if (DeviceMem.time_xyt[index] + p.weight < LLONG_MAX) atomicAdd(&DeviceMem.time_xyt[index], p.weight); // Check for overflow and add atomically //TODO why LLONG_MAX?
                }
              }
            }
          }


					p.weight = 0; // Set the remaining weight to 0, killing the photon
        }
			}
		}

		w=0;

		if(s > 0.0f) {
			// Drop weight (apparently only when the photon is scattered)
			if((((p.x-inclusion_dc[0].x)*(p.x-inclusion_dc[0].x)) +
          ((p.y-inclusion_dc[0].y)*(p.y-inclusion_dc[0].y)) +
          ((p.z-inclusion_dc[0].z)*(p.z-inclusion_dc[0].z))) <=
          (inclusion_dc[0].r*inclusion_dc[0].r)) {
            // Inside inclusion
				    w = __float2uint_rn(inclusion_dc[0].mua*inclusion_dc[0].mutr*__uint2float_rn(p.weight));
            if (p.weight - w >= 0) // Check for underflow
              p.weight -= w;
            else
              p.weight = 0;
				    Spin(&p,inclusion_dc[0].g,&x,&a);
      }
			else {
        // Outside inclusion
				w = __float2uint_rn(layers_dc[p.layer].mua*layers_dc[p.layer].mutr*__uint2float_rn(p.weight));
        if (p.weight - w >= 0) // Check for underflow
          p.weight -= w;
        else
          p.weight = 0;
				Spin(&p,layers_dc[p.layer].g,&x,&a);
			}
		}

		if(!PhotonSurvive(&p,&x,&a)) // Check if photons survives or not
		{
			if(atomicAdd(DeviceMem.num_terminated_photons,1ULL) < (*num_photons_dc-NUM_THREADS)) {
        // Ok to launch another photon
				LaunchPhoton(&p,&x,&a, DeviceMem);//Launch a new photon
      }
			else {
        // No more photons should be launched.
				DeviceMem.thread_active[begin+tx] = 0u; // Set thread to inactive
				ii = NUMSTEPS_GPU;				// Exit main loop
			}

		}
	}//end main for loop!
	__syncthreads();//necessary?

	//save the state of the MC simulation in global memory before exiting
	DeviceMem.p[begin+tx] = p;	//This one is incoherent!!!
	DeviceMem.x[begin+tx] = x; //this one also seems to be coherent

}//end MCd


__global__ void MCd3D(MemStruct DeviceMem)
{

  // Block index
  int bx=blockIdx.x;

  // Thread index
  int tx=threadIdx.x;

  // First element processed by the block
  int begin=NUM_THREADS_PER_BLOCK*bx;

  // 3D bulk matrix width and height
  const unsigned int num_x = __float2uint_rn(4*(*esp_dc)*__int2float_rn(*grid_size_dc));
  const unsigned int num_y = __float2uint_rn(4*(*esp_dc)*__int2float_rn(*grid_size_dc));

  // Size of output images
  const float size_x = __fdividef(det_dc[0].dx*__int2float_rn(det_dc[0].nx),2.);
  const float size_y = __fdividef(det_dc[0].dy*__int2float_rn(det_dc[0].ny),2.);

  // Size of time array
  const int num_x_tdet = det_dc[0].x_temp_numdets;
  const int num_y_tdet = det_dc[0].y_temp_numdets;
  const long num_tbins = det_dc[0].temp_bins;

  // Last Bulk
  const unsigned short last_bulk = *n_bulks_dc+1;

  unsigned long long int x=DeviceMem.x[begin+tx];//coherent
	unsigned int a=DeviceMem.a[begin+tx];//coherent

	float s;	// step length
	unsigned int index; // temporal variable to store indexes to arrays
	unsigned int w; // photon weight
  //bool in_glass = FALSE;

	PhotonStruct p = DeviceMem.p[begin+tx];

	unsigned int new_bulk; // int storing bulk descriptor of current position

	// First, make sure the thread (photon) is active
	unsigned int ii = 0;
	if(!DeviceMem.thread_active[begin+tx]) ii = NUMSTEPS_GPU;

	for(;ii<NUMSTEPS_GPU;ii++) {
    // Main while loop

  	if(bulks_dc[p.bulkpos].mutr!=FLT_MAX) {
			s = -__logf(rand_MWC_oc(&x,&a))*bulks_dc[p.bulkpos].mutr;//sample step length [cm] //HERE AN OPEN_OPEN FUNCTION WOULD BE APPRECIATED
      //s = bulks_dc[p.bulkpos].mutr;

      new_bulk = p.bulkpos;

      //int side_scape = 0;

      //Check for upwards reflection/transmission and move to surface
  		if(p.z+s*p.dz<0.) {
        new_bulk = 0;
        //side_scape = 0;
        s = __fdividef(-p.z,p.dz);
      }

      //Check for downward reflection/transmission and move to surface
  		if(p.z+s*p.dz>(*esp_dc)){
        new_bulk = last_bulk;
        //side_scape = 0;
        s = __fdividef((*esp_dc)-p.z,p.dz);
      }

      // Move photon TODO:here?
      p.x += p.dx*s;
      p.y += p.dy*s;
      p.z += p.dz*s;

      //Update time of flight
      p.tof += (unsigned long)(s/C_CMFS);

      // Retrieve bulk position
      if(new_bulk!=0 && new_bulk!=last_bulk) {
        if (fabsf(p.x)<2*(*esp_dc) && fabsf(p.y)<2*(*esp_dc) && p.z<(*esp_dc)){
          // Inside space of 3D matrix. Calculate index of the current voxel
          // Use round to zero so there are no over sampled voxels (for ex: (max_x,0,0) and (0,1,0) should not map to the same voxel)
          index = __float2uint_rz((p.x+2*(*esp_dc))*__int2float_rn(*grid_size_dc))
                + num_x * (__float2uint_rz((p.y+2*(*esp_dc))*__int2float_rn(*grid_size_dc))
                + num_y * __float2uint_rz((p.z)*__int2float_rn(*grid_size_dc)));
          new_bulk = DeviceMem.bulk_info[index];
        }
        else {
          // Photon scaped to the sides
          // Outside space of 3D matrix, assume inside homogeneous medium (should we assume it is outside the bulk (bulkpos=0)? TODO)
          new_bulk = 1;
        }
      }

    }
    else {
			//s = 100.0f; //temporary, say the step in glass is 100 cm.
      //in_glass=TRUE;
      s = 100.0f;
      new_bulk = MoveToFirstBoundary(&p, p.bulkpos, DeviceMem.bulk_info, s);
      //printf("chau\n");

    }

		//Check for layer transitions and in case, calculate s


    //if(p.z>(*esp_dc)) p.z=(*esp_dc);//needed? TODO
		//if(p.z<0.) p.z=0.;//needed? TODO

    unsigned int reflected;

		if(new_bulk != p.bulkpos) {
      // If changing voxel do Reflect
      reflected = Reflect(&p,new_bulk,&x,&a,2);
    }

    // Accumulate photon hitting density if needed
    if(fabsf(p.x)<2*(*esp_dc) && fabsf(p.y)<2*(*esp_dc) && p.z<(*esp_dc) && *fhd_activated_dc == 1) {
      // Inside space of 3D matrix. Calculate index of the current voxel
      // Use round to zero so there are no over sampled voxels (for ex: (max_x,0,0) and (0,1,0) should not map to the same voxel)
      index = __float2uint_rz((p.x+2*(*esp_dc))*__int2float_rn(*grid_size_dc))
            + num_x * (__float2uint_rz((p.y+2*(*esp_dc))*__int2float_rn(*grid_size_dc))
            + num_y * __float2uint_rz((p.z)*__int2float_rn(*grid_size_dc)));
      if (DeviceMem.fhd[index] + p.weight < LLONG_MAX) atomicAdd(&DeviceMem.fhd[index], p.weight); // Check for overflow and add atomically
      //TODO why LLONG_MAX?
    }


    if(p.bulkpos == 0){
      // Photon is outside bulk and reflected
      if(fabsf(p.x-det_dc[0].x0)<size_x && fabsf(p.y-det_dc[0].y0)<size_y) {
        // Photon is detectable
        // Use round to zero so there are no over sampled pixels (for ex: (max_x,0) and (0,1) should not map to the same pixel)
        index=__float2uint_rz(__fdividef(p.y-det_dc[0].y0+size_y,det_dc[0].dy)) * det_dc[0].nx +
              __float2uint_rz(__fdividef(p.x-det_dc[0].x0+size_x,det_dc[0].dx));
				if ((DeviceMem.Rd_xy[index] + p.weight) < LLONG_MAX) atomicAdd(&DeviceMem.Rd_xy[index], p.weight); // Check for overflow and add atomicall
				}

      if (*do_temp_sim_dc==1u && det_dc[0].temp_rort==0u && p.tof<det_dc[0].max_temp){
        // Save time value in apropiate bin
        for (int xpos = 0; xpos < num_x_tdet; xpos++){
          for (int ypos = 0; ypos < num_y_tdet; ypos++){
            if ((powf((p.x - DeviceMem.tdet_pos_x[xpos]),2) + powf((p.y - DeviceMem.tdet_pos_y[ypos]),2)) < powf((det_dc[0].temp_det_r),2)){
              // Inside time detector
              index = xpos + num_x_tdet*(ypos + num_y_tdet*__float2uint_rz(__fdividef(p.tof, det_dc[0].max_temp)*num_tbins));
              if (DeviceMem.time_xyt[index] + p.weight < LLONG_MAX) atomicAdd(&DeviceMem.time_xyt[index], p.weight); // Check for overflow and add atomically //TODO why LLONG_MAX?
            }
          }
        }
      }

      p.weight = 0u; // Set the remaining weight to 0, effectively killing the photon
      s = 0.0f;
    }

    if(p.bulkpos == last_bulk){
      // Photon is outside bulk and transmitted
      if(fabsf(p.x-det_dc[0].x0)<size_x && fabsf(p.y-det_dc[0].y0)<size_y) {
        // Photon transmitted
        // Use round to zero so there are no over sampled pixels (for ex: (max_x,0) and (0,1) should not map to the same pixel)
        index=__float2uint_rz(__fdividef(p.y-det_dc[0].y0+size_y,det_dc[0].dy)) * det_dc[0].nx +
              __float2uint_rz(__fdividef(p.x-det_dc[0].x0+size_x,det_dc[0].dx));
        if ((DeviceMem.Tt_xy[index] + p.weight) < LLONG_MAX) atomicAdd(&DeviceMem.Tt_xy[index], p.weight); // Check for overflow and add atomically
        }

      if (*do_temp_sim_dc==1u && det_dc[0].temp_rort==1u && p.tof<det_dc[0].max_temp){
        // Save time value in apropiate bin
        for (int xpos = 0; xpos < num_x_tdet; xpos++){
          for (int ypos = 0; ypos < num_y_tdet; ypos++){
            if ((powf((p.x - DeviceMem.tdet_pos_x[xpos]),2) + powf((p.y - DeviceMem.tdet_pos_y[ypos]),2)) < powf((det_dc[0].temp_det_r),2)){
              // Inside time detector
              index = xpos + num_x_tdet*(ypos + num_y_tdet*__float2uint_rz(__fdividef(p.tof, det_dc[0].max_temp)*num_tbins));
              if (DeviceMem.time_xyt[index] + p.weight < LLONG_MAX) atomicAdd(&DeviceMem.time_xyt[index], p.weight); // Check for overflow and add atomically //TODO why LLONG_MAX?
            }
          }
        }
      }

      p.weight = 0u; // Set the remaining weight to 0, effectively killing the photon
      s = 0.0f;
      }

		w=0;

		if(s > 0.0f) {
    	// Drop weight (apparently only when the photon is scattered)
			w = __float2uint_ru(bulks_dc[p.bulkpos].mua*bulks_dc[p.bulkpos].mutr*__uint2float_rn(p.weight));
			if (p.weight - w >= 0 && w > 0) // Check for underflow
        p.weight -= w;
      else
        p.weight = 0u;
			Spin(&p,bulks_dc[p.bulkpos].g,&x,&a);
		}

		if(!PhotonSurvive(&p,&x,&a)) {
      // Check if photons survives or not
			if(atomicAdd(DeviceMem.num_terminated_photons,1ULL) < (*num_photons_dc-NUM_THREADS)) {
        // Ok to launch another photon
				LaunchPhoton(&p,&x,&a, DeviceMem);//Launch a new photon
      }
			else {
{'fraction_non_alphanumeric': 0.14172695006637948, 'fraction_numerical': 0.015985261047441003, 'mean_word_length': 2.7969344717621643, 'pattern_counts': {'":': 0, '<': 111, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3471911', 'n_tokens_mistral': 14525, 'n_tokens_neox': 13696, 'n_words': 3054}