test: dev-align3

common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/SVT/SVTConstants.java

@@ -411,7 +411,7 @@ else if( value == "tube" ) // offset by boxNum to reset row for CCDB table
                         SUPPORTRADIUS[region] = cp.getDouble(ccdbPath+"region/CuSupportInnerRadius", region); // radius to inner side of heatSinkRidge
 
                         for( int m = 0; m < NMODULES; m++ )
-                        {			
+		{			
 					switch( m ) 
 					{
 					case 0: // U = lower / inner

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

@@ -322,6 +322,13 @@ public HitOnTrack(int layer, int sector, StateVec sv, double tRes, double fRes,
             this.filteredResidual = fRes;
             this.smoothedResidual = sRes;
         }
+
+        @Override
+        public String toString() {
+            String str = String.format("HitOnTrack layer=%d sector=%d x=%.3f y=%.3f z=%.3f px=%.3f py=%.3f pz=%.3f residual=%.3f", 
+                                    this.layer, this.sector, this.x, this.y, this.z, this.px, this.py, this.pz, this.residual);
+            return str;
+        }
     }
 
     public void printConfig() {

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

@@ -47,7 +47,7 @@ public double[][] filterCovMat(double[] H, double[][] Carr, double V) {
         }
         return Cinv;
     }
-        
+
     public double[][] smoothCovMat(double[][] C_n_kp1, double[][] C_k, double[][] A , double[][] C_k_kp1) {
         double[][] At = null;
         try {

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

@@ -87,7 +87,7 @@ public void runFitterNoFiltFailSafe() {
     public void runFitter() {
         this.runFitter(sv, mv);
     }
-
+    
     public double getChi2() {
         return this.getChi2(0);
     }
@@ -186,9 +186,10 @@ public void runFitter(AStateVecs sv, AMeasVecs mv) {
             if(Double.isNaN(newchisq) ||
                sv.smoothed().get(0)==null ||
                sv.smoothed().get(0).kappa==0 || 
-               Double.isNaN(sv.smoothed().get(0).kappa)) {
-                this.setFitFailed = true; 
-                break; 
+               Double.isNaN(sv.smoothed().get(0).kappa) || 
+               Double.isNaN(sv.smoothed().get(0).dz)) {
+                this.setFitFailed = true;
+                break;
             }            
             // if chi2 improved and curvature is non-zero, save fit results but continue iterating
             else if(newchisq < this.chi2) {

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

@@ -7,7 +7,6 @@
 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;
@@ -18,11 +17,11 @@
  */
 public class StateVecs extends AStateVecs {
 
-
+     
     @Override
     public boolean getStateVecPosAtMeasSite(StateVec sv, AMeasVecs.MeasVec mv, Swim swim) {
         double[] swimPars = new double[7];
-        
+
         if(mv.surface==null) return false;
 
         int dir = (int) Math.signum(mv.k-sv.k);
@@ -41,27 +40,32 @@ public boolean getStateVecPosAtMeasSite(StateVec sv, AMeasVecs.MeasVec mv, Swim
         if(this.straight) {
 
             Point3D st   = new Point3D(sv.x, sv.y, sv.z); 
-            Vector3D stu = new Vector3D(sv.px,sv.py,sv.pz).asUnit();
+            Vector3D stu = new Vector3D(sv.px*dir,sv.py*dir,sv.pz*dir).asUnit();
 
             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();
                 sv.y = inters.y();
-                sv.z = inters.z();  
+                sv.z = inters.z();                   
                 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()); 
-                } 
+                double l1 = (-(st.x()*stu.x()+st.y()*stu.y())+delta)/(stu.x()*stu.x()+stu.y()*stu.y());
+                double l2 = (-(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) {
+                double l = l1;
+                if(l1>0 && l2<0) l = l1;
+                else if(l1<0 && l2>0) l = l2;
+                else if(l1>0 && l2>0) l = l1 < l2 ? l1 : l2;
+                else return false;
                 Point3D inters = new Point3D(st.x()+l*stu.x(),st.y()+l*stu.y(),st.z()+l*stu.z());
+                if(l<0) System.out.println("arghhhhhhhhhh "+ mv.layer + " " + mv.surface.getSector() + " " + st.toString() + " " + stu.toString() + " " + inters.toString());
                 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();
@@ -241,7 +245,7 @@ public double[][] Q(StateVec vec, AMeasVecs mv) {
         double cosEntranceAngle = this.getLocalDirAtMeasSite(vec, mv.measurements.get(vec.k));
 
         double p = Math.sqrt(vec.px*vec.px + vec.py*vec.py + vec.pz*vec.pz);
-        if(this.straight) p = 1;
+        if(this.straight) p = 100;
 
         // Highland-Lynch-Dahl formula
         double sctRMS = surf.getThetaMS(p, mass, cosEntranceAngle);

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

@@ -53,7 +53,12 @@ public void runFitter(AStateVecs sv, AMeasVecs mv) {
             // chi2
             double newchisq = this.calc_chi2(sv, mv); 
             // if curvature is 0, fit failed
-            if(Double.isNaN(newchisq) || sv.smoothed().get(0)==null) {
+            if(Double.isNaN(newchisq) || 
+               sv.smoothed().get(0)==null ||
+               Double.isNaN(sv.smoothed().get(0).x0) ||
+               Double.isNaN(sv.smoothed().get(0).z0) ||
+               Double.isNaN(sv.smoothed().get(0).tx) ||
+               Double.isNaN(sv.smoothed().get(0).tz)) {
                 this.setFitFailed = true;
                 break;
             }

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

@@ -108,7 +108,7 @@ public double[][] Q(StateVec vec, AMeasVecs mv) {
         Surface surf = mv.measurements.get(vec.k).surface;
         double cosEntranceAngle = this.getLocalDirAtMeasSite(vec, mv.measurements.get(vec.k));
 
-        double p = 1;
+        double p = 100;
 
         // Highland-Lynch-Dahl formula
         double sctRMS = surf.getThetaMS(p, mass, cosEntranceAngle);

reconstruction/cvt/pom.xml

@@ -18,8 +18,8 @@
 
     <dependencies>
 
-      <dependency>
-        <groupId>org.jlab.clas</groupId>
+    <dependency>
+      <groupId>org.jlab.clas</groupId>
         <artifactId>clas-jcsg</artifactId>
         <version>10.0.6-SNAPSHOT</version>
       </dependency>

reconstruction/cvt/src/main/java/org/jlab/rec/cvt/Constants.java

@@ -43,7 +43,7 @@ public static Constants getInstance() {
             }
             return instance;
     }
-    
+
     private static boolean ConstantsLoaded;
 
     // parameters configurable from yaml
@@ -326,8 +326,8 @@ public static double[][] scaleCovMat(double[][] matrix) {
         }
         return scaledMatrix;
     }
-    
-    
+
+
     private static final double D0     = 10;                    // 10 mm
     private static final double DPHI0  = Math.toRadians(10);    // 10 deg
     private static final double DRHO   = 0.01;                  // ~6-7 on kappa, i.e. 150 MeV on pt

reconstruction/cvt/src/main/java/org/jlab/rec/cvt/Geometry.java

@@ -1,5 +1,6 @@
 package org.jlab.rec.cvt;
 
+import cnuphys.magfield.MagneticFields;
 import java.util.ArrayList;
 import java.util.Arrays;
 import java.util.List;
@@ -74,6 +75,10 @@ public synchronized static void initialize(String variation, int run, IndexedTab
 
     private synchronized void load(String variation, int run, IndexedTable svtLorentz, IndexedTable bmtVoltage) {
 
+        if(Constants.getInstance().isCosmics) {
+            MagneticFields.getInstance().getSolenoid().setScaleFactor(1e-7);
+            MagneticFields.getInstance().getTorus().setScaleFactor(1e-7); 
+        }
         // Load target
         ConstantProvider providerTG = GeometryFactory.getConstants(DetectorType.TARGET, run, variation);
         this.zTarget = providerTG.getDouble("/geometry/target/position",0)*10;

reconstruction/cvt/src/main/java/org/jlab/rec/cvt/alignment/AlignmentBankReader.java

@@ -0,0 +1,156 @@
+package org.jlab.rec.cvt.alignment;
+
+import java.util.ArrayList;
+import java.util.List;
+import java.util.Collections;
+
+import org.jlab.detector.base.DetectorType;
+import org.jlab.geom.prim.Line3D;
+import org.jlab.geom.prim.Vector3D;
+import org.jlab.io.base.DataBank;
+import org.jlab.io.base.DataEvent;
+import org.jlab.rec.cvt.banks.RecoBankReader;
+import org.jlab.rec.cvt.bmt.BMTGeometry;
+import org.jlab.rec.cvt.cluster.Cluster;
+import org.jlab.rec.cvt.cross.Cross;
+import org.jlab.rec.cvt.hit.Hit;
+import org.jlab.rec.cvt.hit.Strip;
+import org.jlab.rec.cvt.track.StraightTrack;
+import org.jlab.rec.cvt.track.Track;
+import org.jlab.rec.cvt.trajectory.Helix;
+
+/**
+ *
+ * @author spaul
+ *
+ */
+public class AlignmentBankReader {
+
+    public List<StraightTrack> getCosmics(DataEvent event) {
+
+
+        var SVThits = RecoBankReader.readBSTHitBank(event, "BST::Hits");
+        var BMThits = RecoBankReader.readBMTHitBank(event, "BMT::Hits");
+        if(SVThits!= null) {
+            Collections.sort(SVThits);
+        }
+        if(BMThits!=null) {
+            for(Hit hit : BMThits) {
+                hit.getStrip().calcBMTStripParams(hit.getSector(), hit.getLayer(), swimmer);
+            }
+            Collections.sort(BMThits);
+        }
+
+        var _SVTclusters = RecoBankReader.readBSTClusterBank(event, SVThits, "BSTRec::Clusters");
+        var _BMTclusters = RecoBankReader.readBMTClusterBank(event, BMThits, "BMTRec::Clusters");
+
+
+        var _SVTcrosses = RecoBankReader.readBSTCrossBank(event, _SVTclusters, "BSTRec::Crosses");
+        var _BMTcrosses = RecoBankReader.readBMTCrossBank(event, _BMTclusters, "BMTRec::Crosses");
+        if(_SVTcrosses!=null) {
+            for(Cross cross : _SVTcrosses) {
+                cross.setCluster1(_SVTclusters.get(cross.getCluster1().getId()-1));
+                cross.setCluster2(_SVTclusters.get(cross.getCluster2().getId()-1)); 
+            }
+        }
+        if(_BMTcrosses!=null) {
+            for(Cross cross : _BMTcrosses) {
+                cross.setCluster1(_BMTclusters.get(cross.getCluster1().getId()-1));
+            }
+        }
+
+        List<StraightTrack> tracks = RecoBankReader.readCVTCosmicsBank(event, "CVTRec::Cosmics");
+        if(tracks == null) 
+            return null;
+
+        for(StraightTrack track : tracks) {
+
+            List<Cross> crosses = new ArrayList<>();
+            for(Cross c : track) {
+                if(_SVTcrosses!=null && c.getDetector()==DetectorType.BST) {
+                    for(Cross cross : _SVTcrosses) {
+                        if(c.getId() == cross.getId())
+                            crosses.add(cross);
+                    }
+                }
+                if(_BMTcrosses!=null && c.getDetector()==DetectorType.BMT) {
+                    for(Cross cross : _BMTcrosses) {
+                        if(c.getId() == cross.getId())
+                            crosses.add(cross);
+                    }
+                }
+            }
+            track.clear();
+            track.addAll(crosses);
+        }
+
+        return tracks;
+    }
+
+    public List<Track> getTracks(DataEvent event) {
+
+
+        var SVThits = RecoBankReader.readBSTHitBank(event, "BST::Hits");
+        var BMThits = RecoBankReader.readBMTHitBank(event, "BMT::Hits");
+        if(SVThits!= null) {
+            Collections.sort(SVThits);
+        }
+        if(BMThits!=null) {
+            for(Hit hit : BMThits) {
+                hit.getStrip().calcBMTStripParams(hit.getSector(), hit.getLayer(), swimmer);
+            }
+            Collections.sort(BMThits);
+        }
+
+        var _SVTclusters = RecoBankReader.readBSTClusterBank(event, SVThits, "BSTRec::Clusters");
+        var _BMTclusters = RecoBankReader.readBMTClusterBank(event, BMThits, "BMT::Clusters");
+
+
+        var _SVTcrosses = RecoBankReader.readBSTCrossBank(event, _SVTclusters, "BSTRec::Crosses");
+        var _BMTcrosses = RecoBankReader.readBMTCrossBank(event, _BMTclusters, "BMTRec::Crosses");
+        if(_SVTcrosses!=null) {
+            for(Cross cross : _SVTcrosses) {
+                cross.setCluster1(_SVTclusters.get(cross.getCluster1().getId()-1));
+                cross.setCluster2(_SVTclusters.get(cross.getCluster2().getId()-1)); 
+            }
+        }
+        if(_BMTcrosses!=null) {
+            for(Cross cross : _BMTcrosses) {
+                cross.setCluster1(_BMTclusters.get(cross.getCluster1().getId()-1));
+            }
+        }
+        double xb = 0.0;  /////// FIXME
+        double yb = 0.0;  /////// FIXME
+        _CVTseeds = RecoBankReader.readCVTSeedsBank(event, xb, yb, _SVTcrosses, _BMTcrosses, "CVT::Seeds");
+
+        List<Track> tracks = RecoBankReader.readCVTTracksBank(event, xb, yb, _CVTseeds, _SVTcrosses, _BMTcrosses, "CVTRec::Tracks");
+        if(tracks == null) 
+            return null;
+
+        for(Track track : tracks) {
+
+            List<Cross> crosses = new ArrayList<>();
+            for(Cross c : track) {
+                if(_SVTcrosses!=null && c.getDetector()==DetectorType.BST) {
+                    for(Cross cross : _SVTcrosses) {
+                        if(c.getId() == cross.getId())
+                            crosses.add(cross);
+                    }
+                }
+                if(_BMTcrosses!=null && c.getDetector()==DetectorType.BMT) {
+                    for(Cross cross : _BMTcrosses) {
+                        if(c.getId() == cross.getId())
+                            crosses.add(cross);
+                    }
+                }
+            }
+            track.clear();
+            track.addAll(crosses);
+
+        }
+
+        return tracks;
+    }
+
+
+}

reconstruction/cvt/src/main/java/org/jlab/rec/cvt/alignment/AlignmentBankWriter.java

@@ -0,0 +1,39 @@
+package org.jlab.rec.cvt.alignment;
+
+import java.util.List;
+
+import org.jlab.io.base.DataBank;
+import org.jlab.io.base.DataEvent;
+import org.ejml.simple.SimpleMatrix;
+
+/**
+ *
+ * @author spaul
+ *
+ */
+public class AlignmentBankWriter {
+
+	public void write_Matrix(DataEvent event, String matrixName, List<SimpleMatrix> matrices) {
+
+		//System.out.println("attempting to write matrices");
+		if(event == null)
+			System.out.println("event is null");
+		if(matrixName == null)
+			System.out.println("matrixName is null");
+		if(matrices == null)
+			System.out.println("matrix list is null");
+		DataBank bank = event.createBank("Align::" + matrixName, matrices.size());
+		for(int i = 0; i< matrices.size(); i++) {
+			bank.setShort("rows",i,(short) matrices.get(i).numRows());
+			bank.setShort("columns",i,(short) matrices.get(i).numCols());
+			for(int j = 0; j<matrices.get(i).numRows(); j++) {
+				for(int k = 0; k<matrices.get(i).numCols(); k++) {
+                                    if(Double.isNaN(matrices.get(i).get(j,k))) 
+                                        System.out.println(event.getBank("RUN::config").getInt("event", 0) + " " + matrixName + " " + matrices.get(i).get(j,k));
+					bank.setFloat("element_"+j+"_"+k, i, (float) matrices.get(i).get(j,k));
+				}
+			}
+		}
+		event.appendBank(bank);
+	}	
+}