From 56cd80966b8799442448afbf3a580055ae32544e Mon Sep 17 00:00:00 2001 From: ziegler Date: Thu, 12 Dec 2024 10:41:56 -0500 Subject: [PATCH 01/12] Bug fix in setting the measurements for cosmics. Added missing hemisphere for BST clusters. --- .../rec/cvt/measurement/Measurements.java | 39 ++++++++++++------- 1 file changed, 25 insertions(+), 14 deletions(-) diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java index b00ca82c4c..810ac76713 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java @@ -283,11 +283,14 @@ private void addClusters(DetectorType type, List clusters) { // } private List getClusterSurfaces(DetectorType type, List clusters, Ray ray) { - - List surfaces = this.getClusterSurfaces(type, clusters, ray); - for(Surface surf : surfaces) { - surf.hemisphere = this.getHemisphere(ray, surf); - } + List surfaces = new ArrayList<>(); + for(Cluster cluster : clusters) { + Surface surf = this.getClusterSurface(type, cluster); + if(surf!=null) { + surf.hemisphere = this.getHemisphere(ray, surf); + surfaces.add(surf); + } + } return surfaces; } @@ -341,7 +344,8 @@ private void addMissing(StraightTrack ray) { if(surface == null) continue; surface.hemisphere=hemisphere; surface.passive=true; - if(debug) System.out.println("Generating surface for missing index " + i + " detector " + type.getName() + " layer " + layer + " sector " + surface.getSector()); + if(debug) System.out.println("Generating surface for missing index " + i +" id "+id+ " detector " + type.getName() + " layer " + layer + " sector " + surface.getSector()+" hemisphere "+ + hemisphere); this.add(i, surface); } } @@ -434,15 +438,22 @@ private boolean isCrossed(Ray ray, Surface surface){ } private int getHemisphere(Ray ray, Surface surface){ - if(surface.cylinder==null) - return 0; - List trajs = new ArrayList<>(); - Line3D line = ray.toLine(); - if(surface.cylinder.intersection(line, trajs)>= 1) { - return (int) Math.signum(trajs.get(0).y()); + int h =0; + if(surface.cylinder!=null) { + List trajs = new ArrayList<>(); + Line3D line = ray.toLine(); + if(surface.cylinder.intersection(line, trajs)>= 1) { + h= (int) Math.signum(trajs.get(0).y()); + } } - else - return 0; + if(surface.plane!=null) { + Point3D traj = new Point3D(); + Line3D line = ray.toLine(); + if(surface.plane.intersection(line, traj)>= 1) { + h= (int) Math.signum(traj.y()); + } + } + return h; } private int getHemisphere(Cluster cluster, Seed seed, Surface surface) { if(surface.cylinder==null) From f1f1f029cdc70c74cb6094194090fca05c2bcd54 Mon Sep 17 00:00:00 2001 From: ziegler Date: Fri, 13 Dec 2024 18:03:05 -0500 Subject: [PATCH 02/12] fix in propagation to surface for shifted cylinders. --- .../kalmanfilter/helical/StateVecs.java | 44 ++++++++++++------- .../kalmanfilter/straight/StateVecs.java | 42 ++++++++++-------- 2 files changed, 51 insertions(+), 35 deletions(-) diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java index 067bcf5c5a..7c4d7ee661 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java @@ -1,5 +1,7 @@ package org.jlab.clas.tracking.kalmanfilter.helical; +import java.util.ArrayList; +import java.util.List; import org.jlab.clas.swimtools.Swim; import org.jlab.clas.tracking.kalmanfilter.AMeasVecs; import org.jlab.clas.tracking.kalmanfilter.AMeasVecs.MeasVec; @@ -42,9 +44,8 @@ public boolean getStateVecPosAtMeasSite(StateVec sv, AMeasVecs.MeasVec mv, Swim Point3D st = new Point3D(sv.x, sv.y, sv.z); Vector3D stu = new Vector3D(sv.px,sv.py,sv.pz).asUnit(); - + Line3D toPln = new Line3D(st, stu); if(mv.surface.plane!=null) { - Line3D toPln = new Line3D(st, stu); Point3D inters = new Point3D(); int ints = mv.surface.plane.intersection(toPln, inters); sv.x = inters.x(); @@ -53,21 +54,30 @@ public boolean getStateVecPosAtMeasSite(StateVec sv, AMeasVecs.MeasVec mv, Swim sv.path = inters.distance(st); } else if(mv.surface.cylinder!=null) { - mv.surface.toLocal().apply(st); - mv.surface.toLocal().apply(stu); - double r = mv.surface.cylinder.baseArc().radius(); - double delta = Math.sqrt((st.x()*stu.x()+st.y()*stu.y())*(st.x()*stu.x()+st.y()*stu.y())-(-r*r+st.x()*st.x()+st.y()*st.y())*(stu.x()*stu.x()+stu.y()*stu.y())); - double l = (-(st.x()*stu.x()+st.y()*stu.y())+delta)/(stu.x()*stu.x()+stu.y()*stu.y()); - if(Math.signum(st.y()+l*stu.y())!=mv.hemisphere) { - l = (-(st.x()*stu.x()+st.y()*stu.y())-delta)/(stu.x()*stu.x()+stu.y()*stu.y()); - } - Point3D inters = new Point3D(st.x()+l*stu.x(),st.y()+l*stu.y(),st.z()+l*stu.z()); - mv.surface.toGlobal().apply(inters); - // RDV: should switch to use clas-geometry intersection method, not done now to alwys return a value - sv.x = inters.x(); - sv.y = inters.y(); - sv.z = inters.z(); - sv.path = inters.distance(st); +// mv.surface.toLocal().apply(st); +// mv.surface.toLocal().apply(stu); +// double r = mv.surface.cylinder.baseArc().radius(); +// double delta = Math.sqrt((st.x()*stu.x()+st.y()*stu.y())*(st.x()*stu.x()+st.y()*stu.y())-(-r*r+st.x()*st.x()+st.y()*st.y())*(stu.x()*stu.x()+stu.y()*stu.y())); +// double l = (-(st.x()*stu.x()+st.y()*stu.y())+delta)/(stu.x()*stu.x()+stu.y()*stu.y()); +// if(Math.signum(st.y()+l*stu.y())!=mv.hemisphere) { +// l = (-(st.x()*stu.x()+st.y()*stu.y())-delta)/(stu.x()*stu.x()+stu.y()*stu.y()); +// } +// Point3D inters = new Point3D(st.x()+l*stu.x(),st.y()+l*stu.y(),st.z()+l*stu.z()); +// mv.surface.toGlobal().apply(inters); +// // RDV: should switch to use clas-geometry intersection method, not done now to alwys return a value +// sv.x = inters.x(); +// sv.y = inters.y(); +// sv.z = inters.z(); +// sv.path = inters.distance(st); + List inters = new ArrayList(); + int ints = mv.surface.cylinder.intersection(toPln, inters); + if(ints<1) return false; + + sv.x = inters.get(0).x() ; + sv.y = inters.get(0).y() ; + sv.z = inters.get(0).z() ; + sv.path = inters.get(0).distance(st); + } } else { if(swim==null) { // applicable only to planes parallel to the z -axis diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java index d811b23f9b..08afe43188 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java @@ -1,5 +1,7 @@ package org.jlab.clas.tracking.kalmanfilter.straight; +import java.util.ArrayList; +import java.util.List; import org.jlab.clas.swimtools.Swim; import org.jlab.clas.tracking.kalmanfilter.AMeasVecs; import org.jlab.clas.tracking.kalmanfilter.AMeasVecs.MeasVec; @@ -25,7 +27,7 @@ public boolean getStateVecPosAtMeasSite(StateVec vec, MeasVec mv, Swim swim) { Point3D ref = new Point3D(vec.x0,vec.y0,vec.z0); Vector3D u = new Vector3D(vec.tx, 1, vec.tz).asUnit(); - + Line3D toPln = new Line3D(ref,u); if(mv.k==0) { vec.x = vec.x0-vec.y0*vec.tx; vec.y = 0; @@ -34,7 +36,6 @@ public boolean getStateVecPosAtMeasSite(StateVec vec, MeasVec mv, Swim swim) { } else if(mv.hemisphere!=0) { if(mv.surface.plane!=null) { - Line3D toPln = new Line3D(ref,u); Point3D inters = new Point3D(); int ints = mv.surface.plane.intersection(toPln, inters); vec.x = inters.x() ; @@ -43,22 +44,27 @@ else if(mv.hemisphere!=0) { vec.path = inters.distance(ref); } if(mv.surface.cylinder!=null) { - mv.surface.toLocal().apply(ref); - mv.surface.toLocal().apply(u); - double r = 0.5*(mv.surface.cylinder.baseArc().radius()+mv.surface.cylinder.highArc().radius()); - double delta = Math.sqrt((ref.x()*u.x()+ref.y()*u.y())*(ref.x()*u.x()+ref.y()*u.y()) - -(-r*r+ref.x()*ref.x()+ref.y()*ref.y())*(u.x()*u.x()+u.y()*u.y())); - double l = (-(ref.x()*u.x()+ref.y()*u.y())+delta)/(u.x()*u.x()+u.y()*u.y()); - if(Math.signum(ref.y()+l*u.y())!=mv.hemisphere) { - l = (-(ref.x()*u.x()+ref.y()*u.y())-delta)/(u.x()*u.x()+u.y()*u.y()); - } - - Point3D cylInt = new Point3D(ref.x()+l*u.x(),ref.y()+l*u.y(),ref.z()+l*u.z()); - mv.surface.toGlobal().apply(cylInt); - vec.x = cylInt.x(); - vec.y = cylInt.y(); - vec.z = cylInt.z(); - vec.path = cylInt.distance(ref); + //mv.surface.toLocal().apply(ref); + //mv.surface.toLocal().apply(u); +// double r = 0.5*(mv.surface.cylinder.baseArc().radius()+mv.surface.cylinder.highArc().radius()); +// double delta = Math.sqrt((ref.x()*u.x()+ref.y()*u.y())*(ref.x()*u.x()+ref.y()*u.y()) +// -(-r*r+ref.x()*ref.x()+ref.y()*ref.y())*(u.x()*u.x()+u.y()*u.y())); +// double l = (-(ref.x()*u.x()+ref.y()*u.y())+delta)/(u.x()*u.x()+u.y()*u.y()); +// if(Math.signum(ref.y()+l*u.y())!=mv.hemisphere) { +// l = (-(ref.x()*u.x()+ref.y()*u.y())-delta)/(u.x()*u.x()+u.y()*u.y()); +// } +// +// Point3D cylInt = new Point3D(ref.x()+l*u.x(),ref.y()+l*u.y(),ref.z()+l*u.z()); +// mv.surface.toGlobal().apply(cylInt); + List inters = new ArrayList(); + int ints = mv.surface.cylinder.intersection(toPln, inters); + if(ints<1) return false; + + vec.x = inters.get(0).x() ; + vec.y = inters.get(0).y() ; + vec.z = inters.get(0).z() ; + vec.path = inters.get(0).distance(ref); + } return true; } From 6bf24ee299b21bf205273f07b9b7c78a1b0c8498 Mon Sep 17 00:00:00 2001 From: ziegler Date: Fri, 13 Dec 2024 18:08:29 -0500 Subject: [PATCH 03/12] Fix measurements sorting by Y intersection of the cosmic ray with the surface. --- .../rec/cvt/measurement/Measurements.java | 56 +++++++++++++++---- 1 file changed, 45 insertions(+), 11 deletions(-) diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java index 810ac76713..166dac3f5c 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java @@ -1,7 +1,9 @@ package org.jlab.rec.cvt.measurement; import java.util.ArrayList; +import java.util.LinkedHashMap; import java.util.List; +import java.util.Map; import org.jlab.clas.tracking.kalmanfilter.Material; import org.jlab.clas.tracking.kalmanfilter.Surface; import org.jlab.clas.tracking.objects.Strip; @@ -79,7 +81,7 @@ private void initCosmicSurfaces() { private void add(int index, Surface surface) { if(!(0<=index && index getMeasurements(StraightTrack cosmic) { this.addClusters(cosmic); this.addMissing(cosmic); List surfaces = new ArrayList<>(); + Map yMp = new LinkedHashMap<>(); for(Surface surf : cvtSurfaces) { if(surf!=null) { - if(debug) System.out.println(surf.toString()); if(surf.passive && surf.getIndex()!=0 && !this.isCrossed(cosmic.getRay(), surf)) { - if(debug) System.out.println("Removing surface " + surf.passive + " " + this.isCrossed(cosmic.getRay(), surf)); + if(debug) System.out.println("Removing surface "+surf.toString()+" " + surf.passive + " " + this.isCrossed(cosmic.getRay(), surf)); continue; } - surfaces.add(surf); + double sy = this.getY(cosmic.getRay(), surf); + yMp.put(sy, surf); } } - return surfaces; + Map sortedMap = new LinkedHashMap<>(); + yMp.entrySet() + .stream() + .sorted((entry1, entry2) -> entry2.getKey().compareTo(entry1.getKey())) // Sort in descending order + .forEach(entry -> sortedMap.put(entry.getKey(), entry.getValue())); // Insert sorted entries into new map + + if(debug) sortedMap.forEach((key, value) -> System.out.println(key + ": " + value)); + int i=0; + for(Surface surf : sortedMap.values()) { + surf.setIndex(i++); + surfaces.add(surf); + } + if(debug) sortedMap.forEach((key, value) -> System.out.println(key + ":: " + value)); + return surfaces; } public List getActiveMeasurements(StraightTrack cosmic) { @@ -429,12 +445,33 @@ else if(type==DetectorType.BMT) { } + private double getY(Ray ray, Surface surface) { + double y =0; + if(surface.cylinder!=null) { + List trajs = new ArrayList<>(); + Line3D line = ray.toLine(); + if(surface.cylinder.intersection(line, trajs)>= 1) { + y= trajs.get(0).y(); + } + } + + if(surface.plane!=null) { + Point3D traj = new Point3D(); + Line3D line = ray.toLine(); + if(surface.plane.intersection(line, traj)>= 1) { + y= traj.y(); + } + } + return y; + } + private boolean isCrossed(Ray ray, Surface surface){ if(surface.cylinder==null) return true; List trajs = new ArrayList<>(); Line3D line = ray.toLine(); - return surface.cylinder.intersection(line, trajs) > 1; + return (surface.cylinder.intersection(line, trajs) >= 1 && + surface.hemisphere==(int) Math.signum(trajs.get(0).y())); } private int getHemisphere(Ray ray, Surface surface){ @@ -446,12 +483,9 @@ private int getHemisphere(Ray ray, Surface surface){ h= (int) Math.signum(trajs.get(0).y()); } } + if(surface.plane!=null) { - Point3D traj = new Point3D(); - Line3D line = ray.toLine(); - if(surface.plane.intersection(line, traj)>= 1) { - h= (int) Math.signum(traj.y()); - } + h=(int) surface.hemisphere; } return h; } From a1d3d3755ec4bc5f0d92f650e9671cc9453bbb6e Mon Sep 17 00:00:00 2001 From: ziegler Date: Sat, 4 Jan 2025 09:48:06 -0500 Subject: [PATCH 04/12] removed unused import --- .../tracking/kalmanfilter/helical/StateVecs.java | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java index 7c4d7ee661..14eeb15ce2 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java @@ -9,7 +9,6 @@ import org.jlab.clas.tracking.kalmanfilter.Surface; import org.jlab.clas.tracking.kalmanfilter.Units; import org.jlab.clas.tracking.trackrep.Helix; -import org.jlab.clas.tracking.utilities.MatrixOps; import org.jlab.geom.prim.Line3D; import org.jlab.geom.prim.Point3D; import org.jlab.geom.prim.Vector3D; @@ -72,13 +71,13 @@ else if(mv.surface.cylinder!=null) { List inters = new ArrayList(); int ints = mv.surface.cylinder.intersection(toPln, inters); if(ints<1) return false; + + sv.x = inters.get(0).x() ; + sv.y = inters.get(0).y() ; + sv.z = inters.get(0).z() ; + sv.path = inters.get(0).distance(st); - sv.x = inters.get(0).x() ; - sv.y = inters.get(0).y() ; - sv.z = inters.get(0).z() ; - sv.path = inters.get(0).distance(st); - - } + } } else { if(swim==null) { // applicable only to planes parallel to the z -axis Helix helix = sv.getHelix(xref, yref); From 6aaf945c375fa02631530b5e9fc260f98495bab5 Mon Sep 17 00:00:00 2001 From: ziegler Date: Sat, 4 Jan 2025 09:49:07 -0500 Subject: [PATCH 05/12] Initialized chi^2 --- .../main/java/org/jlab/clas/tracking/kalmanfilter/AKFitter.java | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/AKFitter.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/AKFitter.java index a940a6f9fe..ca09927de3 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/AKFitter.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/AKFitter.java @@ -24,7 +24,7 @@ public abstract class AKFitter { // return variables public boolean setFitFailed = false; - public double chi2; + public double chi2 = Double.POSITIVE_INFINITY; public int NDF; public int NDF0; public int numIter; From 861caf32b7f2d7cf3a888002dbfd96af9a4932ce Mon Sep 17 00:00:00 2001 From: ziegler Date: Sat, 4 Jan 2025 09:55:26 -0500 Subject: [PATCH 06/12] Added y intersect for rays. This is used to sort the surfaces from top to bottom. In the case of a shallow angle ray there can be 2 intersection with a passive surface in the same hemisphere, hence the need to sort by y. --- .../clas/tracking/kalmanfilter/Surface.java | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/Surface.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/Surface.java index 59348bd03d..75bc6fbced 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/Surface.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/Surface.java @@ -42,6 +42,7 @@ public class Surface implements Comparable { public double swimAccuracy; public boolean passive = false; public double hemisphere = 1; + public double rayYinterc = 9999; public double[] unc = new double[2]; public double[] doca = new double[2]; @@ -421,11 +422,18 @@ public void setTransformation(Transformation3D transform) { @Override public int compareTo(Surface o) { - if (this.index > o.index) { - return 1; + if(this.rayYinterc==9999) { + if (this.index > o.index) { + return 1; + } else { + return -1; + } } else { - return -1; - } + if (this.rayYinterc < o.rayYinterc) { + return 1; + } else { + return -1; + } + } } - } From 4fb1e223006274c20cf8d0704c803bdd9e2320d0 Mon Sep 17 00:00:00 2001 From: ziegler Date: Sat, 4 Jan 2025 09:57:47 -0500 Subject: [PATCH 07/12] For cosmic fitting, save the output for the best chi^2; do not kill the track if the next iteration chi^2 is worse. --- .../clas/tracking/kalmanfilter/straight/KFitter.java | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/KFitter.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/KFitter.java index 4c212a492d..003a5008aa 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/KFitter.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/KFitter.java @@ -54,22 +54,24 @@ public void runFitter(AStateVecs sv, AMeasVecs mv) { double newchisq = this.calc_chi2(sv, mv); // if curvature is 0, fit failed if(Double.isNaN(newchisq) || sv.smoothed().get(0)==null) { - this.setFitFailed = true; + this.setFitFailed = true; break; } else if(newchisq < this.chi2) { this.chi2 = newchisq; + finalStateVec = sv.new StateVec(sv.smoothed().get(0)); this.setTrajectory(sv, mv); - setFitFailed = false; + //System.out.println(finalStateVec.toString()); + //setFitFailed = false; } // stop if chi2 got worse else { break; } } - if(!this.setFitFailed) { - finalStateVec = sv.new StateVec(sv.smoothed().get(0)); - } + //if(!this.setFitFailed) { + // finalStateVec = sv.new StateVec(sv.smoothed().get(0)); + //} } @Override From 4528930e8b39123bad375994791e55b0743df9e0 Mon Sep 17 00:00:00 2001 From: ziegler Date: Sat, 4 Jan 2025 10:00:07 -0500 Subject: [PATCH 08/12] Use global fit y-intersection with the surface to assign fitted track intersection point when there are 2 intersection points of the ray with the cylinder. --- .../kalmanfilter/straight/StateVecs.java | 27 ++++++++++++++----- 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java index 08afe43188..fd493879aa 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java @@ -58,12 +58,26 @@ else if(mv.hemisphere!=0) { // mv.surface.toGlobal().apply(cylInt); List inters = new ArrayList(); int ints = mv.surface.cylinder.intersection(toPln, inters); - if(ints<1) return false; - - vec.x = inters.get(0).x() ; - vec.y = inters.get(0).y() ; - vec.z = inters.get(0).z() ; - vec.path = inters.get(0).distance(ref); + if(ints<1) { + if(!mv.surface.passive) { + return false; + } else { + mv.skip=true; + return true; + } + } else { + int index =0; + if(ints>1) { + double dh0 = Math.abs(inters.get(0).y()-mv.surface.rayYinterc); + double dh1 = Math.abs(inters.get(1).y()-mv.surface.rayYinterc); + if(dh1 Date: Sat, 4 Jan 2025 10:03:38 -0500 Subject: [PATCH 09/12] Assign correct hemisphere for cosmics. Use cluster information to determine y-intersection with the surface and hemisphere for active surfaces. Compute the y-intersection with the surface to be used for sorting in the kalman fitter. --- .../rec/cvt/measurement/Measurements.java | 237 ++++++++++++------ 1 file changed, 167 insertions(+), 70 deletions(-) diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java index 166dac3f5c..461152bcc9 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java @@ -1,7 +1,7 @@ package org.jlab.rec.cvt.measurement; import java.util.ArrayList; -import java.util.LinkedHashMap; +import java.util.HashMap; import java.util.List; import java.util.Map; import org.jlab.clas.tracking.kalmanfilter.Material; @@ -64,18 +64,18 @@ private void initTargetSurfaces(double xbeam, double ybeam, boolean beamSpot) { private void initCosmicSurfaces() { cvtSurfaces = new Surface[NSURFACES*2+3]; this.add(MLayer.COSMICPLANE.getIndex(1), this.getCosmicPlane()); - this.add(MLayer.SCHAMBER.getIndex(1), Geometry.getInstance().getScatteringChamber()); - this.add(MLayer.SHIELD.getIndex(1), Geometry.getInstance().getTargetShield()); - this.add(MLayer.INNERSVTCAGE.getIndex(1), Geometry.getInstance().getSVT().getFaradayCageSurfaces(0)); - this.add(MLayer.OUTERSVTCAGE.getIndex(1), Geometry.getInstance().getSVT().getFaradayCageSurfaces(1)); - this.add(MLayer.BMTINNERTUBE.getIndex(1), Geometry.getInstance().getBMT().getInnerTube()); - this.add(MLayer.BMTOUTERTUBE.getIndex(1), Geometry.getInstance().getBMT().getOuterTube()); - this.add(MLayer.SCHAMBER.getIndex(-1), Geometry.getInstance().getScatteringChamber()); + this.add(MLayer.SCHAMBER.getIndex(1), Geometry.getInstance().getScatteringChamber(),1); + this.add(MLayer.SHIELD.getIndex(1), Geometry.getInstance().getTargetShield(),1); + this.add(MLayer.INNERSVTCAGE.getIndex(1), Geometry.getInstance().getSVT().getFaradayCageSurfaces(0),1); + this.add(MLayer.OUTERSVTCAGE.getIndex(1), Geometry.getInstance().getSVT().getFaradayCageSurfaces(1),1); + this.add(MLayer.BMTINNERTUBE.getIndex(1), Geometry.getInstance().getBMT().getInnerTube(),1); + this.add(MLayer.BMTOUTERTUBE.getIndex(1), Geometry.getInstance().getBMT().getOuterTube(),1); + this.add(MLayer.SCHAMBER.getIndex(-1), Geometry.getInstance().getScatteringChamber(),-1); this.add(MLayer.SHIELD.getIndex(-1), Geometry.getInstance().getTargetShield(), -1); this.add(MLayer.INNERSVTCAGE.getIndex(-1), Geometry.getInstance().getSVT().getFaradayCageSurfaces(0), -1); this.add(MLayer.OUTERSVTCAGE.getIndex(-1), Geometry.getInstance().getSVT().getFaradayCageSurfaces(1), -1); - this.add(MLayer.BMTINNERTUBE.getIndex(-1), Geometry.getInstance().getBMT().getInnerTube()); - this.add(MLayer.BMTOUTERTUBE.getIndex(-1), Geometry.getInstance().getBMT().getOuterTube()); + this.add(MLayer.BMTINNERTUBE.getIndex(-1), Geometry.getInstance().getBMT().getInnerTube(),-1); + this.add(MLayer.BMTOUTERTUBE.getIndex(-1), Geometry.getInstance().getBMT().getOuterTube(),-1); } private void add(int index, Surface surface) { @@ -183,6 +183,7 @@ private Surface getCosmicPlane() { cosmic.addMaterial(Geometry.VACUUM); cosmic.setError(1); cosmic.hemisphere = 1; + cosmic.rayYinterc=0; cosmic.passive = true; return cosmic; } @@ -194,7 +195,7 @@ public List getMeasurements(Seed seed) { List surfaces = new ArrayList<>(); for(Surface surf : cvtSurfaces) { if(surf!=null) { - if(debug) System.out.println(surf.toString()); + if(debug) System.out.println("USING SURFACE "+surf.toString()); surfaces.add(surf); } } @@ -207,7 +208,7 @@ public List getActiveMeasurements(Seed seed) { for(Surface surf : surfaces) { if(surf.passive && surf.getIndex()!=0) continue; active.add(surf); - if(debug) System.out.println(surf.toString()); + if(debug) System.out.println("ACTIVE SURFACE "+surf.toString()); } return active; } @@ -231,31 +232,29 @@ public List getMeasurements(StraightTrack cosmic) { this.reset(); this.addClusters(cosmic); this.addMissing(cosmic); - List surfaces = new ArrayList<>(); - Map yMp = new LinkedHashMap<>(); + List surfaces ; + Map sMap = new HashMap<>(); for(Surface surf : cvtSurfaces) { if(surf!=null) { - if(surf.passive && surf.getIndex()!=0 && !this.isCrossed(cosmic.getRay(), surf)) { - if(debug) System.out.println("Removing surface "+surf.toString()+" " + surf.passive + " " + this.isCrossed(cosmic.getRay(), surf)); - continue; + + if(surf.passive && surf.getIndex()!=0) { + surf.hemisphere = this.getHemisphere((int)surf.hemisphere, cosmic.getRay(), surf); + + if(!this.isCrossed(cosmic.getRay(), surf)) { + if(debug) System.out.println("Removing surface "+surf.toString()+" " + surf.passive + " " + this.isCrossed(cosmic.getRay(), surf)); + continue; + } } - double sy = this.getY(cosmic.getRay(), surf); - yMp.put(sy, surf); + + sMap.put(surf.getIndex(), surf); + //surfaces.add(surf); } } - Map sortedMap = new LinkedHashMap<>(); - yMp.entrySet() - .stream() - .sorted((entry1, entry2) -> entry2.getKey().compareTo(entry1.getKey())) // Sort in descending order - .forEach(entry -> sortedMap.put(entry.getKey(), entry.getValue())); // Insert sorted entries into new map - - if(debug) sortedMap.forEach((key, value) -> System.out.println(key + ": " + value)); - int i=0; - for(Surface surf : sortedMap.values()) { - surf.setIndex(i++); - surfaces.add(surf); - } - if(debug) sortedMap.forEach((key, value) -> System.out.println(key + ":: " + value)); + surfaces = new ArrayList<>(sMap.values()); + if(debug) + for(Surface surf : surfaces) { + System.out.println("USED "+surf.toString()); + } return surfaces; } @@ -265,7 +264,6 @@ public List getActiveMeasurements(StraightTrack cosmic) { for(Surface surf : surfaces) { if(surf.passive && surf.getIndex()!=0) continue; active.add(surf); - if(debug) System.out.println(surf.toString()); } return active; } @@ -303,7 +301,7 @@ private List getClusterSurfaces(DetectorType type, List cluste for(Cluster cluster : clusters) { Surface surf = this.getClusterSurface(type, cluster); if(surf!=null) { - surf.hemisphere = this.getHemisphere(ray, surf); + surf.hemisphere = this.getHemisphere(cluster, ray, surf); surfaces.add(surf); } } @@ -358,10 +356,10 @@ private void addMissing(StraightTrack ray) { DetectorType type = MLayer.getDetectorType(id); Surface surface = this.getDetectorSurface(ray, type, layer, hemisphere); if(surface == null) continue; - surface.hemisphere=hemisphere; + surface.hemisphere = this.getHemisphere(hemisphere, ray.getRay(), surface); //resets hemisphere for shallow angle tracks surface.passive=true; if(debug) System.out.println("Generating surface for missing index " + i +" id "+id+ " detector " + type.getName() + " layer " + layer + " sector " + surface.getSector()+" hemisphere "+ - hemisphere); + surface.hemisphere+" y "+surface.rayYinterc); this.add(i, surface); } } @@ -445,50 +443,149 @@ else if(type==DetectorType.BMT) { } - private double getY(Ray ray, Surface surface) { - double y =0; - if(surface.cylinder!=null) { - List trajs = new ArrayList<>(); - Line3D line = ray.toLine(); - if(surface.cylinder.intersection(line, trajs)>= 1) { - y= trajs.get(0).y(); - } + + + private double getY(Cluster cluster, Ray ray, Surface surface) { + if (ray == null || surface == null) { + return Double.POSITIVE_INFINITY; } - - if(surface.plane!=null) { - Point3D traj = new Point3D(); - Line3D line = ray.toLine(); - if(surface.plane.intersection(line, traj)>= 1) { - y= traj.y(); - } + + Line3D line = ray.toLine(); // Shared line for both cylinder and plane + if (surface.cylinder != null) { + return handleCylinderIntersection(cluster, surface.cylinder, line); + } + + if (surface.plane != null) { + return handlePlaneIntersection(line, surface.plane); + } + + return Double.POSITIVE_INFINITY; // Default if no surface type matched + } + + private double handleCylinderIntersection(Cluster cluster, Cylindrical3D cylinder, Line3D line) { + if (cluster.getType() != BMTType.C) { + return cluster.center().y(); // For non-C type clusters, return the center y + } + + List intersections = new ArrayList<>(); + int intersectionCount = cylinder.intersection(line, intersections); + + switch (intersectionCount) { + case 0: + return Double.POSITIVE_INFINITY; + case 1: + return intersections.get(0).y(); + case 2: + // Choose the intersection closest to the cluster center on the z-axis + Point3D closest = getClosestIntersectionToClusterZ(cluster, intersections); + return closest.y(); + default: + return Double.POSITIVE_INFINITY; + } + } + + private double getY(int h, Ray ray, Surface surface) { + if (ray == null || surface == null) { + return Double.POSITIVE_INFINITY; + } + + Line3D line = ray.toLine(); // Shared line for both cylinder and plane + if (surface.cylinder != null) { + return handleCylinderIntersection(h, surface.cylinder, line); + } + + if (surface.plane != null) { + return handlePlaneIntersection(line, surface.plane); } - return y; + + return Double.POSITIVE_INFINITY; // Default if no surface type matched } - private boolean isCrossed(Ray ray, Surface surface){ - if(surface.cylinder==null) - return true; - List trajs = new ArrayList<>(); - Line3D line = ray.toLine(); - return (surface.cylinder.intersection(line, trajs) >= 1 && - surface.hemisphere==(int) Math.signum(trajs.get(0).y())); + private double handleCylinderIntersection(int h, Cylindrical3D cylinder, Line3D line) { + + List intersections = new ArrayList<>(); + int intersectionCount = cylinder.intersection(line, intersections); + + switch (intersectionCount) { + case 0: + return Double.POSITIVE_INFINITY; + case 1: + return intersections.get(0).y(); + case 2: + // Choose + double yFirst = intersections.get(0).y(); + double ySecond =intersections.get(1).y(); + + return (int) Math.signum(yFirst)==h ? yFirst : ySecond; + + default: + return Double.POSITIVE_INFINITY; + } } - private int getHemisphere(Ray ray, Surface surface){ - int h =0; - if(surface.cylinder!=null) { - List trajs = new ArrayList<>(); - Line3D line = ray.toLine(); - if(surface.cylinder.intersection(line, trajs)>= 1) { - h= (int) Math.signum(trajs.get(0).y()); - } + private double handlePlaneIntersection(Line3D line, Plane3D plane) { + Point3D intersection = new Point3D(); + if (plane.intersection(line, intersection) != 0) { + return intersection.y(); } + return Double.POSITIVE_INFINITY; + } + + private Point3D getClosestIntersectionToClusterZ(Cluster cluster, List intersections) { + Point3D first = intersections.get(0); + Point3D second = intersections.get(1); + + double distanceToFirst = Math.abs(cluster.center().z() - first.z()); + double distanceToSecond = Math.abs(cluster.center().z() - second.z()); + + return distanceToFirst < distanceToSecond ? first : second; + } + + private int getHemisphere(Cluster cluster, Ray ray, Surface surface){ + double Y = this.getY(cluster, ray, surface); + if(Y == Double.POSITIVE_INFINITY) return 0; + surface.rayYinterc=Y; + return (int) Math.signum(Y); + + } + + private int getHemisphere(int h, Ray ray, Surface surface){ + double Y = this.getY(h, ray, surface); + if(Y == Double.POSITIVE_INFINITY) return 0; + surface.rayYinterc=Y; + return (int) Math.signum(Y); - if(surface.plane!=null) { - h=(int) surface.hemisphere; + } + + private boolean isCrossed(Ray ray, Surface surface) { + if (surface.cylinder == null) { + return false; + } + + List intersections = new ArrayList<>(); + Line3D line = ray.toLine(); + int intersectionCount = surface.cylinder.intersection(line, intersections); + + if (intersectionCount == 0) { + return false; + } + + // Get the hemisphere sign for comparison + int hemisphereSign = (int) Math.signum(surface.hemisphere); + + switch (intersectionCount) { + case 1: + // Check if the single intersection matches the hemisphere + return hemisphereSign == (int) Math.signum(intersections.get(0).y()); + case 2: + // Check if either of the two intersections match the hemisphere + return hemisphereSign == (int) Math.signum(intersections.get(0).y()) || + hemisphereSign == (int) Math.signum(intersections.get(1).y()); + default: + return false; // If more than 2 intersections, return false (shouldn't happen) } - return h; } + private int getHemisphere(Cluster cluster, Seed seed, Surface surface) { if(surface.cylinder==null) return 0; From bb510112bbe111b4fc9eef4b4cb046b72ee0c2ca Mon Sep 17 00:00:00 2001 From: ziegler Date: Sat, 4 Jan 2025 10:04:10 -0500 Subject: [PATCH 10/12] Cosmic debug statement --- .../java/org/jlab/rec/cvt/services/CVTEngine.java | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/CVTEngine.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/CVTEngine.java index e278b0c9c7..038b736911 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/CVTEngine.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/CVTEngine.java @@ -186,6 +186,18 @@ public int getRun(DataEvent event) { } return run; } + + public void getEventNum(DataEvent event) { + + if (event.hasBank("RUN::config") == false) { + System.err.println("RUN CONDITIONS NOT READ!"); + return ; + } + + DataBank bank = event.getBank("RUN::config"); + System.out.println("EVENT NUMBER "+bank.getInt("event", 0)); + + } public int getPid() { return pid; @@ -314,7 +326,8 @@ public boolean processDataEvent(DataEvent event) { banks.add(RecoBankWriter.fillStraightTracksBank(event, tracks, "CVTRec::Cosmics")); banks.add(RecoBankWriter.fillStraightTracksTrajectoryBank(event, tracks, "CVTRec::Trajectory")); banks.add(RecoBankWriter.fillStraightTrackKFTrajectoryBank(event, tracks, "CVTRec::KFTrajectory")); - } + //if(tracks.isEmpty() && seeds!=null && !seeds.isEmpty()) this.getEventNum(event); + } } else { double[] xyBeam = CVTReconstruction.getBeamSpot(event, beamPos); From 508ee560f791ee5b21caf54898bc09ed6acc99a7 Mon Sep 17 00:00:00 2001 From: ziegler Date: Thu, 9 Jan 2025 19:48:46 -0500 Subject: [PATCH 11/12] Change intersection of ray with surface from y-coordinate to Point3D. Change isCrossed method to use the value of the hemisphere which is calculated as the sign of the y value of the intersection of the ray with the surface if it crosses it or 0 if it does not. --- .../rec/cvt/measurement/Measurements.java | 109 +++++++++--------- 1 file changed, 53 insertions(+), 56 deletions(-) diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java index 461152bcc9..c7f44ff661 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/measurement/Measurements.java @@ -183,7 +183,7 @@ private Surface getCosmicPlane() { cosmic.addMaterial(Geometry.VACUUM); cosmic.setError(1); cosmic.hemisphere = 1; - cosmic.rayYinterc=0; + int inters=plane.intersection(new Line3D(ep1, ep2), cosmic.rayInterc); cosmic.passive = true; return cosmic; } @@ -240,8 +240,8 @@ public List getMeasurements(StraightTrack cosmic) { if(surf.passive && surf.getIndex()!=0) { surf.hemisphere = this.getHemisphere((int)surf.hemisphere, cosmic.getRay(), surf); - if(!this.isCrossed(cosmic.getRay(), surf)) { - if(debug) System.out.println("Removing surface "+surf.toString()+" " + surf.passive + " " + this.isCrossed(cosmic.getRay(), surf)); + if(!this.isCrossed(surf)) { + if(debug) System.out.println("Removing surface "+surf.toString()+" " + surf.passive ); continue; } } @@ -359,7 +359,7 @@ private void addMissing(StraightTrack ray) { surface.hemisphere = this.getHemisphere(hemisphere, ray.getRay(), surface); //resets hemisphere for shallow angle tracks surface.passive=true; if(debug) System.out.println("Generating surface for missing index " + i +" id "+id+ " detector " + type.getName() + " layer " + layer + " sector " + surface.getSector()+" hemisphere "+ - surface.hemisphere+" y "+surface.rayYinterc); + surface.hemisphere+" y "+surface.rayInterc); this.add(i, surface); } } @@ -445,9 +445,9 @@ else if(type==DetectorType.BMT) { - private double getY(Cluster cluster, Ray ray, Surface surface) { + private Point3D getRayInters(Cluster cluster, Ray ray, Surface surface) { if (ray == null || surface == null) { - return Double.POSITIVE_INFINITY; + return new Point3D(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY); } Line3D line = ray.toLine(); // Shared line for both cylinder and plane @@ -459,12 +459,12 @@ private double getY(Cluster cluster, Ray ray, Surface surface) { return handlePlaneIntersection(line, surface.plane); } - return Double.POSITIVE_INFINITY; // Default if no surface type matched + return new Point3D(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY); // Default if no surface type matched } - private double handleCylinderIntersection(Cluster cluster, Cylindrical3D cylinder, Line3D line) { + private Point3D handleCylinderIntersection(Cluster cluster, Cylindrical3D cylinder, Line3D line) { if (cluster.getType() != BMTType.C) { - return cluster.center().y(); // For non-C type clusters, return the center y + return cluster.center(); // For non-C type clusters, return the center } List intersections = new ArrayList<>(); @@ -472,21 +472,21 @@ private double handleCylinderIntersection(Cluster cluster, Cylindrical3D cylinde switch (intersectionCount) { case 0: - return Double.POSITIVE_INFINITY; + return new Point3D(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY); case 1: - return intersections.get(0).y(); + return intersections.get(0); case 2: // Choose the intersection closest to the cluster center on the z-axis Point3D closest = getClosestIntersectionToClusterZ(cluster, intersections); - return closest.y(); + return closest; default: - return Double.POSITIVE_INFINITY; + return new Point3D(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY); } } - private double getY(int h, Ray ray, Surface surface) { + private Point3D getRayInters(int h, Ray ray, Surface surface) { if (ray == null || surface == null) { - return Double.POSITIVE_INFINITY; + return new Point3D(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY); } Line3D line = ray.toLine(); // Shared line for both cylinder and plane @@ -498,37 +498,41 @@ private double getY(int h, Ray ray, Surface surface) { return handlePlaneIntersection(line, surface.plane); } - return Double.POSITIVE_INFINITY; // Default if no surface type matched + return new Point3D(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY); // Default if no surface type matched } - private double handleCylinderIntersection(int h, Cylindrical3D cylinder, Line3D line) { + private Point3D handleCylinderIntersection(int h, Cylindrical3D cylinder, Line3D line) { List intersections = new ArrayList<>(); int intersectionCount = cylinder.intersection(line, intersections); switch (intersectionCount) { case 0: - return Double.POSITIVE_INFINITY; + return new Point3D(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY); case 1: - return intersections.get(0).y(); + return intersections.get(0); case 2: // Choose - double yFirst = intersections.get(0).y(); - double ySecond =intersections.get(1).y(); + Point3D yFirst = intersections.get(0); + Point3D ySecond =intersections.get(1); - return (int) Math.signum(yFirst)==h ? yFirst : ySecond; + if(Math.signum(yFirst.y())==h) { + return yFirst ; + } else { + return ySecond; + } default: - return Double.POSITIVE_INFINITY; + return new Point3D(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY); } } - private double handlePlaneIntersection(Line3D line, Plane3D plane) { + private Point3D handlePlaneIntersection(Line3D line, Plane3D plane) { Point3D intersection = new Point3D(); if (plane.intersection(line, intersection) != 0) { - return intersection.y(); + return intersection; } - return Double.POSITIVE_INFINITY; + return new Point3D(Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY,Double.POSITIVE_INFINITY); } private Point3D getClosestIntersectionToClusterZ(Cluster cluster, List intersections) { @@ -542,50 +546,43 @@ private Point3D getClosestIntersectionToClusterZ(Cluster cluster, List } private int getHemisphere(Cluster cluster, Ray ray, Surface surface){ - double Y = this.getY(cluster, ray, surface); - if(Y == Double.POSITIVE_INFINITY) return 0; - surface.rayYinterc=Y; - return (int) Math.signum(Y); + Point3D Y = this.getRayInters(cluster, ray, surface); + if(Y.y() == Double.POSITIVE_INFINITY) return 0; + surface.rayInterc=Y; + return (int) Math.signum(Y.y()); } private int getHemisphere(int h, Ray ray, Surface surface){ - double Y = this.getY(h, ray, surface); - if(Y == Double.POSITIVE_INFINITY) return 0; - surface.rayYinterc=Y; - return (int) Math.signum(Y); + Point3D Y = this.getRayInters(h, ray, surface); + if(Y.y() == Double.POSITIVE_INFINITY) return 0; + surface.rayInterc=Y; + return (int) Math.signum(Y.y()); } - private boolean isCrossed(Ray ray, Surface surface) { + private boolean isCrossed(Surface surface) { if (surface.cylinder == null) { return false; } - - List intersections = new ArrayList<>(); - Line3D line = ray.toLine(); - int intersectionCount = surface.cylinder.intersection(line, intersections); - - if (intersectionCount == 0) { + //is Crossed is called after getHemisphere which computes the intersection of the ray with the cylinder + // so there is no need to recompute it and it can just use the hemisphere to determine if the intersection was valid + if(surface.hemisphere==0) { return false; + } else { + return true; } - - // Get the hemisphere sign for comparison - int hemisphereSign = (int) Math.signum(surface.hemisphere); - - switch (intersectionCount) { - case 1: - // Check if the single intersection matches the hemisphere - return hemisphereSign == (int) Math.signum(intersections.get(0).y()); - case 2: - // Check if either of the two intersections match the hemisphere - return hemisphereSign == (int) Math.signum(intersections.get(0).y()) || - hemisphereSign == (int) Math.signum(intersections.get(1).y()); - default: - return false; // If more than 2 intersections, return false (shouldn't happen) + + } + //functionality to find out if a surface is crossed for surfaces where the hemisphere has not been set yet (not used at present + private boolean isCrossed(Cluster cluster, Ray ray, Surface surface) { + int h = this.getHemisphere(cluster, ray, surface); + if(h==0) { + return false; + } else { + return true; } } - private int getHemisphere(Cluster cluster, Seed seed, Surface surface) { if(surface.cylinder==null) return 0; From d5699783974fbe62e597f8039daaa50ef5209c95 Mon Sep 17 00:00:00 2001 From: ziegler Date: Thu, 9 Jan 2025 20:35:17 -0500 Subject: [PATCH 12/12] Remove commented out old code. Fix intersection of ray with cylinder for tracks from the target. Intersection initialized as positive infinity. --- .../clas/tracking/kalmanfilter/Surface.java | 6 +++--- .../kalmanfilter/helical/StateVecs.java | 17 +--------------- .../kalmanfilter/straight/KFitter.java | 3 ++- .../kalmanfilter/straight/StateVecs.java | 20 ++++--------------- 4 files changed, 10 insertions(+), 36 deletions(-) diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/Surface.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/Surface.java index 75bc6fbced..d17b281cfa 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/Surface.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/Surface.java @@ -42,7 +42,7 @@ public class Surface implements Comparable { public double swimAccuracy; public boolean passive = false; public double hemisphere = 1; - public double rayYinterc = 9999; + public Point3D rayInterc = new Point3D(Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY, Double.POSITIVE_INFINITY); public double[] unc = new double[2]; public double[] doca = new double[2]; @@ -422,14 +422,14 @@ public void setTransformation(Transformation3D transform) { @Override public int compareTo(Surface o) { - if(this.rayYinterc==9999) { + if(this.rayInterc.y()==Double.POSITIVE_INFINITY) { if (this.index > o.index) { return 1; } else { return -1; } } else { - if (this.rayYinterc < o.rayYinterc) { + if (this.rayInterc.y() < o.rayInterc.y()) { return 1; } else { return -1; diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java index 14eeb15ce2..c9ea7aa7eb 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java @@ -53,23 +53,8 @@ public boolean getStateVecPosAtMeasSite(StateVec sv, AMeasVecs.MeasVec mv, Swim sv.path = inters.distance(st); } else if(mv.surface.cylinder!=null) { -// mv.surface.toLocal().apply(st); -// mv.surface.toLocal().apply(stu); -// double r = mv.surface.cylinder.baseArc().radius(); -// double delta = Math.sqrt((st.x()*stu.x()+st.y()*stu.y())*(st.x()*stu.x()+st.y()*stu.y())-(-r*r+st.x()*st.x()+st.y()*st.y())*(stu.x()*stu.x()+stu.y()*stu.y())); -// double l = (-(st.x()*stu.x()+st.y()*stu.y())+delta)/(stu.x()*stu.x()+stu.y()*stu.y()); -// if(Math.signum(st.y()+l*stu.y())!=mv.hemisphere) { -// l = (-(st.x()*stu.x()+st.y()*stu.y())-delta)/(stu.x()*stu.x()+stu.y()*stu.y()); -// } -// Point3D inters = new Point3D(st.x()+l*stu.x(),st.y()+l*stu.y(),st.z()+l*stu.z()); -// mv.surface.toGlobal().apply(inters); -// // RDV: should switch to use clas-geometry intersection method, not done now to alwys return a value -// sv.x = inters.x(); -// sv.y = inters.y(); -// sv.z = inters.z(); -// sv.path = inters.distance(st); List inters = new ArrayList(); - int ints = mv.surface.cylinder.intersection(toPln, inters); + int ints = mv.surface.cylinder.intersectionRay(toPln, inters); if(ints<1) return false; sv.x = inters.get(0).x() ; diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/KFitter.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/KFitter.java index 003a5008aa..5a2b76ab18 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/KFitter.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/KFitter.java @@ -62,7 +62,8 @@ else if(newchisq < this.chi2) { finalStateVec = sv.new StateVec(sv.smoothed().get(0)); this.setTrajectory(sv, mv); //System.out.println(finalStateVec.toString()); - //setFitFailed = false; + //setFitFailed = false; //allow for new iteration to be worse but save the track from the best iteration. + // There are instances where the next iteration is worse at the 1 per 10000 level. This would loose the track... } // stop if chi2 got worse else { diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java index fd493879aa..5cf50fae3b 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java +++ b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java @@ -44,18 +44,6 @@ else if(mv.hemisphere!=0) { vec.path = inters.distance(ref); } if(mv.surface.cylinder!=null) { - //mv.surface.toLocal().apply(ref); - //mv.surface.toLocal().apply(u); -// double r = 0.5*(mv.surface.cylinder.baseArc().radius()+mv.surface.cylinder.highArc().radius()); -// double delta = Math.sqrt((ref.x()*u.x()+ref.y()*u.y())*(ref.x()*u.x()+ref.y()*u.y()) -// -(-r*r+ref.x()*ref.x()+ref.y()*ref.y())*(u.x()*u.x()+u.y()*u.y())); -// double l = (-(ref.x()*u.x()+ref.y()*u.y())+delta)/(u.x()*u.x()+u.y()*u.y()); -// if(Math.signum(ref.y()+l*u.y())!=mv.hemisphere) { -// l = (-(ref.x()*u.x()+ref.y()*u.y())-delta)/(u.x()*u.x()+u.y()*u.y()); -// } -// -// Point3D cylInt = new Point3D(ref.x()+l*u.x(),ref.y()+l*u.y(),ref.z()+l*u.z()); -// mv.surface.toGlobal().apply(cylInt); List inters = new ArrayList(); int ints = mv.surface.cylinder.intersection(toPln, inters); if(ints<1) { @@ -67,11 +55,11 @@ else if(mv.hemisphere!=0) { } } else { int index =0; - if(ints>1) { - double dh0 = Math.abs(inters.get(0).y()-mv.surface.rayYinterc); - double dh1 = Math.abs(inters.get(1).y()-mv.surface.rayYinterc); + if(ints>1) {//pick the intersection that is closest to the global fit ray intersection + double dh0 = Math.abs(inters.get(0).y()-mv.surface.rayInterc.y()); + double dh1 = Math.abs(inters.get(1).y()-mv.surface.rayInterc.y()); if(dh1