From 1b00b0b836a4368daa562e368d938cb16047cff6 Mon Sep 17 00:00:00 2001 From: Andreu Date: Wed, 19 May 2021 22:25:57 -0400 Subject: [PATCH] Adding insensitive blocking layers at the top and bottom of the detector. The new parameters BLOCKING_LAYER_BOTTOM and BLOCKING_LAYER_TOP define two blocking layers in the detector slab (8-micron-thick by default, following Zhou et al., Med. Phys. 34, 1098-1109 (2007)) where interactions can happen (and fluorescence generated) but the deposited energy does not contribute to the image. The insensitive top layer causes a measurable drop in DQE(0), but it does not affect MTF as implemented. Pixel values will be reduced (less energy detected per history). Changes marked with keyword: !!BLOCKING_LAYER!! --- MC-GPU_kernel_v1.5b.cu | 22 ++++++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/MC-GPU_kernel_v1.5b.cu b/MC-GPU_kernel_v1.5b.cu index d712606..f9f7b97 100644 --- a/MC-GPU_kernel_v1.5b.cu +++ b/MC-GPU_kernel_v1.5b.cu @@ -37,8 +37,21 @@ // //////////////////////////////////////////////////////////////////////////////// +// ** This software is described in the following reference (please cite it in yuor papers): +// Andreu Badal, Diksha Sharma, Christian G.Graff, Rongping Zeng, and Aldo Badano, Mammography and breast +// tomosynthesis simulator for virtual clinical trials, Computer Physics Communications 261, p. 107779 (2021) +// https://doi.org/10.1016/j.cpc.2020.107779 +// ** Update May 2021 ** !!BLOCKING_LAYER!! +// Enabling blocking (or dead) layers at the top and bottom of the detector slab. +// Interactions in these layers will not contribute to the pixel value, but their fluorescence will +// be tracked and might be detected somewhere else. The insensitive top layer causes a measurable +// drop in DQE(0), but it does not affect MTF as implemented. Pixel values will be reduced (less +// energy detected per history. [Reference: Zhou et al., Med. Phys. 34, 1098-1109 (2007)] +#define BLOCKING_LAYER_TOP 0.0008f // [cm] = 8 micron. Thickness layer closer to source. !!BLOCKING_LAYER!! +#define BLOCKING_LAYER_BOTTOM 0.0008f // [cm] = 8 micron. Thickness layer further from source. !!BLOCKING_LAYER!! + //////////////////////////////////////////////////////////////////////////////// //! Initialize the image array, ie, set all pixels to zero //! Essentially, this function has the same effect as the command: @@ -1495,10 +1508,14 @@ __device__ inline void tally_image(float* energy, float3* position, float3* dire } } - // -- Particle enters the detector! Tally the particle energy in the corresponding pixel (in tenths of meV): + // -- Particle enters the detector! Tally the particle energy in the corresponding pixel (in integer SCALE_eV fractions of eV): // Using a CUDA atomic function (not available for global floats yet) to read and increase the pixel value in a single instruction, blocking interferences from other threads. // The offset for the primaries or scatter images are calculated considering that: // scatter_state=0 for non-scattered, =1 for Compton, =2 for Rayleigh, and =3 for multiple scatter. + + // - Do not count the energy deposited inside the top or bottom blocking (dead) layers of the detector (fluorescence still generated). !!BLOCKING_LAYER!! + if ((position->y > BLOCKING_LAYER_BOTTOM) && (position->y < (detector_data_SHARED->scintillator_thickness-BLOCKING_LAYER_TOP))) // !!BLOCKING_LAYER!! + atomicAdd(( image + // Pointer to beginning of image array (int)(*scatter_state) * detector_data_SHARED->total_num_pixels + // Offset to corresponding scatter image (pixel_coord_x + pixel_coord_z*(detector_data_SHARED->num_pixels.x)) ), // Offset to the corresponding pixel @@ -1520,7 +1537,8 @@ __device__ inline void tally_image(float* energy, float3* position, float3* dire dist_detector = -detector_data_SHARED->fluorescence_MFP*logf(ranecu(seed)); // -- Tally fluorescence energy in the corresponding pixel, unless escaped: position->y = position->y + dist_detector*direction->y; - if ((position->y>0.0f) && (position->yscintillator_thickness)) + + if ((position->y > BLOCKING_LAYER_BOTTOM) && (position->y < (detector_data_SHARED->scintillator_thickness-BLOCKING_LAYER_TOP))) // !!BLOCKING_LAYER!! { position->x = position->x + dist_detector*direction->x; pixel_coord_x = __float2int_rd((position->x - detector_data_SHARED->offset.x + 0.5f*detector_data_SHARED->width_X) * detector_data_SHARED->inv_pixel_size_X); // CUDA intrinsic function converts float to integer rounding down (to minus inf)