diff --git a/common-tools/clas-detector/src/main/java/org/jlab/detector/pulse/ModeAHDC.java b/common-tools/clas-detector/src/main/java/org/jlab/detector/pulse/ModeAHDC.java index 7cd68df1b1..a1beb05060 100644 --- a/common-tools/clas-detector/src/main/java/org/jlab/detector/pulse/ModeAHDC.java +++ b/common-tools/clas-detector/src/main/java/org/jlab/detector/pulse/ModeAHDC.java @@ -10,23 +10,23 @@ /** * A new extraction method dedicated to the AHDC signal waveform - * + * * Some blocks of code are inspired by MVTFitter.java and Bonus12 (`createBonusBank()`) - * + * * To do list: * - read pedestal from the DB when a baseline cannot be computed * - change the definition of ADC to be the max of the peak? * - fit waveforms to define arrival time (and charge) - * + * * @author ftouchte, pilleux */ public class ModeAHDC extends HipoExtractor { - + //Parameters, to be read from DB? //Saturation threshold should be 4095 (2^12-1) //But in practice waveforms saturate below it - private final short ADC_LIMIT = 3500; + private final short ADC_LIMIT = 3500; //number of consecutive samples exceeding threshold to consider a saturation plateau private final int consecutiveSaturatingSamples = 3; //Sampling time in ns @@ -44,18 +44,7 @@ public class ModeAHDC extends HipoExtractor { private final float defaultBaseline = 300; // dream clock for fine timestamp correction private final float dream_clock = 8.0f; - - //Waveform and corresponding pulse - //This is the CURRENT pulse, it is initialized - //at every extraction. This is not the cleanest - //but it avoids feeding the same arguments to all methods - private short[] samples; - private Pulse pulse; - private long time_ZS; - private int binMax; - private int numberOfBins; - private int effectiveNumberOfBins; - + //Waveform types: //wfType 6 ⇒ too small (nsamples <= 10) //wfType 5 ⇒ decreasing baseline (or leadingEdgeTime fails) @@ -65,83 +54,91 @@ public class ModeAHDC extends HipoExtractor { //wfType 1 ⇒ saturing //wfType 0 ⇒ OK + private static class AHDCPulse extends Pulse { + private long time_ZS; + private int binMax; + private int numberOfBins; + private int effectiveNumberOfBins; + private short[] samples; + } + /** * */ - public void preProcess() { - assignValidType(0); - this.numberOfBins = samples.length; + public void preProcess(AHDCPulse p) { + assignValidType(0, p); + p.numberOfBins = p.samples.length; // check if the signal is too short int countNotZeros = 0; boolean FirstNonZeroFromTheTailReached = false; - for (int i = samples.length-1; i >= 0; i--) { - if (samples[i] > 0) { + for (int i = p.samples.length-1; i >= 0; i--) { + if (p.samples[i] > 0) { countNotZeros++; if (!FirstNonZeroFromTheTailReached) { FirstNonZeroFromTheTailReached = true; - this.effectiveNumberOfBins = i+1; // this value can be identified as the effective number of bins + p.effectiveNumberOfBins = i+1; // this value can be identified as the effective number of bins } } } - if (countNotZeros <= 10) assignInvalidType(); + if (countNotZeros <= 10) assignInvalidType(p); } /** - * This method computes the waveform baseline - * from the average of the first few samples - */ - public void baselineComputation() + * This method computes the waveform baseline + * from the average of the first few samples + */ + public void baselineComputation(AHDCPulse p) { //Default baseline float adcOffset = -99; - - //Check if the wf has sufficient length, if not it is invalid - if (numberOfBins >= this.nBaselineSamples+1){ + //Check if the wf has sufficient length, if not it is invalid + if (p.numberOfBins >= this.nBaselineSamples+1){ + //The baseline is the average of the first few samples for(int i=0; ithis.pulse.wftype) this.pulse.wftype = (short) type; + if(type>p.wftype) p.wftype = (short) type; } /** - * This method assigns an invalid type (-1) to the waveform - */ - public void assignInvalidType(){ - this.pulse.wftype = 6; + * This method assigns an invalid type (-1) to the waveform + */ + public void assignInvalidType(AHDCPulse p){ + p.wftype = 6; } /** - * This method subtracts the baseline, - * computes the max ADC and corresponding time - * and the integral of the wf - * - * @return an int for status - * - */ - public int waveformADCProcessing() + * This method subtracts the baseline, + * computes the max ADC and corresponding time + * and the integral of the wf + * + * @return an int for status + * + */ + public int waveformADCProcessing(AHDCPulse p) { - //int numberOfBins = this.samples.length; + //int numberOfBins = this.samples.length; //Initialize to dummy values - this.pulse.adcMax = -99; - this.binMax = -99; + p.adcMax = -99; + p.binMax = -99; //For invalid wf, dummy values kept //if(this.pulse.wftype == 6) return 0; @@ -154,20 +151,20 @@ public int waveformADCProcessing() int endPlateau = -99; //Looping through samples - for (int bin = 0; bin < numberOfBins; bin++){ + for (int bin = 0; bin < p.numberOfBins; bin++){ //Baseline subtraction - this.samples[bin] = (short) Math.max(this.samples[bin] - this.pulse.pedestal, 0); + p.samples[bin] = (short) Math.max(p.samples[bin] - p.pedestal, 0); //Look for maximum - if ((this.pulse.adcMax < this.samples[bin]) && (bin >= nBaselineSamples)) { // the max should not be in the baseline - this.pulse.adcMax = this.samples[bin]; - this.binMax = bin; + if ((p.adcMax < p.samples[bin]) && (bin >= nBaselineSamples)) { // the max should not be in the baseline + p.adcMax = p.samples[bin]; + p.binMax = bin; } //Integration - this.pulse.integral += this.samples[bin]; + p.integral += p.samples[bin]; //Check for saturation: if the sample exceeds the ths //and if it belongs to a group of saturating samples - if(this.samples[bin]>(this.ADC_LIMIT-this.pulse.pedestal)){ + if(p.samples[bin]>(this.ADC_LIMIT-p.pedestal)){ //if this is the first sample saturating, start a new plateau if(nsaturating==0) startPlateau = bin; //else continue the plateau @@ -176,83 +173,83 @@ public int waveformADCProcessing() //we store the maximum length of consecutive saturating samples if(nsaturating>maxNsaturating) maxNsaturating = nsaturating; } - else nsaturating = 0;//this sample is not in a row of saturating samples + else nsaturating = 0;//this sample is not in a row of saturating samples } - - //Saturating samples: if the number of consecutive saturating samples + + //Saturating samples: if the number of consecutive saturating samples //is greater than the limit, we consider the wf saturated if(maxNsaturating>=this.consecutiveSaturatingSamples){ - this.assignValidType(1); + this.assignValidType(1, p); //The time is taken as the middle of the plateau - this.binMax = (startPlateau + endPlateau)/2; + p.binMax = (startPlateau + endPlateau)/2; } - - //Define the pulse time as the peak time - this.pulse.time = (this.binMax + this.time_ZS)*this.samplingTime; + + //Define the pulse time as the peak time + p.time = (p.binMax + p.time_ZS)*this.samplingTime; //If there are 2*npts+1 points around the peak //(if the peak is not at an edge of the window) //Then the peak ADC value is revisited to be the average of these //TO DO: just use the adcmax??? int npts = 1; - if (this.pulse.adcMax == 0) { assignValidType(5);} // classification before the fit - if ((this.binMax - npts >= 0) && (this.binMax + npts <= numberOfBins - 1)){ - this.pulse.adcMax = 0; - for (int bin = this.binMax - npts; bin <= this.binMax + npts; bin++) this.pulse.adcMax += this.samples[bin]; - this.pulse.adcMax = this.pulse.adcMax/(2*npts+1); + if (p.adcMax == 0) { assignValidType(5, p);} // classification before the fit + if ((p.binMax - npts >= 0) && (p.binMax + npts <= p.numberOfBins - 1)){ + p.adcMax = 0; + for (int bin = p.binMax - npts; bin <= p.binMax + npts; bin++) p.adcMax += p.samples[bin]; + p.adcMax = p.adcMax/(2*npts+1); } return 0; - } + } /** - * This method computes the leading edge time - * and time over threshold - * over a constant fraction of the peak value - * - * @return an int for status - * - */ - public int waveformCFAprocessing(){ + * This method computes the leading edge time + * and time over threshold + * over a constant fraction of the peak value + * + * @return an int for status + * + */ + public int waveformCFAprocessing(AHDCPulse p){ - //Initialization - this.pulse.leadingEdgeTime = -9999; - this.pulse.trailingEdgeTime = -9999; - this.pulse.timeOverThreshold = -9999; + //Initialization + p.leadingEdgeTime = -9999; + p.trailingEdgeTime = -9999; + p.timeOverThreshold = -9999; //Set the CFA threshold - float threshold = this.amplitudeFractionCFA*this.pulse.adcMax; + float threshold = this.amplitudeFractionCFA*p.adcMax; //Crossing the threshold before the peak int binRise = -99; - for (int bin = 0; bin < binMax - 1; bin++){ - if (this.samples[bin] < threshold && this.samples[bin+1] >= threshold) + for (int bin = 0; bin < p.binMax - 1; bin++){ + if (p.samples[bin] < threshold && p.samples[bin+1] >= threshold) binRise = bin; //Here we keep only the last time the signal crosses the threshold } //If the waveform does not cross the ths before the peak //it was early or remant from a previous event and we cannot define a leading time //we set the time at zero - if(binRise==-99) { - this.assignValidType(5); - this.pulse.leadingEdgeTime = (0 + this.time_ZS)*this.samplingTime; + if(binRise==-99) { + this.assignValidType(5, p); + p.leadingEdgeTime = (0 + p.time_ZS)*this.samplingTime; } else { float slopeRise = 0; //linear interpolation //threshold = leadingtime*(ADC1-ADC0)+ADC0 - if (binRise + 1 <= numberOfBins-1) - slopeRise = this.samples[binRise+1] - this.samples[binRise]; - float fittedBinRise = (slopeRise == 0) ? binRise : binRise + (threshold - this.samples[binRise])/slopeRise; - this.pulse.leadingEdgeTime = (fittedBinRise + this.time_ZS)*this.samplingTime; + if (binRise + 1 <= p.numberOfBins-1) + slopeRise = p.samples[binRise+1] - p.samples[binRise]; + float fittedBinRise = (slopeRise == 0) ? binRise : binRise + (threshold - p.samples[binRise])/slopeRise; + p.leadingEdgeTime = (fittedBinRise + p.time_ZS)*this.samplingTime; } //Crossing the threshold back down after the peak int binFall = -99; - for (int bin = binMax; bin < numberOfBins-1; bin++){ - if (this.samples[bin] > threshold - && this.samples[bin+1] <= threshold - && this.samples[bin]>0 - && this.samples[bin+1]>0){ + for (int bin = p.binMax; bin < p.numberOfBins-1; bin++){ + if (p.samples[bin] > threshold + && p.samples[bin+1] <= threshold + && p.samples[bin]>0 + && p.samples[bin+1]>0){ binFall = bin+1; break; //We keep only the first time the signal crosses back the threshold } @@ -261,87 +258,85 @@ public int waveformCFAprocessing(){ //it was late and falling edge cannot be defined //the time correspond to the end of the signal if(binFall==-99) { - this.assignValidType(2); - this.pulse.trailingEdgeTime = (this.effectiveNumberOfBins -1 + this.time_ZS)*this.samplingTime; + this.assignValidType(2, p); + p.trailingEdgeTime = (p.effectiveNumberOfBins -1 + p.time_ZS)*this.samplingTime; } else { - float slopeFall = 0; - if (binFall - 1 >= 0) - slopeFall = this.samples[binFall] - this.samples[binFall-1]; - float fittedBinFall = (slopeFall == 0) ? binFall : binFall-1 + (threshold - this.samples[binFall-1])/slopeFall; - this.pulse.trailingEdgeTime = (fittedBinFall + this.time_ZS)*this.samplingTime; + float slopeFall = 0; + if (binFall - 1 >= 0) + slopeFall = p.samples[binFall] - p.samples[binFall-1]; + float fittedBinFall = (slopeFall == 0) ? binFall : binFall-1 + (threshold - p.samples[binFall-1])/slopeFall; + p.trailingEdgeTime = (fittedBinFall + p.time_ZS)*this.samplingTime; } - - this.pulse.timeOverThreshold = this.pulse.trailingEdgeTime - this.pulse.leadingEdgeTime; + + p.timeOverThreshold = p.trailingEdgeTime - p.leadingEdgeTime; //if ((this.pulse.timeOverThreshold < 300) || (this.pulse.timeOverThreshold > 750)) this.assignValidType(4); // bad tot - if (this.pulse.timeOverThreshold < 300) this.assignValidType(4); // tot too small + if (p.timeOverThreshold < 300) this.assignValidType(4, p); // tot too small return 0; - } + } - /** - * This methods extracts a time using the Constant Fraction Discriminator (CFD) algorithm - * as described in https://commons.wikimedia.org/wiki/File:CFD_Diagram1.jpg for example - * - * @return an int for status - * - */ - public int computeTimeUsingConstantFractionDiscriminator(){ - //Dummy values for hits for which the leading edge is not defined - this.pulse.constantFractionTime = -99; - if(this.pulse.wftype>=4) return 0; - - //int numberOfBins = this.samples.length; - float[] signal = new float[numberOfBins]; - - // signal generation - for (int bin = 0; bin < numberOfBins; bin++){ - signal[bin] = (1 - this.fractionCFD)*this.samples[bin]; // we fill it with a fraction of the original signal - if (bin < numberOfBins - this.binDelayCFD) - signal[bin] += -1*this.fractionCFD*this.samples[bin + this.binDelayCFD]; // we advance and invert a complementary fraction of the original signal and superimpose it to the previous signal - } - // determine the two humps - int binHumpSup = 0; - int binHumpInf = 0; - for (int bin = 0; bin < numberOfBins; bin++){ - if (signal[bin] > signal[binHumpSup]) - binHumpSup = bin; - } - for (int bin = 0; bin < binHumpSup; bin++){ // this loop has been added to be sure : binHumpInf < binHumpSup - if (signal[bin] < signal[binHumpInf]) - binHumpInf = bin; - } - // research for zero - int binZero = 0; - for (int bin = binHumpInf; bin <= binHumpSup; bin++){ - if (signal[bin] < 0) - binZero = bin; // last pass below zero - } // at this stage : binZero < constantFractionTime/samplingTime <= binZero + 1 // constantFractionTime is determined by assuming a linear fit between binZero and binZero + 1 - float slopeCFD = 0; - if (binZero + 1 <= numberOfBins) - slopeCFD = signal[binZero+1] - signal[binZero]; - float fittedBinZero = (slopeCFD == 0) ? binZero : binZero + (0 - signal[binZero])/slopeCFD; - this.pulse.constantFractionTime = (fittedBinZero + this.time_ZS)*this.samplingTime; - - return 0; + /** + * This methods extracts a time using the Constant Fraction Discriminator (CFD) algorithm + * as described in https://commons.wikimedia.org/wiki/File:CFD_Diagram1.jpg for example + * + * @return an int for status + * + */ + public int computeTimeUsingConstantFractionDiscriminator(AHDCPulse p){ + //Dummy values for hits for which the leading edge is not defined + p.constantFractionTime = -99; + if(p.wftype>=4) return 0; + + //int numberOfBins = this.samples.length; + float[] signal = new float[p.numberOfBins]; + + // signal generation + for (int bin = 0; bin < p.numberOfBins; bin++){ + signal[bin] = (1 - this.fractionCFD)*p.samples[bin]; // we fill it with a fraction of the original signal + if (bin < p.numberOfBins - this.binDelayCFD) + signal[bin] += -1*this.fractionCFD*p.samples[bin + this.binDelayCFD]; // we advance and invert a complementary fraction of the original signal and superimpose it to the previous signal } - - /** - * Apply fine timestamp correction - * - * adapted from decode/MVTFitter.java - * - * @param timestamp for fine time correction - * @param fineTimeStampResolution correspond to the dream clock (usually equals to 8; but to be checked!) - */ - - private void fineTimeStampCorrection(long timestamp, float fineTimeStampResolution) { - long fineTimeStamp = timestamp & 0x00000007; // keep and convert last 3 bits of binary timestamp - this.pulse.leadingEdgeTime += (double) ((fineTimeStamp+0.5) * fineTimeStampResolution); //fineTimeStampCorrection, only on the leadingEdgeTime for the moment (we don't use timeMax or constantFractionTime in the reconstruction yet) - } - - + // determine the two humps + int binHumpSup = 0; + int binHumpInf = 0; + for (int bin = 0; bin < p.numberOfBins; bin++){ + if (signal[bin] > signal[binHumpSup]) + binHumpSup = bin; + } + for (int bin = 0; bin < binHumpSup; bin++){ // this loop has been added to be sure : binHumpInf < binHumpSup + if (signal[bin] < signal[binHumpInf]) + binHumpInf = bin; + } + // research for zero + int binZero = 0; + for (int bin = binHumpInf; bin <= binHumpSup; bin++){ + if (signal[bin] < 0) + binZero = bin; // last pass below zero + } // at this stage : binZero < constantFractionTime/samplingTime <= binZero + 1 // constantFractionTime is determined by assuming a linear fit between binZero and binZero + 1 + float slopeCFD = 0; + if (binZero + 1 <= p.numberOfBins) + slopeCFD = signal[binZero+1] - signal[binZero]; + float fittedBinZero = (slopeCFD == 0) ? binZero : binZero + (0 - signal[binZero])/slopeCFD; + p.constantFractionTime = (fittedBinZero + p.time_ZS)*this.samplingTime; + + return 0; + } + + /** + * Apply fine timestamp correction + * + * adapted from decode/MVTFitter.java + * + * @param timestamp for fine time correction + * @param fineTimeStampResolution correspond to the dream clock (usually equals to 8; but to be checked!) + */ + + private void fineTimeStampCorrection(long timestamp, float fineTimeStampResolution, AHDCPulse p) { + long fineTimeStamp = timestamp & 0x00000007; // keep and convert last 3 bits of binary timestamp + p.leadingEdgeTime += (double) ((fineTimeStamp+0.5) * fineTimeStampResolution); //fineTimeStampCorrection, only on the leadingEdgeTime for the moment (we don't use timeMax or constantFractionTime in the reconstruction yet) + } /** * This method extracts relevant information from the waveform @@ -356,39 +351,39 @@ private void fineTimeStampCorrection(long timestamp, float fineTimeStampResoluti @Override public List extract(NamedEntry pars, int id, long timestamp, long time_ZS, short... samples){ - //Initialize everything each time a new wf is read + //Initialize everything each time a new wf is read //and new adc information is extracted. //This is messy but avoids feeding the same arguments to all methods - this.pulse = new Pulse(); - this.pulse.id = id; - this.pulse.timestamp = timestamp; - this.time_ZS = time_ZS; - this.samples = samples; - this.binMax = -99; + AHDCPulse pulse = new AHDCPulse(); + pulse.id = id; + pulse.timestamp = timestamp; + pulse.time_ZS = time_ZS; + pulse.samples = samples; + pulse.binMax = -99; - List output = new ArrayList<>(); - // Pre Process - this.preProcess(); + this.preProcess(pulse); + //Baseline computation - this.baselineComputation(); - + this.baselineComputation(pulse); + //Get the ADC information from the pulse (peak ADC and time, integral) - this.waveformADCProcessing(); + this.waveformADCProcessing(pulse); //Get the time overr threshold - this.waveformCFAprocessing(); + this.waveformCFAprocessing(pulse); //Get the CFD time - this.computeTimeUsingConstantFractionDiscriminator(); - - // Fine timestamp correction on leadingEdgeTime - this.fineTimeStampCorrection(timestamp, dream_clock); + this.computeTimeUsingConstantFractionDiscriminator(pulse); - output.add(this.pulse); + // Fine timestamp correction on leadingEdgeTime + this.fineTimeStampCorrection(timestamp, dream_clock, pulse); + + List output = new ArrayList<>(); + output.add(pulse); return output; } - + @Override public void update(int n, IndexedTable it, DataEvent event, String wfBankName, String adcBankName) { DataBank wf = event.getBank(wfBankName); @@ -412,26 +407,26 @@ public void update(int n, IndexedTable it, DataEvent event, String wfBankName, S } } } - + @Override protected void update(int n, IndexedTable it, Bank wfBank, Bank adcBank) { - if (wfBank.getRows() > 0) { - List pulses = getPulses(n, it, wfBank); - adcBank.reset(); - adcBank.setRows(pulses!=null ? pulses.size() : 0); - if (pulses!=null && !pulses.isEmpty()) { - for (int i=0; i 0) { + List pulses = getPulses(n, it, wfBank); + adcBank.reset(); + adcBank.setRows(pulses!=null ? pulses.size() : 0); + if (pulses!=null && !pulses.isEmpty()) { + for (int i=0; i pair = trackAIResults.get(i); - bank.setInt("track_id", i, pair.getKey()); + bank.setInt("trackid", i, pair.getKey()); bank.setInt("matched_atof_hit_id", i, pair.getValue()); } event.appendBank(bank); diff --git a/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java index 4f8dda448f..b7cbb1b01c 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java @@ -6,7 +6,6 @@ import org.jlab.io.base.DataEvent; import org.jlab.io.hipo.HipoDataSource; import org.jlab.io.hipo.HipoDataSync; -import org.jlab.jnp.hipo4.data.SchemaFactory; import org.jlab.rec.ahdc.AI.*; import org.jlab.rec.ahdc.Banks.RecoBankWriter; import org.jlab.rec.ahdc.Cluster.Cluster; @@ -15,7 +14,6 @@ import org.jlab.rec.ahdc.HelixFit.HelixFitJava; import org.jlab.rec.ahdc.Hit.Hit; import org.jlab.rec.ahdc.Hit.HitReader; -import org.jlab.rec.ahdc.Hit.TrueHit; import org.jlab.rec.ahdc.HoughTransform.HoughTransform; import org.jlab.rec.ahdc.KalmanFilter.KalmanFilter; import org.jlab.rec.ahdc.KalmanFilter.MaterialMap; @@ -24,7 +22,6 @@ import org.jlab.rec.ahdc.Track.Track; import org.jlab.rec.ahdc.ModeTrackFinding; import java.io.File; -import java.lang.reflect.Array; import java.util.*; import java.util.logging.Logger; diff --git a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java index c4483be4eb..6a54552d89 100644 --- a/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java +++ b/reconstruction/alert/src/main/java/org/jlab/service/alert/ALERTEngine.java @@ -3,7 +3,6 @@ import java.util.ArrayList; import java.util.concurrent.atomic.AtomicInteger; import java.io.File; -import java.util.*; import org.jlab.detector.calib.utils.DatabaseConstantProvider; import org.jlab.geom.base.Detector; @@ -16,10 +15,6 @@ import org.jlab.clas.reco.ReconstructionEngine; import org.jlab.clas.swimtools.Swim; -import org.jlab.rec.ahdc.AI.InterCluster; -import org.jlab.rec.ahdc.AI.PreClustering; -import org.jlab.rec.ahdc.Hit.Hit; -import org.jlab.rec.ahdc.PreCluster.PreCluster; import org.jlab.rec.alert.TrackMatchingAI.ModelTrackMatching; import org.jlab.rec.alert.banks.RecoBankWriter; @@ -153,12 +148,12 @@ public boolean processDataEvent(DataEvent event) { ArrayList> matched_ATOF_hit_id = new ArrayList<>(); for (int i = 0; i < bank_AHDCtracks.rows(); i++) { - int track_id = bank_AHDCtracks.getInt("track_id", i); + int track_id = bank_AHDCtracks.getInt("trackid", i); // Get all interclusters for this track ArrayList> interClusters = new ArrayList<>(); for (int j = 0; j < bank_AHDCInterclusters.rows(); j++) { - int intercluster_track_id = bank_AHDCInterclusters.getInt("track_id", j); + int intercluster_track_id = bank_AHDCInterclusters.getInt("trackid", j); if (intercluster_track_id == track_id) { float x = bank_AHDCInterclusters.getFloat("x", j); float y = bank_AHDCInterclusters.getFloat("y", j);