Skip to content

fix some bugs in cvt reconstructions #1392

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
tongtongcao wants to merge 1 commit into
base: development
Choose a base branch
from
Open

fix some bugs in cvt reconstructions #1392

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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 reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/RecUtilities.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -495,7 +495,7 @@ public boolean reFitCircle(Seed seed, int iter, double xb, double yb) {
if(fitStatus) {
CircleFitPars pars = circlefit.getFit();
seed.getHelix().setCurvature(pars.rho());
seed.getHelix().setDCA(-pars.doca());
seed.getHelix().setDCA(pars.doca());
seed.getHelix().setPhiAtDCA(pars.phi());
seed.update_Crosses(xb,yb);
}
Expand Down
7 changes: 6 additions & 1 deletion reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/TracksFromTargetRec.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -132,7 +132,7 @@ public List<Seed> getSeeds(List<ArrayList<Cluster>> clusters, List<ArrayList<Cro
if(Constants.getInstance().seedingDebugMode) {
System.out.println("Before overlap remover");
System.out.println(s.toString());
}
}
s.setKey(s.new Key(s));
}
if(Constants.getInstance().removeOverlappingSeeds)
Expand Down Expand Up @@ -171,6 +171,11 @@ public List<Seed> getSeeds(List<ArrayList<Cluster>> clusters, List<ArrayList<Cro
}
}

// To match with d0 defition in tracking
for(Seed s : seeds) {
s.getHelix().setDCA(-s.getHelix().getDCA());
}

this.CVTseeds = seeds;
// Got seeds;
return seeds;
Expand Down
189 changes: 131 additions & 58 deletions reconstruction/cvt/src/main/java/org/jlab/rec/cvt/track/TrackSeederXY.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ public TrackSeederXY(double xb, double yb) {
sortedCrosses = new ArrayList<>();
for(int i =0; i<NBINS; i++) {
sortedCrosses.add(i, new ArrayList<>() );
for(int l =0; l<3; l++) {
for(int l =0; l<6; l++) {
sortedCrosses.get(i).add(l,new ArrayList<>() );
}
}
Expand Down Expand Up @@ -100,70 +100,144 @@ public void findSeedCrossList(List<Cross> crosses) {


/*
Scans overphase space to find groups of BMT crosses
Scans overphase space to find groups of SVT and BMT-Z crosses
*/
private void findSeedCrossesFixedBin(List<Cross> crosses, double phiShift) {
for(int b =0; b<NBINS; b++) {
for(int l =0; l<3; l++) {

// Clear previous contents
for (int b = 0; b < NBINS; b++) {
for (int l = 0; l < 6; l++) {
sortedCrosses.get(b).get(l).clear();
}
}
int[][] LPhi = new int[NBINS][3];
for (int i = 0; i < crosses.size(); i++) {
double phi = Math.toDegrees(crosses.get(i).getPoint().toVector3D().phi());

phi += phiShift;
if (phi < 0) {
phi += 360;
// Number of crosses in each bin/layer
int[][] LPhi = new int[NBINS][6];

// ------------------------------------------------------------
// 1. Sort crosses into phi bins and detector layers
// ------------------------------------------------------------
for (Cross cross : crosses) {

double phi = Math.toDegrees(
cross.getPoint().toVector3D().phi()
);

// Shift phi and wrap into [0, 360)
phi = (phi + phiShift + 360.0) % 360.0;

int binIdx = (int) (phi / (360.0 / NBINS));

// Protect against numerical rounding
if (binIdx >= NBINS) {
binIdx = NBINS - 1;
}

int binIdx = (int) (phi / (360./NBINS) );

if(binIdx>NBINS-1)
binIdx = NBINS-1;
sortedCrosses.get(binIdx).get(crosses.get(i).getRegion() - 1).add(crosses.get(i));
LPhi[binIdx][crosses.get(i).getRegion() - 1]++;
int layer = -1;

if (cross.getDetector() == DetectorType.BST) {

// BST region 1,2,3 -> layer 0,1,2
layer = cross.getRegion() - 1;

} else if (cross.getDetector() == DetectorType.BMT) {

// BMT region 1,2,3 -> layer 3,4,5
layer = cross.getRegion() + 2;
}

// Ignore unexpected detector/region
if (layer < 0 || layer >= 6) {
continue;
}

sortedCrosses.get(binIdx).get(layer).add(cross);
LPhi[binIdx][layer]++;
}


// ------------------------------------------------------------
// 2. For each phi bin, generate ALL combinations
// ------------------------------------------------------------
for (int b = 0; b < NBINS; b++) {
int max_layers =0;
for (int la = 0; la < 3; la++) {
if(LPhi[b][la]>0)
max_layers++;

// Find occupied layers
List<Integer> occupiedLayers = new ArrayList<>();

for (int la = 0; la < 6; la++) {
if (LPhi[b][la] > 0) {
occupiedLayers.add(la);
}
}

// Need at least two different layers
if (occupiedLayers.size() < 2) {
continue;
}
if (sortedCrosses.get(b) != null && max_layers >= 2) {
double SumLyr=0;
while(LPhi[b][0]+LPhi[b][1]+ LPhi[b][2]>=max_layers) {
if(SumLyr!=LPhi[b][0]+LPhi[b][1]+ LPhi[b][2]) {
SumLyr = LPhi[b][0]+LPhi[b][1]+ LPhi[b][2];
}
ArrayList<Cross> hits = new ArrayList<>();
for (int la = 0; la < 3; la++) {
if (sortedCrosses.get(b).get(la) != null && LPhi[b][la]>0) {
if (sortedCrosses.get(b).get(la).get(LPhi[b][la]-1) != null
&& sortedCrosses.get(b).get(la).size()>0) {
hits.add(sortedCrosses.get(b).get(la).get(LPhi[b][la]-1));

if(LPhi[b][la]>1)
LPhi[b][la]--;
if(SumLyr==max_layers)
LPhi[b][la]=0;
}
}
}

if (hits.size() >= 2) {
double seedIdx=0;
int s = hits.size();
int index = (int) Math.pow(2,s);
for(Cross c : hits) {
seedIdx +=c.getId()*Math.pow(10, index);
index-=4;
}
seedMap.put(seedIdx, hits);
}

// Cross combos for any 2 more layers
for (int nLayers = 2; nLayers <= occupiedLayers.size(); nLayers++) {
generateLayerCombinations(b, occupiedLayers, 0, nLayers, new ArrayList<>());
}
}
}

/*
* Generate combinations of occupied layers.
*
*/
private void generateLayerCombinations(int binIdx, List<Integer> occupiedLayers, int start, int nLayers, ArrayList<Integer> selectedLayers) {

// ------------------------------------------------------------
// Select the required number of layers
// ------------------------------------------------------------
if (selectedLayers.size() == nLayers) {
ArrayList<Cross> currentCombination = new ArrayList<>();
generateCrossCombinations(binIdx, selectedLayers, 0, currentCombination);
return;
}

// ------------------------------------------------------------
// Select next layer
// ------------------------------------------------------------
for (int i = start; i < occupiedLayers.size(); i++) {
selectedLayers.add(occupiedLayers.get(i));
generateLayerCombinations(binIdx, occupiedLayers, i + 1, nLayers, selectedLayers);
selectedLayers.remove(selectedLayers.size() - 1);
}
}

/*
* Generate Cartesian products of Crosses from the selected layers.
*
*/
private void generateCrossCombinations(int binIdx, List<Integer> selectedLayers, int layerIndex, ArrayList<Cross> currentCombination) {
// ------------------------------------------------------------
// All selected layers have been processed
// ------------------------------------------------------------
if (layerIndex == selectedLayers.size()) {
if (currentCombination.size() >= 2) {
ArrayList<Cross> crosses = new ArrayList<>(currentCombination);
long seedIdx = 0;
for (Cross c : crosses) {
seedIdx = seedIdx * 10000L + c.getId();
}

seedMap.put((double) seedIdx, crosses);
}
return;
}

// Current detector layer
int layer = selectedLayers.get(layerIndex);
List<Cross> layerCrosses = sortedCrosses.get(binIdx).get(layer);

// ------------------------------------------------------------
// Try EVERY Cross in this layer
// ------------------------------------------------------------
for (Cross cross : layerCrosses) {
currentCombination.add(cross);
generateCrossCombinations(binIdx, selectedLayers, layerIndex + 1,currentCombination);
currentCombination.remove(currentCombination.size() - 1);
}
}

Expand Down Expand Up @@ -237,15 +311,14 @@ private int countType(List<Cross> cand, DetectorType dt) {
}
private double calcResi(double rho, double d0, double phi0, double xc, double yc) {
double r = Math.sqrt(xc*xc+yc*yc);
double par = 1. - ((r * r - d0 * d0) * rho * rho) / (2. * (1. + d0 * Math.abs(rho)));
double par = (2. * (1. + d0 * rho) - (r * r - d0 * d0) * rho * rho) / (2. * Math.abs(1. + d0 * rho));

double newPathLength = Math.abs(Math.acos(par) / rho);
int charge = (int) Math.signum(rho);
double alpha = -newPathLength * rho;

double x = d0 * charge * Math.sin(phi0) + (charge / Math.abs(rho))
* (Math.sin(phi0) - Math.cos(alpha) * Math.sin(phi0) - Math.sin(alpha) * Math.cos(phi0));
double y = -d0 * charge * Math.cos(phi0) - (charge / Math.abs(rho))
* (Math.cos(phi0) + Math.sin(alpha) * Math.sin(phi0) - Math.cos(alpha) * Math.cos(phi0));

double x = d0 * Math.sin(phi0) + 1./rho * (Math.sin(phi0) - Math.sin(phi0 + alpha));
double y = -d0 * Math.cos(phi0) - 1./rho * (Math.cos(phi0) - Math.cos(phi0 + alpha));

double res2 = ((x-xc)*(x-xc)+(y-yc)*(y-yc));
return Math.sqrt(res2);
Expand Down
35 changes: 6 additions & 29 deletions reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/Helix.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -178,15 +178,6 @@ public double radius() {
return 1. / C;
}

// (x,y) coordinates of the circle center
public double xCen() {
return (radius() - this.getDCA()) * Math.sin(this.getPhiAtDCA());
}

public double yCen() {
return (-radius() + this.getDCA()) * Math.cos(this.getPhiAtDCA());
}

// (x,y) coordinates of the dca
public double xDCA() {
return -this.getDCA() * Math.sin(this.getPhiAtDCA());
Expand All @@ -198,7 +189,7 @@ public double yDCA() {

public Point3D getVertex() {
return new Point3D(this.xDCA()+this.getXb(),this.yDCA()+this.getYb(),this.getZ0());
}
}

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

return new Vector3D(px,py,pz);
}

public double getArcLengthDCA(Point3D refpoint) {
//insure that the refpoint is on the helix
if (refpoint == null) {
return 0;
}
double refX = radius() * Math.cos(refpoint.toVector3D().phi());
double refY = radius() * Math.sin(refpoint.toVector3D().phi());
double arclen = arcLength(xCen(), yCen(), radius(), xCen(), yCen(), refX, refY);
return arclen;
}
}

// this method finds the arclength between 2 points in a circle
// this private method is used to get the pathlength from a point on the helical track to the distance of closest approach
Expand Down Expand Up @@ -272,16 +252,13 @@ public Point3D getPointAtRadius(double r) {
return new Point3D(x, y, z);
}


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

double alpha = -newPathLength * omega;

double x = d0 * charge * Math.sin(phi0) + (charge / Math.abs(omega))
* (Math.sin(phi0) - Math.cos(alpha) * Math.sin(phi0) - Math.sin(alpha) * Math.cos(phi0))+xb;
double y = -d0 * charge * Math.cos(phi0) - (charge / Math.abs(omega))
* (Math.cos(phi0) + Math.sin(alpha) * Math.sin(phi0) - Math.cos(alpha) * Math.cos(phi0))+yb;
double x = d0 * Math.sin(phi0) + 1./omega * (Math.sin(phi0) - Math.sin(phi0 + alpha));
double y = -d0 * Math.cos(phi0) - 1./omega * (Math.cos(phi0) - Math.cos(phi0 + alpha));
double z = z0 + newPathLength * tandip;

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

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

double alpha = newPathLength * omega;
Expand Down
Loading