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
wants to merge 12 commits into
base: development
Choose a base branch
from
+255 −91
Open

Iss388 cvt cosmics #421

Show file tree
Hide file tree
Changes from 10 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 number Diff line number Diff 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 number Diff line number Diff 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 double rayYinterc = 9999;
Copy link
Collaborator

@raffaelladevita raffaelladevita Jan 9, 2025

Choose a reason for hiding this comment

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

let's make it inf or nan

Copy link
Collaborator

@raffaelladevita raffaelladevita Jan 9, 2025

Choose a reason for hiding this comment

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

Also, to be more general, it could make sense to add a Point3D to store the full information of the trajectory point instead of only y


public double[] unc = new double[2];
public double[] doca = new double[2];
Expand Down Expand 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.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;
}
}
}

}
47 changes: 28 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 number Diff line number Diff 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 Down Expand 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,31 @@ 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);
Copy link
Collaborator

@raffaelladevita raffaelladevita Jan 9, 2025

Choose a reason for hiding this comment

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

let's remove the commented-out lines...they will still be in the history

// 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.intersection(toPln, inters);
if(ints<1) return false;

sv.x = inters.get(0).x() ;
Copy link
Collaborator

@raffaelladevita raffaelladevita Jan 9, 2025

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
12 changes: 7 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 number Diff line number Diff line change
Expand Up @@ -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;
Copy link
Collaborator

@raffaelladevita raffaelladevita Jan 9, 2025

Choose a reason for hiding this comment

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

I assume setFitFailed is not really used elsewhere but it may be better to still set it to false
(and remove the commented text)

}
// 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
55 changes: 37 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 number Diff line number Diff 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,41 @@ 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());
//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<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) {
double dh0 = Math.abs(inters.get(0).y()-mv.surface.rayYinterc);
Copy link
Collaborator

@raffaelladevita raffaelladevita Jan 9, 2025

Choose a reason for hiding this comment

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

Shouldn't the right intersection be the one with the largest y among the ones with y<rayInterc?

double dh1 = Math.abs(inters.get(1).y()-mv.surface.rayYinterc);
if(dh1<dh0)
index=1;
}
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 Down Expand Up @@ -172,7 +192,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
Loading