From b8e7b443242c716c12006442c2738e09ed77c0c9 Mon Sep 17 00:00:00 2001
From: Stan Seibert <stan@mtrr.org>
Date: Thu, 25 Aug 2011 20:57:36 -0400
Subject: A new PDF evaluation method that does not require storage
 proportional to [number of bins] * [number of PMTs].  Instead we accumulate
 information as the Monte Carlo runs in order to evaluate the PDFs only at the
 points required for the likelihood calculation.

This new interface has been propagated all the way up from the GPU class
through the Simulation class to the Likelihood class.  We have preserved
the full binned histogram implementation in case we need it in the future.
---
 src/daq.cu | 98 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 98 insertions(+)

(limited to 'src')

diff --git a/src/daq.cu b/src/daq.cu
index a34fbbe..cb53d56 100644
--- a/src/daq.cu
+++ b/src/daq.cu
@@ -106,5 +106,103 @@ __global__ void bin_hits(int nchannels,
 		pdf[bin] += 1;
 	}
 }
+	
+__global__ void accumulate_pdf_eval(int time_only,
+				    int nchannels,
+				    unsigned int *event_hit,
+				    float *event_time,
+				    float *event_charge,
+				    float *mc_time,
+				    unsigned int *mc_charge, // quirk of DAQ!
+				    unsigned int *hitcount, 
+				    unsigned int *bincount,
+				    float min_twidth, float tmin, float tmax,
+				    float min_qwidth, float qmin, float qmax,
+				    float *nearest_mc,
+				    int min_bin_content)
+{
+	int id = threadIdx.x + blockDim.x * blockIdx.x;
+	
+	if (id >= nchannels)
+		return;
+	
+	// Was this channel hit in the Monte Carlo?
+	float channel_mc_time = mc_time[id];
+	if (channel_mc_time >= 1e8)
+		return; // Nothing else to do
+	
+	// Is this channel inside the range of the PDF?
+	float distance;
+	int channel_bincount = bincount[id];
+	if (time_only) {
+		if (channel_mc_time < tmin || channel_mc_time > tmax)
+			return;  // Nothing else to do
+		
+		// This MC information is contained in the PDF
+		hitcount[id] += 1;
+	  
+		// Was this channel hit in the event-of-interest?
+		int channel_event_hit = event_hit[id];
+		if (!channel_event_hit)
+			return; // No need to update PDF value for unhit channel
+		
+		// Are we inside the minimum size bin?
+		float channel_event_time = event_time[id];
+		distance = fabsf(channel_mc_time - channel_event_time);
+		if (distance < min_twidth/2.0)
+		  bincount[id] = channel_bincount + 1;
+
+	} else { // time and charge PDF
+		float channel_mc_charge = mc_charge[id]; // int->float conversion because DAQ just returns an integer
+	
+		if (channel_mc_time < tmin || channel_mc_time > tmax ||
+		    channel_mc_charge < qmin || channel_mc_charge > qmax)
+			return;  // Nothing else to do
+		
+		// This MC information is contained in the PDF
+		hitcount[id] += 1;
+	  
+		// Was this channel hit in the event-of-interest?
+		int channel_event_hit = event_hit[id];
+		if (!channel_event_hit)
+			return; // No need to update PDF value for unhit channel
+		
+		// Are we inside the minimum size bin?
+		float channel_event_time = event_time[id];
+		float channel_event_charge = event_charge[id];
+		float normalized_time_distance = fabsf(channel_event_time - channel_mc_time)/min_twidth/2.0;
+		float normalized_charge_distance = fabsf(channel_event_charge - channel_mc_charge)/min_qwidth/2.0;
+		distance = sqrt(normalized_time_distance*normalized_time_distance + normalized_charge_distance*normalized_charge_distance);
+
+		if (distance < 1.0f)
+		  bincount[id] = channel_bincount + 1;
+	}
+
+	// Do we need to keep updating the nearest_mc list?
+	if (channel_bincount + 1 >= min_bin_content)
+		return; // No need to perform insertion into nearest_mc because bincount is a better estimate of the PDF value
+
+	// insertion sort the distance into the array of nearest MC points
+	int offset = min_bin_content * id;
+	int entry = min_bin_content - 1;
+	
+	// If last entry less than new entry, nothing to update
+	if (distance > nearest_mc[offset + entry])
+		return;
+
+	// Find where to insert the new entry while sliding the rest
+	// to the right
+	entry--;
+	while (entry >= 0) {
+		if (nearest_mc[offset+entry] >= distance)
+			nearest_mc[offset+entry+1] = nearest_mc[offset+entry];
+		else
+			break;
+		entry--;
+	}
+
+	nearest_mc[offset+entry+1] = distance;
+}
+
 
 } // extern "C"
-- 
cgit