visually tested optics code. added models of the inner and outer meshes for the 12" hamamatsu and sno pmts. ratdb.py is able to parse ratdb files. chromaticity.py provides a function to map wavelength -> rgb color. lbne detector model now includes an outer black cylinder and pmts with a glass layer and photocathode/reflective surfaces.

4 files changed, 167 insertions, 170 deletions

diff --git a/src/intersect.h b/src/intersect.h
index e6566f3..bfde782 100644
--- a/src/intersect.h
+++ b/src/intersect.h
@@ -44,7 +44,7 @@ __device__ bool intersect_triangle(const float3 &origin, const float3 &direction
 		    (m.a01*m.a20 - m.a00*m.a21)*b.y +
 		    (m.a00*m.a11 - m.a01*m.a10)*b.z)/determinant;
 
-	if (u3 < 0.0f || (1.0f-u1-u2) < 0.0f)
+	if (u3 <= 0.0f || (1.0f-u1-u2) < 0.0f)
 		return false;
 
 	distance = u3;
@@ -84,32 +84,34 @@ __device__ bool intersect_box(const float3 &origin, const float3 &direction, con
 {
 	float kmin, kmax, kymin, kymax, kzmin, kzmax;
 
-	if (direction.x >= 0.0f)
+	float divx = 1.0f/direction.x;
+	if (divx >= 0.0f)
 	{
-		kmin = (lower_bound.x - origin.x)/direction.x;
-		kmax = (upper_bound.x - origin.x)/direction.x;
+		kmin = (lower_bound.x - origin.x)*divx;
+		kmax = (upper_bound.x - origin.x)*divx;
 	}
 	else
 	{
-		kmin = (upper_bound.x - origin.x)/direction.x;
-		kmax = (lower_bound.x - origin.x)/direction.x;
+		kmin = (upper_bound.x - origin.x)*divx;
+		kmax = (lower_bound.x - origin.x)*divx;
 	}
 
-	if (kmax < kmin)
+	if (kmax < 0.0f)
 		return false;
 
-	if (direction.y >= 0.0f)
+	float divy = 1.0f/direction.y;
+	if (divy >= 0.0f)
 	{
-		kymin = (lower_bound.y - origin.y)/direction.y;
-		kymax = (upper_bound.y - origin.y)/direction.y;
+		kymin = (lower_bound.y - origin.y)*divy;
+		kymax = (upper_bound.y - origin.y)*divy;
 	}
 	else
 	{
-		kymin = (upper_bound.y - origin.y)/direction.y;
-		kymax = (lower_bound.y - origin.y)/direction.y;
+		kymin = (upper_bound.y - origin.y)*divy;
+		kymax = (lower_bound.y - origin.y)*divy;
 	}
 
-	if (kymax < kymin)
+	if (kymax < 0.0f)
 		return false;
 
 	if (kymin > kmin)
@@ -121,18 +123,19 @@ __device__ bool intersect_box(const float3 &origin, const float3 &direction, con
 	if (kmin > kmax)
 		return false;
 
-	if (direction.z >= 0.0f)
+	float divz = 1.0f/direction.z;
+	if (divz >= 0.0f)
 	{
-		kzmin = (lower_bound.z - origin.z)/direction.z;
-		kzmax = (upper_bound.z - origin.z)/direction.z;
+		kzmin = (lower_bound.z - origin.z)*divz;
+		kzmax = (upper_bound.z - origin.z)*divz;
 	}
 	else
 	{
-		kzmin = (upper_bound.z - origin.z)/direction.z;
-		kzmax = (lower_bound.z - origin.z)/direction.z;
+		kzmin = (upper_bound.z - origin.z)*divz;
+		kzmax = (lower_bound.z - origin.z)*divz;
 	}
 
-	if (kzmax < kzmin)
+	if (kzmax < 0.0f)
 		return false;
 
 	if (kzmin > kmin)
@@ -144,9 +147,6 @@ __device__ bool intersect_box(const float3 &origin, const float3 &direction, con
 	if (kmin > kmax)
 		return false;
 
-	if (kmax < 0.0f)
-		return false;
-
 	return true;
 }
 
diff --git a/src/kernel.cu b/src/kernel.cu
index b2e0903..a2d2d71 100644
--- a/src/kernel.cu
+++ b/src/kernel.cu
@@ -11,7 +11,15 @@
 #include "random.h"
 
 #define STACK_SIZE 500
-#define MAX_DEPTH 5
+
+enum
+{
+	DEBUG = -2,
+	INIT = -1,
+	NO_HIT,
+	BULK_ABSORB,
+	SURFACE_ABSORB
+};
 
 /* flattened triangle mesh */
 texture<float4, 1, cudaReadModeElementType> vertices;
@@ -50,11 +58,11 @@ __device__ bool intersect_node(const float3 &origin, const float3 &direction, co
    direction `direction` and the global mesh texture. If the ray intersects
    the texture return the index of the triangle which the ray intersected,
    else return -1. */
-__device__ int intersect_mesh(const float3 &origin, const float3& direction, const int start_node, const int first_node, float &distance)
+__device__ int intersect_mesh(const float3 &origin, const float3& direction, const int start_node, const int first_node, float &min_distance, int last_hit_triangle = -1)
 {
-	int triangle_idx = -1;
+	int triangle_index = -1;
 
-	float min_distance;
+	float distance;
 
 	if (!intersect_node(origin, direction, start_node))
 		return -1;
@@ -76,8 +84,14 @@ __device__ int intersect_mesh(const float3 &origin, const float3& direction, con
 		if (*node >= first_node)
 		{
 			for (i=0; i < length; i++)
+			{
 				if (intersect_node(origin, direction, index+i))
-					*node++ = index+i;
+				{
+					//*node++ = index+i;
+					*node = index+i;
+					node++;
+				}
+			}
 
 			if (node == head)
 				break;
@@ -88,6 +102,9 @@ __device__ int intersect_mesh(const float3 &origin, const float3& direction, con
 		{
 			for (i=0; i < length; i++)
 			{
+				if (last_hit_triangle == index+i)
+					continue;
+
 				uint4 triangle_data = triangles[index+i];
 
 				float3 v0 = make_float3(tex1Dfetch(vertices, triangle_data.x));
@@ -96,16 +113,16 @@ __device__ int intersect_mesh(const float3 &origin, const float3& direction, con
 
 				if (intersect_triangle(origin, direction, v0, v1, v2, distance))
 				{
-					if (triangle_idx == -1)
+					if (triangle_index == -1)
 					{
-						triangle_idx = index + i;
+						triangle_index = index + i;
 						min_distance = distance;
 						continue;
 					}
 
 					if (distance < min_distance)
 					{
-						triangle_idx = index + i;
+						triangle_index = index + i;
 						min_distance = distance;
 					}
 				}
@@ -118,130 +135,125 @@ __device__ int intersect_mesh(const float3 &origin, const float3& direction, con
 	} // while loop
 	while (node != head);
 
-	return triangle_idx;
+	return triangle_index;
 }
 
-__device__ curandState rng_states[256*512];
+__device__ curandState rng_states[1000000];
 
 extern "C"
 {
 
 __global__ void set_pointer(uint4 *triangle_ptr)
 {
-  triangles = triangle_ptr;
+	triangles = triangle_ptr;
 }
 
 /* Initialize random number states */
-__global__ void init_rng(unsigned long long seed, unsigned long long offset)
+__global__ void init_rng(int nthreads, unsigned long long seed, unsigned long long offset)
 {
 	int id = blockIdx.x*blockDim.x + threadIdx.x;
+
+	if (id >= nthreads)
+		return;
+
 	curand_init(seed, id, offset, rng_states+id);
 }
 
 /* Translate `points` by the vector `v` */
 __global__ void translate(int nthreads, float3 *points, float3 v)
 {
-	int idx = blockIdx.x*blockDim.x + threadIdx.x;
+	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
-	if (idx >= nthreads)
+	if (id >= nthreads)
 		return;
 
-	*(points+idx) += v;
+	points[id] += v;
 }
 
 /* Rotate `points` through an angle `phi` counter-clockwise about the
    axis `axis` (when looking towards +infinity). */
 __global__ void rotate(int nthreads, float3 *points, float phi, float3 axis)
 {
-	int idx = blockIdx.x*blockDim.x + threadIdx.x;
+	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
-	if (idx >= nthreads)
+	if (id >= nthreads)
 		return;
 
-	*(points+idx) = rotate(*(points+idx), phi, axis);
+	points[id] = rotate(points[id], phi, axis);
 }
 
-/* Trace the rays starting at `origins` traveling in the direction `directions`
+/* Trace the rays starting at `positions` traveling in the direction `directions`
    to their intersection with the global mesh. If the ray intersects the mesh
    set the pixel associated with the ray to a 32 bit color whose brightness is
    determined by the cosine of the angle between the ray and the normal of the
    triangle it intersected, else set the pixel to 0. */
-__global__ void ray_trace(int nthreads, float3 *origins, float3 *directions, int start_node, int first_node, int *pixels)
+__global__ void ray_trace(int nthreads, float3 *positions, float3 *directions, int start_node, int first_node, int *pixels)
 {
-	int idx = blockIdx.x*blockDim.x + threadIdx.x;
+	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
-	if (idx >= nthreads)
+	if (id >= nthreads)
 		return;
 
-	float3 origin = *(origins+idx);
-	float3 direction = *(directions+idx);
+	float3 position = positions[id];
+	float3 direction = directions[id];
 	direction /= norm(direction);
 
 	float distance;
 
-	int intersection_idx = intersect_mesh(origin, direction, start_node, first_node, distance);
+	int triangle_index = intersect_mesh(position, direction, start_node, first_node, distance);
 
-	if (intersection_idx == -1)
+	if (triangle_index == -1)
 	{
-		*(pixels+idx) = 0;
+		pixels[id] = 0;
 	}
 	else
 	{
-		uint4 triangle_data = triangles[intersection_idx];
+		uint4 triangle_data = triangles[triangle_index];
 
 		float3 v0 = make_float3(tex1Dfetch(vertices, triangle_data.x));
 		float3 v1 = make_float3(tex1Dfetch(vertices, triangle_data.y));
 		float3 v2 = make_float3(tex1Dfetch(vertices, triangle_data.z));
 
-		*(pixels+idx) = get_color(direction, v0, v1, v2, triangle_data.w);
+		pixels[id] = get_color(direction, v0, v1, v2, triangle_data.w);
 	}
+
 } // ray_trace
 
-/* Propagate the photons starting at `positions` traveling in the direction
-   `directions` to their intersection with the global mesh. If the ray
-   intersects the mesh set the hit_solid array value associated with the
-   photon to the triangle index of the triangle the photon intersected, else
-   set the hit_solid array value to -1. */
-__global__ void propagate(int nthreads, float3 *positions, float3 *directions, float *wavelengths, float3 *polarizations, float *times, int *states, int *last_hit_triangles, int start_node, int first_node)
+__global__ void propagate(int nthreads, float3 *positions, float3 *directions, float *wavelengths, float3 *polarizations, float *times, int *states, int *last_hit_triangles, int start_node, int first_node, int max_steps)
 {
 	int id = blockIdx.x*blockDim.x + threadIdx.x;
 
-	//states[id] = interp_property(materials[0].refractive_index
+	if (id >= nthreads)
+		return;
 
-	//return;
+	int state = states[id];
 
-	if (id >= nthreads)
+	if (state != INIT)
 		return;
 
+	curandState rng = rng_states[id];
 	float3 position = positions[id];
 	float3 direction = directions[id];
-	direction /= norm(direction);
 	float3 polarization = polarizations[id];
 	float wavelength = wavelengths[id];
 	float time = times[id];
+	int last_hit_triangle = last_hit_triangles[id];
 
-	int last_hit_triangle;
-
-	curandState rng = rng_states[id];
-
-	unsigned int depth=0;
-	while (true)
+	int steps = 0;
+	while (steps < max_steps)
 	{
-		depth++;
+		steps++;
 
-		if (depth > MAX_DEPTH)
-		{
-			states[id] = MAX_DEPTH_REACHED;
-			break;
-		}
+		direction /= norm(direction);
+		polarization /= norm(polarization);
 
 		float distance;
-		
-	        last_hit_triangle = intersect_mesh(position, direction, start_node, first_node, distance);
-		
+
+	        last_hit_triangle = intersect_mesh(position, direction, start_node, first_node, distance, last_hit_triangle);
+
 		if (last_hit_triangle == -1)
 		{
-			states[id] = NO_HIT;
+			state = NO_HIT;
 			break;
 		}
 
@@ -250,17 +262,11 @@ __global__ void propagate(int nthreads, float3 *positions, float3 *directions, f
 		float3 v0 = make_float3(tex1Dfetch(vertices, triangle_data.x));
 		float3 v1 = make_float3(tex1Dfetch(vertices, triangle_data.y));
 		float3 v2 = make_float3(tex1Dfetch(vertices, triangle_data.z));
-		
-		int material_in_index = 0xFF & (triangle_data.w >> 24);
-		int material_out_index = 0xFF & (triangle_data.w >> 16);
-		int surface_index = 0xFF & (triangle_data.w >> 8);
 
-		if (material_in_index < 0 || material_out_index < 0 || surface_index < 0)
-		{
-			states[id] = -2;
-			break;
-		}
-		
+		int material_in_index = tex1Dfetch(material1_lookup, last_hit_triangle);
+		int material_out_index = tex1Dfetch(material2_lookup, last_hit_triangle);
+		int surface_index = tex1Dfetch(surface_lookup, last_hit_triangle);
+
 		float3 surface_normal = cross(v1-v0,v2-v0);
 		surface_normal /= norm(surface_normal);
 		
@@ -292,9 +298,12 @@ __global__ void propagate(int nthreads, float3 *positions, float3 *directions, f
 		{
 			if (absorption_distance <= distance)
 			{
-				time = absorption_distance/(SPEED_OF_LIGHT/refractive_index1);
-				position = position + absorption_distance*direction;
-				states[id] = BULK_ABSORB;
+				time += absorption_distance/(SPEED_OF_LIGHT/refractive_index1);
+				position += absorption_distance*direction;
+				state = BULK_ABSORB;
+
+				last_hit_triangle = -1;
+
 				break;
 			} // photon is absorbed in material1
 		}
@@ -302,37 +311,48 @@ __global__ void propagate(int nthreads, float3 *positions, float3 *directions, f
 		{
 			if (scattering_distance <= distance)
 			{
-				time = scattering_distance/(SPEED_OF_LIGHT/refractive_index1);
-				position = position + scattering_distance*direction;
+				time += scattering_distance/(SPEED_OF_LIGHT/refractive_index1);
+				position += scattering_distance*direction;
 
-				float x,y;
+				float theta,y;
 				while (true)
 				{
 					y = curand_uniform(&rng);
-					x = uniform(0, 2*PI);
+					theta = uniform(&rng, 0, 2*PI);
 
-					if (y < powf(cosf(x),2))
+					if (y < powf(cosf(theta),2))
 						break;
 				}
 				
-				float phi = uniform(0, 2*PI);
-
-				float3 old_direction = direction;
-				float3 old_polarization = polarization;
+				float phi = uniform(&rng, 0, 2*PI);
 
-				direction = rotate(old_direction, x, cross(old_polarization, old_direction));
-				polarization = rotate(old_polarization, x, cross(old_polarization, old_direction));
+				float3 b = cross(polarization, direction);
+				float3 c = polarization;
 
-				direction = rotate(direction, phi, old_polarization);
-				polarization = rotate(polarization, phi, old_polarization);
+				direction = rotate(direction, theta, b);
+				direction = rotate(direction, phi, c);
+				polarization = rotate(polarization, theta, b);
+				polarization = rotate(polarization, phi, c);
 
-				direction /= norm(direction);
-				polarization /= norm(polarization);
+				last_hit_triangle = -1;
 
 				continue;
 			} // photon is scattered in material1
-		}
-		
+
+		} // if scattering_distance < absorption_distance
+
+		position += distance*direction;
+		time += distance/(SPEED_OF_LIGHT/refractive_index1);
+
+		// p is normal to the plane of incidence
+		float3 p = cross(direction, surface_normal);
+		p /= norm(p);
+
+		float normal_coefficient = dot(polarization, p);
+		float normal_probability = normal_coefficient*normal_coefficient;
+
+		float incident_angle = acosf(dot(surface_normal,-direction));
+
 		if (surface_index != -1)
 		{
 			Surface surface = surfaces[surface_index];
@@ -346,95 +366,81 @@ __global__ void propagate(int nthreads, float3 *positions, float3 *directions, f
 			reflection_diffuse /= sum;
 			reflection_specular /= sum;
 
-			float p = curand_uniform(&rng);
+			float uniform_sample = curand_uniform(&rng);
 
-			if (p < absorption)
+			if (uniform_sample < absorption)
 			{
 				// absorb
-				states[id] = SURFACE_ABSORB;
+				state = SURFACE_ABSORB;
 				break;
 			}
-			else if (p < absorption + reflection_diffuse)
+			else if (uniform_sample < absorption + reflection_diffuse)
 			{
 				// diffusely reflect
-				while (true)
-				{
-					direction = uniform_sphere(&rng);
+				direction = uniform_sphere(&rng);
 
-					if (dot(direction, surface_normal) > 0.0f)
-					{
-						// randomize polarization
-						polarization = cross(uniform_sphere(&rng), direction);
-						break;
-					}
-				}
+				if (dot(direction, surface_normal) < 0.0f)
+					direction = -direction;
+
+				// randomize polarization ?
+				polarization = cross(uniform_sphere(&rng), direction);
+				polarization /= norm(polarization);
 
 				continue;
 			}
 			else
 			{
 				// specularly reflect
-				direction = rotate(direction, -PI/2.0f, cross(direction, surface_normal));
 
-				// polarization ?
+				direction = rotate(surface_normal, incident_angle, p);
 
 				continue;
 			}
-		} // surface
-		
-		// compute reflection/transmission probability
-
-		// p is normal to the plane of incidence
-		float3 p = cross(direction, surface_normal);
 
-		float normal_coefficient = dot(polarization, p);
-		float normal_probability = normal_coefficient*normal_coefficient;
+		} // surface
 
-		float incident_angle = acosf(dot(surface_normal,-direction));
 		float refracted_angle = asinf(sinf(incident_angle)*refractive_index1/refractive_index2);
 
 		if (curand_uniform(&rng) < normal_probability)
 		{
-			// photon polarization normal to plane of incidence
+
 			float reflection_coefficient = -sinf(incident_angle-refracted_angle)/sinf(incident_angle+refracted_angle);
 			float reflection_probability = reflection_coefficient*reflection_coefficient;
 
-			polarization = p;
-
-			if (curand_uniform(&rng) < reflection_probability)
-			{
-				// reflect off surface
-				direction = rotate(surface_normal, PI/2.0f, p);
-			}
+			if ((curand_uniform(&rng) < reflection_probability) || isnan(refracted_angle))
+				direction = rotate(surface_normal, incident_angle, p);
 			else
-			{
-				// transmit
 				direction = rotate(surface_normal, PI-refracted_angle, p);
-			}
+
+			polarization = p;
+
 			continue;
-		}
+		}  // photon polarization normal to plane of incidence
 		else
 		{
-			// photon polarization parallel to surface
 			float reflection_coefficient = tanf(incident_angle-refracted_angle)/tanf(incident_angle+refracted_angle);
 			float reflection_probability = reflection_coefficient*reflection_coefficient;
 
-			if (curand_uniform(&rng) < reflection_probability)
-			{
-				// reflect off surface
-				direction = rotate(surface_normal, PI/2.0f, p);
-				polarization = cross(p, direction);
-			}
+			if ((curand_uniform(&rng) < reflection_probability) || isnan(refracted_angle))
+				direction = rotate(surface_normal, incident_angle, p);
 			else
-			{
-				// transmit
 				direction = rotate(surface_normal, PI-refracted_angle, p);
-				polarization = cross(p, direction);
-			}
-			continue;
-		}
 
-	} // while(true)
+			polarization = cross(p, direction);
+
+			continue;
+		} // photon polarization parallel to surface
+
+	} // while(nsteps < max_steps)
+
+	states[id] = state;
+	rng_states[id] = rng;
+	positions[id] = position;
+	directions[id] = direction;
+	polarizations[id] = polarization;
+	wavelengths[id] = wavelength;
+	times[id] = time;
+	last_hit_triangles[id] = last_hit_triangle;
 
 } // propagate
 
diff --git a/src/materials.h b/src/materials.h
index 47b7d22..77c9f43 100644
--- a/src/materials.h
+++ b/src/materials.h
@@ -6,15 +6,6 @@ __device__ float max_wavelength;
 __device__ float wavelength_step;
 __device__ unsigned int wavelength_size;
 
-enum
-{
-	INIT = -1,
-	MAX_DEPTH_REACHED,
-	NO_HIT,
-	BULK_ABSORB,
-	SURFACE_ABSORB
-};
-
 struct Material
 {
 	float *refractive_index;
@@ -40,10 +31,9 @@ __device__ float interp_property(const float &x, const float *fp)
 	if (x > max_wavelength)
 		return fp[wavelength_size-1];
 
-	unsigned int jl = (x-min_wavelength)/wavelength_step;
-	float xl = min_wavelength + jl*wavelength_step;
+	int jl = (x-min_wavelength)/wavelength_step;
 
-	return xl + (x-xl)*(fp[jl+1]-fp[jl])/wavelength_step;
+	return fp[jl] + (x-(min_wavelength + jl*wavelength_step))*(fp[jl+1]-fp[jl])/wavelength_step;
 }
 
 extern "C"
diff --git a/src/physical_constants.h b/src/physical_constants.h
index 2ff87cd..d591e93 100644
--- a/src/physical_constants.h
+++ b/src/physical_constants.h
@@ -3,5 +3,6 @@
 
 #define SPEED_OF_LIGHT 2.99792458e8
 #define PI 3.141592653589793
+#define EPSILON 1.0e-3
 
 #endif