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 c723fb2a3f..cfdce11c30 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 @@ -708,7 +708,9 @@ private boolean filter(int k, boolean forward, double annealingFactor) { return 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 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); @@ -857,7 +859,9 @@ public Matrix filterCovMat(double[] H, Matrix Ci, double V) { 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; @@ -880,9 +884,9 @@ 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 += svc.deltaPath; + path += (forward ? 1 : -1) * svc.deltaPath; svc.setPathLength(path); - + double V0 = mv.measurements.get(0).surface.unc[0]; Point3D point = new Point3D(svc.x, svc.y, mv.measurements.get(0).surface.measPoint.z()); @@ -922,7 +926,7 @@ private void calcFinalChisq(int sector, boolean nofilter) { double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0]); svc = sv.transported(forward).get(k1 + 1); - path += svc.deltaPath; + path += (forward ? 1 : -1) * svc.deltaPath; svc.setPathLength(path); svc.setProjector(mv.measurements.get(k1 + 1).surface.wireLine[0].origin().x()); svc.setProjectorDoca(h); @@ -975,7 +979,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 += svc.deltaPath; + path += (forward ? 1 : -1) * svc.deltaPath; svc.setPathLength(path); Point3D point = new Point3D(svc.x, svc.y, mv.measurements.get(0).surface.measPoint.z()); @@ -1047,7 +1051,7 @@ private void calcFinalChisqDAF(int sector, boolean nofilter) { } svc = sv.transported(forward).get(k1 + 1); - path += svc.deltaPath; + path += (forward ? 1 : -1) * 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()); @@ -1116,6 +1120,10 @@ private void calcFinalChisqDAF(int sector, boolean nofilter) { 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 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 232bcaa7d0..fb563fcc8b 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 @@ -406,6 +406,8 @@ 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; @@ -426,11 +428,11 @@ private void calcFinalChisq(int sector, boolean nofilter) { kfStateVecsAlongTrajectory = new ArrayList<>(); if (sVec != null && sVec.CM != null) { - boolean forward = false; + boolean forward = false; sv.transport(sector, k, 0, sVec, mv, this.getSwimmer(), forward); StateVec svc = sv.transported(forward).get(0); - path += svc.deltaPath; + path += (forward ? 1 : -1) * svc.deltaPath; svc.setPathLength(path); double V0 = mv.measurements.get(0).surface.unc[0]; @@ -473,7 +475,7 @@ private void calcFinalChisq(int sector, boolean nofilter) { double h = mv.hDoca(point, mv.measurements.get(k1 + 1).surface.wireLine[0]); svc = sv.transported(forward).get(k1+1); - path += svc.deltaPath; + path += (forward ? 1 : -1) * svc.deltaPath; svc.setPathLength(path); svc.setProjector(mv.measurements.get(k1 + 1).surface.wireLine[0].origin().x()); svc.setProjectorDoca(h); @@ -504,6 +506,10 @@ 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 f69dd0bfc9..269e966c21 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 @@ -53,7 +53,7 @@ public void initFromHB(StateVec initSV, double beta) { this.trackTrajS.clear(); this.trackTrajT.put(0, new StateVec(initSV)); } - + /** * * @param sector @@ -150,6 +150,104 @@ public Matrix transport(int sector, int i, double Zf, StateVec iVec, AMeasVecs m return fVec.CM; + } + + /** + * + * @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.FINE, " 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; + } /** 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 d020322477..25473ca915 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 @@ -94,7 +94,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; @@ -105,10 +105,10 @@ 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/time_to_distance/t2d_pressure"; + public static final String T2DPRESSURE = "/calibration/dc/v2/t2d_pressure"; public static final String PRESSURE = "/hall/weather/pressure"; - 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 T2DPRESSUREREF = "/calibration/dc/v2/ref_pressure"; + public static final String T0CORRECTION = "/calibration/dc/v2/t0"; 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"; @@ -296,13 +296,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; @@ -429,7 +429,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 faa907a8cf..102e9a4b47 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 @@ -791,8 +791,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 36e802de15..19f27b1c50 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 @@ -385,27 +385,28 @@ 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); +// 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); } diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/timetodistance/T2DFunctions.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/timetodistance/T2DFunctions.java index 9a863dbcce..4892d5ceb8 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/timetodistance/T2DFunctions.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/timetodistance/T2DFunctions.java @@ -321,7 +321,7 @@ public static double polyFcnMac(double x, double alpha, double bfield, double v_ double time = 0; // alpha correction double cos30minusalpha=Math.cos(Math.toRadians(30.-alpha)); - double dmaxalpha = dmax*cos30minusalpha; + double dmaxalpha = dmax*cos30minusalpha; double xhatalpha = x/dmaxalpha; // rcapital is an intermediate parameter double rcapital = R*dmax; 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 1ae2e32755..1309a42012 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,10 +1,21 @@ package org.jlab.rec.dc.timetodistance; -import java.math.RoundingMode; -import java.text.DecimalFormat; +import java.util.ArrayList; +import java.util.HashMap; +import java.util.List; +import java.util.Map; import java.util.logging.Level; import java.util.logging.Logger; +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; @@ -17,8 +28,8 @@ public TableLoader() { private static boolean T2DLOADED = false; - private static final int NBINST=2000; - + 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; @@ -28,39 +39,185 @@ 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][NBINST]; // sector slyr alpha Bfield time bins [s][r][ibfield][icosalpha][tbin] + 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 int maxTBin = -1; - //public static double[] distbetaValues = new double[]{0.16, 0.16, 0.08, 0.08, 0.08, 0.08}; - /* - * - */ + 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(){ - TimeToDistanceEstimator tde = new TimeToDistanceEstimator(); - for(int s = 0; s<1; s++ ){ // loop over sectors - for(int r = 4; r<5; r++ ){ //loop over slys - for(int ibfield =0; ibfield<1; ibfield++) { - for (int tb = 250; tb< 300; tb++) { - LOGGER.log(Level.FINE, " NEW TIME BIN "); - for(int icosalpha =0; icosalpha 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) { @@ -93,148 +250,135 @@ private static synchronized void FillAlpha() { AlphaBounds[0][0] = 0; AlphaBounds[5][1] = 30; } - public static synchronized void Fill(IndexedTable t2dPressure, IndexedTable t2dPressRef, IndexedTable pressure) { - if (T2DLOADED) return; - - + public static boolean useP = true; + public static synchronized void getConstants(IndexedTable t2dPressure, IndexedTable t2dPressRef, IndexedTable pressure){ double p_ref = t2dPressRef.getDoubleValue("pressure", 0,0,0); double p = pressure.getDoubleValue("value", 0,0,3); double dp = p - p_ref; + double dp2scale = 0; + double dpscale = 1; + + if(!useP) + dpscale = 0; for(int s = 0; s<6; s++ ){ // loop over sectors for(int r = 0; r<6; r++ ){ //loop over slys // Fill constants FracDmaxAtMinVel[s][r] = t2dPressure.getDoubleValue("c1_a0", s+1,r+1,0) - +t2dPressure.getDoubleValue("c1_a1", s+1,r+1,0)*dp; + +t2dPressure.getDoubleValue("c1_a1", s+1,r+1,0)*dp*dpscale; v0[s][r] = t2dPressure.getDoubleValue("v0_a0", s+1,r+1,0) - +t2dPressure.getDoubleValue("v0_a1", s+1,r+1,0)*dp - +t2dPressure.getDoubleValue("v0_a2", s+1,r+1,0)*dp*dp; + +t2dPressure.getDoubleValue("v0_a1", s+1,r+1,0)*dp*dpscale + +t2dPressure.getDoubleValue("v0_a2", s+1,r+1,0)*dp*dp*dp2scale; vmid[s][r] = t2dPressure.getDoubleValue("vmid_a0", s+1,r+1,0) - +t2dPressure.getDoubleValue("vmid_a1", s+1,r+1,0)*dp - +t2dPressure.getDoubleValue("vmid_a2", s+1,r+1,0)*dp*dp; - delta_bfield_coefficient[s][r] = t2dPressure.getDoubleValue("delta_bfield_a0", s+1,r+1,0) - +t2dPressure.getDoubleValue("delta_bfield_a1", s+1,r+1,0)*dp - +t2dPressure.getDoubleValue("delta_bfield_a2", s+1,r+1,0)*dp*dp; - b1[s][r] = t2dPressure.getDoubleValue("b1_a0", s+1,r+1,0) - +t2dPressure.getDoubleValue("b1_a1", s+1,r+1,0)*dp - +t2dPressure.getDoubleValue("b1_a2", s+1,r+1,0)*dp*dp; - b2[s][r] = t2dPressure.getDoubleValue("b2_a0", s+1,r+1,0) - +t2dPressure.getDoubleValue("b2_a1", s+1,r+1,0)*dp - +t2dPressure.getDoubleValue("b2_a2", s+1,r+1,0)*dp*dp; - b3[s][r] = t2dPressure.getDoubleValue("b3_a0", s+1,r+1,0) - +t2dPressure.getDoubleValue("b3_a1", s+1,r+1,0)*dp - +t2dPressure.getDoubleValue("b3_a2", s+1,r+1,0)*dp*dp; - b4[s][r] = t2dPressure.getDoubleValue("b4_a0", s+1,r+1,0) - +t2dPressure.getDoubleValue("b4_a1", s+1,r+1,0)*dp - +t2dPressure.getDoubleValue("b4_a2", s+1,r+1,0)*dp*dp; + +t2dPressure.getDoubleValue("vmid_a1", s+1,r+1,0)*dp*dpscale + +t2dPressure.getDoubleValue("vmid_a2", s+1,r+1,0)*dp*dp*dp2scale; + distbeta[s][r] = t2dPressure.getDoubleValue("distbeta_a0", s+1,r+1,0) + +t2dPressure.getDoubleValue("distbeta_a1", s+1,r+1,0)*dp*dpscale + +t2dPressure.getDoubleValue("distbeta_a2", s+1,r+1,0)*dp*dp*dp2scale; + if(r>1 && r<4) { + delta_bfield_coefficient[s][r] = t2dPressure.getDoubleValue("delta_bfield_a0", s+1,r+1,0) + +t2dPressure.getDoubleValue("delta_bfield_a1", s+1,r+1,0)*dp*dpscale + +t2dPressure.getDoubleValue("delta_bfield_a2", s+1,r+1,0)*dp*dp*dp2scale + +t2dPressure.getDoubleValue("delta_bfield_a1", s+1,r+1,0)*dp*dpscale + +t2dPressure.getDoubleValue("delta_bfield_a2", s+1,r+1,0)*dp*dp*dp2scale; + b1[s][r] = t2dPressure.getDoubleValue("b1_a0", s+1,r+1,0) + +t2dPressure.getDoubleValue("b1_a1", s+1,r+1,0)*dp*dpscale + +t2dPressure.getDoubleValue("b1_a2", s+1,r+1,0)*dp*dp*dp2scale; + b2[s][r] = t2dPressure.getDoubleValue("b2_a0", s+1,r+1,0) + +t2dPressure.getDoubleValue("b2_a1", s+1,r+1,0)*dp*dpscale + +t2dPressure.getDoubleValue("b2_a2", s+1,r+1,0)*dp*dp*dp2scale; + b3[s][r] = t2dPressure.getDoubleValue("b3_a0", s+1,r+1,0) + +t2dPressure.getDoubleValue("b3_a1", s+1,r+1,0)*dp*dpscale + +t2dPressure.getDoubleValue("b3_a2", s+1,r+1,0)*dp*dp*dp2scale; + b4[s][r] = t2dPressure.getDoubleValue("b4_a0", s+1,r+1,0) + +t2dPressure.getDoubleValue("b4_a1", s+1,r+1,0)*dp*dpscale + +t2dPressure.getDoubleValue("b4_a2", s+1,r+1,0)*dp*dp*dp2scale; + } Tmax[s][r] = t2dPressure.getDoubleValue("tmax_a0", s+1,r+1,0) - +t2dPressure.getDoubleValue("tmax_a1", s+1,r+1,0)*dp - +t2dPressure.getDoubleValue("tmax_a2", s+1,r+1,0)*dp*dp; - + +t2dPressure.getDoubleValue("tmax_a1", s+1,r+1,0)*dp*dpscale + +t2dPressure.getDoubleValue("tmax_a2", s+1,r+1,0)*dp*dp*dp2scale; } } - Fill(); } - public static synchronized void Fill(IndexedTable tab) { - //CCDBTables 0 = "/calibration/dc/signal_generation/doca_resolution"; - //CCDBTables 1 = "/calibration/dc/time_to_distance/t2d"; - //CCDBTables 2 = "/calibration/dc/time_corrections/T0_correction"; - if (T2DLOADED) return; - + public static synchronized void FillTable() { + double stepSize = 0.00010; for(int s = 0; s<6; s++ ){ // loop over sectors - for(int r = 0; r<6; r++ ){ //loop over slys - // Fill constants - delta_T0[s][r] = tab.getDoubleValue("delta_T0", s+1,r+1,0); - FracDmaxAtMinVel[s][r] = tab.getDoubleValue("c1", s+1,r+1,0);//use same table. names strings - deltanm[s][r] = tab.getDoubleValue("deltanm", s+1,r+1,0); - v0[s][r] = tab.getDoubleValue("v0", s+1,r+1,0); - vmid[s][r] = tab.getDoubleValue("c2", s+1,r+1,0); - delta_bfield_coefficient[s][r] = tab.getDoubleValue("delta_bfield_coefficient", s+1,r+1,0); - b1[s][r] = tab.getDoubleValue("b1", s+1,r+1,0); - b2[s][r] = tab.getDoubleValue("b2", s+1,r+1,0); - b3[s][r] = tab.getDoubleValue("b3", s+1,r+1,0); - b4[s][r] = tab.getDoubleValue("b4", s+1,r+1,0); - Tmax[s][r] = tab.getDoubleValue("tmax", s+1,r+1,0); - // end fill constants + //System.out.println("sector "+(s+1)+" sly "+(r+1)+" v0 "+v0[s][r]+" vmid "+vmid[s][r]+" R "+FracDmaxAtMinVel[s][r]); + double dmax = 2.*Constants.getInstance().wpdist[r]; + //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"); + continue; + } + if(tbin>maxTBin) + maxTBin = tbin; + + if(timebfielddmax) { + DISTFROMTIME[s][r][ibfield][icosalpha][ibeta][tbin]=dmax; + idist=nxmax; + } + } + } + } + } } } - Fill(); - } - public static synchronized void Fill() { + TableLoader.fillMissingTableBins(); + } + public static synchronized void Fill(IndexedTable t2dPressure, IndexedTable t2dPressRef, IndexedTable pressure) { + //CCDBTables 0 = "/calibration/dc/signal_generation/doca_resolution"; //CCDBTables 1 = "/calibration/dc/time_to_distance/t2d"; //CCDBTables 2 = "/calibration/dc/time_corrections/T0_correction"; - - double stepSize = 0.0010; + if (T2DLOADED) return; FillAlpha(); + getConstants(t2dPressure, t2dPressRef, pressure); + FillTable(); - for(int s = 0; s<6; s++ ){ // loop over sectors - - for(int r = 0; r<6; r++ ){ //loop over slys - - // constants filled - //LOGGER.log(Level.FINE, v0[s][r]+" "+vmid[s][r]+" "+FracDmaxAtMinVel[s][r]); - double dmax = 2.*Constants.getInstance().wpdist[r]; - //double tmax = CCDBConstants.getTMAXSUPERLAYER()[s][r]; - for(int ibfield =0; ibfieldNBINST-1) { - //System.err.println("Problem with tbin"); - continue; - } - if(tbin>maxTBin) - maxTBin = tbin; - //if(tbin>maxBinIdxT[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(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; - } */ - } - } - } - } - } - TableLoader.fillMissingTableBins(); - //TableLoader.test(); + System.out.println(" T2D TABLE FILLED....."); T2DLOADED = true; - } + } - private static synchronized void fillMissingTableBins() { + private static void fillMissingTableBins() { for(int s = 0; s<6; s++ ){ // loop over sectors @@ -244,17 +388,21 @@ private static synchronized void fillMissingTableBins() { for(int icosalpha =0; icosalphadmax) x=dmax; - - if(Constants.getInstance().getT2D()==0) { - - return T2DFunctions.ExpoFcn(x, alpha, bfield, v0[s][r], deltanm[s][r], 0.615, - tmax, dmax, delBf, Bb1, Bb2, Bb3, Bb4, superlayer) + delta_T0[s][r]; - } else { - return T2DFunctions.polyFcnMac(x, alpha, bfield, v0[s][r], vmid[s][r], FracDmaxAtMinVel[s][r], + return polyFcnMac(x, alpha, bfield, v0[s][r], vmid[s][r], FracDmaxAtMinVel[s][r], tmax, dmax, delBf, Bb1, Bb2, Bb3, Bb4, superlayer) ; - } + + } + + public static synchronized double getDeltaTimeBeta(double x, double beta, double distbeta, double v_0) { + + double value = (0.5*Math.pow(beta*beta*distbeta,3)*x/(Math.pow(beta*beta*distbeta,3)+x*x*x))/v_0; + + return value; } + public static double[][] delta_T0 = new double[6][6]; public static double[][] delta_bfield_coefficient = new double[6][6]; - public static double[][] deltanm = new double[6][6]; + public static double[][] distbeta = new double[6][6]; public static double[][] vmid = new double[6][6]; public static double[][] v0 = new double[6][6]; public static double[][] b1 = new double[6][6]; @@ -297,6 +447,31 @@ public static synchronized double calc_Time(double x, double alpha, double bfiel public static double[][] b3 = new double[6][6]; public static double[][] b4 = new double[6][6]; public static double[][] Tmax = new double[6][6]; - public static double[][] FracDmaxAtMinVel = new double[6][6]; // fraction of dmax corresponding to the point in the cell where the velocity is minimal + public static double[][] FracDmaxAtMinVel = new double[6][6]; + private static class FitLine extends Func1D { + + public FitLine(String name, double min, double max, + int s, int r, int ibeta, double alpha, double bfield) { + super(name, min, max); + this.s = s; + this.r = r; + this.ibeta = ibeta; + this.alpha = alpha; + this.bfield = bfield; + } + private final double alpha; + private final double bfield; + private final int ibeta; + private final int s; + private final int r; + + @Override + public double evaluate(double x) { + double timebfield = calc_Time( x, alpha, bfield, s+1, r+1) ; + double deltatime_beta = getDeltaTimeBeta(x,betaValues[ibeta],distbeta[s][r],v0[s][r]); + timebfield+=deltatime_beta; + return timebfield; + } + } } diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/timetodistance/TimeToDistanceEstimator.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/timetodistance/TimeToDistanceEstimator.java index 141a631ed3..998541d165 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/timetodistance/TimeToDistanceEstimator.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/timetodistance/TimeToDistanceEstimator.java @@ -1,17 +1,19 @@ package org.jlab.rec.dc.timetodistance; -import java.math.RoundingMode; -import java.text.DecimalFormat; import java.util.logging.Level; import java.util.logging.Logger; +import org.jlab.rec.dc.Constants; +import static org.jlab.rec.dc.timetodistance.TableLoader.BfieldValues; +import static org.jlab.rec.dc.timetodistance.TableLoader.calc_Time; public class TimeToDistanceEstimator { + public TimeToDistanceEstimator() { + // TODO Auto-generated constructor stub } - - private static final Logger LOGGER = Logger.getLogger(TimeToDistanceEstimator.class.getName()); + public boolean t2DPrecisionImprov=false; /** * * @param x value on grid @@ -27,104 +29,224 @@ private double interpolateLinear(double x0, double xa, double xb, double ya, dou x=xb; if(xTableLoader.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 binlowB = this.getBIdx(B); - int binhighB = binlowB + 1; - - if(binhighB > TableLoader.maxBinIdxB) { - binhighB = TableLoader.maxBinIdxB; + 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; + } } - - 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 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; + } } - //if(binhighAlpha==binlowAlpha) { - // binlowAlpha=binhighAlpha-1; - //} - - double alpha1 = this.getAlphaFromAlphaIdx(binlowAlpha); - double alpha2 = this.getAlphaFromAlphaIdx(binhighAlpha); + double alphaValueLo = this.getAlphaFromAlphaIdx(alphaLo); + double alphaValueHi = this.getAlphaFromAlphaIdx(alphaHi); // interpolate in B: - 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); + 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]); - // 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); + //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); + + 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"); + + double dmax = 2.*Constants.getInstance().wpdist[SlyrIdx]; +// st +=(" x "+x+" dmax "+dmax+"\n"); + + 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; + public int getTimeIdx(double t1) { + int idx=(int) Math.floor(t1 / TableLoader.timeBinWidth); + if(idx<0) idx=0; + return idx; } + /** * * @param b1 bfield value in T * @return B field bin */ - public int getBIdx(double b1) { - -// int binIdx = (int) ((1+b1)*2) -2; -// if(binIdx<0) { -// binIdx = TableLoader.minBinIdxB; -// } -// if(binIdx>TableLoader.maxBinIdxB) { -// binIdx = TableLoader.maxBinIdxB; -// } + public int getBIdx(double b1) { + int binIdx = (int) Math.floor(b1*b1); 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; } /** @@ -220,34 +311,18 @@ private int getAlphaIdx(double alpha) { } return binIdx; } - - private int getTimeNextIdx(double t, int SecIdx, int SlyrIdx, int binlowB, int binlowAlpha) { - int binlowT = this.getTimeIdx(t, SecIdx, SlyrIdx, binlowB, binlowAlpha); - int binhighT = binlowT + 1; - - if(binhighT>TableLoader.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; + 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] && beta 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); @@ -927,6 +925,11 @@ private List findStraightTracks(CrossList crossList, DCGeant4Factory DcDe cand.set_FitConvergenceStatus(kFZRef.ConvStatus); cand.set_Id(cands.size() + 1); cand.set_CovMat(kFZRef.finalStateVec.CM); + + 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 @@ -1071,7 +1074,6 @@ private List findCurvedTracks(CrossList crossList, DCGeant4Factory DcDete kFZRef.init(measSurfaces, initSV); kFZRef.runFitter(); - List kfStateVecsAlongTrajectory = setKFStateVecsAlongTrajectory(kFZRef); if (kFZRef.finalStateVec == null) { continue; @@ -1093,15 +1095,21 @@ private List findCurvedTracks(CrossList crossList, DCGeant4Factory DcDete cand.set_FitNDF(kFZRef.NDF); cand.set_FitConvergenceStatus(kFZRef.ConvStatus); - cand.set_CovMat(kFZRef.finalStateVec.CM); - cand.setStateVecs(kfStateVecsAlongTrajectory); - + cand.set_CovMat(kFZRef.finalStateVec.CM); + cand.setFinalStateVec(fitStateVec); cand.set_Id(cands.size() + 1); this.setTrackPars(cand, traj, trjFind, fitStateVec, fitStateVec.getZ(), DcDetector, dcSwim); + + 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); @@ -1117,7 +1125,7 @@ private List findCurvedTracks(CrossList crossList, DCGeant4Factory DcDete return cands; } - public List setKFStateVecsAlongTrajectory(KFitter kFZRef) { + public List setKFStateVecsAlongTrajectory(KFitter kFZRef, double deltaPathToVtx) { List kfStateVecsAlongTrajectory = new ArrayList<>(); for(int i = 0; i < kFZRef.kfStateVecsAlongTrajectory.size(); i++) { @@ -1125,7 +1133,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()); + sv.setPathLength(svc.getPathLength() + deltaPathToVtx); // Transition for the starting point from the final point at the last layer to vertex sv.setProjector(svc.getProjector()); sv.setProjectorDoca(svc.getProjectorDoca()); kfStateVecsAlongTrajectory.add(sv); @@ -1134,7 +1142,7 @@ public List setKFStateVecsAlongTrajectory(K return kfStateVecsAlongTrajectory; } - public List setKFStateVecsAlongTrajectory(KFitterStraight kFZRef) { + public List setKFStateVecsAlongTrajectory(KFitterStraight kFZRef, double deltaPathToVtx) { List kfStateVecsAlongTrajectory = new ArrayList<>(); for(int i = 0; i < kFZRef.kfStateVecsAlongTrajectory.size(); i++) { @@ -1142,7 +1150,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()); + sv.setPathLength(svc.getPathLength() + deltaPathToVtx); // Transition for the starting point from the final point at the last layer to vertex 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 2bd9d1e9f9..5ecf48569d 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 @@ -90,13 +90,13 @@ public boolean processDataEvent(DataEvent event) { Swim dcSwim = new Swim(); // fill T2D table - if(Constants.getInstance().getT2D()==0) { - TableLoader.Fill(this.getConstantsManager().getConstants(run, Constants.TIME2DIST)); - } else { +// if(Constants.getInstance().getT2D()==0) { +// TableLoader.Fill(this.getConstantsManager().getConstants(run, Constants.TIME2DIST)); +// } else { TableLoader.Fill(this.getConstantsManager().getConstants(run, Constants.T2DPRESSURE), this.getConstantsManager().getConstants(run, Constants.T2DPRESSUREREF), this.getConstantsManager().getConstants(run, Constants.PRESSURE)); - } +// } ClusterFitter cf = new ClusterFitter(); ClusterCleanerUtilities ct = new ClusterCleanerUtilities(); @@ -235,7 +235,6 @@ public boolean processDataEvent(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) { @@ -253,15 +252,18 @@ public boolean processDataEvent(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; } - + // get CovMat at vertex - Point3D VTCS = crosses.get(0).getCoordsInSector(TrackArray1.get_Vtx0().x(), TrackArray1.get_Vtx0().y(), TrackArray1.get_Vtx0().z()); + Point3D VTCS = crosses.get(0).getCoordsInTiltedSector(TrackArray1.get_Vtx0().x(), TrackArray1.get_Vtx0().y(), TrackArray1.get_Vtx0().z()); TrackArray1.set_CovMat(kFZRef.propagateToVtx(crosses.get(0).get_Sector(), VTCS.z())); + + 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); @@ -277,8 +279,6 @@ public boolean processDataEvent(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 @@ -296,15 +296,18 @@ public boolean processDataEvent(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; } // get CovMat at vertex - Point3D VTCS = crosses.get(0).getCoordsInSector(TrackArray1.get_Vtx0().x(), TrackArray1.get_Vtx0().y(), TrackArray1.get_Vtx0().z()); + Point3D VTCS = crosses.get(0).getCoordsInTiltedSector(TrackArray1.get_Vtx0().x(), TrackArray1.get_Vtx0().y(), TrackArray1.get_Vtx0().z()); TrackArray1.set_CovMat(kFZRef.propagateToVtx(crosses.get(0).get_Sector(), VTCS.z())); + + 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); @@ -361,7 +364,7 @@ public boolean processDataEvent(DataEvent event) { return true; } - public List setKFStateVecsAlongTrajectory(KFitter kFZRef) { + public List setKFStateVecsAlongTrajectory(KFitter kFZRef, double deltaPathToVtx) { List kfStateVecsAlongTrajectory = new ArrayList<>(); for(int i = 0; i < kFZRef.kfStateVecsAlongTrajectory.size(); i++) { @@ -369,7 +372,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()); + sv.setPathLength(svc.getPathLength() + deltaPathToVtx); // Transition for the starting point from the final point at the last layer to vertex sv.setProjector(svc.getProjector()); sv.setProjectorDoca(svc.getProjectorDoca()); sv.setDAFWeight(svc.getFinalDAFWeight()); @@ -380,7 +383,7 @@ public List setKFStateVecsAlongTrajectory(K return kfStateVecsAlongTrajectory; } - public List setKFStateVecsAlongTrajectory(KFitterStraight kFZRef) { + public List setKFStateVecsAlongTrajectory(KFitterStraight kFZRef, double deltaPathToVtx) { List kfStateVecsAlongTrajectory = new ArrayList<>(); for(int i = 0; i < kFZRef.kfStateVecsAlongTrajectory.size(); i++) { @@ -388,7 +391,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()); + sv.setPathLength(svc.getPathLength() + deltaPathToVtx); // Transition for the starting point from the final point at the last layer to vertex sv.setProjector(svc.getProjector()); sv.setProjectorDoca(svc.getProjectorDoca()); kfStateVecsAlongTrajectory.add(sv);