Skip to content

Commit 28065a9

Browse files
tongtongcaotongtongcao
committed
fix some bugs
1 parent 7174be6 commit 28065a9

4 files changed

Lines changed: 144 additions & 89 deletions

File tree

‎reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/RecUtilities.java‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -495,7 +495,7 @@ public boolean reFitCircle(Seed seed, int iter, double xb, double yb) {
495495
if(fitStatus) {
496496
CircleFitPars pars = circlefit.getFit();
497497
seed.getHelix().setCurvature(pars.rho());
498-
seed.getHelix().setDCA(-pars.doca());
498+
seed.getHelix().setDCA(pars.doca());
499499
seed.getHelix().setPhiAtDCA(pars.phi());
500500
seed.update_Crosses(xb,yb);
501501
}

‎reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/TracksFromTargetRec.java‎

Lines changed: 6 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -132,7 +132,7 @@ public List<Seed> getSeeds(List<ArrayList<Cluster>> clusters, List<ArrayList<Cro
132132
if(Constants.getInstance().seedingDebugMode) {
133133
System.out.println("Before overlap remover");
134134
System.out.println(s.toString());
135-
}
135+
}
136136
s.setKey(s.new Key(s));
137137
}
138138
if(Constants.getInstance().removeOverlappingSeeds)
@@ -171,6 +171,11 @@ public List<Seed> getSeeds(List<ArrayList<Cluster>> clusters, List<ArrayList<Cro
171171
}
172172
}
173173

174+
// To match with d0 defition in tracking
175+
for(Seed s : seeds) {
176+
s.getHelix().setDCA(-s.getHelix().getDCA());
177+
}
178+
174179
this.CVTseeds = seeds;
175180
// Got seeds;
176181
return seeds;

‎reconstruction/cvt/src/main/java/org/jlab/rec/cvt/track/TrackSeederXY.java‎

Lines changed: 131 additions & 58 deletions
Original file line numberDiff line numberDiff line change
@@ -33,7 +33,7 @@ public TrackSeederXY(double xb, double yb) {
3333
sortedCrosses = new ArrayList<>();
3434
for(int i =0; i<NBINS; i++) {
3535
sortedCrosses.add(i, new ArrayList<>() );
36-
for(int l =0; l<3; l++) {
36+
for(int l =0; l<6; l++) {
3737
sortedCrosses.get(i).add(l,new ArrayList<>() );
3838
}
3939
}
@@ -100,70 +100,144 @@ public void findSeedCrossList(List<Cross> crosses) {
100100

101101

102102
/*
103-
Scans overphase space to find groups of BMT crosses
103+
Scans overphase space to find groups of SVT and BMT-Z crosses
104104
*/
105105
private void findSeedCrossesFixedBin(List<Cross> crosses, double phiShift) {
106-
for(int b =0; b<NBINS; b++) {
107-
for(int l =0; l<3; l++) {
106+
107+
// Clear previous contents
108+
for (int b = 0; b < NBINS; b++) {
109+
for (int l = 0; l < 6; l++) {
108110
sortedCrosses.get(b).get(l).clear();
109111
}
110112
}
111-
int[][] LPhi = new int[NBINS][3];
112-
for (int i = 0; i < crosses.size(); i++) {
113-
double phi = Math.toDegrees(crosses.get(i).getPoint().toVector3D().phi());
114113

115-
phi += phiShift;
116-
if (phi < 0) {
117-
phi += 360;
114+
// Number of crosses in each bin/layer
115+
int[][] LPhi = new int[NBINS][6];
116+
117+
// ------------------------------------------------------------
118+
// 1. Sort crosses into phi bins and detector layers
119+
// ------------------------------------------------------------
120+
for (Cross cross : crosses) {
121+
122+
double phi = Math.toDegrees(
123+
cross.getPoint().toVector3D().phi()
124+
);
125+
126+
// Shift phi and wrap into [0, 360)
127+
phi = (phi + phiShift + 360.0) % 360.0;
128+
129+
int binIdx = (int) (phi / (360.0 / NBINS));
130+
131+
// Protect against numerical rounding
132+
if (binIdx >= NBINS) {
133+
binIdx = NBINS - 1;
118134
}
119135

120-
int binIdx = (int) (phi / (360./NBINS) );
121-
122-
if(binIdx>NBINS-1)
123-
binIdx = NBINS-1;
124-
sortedCrosses.get(binIdx).get(crosses.get(i).getRegion() - 1).add(crosses.get(i));
125-
LPhi[binIdx][crosses.get(i).getRegion() - 1]++;
136+
int layer = -1;
137+
138+
if (cross.getDetector() == DetectorType.BST) {
139+
140+
// BST region 1,2,3 -> layer 0,1,2
141+
layer = cross.getRegion() - 1;
142+
143+
} else if (cross.getDetector() == DetectorType.BMT) {
144+
145+
// BMT region 1,2,3 -> layer 3,4,5
146+
layer = cross.getRegion() + 2;
147+
}
148+
149+
// Ignore unexpected detector/region
150+
if (layer < 0 || layer >= 6) {
151+
continue;
152+
}
153+
154+
sortedCrosses.get(binIdx).get(layer).add(cross);
155+
LPhi[binIdx][layer]++;
126156
}
127-
157+
158+
// ------------------------------------------------------------
159+
// 2. For each phi bin, generate ALL combinations
160+
// ------------------------------------------------------------
128161
for (int b = 0; b < NBINS; b++) {
129-
int max_layers =0;
130-
for (int la = 0; la < 3; la++) {
131-
if(LPhi[b][la]>0)
132-
max_layers++;
162+
163+
// Find occupied layers
164+
List<Integer> occupiedLayers = new ArrayList<>();
165+
166+
for (int la = 0; la < 6; la++) {
167+
if (LPhi[b][la] > 0) {
168+
occupiedLayers.add(la);
169+
}
170+
}
171+
172+
// Need at least two different layers
173+
if (occupiedLayers.size() < 2) {
174+
continue;
133175
}
134-
if (sortedCrosses.get(b) != null && max_layers >= 2) {
135-
double SumLyr=0;
136-
while(LPhi[b][0]+LPhi[b][1]+ LPhi[b][2]>=max_layers) {
137-
if(SumLyr!=LPhi[b][0]+LPhi[b][1]+ LPhi[b][2]) {
138-
SumLyr = LPhi[b][0]+LPhi[b][1]+ LPhi[b][2];
139-
}
140-
ArrayList<Cross> hits = new ArrayList<>();
141-
for (int la = 0; la < 3; la++) {
142-
if (sortedCrosses.get(b).get(la) != null && LPhi[b][la]>0) {
143-
if (sortedCrosses.get(b).get(la).get(LPhi[b][la]-1) != null
144-
&& sortedCrosses.get(b).get(la).size()>0) {
145-
hits.add(sortedCrosses.get(b).get(la).get(LPhi[b][la]-1));
146-
147-
if(LPhi[b][la]>1)
148-
LPhi[b][la]--;
149-
if(SumLyr==max_layers)
150-
LPhi[b][la]=0;
151-
}
152-
}
153-
}
154-
155-
if (hits.size() >= 2) {
156-
double seedIdx=0;
157-
int s = hits.size();
158-
int index = (int) Math.pow(2,s);
159-
for(Cross c : hits) {
160-
seedIdx +=c.getId()*Math.pow(10, index);
161-
index-=4;
162-
}
163-
seedMap.put(seedIdx, hits);
164-
}
176+
177+
// Cross combos for any 2 more layers
178+
for (int nLayers = 2; nLayers <= occupiedLayers.size(); nLayers++) {
179+
generateLayerCombinations(b, occupiedLayers, 0, nLayers, new ArrayList<>());
180+
}
181+
}
182+
}
183+
184+
/*
185+
* Generate combinations of occupied layers.
186+
*
187+
*/
188+
private void generateLayerCombinations(int binIdx, List<Integer> occupiedLayers, int start, int nLayers, ArrayList<Integer> selectedLayers) {
189+
190+
// ------------------------------------------------------------
191+
// Select the required number of layers
192+
// ------------------------------------------------------------
193+
if (selectedLayers.size() == nLayers) {
194+
ArrayList<Cross> currentCombination = new ArrayList<>();
195+
generateCrossCombinations(binIdx, selectedLayers, 0, currentCombination);
196+
return;
197+
}
198+
199+
// ------------------------------------------------------------
200+
// Select next layer
201+
// ------------------------------------------------------------
202+
for (int i = start; i < occupiedLayers.size(); i++) {
203+
selectedLayers.add(occupiedLayers.get(i));
204+
generateLayerCombinations(binIdx, occupiedLayers, i + 1, nLayers, selectedLayers);
205+
selectedLayers.remove(selectedLayers.size() - 1);
206+
}
207+
}
208+
209+
/*
210+
* Generate Cartesian products of Crosses from the selected layers.
211+
*
212+
*/
213+
private void generateCrossCombinations(int binIdx, List<Integer> selectedLayers, int layerIndex, ArrayList<Cross> currentCombination) {
214+
// ------------------------------------------------------------
215+
// All selected layers have been processed
216+
// ------------------------------------------------------------
217+
if (layerIndex == selectedLayers.size()) {
218+
if (currentCombination.size() >= 2) {
219+
ArrayList<Cross> crosses = new ArrayList<>(currentCombination);
220+
long seedIdx = 0;
221+
for (Cross c : crosses) {
222+
seedIdx = seedIdx * 10000L + c.getId();
165223
}
224+
225+
seedMap.put((double) seedIdx, crosses);
166226
}
227+
return;
228+
}
229+
230+
// Current detector layer
231+
int layer = selectedLayers.get(layerIndex);
232+
List<Cross> layerCrosses = sortedCrosses.get(binIdx).get(layer);
233+
234+
// ------------------------------------------------------------
235+
// Try EVERY Cross in this layer
236+
// ------------------------------------------------------------
237+
for (Cross cross : layerCrosses) {
238+
currentCombination.add(cross);
239+
generateCrossCombinations(binIdx, selectedLayers, layerIndex + 1,currentCombination);
240+
currentCombination.remove(currentCombination.size() - 1);
167241
}
168242
}
169243

@@ -237,15 +311,14 @@ private int countType(List<Cross> cand, DetectorType dt) {
237311
}
238312
private double calcResi(double rho, double d0, double phi0, double xc, double yc) {
239313
double r = Math.sqrt(xc*xc+yc*yc);
240-
double par = 1. - ((r * r - d0 * d0) * rho * rho) / (2. * (1. + d0 * Math.abs(rho)));
314+
double par = (2. * (1. + d0 * rho) - (r * r - d0 * d0) * rho * rho) / (2. * Math.abs(1. + d0 * rho));
315+
241316
double newPathLength = Math.abs(Math.acos(par) / rho);
242317
int charge = (int) Math.signum(rho);
243318
double alpha = -newPathLength * rho;
244-
245-
double x = d0 * charge * Math.sin(phi0) + (charge / Math.abs(rho))
246-
* (Math.sin(phi0) - Math.cos(alpha) * Math.sin(phi0) - Math.sin(alpha) * Math.cos(phi0));
247-
double y = -d0 * charge * Math.cos(phi0) - (charge / Math.abs(rho))
248-
* (Math.cos(phi0) + Math.sin(alpha) * Math.sin(phi0) - Math.cos(alpha) * Math.cos(phi0));
319+
320+
double x = d0 * Math.sin(phi0) + 1./rho * (Math.sin(phi0) - Math.sin(phi0 + alpha));
321+
double y = -d0 * Math.cos(phi0) - 1./rho * (Math.cos(phi0) - Math.cos(phi0 + alpha));
249322

250323
double res2 = ((x-xc)*(x-xc)+(y-yc)*(y-yc));
251324
return Math.sqrt(res2);

‎reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/Helix.java‎

Lines changed: 6 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -178,15 +178,6 @@ public double radius() {
178178
return 1. / C;
179179
}
180180

181-
// (x,y) coordinates of the circle center
182-
public double xCen() {
183-
return (radius() - this.getDCA()) * Math.sin(this.getPhiAtDCA());
184-
}
185-
186-
public double yCen() {
187-
return (-radius() + this.getDCA()) * Math.cos(this.getPhiAtDCA());
188-
}
189-
190181
// (x,y) coordinates of the dca
191182
public double xDCA() {
192183
return -this.getDCA() * Math.sin(this.getPhiAtDCA());
@@ -198,7 +189,7 @@ public double yDCA() {
198189

199190
public Point3D getVertex() {
200191
return new Point3D(this.xDCA()+this.getXb(),this.yDCA()+this.getYb(),this.getZ0());
201-
}
192+
}
202193

203194
public double getPt(double solenoidMag) {
204195
double pt = Constants.LIGHTVEL * this.radius() * solenoidMag;
@@ -214,18 +205,7 @@ public Vector3D getPXYZ(double solenoidMag) {
214205
double py = pt*Math.sin(this.getPhiAtDCA());
215206

216207
return new Vector3D(px,py,pz);
217-
}
218-
219-
public double getArcLengthDCA(Point3D refpoint) {
220-
//insure that the refpoint is on the helix
221-
if (refpoint == null) {
222-
return 0;
223-
}
224-
double refX = radius() * Math.cos(refpoint.toVector3D().phi());
225-
double refY = radius() * Math.sin(refpoint.toVector3D().phi());
226-
double arclen = arcLength(xCen(), yCen(), radius(), xCen(), yCen(), refX, refY);
227-
return arclen;
228-
}
208+
}
229209

230210
// this method finds the arclength between 2 points in a circle
231211
// 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) {
272252
return new Point3D(x, y, z);
273253
}
274254

275-
276-
double par = 1. - ((r * r - d0 * d0) * omega * omega) / (2. * (1. + d0 * Math.abs(omega)));
255+
double par = (2. * (1. + d0 * omega) - (r * r - d0 * d0) * omega * omega) / (2. * Math.abs(1. + d0 * omega));
277256
double newPathLength = Math.abs(Math.acos(par) / omega);
278257

279258
double alpha = -newPathLength * omega;
280259

281-
double x = d0 * charge * Math.sin(phi0) + (charge / Math.abs(omega))
282-
* (Math.sin(phi0) - Math.cos(alpha) * Math.sin(phi0) - Math.sin(alpha) * Math.cos(phi0))+xb;
283-
double y = -d0 * charge * Math.cos(phi0) - (charge / Math.abs(omega))
284-
* (Math.cos(phi0) + Math.sin(alpha) * Math.sin(phi0) - Math.cos(alpha) * Math.cos(phi0))+yb;
260+
double x = d0 * Math.sin(phi0) + 1./omega * (Math.sin(phi0) - Math.sin(phi0 + alpha));
261+
double y = -d0 * Math.cos(phi0) - 1./omega * (Math.cos(phi0) - Math.cos(phi0 + alpha));
285262
double z = z0 + newPathLength * tandip;
286263

287264
return new Point3D(x, y, z);
@@ -301,7 +278,7 @@ public Vector3D getTrackDirectionAtRadius(double r) {
301278
return new Vector3D(ux, uy, uz);
302279
}
303280

304-
double par = 1. - ((r * r - d0 * d0) * omega * omega) / (2. * (1. + d0 * Math.abs(omega)));
281+
double par = (2. * (1. + d0 * omega) - (r * r - d0 * d0) * omega * omega) / (2. * Math.abs(1. + d0 * omega));
305282
double newPathLength = Math.abs(Math.acos(par) / omega);
306283

307284
double alpha = newPathLength * omega;

0 commit comments

Comments
 (0)