fix some bugs in cvt reconstructions

common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/AStateVecs.java

@@ -509,7 +509,7 @@ public Helix getHelix(double xref, double yref) {
             double z0       = vec.z0 + vec.dz;
             //double omega    = turningSign * KFitter.polarity / R;  
             double omega    = Math.signum(vec.kappa) / R;  
-            double d0       = -vec.d_rho;
+            double d0       = vec.d_rho;
 
             Helix helix = new Helix(d0, phi0, omega, z0, tanDip, turningSign, bfield, xref, yref, units);
 

common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java

@@ -297,13 +297,13 @@ public void init(Helix helix, double[][] cov, double xref, double yref, double z
         initSV.kappa = initSV.alpha * helix.getOmega();
         // convert from the helix class representation to KF
         // kappa = alpha/omega = - q / pt
-        // drho  = - d0
+        // drho  =  d0
         // phi0  = (phi0 - 90)
 //        initSV.kappa = initSV.alpha * trk.getOmega();
 //        initSV.phi0  = trk.getPhi0()-Math.PI/2;
 //        initSV.dz    = trk.getZ();
 //        initSV.tanL  = trk.getTanL();
-//        initSV.d_rho = -trk.getD0();
+//        initSV.d_rho = trk.getD0();
 
         // recalculate helix parameters from vertex and momentum using chosen pivot
         initSV.updateHelix();

common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/trackrep/Helix.java

@@ -78,9 +78,9 @@ public Helix(double x0, double y0, double z0, double px0, double py0, double pz0
         double S = Math.sin(_phi0);
         double C = Math.cos(_phi0);
         if(Math.abs(S)>=Math.abs(C)) {
-            _d0 = -(x0-xb)/S;
+            _d0 = (x0-xb)/S;
         } else {
-            _d0 = (y0-yb)/C;
+            _d0 = -(y0-yb)/C;
         }
         _xb = xb;
         _yb = yb;
@@ -110,10 +110,10 @@ public void reset(double d0, double phi0, double omega, double z0, double tanL,
 
     public final void update() {
         setR(1./Math.abs(getOmega()));
-        _xd = -getD0()*getSinphi0()+_xb;
-        _yd =  getD0()*getCosphi0()+_yb;
-        _xc = -(-(double)KFitter.polarity*_turningSign*_R + _d0)*getSinphi0()+_xb; 
-        _yc =  (-(double)KFitter.polarity*_turningSign*_R + _d0)*getCosphi0()+_yb;
+        _xd = getD0()*getSinphi0()+_xb;
+        _yd =  -getD0()*getCosphi0()+_yb;
+        _xc = ((double)KFitter.polarity*_turningSign*_R + _d0)*getSinphi0()+_xb; 
+        _yc =  -((double)KFitter.polarity*_turningSign*_R + _d0)*getCosphi0()+_yb;
         _x  = getX(0);
         _y  = getY(0);
         _z  = getZ(0);

reconstruction/cvt/src/main/java/org/jlab/rec/cvt/banks/RecoBankWriter.java

@@ -519,7 +519,7 @@ public static DataBank fillTrackBank(DataEvent event, List<Track> trkcands, Stri
             bank.setFloat("phi0", i,  (float) helix.getPhiAtDCA());
             bank.setFloat("tandip", i,  (float) helix.getTanDip());
             bank.setFloat("z0", i,  (float) (helix.getZ0()/10.));
-            bank.setFloat("d0", i,  (float) (helix.getDCA()/10.));
+            bank.setFloat("d0", i,  (float) -(helix.getDCA()/10.));
             bank.setFloat("xb", i,  (float) (helix.getXb()/10.0));
             bank.setFloat("yb", i,  (float) (helix.getYb()/10.0));
             // this is the format of the covariance matrix for helical tracks
@@ -532,8 +532,8 @@ public static DataBank fillTrackBank(DataEvent event, List<Track> trkcands, Stri
             double[][] covmatrix = helix.getCovMatrix();
             if (covmatrix != null) {
                 bank.setFloat("cov_d02", i,  (float) covmatrix[0][0]/10/10 );
-                bank.setFloat("cov_d0phi0", i,  (float) covmatrix[0][1]/10 );
-                bank.setFloat("cov_d0rho", i,  (float) covmatrix[0][2] );
+                bank.setFloat("cov_d0phi0", i,  (float) -covmatrix[0][1]/10 );
+                bank.setFloat("cov_d0rho", i,  (float) -covmatrix[0][2] );
                 bank.setFloat("cov_phi02", i,  (float) covmatrix[1][1] );
                 bank.setFloat("cov_phi0rho", i,  (float) covmatrix[1][2]*10 );
                 bank.setFloat("cov_rho2", i,  (float) covmatrix[2][2]*10*10 );
@@ -592,7 +592,7 @@ public static DataBank fillUTrackBank(DataEvent event, List<Track> trkcands, Str
             bank.setFloat("phi0", i,  (float) helix.getPhiAtDCA());
             bank.setFloat("tandip", i,  (float) helix.getTanDip());
             bank.setFloat("z0", i,  (float) (helix.getZ0()/10.));
-            bank.setFloat("d0", i,  (float) (helix.getDCA()/10.));
+            bank.setFloat("d0", i,  (float) -(helix.getDCA()/10.));
             bank.setFloat("xb", i,  (float) (helix.getXb()/10.0));
             bank.setFloat("yb", i,  (float) (helix.getYb()/10.0));
             // this is the format of the covariance matrix for helical tracks
@@ -605,8 +605,8 @@ public static DataBank fillUTrackBank(DataEvent event, List<Track> trkcands, Str
             double[][] covmatrix = helix.getCovMatrix();
             if (covmatrix != null) {
                 bank.setFloat("cov_d02", i,  (float) covmatrix[0][0]/10/10 );
-                bank.setFloat("cov_d0phi0", i,  (float) covmatrix[0][1]/10 );
-                bank.setFloat("cov_d0rho", i,  (float) covmatrix[0][2] );
+                bank.setFloat("cov_d0phi0", i,  -(float) covmatrix[0][1]/10 );
+                bank.setFloat("cov_d0rho", i,  -(float) covmatrix[0][2] );
                 bank.setFloat("cov_phi02", i,  (float) covmatrix[1][1] );
                 bank.setFloat("cov_phi0rho", i,  (float) covmatrix[1][2]*10 );
                 bank.setFloat("cov_rho2", i,  (float) covmatrix[2][2]*10*10 );

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

@@ -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);
             }

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

@@ -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) 

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

@@ -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<>() );
             }
         }
@@ -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);
         }
     }
 
@@ -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);