diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitter.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitter.java index 2ecff1db9e..90f47e65ef 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitter.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitter.java @@ -718,8 +718,6 @@ private boolean filter(int k, boolean forward, double annealingFactor) { } } - // Since no vertex inforamtion, the starting point for path length is the final point at the last layer. - // After vertex information is obtained, transition for the starting point from the final point to vertex will be taken. private boolean filter(int k, boolean forward) { StateVec sVec = sv.transported(forward).get(k); org.jlab.clas.tracking.kalmanfilter.AMeasVecs.MeasVec mVec = mv.measurements.get(k); @@ -895,7 +893,7 @@ public void calcFinalChisq(int sector, boolean nofilter) { sv.transport(sector, k, 0, sVec, mv, this.getSwimmer(), forward); StateVec svc = sv.transported(forward).get(0); - path += (forward ? 1 : -1) * svc.deltaPath; + path += svc.deltaPath; svc.setPathLength(path); double V0 = mv.measurements.get(0).surface.unc[0]; @@ -941,7 +939,7 @@ public void calcFinalChisq(int sector, boolean nofilter) { double h = mv.hDoca(line, mv.measurements.get(k1 + 1).surface.wireLine[0]); svc = sv.transported(forward).get(k1 + 1); - path += (forward ? 1 : -1) * svc.deltaPath; + path += svc.deltaPath; svc.setPathLength(path); svc.setProjector(mv.measurements.get(k1 + 1).surface.wireLine[0].origin().x()); svc.setProjectorDoca(h); @@ -994,7 +992,7 @@ private void calcFinalChisqDAF(int sector, boolean nofilter) { sv.transport(sector, k, 0, sVec, mv, this.getSwimmer(), forward); StateVec svc = sv.transported(forward).get(0); - path += (forward ? 1 : -1) * svc.deltaPath; + path += svc.deltaPath; svc.setPathLength(path); Point3D point = new Point3D(svc.x, svc.y, mv.measurements.get(0).surface.measPoint.z()); @@ -1068,7 +1066,7 @@ private void calcFinalChisqDAF(int sector, boolean nofilter) { } svc = sv.transported(forward).get(k1 + 1); - path += (forward ? 1 : -1) * svc.deltaPath; + path += svc.deltaPath; svc.setPathLength(path); point = new Point3D(sv.transported(forward).get(k1 + 1).x, sv.transported(forward).get(k1 + 1).y, mv.measurements.get(k1 + 1).surface.measPoint.z()); @@ -1140,10 +1138,6 @@ public Matrix propagateToVtx(int sector, double Zf) { return sv.transport(sector, finalStateVec.k, Zf, finalStateVec, mv, this.getSwimmer()); } - public double getDeltaPathToVtx(int sector, double Zf) { - return sv.getDeltaPath(sector, finalStateVec.k, Zf, finalStateVec, mv, this.getSwimmer()); - } - //Todo: apply the common funciton to replace current function above @Override public void runFitter(AStateVecs sv, AMeasVecs mv) { diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitterStraight.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitterStraight.java index 60943ac4c8..6435c995e7 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitterStraight.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/KFitterStraight.java @@ -411,8 +411,6 @@ private void calcFinalChisq(int sector) { calcFinalChisq(sector, false); } - // Since no vertex inforamtion, the starting point for path length is the final point at the last layer. - // After vertex information is obtained, transition for the starting point from the final point to vertex will be taken. private void calcFinalChisq(int sector, boolean nofilter) { int k = svzLength - 1; this.chi2 = 0; @@ -437,7 +435,7 @@ private void calcFinalChisq(int sector, boolean nofilter) { sv.transport(sector, k, 0, sVec, mv, this.getSwimmer(), forward); StateVec svc = sv.transported(forward).get(0); - path += (forward ? 1 : -1) * svc.deltaPath; + path += svc.deltaPath; svc.setPathLength(path); double V0 = mv.measurements.get(0).surface.unc[0]; @@ -484,7 +482,7 @@ private void calcFinalChisq(int sector, boolean nofilter) { double h = mv.hDoca(line, mv.measurements.get(k1 + 1).surface.wireLine[0]); svc = sv.transported(forward).get(k1+1); - path += (forward ? 1 : -1) * svc.deltaPath; + path += svc.deltaPath; svc.setPathLength(path); svc.setProjector(mv.measurements.get(k1 + 1).surface.wireLine[0].origin().x()); svc.setProjectorDoca(h); @@ -515,10 +513,6 @@ public Matrix propagateToVtx(int sector, double Zf) { return sv.transport(sector, finalStateVec.k, Zf, finalStateVec, mv, this.getSwimmer()); } - public double getDeltaPathToVtx(int sector, double Zf) { - return sv.getDeltaPath(sector, finalStateVec.k, Zf, finalStateVec, mv, this.getSwimmer()); - } - @Override public void runFitter(AStateVecs sv, AMeasVecs mv) { throw new UnsupportedOperationException("Not supported yet."); diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/StateVecs.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/StateVecs.java index 1fe0751241..738b078e58 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/StateVecs.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/zReference/StateVecs.java @@ -152,104 +152,6 @@ public Matrix transport(int sector, int i, double Zf, StateVec iVec, AMeasVecs m } - /** - * - * @param sector - * @param i initial state vector index - * @param Zf - * @param iVec state vector at the initial index - * @param mv measurements - */ - public double getDeltaPath(int sector, int i, double Zf, StateVec iVec, AMeasVecs mv, Swim swimmer) { // s = signed step-size - - double stepSize = 1.0; - StateVec fVec = new StateVec(0); - fVec.x = iVec.x; - fVec.y = iVec.y; - fVec.z = iVec.z; - fVec.tx = iVec.tx; - fVec.ty = iVec.ty; - fVec.Q = iVec.Q; - fVec.B = iVec.B; - Matrix5x5.copy(iVec.CM, fVec.CM); - - double s = 0; - double zInit = mv.measurements.get(i).surface.measPoint.z(); - double BatMeas = iVec.B; - - double z = zInit; - - while (Math.signum(Zf - zInit) * z < Math.signum(Zf - zInit) * Zf) { - z = fVec.z; - if (z == Zf) { - break; - } - - double x = fVec.x; - double y = fVec.y; - double tx = fVec.tx; - double ty = fVec.ty; - double Q = fVec.Q; - double dPath = fVec.deltaPath; - Matrix cMat = new Matrix(); - Matrix5x5.copy(fVec.CM, cMat); - s = Math.signum(Zf - zInit) * stepSize; - - // LOGGER.log(Level.FINEST, " from "+(float)Z[i]+" to "+(float)Z[f]+" at "+(float)z+" By is "+bf[1]+" B is "+Math.sqrt(bf[0]*bf[0]+bf[1]*bf[1]+bf[2]*bf[2])/Bmax+" stepSize is "+s); - if (Math.signum(Zf - zInit) * (z + s) > Math.signum(Zf - zInit) * Zf) { - s = Math.signum(Zf - zInit) * Math.abs(Zf - z); - } - - //rk.RK4transport(sector, Q, x, y, z, tx, ty, s, swimmer, cMat, fVec, dPath); - rk.RK4transport(sector, s, swimmer, cMat, fVec); - - // Q process noise matrix estimate - double p = Math.abs(1. / iVec.Q); - - double X0 = this.getX0(mv.measurements.get(i).surface, z, Z); - double t_ov_X0 = Math.abs(s) / X0;//path length in radiation length units = t/X0 [true path length/ X0] ; Ar radiation length = 14 cm - - double beta = this.beta; - if (beta > 1.0 || beta <= 0) { - beta = 1.0; - } - - double sctRMS = 0; - - ////// Todo: Modify multi-scattering or remove it; After update, some parameters, like iteration termintion chonditions, may need to be updated. - // Speed of light should be 1 - // From one measurement site to another, F and Q should be calculated separaetely with multiple steps; and then C' = FTCF + Q - if (Math.abs(s) > 0) { - sctRMS = ((0.0136) / (beta * PhysicsConstants.speedOfLight() * p)) * Math.sqrt(t_ov_X0) - * (1 + 0.038 * Math.log(t_ov_X0)); - } - - double cov_txtx = (1 + tx * tx) * (1 + tx * tx + ty * ty) * sctRMS * sctRMS; - double cov_tyty = (1 + ty * ty) * (1 + tx * tx + ty * ty) * sctRMS * sctRMS; - double cov_txty = tx * ty * (1 + tx * tx + ty * ty) * sctRMS * sctRMS; - - fMS.set( - 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, - 0, 0, cov_txtx, cov_txty, 0, - 0, 0, cov_txty, cov_tyty, 0, - 0, 0, 0, 0, 0 - ); - - Matrix5x5.copy(fVec.CM, copyMatrix); - Matrix5x5.add(copyMatrix, fMS, fVec.CM); - - if (Math.abs(fVec.B - BatMeas) < 0.0001) { - stepSize *= 2; - } - - BatMeas = fVec.B; - } - - return fVec.deltaPath; - - } - /** * * @param sector diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/Constants.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/Constants.java index 82473ccf2b..4cb45072ea 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/Constants.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/Constants.java @@ -96,7 +96,7 @@ public static Constants getInstance() { public int SECTORSELECT = 0; public int NSUPERLAYERTRACKING = 5; private boolean USETSTART = true; - //private boolean USETIMETBETA = false; + private boolean USETIMETBETA = false; private boolean CHECKBETA = false; private int T2D = 1; // 1=polynomial, 0=exponential private boolean USEDOUBLETS = false; @@ -107,11 +107,11 @@ public static Constants getInstance() { public static final String TT = "/daq/tt/dc"; public static final String DOCARES = "/calibration/dc/signal_generation/doca_resolution"; public static final String TIME2DIST = "/calibration/dc/time_to_distance/time2dist"; - public static final String T2DPRESSURE = "/calibration/dc/v2/t2d_pressure"; + public static final String T2DPRESSURE = "/calibration/dc/time_to_distance/t2d_pressure"; public static final String PRESSURE = "/hall/weather/pressure"; - public static final String T2DPRESSUREREF = "/calibration/dc/v2/ref_pressure"; - public static final String T0CORRECTION = "/calibration/dc/v2/t0"; - public static final String TDCTCUTS = "/calibration/dc/v2/tdc_cuts"; + public static final String T2DPRESSUREREF = "/calibration/dc/time_to_distance/ref_pressure"; + public static final String T0CORRECTION = "/calibration/dc/time_corrections/T0Corrections"; + public static final String TDCTCUTS = "/calibration/dc/time_corrections/tdctimingcuts"; public static final String WIRESTAT = "/calibration/dc/tracking/wire_status"; public static final String TIMEJITTER = "/calibration/dc/time_jitter"; public static final String BEAMPOS = "/geometry/beam/position"; @@ -298,13 +298,13 @@ public void setUSETSTART(boolean usetstart) { USETSTART = usetstart; } -// public boolean useUSETIMETBETA() { -// return USETIMETBETA; -// } -// -// public void setUSETIMETBETA(boolean usetimebeta) { -// USETIMETBETA = usetimebeta; -// } + public boolean useUSETIMETBETA() { + return USETIMETBETA; + } + + public void setUSETIMETBETA(boolean usetimebeta) { + USETIMETBETA = usetimebeta; + } public double getWIREDIST() { return WIREDIST; @@ -424,7 +424,7 @@ public void printConfig(String engine) { LOGGER.log(Level.INFO, "["+engine+"] run with wire ministagger = " + MINISTAGGERSTATUS.getName()); LOGGER.log(Level.INFO, "["+engine+"] run with wire feedthroughs = " + FEEDTHROUGHSSTATUS.getName()); LOGGER.log(Level.INFO, "["+engine+"] run with wire distortions = " + ENDPLATESBOWING); - //LOGGER.log(Level.INFO, "["+engine+"] run with with time Beta correction (is false for doca Beta correction) = " + USETIMETBETA); + LOGGER.log(Level.INFO, "["+engine+"] run with with time Beta correction (is false for doca Beta correction) = " + USETIMETBETA); LOGGER.log(Level.INFO, "["+engine+"] run with with Beta cut = " + CHECKBETA); LOGGER.log(Level.INFO, "["+engine+"] run with time to distance function set to exponential/polynomial (0/1) = " + T2D); LOGGER.log(Level.INFO, "["+engine+"] run with with hit doublets recovery = " + USEDOUBLETS); diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/HitReader.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/HitReader.java index e9e88f56b8..f33e6e5e3d 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/HitReader.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/HitReader.java @@ -13,6 +13,7 @@ import org.jlab.rec.dc.hit.Hit; import org.jlab.rec.dc.timetodistance.TimeToDistanceEstimator; +import org.jlab.clas.swimtools.Swimmer; import org.jlab.detector.banks.RawBank.OrderType; import org.jlab.detector.banks.RawDataBank; import org.jlab.detector.calib.utils.ConstantsManager; @@ -54,6 +55,8 @@ public class HitReader { private List _HBHits; //hit-based tracking hit information private List _TBHits; //time-based tracking hit information + private final double timeBuf = 25.0; + private static final Logger LOGGER = Logger.getLogger(HitReader.class.getName()); public HitReader(Banks names, DCGeant4Factory detector) { @@ -204,15 +207,49 @@ public void fetch_DCHits(DataEvent event) { double timeCutMin = 0; double timeCutMax = 0; + double timeCutLC = 0; int region = ((superlayer + 1) / 2); - long hash = IndexedTable.DEFAULT_GENERATOR.hashCode(0,region,wire); - timeCutMin = tdccuts.getIntValueByHash("min", hash); - timeCutMax = tdccuts.getIntValueByHash("max", hash); + switch (region) { + case 1: + timeCutMin = tdccuts.getIntValue("MinEdge", 0, region, 0); + timeCutMax = tdccuts.getIntValue("MaxEdge", 0, region, 0); + break; + case 2: + if (wire <= 56) { + timeCutLC = tdccuts.getIntValue("LinearCoeff", 0, region, 1); + timeCutMin = tdccuts.getIntValue("MinEdge", 0, region, 1); + timeCutMax = tdccuts.getIntValue("MaxEdge", 0, region, 1); + } + if (wire > 56) { + timeCutLC = tdccuts.getIntValue("LinearCoeff", 0, region, 56); + timeCutMin = tdccuts.getIntValue("MinEdge", 0, region, 56); + timeCutMax = tdccuts.getIntValue("MaxEdge", 0, region, 56); + } + break; + case 3: + timeCutMin = tdccuts.getIntValue("MinEdge", 0, region, 0); + timeCutMax = tdccuts.getIntValue("MaxEdge", 0, region, 0)+timeBuf; + break; + } boolean passTimingCut = false; - if (tdc > timeCutMin && tdc < timeCutMax) + if (region == 1 && tdc > timeCutMin && tdc < timeCutMax) + passTimingCut = true; + if (region == 2) { + double Bscale = Swimmer.getTorScale() * Swimmer.getTorScale(); + if (wire >= 56) { + if (tdc > timeCutMin && + tdc < timeCutMax + timeCutLC * (double) (112 - wire / 56) * Bscale) + passTimingCut = true; + } else { + if (tdc > timeCutMin && + tdc < timeCutMax + timeCutLC * (double) (56 - wire / 56) * Bscale) + passTimingCut = true; + } + } + if (region == 3 && tdc > timeCutMin && tdc < timeCutMax) passTimingCut = true; if (passTimingCut) { // cut on spurious hits @@ -657,8 +694,8 @@ private double[] getT0(int sector, int superlayer, int cable = this.getCableID1to6(layer, wire); int slot = this.getSlotID1to7(wire); - double t0 = t0Table.getDoubleValue("t0correction", sector, superlayer, slot, cable); - double t0E = t0Table.getDoubleValue("t0error", sector, superlayer, slot, cable); + double t0 = t0Table.getDoubleValue("T0Correction", sector, superlayer, slot, cable); + double t0E = t0Table.getDoubleValue("T0Error", sector, superlayer, slot, cable); T0Corr[0] = t0; T0Corr[1] = t0E; diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/hit/FittedHit.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/hit/FittedHit.java index c9dec90771..edc6867e4e 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/hit/FittedHit.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/hit/FittedHit.java @@ -383,28 +383,27 @@ public void set_TimeToDistance(DataEvent event, double trkAngle, double B, Index double deltatime_beta = 0; double deltadoca_beta = 0; -// if(Constants.getInstance().useUSETIMETBETA()==true) { -// deltatime_beta = calcDeltaTimeBetaTFCN(this.get_Time(), tab, beta); -// } -// -// if(event.hasBank("MC::Particle")==false) { -// distance = tde.interpolateOnGrid(B, Math.toDegrees(ralpha), -// this.getCorrectedTime(this.get_Time(), deltatime_beta), -// secIdx, slIdx); -// } else { -// distance = tde.interpolateOnGrid(B, Math.toDegrees(ralpha), -// this.getCorrectedTime(this.get_Time(), 0), -// secIdx, slIdx) ; -// } -// //get deltadoca -// if(Constants.getInstance().useUSETIMETBETA()==false) { -// deltadoca_beta = calcDeltaDocaBeta(distance, tab, beta); -// } -// -// distance -=deltadoca_beta; - //this.set_DeltaTimeBeta(deltatime_beta); - //this.set_DeltaDocaBeta(deltadoca_beta); - distance = tde.interpolateOnGrid(B, Math.toDegrees(ralpha), beta, this.get_Time(), secIdx, slIdx); + if(Constants.getInstance().useUSETIMETBETA()==true) { + deltatime_beta = calcDeltaTimeBetaTFCN(this.get_Time(), tab, beta); + } + + if(event.hasBank("MC::Particle")==false) { + distance = tde.interpolateOnGrid(B, Math.toDegrees(ralpha), + this.getCorrectedTime(this.get_Time(), deltatime_beta), + secIdx, slIdx); + } else { + distance = tde.interpolateOnGrid(B, Math.toDegrees(ralpha), + this.getCorrectedTime(this.get_Time(), 0), + secIdx, slIdx) ; + } + //get deltadoca + if(Constants.getInstance().useUSETIMETBETA()==false) { + deltadoca_beta = calcDeltaDocaBeta(distance, tab, beta); + } + + distance -=deltadoca_beta; + this.set_DeltaTimeBeta(deltatime_beta); + this.set_DeltaDocaBeta(deltadoca_beta); } diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/timetodistance/TableLoader.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/timetodistance/TableLoader.java index d21bd960c0..806033b5f6 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/timetodistance/TableLoader.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/timetodistance/TableLoader.java @@ -1,19 +1,7 @@ package org.jlab.rec.dc.timetodistance; -import java.util.ArrayList; -import java.util.HashMap; -import java.util.List; -import java.util.Map; -import org.jlab.detector.base.DetectorType; -import org.jlab.detector.base.GeometryFactory; -import org.jlab.geom.base.ConstantProvider; -import org.jlab.groot.data.H1F; -import org.jlab.groot.data.H2F; -import org.jlab.groot.math.Func1D; -import org.jlab.groot.ui.TCanvas; import org.jlab.rec.dc.Constants; import static org.jlab.rec.dc.timetodistance.T2DFunctions.polyFcnMac; -import org.jlab.service.dc.DCEngine; import org.jlab.utils.groups.IndexedTable; @@ -25,7 +13,6 @@ public TableLoader() { private static boolean T2DLOADED = false; public static final int NBINST=2000; - public static final double[] betaValues = new double[] {0.6, 0.7, 0.8, 0.9, 1.0}; public static final double[] BfieldValues = new double[]{0.0000, 1.0000, 1.4142, 1.7321, 2.0000, 2.2361, 2.4495, 2.6458}; public static int minBinIdxB = 0; public static int maxBinIdxB = BfieldValues.length-1; @@ -35,187 +22,10 @@ public TableLoader() { private static final double[][] AlphaBounds = new double[6][2]; public static int minBinIdxT = 0; public static final int[][][][] maxBinIdxT = new int[6][6][8][6]; - public static double[][][][][][] DISTFROMTIME = new double[6][6][maxBinIdxB+1][maxBinIdxAlpha+1][betaValues.length][NBINST]; // sector slyr alpha Bfield time bins [s][r][ibfield][icosalpha][tbin] - public static int timeBinWidth = 2; //ns + public static double[][][][][] DISTFROMTIME = new double[6][6][maxBinIdxB+1][maxBinIdxAlpha+1][NBINST]; // sector slyr alpha Bfield time bins [s][r][ibfield][icosalpha][tbin] public static int maxTBin = -1; - public static void main(String[] args) { - DCEngine dce = new DCEngine("test"); - dce.setVariation("default"); - dce.LoadTables(); - ConstantProvider provider = GeometryFactory.getConstants(DetectorType.DC, 11, "default"); - for(int l=0; l<6; l++) { - Constants.getInstance().wpdist[l] = provider.getDouble("/geometry/dc/superlayer/wpdist", l); - } - int run = 18331; - TableLoader.Fill(dce.getConstantsManager().getConstants(run, Constants.T2DPRESSURE), - dce.getConstantsManager().getConstants(run, Constants.T2DPRESSUREREF), - dce.getConstantsManager().getConstants(run, Constants.PRESSURE)); - TableLoader.Fill(dce.getConstantsManager().getConstants(run, Constants.T2DPRESSURE), - dce.getConstantsManager().getConstants(run, Constants.T2DPRESSUREREF), - dce.getConstantsManager().getConstants(run, Constants.PRESSURE)); - test(); - } - - public static void test(){ - - plotFcns(); - - TimeToDistanceEstimator tde = new TimeToDistanceEstimator(); - tde.t2DPrecisionImprov=true; - double[] beta = new double[]{0.65, 0.66, 0.67, 0.68, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0}; - double[] Bfield = new double[]{0.0, 0.5, 1.0, 1.2, 1.3, 1.4, 1.5, 2.0, 2.1, 2.15, 2.2, 2.25, 2.26, 2.27, 2.28, 2.29}; - double[] htmax = new double[]{200,210,1050,1200,720,800}; - double[] hdmax = new double[] {0.77, 0.8, 1.1, 1.2, 1.8, 1.84}; - - List hts = new ArrayList<>(); //histo to check table and interp. from time to idistance (by interpolation) - //to calculated time (from dist.to t) in seconds; as a function of time - List hds = new ArrayList<>(); //histo to check table and interp. from distance to calculated time (from dist.to t) - //to idistance (by interpolation) in microns; as a function of distance - List hdsSL3 = new ArrayList<>(); //histo to check table and interp. from distance to calculated time (from dist.to t) - //for various values of B-field - List hdsSL4 = new ArrayList<>(); //histo to check table and interp. from distance to calculated time (from dist.to t) - //for various values of B-field - List hd2s = new ArrayList<>();// as s function of distance - List ht2d = new ArrayList<>();// time to distance from interpolation - for(int r = 0; r<6; r++ ){ //loop over slys - hts.add(new H1F("ht"+(r+1), "time resolution (ns)", "Counts/0.1 ns", 400, -20.0,20.0)); - hds.add(new H1F("hd"+(r+1), "doca resolution (um)", "Counts/0.1 um", 400, -20.0,20.0)); - //public H2F(String name, String title, int bx, double xmin, double xmax, int by, double ymin, double ymax) - hd2s.add(new H2F("hd2"+(r+1), "doca resolution (um) vs doca (cm)", (int)(hdmax[r]*100), 0, hdmax[r], 400, -20.0,20.0)); - ht2d.add(new H2F("ht2d"+(r+1), "doca (cm) vs time(time (ns)", (int)htmax[r], 0, htmax[r], (int)(hdmax[r]*100), 0, hdmax[r])); - } - for (int b =0; b3) maxBidx=1; - for(int ibfield =0; ibfield fmap = new HashMap<>(); - for (int i = 0; i < 6; i++) { - CanD2T.cd(i); - double[] pars = new double[11]; - pars[0] = TableLoader.v0[0][i]; - pars[1] = TableLoader.vmid[0][i]; - pars[2] = TableLoader.FracDmaxAtMinVel[0][i]; - pars[3] = TableLoader.Tmax[0][i]; - pars[4] = TableLoader.distbeta[0][i]; - pars[5] = TableLoader.delta_bfield_coefficient[0][i]; - pars[6] = TableLoader.b1[0][i]; - pars[7] = TableLoader.b2[0][i]; - pars[8] = TableLoader.b3[0][i]; - pars[9] = TableLoader.b4[0][i]; - pars[10] = 2.*Constants.getInstance().wpdist[i];//fix dmax - - for(int j =0; j0 && (i<2 || i>3) ) continue; - String name = "f"; - name+=i; - name+="."; - name+=j; - name+="."; - name+=k; - fmap.put(name, new FitLine(name, 0, pars[10], 0, i, 3, alpha, BfieldValues[k])); - fmap.get(name).setLineWidth(2); - fmap.get(name).setLineColor(k+1); - fmap.get(name).setRange(0, pars[10]); - fmap.get(name).getAttributes().setTitleX("doca (cm)"); - fmap.get(name).getAttributes().setTitleY("time (ns)"); - CanD2T.draw(fmap.get(name), "same"); - } - } - } - } - private static int getAlphaBin(double Alpha) { int bin = 0; for(int b =0; b<6; b++) { @@ -270,9 +80,6 @@ private static void getConstants(IndexedTable t2dPressure, IndexedTable t2dPress vmid[s][r] = t2dPressure.getDoubleValueByHash("vmid_a0", hash) +t2dPressure.getDoubleValueByHash("vmid_a1", hash)*dp*dpscale +t2dPressure.getDoubleValueByHash("vmid_a2", hash)*dp*dp*dp2scale; - distbeta[s][r] = t2dPressure.getDoubleValueByHash("distbeta_a0", hash) - +t2dPressure.getDoubleValueByHash("distbeta_a1", hash)*dp*dpscale - +t2dPressure.getDoubleValueByHash("distbeta_a2", hash)*dp*dp*dp2scale; if(r>1 && r<4) { delta_bfield_coefficient[s][r] = t2dPressure.getDoubleValueByHash("delta_bfield_a0", hash) +t2dPressure.getDoubleValueByHash("delta_bfield_a1", hash)*dp*dpscale @@ -306,25 +113,16 @@ public static void FillTable() { //double tmax = CCDBConstants.getTMAXSUPERLAYER()[s][r]; for(int ibfield =0; ibfieldmaxTime) - maxTime=timebfield; - //System.out.println("T "+timebfield+" maxT "+maxTime+" x "+x); int tbin = (int) Math.floor(timebfield/2); if(tbin<0 || tbin>NBINST-1) { //System.err.println("Problem with tbin"); @@ -332,29 +130,32 @@ public static void FillTable() { } if(tbin>maxTBin) maxTBin = tbin; - - if(timebfieldmaxBinIdxT[s][r][ibfield][icosalpha]) { + //maxBinIdxT[s][r][ibfield][icosalpha] = NBINST; + //} //LOGGER.log(Level.FINE, "tbin "+tbin+" tmax "+tmax+ "s "+s+" sl "+r ); + if(DISTFROMTIME[s][r][ibfield][icosalpha][tbin]==0) { + // firstbin = bi + // bincount = 0; + DISTFROMTIME[s][r][ibfield][icosalpha][tbin]=x; + } else { + // test for getting center of the bin (to be validated): + //double prevTime = calc_Time(x-stepSize, alpha, bfield, s+1, r+1); + //if(x>DISTFROMTIME[s][r][ibfield][icosalpha][tbin] + // && Math.abs((double)(2.*tbin+1)-timebfield)<=Math.abs((double)(2.*tbin+1)-prevTime)) { + // DISTFROMTIME[s][r][ibfield][icosalpha][tbin]=x; + //} + // bincount++; + DISTFROMTIME[s][r][ibfield][icosalpha][tbin]+=stepSize; } - } - if(DISTFROMTIME[s][r][ibfield][icosalpha][ibeta][tbin]>dmax) { - DISTFROMTIME[s][r][ibfield][icosalpha][ibeta][tbin]=dmax; - idist=nxmax; - } + /* if(timebfield>timebfield_max) { + DISTFROMTIME[s][r][ibfield][icosalpha][tbin]=x-stepSize*0.5; + if(DISTFROMTIME[s][r][ibfield][icosalpha][tbin]>dmax) + DISTFROMTIME[s][r][ibfield][icosalpha][tbin] = dmax; + } */ } } } - } } } @@ -375,7 +176,7 @@ public static synchronized void Fill(IndexedTable t2dPressure, IndexedTable t2dP T2DLOADED = true; } - private static void fillMissingTableBins() { + private static synchronized void fillMissingTableBins() { for(int s = 0; s<6; s++ ){ // loop over sectors @@ -385,21 +186,17 @@ private static void fillMissingTableBins() { for(int icosalpha =0; icosalphaTableLoader.NBINST-1) { - tHi=TableLoader.NBINST-1; - } - } - double timeLo=2*tLo+1; - double timeHi=2*tHi+1; - //get the beta bins - int binBeta = this.getBetaIdx(beta); - int betaLo = 0; - int betaHi = 0; - double betaCent = TableLoader.betaValues[binBeta]; - if(betaTableLoader.betaValues.length-1) { - betaHi=TableLoader.betaValues.length-1; - } - } - double betaValueLo = TableLoader.betaValues[betaLo]; - double betaValueHigh = TableLoader.betaValues[betaHi]; - - //get the Bfield bins double B = Math.abs(Bf); - int binB = this.getBIdx(B); - double BfCen = BfieldValues[binB]; - int BfLo=0; - int BfHi=0; - if(BTableLoader.BfieldValues.length-1) { - BfHi=TableLoader.BfieldValues.length-1; + + int binlowB = this.getBIdx(B); + int binhighB = binlowB + 1; + + if(binhighB > TableLoader.maxBinIdxB) { + binhighB = TableLoader.maxBinIdxB; } - } - double BLo = BfieldValues[BfLo]; - double BHi = BfieldValues[BfHi]; - // get the alpha bins - int alphaBin = this.getAlphaIdx(alpha); - int alphaLo = 0; - int alphaHi = 0; - double alphaCenValue = this.getAlphaFromAlphaIdx(alphaBin); - if(alpha TableLoader.maxBinIdxAlpha) { - alphaHi = TableLoader.maxBinIdxAlpha; + + double B1 = TableLoader.BfieldValues[binlowB]; + double B2 = TableLoader.BfieldValues[binhighB]; + + // for alpha ranges + int binlowAlpha = this.getAlphaIdx(alpha); + int binhighAlpha = binlowAlpha + 1; + + if(binhighAlpha > TableLoader.maxBinIdxAlpha) { + binhighAlpha = TableLoader.maxBinIdxAlpha; } - } - double alphaValueLo = this.getAlphaFromAlphaIdx(alphaLo); - double alphaValueHi = this.getAlphaFromAlphaIdx(alphaHi); + //if(binhighAlpha==binlowAlpha) { + // binlowAlpha=binhighAlpha-1; + //} + + double alpha1 = this.getAlphaFromAlphaIdx(binlowAlpha); + double alpha2 = this.getAlphaFromAlphaIdx(binhighAlpha); // interpolate in B: - double Bfc = B*B; - double Bfa = BLo*BLo; - double Bfb = BHi*BHi; - double f_B_alpha1_beta1_t1 = interpolateLinear(Bfc, Bfa, Bfb, - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfLo][alphaLo][betaLo][tLo], - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfHi][alphaLo][betaLo][tLo]); - double f_B_alpha2_beta1_t1 = interpolateLinear(Bfc, Bfa, Bfb, - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfLo][alphaHi][betaLo][tLo], - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfHi][alphaHi][betaLo][tLo]); - double f_B_alpha1_beta1_t2 = interpolateLinear(Bfc, Bfa, Bfb, - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfLo][alphaLo][betaLo][tHi], - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfHi][alphaLo][betaLo][tHi]); - double f_B_alpha2_beta1_t2 = interpolateLinear(Bfc, Bfa, Bfb, - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfLo][alphaHi][betaLo][tHi], - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfHi][alphaHi][betaLo][tHi]); - double f_B_alpha1_beta2_t1 = interpolateLinear(Bfc, Bfa, Bfb, - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfLo][alphaLo][betaHi][tLo], - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfHi][alphaLo][betaHi][tLo]); - double f_B_alpha2_beta2_t1 = interpolateLinear(Bfc, Bfa, Bfb, - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfLo][alphaHi][betaHi][tLo], - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfHi][alphaHi][betaHi][tLo]); - double f_B_alpha1_beta2_t2 = interpolateLinear(Bfc, Bfa, Bfb, - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfLo][alphaLo][betaHi][tHi], - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfHi][alphaLo][betaHi][tHi]); - double f_B_alpha2_beta2_t2 = interpolateLinear(Bfc, Bfa, Bfb, - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfLo][alphaHi][betaHi][tHi], - TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][BfHi][alphaHi][betaHi][tHi]); + double f_B_alpha1_t1 = interpolateLinear(B*B, B1*B1, B2*B2, + TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binlowB][binlowAlpha][this.getTimeIdx(t, SecIdx, SlyrIdx, binlowB, binlowAlpha)], + TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binhighB][binlowAlpha][this.getTimeIdx(t, SecIdx, SlyrIdx, binhighB, binlowAlpha)]); + double f_B_alpha2_t1 = interpolateLinear(B*B, B1*B1, B2*B2, + TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binlowB][binhighAlpha][this.getTimeIdx(t, SecIdx, SlyrIdx, binlowB, binhighAlpha)], + TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binhighB][binhighAlpha][this.getTimeIdx(t, SecIdx, SlyrIdx, binhighB, binhighAlpha)]); + double f_B_alpha1_t2 = interpolateLinear(B*B, B1*B1, B2*B2, + TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binlowB][binlowAlpha][this.getTimeNextIdx(t, SecIdx, SlyrIdx, binlowB, binlowAlpha)], + TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binhighB][binlowAlpha][this.getTimeNextIdx(t, SecIdx, SlyrIdx, binhighB, binlowAlpha)]); + double f_B_alpha2_t2 = interpolateLinear(B*B, B1*B1, B2*B2, + TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binlowB][binhighAlpha][this.getTimeNextIdx(t, SecIdx, SlyrIdx, binlowB, binhighAlpha)], + TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binhighB][binhighAlpha][this.getTimeNextIdx(t, SecIdx, SlyrIdx, binhighB, binhighAlpha)]); + // interpolate in d for 2 values of alpha: + double f_B_alpha1_t = interpolateLinear(t, this.getTimeIdx(t, SecIdx, SlyrIdx, binlowB, binlowAlpha)*2., this.getTimeNextIdx(t, SecIdx, SlyrIdx, binhighB, binlowAlpha)*2., f_B_alpha1_t1, f_B_alpha1_t2); + double f_B_alpha2_t = interpolateLinear(t, this.getTimeIdx(t, SecIdx, SlyrIdx, binlowB, binhighAlpha)*2., this.getTimeNextIdx(t, SecIdx, SlyrIdx, binhighB, binhighAlpha)*2., f_B_alpha2_t1, f_B_alpha2_t2); + //LOGGER.log(Level.FINE, TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binlowB][binlowAlpha][this.getTimeIdx(t, SecIdx, SlyrIdx, binlowB, binlowAlpha)]); + //LOGGER.log(Level.FINE, SlyrIdx+" binlowB "+binlowB+" binlowAlpha "+binlowAlpha+" t "+this.getTimeIdx(t, SecIdx, SlyrIdx, binlowB, binlowAlpha)+" time "+t); + //LOGGER.log(Level.FINE, TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binlowB][binhighAlpha][this.getTimeIdx(t, SecIdx, SlyrIdx, binlowB, binhighAlpha)]); + //LOGGER.log(Level.FINE, SlyrIdx+" binlowB "+binlowB+" binhighAlpha "+binhighAlpha+" t "+this.getTimeIdx(t, SecIdx, SlyrIdx, binlowB, binhighAlpha)+" time "+t); + //LOGGER.log(Level.FINE, TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binhighB][binlowAlpha][this.getTimeIdx(t, SecIdx, SlyrIdx, binhighB, binlowAlpha)]); + //LOGGER.log(Level.FINE, SlyrIdx+" binhighB "+binhighB+" binlowAlpha "+binlowAlpha+" t "+this.getTimeIdx(t, SecIdx, SlyrIdx, binhighB, binlowAlpha)+" time "+t); + //LOGGER.log(Level.FINE, TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binhighB][binhighAlpha][this.getTimeIdx(t, SecIdx, SlyrIdx, binhighB, binhighAlpha)]); + //LOGGER.log(Level.FINE, SlyrIdx+" binhighB "+binhighB+" binhighAlpha "+binhighAlpha+" t "+this.getTimeIdx(t, SecIdx, SlyrIdx, binhighB, binhighAlpha)+" time "+t); + //LOGGER.log(Level.FINE, " f_B_alpha1_t1 "+f_B_alpha1_t1+" f_B_alpha2_t1 "+f_B_alpha2_t1 + // +" f_B_alpha1_t2 "+f_B_alpha1_t2+" f_B_alpha2_t2 "+f_B_alpha2_t2 + // +" f_B_alpha1_t "+f_B_alpha1_t+" f_B_alpha2_t "+f_B_alpha2_t); - //interpolate in alpha - double f_B_alpha_beta1_t1 = interpolateLinear(alpha, alphaValueLo, alphaValueHi, f_B_alpha1_beta1_t1, f_B_alpha2_beta1_t1); - double f_B_alpha_beta2_t1 = interpolateLinear(alpha, alphaValueLo, alphaValueHi, f_B_alpha1_beta2_t1, f_B_alpha2_beta2_t1); - double f_B_alpha_beta1_t2 = interpolateLinear(alpha, alphaValueLo, alphaValueHi, f_B_alpha1_beta1_t2, f_B_alpha2_beta1_t2); - double f_B_alpha_beta2_t2 = interpolateLinear(alpha, alphaValueLo, alphaValueHi, f_B_alpha1_beta2_t2, f_B_alpha2_beta2_t2); - //interpolate in beta - double f_B_alpha_beta_t1 = interpolateLinear(beta, betaValueLo, betaValueHigh,f_B_alpha_beta1_t1,f_B_alpha_beta2_t1); - double f_B_alpha_beta_t2 = interpolateLinear(beta, betaValueLo, betaValueHigh,f_B_alpha_beta1_t2,f_B_alpha_beta2_t2); - //interpolate in time - double f_B_alpha_beta_t = interpolateLinear(t, timeLo, timeHi, f_B_alpha_beta_t1, f_B_alpha_beta_t2); + // interpolate in alpha: (cos30-cosA) + double f_B_alpha_t = interpolateLinear(Math.cos(Math.toRadians(30.))-Math.cos(Math.toRadians(alpha)), + Math.cos(Math.toRadians(30.))-Math.cos(Math.toRadians(alpha1)), + Math.cos(Math.toRadians(30.))-Math.cos(Math.toRadians(alpha2)), f_B_alpha1_t, f_B_alpha2_t); - double x = f_B_alpha_beta_t; -// String st =new String(); -// st +="time "+t+" beta "+beta+" BETA: ["+betaValueLo+" "+betaCent+" "+betaValueHigh+"]"+" alpha "+alpha+"\n"; -// st +=(" f_B_alpha1_beta1_t1 "+f_B_alpha1_beta1_t1+" f_B_alpha2_beta1_t1 "+f_B_alpha2_beta1_t1+"\n"); -// st +=(" f_B_alpha1_beta2_t1 "+f_B_alpha1_beta2_t1+" f_B_alpha2_beta2_t1 "+f_B_alpha2_beta2_t1+"\n"); -// st +=(" f_B_alpha1_beta1_t2 "+f_B_alpha1_beta1_t2+" f_B_alpha2_beta1_t2 "+f_B_alpha2_beta1_t2+"\n"); -// st +=(" f_B_alpha1_beta2_t2 "+f_B_alpha1_beta2_t2+" f_B_alpha2_beta2_t2 "+f_B_alpha2_beta2_t2+"\n"); -// st +=(" f_B_alpha_beta1_t1 "+f_B_alpha_beta1_t1+" f_B_alpha_beta2_t1 "+f_B_alpha_beta2_t1+"\n"); -// st +=(" f_B_alpha_beta1_t2 "+f_B_alpha_beta1_t2+" f_B_alpha_beta2_t2 "+f_B_alpha_beta2_t2+"\n"); -// st +=(" f_B_alpha_beta_t1 "+f_B_alpha_beta_t1+" f_B_alpha_beta_t2 "+f_B_alpha_beta_t2+"\n"); -// st +=(" f_B_alpha_t "+f_B_alpha_beta_t+"\n"); + return f_B_alpha_t; - double dmax = 2.*Constants.getInstance().wpdist[SlyrIdx]; -// st +=(" x "+x+" dmax "+dmax+"\n"); + /* + // interpolate in B: + double f_B_alpha1_t = interpolateLinear(B*B, B1*B1, B2*B2, + TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binlowB][binlowAlpha][this.getTimeIdx(t, SecIdx, SlyrIdx, binlowB, binlowAlpha)], + TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binhighB][binlowAlpha][this.getTimeIdx(t, SecIdx, SlyrIdx, binhighB, binlowAlpha)]); + double f_B_alpha2_t = interpolateLinear(B*B, B1*B1, B2*B2, + TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binlowB][binhighAlpha][this.getTimeIdx(t, SecIdx, SlyrIdx, binlowB, binhighAlpha)], + TableLoader.DISTFROMTIME[SecIdx][SlyrIdx][binhighB][binhighAlpha][this.getTimeIdx(t, SecIdx, SlyrIdx, binhighB, binhighAlpha)]); - if(x>dmax) { -// setDebug(st); - return dmax; - } - if(!this.t2DPrecisionImprov) return x; - //Reolution improvement to compensate for non-linearity accross bin not accounted for in interpolation - double calctime = calc_Time( x, alpha, B, SecIdx+1, SlyrIdx+1) ; - double deltatime_beta = TableLoader.getDeltaTimeBeta(x,beta,TableLoader.distbeta[SecIdx][SlyrIdx],TableLoader.v0[SecIdx][SlyrIdx]); - calctime+=deltatime_beta; -// st +=(" t "+t+" calctime "+calctime+" deltatime_beta "+deltatime_beta+"\n"); - if(calctime>t) { - double tref=0; - for(int i = 1; i<100; i++) { - x-=0.0001*i; - if(x<0) return x; - calctime = calc_Time( x, alpha, B, SecIdx+1, SlyrIdx+1) ; - deltatime_beta = TableLoader.getDeltaTimeBeta(x,beta,TableLoader.distbeta[SecIdx][SlyrIdx],TableLoader.v0[SecIdx][SlyrIdx]); - calctime+=deltatime_beta; -// st +=(i+"] x "+x+" t "+t+" calct "+calctime+"\n"); - if(calctimecalctime) { - double tref=0; - for(int i = 1; i<100; i++) { - x+=0.0001*i; - calctime = calc_Time( x, alpha, B, SecIdx+1, SlyrIdx+1) ; - deltatime_beta = TableLoader.getDeltaTimeBeta(x,beta,TableLoader.distbeta[SecIdx][SlyrIdx],TableLoader.v0[SecIdx][SlyrIdx]); - calctime+=deltatime_beta; -// st +=(i+"] x "+x+" t "+t+" calct "+calctime+"\n"); - if(x>dmax) { -// setDebug(st); - return dmax; - } - if(tTableLoader.maxTBin) { + binIdx = TableLoader.maxTBin ; + } + return binIdx; + } /** * * @param b1 bfield value in T * @return B field bin */ public int getBIdx(double b1) { - int binIdx = (int) Math.floor(b1*b1); + +// int binIdx = (int) ((1+b1)*2) -2; +// if(binIdx<0) { +// binIdx = TableLoader.minBinIdxB; +// } +// if(binIdx>TableLoader.maxBinIdxB) { +// binIdx = TableLoader.maxBinIdxB; +// } int maxBinIdxB = TableLoader.BfieldValues.length-1; + DecimalFormat df = new DecimalFormat("#"); + df.setRoundingMode(RoundingMode.CEILING); + int binIdx =0; + try{ + binIdx = Integer.parseInt(df.format(b1*b1) ) -1; + } catch (NumberFormatException e) { + LOGGER.log(Level.WARNING, " field bin error "+b1+" "); + } if(binIdx<0) { binIdx = 0; } if(binIdx>maxBinIdxB) binIdx = maxBinIdxB; - return binIdx; } /** @@ -310,17 +221,33 @@ private int getAlphaIdx(double alpha) { return binIdx; } - private int getBetaIdx(double beta) { - if(beta>=1.0) return TableLoader.betaValues.length-1; - int value = TableLoader.betaValues.length-1; - for(int i = 0; i=TableLoader.betaValues[i] && betaTableLoader.maxBinIdxT[SecIdx][SlyrIdx][binlowB][binlowAlpha]) { + binhighT=TableLoader.maxBinIdxT[SecIdx][SlyrIdx][binlowB][binlowAlpha]; } + return binhighT; + } + + /** + * @param slyIdx superlayer index + * @param time + * @return test doca corr + */ + public double addDOCACorr(double time, int slyIdx) { + double dDoca = 0; + if(slyIdx+1 == 5 || slyIdx+1 ==6) { + if(time>600) { + dDoca = 0.15; + } else { + dDoca = (7.6e-3 - 2.4e-4*time +9.8e-3*time*time - 3.8e-6*time*time*time)*5.5410595e-05; + } + //LOGGER.log(Level.FINE, "time "+time +" added doca "+(float)dDoca); } - return value; + return dDoca; } +} - -} diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java index f90b26236a..28edab4186 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/TrackCandListFinder.java @@ -290,15 +290,9 @@ public List getTrackCandsAI(CrossList crossList, DCGeant4Factory DcDetect stateVec.getZ(), DcDetector, dcSwim); - if(cand.get_Vtx0() != null){ - Point3D VTCS = cand.get(cand.size()-1).getCoordsInTiltedSector(cand.get_Vtx0().x(), cand.get_Vtx0().y(), cand.get_Vtx0().z()); - finalSV.k = measSurfaces.size()-1; - kFZRef.finalStateVec = finalSV; - double deltaPathToVtx = kFZRef.getDeltaPathToVtx(sector, VTCS.z()); - - List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef, deltaPathToVtx); + List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef); cand.setStateVecs(kfStateVecsAlongTrajectory); - } + if (kFZRef.chi2 < Constants.MAXCHI2) cands.add(cand); } @@ -940,6 +934,8 @@ private List findStraightTracks(CrossList crossList, DCGeant4Factory DcDe continue; } + List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef); + StateVec fn = new StateVec(); if (!kFZRef.setFitFailed && kFZRef.finalStateVec != null) { fn.set(kFZRef.finalStateVec.x, kFZRef.finalStateVec.y, kFZRef.finalStateVec.tx, kFZRef.finalStateVec.ty); @@ -958,14 +954,7 @@ private List findStraightTracks(CrossList crossList, DCGeant4Factory DcDe cand.set_FitConvergenceStatus(kFZRef.ConvStatus); cand.set_Id(cands.size() + 1); cand.set_CovMat(kFZRef.finalStateVec.CM); - - if (cand.get_Vtx0() != null) { - Point3D VTCS = cand.get(cand.size()-1).getCoordsInTiltedSector(cand.get_Vtx0().x(), cand.get_Vtx0().y(), cand.get_Vtx0().z()); - double deltaPathToVtx = kFZRef.getDeltaPathToVtx(cand.get(cand.size()-1).get_Sector(), VTCS.z()); - - List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef, deltaPathToVtx); cand.setStateVecs(kfStateVecsAlongTrajectory); - } // add candidate to list of tracks cands.add(cand); @@ -1100,6 +1089,7 @@ private List findCurvedTracks(CrossList crossList, DCGeant4Factory DcDete kFZRef.init(measSurfaces, initSV); kFZRef.runFitter(); + List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef); if (kFZRef.finalStateVec != null) { if (kFZRef.chi2 < Constants.MAXCHI2) { @@ -1120,6 +1110,7 @@ private List findCurvedTracks(CrossList crossList, DCGeant4Factory DcDete cand.set_FitConvergenceStatus(kFZRef.ConvStatus); cand.set_CovMat(kFZRef.finalStateVec.CM); + cand.setStateVecs(kfStateVecsAlongTrajectory); cand.setFinalStateVec(fitStateVec); cand.set_Id(cands.size() + 1); @@ -1128,14 +1119,6 @@ private List findCurvedTracks(CrossList crossList, DCGeant4Factory DcDete fitStateVec.getZ(), DcDetector, dcSwim); - if (cand.get_Vtx0() != null) { - Point3D VTCS = cand.get(cand.size()-1).getCoordsInTiltedSector(cand.get_Vtx0().x(), cand.get_Vtx0().y(), cand.get_Vtx0().z()); - double deltaPathToVtx = kFZRef.getDeltaPathToVtx(cand.get(cand.size()-1).get_Sector(), VTCS.z()); - - List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef, deltaPathToVtx); - cand.setStateVecs(kfStateVecsAlongTrajectory); - } - // add candidate to list of tracks if (cand.fit_Successful = true) { cands.add(cand); @@ -1151,7 +1134,7 @@ private List findCurvedTracks(CrossList crossList, DCGeant4Factory DcDete return cands; } - public List setKFStateVecsAlongTrajectory(KFitter kFZRef, double deltaPathToVtx) { + public List setKFStateVecsAlongTrajectory(KFitter kFZRef) { List kfStateVecsAlongTrajectory = new ArrayList<>(); for(int i = 0; i < kFZRef.kfStateVecsAlongTrajectory.size(); i++) { @@ -1159,7 +1142,7 @@ public List setKFStateVecsAlongTrajectory(K org.jlab.rec.dc.trajectory.StateVec sv = new org.jlab.rec.dc.trajectory.StateVec(svc.x, svc.y, svc.tx, svc.ty); sv.setZ(svc.z); sv.setB(svc.B); - sv.setPathLength(svc.getPathLength() + deltaPathToVtx); // Transition for the starting point from the final point at the last layer to vertex + sv.setPathLength(svc.getPathLength()); sv.setProjector(svc.getProjector()); sv.setProjectorDoca(svc.getProjectorDoca()); kfStateVecsAlongTrajectory.add(sv); @@ -1168,7 +1151,7 @@ public List setKFStateVecsAlongTrajectory(K return kfStateVecsAlongTrajectory; } - public List setKFStateVecsAlongTrajectory(KFitterStraight kFZRef, double deltaPathToVtx) { + public List setKFStateVecsAlongTrajectory(KFitterStraight kFZRef) { List kfStateVecsAlongTrajectory = new ArrayList<>(); for(int i = 0; i < kFZRef.kfStateVecsAlongTrajectory.size(); i++) { @@ -1176,7 +1159,7 @@ public List setKFStateVecsAlongTrajectory(K org.jlab.rec.dc.trajectory.StateVec sv = new org.jlab.rec.dc.trajectory.StateVec(svc.x, svc.y, svc.tx, svc.ty); sv.setZ(svc.z); sv.setB(svc.B); - sv.setPathLength(svc.getPathLength() + deltaPathToVtx); // Transition for the starting point from the final point at the last layer to vertex + sv.setPathLength(svc.getPathLength()); sv.setProjector(svc.getProjector()); sv.setProjectorDoca(svc.getProjectorDoca()); kfStateVecsAlongTrajectory.add(sv); diff --git a/reconstruction/dc/src/main/java/org/jlab/service/dc/DCTBEngine.java b/reconstruction/dc/src/main/java/org/jlab/service/dc/DCTBEngine.java index 4bf0281c17..8baf5f2ca8 100644 --- a/reconstruction/dc/src/main/java/org/jlab/service/dc/DCTBEngine.java +++ b/reconstruction/dc/src/main/java/org/jlab/service/dc/DCTBEngine.java @@ -258,6 +258,7 @@ public boolean processDataEventUser(DataEvent event) { getInitState(TrackArray1, measSurfaces.get(0).measPoint.z(), initSV, kFZRef, dcSwim, new float[3]); kFZRef.initFromHB(measSurfaces, initSV, TrackArray1.get(0).get(0).get(0).get_Beta()); kFZRef.runFitter(); + List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef); StateVec fn = new StateVec(); if (kFZRef.setFitFailed==false && kFZRef.finalStateVec!=null) { @@ -275,6 +276,7 @@ public boolean processDataEventUser(DataEvent event) { TrackArray1.set_FitChi2(kFZRef.chi2); TrackArray1.set_FitNDF(kFZRef.NDF); + TrackArray1.setStateVecs(kfStateVecsAlongTrajectory); TrackArray1.set_FitConvergenceStatus(kFZRef.ConvStatus); if (TrackArray1.get_Vtx0().toVector3D().mag() > 500) { continue; @@ -285,10 +287,6 @@ public boolean processDataEventUser(DataEvent event) { TrackArray1.set_CovMat(kFZRef.propagateToVtx(TrackArray1.get(TrackArray1.size()-1).get_Sector(), VTCS.z())); TrackArray1.transCMToGlobal(); - double deltaPathToVtx = kFZRef.getDeltaPathToVtx(TrackArray1.get(TrackArray1.size()-1).get_Sector(), VTCS.z()); - List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef, deltaPathToVtx); - TrackArray1.setStateVecs(kfStateVecsAlongTrajectory); - if (TrackArray1.isGood()) { trkcands.add(TrackArray1); } @@ -303,6 +301,8 @@ public boolean processDataEventUser(DataEvent event) { kFZRef.initFromHB(measSurfaces, initSV, TrackArray1.get(0).get(0).get(0).get_Beta(), useDAF); kFZRef.runFitter(useDAF); + List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef); + StateVec fn = new StateVec(); if (kFZRef.setFitFailed==false && kFZRef.finalStateVec!=null) { // set the state vector at the last measurement site @@ -320,6 +320,7 @@ public boolean processDataEventUser(DataEvent event) { TrackArray1.set_FitChi2(kFZRef.chi2); TrackArray1.set_FitNDF(kFZRef.NDF); TrackArray1.set_NDFDAF(kFZRef.getNDFDAF()); + TrackArray1.setStateVecs(kfStateVecsAlongTrajectory); TrackArray1.set_FitConvergenceStatus(kFZRef.ConvStatus); if (TrackArray1.get_Vtx0().toVector3D().mag() > 500) { continue; @@ -330,10 +331,6 @@ public boolean processDataEventUser(DataEvent event) { TrackArray1.set_CovMat(kFZRef.propagateToVtx(TrackArray1.get(TrackArray1.size()-1).get_Sector(), VTCS.z())); TrackArray1.transCMToGlobal(); - double deltaPathToVtx = kFZRef.getDeltaPathToVtx(TrackArray1.get(TrackArray1.size()-1).get_Sector(), VTCS.z()); - List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef, deltaPathToVtx); - TrackArray1.setStateVecs(kfStateVecsAlongTrajectory); - if (TrackArray1.isGood()) { trkcands.add(TrackArray1); } @@ -389,7 +386,7 @@ public boolean processDataEventUser(DataEvent event) { return true; } - public List setKFStateVecsAlongTrajectory(KFitter kFZRef, double deltaPathToVtx) { + public List setKFStateVecsAlongTrajectory(KFitter kFZRef) { List kfStateVecsAlongTrajectory = new ArrayList<>(); for(int i = 0; i < kFZRef.kfStateVecsAlongTrajectory.size(); i++) { @@ -397,7 +394,7 @@ public List setKFStateVecsAlongTrajectory(K org.jlab.rec.dc.trajectory.StateVec sv = new org.jlab.rec.dc.trajectory.StateVec(svc.x, svc.y, svc.tx, svc.ty); sv.setZ(svc.z); sv.setB(svc.B); - sv.setPathLength(svc.getPathLength() + deltaPathToVtx); // Transition for the starting point from the final point at the last layer to vertex + sv.setPathLength(svc.getPathLength()); sv.setProjector(svc.getProjector()); sv.setProjectorDoca(svc.getProjectorDoca()); sv.setDAFWeight(svc.getFinalDAFWeight()); @@ -408,7 +405,7 @@ public List setKFStateVecsAlongTrajectory(K return kfStateVecsAlongTrajectory; } - public List setKFStateVecsAlongTrajectory(KFitterStraight kFZRef, double deltaPathToVtx) { + public List setKFStateVecsAlongTrajectory(KFitterStraight kFZRef) { List kfStateVecsAlongTrajectory = new ArrayList<>(); for(int i = 0; i < kFZRef.kfStateVecsAlongTrajectory.size(); i++) { @@ -416,7 +413,7 @@ public List setKFStateVecsAlongTrajectory(K org.jlab.rec.dc.trajectory.StateVec sv = new org.jlab.rec.dc.trajectory.StateVec(svc.x, svc.y, svc.tx, svc.ty); sv.setZ(svc.z); sv.setB(svc.B); - sv.setPathLength(svc.getPathLength() + deltaPathToVtx); // Transition for the starting point from the final point at the last layer to vertex + sv.setPathLength(svc.getPathLength()); sv.setProjector(svc.getProjector()); sv.setProjectorDoca(svc.getProjectorDoca()); kfStateVecsAlongTrajectory.add(sv);