Kernel estimation implementation for calculating PDFs at each PMT.

A little rough around the edges, and still needs some development work.

1 files changed, 139 insertions, 0 deletions

diff --git a/chroma/cuda/pdf.cu b/chroma/cuda/pdf.cu
index c9c2ea4..9f547d0 100644
--- a/chroma/cuda/pdf.cu
+++ b/chroma/cuda/pdf.cu
@@ -122,4 +122,143 @@ accumulate_pdf_eval(int time_only, int nchannels, unsigned int *event_hit,
     nearest_mc[offset+entry+1] = distance;
 }
 
+
+__global__ void
+accumulate_moments(int time_only, int nchannels,
+		   float *mc_time,
+		   unsigned int *mc_charge, // quirk of DAQ!
+		   float tmin, float tmax,
+		   float qmin, float qmax,
+		   unsigned int *mom0,
+		   float *t_mom1,
+		   float *t_mom2,
+		   float *q_mom1,
+		   float *q_mom2)
+
+{
+    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];
+
+    // Is this channel inside the range of the PDF?
+    if (time_only) {
+	if (channel_mc_time < tmin || channel_mc_time > tmax)
+	    return;  // Nothing else to do
+	
+	mom0[id] += 1;
+	t_mom1[id] += channel_mc_time;
+	t_mom2[id] += channel_mc_time*channel_mc_time;
+    }
+    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
+
+	mom0[id] += 1;
+	t_mom1[id] += channel_mc_time;
+	t_mom2[id] += channel_mc_time*channel_mc_time;
+	q_mom1[id] += channel_mc_charge;
+	q_mom2[id] += channel_mc_charge*channel_mc_charge;
+    }
+}
+
+static const float invroot2 = 0.70710678118654746f; // 1/sqrt(2)
+static const float rootPiBy2 = 1.2533141373155001f; // sqrt(M_PI/2)
+
+__global__ void
+accumulate_kernel_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!
+		       float tmin, float tmax,
+		       float qmin, float qmax,
+		       float *inv_time_bandwidths,
+		       float *inv_charge_bandwidths,
+		       unsigned int *hitcount,
+		       float *pdf_values)
+{
+    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];
+	
+    // Is this channel inside the range of the PDF?
+    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
+
+	// Kernel argument
+	float channel_event_time = event_time[id];
+	float inv_bandwidth = inv_time_bandwidths[id];
+	float arg = (channel_mc_time - channel_event_time) * inv_bandwidth;
+
+	// evaluate 1D Gaussian normalized within time window
+	float term = expf(-0.5f * arg * arg) * inv_bandwidth;
+
+	float norm = tmax - tmin;
+	if (inv_bandwidth > 0.0f) {
+	  float loarg = (tmin - channel_mc_time)*inv_bandwidth*invroot2;
+	  float hiarg = (tmax - channel_mc_time)*inv_bandwidth*invroot2;
+	  norm = (erff(hiarg) - erff(loarg)) * rootPiBy2;
+	}
+	pdf_values[id] += term / norm;
+    }
+    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
+
+
+	// Kernel argument: time dim
+	float channel_event_obs = event_time[id];
+	float inv_bandwidth = inv_time_bandwidths[id];
+	float arg = (channel_mc_time - channel_event_obs) * inv_bandwidth;
+
+	float loarg = (tmin - channel_mc_time)*inv_bandwidth*invroot2;
+	float hiarg = (tmax - channel_mc_time)*inv_bandwidth*invroot2;
+	float norm = (erff(hiarg) - erff(loarg)) * rootPiBy2;
+
+	// Kernel argument: charge dim
+	channel_event_obs = event_time[id];
+	inv_bandwidth = inv_time_bandwidths[id];
+	arg *= (channel_mc_charge - channel_event_obs) * inv_bandwidth;
+
+	loarg = (tmin - channel_mc_charge)*inv_bandwidth*invroot2;
+	hiarg = (tmax - channel_mc_charge)*inv_bandwidth*invroot2;
+	norm *= (erff(hiarg) - erff(loarg)) * rootPiBy2;
+
+	// evaluate 2D Gaussian normalized within time window
+	float term = expf(-0.5f * arg * arg);
+
+	pdf_values[id] += term / norm;
+    }
+}
+
+
 } // extern "C"