diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/RecUtilities.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/RecUtilities.java index 0204c84d20..9398d10b66 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/RecUtilities.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/RecUtilities.java @@ -495,7 +495,7 @@ public boolean reFitCircle(Seed seed, int iter, double xb, double yb) { if(fitStatus) { CircleFitPars pars = circlefit.getFit(); seed.getHelix().setCurvature(pars.rho()); - seed.getHelix().setDCA(-pars.doca()); + seed.getHelix().setDCA(pars.doca()); seed.getHelix().setPhiAtDCA(pars.phi()); seed.update_Crosses(xb,yb); } diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/TracksFromTargetRec.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/TracksFromTargetRec.java index d9ffecf6c5..b2854ae234 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/TracksFromTargetRec.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/TracksFromTargetRec.java @@ -132,7 +132,7 @@ public List getSeeds(List> clusters, List getSeeds(List> clusters, List(); for(int i =0; i() ); - for(int l =0; l<3; l++) { + for(int l =0; l<6; l++) { sortedCrosses.get(i).add(l,new ArrayList<>() ); } } @@ -100,70 +100,144 @@ public void findSeedCrossList(List crosses) { /* - Scans overphase space to find groups of BMT crosses + Scans overphase space to find groups of SVT and BMT-Z crosses */ private void findSeedCrossesFixedBin(List crosses, double phiShift) { - for(int b =0; b= NBINS) { + binIdx = NBINS - 1; } - int binIdx = (int) (phi / (360./NBINS) ); - - if(binIdx>NBINS-1) - binIdx = NBINS-1; - sortedCrosses.get(binIdx).get(crosses.get(i).getRegion() - 1).add(crosses.get(i)); - LPhi[binIdx][crosses.get(i).getRegion() - 1]++; + int layer = -1; + + if (cross.getDetector() == DetectorType.BST) { + + // BST region 1,2,3 -> layer 0,1,2 + layer = cross.getRegion() - 1; + + } else if (cross.getDetector() == DetectorType.BMT) { + + // BMT region 1,2,3 -> layer 3,4,5 + layer = cross.getRegion() + 2; + } + + // Ignore unexpected detector/region + if (layer < 0 || layer >= 6) { + continue; + } + + sortedCrosses.get(binIdx).get(layer).add(cross); + LPhi[binIdx][layer]++; } - + + // ------------------------------------------------------------ + // 2. For each phi bin, generate ALL combinations + // ------------------------------------------------------------ for (int b = 0; b < NBINS; b++) { - int max_layers =0; - for (int la = 0; la < 3; la++) { - if(LPhi[b][la]>0) - max_layers++; + + // Find occupied layers + List occupiedLayers = new ArrayList<>(); + + for (int la = 0; la < 6; la++) { + if (LPhi[b][la] > 0) { + occupiedLayers.add(la); + } + } + + // Need at least two different layers + if (occupiedLayers.size() < 2) { + continue; } - if (sortedCrosses.get(b) != null && max_layers >= 2) { - double SumLyr=0; - while(LPhi[b][0]+LPhi[b][1]+ LPhi[b][2]>=max_layers) { - if(SumLyr!=LPhi[b][0]+LPhi[b][1]+ LPhi[b][2]) { - SumLyr = LPhi[b][0]+LPhi[b][1]+ LPhi[b][2]; - } - ArrayList hits = new ArrayList<>(); - for (int la = 0; la < 3; la++) { - if (sortedCrosses.get(b).get(la) != null && LPhi[b][la]>0) { - if (sortedCrosses.get(b).get(la).get(LPhi[b][la]-1) != null - && sortedCrosses.get(b).get(la).size()>0) { - hits.add(sortedCrosses.get(b).get(la).get(LPhi[b][la]-1)); - - if(LPhi[b][la]>1) - LPhi[b][la]--; - if(SumLyr==max_layers) - LPhi[b][la]=0; - } - } - } - - if (hits.size() >= 2) { - double seedIdx=0; - int s = hits.size(); - int index = (int) Math.pow(2,s); - for(Cross c : hits) { - seedIdx +=c.getId()*Math.pow(10, index); - index-=4; - } - seedMap.put(seedIdx, hits); - } + + // Cross combos for any 2 more layers + for (int nLayers = 2; nLayers <= occupiedLayers.size(); nLayers++) { + generateLayerCombinations(b, occupiedLayers, 0, nLayers, new ArrayList<>()); + } + } + } + + /* + * Generate combinations of occupied layers. + * + */ + private void generateLayerCombinations(int binIdx, List occupiedLayers, int start, int nLayers, ArrayList selectedLayers) { + + // ------------------------------------------------------------ + // Select the required number of layers + // ------------------------------------------------------------ + if (selectedLayers.size() == nLayers) { + ArrayList currentCombination = new ArrayList<>(); + generateCrossCombinations(binIdx, selectedLayers, 0, currentCombination); + return; + } + + // ------------------------------------------------------------ + // Select next layer + // ------------------------------------------------------------ + for (int i = start; i < occupiedLayers.size(); i++) { + selectedLayers.add(occupiedLayers.get(i)); + generateLayerCombinations(binIdx, occupiedLayers, i + 1, nLayers, selectedLayers); + selectedLayers.remove(selectedLayers.size() - 1); + } + } + + /* + * Generate Cartesian products of Crosses from the selected layers. + * + */ + private void generateCrossCombinations(int binIdx, List selectedLayers, int layerIndex, ArrayList currentCombination) { + // ------------------------------------------------------------ + // All selected layers have been processed + // ------------------------------------------------------------ + if (layerIndex == selectedLayers.size()) { + if (currentCombination.size() >= 2) { + ArrayList crosses = new ArrayList<>(currentCombination); + long seedIdx = 0; + for (Cross c : crosses) { + seedIdx = seedIdx * 10000L + c.getId(); } + + seedMap.put((double) seedIdx, crosses); } + return; + } + + // Current detector layer + int layer = selectedLayers.get(layerIndex); + List layerCrosses = sortedCrosses.get(binIdx).get(layer); + + // ------------------------------------------------------------ + // Try EVERY Cross in this layer + // ------------------------------------------------------------ + for (Cross cross : layerCrosses) { + currentCombination.add(cross); + generateCrossCombinations(binIdx, selectedLayers, layerIndex + 1,currentCombination); + currentCombination.remove(currentCombination.size() - 1); } } @@ -237,15 +311,14 @@ private int countType(List cand, DetectorType dt) { } private double calcResi(double rho, double d0, double phi0, double xc, double yc) { double r = Math.sqrt(xc*xc+yc*yc); - double par = 1. - ((r * r - d0 * d0) * rho * rho) / (2. * (1. + d0 * Math.abs(rho))); + double par = (2. * (1. + d0 * rho) - (r * r - d0 * d0) * rho * rho) / (2. * Math.abs(1. + d0 * rho)); + double newPathLength = Math.abs(Math.acos(par) / rho); int charge = (int) Math.signum(rho); double alpha = -newPathLength * rho; - - double x = d0 * charge * Math.sin(phi0) + (charge / Math.abs(rho)) - * (Math.sin(phi0) - Math.cos(alpha) * Math.sin(phi0) - Math.sin(alpha) * Math.cos(phi0)); - double y = -d0 * charge * Math.cos(phi0) - (charge / Math.abs(rho)) - * (Math.cos(phi0) + Math.sin(alpha) * Math.sin(phi0) - Math.cos(alpha) * Math.cos(phi0)); + + double x = d0 * Math.sin(phi0) + 1./rho * (Math.sin(phi0) - Math.sin(phi0 + alpha)); + double y = -d0 * Math.cos(phi0) - 1./rho * (Math.cos(phi0) - Math.cos(phi0 + alpha)); double res2 = ((x-xc)*(x-xc)+(y-yc)*(y-yc)); return Math.sqrt(res2); diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/Helix.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/Helix.java index 107e5e7e2d..7cd7a61bdc 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/Helix.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/Helix.java @@ -178,15 +178,6 @@ public double radius() { return 1. / C; } - // (x,y) coordinates of the circle center - public double xCen() { - return (radius() - this.getDCA()) * Math.sin(this.getPhiAtDCA()); - } - - public double yCen() { - return (-radius() + this.getDCA()) * Math.cos(this.getPhiAtDCA()); - } - // (x,y) coordinates of the dca public double xDCA() { return -this.getDCA() * Math.sin(this.getPhiAtDCA()); @@ -198,7 +189,7 @@ public double yDCA() { public Point3D getVertex() { return new Point3D(this.xDCA()+this.getXb(),this.yDCA()+this.getYb(),this.getZ0()); - } + } public double getPt(double solenoidMag) { double pt = Constants.LIGHTVEL * this.radius() * solenoidMag; @@ -214,18 +205,7 @@ public Vector3D getPXYZ(double solenoidMag) { double py = pt*Math.sin(this.getPhiAtDCA()); return new Vector3D(px,py,pz); - } - - public double getArcLengthDCA(Point3D refpoint) { - //insure that the refpoint is on the helix - if (refpoint == null) { - return 0; - } - double refX = radius() * Math.cos(refpoint.toVector3D().phi()); - double refY = radius() * Math.sin(refpoint.toVector3D().phi()); - double arclen = arcLength(xCen(), yCen(), radius(), xCen(), yCen(), refX, refY); - return arclen; - } + } // this method finds the arclength between 2 points in a circle // this private method is used to get the pathlength from a point on the helical track to the distance of closest approach @@ -272,16 +252,13 @@ public Point3D getPointAtRadius(double r) { return new Point3D(x, y, z); } - - double par = 1. - ((r * r - d0 * d0) * omega * omega) / (2. * (1. + d0 * Math.abs(omega))); + double par = (2. * (1. + d0 * omega) - (r * r - d0 * d0) * omega * omega) / (2. * Math.abs(1. + d0 * omega)); double newPathLength = Math.abs(Math.acos(par) / omega); double alpha = -newPathLength * omega; - double x = d0 * charge * Math.sin(phi0) + (charge / Math.abs(omega)) - * (Math.sin(phi0) - Math.cos(alpha) * Math.sin(phi0) - Math.sin(alpha) * Math.cos(phi0))+xb; - double y = -d0 * charge * Math.cos(phi0) - (charge / Math.abs(omega)) - * (Math.cos(phi0) + Math.sin(alpha) * Math.sin(phi0) - Math.cos(alpha) * Math.cos(phi0))+yb; + double x = d0 * Math.sin(phi0) + 1./omega * (Math.sin(phi0) - Math.sin(phi0 + alpha)); + double y = -d0 * Math.cos(phi0) - 1./omega * (Math.cos(phi0) - Math.cos(phi0 + alpha)); double z = z0 + newPathLength * tandip; return new Point3D(x, y, z); @@ -301,7 +278,7 @@ public Vector3D getTrackDirectionAtRadius(double r) { return new Vector3D(ux, uy, uz); } - double par = 1. - ((r * r - d0 * d0) * omega * omega) / (2. * (1. + d0 * Math.abs(omega))); + double par = (2. * (1. + d0 * omega) - (r * r - d0 * d0) * omega * omega) / (2. * Math.abs(1. + d0 * omega)); double newPathLength = Math.abs(Math.acos(par) / omega); double alpha = newPathLength * omega;