diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/utilities/ProbChi2perNDF.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/utilities/ProbChi2perNDF.java new file mode 100644 index 0000000000..0a67deb03a --- /dev/null +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/utilities/ProbChi2perNDF.java @@ -0,0 +1,23 @@ +package org.jlab.clas.tracking.utilities; + +import org.apache.commons.math3.special.Gamma; + +/** Chi-square probability utilities */ +public class ProbChi2perNDF { + + /** + * Returns the upper-tail probability (survival function, or p-value) of the + * chi-square distribution: P(X >= chi2) for X ~ χ²(ndf). + * @param chi2 the chi-square statistic + * @param ndf the number of degrees of freedom + * @return the probability of observing a chi-square value at least as large as {@code chi2}, in the range [0, 1] + */ + public static double prob(double chi2, int ndf) { + if(ndf <= 0) + return 0.0; + if(chi2 <= 0.0) { // edge cases follow ROOT's `TMath::Prob` conventions + return chi2 < 0.0 ? 0.0 : 1.0; + } + return Gamma.regularizedGammaQ(0.5 * ndf, 0.5 * chi2); + } +} diff --git a/common-tools/coat-libs/pom.xml b/common-tools/coat-libs/pom.xml index 86ae62720b..136ce56b3c 100644 --- a/common-tools/coat-libs/pom.xml +++ b/common-tools/coat-libs/pom.xml @@ -36,12 +36,6 @@ j4ml-clas12 - - org.jlab.plugins - trackfitter - 1.0 - - org.jlab.clas clas-io diff --git a/reconstruction/cvt/pom.xml b/reconstruction/cvt/pom.xml index 0196338c09..4f49cd1f7e 100644 --- a/reconstruction/cvt/pom.xml +++ b/reconstruction/cvt/pom.xml @@ -85,11 +85,6 @@ test jar - - org.jlab.plugins - trackfitter - 1.0 - junit junit diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/Constants.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/Constants.java index 51303f0636..c3bacfd08e 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/Constants.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/Constants.java @@ -269,7 +269,7 @@ public static double getRCUT() { return RCUT; } - public static void setRCUT(double r) { + private static void setRCUT(double r) { RCUT = r; } @@ -281,7 +281,7 @@ public static double getZRANGE() { return ZRANGE; } - public static void setZRANGE(double zr) { + private static void setZRANGE(double zr) { ZRANGE = zr; } diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/fit/CosmicFitter.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/fit/CosmicFitter.java index e8032d2aa4..ab0cdd9bf5 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/fit/CosmicFitter.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/fit/CosmicFitter.java @@ -6,7 +6,7 @@ import org.jlab.geom.prim.Point3D; import org.jlab.geom.prim.Vector3D; import org.jlab.rec.cvt.trajectory.Ray; -import trackfitter.fitter.utilities.ProbChi2perNDF; +import org.jlab.clas.tracking.utilities.ProbChi2perNDF; /** * A fitter which does sequential fit (for x, y coordinates) and then (for r, z @@ -21,7 +21,6 @@ public class CosmicFitter { private CosmicFitPars _rayfitoutput; private LineFitPars _linefitparsRhoZ; - private LineFitPars _linefitparsYX; /** * Status of the HelicalTrackFit */ @@ -164,7 +163,7 @@ public FitStatus fit2(List X, List Y, List Z, List(X.size()), X.size()); + _linefitYX.fitStatus(Y, X, errRt, new ArrayList<>(X.size()), X.size()); _linefitRhoZ = new LineFitter(); boolean linefitstatusOK = _linefitRhoZ.fitStatus(Rho, Z, errRho, ErrZ, Z.size()); @@ -174,7 +173,6 @@ public FitStatus fit2(List X, List Y, List Z, List Constants.CHI2CUT * (this.getNDF() + 5)) return false; if(this.getNDF() < Constants.NDFCUT) return false; - if(this.getHelix().radius() < Constants.getInstance().getRCUT()) + if(this.getHelix().radius() < Constants.getRCUT()) return false; - if(Math.abs(Geometry.getInstance().getTargetZOffset()-this.getHelix().getZ0()) > Geometry.getInstance().getTargetHalfLength()+Constants.getInstance().getZRANGE()) + if(Math.abs(Geometry.getInstance().getTargetZOffset()-this.getHelix().getZ0()) > Geometry.getInstance().getTargetHalfLength()+Constants.getZRANGE()) return false; else return true; diff --git a/reconstruction/dc/pom.xml b/reconstruction/dc/pom.xml index 142352eb2f..cb0bbee9c9 100644 --- a/reconstruction/dc/pom.xml +++ b/reconstruction/dc/pom.xml @@ -97,11 +97,6 @@ snr 14.2.1-SNAPSHOT - - org.jlab.plugins - trackfitter - 1.0 - org.jlab.clas swim-tools diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/RecoBankWriter.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/RecoBankWriter.java index d018dd3c8d..a542d0dfed 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/RecoBankWriter.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/banks/RecoBankWriter.java @@ -13,7 +13,7 @@ import org.jlab.rec.dc.segment.Segment; import org.jlab.rec.dc.track.Track; -import trackfitter.fitter.utilities.*; +import org.jlab.clas.tracking.utilities.ProbChi2perNDF; /** * A class to fill the reconstructed DC banks @@ -36,8 +36,8 @@ public RecoBankWriter(Banks names) { public void updateListsWithClusterInfo(List fhits, List clusters) { - ArrayList rmHits = new ArrayList(); - ArrayList addHits = new ArrayList(); + List rmHits = new ArrayList<>(); + List addHits = new ArrayList<>(); for (int i = 0; i < clusters.size(); i++) { clusters.get(i).set_Id(i + 1); for (int j = 0; j < clusters.get(i).size(); j++) { @@ -204,8 +204,6 @@ public DataBank fillHBClustersBank(DataEvent event, List cluslist DataBank bank = event.createBank(name, cluslist.size()); int[] hitIdxArray = new int[12]; - if(cluslist==null) - return bank; for (int i = 0; i < cluslist.size(); i++) { if (cluslist.get(i) == null || cluslist.get(i).get_Id() == -1) { continue; diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/cross/CrossListFinder.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/cross/CrossListFinder.java index 5399b2dadb..968139284b 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/cross/CrossListFinder.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/cross/CrossListFinder.java @@ -15,7 +15,7 @@ import org.jlab.rec.dc.segment.Segment; import org.jlab.rec.dc.timetodistance.TimeToDistanceEstimator; import org.jlab.utils.groups.IndexedTable; -import trackfitter.fitter.LineFitter; +import org.jlab.rec.dc.track.fit.basefit.LineFitter; /** * A class with methods used to find lists of crosses. This is the Pattern Recognition step used in track seeding, to 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 547f629302..f90b26236a 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 @@ -8,7 +8,6 @@ import java.util.Map; import java.util.logging.Level; import java.util.logging.Logger; -import org.jlab.clas.clas.math.FastMath; import org.jlab.clas.swimtools.Swim; import org.jlab.detector.geant4.v2.DCGeant4Factory; import org.jlab.geom.prim.Line3D; @@ -24,8 +23,8 @@ import org.jlab.rec.dc.trajectory.Trajectory; import org.jlab.rec.dc.trajectory.TrajectoryFinder; -import trackfitter.fitter.LineFitPars; -import trackfitter.fitter.LineFitter; +import org.jlab.rec.dc.track.fit.basefit.LineFitter; +import org.jlab.rec.dc.track.fit.basefit.LineFitPars; import org.jlab.clas.tracking.kalmanfilter.Surface; import org.jlab.clas.tracking.kalmanfilter.zReference.KFitter; @@ -95,47 +94,47 @@ private boolean PassNSuperlayerTracking(List crossesInTrk, Track cand) { return pass; } - private double getHitBasedFitChi2ToCrosses(int sector, double x1, double y1, double z1, - double x2, double y2, double z2, double x3, - double y3, double z3, double p, int q, double x, - double y, double z, double tanThX, double tanThY, Swim dcSwim) { - double pz = p / Math.sqrt(tanThX * tanThX + tanThY * tanThY + 1); - - dcSwim.SetSwimParameters(x, y, z, - -pz * tanThX, -pz * tanThY, -pz, - -q); - double chi2 = 0; // assume err =1 on points - double[] R = dcSwim.SwimToPlaneTiltSecSys(sector, z3); - if (R == null) { - return Double.POSITIVE_INFINITY; - } - - //chi2 += (R[0] - x3) * (R[0] - x3) + (R[1] - y3) * (R[1] - y3); - chi2 += (Math.sqrt(R[0] * R[0] + R[1] * R[1]) - Math.sqrt(x3 * x3 + y3 * y3)) * (Math.sqrt(R[0] * R[0] + R[1] * R[1]) - Math.sqrt(x3 * x3 + y3 * y3)); - dcSwim.SetSwimParameters(R[0], R[1], R[2], - R[3], R[4], R[5], - -q); - R = dcSwim.SwimToPlaneTiltSecSys(sector, z2); - if (R == null) { - return Double.POSITIVE_INFINITY; - } - - dcSwim.SetSwimParameters(R[0], R[1], R[2], - R[3], R[4], R[5], - -q); - chi2 += (R[0] - x2) * (R[0] - x2) + (R[1] - y2) * (R[1] - y2); - dcSwim.SetSwimParameters(R[0], R[1], R[2], - R[3], R[4], R[5], - -q); - R = dcSwim.SwimToPlaneTiltSecSys(sector, z1); - if (R == null) { - return Double.POSITIVE_INFINITY; - } - - chi2 += (R[0] - x1) * (R[0] - x1) + (R[1] - y1) * (R[1] - y1); - - return chi2; - } + // private double getHitBasedFitChi2ToCrosses(int sector, double x1, double y1, double z1, + // double x2, double y2, double z2, double x3, + // double y3, double z3, double p, int q, double x, + // double y, double z, double tanThX, double tanThY, Swim dcSwim) { + // double pz = p / Math.sqrt(tanThX * tanThX + tanThY * tanThY + 1); + // + // dcSwim.SetSwimParameters(x, y, z, + // -pz * tanThX, -pz * tanThY, -pz, + // -q); + // double chi2 = 0; // assume err =1 on points + // double[] R = dcSwim.SwimToPlaneTiltSecSys(sector, z3); + // if (R == null) { + // return Double.POSITIVE_INFINITY; + // } + // + // //chi2 += (R[0] - x3) * (R[0] - x3) + (R[1] - y3) * (R[1] - y3); + // chi2 += (Math.sqrt(R[0] * R[0] + R[1] * R[1]) - Math.sqrt(x3 * x3 + y3 * y3)) * (Math.sqrt(R[0] * R[0] + R[1] * R[1]) - Math.sqrt(x3 * x3 + y3 * y3)); + // dcSwim.SetSwimParameters(R[0], R[1], R[2], + // R[3], R[4], R[5], + // -q); + // R = dcSwim.SwimToPlaneTiltSecSys(sector, z2); + // if (R == null) { + // return Double.POSITIVE_INFINITY; + // } + // + // dcSwim.SetSwimParameters(R[0], R[1], R[2], + // R[3], R[4], R[5], + // -q); + // chi2 += (R[0] - x2) * (R[0] - x2) + (R[1] - y2) * (R[1] - y2); + // dcSwim.SetSwimParameters(R[0], R[1], R[2], + // R[3], R[4], R[5], + // -q); + // R = dcSwim.SwimToPlaneTiltSecSys(sector, z1); + // if (R == null) { + // return Double.POSITIVE_INFINITY; + // } + // + // chi2 += (R[0] - x1) * (R[0] - x1) + (R[1] - y1) * (R[1] - y1); + // + // return chi2; + // } /* ** a: parameter a for fitting function f(x) = a*z^2 + b*z + c in xz plane @@ -208,7 +207,7 @@ private int calcInitTrkQ(double a, double TORSCALE) { } public List getTrackCandsAI(CrossList crossList, DCGeant4Factory DcDetector, Swim dcSwim, HBTrackStateEstimator hbTSEstimator){ - List cands = new ArrayList(); + List cands = new ArrayList<>(); for (List aCrossList : crossList) { Track cand = new Track(); @@ -752,24 +751,6 @@ else if(t1.get_FitChi2()/t1.get_FitNDF()==t2.get_FitChi2()/t2.get_FitNDF() && i> trkcands.addAll(selectedTracks); } - /** - * @param selectedTracks the list of selected tracks - * @param selectedTrk the selected track - * @return a boolean indicating if the track is in the list - */ - private boolean ListContainsTrack(List selectedTracks, Track selectedTrk) { - boolean isInList = false; - for (Track trk : selectedTracks) { - if (trk == null) { - continue; - } - if (trk.get_Id() == selectedTrk.get_Id()) { - isInList = true; - } - } - return isInList; - } - /** * @param track the track * @param trkcands the list of candidates @@ -880,31 +861,31 @@ public void setHitDoubletsInfo(Segment seg) { } } - private double calcCurvSign(Track cand) { - double P0x = cand.get(0).get_Point().z(); - double P1x = cand.get(1).get_Point().z(); - double P2x = cand.get(2).get_Point().z(); - double P0y = cand.get(0).get_Point().x(); - double P1y = cand.get(1).get_Point().x(); - double P2y = cand.get(2).get_Point().x(); - - if (Math.abs(P1x - P0x) < 1.0e-18 || Math.abs(P2x - P1x) < 1.0e-18) { - return 0.0; - } - - // Find the intersection of the lines joining the innermost to middle and middle to outermost point - double ma = (P1y - P0y) / (P1x - P0x); - double mb = (P2y - P1y) / (P2x - P1x); - - if (Math.abs(mb - ma) < 1.0e-18) { - return 0.0; - } - - double xcen = 0.5 * (ma * mb * (P0y - P2y) + mb * (P0x + P1x) - ma * (P1x + P2x)) / (mb - ma); - double ycen = (-1. / mb) * (xcen - 0.5 * (P1x + P2x)) + 0.5 * (P1y + P2y); - - return Math.signum(ycen); - } + // private double calcCurvSign(Track cand) { + // double P0x = cand.get(0).get_Point().z(); + // double P1x = cand.get(1).get_Point().z(); + // double P2x = cand.get(2).get_Point().z(); + // double P0y = cand.get(0).get_Point().x(); + // double P1y = cand.get(1).get_Point().x(); + // double P2y = cand.get(2).get_Point().x(); + // + // if (Math.abs(P1x - P0x) < 1.0e-18 || Math.abs(P2x - P1x) < 1.0e-18) { + // return 0.0; + // } + // + // // Find the intersection of the lines joining the innermost to middle and middle to outermost point + // double ma = (P1y - P0y) / (P1x - P0x); + // double mb = (P2y - P1y) / (P2x - P1x); + // + // if (Math.abs(mb - ma) < 1.0e-18) { + // return 0.0; + // } + // + // double xcen = 0.5 * (ma * mb * (P0y - P2y) + mb * (P0x + P1x) - ma * (P1x + P2x)) / (mb - ma); + // double ycen = (-1. / mb) * (xcen - 0.5 * (P1x + P2x)) + 0.5 * (P1y + P2y); + // + // return Math.signum(ycen); + // } private List findStraightTracks(CrossList crossList, DCGeant4Factory DcDetector, double TORSCALE, Swim dcSwim) { diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/basefit/LineFitPars.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/basefit/LineFitPars.java index 46d078eb01..d6414be29d 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/basefit/LineFitPars.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/basefit/LineFitPars.java @@ -1,6 +1,6 @@ package org.jlab.rec.dc.track.fit.basefit; -import trackfitter.fitter.utilities.ProbChi2perNDF; +import org.jlab.clas.tracking.utilities.ProbChi2perNDF; /** * The fit parameters of a line fit returned by LineFitter diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/basefit/LineFitter.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/basefit/LineFitter.java index 443db3eb7d..28f3aa6c40 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/basefit/LineFitter.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/track/fit/basefit/LineFitter.java @@ -1,6 +1,5 @@ package org.jlab.rec.dc.track.fit.basefit; -import java.util.ArrayList; import java.util.List; /** A least square fitting method @@ -15,9 +14,9 @@ public class LineFitter { // the constructor public LineFitter() { } - private final List w = new ArrayList(); + // fit status - public boolean fitStatus(List x, List y, List sigma_x, List sigma_y, int nbpoints) { + public boolean fitStatus(double[] x, double[] y, double[] sigma_x, double[] sigma_y, int nbpoints) { boolean fitStat = false; if (nbpoints>=2) { // must have enough points to do the fit @@ -25,22 +24,21 @@ public boolean fitStatus(List x, List y, List sigma_x, L // initialize weight-sum and moments double Sw, Sx, Sy, Sxx, Sxy; Sw = Sx = Sy = Sxx = Sxy =0.; + double[] w = new double[nbpoints]; - w.clear(); - ((ArrayList) w).ensureCapacity(nbpoints); for (int i = 0; i0 @@ -61,8 +59,8 @@ public boolean fitStatus(List x, List y, List sigma_x, L double chi_2 = 0.; double pointchi_2[] = new double[nbpoints]; //individual chi2 for each fitted point for (int j = 0; j x, List y, List sigma_x, L // if there is a fit return true return fitStat; } + + // fit status (overload) + public boolean fitStatus(List x, List y, List sigma_x, List sigma_y, int nbpoints) { + return fitStatus( + x.stream().mapToDouble(Double::doubleValue).toArray(), + y.stream().mapToDouble(Double::doubleValue).toArray(), + sigma_x.stream().mapToDouble(Double::doubleValue).toArray(), + sigma_y.stream().mapToDouble(Double::doubleValue).toArray(), + nbpoints); + } + // return the fit result public LineFitPars getFit() { return _linefitresult; diff --git a/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/TrajectoryFinder.java b/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/TrajectoryFinder.java index 9fa71474d4..c153425457 100644 --- a/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/TrajectoryFinder.java +++ b/reconstruction/dc/src/main/java/org/jlab/rec/dc/trajectory/TrajectoryFinder.java @@ -7,7 +7,7 @@ import org.jlab.geom.prim.Vector3D; import org.jlab.rec.dc.Constants; import org.jlab.rec.dc.cross.Cross; -import trackfitter.fitter.LineFitter; +import org.jlab.rec.dc.track.fit.basefit.LineFitter; import Jama.Matrix; /** @@ -175,7 +175,7 @@ public void fitTrajectory(List candCrossList) { } lineFit = new LineFitter(); - boolean linefitstatusOK = lineFit.fitStatus(z, theta_x, new double[3], theta_x_err, 3); + lineFit.fitStatus(z, theta_x, new double[3], theta_x_err, 3); TrajChisqProbFitXZ = lineFit.getFit().getProb(); x_fitCoeff = quadraticLRFit(z, x, x_err);