Individual pixel gain in spe

examples/viewer.xml

@@ -25,6 +25,7 @@
     <fact.extraction.SinglePulseExtraction
         dataKey="DataCalibrated"
         outputKey="PhotonArrivals"
+        url="${integralGainFile}"
     />
 
     <fact.utils.ConvertSinglePulses2Timeseries

src/main/java/fact/extraction/SinglePulseExtraction.java

@@ -8,6 +8,12 @@
 import java.util.Arrays;
 import fact.features.singlePulse.timeSeriesExtraction.SinglePulseExtractor;
 import fact.features.singlePulse.timeSeriesExtraction.ElementWise;
+import fact.features.singlePulse.timeSeriesExtraction.TemplatePulse;
+import org.slf4j.Logger;
+import org.slf4j.LoggerFactory;
+import stream.io.SourceURL;
+import stream.io.CsvStream;
+
 
 /**
  * Extracts a list of arrival slice positions of photons for each pixel
@@ -34,6 +40,9 @@
  * and some modifications from Jens Buss
 */
 public class SinglePulseExtraction implements Processor {
+
+    static Logger log = LoggerFactory.getLogger(SinglePulseExtraction.class);
+
     @Parameter(required=true, description="")
     private String dataKey = null;
 
@@ -83,6 +92,16 @@ public class SinglePulseExtraction implements Processor {
     private int npix = Constants.NUMBEROFPIXEL;
     private int roi  = 300;
 
+    @Parameter(
+        required = false,
+        description = "The url to the inputfiles for the gain calibration constants",
+        defaultValue="classpath:/default/gain_sorted_20131127.csv"
+    )
+    protected SourceURL url = null;
+
+    protected double[] integralGains = null;
+    protected double factSinglePePulseIntegral;
+
     @Override
     public Data process(Data input) {
 
@@ -114,7 +133,10 @@ public Data process(Data input) {
             );
 
             double[] pixelTimeSeries = ElementWise.multiply(
-                pixelTimeSeriesInMv, 1.0/config.factSinglePeAmplitudeInMv);
+                pixelTimeSeriesInMv,
+                factSinglePePulseIntegral /
+                    integralGains[pix]
+            );
 
             SinglePulseExtractor.Result result = spe.extractFromTimeSeries(
                 pixelTimeSeries);
@@ -148,6 +170,42 @@ private void addStartSliceOffset(int[][] arr) {
         }
     }
 
+    public double[] loadIntegralGainFile(SourceURL inputUrl, Logger log) {
+        factSinglePePulseIntegral = TemplatePulse.factSinglePePulseIntegral();
+
+        double[] integralGains = new double[npix];
+        Data integralGainData = null;
+        try {
+            CsvStream stream = new CsvStream(inputUrl, " ");
+            stream.setHeader(false);
+            stream.init();
+            integralGainData = stream.readNext();
+
+            for (int i = 0 ; i < npix ; i++){
+                String key = "column:" + (i);
+                integralGains[i] = (Double) integralGainData.get(key);
+            }
+            return integralGains;
+
+        } catch (Exception e) {
+            log.error("Failed to load integral Gain data: {}", e.getMessage());
+            e.printStackTrace();
+            return null;
+        }
+    }
+
+    public void setUrl(SourceURL url) {
+        try {
+            integralGains = loadIntegralGainFile(url,log);
+        } catch (Exception e) {
+            throw new RuntimeException(e.getMessage());
+        }
+        this.url = url;
+    }
+
+    public SourceURL getUrl() {
+        return url;
+    }
 
 
     public void setDataKey(String dataKey) {

src/main/java/fact/features/singlePulse/timeSeriesExtraction/TemplatePulse.java

@@ -42,21 +42,21 @@ public class TemplatePulse {
     *                      |        |                     /
     *                      0   lengthInSlices
     */
-    public static double[] factSinglePePulse(int lengthInSlices) {
-
-        final double periodeSliceInNs = 0.5;
+    public static double[] factSinglePePulse(int lengthInSlices, double samplingPeriodInNs) {
 
         final double[] time = new double[lengthInSlices];
-        for (int i = 0; i < time.length; i++) {time[i] = i*periodeSliceInNs;}
+        for (int i = 0; i < time.length; i++) {
+            time[i] = i * samplingPeriodInNs;
+        }
 
         double[] template = new double[lengthInSlices];
         for (int i = 0; i < time.length; i++) {
 
-            final double amplitude = 
+            final double amplitude =
                 1.626*
                 (1.0-Math.exp(-0.3803*time[i]))
                 *Math.exp(-0.0649*time[i]);
-            
+
             if(amplitude < 0.0) {
                 template[i] = 0.0;
             }else{
@@ -65,4 +65,70 @@ public static double[] factSinglePePulse(int lengthInSlices) {
         }
         return template;
     }
+
+    public static double[] factSinglePePulse(int lengthInSlices) {
+
+        final double samplingPeriodInNs = 0.5;
+        return factSinglePePulse(lengthInSlices, samplingPeriodInNs);
+
+    }
+
+    public static double factSinglePePulseIntegral() {
+
+        final double half_height = 0.5; // this is fix forever :-D
+        final double samplingPeriodInNs = 1e-3;
+        final int integrationWindowLengthInSlices = 30;
+        final double sliceInNs = 0.5;
+        final int factor = (int)(sliceInNs / samplingPeriodInNs);
+        final double pulseLengthInNs = integrationWindowLengthInSlices/sliceInNs + 3.;
+        final int templateLength = (int) (pulseLengthInNs / samplingPeriodInNs);
+        double[] template = factSinglePePulse(templateLength, samplingPeriodInNs);
+
+        // resample in 500ps steps, i.e. take only every 500th sample
+        double[] sums = new double[factor];
+        double[] resampled_pulse;
+        for (int start_offset = 0; start_offset < factor; start_offset++){
+            resampled_pulse = new double[templateLength/factor];
+            for (int i=0; i<resampled_pulse.length; i++){
+                resampled_pulse[i] = template[start_offset+(i*factor)];
+            }
+
+            sums[start_offset] = integrate(
+                resampled_pulse,
+                findHalfHeightPosition(resampled_pulse, 1.),
+                integrationWindowLengthInSlices
+            );
+        }
+
+        return mean(sums);
+
+    }
+
+    public static int findHalfHeightPosition(double[] samples, double height) {
+        for (int i= 1; i < samples.length; i++){
+            if (samples[i] > height/2.) {
+                return i - 1;
+            }
+        }
+        return samples.length;
+    }
+
+    public static double integrate(double[] samples, int start, int length) {
+        double sum = 0;
+        for (int i=0; i < length; i++) {
+            sum += samples[start+i];
+        }
+        return sum;
+    }
+
+    public static double mean(double[] foo) {
+        if (foo.length == 0) { return Double.NaN; }
+
+        double m = 0.;
+        for (double a: foo){
+            m += a;
+        }
+        return m / foo.length;
+    }
+
 }

src/test/java/fact/utils/TemplatePulseTest.java

@@ -0,0 +1,58 @@
+package fact.utils;
+
+import fact.features.singlePulse.timeSeriesExtraction.TemplatePulse;
+import junit.framework.Assert;
+import org.junit.Test;
+
+public class TemplatePulseTest {
+
+    @Test
+    public void testIfResultNear24() {
+
+        double result = TemplatePulse.factSinglePePulseIntegral();
+        System.out.println(result);
+        Assert.assertTrue(
+            Math.abs(result - 24.37) < 0.01
+        );
+    }
+
+    @Test
+    public void testIntegrateEmptySamples() {
+        Throwable e = null;
+        try {
+            double[] samples = new double[]{};
+            TemplatePulse.integrate(samples, 0, 10);
+        }catch(Throwable ex) {
+            e = ex;
+        }
+        Assert.assertTrue(e instanceof IndexOutOfBoundsException);
+    }
+
+    @Test
+    public void testIntegrate() {
+        double[] samples = new double[]{1,2,3};
+        Assert.assertEquals(
+            6.,
+            TemplatePulse.integrate(samples, 0, 3)
+        );
+    }
+
+    @Test
+    public void testMean() {
+        double[] foo = new double[]{1,2,3};
+        Assert.assertEquals(
+            2.,
+            TemplatePulse.mean(foo)
+        );
+    }
+
+    @Test
+    public void testfindHalfHeightPosition() {
+        double[] foo = new double[]{0.1, 0.3, 0.7, 1, 1, 0.5, 0.3};
+        Assert.assertEquals(
+            1,
+            TemplatePulse.findHalfHeightPosition(foo, 1.)
+        );
+    }
+
+}