Author: [log in to unmask]
Date: Sun May 24 14:17:59 2015
New Revision: 3017

Log:
separate class for pulse fitting

Added:
    java/trunk/users/src/main/java/org/hps/users/baltzell/EcalPulseFitter.java

Added: java/trunk/users/src/main/java/org/hps/users/baltzell/EcalPulseFitter.java
 =============================================================================
--- java/trunk/users/src/main/java/org/hps/users/baltzell/EcalPulseFitter.java	(added)
+++ java/trunk/users/src/main/java/org/hps/users/baltzell/EcalPulseFitter.java	Sun May 24 14:17:59 2015
@@ -0,0 +1,156 @@
+package org.hps.users.baltzell;
+
+import hep.aida.IAnalysisFactory;
+import hep.aida.IDataPointSet;
+import hep.aida.IFitFactory;
+import hep.aida.IFitResult;
+import hep.aida.IFitter;
+import hep.aida.IFunction;
+import hep.aida.IFunctionFactory;
+
+import java.io.FileWriter;
+import java.io.IOException;
+import java.util.logging.Level;
+import java.util.logging.Logger;
+
+import org.lcsim.util.aida.AIDA;
+
+public class EcalPulseFitter {
+
+    AIDA aida = AIDA.defaultInstance();
+    IAnalysisFactory analysisFactory = aida.analysisFactory();
+    IFunctionFactory functionFactory = analysisFactory.createFunctionFactory(null);
+    IFitFactory fitFactory = analysisFactory.createFitFactory();
+    IFitter fitter=fitFactory.createFitter();
+    IFunction fitFunction=new Ecal3PoleFunction();
+    IDataPointSet fitData=aida.analysisFactory().createDataPointSetFactory(null).create("ADC DataPointSet", 2);
+
+    private final int nsPerSample=4;
+    public int debug = 0;
+    public String fitFileName = null;
+    public FileWriter fitFileWriter = null;
+    public boolean fixShapeParameter = false;
+    final private double threePoleWidth=2.478;
+    
+    
+    public IFitResult fitPulse(short samples[],int threshCross,double maxADC,double noise) {
+       
+        if (debug>0) System.err.println("FITTING.....................................................");
+     
+        final int threshRange[]={6,25};
+        final int fitRange[]={-10,15};
+        final int pedRange[]={-10,-5};
+       
+        // don't bother with pulses far from trigger:
+        if (threshCross < threshRange[0] || threshCross > threshRange[1]) return null;
+    
+        // calculate pedestal for initializing fit parameters:
+        int nped=0;
+        double ped=0;
+        for (int ii=threshCross+pedRange[0]; ii<threshCross+pedRange[1]; ii++) 
+        {
+            if (ii<0) continue;
+            if (ii>=samples.length) break;
+            ped += samples[ii];
+            nped++;
+        }
+        
+        // don't bother trying to fit:
+        if (nped==0) return null;
+        ped /= nped;
+       
+        // choose points to fit and get starting value for pulse integral:
+        fitData.clear();
+        int nFitPoints=0;
+        int sumADC=0;
+//        int maxADC=0;
+        for (int ii=threshCross+fitRange[0]; ii<threshCross+fitRange[1]; ii++) 
+        {
+            if (ii<0) continue;
+            if (ii>=samples.length) break;
+            if (debug>0) System.err.print(ii+":"+samples[ii]+" ");
+//            if (samples[ii] > maxADC) maxADC=samples[ii];
+            sumADC += samples[ii];
+            fitData.addPoint();
+            fitData.point(nFitPoints).coordinate(0).setValue(ii);
+            fitData.point(nFitPoints).coordinate(1).setValue(samples[ii]);
+            fitData.point(nFitPoints).coordinate(1).setErrorMinus(noise);
+            fitData.point(nFitPoints).coordinate(1).setErrorPlus(noise);
+            nFitPoints++;
+        }
+        if (debug>0) System.err.print("\n");
+       
+        // don't bother trying to fit:
+        if (nFitPoints < 10) return null;
+        if (maxADC < ped) return null;
+        
+        if (debug>0) System.err.println("------- "+ped+" "+threshCross+" "+maxADC);
+
+        final double pulseIntegral = sumADC-ped*nFitPoints;
+ 
+        // initialize parameters:
+        fitFunction.setParameter("pedestal",ped);
+        fitFunction.setParameter("time0",(double)threshCross-2);
+        fitFunction.setParameter("integral",pulseIntegral>0 ? pulseIntegral : 2);
+        fitFunction.setParameter("width",threePoleWidth);
+
+        // constrain parameters:
+        fitter.fitParameterSettings("integral").setBounds(0,999999);
+        fitter.fitParameterSettings("time0").setBounds(1,30);
+        fitter.fitParameterSettings("width").setBounds(0.1,5);
+        if (fixShapeParameter) fitter.fitParameterSettings("width").setFixed(true);
+
+        ((Ecal3PoleFunction)fitFunction).setDebug(debug>1);
+
+        if (debug>0) {
+            System.err.println(String.format("A= %8.2f",fitFunction.parameter("integral")));
+            System.err.println(String.format("T= %8.2f",fitFunction.parameter("time0")*4));
+            System.err.println(String.format("P= %8.2f",fitFunction.parameter("pedestal")));
+            System.err.println(String.format("S= %8.2f",fitFunction.parameter("width")));
+        } else Logger.getLogger("org.freehep.math.minuit").setLevel(Level.OFF);
+
+        IFitResult fitResult = fitter.fit(fitData,fitFunction);
+        writeFit(samples,fitResult);
+       
+        if (debug>0) {
+            final double P = fitResult.fittedParameter("pedestal");
+            final double I = fitResult.fittedParameter("integral");
+            final double T = fitResult.fittedParameter("time0")*nsPerSample;
+            final double S = fitResult.fittedParameter("width");
+            final double Q = fitResult.quality();
+            final double E[] = fitResult.errors();
+
+            if (debug>0) {
+               System.err.println(";;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;");
+               System.err.println(String.format("I = %8.2f %8.2f",I,sumADC-P*30));
+               System.err.println(String.format("T = %8.2f %8.2f",T,threshCross));
+               System.err.println(String.format("P = %8.2f %8.2f",P,ped));
+               System.err.println(String.format("S = %8.2f %8.2f",S,threePoleWidth));
+               System.err.println(String.format("M = %8.2f",maxADC-P));
+            }
+        }
+        
+        return fitResult;
+    }
+
+    
+    public void writeFit(short samples[],IFitResult fit) {
+        if (fitFileName == null) return;
+        if (fitFileWriter == null) {
+            try { fitFileWriter=new FileWriter(fitFileName); }
+            catch (IOException ee) { throw new RuntimeException("Error opening file "+fitFileName,ee); }
+        }
+        try {
+            for (final short ss : samples) fitFileWriter.write(String.format("%6d ",ss));
+            fitFileWriter.write(String.format("%8.3f",fit.fittedParameter("integral")));
+            fitFileWriter.write(String.format("%8.3f",fit.fittedParameter("time0")));
+            fitFileWriter.write(String.format("%8.3f",fit.fittedParameter("pedestal")));
+            fitFileWriter.write(String.format("%8.3f",fit.fittedParameter("width")));
+            fitFileWriter.write(String.format("%8.3f",fit.quality()));
+            fitFileWriter.write("\n");
+        } catch (IOException ee) {
+            throw new RuntimeException("Error writing file "+fitFileName,ee);
+        } 
+    }
+    
+}