Skip to content

Commit fd4f9d9

Browse files
zyzhangjlabZeyu Zhangftouchte
authored
Fix Nhits issue & fill track bank (#1137)
Co-authored-by: Zeyu Zhang <[email protected]> Co-authored-by: Felix Touchte Codjo <[email protected]>
1 parent bd8fe01 commit fd4f9d9

4 files changed

Lines changed: 127 additions & 14 deletions

File tree

‎reconstruction/alert/src/main/java/org/jlab/rec/ahdc/HelixFit/HelixFitJava.java‎

Lines changed: 81 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -997,7 +997,81 @@ public HelixFitObject HelixFit(int PointNum, double szPos[][], int fit_track_to_
997997
return h;
998998
}
999999

1000-
// Public weighted wrapper (similar to the original HelixFit)
1000+
// Compute the path length along the helix described by h,
1001+
// for the portion where the radial distance to the beamline
1002+
// (sqrt(x^2 + y^2)) is between rMin and rMax [mm].
1003+
private double computePathFromHelix(HelixFitObject h, double rMin, double rMax) {
1004+
1005+
// Helix parameters from the fit
1006+
double A = h.get_A(); // circle center X in XY-plane
1007+
double B = h.get_B(); // circle center Y in XY-plane
1008+
double X0 = h.get_X0(); // point of closest approach in XY
1009+
double Y0 = h.get_Y0();
1010+
double thDeg = h.get_Theta(); // already stored in degrees in the wrapper
1011+
double theta = Math.toRadians(thDeg); // convert to radians for sin()
1012+
1013+
// Circle radius in XY-plane
1014+
double R = Math.sqrt((X0 - A)*(X0 - A) + (Y0 - B)*(Y0 - B));
1015+
if (R < 1e-6) {
1016+
// Degenerate circle
1017+
return 0.0;
1018+
}
1019+
1020+
// Distance from origin to circle center
1021+
double C = Math.sqrt(A*A + B*B);
1022+
1023+
// Minimal / maximal distance from the circle to the origin
1024+
double Dmin = Math.abs(C - R);
1025+
double Dmax = C + R;
1026+
1027+
// If there is no overlap between [rMin, rMax] and [Dmin, Dmax], path = 0
1028+
if (rMax <= Dmin || rMin >= Dmax) {
1029+
return 0.0;
1030+
}
1031+
1032+
// Clip the requested radial window to what the circle actually covers
1033+
double rIn = Math.max(rMin, Dmin);
1034+
double rOut = Math.min(rMax, Dmax);
1035+
if (rOut <= rIn) {
1036+
return 0.0;
1037+
}
1038+
1039+
// Helper constants for r^2(psi) = C0 + 2 R C cos(psi)
1040+
double C0 = A*A + B*B + R*R;
1041+
1042+
// Helper to compute psi(r) in [0,π] from r:
1043+
// cos(psi) = (r^2 - C0)/(2 R C)
1044+
java.util.function.DoubleFunction<Double> psiOfR = (double r) -> {
1045+
double num = r*r - C0;
1046+
double den = 2.0 * R * C;
1047+
if (Math.abs(den) < 1e-12) {
1048+
return 0.0;
1049+
}
1050+
double u = num / den;
1051+
// Clamp to [-1,1] to avoid NaN from acos
1052+
if (u > 1.0) u = 1.0;
1053+
if (u < -1.0) u = -1.0;
1054+
return Math.acos(u); // in [0, π]
1055+
};
1056+
1057+
double psiIn = psiOfR.apply(rIn);
1058+
double psiOut = psiOfR.apply(rOut);
1059+
1060+
double dPsi = Math.abs(psiOut - psiIn);
1061+
1062+
// ds/dpsi = R / sin(theta)
1063+
double sinTheta = Math.sin(theta);
1064+
if (Math.abs(sinTheta) < 1e-6) {
1065+
// Track almost parallel to the beam; this formula is ill-defined.
1066+
// Return 0 for safety or some large sentinel.
1067+
return 0.0;
1068+
}
1069+
1070+
double path = (R / sinTheta) * dPsi;
1071+
1072+
return path;
1073+
}
1074+
// Public weighted wrapper (similar to the original HelixFit)
10011075
public HelixFitObject HelixFitWeighted(int PointNum, double[][] szPos,
10021076
double[] weights, int fit_track_to_beamline) {
10031077

@@ -1188,14 +1262,17 @@ public HelixFitObject helix_fit_with_doca_selection(List<DocaCluster> docaCluste
11881262
}
11891263

11901264
// weights = null => all weights = 1.0 inside helix_fit_weighted
1265+
// Use the DOCa-based bestChi2 as helix chi2
11911266

11921267
double PI=Math.acos(0.0)*2;
11931268
//double Rho=0,Phi=0,Theta=0,X0=0,Y0=0,DCA=0,Chi2=0;
11941269
double Phi_deg;
11951270
double Theta_deg;
11961271

1197-
HelixFitObject h = helix_fit_weighted(nPoints, szPosSel, null, fit_track_to_beamline);;
1198-
1272+
HelixFitObject h = helix_fit_weighted(nPoints, szPosSel, null, fit_track_to_beamline);
1273+
h.set_Chi2(bestChi2);
1274+
double path = computePathFromHelix(h, 30.0, 70.0);
1275+
h.set_path(path);
11991276
Phi_deg=Math.toDegrees(h.get_Phi());
12001277
if(Phi_deg >= 180){
12011278
Phi_deg -= 360;
@@ -1211,6 +1288,4 @@ public HelixFitObject helix_fit_with_doca_selection(List<DocaCluster> docaCluste
12111288

12121289
}
12131290

1214-
}
1215-
1216-
1291+
}

‎reconstruction/alert/src/main/java/org/jlab/rec/ahdc/HelixFit/HelixFitObject.java‎

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,8 @@ public class HelixFitObject {
1616
private double _DCA;
1717
private double _Chi2;
1818
private double _magfield;
19-
19+
private double _path;
20+
2021
public HelixFitObject(){
2122
//default constructor
2223
}
@@ -33,6 +34,7 @@ public HelixFitObject(double Rho, double A, double B, double Phi, double Theta,
3334
_DCA = DCA;
3435
_Chi2 = Chi2;
3536
_magfield = 50;
37+
_path = 0;
3638
}
3739
public double get_Rho(){
3840
return _Rho;
@@ -82,6 +84,9 @@ public double get_DCA(){
8284
public double get_Chi2(){
8385
return _Chi2;
8486
}
87+
public void set_Chi2(double Chi2){
88+
_Chi2 = Chi2;
89+
}
8590
public double get_Mom(){
8691
return 0.3*_magfield*Math.abs(_Rho)/(10*Math.sin(Math.toRadians(_Theta)));
8792
}
@@ -103,5 +108,10 @@ public double get_dEdx(){
103108
public void set_magfield(double magfield){
104109
_magfield = magfield;
105110
}
106-
111+
public double get_path(){
112+
return _path;
113+
}
114+
public void set_path(double path){
115+
_path = path;
116+
}
107117
}

‎reconstruction/alert/src/main/java/org/jlab/rec/ahdc/Track/Track.java‎

Lines changed: 31 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -73,6 +73,8 @@ public void setPositionAndMomentum(HelixFitObject helixFitObject) {
7373
this.px0 = helixFitObject.get_px();
7474
this.py0 = helixFitObject.get_py();
7575
this.pz0 = helixFitObject.get_pz();
76+
this.chi2 = helixFitObject.get_Chi2();
77+
this.path = helixFitObject.get_path();
7678
}
7779

7880
public void setPositionAndMomentumVec(double[] x) {
@@ -176,15 +178,40 @@ public void set_trackId(int _trackId) {
176178
public void set_chi2(double _chi2) { chi2 = _chi2;}
177179
public void set_sum_residuals(double _sum_residuals) { sum_residuals = _sum_residuals;}
178180
public int get_trackId() {return trackId;}
179-
public int get_n_hits() {return n_hits;}
180-
public int get_sum_adc() {return sum_adc;}
181+
public int get_n_hits() {
182+
if (hits == null) {
183+
return 0;
184+
}
185+
return hits.size();
186+
}
187+
public int get_sum_adc() {
188+
if (hits == null || hits.isEmpty()) {
189+
return 0;
190+
}
191+
int sum = 0;
192+
for (Hit h : hits) {
193+
sum += (int) Math.round(h.getADC());
194+
}
195+
return sum;
196+
}
181197
public double get_chi2() {return chi2;}
182198
public double get_sum_residuals() {return sum_residuals;}
183199
// AHDC::track
184200
public void set_dEdx(double _dEdx) { dEdx = _dEdx;}
185201
public void set_p_drift(double _p_drift) { p_drift = _p_drift;}
186202
public void set_path(double _path) { path = _path;}
187-
public double get_dEdx() {return dEdx;}
203+
public double get_dEdx() {
204+
if (path <= 0) {
205+
dEdx = 0;
206+
}else {
207+
int sum = 0;
208+
for (Hit h : hits) {
209+
sum += (int) Math.round(h.getADC());
210+
}
211+
dEdx = sum/path;
212+
}
213+
return dEdx;
214+
}
188215
public double get_p_drift() {return p_drift;}
189216
public double get_path() {return path;}
190217

@@ -196,4 +223,4 @@ public void set_trackId(int _trackId) {
196223
public int get_predicted_ATOF_layer() {return predicted_ATOF_layer;}
197224
public int get_predicted_ATOF_wedge() {return predicted_ATOF_wedge;}
198225

199-
}
226+
}

‎reconstruction/alert/src/main/java/org/jlab/service/ahdc/AHDCEngine.java‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -218,13 +218,14 @@ else if (modeTrackFinding == ModeTrackFinding.CV_Hough) {
218218
// V) Global fit
219219
int trackid = 0;
220220
ArrayList<DocaCluster> all_docaClusters = new ArrayList<>();
221+
AHDC_Tracks.removeIf(track -> track.get_Clusters().size() < 3);
221222
for (Track track : AHDC_Tracks) {
222223
trackid++;
223224
track.set_trackId(trackid);
224225
List<Cluster> originalClusters = track.get_Clusters();
225226
ArrayList<DocaCluster> docaClusters = DocaClusterRefiner.buildRefinedClusters(originalClusters);
226227
all_docaClusters.addAll(docaClusters);
227-
if (docaClusters == null || docaClusters.size() < 3) {
228+
if (docaClusters == null || docaClusters.size() < 3 || originalClusters == null || originalClusters.size() < 3) {
228229
// not enough points, skip helix fit
229230
continue;
230231
}
@@ -314,4 +315,4 @@ public static void main(String[] args) {
314315

315316
System.out.println("finished " + (System.nanoTime() - starttime) * Math.pow(10, -9));
316317
}
317-
}
318+
}

0 commit comments

Comments
 (0)