Skip to content

Iss388 cvt cosmics#421

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
zieglerv wants to merge 12 commits into
base:development
Choose a base branch
from
Open

Iss388 cvt cosmics #421

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
56cd809
Bug fix in setting the measurements for cosmics. Added missing hemis…
Dec 12, 2024
f1f1f02
fix in propagation to surface for shifted cylinders.
Dec 13, 2024
6bf24ee
Fix measurements sorting by Y intersection of the cosmic ray with the…
Dec 13, 2024
a1d3d37
removed unused import
Jan 4, 2025
6aaf945
Initialized chi^2
Jan 4, 2025
861caf3
Added y intersect for rays. This is used to sort the surfaces from t…
Jan 4, 2025
4fb1e22
For cosmic fitting, save the output for the best chi^2; do not kill t…
Jan 4, 2025
4528930
Use global fit y-intersection with the surface to assign fitted track…
Jan 4, 2025
92952b4
Assign correct hemisphere for cosmics.
Jan 4, 2025
bb51011
Cosmic debug statement
Jan 4, 2025
508ee56
Change intersection of ray with surface from y-coordinate to Point3D.
Jan 10, 2025
d569978
Remove commented out old code.
Jan 10, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/AKFitter.java
View file
Open in desktop
Original file line numberDiff line numberDiff line change
Expand Up@@ -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;
Expand Down
18 changes: 13 additions & 5 deletions common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/Surface.java
View file
Open in desktop
Original file line numberDiff line numberDiff line change
Expand Up@@ -42,6 +42,7 @@ public class Surface implements Comparable<Surface> {
public double swimAccuracy;
public boolean passive = false;
public double hemisphere = 1;
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];
Expand DownExpand Up@@ -421,11 +422,18 @@ public void setTransformation(Transformation3D transform) {

@Override
public int compareTo(Surface o) {
if (this.index > o.index) {
return 1;
if(this.rayInterc.y()==Double.POSITIVE_INFINITY) {
if (this.index > o.index) {
return 1;
} else {
return -1;
}
} else {
return -1;
}
if (this.rayInterc.y() < o.rayInterc.y()) {
return 1;
} else {
return -1;
}
}
}

}
32 changes: 13 additions & 19 deletions ...ls/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java
View file
Open in desktop
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,14 @@
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;
import org.jlab.clas.tracking.kalmanfilter.AStateVecs;
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;
Expand DownExpand Up@@ -42,9 +43,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();
Expand All@@ -53,22 +53,16 @@ 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<Point3D> inters = new ArrayList();
int ints = mv.surface.cylinder.intersectionRay(toPln, inters);
if(ints<1) return false;

sv.x = inters.get(0).x() ;

@raffaelladevitaraffaelladevitaJan 9, 2025

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

is the first intersection always the right one?

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);
Expand Down
13 changes: 8 additions & 5 deletions ...ols/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/KFitter.java
View file
Open in desktop
Original file line numberDiff line numberDiff line change
Expand Up@@ -54,22 +54,25 @@ 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; //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 {
break;
}
}
if(!this.setFitFailed) {
finalStateVec = sv.new StateVec(sv.smoothed().get(0));
}
//if(!this.setFitFailed) {
// finalStateVec = sv.new StateVec(sv.smoothed().get(0));
//}
}

@Override
Expand Down
43 changes: 25 additions & 18 deletions ...s/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java
View file
Open in desktop
Original file line numberDiff line numberDiff line change
@@ -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;
Expand All@@ -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;
Expand All@@ -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() ;
Expand All@@ -43,22 +44,29 @@ 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());
List<Point3D> inters = new ArrayList();
int ints = mv.surface.cylinder.intersection(toPln, inters);
if(ints<1) {
if(!mv.surface.passive) {
return false;
} else {
mv.skip=true;
return true;
}

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);
} else {
int index =0;
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<dh0)
index=1; //if the second intersection is closest use its index
}
vec.x = inters.get(index).x() ;
vec.y = inters.get(index).y() ;
vec.z = inters.get(index).z() ;
vec.path = inters.get(index).distance(ref);
}

}
return true;
}
Expand DownExpand Up@@ -172,7 +180,6 @@ public void init(double x0, double z0, double tx, double tz, Units units, double
this.trackTrajB.clear();
this.trackTrajS.clear();
this.trackTrajT.put(0, new StateVec(initSV));

}

@Override
Expand Down
Loading