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.

1 files changed, 98 insertions, 0 deletions

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"