Author: [log in to unmask]
Date: Wed May 27 14:10:48 2015
New Revision: 3040
Log:
~final pulse fitting
Added:
java/trunk/steering-files/src/main/resources/org/hps/steering/users/baltzell/EcalReconPulseFit.lcsim
Modified:
java/trunk/users/src/main/java/org/hps/users/baltzell/Ecal3PoleFunction.java
java/trunk/users/src/main/java/org/hps/users/baltzell/EcalPulseFitter.java
java/trunk/users/src/main/java/org/hps/users/baltzell/EcalRawConverter.java
java/trunk/users/src/main/java/org/hps/users/baltzell/EcalRawConverterDriver.java
Added: java/trunk/steering-files/src/main/resources/org/hps/steering/users/baltzell/EcalReconPulseFit.lcsim
=============================================================================
--- java/trunk/steering-files/src/main/resources/org/hps/steering/users/baltzell/EcalReconPulseFit.lcsim (added)
+++ java/trunk/steering-files/src/main/resources/org/hps/steering/users/baltzell/EcalReconPulseFit.lcsim Wed May 27 14:10:48 2015
@@ -0,0 +1,52 @@
+<lcsim xmlns:xs="http://www.w3.org/2001/XMLSchema-instance"
+ xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/lcsim/1.0/lcsim.xsd">
+ <execute>
+ <driver name="EcalRunningPedestal"/>
+ <driver name="EcalRawConverter" />
+ <driver name="ReconClusterer" />
+ <driver name="LCIOWriter" />
+ <!--<driver name="AidaSaveDriver" />-->
+ <driver name="CleanupDriver" />
+ </execute>
+ <drivers>
+ <driver name="EcalRunningPedestal" type="org.hps.recon.ecal.EcalRunningPedestalDriver">
+ <minLookbackEvents>10</minLookbackEvents>
+ <maxLookbackEvents>50</maxLookbackEvents>
+ </driver>
+ <driver name="EcalRawConverter" type="org.hps.users.baltzell.EcalRawConverterDriver">
+ <ecalCollectionName>EcalCalHits</ecalCollectionName>
+ <use2014Gain>false</use2014Gain>
+ <useTimestamps>false</useTimestamps>
+ <useTruthTime>false</useTruthTime>
+ <useRunningPedestal>true</useRunningPedestal>
+ <useTimeWalkCorrection>false</useTimeWalkCorrection>
+ <emulateFirmware>true</emulateFirmware>
+ <emulateMode7>true</emulateMode7>
+ <leadingEdgeThreshold>12</leadingEdgeThreshold>
+ <nsa>20</nsa>
+ <nsb>100</nsb>
+ <useFit>true</useFit>
+ <fixShapeParameter>true</fixShapeParameter>
+ <debug>0</debug>
+ <!--<fitFileName>ecalPulseFits.txt</fitFileName>-->
+ </driver>
+ <driver name="ReconClusterer" type="org.hps.recon.ecal.cluster.ReconClusterDriver">
+ <logLevel>WARNING</logLevel>
+ <outputClusterCollectionName>EcalClusters</outputClusterCollectionName>
+ <hitEnergyThreshold>0.01</hitEnergyThreshold>
+ <seedEnergyThreshold>0.050</seedEnergyThreshold>
+ <clusterEnergyThreshold>0.100</clusterEnergyThreshold>
+ <minTime>0.0</minTime>
+ <timeWindow>25.0</timeWindow>
+ <useTimeCut>true</useTimeCut>
+ <writeRejectedHitCollection>false</writeRejectedHitCollection>
+ </driver>
+ <driver name="LCIOWriter" type="org.lcsim.util.loop.LCIODriver">
+ <outputFilePath>${outputFile}.slcio</outputFilePath>
+ </driver>
+ <driver name="AidaSaveDriver" type="org.lcsim.job.AidaSaveDriver">
+ <outputFileName>ecalPulseFits.root</outputFileName>
+ </driver>
+ <driver name="CleanupDriver" type="org.lcsim.recon.tracking.digitization.sisim.config.ReadoutCleanupDriver" />
+ </drivers>
+</lcsim>
Modified: java/trunk/users/src/main/java/org/hps/users/baltzell/Ecal3PoleFunction.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/baltzell/Ecal3PoleFunction.java (original)
+++ java/trunk/users/src/main/java/org/hps/users/baltzell/Ecal3PoleFunction.java Wed May 27 14:10:48 2015
@@ -2,9 +2,17 @@
import hep.aida.ref.function.AbstractIFunction;
+/*
+ *
+ * "3-Pole Function:" x**2 * exp(-x/tau)
+ *
+ * Here x is time, and we have pedestal and time offsets, and it's normalized:
+ * PEDESTAL + INTEGRAL / WIDTH^3 / 2 * (TIME-TIME0)**2 * exp(-(TIME-TIME0)/WIDTH)
+ *
+ */
public class Ecal3PoleFunction extends AbstractIFunction {
-
+
protected double pedestal=0;
protected double time0=0;
protected double integral=0;
@@ -21,8 +29,9 @@
init(title);
}
- void setDebug(boolean debug) { this.debug=debug; }
-
+ public void setDebug(boolean debug) { this.debug=debug; }
+ public double maximum() { return value(new double[]{time0+2*width}); }
+
@Override
public double value(double[] v) {
final double time = v[0];
@@ -36,20 +45,20 @@
@Override
public void setParameters(double[] pars) throws IllegalArgumentException {
super.setParameters(pars);
- pedestal=pars[0];
- time0=pars[1];
- integral=pars[2];
- width=pars[3];
+ pedestal = pars[0];
+ time0 = pars[1];
+ integral = pars[2];
+ width = pars[3];
if (debug) System.err.println(String.format("%8.2f %8.2f %8.2f %8.2f",pars[0],pars[1],pars[2],pars[3]));
}
@Override
public void setParameter(String key,double value) throws IllegalArgumentException{
super.setParameter(key,value);
- if (key.equals("pedestal")) pedestal=value;
- else if (key.equals("time0")) time0=value;
- else if (key.equals("integral")) integral=value;
- else if (key.equals("width")) width=value;
+ if (key.equals("pedestal")) pedestal = value;
+ else if (key.equals("time0")) time0 = value;
+ else if (key.equals("integral")) integral = value;
+ else if (key.equals("width")) width = value;
}
}
Modified: java/trunk/users/src/main/java/org/hps/users/baltzell/EcalPulseFitter.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/baltzell/EcalPulseFitter.java (original)
+++ java/trunk/users/src/main/java/org/hps/users/baltzell/EcalPulseFitter.java Wed May 27 14:10:48 2015
@@ -13,34 +13,110 @@
import java.util.logging.Level;
import java.util.logging.Logger;
+import org.hps.conditions.database.DatabaseConditionsManager;
+import org.hps.conditions.ecal.EcalChannelConstants;
+import org.hps.conditions.ecal.EcalConditions;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.geometry.Detector;
import org.lcsim.util.aida.AIDA;
public class EcalPulseFitter {
+
+ public int debug = 0;
+ public String fitFileName = null;
+ public boolean fixShapeParameter = true;
+
+ private final static int nsPerSample=4;
+ private final boolean doTimeFit = false;
+ private FileWriter fitFileWriter = null;
+ private EcalConditions ecalConditions = null;
+
+ // don't bother fitting pulses with threshold crossing outside this sample range:
+ private final static int threshRange[]={7,18}; // (28 ns <--> 72 ns)
+
+ // fit sample range relative to threshold crossing:
+ private final static int fitRange[]={-10,15};
+
+ // sample range relative to threshold crossing used to initialize pedestal fit parameter:
+ private final static int pedRange[]={-10,-5};
+
+ private final static double threePoleWidths[]={
+ 2.44057,2.40057,2.53389,2.32095,2.44267,2.43631,2.40292,2.46911,2.36032,2.42132,2.41473,2.43864,2.37710,2.43278,2.49126,2.41739,2.40560,
+ 2.42042,2.51278,2.46731,2.42647,2.35636,2.55308,2.47517,2.38748,2.48324,2.53726,2.54780,2.50771,2.45188,2.35011,2.35110,2.37575,2.37649,
+ 2.41037,2.38817,2.52422,2.47389,2.40672,2.59415,2.35626,2.44359,2.54118,2.42893,2.43100,2.32748,2.34831,2.30641,2.40003,2.38997,2.41643,
+ 2.45717,2.44830,2.48149,2.44339,2.33572,2.42276,2.43731,2.44205,2.37290,2.54975,2.44517,2.39132,2.48821,2.45687,2.39768,2.61674,2.52142,
+ 2.46620,2.46802,2.54764,2.45238,2.48421,2.44521,2.47204,2.50415,2.47346,2.51053,2.44811,2.54895,2.44839,2.51419,2.45706,2.38596,2.46743,
+ 2.49971,2.37904,2.46478,2.43015,2.52088,2.55846,2.50522,2.30973,2.41450,2.40907,2.31343,2.48668,2.46800,2.43559,2.42113,2.47494,2.36971,
+ 2.44022,2.44193,2.42661,2.41918,2.41102,2.51565,2.34875,2.36750,2.37355,2.48763,2.36659,2.42424,2.44321,2.46651,2.52053,2.46172,2.41316,
+ 2.41670,2.41661,2.53472,2.49345,2.46496,2.40006,2.41084,2.44442,2.49113,2.42821,2.56839,2.45630,2.46344,2.42585,2.43730,2.44339,2.44315,
+ 2.48069,2.39892,2.36987,2.40482,2.40103,2.36307,2.46118,2.47136,2.39327,2.46267,2.30668,2.45917,2.46903,2.51727,2.27056,2.46365,2.35228,
+ 2.52188,2.44052,2.41457,2.25591,2.51210,2.43585,2.57410,2.46877,2.37902,2.46975,2.45628,2.44264,2.36588,2.48233,2.46351,2.37512,2.44228,
+ 2.52664,2.51927,2.37348,2.52284,2.52505,2.55755,2.51522,2.44848,2.51430,2.41387,2.42311,2.42880,2.46546,2.40000,2.45992,2.46841,2.43101,
+ 2.22369,2.43262,2.41324,2.45465,2.40752,2.50491,2.45191,2.40109,2.49099,2.37011,2.70111,2.46843,2.44336,2.44675,2.42308,2.51800,2.45660,
+ 2.43406,2.43297,2.37133,2.50718,2.43035,2.42818,2.47709,2.40578,2.39228,2.47350,2.51876,2.39584,2.44114,2.49120,2.48185,2.59246,2.41612,
+ 2.42914,2.51496,2.38785,2.53339,2.30056,2.25807,2.40500,2.53300,2.37439,2.48418,2.33764,2.47934,2.33080,2.48834,2.48318,2.42860,2.28253,
+ 2.44013,2.38272,2.23632,2.37630,2.34564,2.34538,2.44281,2.33952,2.33574,2.40939,2.44992,2.43986,2.45022,2.39336,2.43253,2.38136,2.41941,
+ 2.56685,2.54684,2.43060,2.33258,2.43024,2.48895,2.51698,2.38837,2.45276,2.38518,2.35826,2.38321,2.31616,2.46480,2.53753,2.28027,2.50727,
+ 2.35369,2.29179,2.06776,2.49429,2.35639,2.42628,2.55657,2.51859,2.49579,2.45617,2.41293,2.47083,2.51653,2.47478,2.46354,2.48828,2.34846,
+ 2.34419,2.43172,2.46589,2.45462,2.52146,2.35946,2.38810,2.46027,2.53848,2.48102,2.42701,2.51750,2.57911,2.45136,2.48329,2.39329,2.48209,
+ 2.35948,2.49431,2.36117,2.26121,2.40899,2.51383,2.41642,2.52102,2.39702,2.43949,2.40867,2.38560,2.46874,2.22937,2.40741,2.36007,2.38784,
+ 2.41052,2.42673,2.39476,2.40239,2.46902,2.55278,2.44661,2.54734,2.41863,2.42451,2.40056,2.44307,2.35066,2.46254,2.43270,2.50439,2.49733,
+ 2.55563,2.48589,2.45219,2.44774,2.41116,2.49081,2.43893,2.43731,2.55774,2.42404,2.37911,2.45140,2.47196,2.50533,2.38292,2.47868,2.48262,
+ 2.39845,2.42290,2.54415,2.43109,2.54306,2.49385,2.35113,2.43233,2.50552,2.49532,2.42667,2.53772,2.41398,2.41968,2.59536,2.52077,2.42026,
+ 2.51269,2.42584,2.54930,2.63547,2.42869,2.43348,2.42402,2.38768,2.41903,2.50153,2.45315,2.45472,2.41113,2.44795,2.50849,2.41817,2.44166,
+ 2.61143,2.38567,2.40737,2.39843,2.43711,2.35324,2.48315,2.38388,2.45665,2.30220,2.41293,2.36640,2.38185,2.40913,2.49233,2.46818,2.51527,
+ 2.54177,2.37051,2.37318,2.38405,2.53258,2.58176,2.53786,2.52385,2.40308,2.45279,2.39370,2.50056,2.41076,2.40789,2.47444,2.52857,2.40741,
+ 2.38524,2.39403,2.42442,2.56829,2.56063,2.58340,2.42897,2.58796,2.56262,2.53480,2.43707,2.38258,2.46584,2.40596,2.42590,2.48281,2.47800
+ };
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();
+ IFunction fitFcn3Pole=new Ecal3PoleFunction();
+ IFunction fitFcnLinear=functionFactory.createFunctionFromScript("fitFcnLinear",1,"p0+x[0]*p1","p0,p1","",null);
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) {
+ public IFitResult fitLeadingEdge(RawTrackerHit hit,int threshCross) {
+
+ final short samples[] = hit.getADCValues();
+ final double noise = findChannel(hit.getCellID()).getCalibration().getNoise();
+
+ int imaxADC=0;
+ for (int ii=threshCross; ii<samples.length-1; ii++)
+ {
+ if (ii<0) continue;
+ if (ii>=samples.length) break;
+ if (samples[ii+1] < samples[ii]) {
+ imaxADC=ii;
+ break;
+ }
+ }
+ fitData.clear();
+ int nFitPoints=0;
+ for (int ii=threshCross-1; ii<imaxADC-1; 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 (nFitPoints<2) return null;
+
+ return fitter.fit(fitData,fitFcnLinear);
+ }
+
+ public IFitResult[] fitPulse(RawTrackerHit hit,int threshCross,double maxADC) {
if (debug>0) System.err.println("FITTING.....................................................");
- final int threshRange[]={6,25};
- final int fitRange[]={-10,15};
- final int pedRange[]={-10,-5};
-
+ final short samples[] = hit.getADCValues();
+ final double noise = findChannel(hit.getCellID()).getCalibration().getNoise();
+ final int cid = ecalConditions.getChannelCollection().findGeometric(hit.getCellID()).getChannelId();
+
// don't bother with pulses far from trigger:
if (threshCross < threshRange[0] || threshCross > threshRange[1]) return null;
@@ -62,6 +138,7 @@
// choose points to fit and get starting value for pulse integral:
fitData.clear();
int nFitPoints=0;
+ int iPeakADC=-1;
int sumADC=0;
// int maxADC=0;
for (int ii=threshCross+fitRange[0]; ii<threshCross+fitRange[1]; ii++)
@@ -70,6 +147,9 @@
if (ii>=samples.length) break;
if (debug>0) System.err.print(ii+":"+samples[ii]+" ");
// if (samples[ii] > maxADC) maxADC=samples[ii];
+ if (iPeakADC<0 && ii<samples.length-1 && samples[ii+1] < samples[ii]) {
+ iPeakADC=ii;
+ }
sumADC += samples[ii];
fitData.addPoint();
fitData.point(nFitPoints).coordinate(0).setValue(ii);
@@ -89,28 +169,29 @@
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);
+ fitFcn3Pole.setParameter("pedestal",ped);
+ fitFcn3Pole.setParameter("time0",(double)threshCross-2);
+ fitFcn3Pole.setParameter("integral",pulseIntegral>0 ? pulseIntegral : 2);
+ fitFcn3Pole.setParameter("width",threePoleWidths[cid-1]);
// constrain parameters:
fitter.fitParameterSettings("integral").setBounds(0,999999);
fitter.fitParameterSettings("time0").setBounds(1,30);
+// fitter.fitParameterSettings("time0").setBounds(threshRange[0]-5,threshRange[1]+20);
fitter.fitParameterSettings("width").setBounds(0.1,5);
if (fixShapeParameter) fitter.fitParameterSettings("width").setFixed(true);
- ((Ecal3PoleFunction)fitFunction).setDebug(debug>1);
+ ((Ecal3PoleFunction)fitFcn3Pole).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")));
+ System.err.println(String.format("A= %8.2f",fitFcn3Pole.parameter("integral")));
+ System.err.println(String.format("T= %8.2f",fitFcn3Pole.parameter("time0")*4));
+ System.err.println(String.format("P= %8.2f",fitFcn3Pole.parameter("pedestal")));
+ System.err.println(String.format("S= %8.2f",fitFcn3Pole.parameter("width")));
} else Logger.getLogger("org.freehep.math.minuit").setLevel(Level.OFF);
- IFitResult fitResult = fitter.fit(fitData,fitFunction);
- writeFit(samples,fitResult);
+ IFitResult fitResult = fitter.fit(fitData,fitFcn3Pole);
+ writeFit(samples,fitResult,cid);
if (debug>0) {
final double P = fitResult.fittedParameter("pedestal");
@@ -119,22 +200,31 @@
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));
+ 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,threePoleWidths[cid-1]));
+ System.err.println(String.format("M = %8.2f",maxADC-P));
+ System.err.println(String.format("Q = %8.2f",Q));
+ }
+
+ // now fit just leading edge:
+ IFitResult timeResult = null;
+ if (doTimeFit) {
+ for (int ii=fitData.size()-1; ii>=0; ii--) {
+ if (fitData.point(ii).coordinate(0).value() > iPeakADC+1) fitData.removePoint(ii);
+ else break;
}
- }
-
- return fitResult;
- }
-
-
- public void writeFit(short samples[],IFitResult fit) {
+ fitter.fit(fitData,fitFcn3Pole);
+ }
+
+ IFitResult fitResults[] = {fitResult,timeResult};
+ return fitResults;
+ }
+
+
+ public void writeFit(short samples[],IFitResult fit,final int cid) {
if (fitFileName == null) return;
if (fitFileWriter == null) {
try { fitFileWriter=new FileWriter(fitFileName); }
@@ -147,10 +237,31 @@
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(String.format("%4d",cid));
fitFileWriter.write("\n");
} catch (IOException ee) {
throw new RuntimeException("Error writing file "+fitFileName,ee);
}
}
-
-}
+
+ public double getSmartPedestal(short samples[]) {
+ double ped = 99999;
+ final int nSamples = 6;
+ for (int ii=0; ii<samples.length-nSamples; ii++) {
+ double aped = 0;
+ for (int jj=ii; jj<ii+nSamples; jj++) aped += samples[jj];
+ aped /= nSamples;
+ if (aped < ped) ped = aped;
+ }
+ return ped;
+ }
+ public void setDetector(Detector detector) {
+ // ECAL combined conditions object.
+ ecalConditions = DatabaseConditionsManager.getInstance().getEcalConditions();
+ }
+ public EcalChannelConstants findChannel(long cellID) {
+ return ecalConditions.getChannelConstants(ecalConditions.getChannelCollection().findGeometric(cellID));
+ }
+
+
+}
Modified: java/trunk/users/src/main/java/org/hps/users/baltzell/EcalRawConverter.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/baltzell/EcalRawConverter.java (original)
+++ java/trunk/users/src/main/java/org/hps/users/baltzell/EcalRawConverter.java Wed May 27 14:10:48 2015
@@ -1,20 +1,12 @@
package org.hps.users.baltzell;
-import org.hps.users.baltzell.Ecal3PoleFunction;
-
-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 hep.aida.IHistogram1D;
+//import hep.aida.IHistogram1D;
+//import hep.aida.IHistogram2D;
import java.awt.event.ActionEvent;
import java.awt.event.ActionListener;
import java.io.FileWriter;
-import java.io.IOException;
import java.util.ArrayList;
import java.util.Map;
import java.util.logging.Level;
@@ -49,7 +41,9 @@
* @author Nathan Baltzell <[log in to unmask]>
* @author Holly Szumila <[log in to unmask]>
*
- * baltzell: New in 2015: (default behavior is still unchanged)
+ * baltzell: New in May, 2015: Pulse Fitting.
+ *
+ * baltzell: New in Early-2015: (default behavior is still unchanged)
*
* Implemented conversion of Mode-1 to Mode-3.
*
@@ -67,18 +61,22 @@
*/
public class EcalRawConverter {
-
- // wtf.
+ EcalPulseFitter pulseFitter = new 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 data=aida.analysisFactory().createDataPointSetFactory(null).create("ADC DataPointSet", 2);
-
- IHistogram1D hWidth = aida.histogram1D("hShape",500,0,5);
- IHistogram1D hQuality = aida.histogram1D("hChiSquare",500,0,5);
+ IHistogram1D hWidth = aida.histogram1D("hWidth",500,0,5);
+ IHistogram2D hWidthE = aida.histogram2D("hWidthE",500,0,5,100,0,5000.0);
+ IHistogram1D hChiSquare = aida.histogram1D("hChiSquare",500,0,5);
+ IHistogram2D hWidthChan = aida.histogram2D("hWidthChan",500,0,5,442,0.5,442.5);
+ IHistogram2D hSamplesChanW = aida.histogram2D("hSamplesChanW",442,0.5,442.5,50,-0.5,49.5);
+ IHistogram2D hSamplesChan = aida.histogram2D("hSamplesChan",442,0.5,442.5,50,-0.5,49.5);
+ IHistogram2D hSamplesChanLoW = aida.histogram2D("hSamplesChanLoW",442,0.5,442.5,50,-0.5,49.5);
+ IHistogram2D hSamplesChanLo = aida.histogram2D("hSamplesChanLo",442,0.5,442.5,50,-0.5,49.5);
+ IHistogram2D hSamplesChanHiW = aida.histogram2D("hSamplesChanHiW",442,0.5,442.5,50,-0.5,49.5);
+ IHistogram2D hSamplesChanHi = aida.histogram2D("hSamplesChanHi",442,0.5,442.5,50,-0.5,49.5);
+ IHistogram2D hWidthEnergy153 = aida.histogram2D("hWidthEnergy153",500,0,5,500,0,5000);
+ */
private boolean useTimeWalkCorrection = false;
private boolean useRunningPedestal = false;
@@ -93,6 +91,11 @@
public String fitFileName = null;
public FileWriter fitFileWriter = null;
public boolean fixShapeParameter = false;
+ //final private double threePoleWidth=2.478;
+
+ public void setUseFit(boolean useFit) { this.useFit=useFit; }
+ public void setFitFileName(String name) { pulseFitter.fitFileName=name; }
+ public void setFixShapeParameter(boolean fix) { pulseFitter.fixShapeParameter=fix; }
/*
* The time for one FADC sample (units = ns).
@@ -322,89 +325,6 @@
return (lastSample-firstSample+1)*getSingleSamplePedestal(event,cellID);
}
- public IFitResult fitPulse(short samples[],int thresholdCrossing,double maxADC,double noise) {
-
- if (thresholdCrossing < 9 || thresholdCrossing > 15) return null;
-
- if (debug>0) System.err.println("FITTING.....................................................");
-
- data.clear();
- int nped=0;
- double ped=0;
- for (int ii=thresholdCrossing-9; ii<thresholdCrossing-5; ii++) {
- if (ii<0 || ii>=samples.length) break;
- ped += samples[ii];
- nped++;
- }
- ped /= nped;
- if (nped==0) return null;
- int npts=0;
- int sumADC=0;
- for (int ii=thresholdCrossing-10; ii<thresholdCrossing+15; ii++) {
- if (ii<0) continue;
- if (ii>=samples.length) break;
- if (debug>0) System.err.print(ii+":"+samples[ii]+" ");
- data.addPoint();
- data.point(npts).coordinate(0).setValue(ii);
- data.point(npts).coordinate(1).setValue(samples[ii]);
- data.point(npts).coordinate(1).setErrorMinus(noise);
- data.point(npts).coordinate(1).setErrorPlus(noise);
- sumADC += samples[ii];
- npts++;
- }
- if (debug>0) System.err.print("\n");
- if (npts>=10) {
- if (debug>0) System.err.println("------- "+ped+" "+thresholdCrossing+" "+maxADC);
- if (maxADC-ped < 0) return null;
- final double pulseIntegral = sumADC-ped*npts;
- fitFunction.setParameter("pedestal",ped);
- fitFunction.setParameter("time0",(double)thresholdCrossing-2);
- fitFunction.setParameter("integral",pulseIntegral>0?pulseIntegral:10);
- fitFunction.setParameter("width",2.478);
- ((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")));
- }
-
- fitter.fitParameterSettings("integral").setBounds(0,999999);
- fitter.fitParameterSettings("time0").setBounds(5,40);
- fitter.fitParameterSettings("width").setBounds(0.1,5);
- if (fixShapeParameter) fitter.fitParameterSettings("width").setFixed(true);
-
- if (debug<1) Logger.getLogger("org.freehep.math.minuit").setLevel(Level.OFF);
- IFitResult fitResult = fitter.fit(data,fitFunction);
- writeFit(samples,fitResult);
- return fitResult;
- }
- return null;
- }
-
- 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);
- }
- }
-
-
-
/*
* Emulate the FADC250 firmware in conversion of Mode-1 waveform to a Mode-3/7 pulse,
* given a time for threshold crossing.
@@ -490,9 +410,13 @@
final int a0 = samples[t0];
final int a1 = samples[t1];
final double slope = (a1 - a0); // units = ADC/sample
- final double yint = a1 - slope * t1; // units = ADC
+ final double yint = a1 - slope * t1; // units = ADC
+
+ // mode-7 time, time at half max:
pulseTime = ((halfMax - a0)/(a1-a0) + t0)* nsPerSample;
-
+
+ // mode-8 time, time at pedestal intersection:
+ //pulseTime = (minADC-yint)/slope * nsPerSample;
}
}
@@ -500,40 +424,50 @@
if (useFit)
{
- final double noise=findChannel(hit.getCellID()).getCalibration().getNoise();
- IFitResult fitResult=fitPulse(samples,thresholdCrossing,maxADC,noise);
-
- if (fitResult != null) {
-
- final double P=fitResult.fittedParameter("pedestal");
- final double A=fitResult.fittedParameter("integral");
- final double T=fitResult.fittedParameter("time0");
- final double S=fitResult.fittedParameter("width");
- final double Q=fitResult.quality();
-
- // finite integral from T to T+NSAMP:
- //final int NSAMP=30;
- //final double I = (A/2/S) * S*S*S* ( 2 - Math.exp(-NSAMP/S) * (Math.pow((S+NSAMP)/S,2) + 1));
-
- // infinite integral:
- //final double I = A*S*S;
+ IFitResult fitResults[] = pulseFitter.fitPulse(hit,thresholdCrossing,maxADC);
+ if (fitResults!=null && fitResults[0]!=null) {
- // after normalizing the function:
- final double I = A;
+ IFitResult fitResult=fitResults[0];
+ IFitResult timeResult=fitResults[1]==null?fitResults[0]:fitResults[1];
+
+ pulseTime = timeResult.fittedParameter("time0")*nsPerSample;
+ sumADC = fitResult.fittedParameter("integral");
+ minADC = fitResult.fittedParameter("pedestal");
+ maxADC = ((Ecal3PoleFunction)fitResult.fittedFunction()).maximum();
- if (debug>0) {
- System.err.println(";;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;");
- System.err.println(String.format("A= %8.2f %8.2f",I,sumADC-P*30));
- System.err.println(String.format("T= %8.2f %8.2f",T*4,pulseTime));
- System.err.println(String.format("P= %8.2f %8.2f",P,minADC));
- System.err.println(String.format("S= %8.2f",S));
- System.err.println(String.format("M= %8.2f %8.2f",A,(maxADC-P)));
- }
-
- hWidth.fill(S);
- hQuality.fill(Q);
- pulseTime = T*4;
- sumADC = I;
+ /*
+ final double width = fitResult.fittedParameter("width");
+ //final double E[] = fitResult.errors();
+ final int cid = getChannelId(hit.getCellID());
+ hWidth.fill(width);
+ hWidthE.fill(width,fitResult.fittedParameter("integral"));
+ if (cid==153) hWidthEnergy153.fill(width,fitResult.fittedParameter("integral"));
+ hChiSquare.fill(fitResult.quality());
+ hWidthChan.fill(width,cid);
+ for (int ii=0; ii<samples.length; ii++) {
+ hSamplesChan.fill(cid,ii);
+ hSamplesChanW.fill(cid,ii,samples[ii]);
+ if (fitResult.fittedParameter("integral") > 2000)
+ {
+ hSamplesChanHi.fill(cid,ii);
+ hSamplesChanHiW.fill(cid,ii,samples[ii]);
+ }
+ else if (fitResult.fittedParameter("integral") < 1000)
+ {
+ hSamplesChanLo.fill(cid,ii);
+ hSamplesChanLoW.fill(cid,ii,samples[ii]);
+ }
+ }
+ */
+
+ /*
+ timeResult = pulseFitter.fitLeadingEdge(hit,thresholdCrossing);
+ if (timeResult != null) {
+ final double p0 = timeResult.fittedParameter("p0");
+ final double p1 = timeResult.fittedParameter("p1");
+ pulseTime = (minADC-p0)/p1 * nsPerSample;
+ }
+ */
}
}
return new double []{pulseTime,sumADC,minADC,maxADC};
@@ -711,6 +645,8 @@
public void setDetector(Detector detector) {
// ECAL combined conditions object.
ecalConditions = DatabaseConditionsManager.getInstance().getEcalConditions();
+
+ pulseFitter.setDetector(detector);
}
/**
@@ -723,4 +659,8 @@
return ecalConditions.getChannelConstants(ecalConditions.getChannelCollection().findGeometric(cellID));
}
+ public int getChannelId(long cellID) {
+ return ecalConditions.getChannelCollection().findGeometric(cellID).getChannelId();
+ }
+
}
Modified: java/trunk/users/src/main/java/org/hps/users/baltzell/EcalRawConverterDriver.java
=============================================================================
--- java/trunk/users/src/main/java/org/hps/users/baltzell/EcalRawConverterDriver.java (original)
+++ java/trunk/users/src/main/java/org/hps/users/baltzell/EcalRawConverterDriver.java Wed May 27 14:10:48 2015
@@ -57,9 +57,9 @@
converter = new EcalRawConverter();
}
- public void setUseFit(boolean useFit) { converter.useFit=useFit; }
- public void setFitFileName(String name) { converter.fitFileName=name; }
- public void setFixShapeParameter(boolean fix) { converter.fixShapeParameter=fix; }
+ public void setUseFit(boolean useFit) { converter.setUseFit(useFit); }
+ public void setFitFileName(String name) { converter.setFitFileName(name); }
+ public void setFixShapeParameter(boolean fix) { converter.setFixShapeParameter(fix); }
/**
