diff --git a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/SVT/SVTConstants.java b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/SVT/SVTConstants.java index dd7c9275a7..336ab0917f 100644 --- a/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/SVT/SVTConstants.java +++ b/common-tools/clas-jcsg/src/main/java/org/jlab/detector/geant4/v2/SVT/SVTConstants.java @@ -400,7 +400,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 diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/AKFitter.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/AKFitter.java index a940a6f9fe..a495bb057c 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/AKFitter.java +++ b/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() { diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/KFCovMatOps.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/KFCovMatOps.java index 5e7b4e84eb..9e90af583c 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/KFCovMatOps.java +++ b/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 { diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/KFitter.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/KFitter.java index f50f453521..abc291acdb 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/KFitter.java +++ b/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) { diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java index 067bcf5c5a..b95940ff5b 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/helical/StateVecs.java +++ b/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,7 +40,7 @@ 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); @@ -49,7 +48,7 @@ public boolean getStateVecPosAtMeasSite(StateVec sv, AMeasVecs.MeasVec mv, Swim 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) { @@ -57,11 +56,16 @@ else if(mv.surface.cylinder!=null) { 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); diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/KFitter.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/KFitter.java index 4c212a492d..cd95f6c36f 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/KFitter.java +++ b/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; } diff --git a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java b/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java index d811b23f9b..efa42dfb7e 100644 --- a/common-tools/clas-tracking/src/main/java/org/jlab/clas/tracking/kalmanfilter/straight/StateVecs.java +++ b/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); diff --git a/reconstruction/cvt/pom.xml b/reconstruction/cvt/pom.xml index 8fa13b5e34..bb3522184e 100644 --- a/reconstruction/cvt/pom.xml +++ b/reconstruction/cvt/pom.xml @@ -18,8 +18,8 @@ - - org.jlab.clas + + org.jlab.clas clas-jcsg 11.1.0-SNAPSHOT diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/Constants.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/Constants.java index 51303f0636..d77a2d71fc 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/Constants.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/Constants.java @@ -41,7 +41,7 @@ public static Constants getInstance() { } return instance; } - + private static boolean ConstantsLoaded; // parameters configurable from yaml @@ -316,8 +316,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 diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/alignment/AlignmentBankReader.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/alignment/AlignmentBankReader.java new file mode 100644 index 0000000000..69c9e8fae2 --- /dev/null +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/alignment/AlignmentBankReader.java @@ -0,0 +1,162 @@ +package org.jlab.rec.cvt.alignment; + +import java.util.ArrayList; +import java.util.List; +import java.util.Collections; +import java.util.Map; + +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 { + + private Map _SVTcrosses; + private Map _SVTclusters; + private Map _BMTcrosses; + private Map _BMTclusters; + + public List getCosmics(DataEvent event) { + + + List SVThits = RecoBankReader.readBSTHitBank(event, "BSTRec::Hits"); + List BMThits = RecoBankReader.readBMTHitBank(event, "BMTRec::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); + } + + _SVTclusters = RecoBankReader.readBSTClusterBank(event, SVThits, "BSTRec::Clusters"); + _BMTclusters = RecoBankReader.readBMTClusterBank(event, BMThits, "BMTRec::Clusters"); + + + _SVTcrosses = RecoBankReader.readBSTCrossBank(event, _SVTclusters, "BSTRec::Crosses"); + _BMTcrosses = RecoBankReader.readBMTCrossBank(event, _BMTclusters, "BMTRec::Crosses"); + if(_SVTcrosses!=null) { + for(Cross cross : _SVTcrosses.values()) { + cross.setCluster1(_SVTclusters.get(cross.getCluster1().getId()-1)); + cross.setCluster2(_SVTclusters.get(cross.getCluster2().getId()-1)); + } + } + if(_BMTcrosses!=null) { + for(Cross cross : _BMTcrosses.values()) { + cross.setCluster1(_BMTclusters.get(cross.getCluster1().getId()-1)); + } + } + + List tracks = RecoBankReader.readCVTCosmicsBank(event, "CVTRec::Cosmics"); + if(tracks == null) + return null; + + for(StraightTrack track : tracks) { + + List crosses = new ArrayList<>(); + for(Cross c : track) { + if(_SVTcrosses!=null && c.getDetector()==DetectorType.BST) { + for(Cross cross : _SVTcrosses.values()) { + if(c.getId() == cross.getId()) + crosses.add(cross); + } + } + if(_BMTcrosses!=null && c.getDetector()==DetectorType.BMT) { + for(Cross cross : _BMTcrosses.values()) { + if(c.getId() == cross.getId()) + crosses.add(cross); + } + } + } + track.clear(); + track.addAll(crosses); + } + + return tracks; + } + + public List getTracks(DataEvent event) { + + + List SVThits = RecoBankReader.readBSTHitBank(event, "BSTRec::Hits"); + List BMThits = RecoBankReader.readBMTHitBank(event, "BMTRec::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); + } + + _SVTclusters = RecoBankReader.readBSTClusterBank(event, SVThits, "BSTRec::Clusters"); + _BMTclusters = RecoBankReader.readBMTClusterBank(event, BMThits, "BMT::Clusters"); + + + _SVTcrosses = RecoBankReader.readBSTCrossBank(event, _SVTclusters, "BSTRec::Crosses"); + _BMTcrosses = RecoBankReader.readBMTCrossBank(event, _BMTclusters, "BMTRec::Crosses"); + if(_SVTcrosses!=null) { + for(Cross cross : _SVTcrosses.values()) { + cross.setCluster1(_SVTclusters.get(cross.getCluster1().getId()-1)); + cross.setCluster2(_SVTclusters.get(cross.getCluster2().getId()-1)); + } + } + if(_BMTcrosses!=null) { + for(Cross cross : _BMTcrosses.values()) { + 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 tracks = RecoBankReader.readCVTTracksBank(event, xb, yb, _CVTseeds, _SVTcrosses, _BMTcrosses, "CVTRec::Tracks"); + if(tracks == null) + return null; + + for(Track track : tracks) { + + List crosses = new ArrayList<>(); + for(Cross c : track) { + if(_SVTcrosses!=null && c.getDetector()==DetectorType.BST) { + for(Cross cross : _SVTcrosses.values()) { + if(c.getId() == cross.getId()) + crosses.add(cross); + } + } + if(_BMTcrosses!=null && c.getDetector()==DetectorType.BMT) { + for(Cross cross : _BMTcrosses.values()) { + if(c.getId() == cross.getId()) + crosses.add(cross); + } + } + } + track.clear(); + track.addAll(crosses); + + } + + return tracks; + } + + +} diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/alignment/AlignmentBankWriter.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/alignment/AlignmentBankWriter.java new file mode 100644 index 0000000000..428d0b81b1 --- /dev/null +++ b/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 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 tracks; + if(isCosmics) { + tracks = reader.getCosmics(event); + } + else { + tracks = reader.getTracks(event); + } + if (tracks==null || tracks.size() > 2) return true; + + + List Is = new ArrayList<>(); + List As = new ArrayList<>(); + List Bs = new ArrayList<>(); + List Vs = new ArrayList<>(); + List ms = new ArrayList<>(); + List cs = new ArrayList<>(); + List qs = new ArrayList<>(); + List trackIDs = new ArrayList<>(); + + + tracksLoop : for (Trajectory track : tracks) { + + if(Math.abs(getDoca(track))>maxDocaCut) + continue; + + int nCrossSVT = 0, nCrossBMT = 0; + int countBMTZ = 0, countBMTC = 0; + for(Cross c : track) { + if (c.getDetector() == DetectorType.BST && !isBMTonly) { + nCrossSVT++; + gCountSVT++; + } + if (c.getDetector() == DetectorType.BMT && !isSVTonly) { + if (c.getCluster1().getType() != BMTType.C || !skipBMTC) { + nCrossBMT++; + } + if (c.getCluster1().getType() == BMTType.C) { + countBMTC++; + gCountBMTC++; + } + if (c.getCluster1().getType() == BMTType.Z) { + gCountBMTZ++; + countBMTZ++; + } + //System.out.println(c.getSector()+" "+c.getRegion() + " " + c.getCluster1().getCentroid()+" " + c.getId()); + } + if (nCrossBMT > 12) { + System.out.println("Too many BMT crosses!"); + return false; + } + } + + + if(nCrossSVT*20.001 && !curvedTracks) { + continue; + } + } + + int paramsFromBeamspot = (isCosmics || ! includeBeamspot ? 0:1); + int cols = nAlignVars*((svtTopBottomSep ? 2*nCrossSVT : nCrossSVT) + nCrossBMT + paramsFromBeamspot); + int rows = 2*nCrossSVT+nCrossBMT + paramsFromBeamspot; + + SimpleMatrix A = new SimpleMatrix(rows, cols);//not sure why there aren't 6 columns + SimpleMatrix B = new SimpleMatrix(rows, 4); + SimpleMatrix V = new SimpleMatrix(rows,rows); + SimpleMatrix m = new SimpleMatrix(rows,1); + SimpleMatrix c = new SimpleMatrix(rows,1); + SimpleMatrix I = new SimpleMatrix(rows,1); + SimpleMatrix q = new SimpleMatrix(4, 1); //track parameters, for plotting kinematic dependence. Not used in KFA. + + if (track.getHelix() == null) { + track.setHelix(createHelixFromRay(track.getRay())); + } + + q.set(0, 0, track.getHelix().getDCA()); + q.set(1, 0, track.getHelix().getPhiAtDCA()); + q.set(2, 0, track.getHelix().getZ0()); + q.set(3, 0, track.getHelix().getTanDip()); + + if(debug) { + System.out.println("track parameters"); + q.print(); + } + + int i = 0; + if(!curvedTracks) { + for (Cross cross : track) { + if (useNewFillMatrices) { + if (cross.getDetector() == DetectorType.BST) { + Cluster cl1 = cross.getCluster1(); + boolean ok = fillMatricesNew(i, ray, cl1, A, B, V, m, c, I, debug, false); + i++; + if (!ok) { //reject track if there's a cluster with really bad values. + if (debug) { + System.out.println("rejecting track due to problem in an SVT layer"); + } + continue tracksLoop; + } + Cluster cl2 = cross.getCluster2(); + ok = fillMatricesNew(i, ray, cl2, A, B, V, m, c, I, debug, false); + i++; + if (!ok) { //reject track if there's a cluster with really bad values. + if (debug) { + System.out.println("rejecting track due to problem in an SVT layer"); + } + continue tracksLoop; + } + } else { + Cluster cl1 = cross.getCluster1(); + boolean ok = true; + if (cl1.getType() == BMTType.Z || !skipBMTC) { + ok = fillMatricesNew(i, ray, cl1, A, B, V, m, c, I, this.debug, false); + } + i++; + if (!ok) { //reject track if there's a cluster with really bad values. + if (debug) { + System.out.println("rejecting track due to problem in a BMT" + cl1.getType().name() + " layer"); + } + continue tracksLoop; + } + //} + } + continue; + } + + } + if (!isCosmics && includeBeamspot) { + + //pseudo cluster for the beamspot + Cluster cl1 = new Cluster(null, null, 0, 0, 0); + cl1.setLine(new Line3D(track.getHelix().getXb(), track.getHelix().getYb(), -100, track.getHelix().getXb(), track.getHelix().getYb(), 100)); + + Vector3D n = ray.getDirVec(); + Vector3D l = new Vector3D(0, 0, 1); + cl1.setN(n); + cl1.setL(l); + cl1.setS(n.cross(l)); + cl1.setResolution(0.6); + + fillMatricesNew(i, ray, cl1, A, B, V, m, c, I, this.debug, true); + + } + } + else { + Helix helix = track.getHelix(); + //curved tracks + for (Cross cross : track) { + if (useNewFillMatrices) { + if (cross.getDetector() == DetectorType.BST) { + Cluster cl1 = cross.getCluster1(); + boolean ok = fillMatricesNew(i, helix, cl1, A, B, V, m, c, I, debug, false); + i++; + if (!ok) { //reject track if there's a cluster with really bad values. + if (debug) { + System.out.println("rejecting track due to problem in an SVT layer"); + } + continue tracksLoop; + } + Cluster cl2 = cross.getCluster2(); + ok = fillMatricesNew(i, helix, cl2, A, B, V, m, c, I, debug, false); + i++; + if (!ok) { //reject track if there's a cluster with really bad values. + if (debug) { + System.out.println("rejecting track due to problem in an SVT layer"); + } + continue tracksLoop; + } + } else { + Cluster cl1 = cross.getCluster1(); + boolean ok = true; + if (cl1.getType() == BMTType.Z || !skipBMTC) { + ok = fillMatricesNew(i, helix, cl1, A, B, V, m, c, I, this.debug, false); + } + i++; + if (!ok) { //reject track if there's a cluster with really bad values. + if (debug) { + System.out.println("rejecting track due to problem in a BMT" + cl1.getType().name() + " layer"); + } + continue tracksLoop; + } + //} + } + } + + } + if(!isCosmics && includeBeamspot) { + + //pseudo cluster for the beamspot + Cluster cl1 = new Cluster(null, null, 0, 0, 0); + cl1.setLine(new Line3D(track.getHelix().getXb(), track.getHelix().getYb(), -100, track.getHelix().getXb(), track.getHelix().getYb(), 100)); + + Vector3D n = ray.getDirVec(); + Vector3D l = new Vector3D(0,0,1); + cl1.setN(n); + cl1.setL(l); + cl1.setS(n.cross(l)); + cl1.setResolution(0.6); + + fillMatricesNew(i, helix, cl1, A,B,V,m,c,I, this.debug, true); + + + } + } + As.add(A); + Bs.add(B); + Vs.add(V); + ms.add(m); + cs.add(c); + Is.add(I); + qs.add(q); + trackIDs.add(track.getId()); + } + + //only include events that have tracks that will be used in alignment + if(!As.isEmpty()) { + AlignmentBankWriter writer = new AlignmentBankWriter(); + writer.write_Matrix(event, "I", Is); + writer.write_Matrix(event, "A", As); + writer.write_Matrix(event, "B", Bs); + writer.write_Matrix(event, "V", Vs); + writer.write_Matrix(event, "m", ms); + writer.write_Matrix(event, "c", cs); + writer.write_Matrix(event, "q", qs); + fillMisc(event,runNum,eventNum,trackIDs,As,Bs,Vs,ms,cs,Is); + } + return true; + + } + private Helix createHelixFromRay(Ray ray) { + Vector3D u = ray.getDirVec(); + Vector3D xref = ray.getRefPoint().toVector3D(); + double phi = Math.atan2(u.y(),u.x()); + Vector3D uT = new Vector3D(Math.cos(phi), Math.sin(phi),0); + Vector3D mscphi = new Vector3D(-Math.sin(phi), Math.cos(phi),0); + double cosdip = Math.hypot(u.x(), u.y()); + double d = mscphi.dot(xref); + double curvature = 0; + double Z0 = xref.z()-u.z()*xref.dot(uT)/u.dot(uT); + double tandip = u.z()/Math.hypot(u.x(), u.y()); + return new Helix(d, phi, curvature, Z0, tandip, 0,0); + } + + + private Ray getRay(Helix h) { + + double d = h.getDCA(); + double z = h.getZ0(); + double phi = h.getPhiAtDCA(); + double td = h.getTanDip(); + double cd = 1/Math.hypot(td, 1); + double sd = td*cd; + double xb = h.getXb(); + double yb = h.getYb(); + //Vector3D u = new Vector3D(-cd*Math.sin(phi), cd*Math.cos(phi), sd); + //Point3D x = new Point3D(d*Math.cos(phi),d*Math.sin(phi), z); + Vector3D u = new Vector3D(cd*Math.cos(phi), cd*Math.sin(phi), sd); + + + Point3D x = new Point3D(-d*Math.sin(phi)+xb,d*Math.cos(phi)+yb, z); + //Point3D x = new Point3D(-d*Math.sin(phi),d*Math.cos(phi), z); + + //System.out.println("xb yb from db" + xb + yb); + //Point3D x = new Point3D(-d*Math.sin(phi),d*Math.cos(phi), z); + //if(u.y() <0) + // u = u.multiply(-1); + //x = x.toVector3D().add(u.multiply(-x.y()/u.y())).toPoint3D(); + Ray ray = new Ray(x, u); + //System.out.println("doca " + d); + //System.out.println("td " + td); + + return ray; + } + + + + private double getDoca(Trajectory track) { + if(track instanceof StraightTrack) { + Ray ray = track.getRay(); + double intercept = ray.getYXInterc(); + double slope = ray.getYXSlope(); + return Math.abs(intercept)/Math.hypot(1, slope); + } else return track.getHelix().getDCA(); + } + + private void fillMisc(DataEvent event, int runNum, int eventNum, List trackIDs, + List As, List Bs, List Vs, + List ms, List cs, List is) { + DataBank bank = event.createBank("Align::misc", trackIDs.size()); + for(int i = 0; i spMax) { //this can only happen if the angle between the track and the normal is small + //System.out.println("rejecting track"); + return false; + } + int index = nAlignables - 1; + + //System.out.println("i = " + i + "; rows = " + A.getRowDimension() + "; cols = " + + A.getColumnDimension()); + Vector3D dmdr = sp.cross(xref).sub(sp.cross(u).multiply(n.dot(xref.clone().sub(e)) / udotn)); + if (orderTx >= 0) { + A.set(i, i * nAlignVars + orderTx, sp.x()); + } + if (orderTy >= 0) { + A.set(i, i * nAlignVars + orderTy, sp.y()); + } + if (orderTz >= 0) { + A.set(i, i * nAlignVars + orderTz, sp.z()); + } + if (orderRx >= 0) { + A.set(i, i * nAlignVars + orderRx, -dmdr.x()); + } + if (orderRy >= 0) { + A.set(i, i * nAlignVars + orderRy, -dmdr.y()); + } + if (orderRz >= 0) { + A.set(i, i * nAlignVars + orderRz, -dmdr.z()); + } + + I.set(i, 0, index); + + Vector3D dmdu = sp.multiply(e.clone().sub(xref).dot(n) / udotn); + if (!this.useDocaPhiZTandip) { + B.set(i, 0, sp.x()); + B.set(i, 1, sp.z()); + B.set(i, 2, dmdu.x()); + B.set(i, 3, dmdu.z()); + } else { + + double phi = Math.atan2(u.y(), u.x()); + Vector3D csphi = new Vector3D(Math.cos(phi), Math.sin(phi), 0); + Vector3D mscphi = new Vector3D(-Math.sin(phi), Math.cos(phi), 0); + double cosdip = Math.hypot(u.x(), u.y()); + double d = mscphi.dot(xref); + B.set(i, 0, -sp.dot(mscphi)); + B.set(i, 1, -sp.dot(mscphi) * n.dot(e.clone().sub(xref)) * cosdip / udotn + sp.dot(csphi) * d); + B.set(i, 2, -sp.z()); + B.set(i, 3, -sp.z() * n.dot(e.clone().sub(xref)) / udotn); + + } + //dm.set(i,0, s.dot(e.minus(extrap))); + + double ci = s.dot(extrap); + double mi = s.dot(e); + + c.set(i, 0, ci); + m.set(i, 0, mi); + + return true; + } + + + private boolean fillMatricesNew(int i, Helix helix, Cluster cl, SimpleMatrix A, SimpleMatrix B, SimpleMatrix V, + SimpleMatrix m, SimpleMatrix c, SimpleMatrix I, boolean debug, boolean isBeamspot) { + Vector3D u= helix.getTrackDirectionAtRadius(cl.getRadius()); + Point3D xref = helix.getPointAtRadius(cl.getRadius()); + Ray ray = new Ray(xref, u); + return fillMatricesNew(i, ray, cl, A, B, V, m, c, I, debug, isBeamspot); + } + + /** + * generic method that uses any type of cluster. + * @param i + * @param ray + * @param cl + * @param A + * @param B + * @param V + * @param m + * @param c + * @param I + * @param string + * @return + */ + private boolean fillMatricesNew(int i, Ray ray, Cluster cl, SimpleMatrix A, SimpleMatrix B, SimpleMatrix V, + SimpleMatrix m, SimpleMatrix c, SimpleMatrix I, boolean debug, boolean isBeamspot) { + int layer = cl.getLayer(); + int sector = cl.getSector(); + + Vector3D l; + Vector3D s; + Vector3D n; + + DetectorType detector = cl.getDetector(); + BMTType bmtType = cl.getType(); + + if (debug) { + System.out.println("\n\nNew method " + detector + " layer " + layer + " sector " + sector); + + } + Vector3D xref = ray.getRefPoint().toVector3D(); + Vector3D u = ray.getDirVec(); + + Vector3D e = null; + Vector3D extrap = null; + if (detector == DetectorType.BST || (detector == DetectorType.BMT && bmtType == BMTType.Z) || isBeamspot) { + l = cl.getL(); + s = cl.getS(); + n = cl.getN(); + + e = cl.getLine().midpoint().toVector3D(); + + double udotn = u.dot(n); + extrap = xref.clone().add(u.multiply(n.dot(e.clone().sub(xref)) / udotn)); + } + else { // BMTC + Vector3D a = cl.getArc().normal(); + + if (debug) { + System.out.println("a: " + a); + } + Vector3D cc = cl.getArc().center().toVector3D(); + Vector3D uT = perp(u, a); + + Vector3D tmp1 = perp(xref.clone().sub(cc), a); + + Vector3D endpoint = cl.getArc().origin().toVector3D(); + + double R = perp(endpoint.clone().sub(cc), a).mag(); + if (debug) { + System.out.println("center: " + cc.toStringBrief()); + System.out.println("R: " + R); + } + double AA = uT.mag2(); + + double BB = 2 * tmp1.dot(uT); + double CC = tmp1.mag2() - R * R; + double lambda_plus = (-BB + Math.sqrt(BB * BB - 4 * AA * CC)) / (2 * AA); + double lambda_minus = (-BB - Math.sqrt(BB * BB - 4 * AA * CC)) / (2 * AA); + Vector3D extrap_plus = xref.clone().add(u.multiply(lambda_plus)); + Vector3D extrap_minus = xref.clone().add(u.multiply(lambda_minus)); + + if (debug) { + System.out.println("extrap is on cylinder: this should be zero: " + (perp(extrap_plus.clone().sub(cc), a).mag() - R)); + } + + //choose the extrapolated point that is closer in z to the measured cluster. + if (Math.abs(extrap_plus.clone().sub(cc).z()) < Math.abs(extrap_minus.clone().sub(cc).z())) { + extrap = extrap_plus; + } else { + extrap = extrap_minus; + } + e = extrap.clone().add(endpoint.clone().sub(extrap).projection(a)); + s = a; + n = perp(extrap.clone().sub(cc), a).asUnit(); + l = s.cross(n); + } + + n = n.sub(l.multiply(n.dot(l))).asUnit(); + s = s.sub(l.multiply(s.dot(l))).asUnit(); + + double udotn = u.dot(n); + if (Math.abs(udotn) < minCosIncident) { + if (debug) { + System.out.println("rejecting track: abs(udotn)<" + minCosIncident); + System.out.println("u = " + u.toString()); + System.out.println("n = " + n.toString()); + } + return false; + } + double sdotu = s.dot(u); + + if (debug) { + System.out.println("e: " + e.toString()); + } + + if (detector == DetectorType.BMT && bmtType == BMTType.Z && debug) { + Vector3D diff = xref.clone(); + double check = l.cross(u).dot(diff); + System.out.println("distance between track and strip, phi,r: " + check + " " + u.phi() + " " + e.mag()); + } + + double resolution = cl.getResolution(); + + V.set(i, i, Math.pow(resolution, 2)); + if (debug) { + System.out.println("resolution " + resolution); + } + + Vector3D sp = s.clone().sub(n.multiply(sdotu / udotn)); + if (sp.mag() > spMax) { //this can only happen if the angle between the track and the normal is small + if (debug) { + System.out.println("rejecting track: sp.magnitude() > " + spMax); + } + return false; + } + + Vector3D dmdr = sp.cross(xref).sub(sp.cross(u).multiply(n.dot(xref.clone().sub(e)) / udotn)); + + if (orderTx >= 0) { + A.set(i, (svtTopBottomSep ? i : i / 2) * nAlignVars + orderTx, sp.x()); + } + if (orderTy >= 0) { + A.set(i, (svtTopBottomSep ? i : i / 2) * nAlignVars + orderTy, sp.y()); + } + if (orderTz >= 0) { + A.set(i, (svtTopBottomSep ? i : i / 2) * nAlignVars + orderTz, sp.z()); + } + if (orderRx >= 0) { + A.set(i, (svtTopBottomSep ? i : i / 2) * nAlignVars + orderRx, -dmdr.x()); + } + if (orderRy >= 0) { + A.set(i, (svtTopBottomSep ? i : i / 2) * nAlignVars + orderRy, -dmdr.y()); + } + if (orderRz >= 0) { + A.set(i, (svtTopBottomSep ? i : i / 2) * nAlignVars + orderRz, -dmdr.z()); + } + + if (detector == DetectorType.BST) { + I.set(i, 0, getIndexSVT(layer - 1, sector - 1)); + } else if (detector == DetectorType.BMT) { + I.set(i, 0, getIndexBMT(layer - 1, sector - 1)); + } else { + I.set(i, 0, INDEXBEAMLINE); + } + Vector3D dmdu = sp.multiply(e.clone().sub(xref).dot(n) / udotn); + if (!this.useDocaPhiZTandip) { + B.set(i, 0, -sp.x()); + B.set(i, 1, -sp.z()); + B.set(i, 2, -dmdu.x()); + B.set(i, 3, -dmdu.z()); + } else { + + double phi = Math.atan2(u.y(), u.x()); + /*if(detector == DetectorType.BMT && bmtType==BMTType.C) { + System.out.println("phi is " + phi); + }*/ + Vector3D csphi = new Vector3D(Math.cos(phi), Math.sin(phi), 0); + Vector3D mscphi = new Vector3D(-Math.sin(phi), Math.cos(phi), 0); + double cosdip = Math.hypot(u.x(), u.y()); + double d = mscphi.dot(xref); + B.set(i, 0, -sp.dot(mscphi)); + B.set(i, 1, -sp.dot(mscphi) * n.dot(e.clone().sub(xref)) * cosdip / udotn + sp.dot(csphi) * d); + B.set(i, 2, -sp.z()); + B.set(i, 3, -sp.z() * n.dot(e.clone().sub(xref)) / udotn); + + } + //dm.set(i,0, s.dot(e.minus(extrap))); + + double ci = s.dot(extrap); + double mi = s.dot(e); + + c.set(i, 0, ci); + m.set(i, 0, mi); + if (debug) { + System.out.println("n.(e-xref): " + n.dot(e.clone().sub(xref))); + System.out.println("u.n: " + udotn); + System.out.println("s: " + s.toString()); + System.out.println("sp: " + sp.toString()); + System.out.println("extrap: " + extrap.toString()); + System.out.println("n: " + n.toString()); + System.out.println("e: " + e.toString()); + System.out.println("xref: " + xref.toString()); + System.out.println("u: " + u.toString()); + System.out.println("m: " + mi); + System.out.println("c: " + ci); + } + if (Math.abs(ci - mi) > maxResidualCutSVT && detector == DetectorType.BST + || Math.abs(ci - mi) > maxResidualCutBMTZ && detector == DetectorType.BMT && bmtType == BMTType.Z + || Math.abs(ci - mi) > maxResidualCutBMTC && detector == DetectorType.BMT && bmtType == BMTType.C) { + if (debug) { + System.out.println("rejecting track: Math.abs(ci-mi)>maxResidualCut"); + } + return false; + } + return true; + } + + + Vector3D perp(Vector3D v, Vector3D a) { + return v.clone().sub(a.multiply(v.dot(a))); + } + + private int getIndexBMT(int layer, int sector) { + if (layer < 0 || sector < 0) { + return -1; + } + return NSVTSENSORS+layer*3+sector; + } + + private int getIndexSVT(int layer, int sector){ + int index = -1; + int region = layer/2; + switch (region) { + case 0: + index = sector; + break; + case 1: + index = SVTGeometry.NSECTORS[0] + sector; + break; + case 2: + index = SVTGeometry.NSECTORS[0] + + SVTGeometry.NSECTORS[2] + sector; + break; + default: + break; + } + if(svtTopBottomSep && layer%2==1) { + index += NSVTSENSORS/2; + } + return index; + + } + + + @Override + public boolean init() { + if(this.getEngineConfiguration() == null || "null".equals(this.getEngineConfiguration())) { + return true; //prevents init from being run twice. + } + + if (this.getEngineConfigString("svtOnly")!=null) { + this.isSVTonly= Boolean.valueOf(this.getEngineConfigString("svtOnly")); + } + System.out.println("["+this.getName()+"] align SVT only set to " + this.isSVTonly); + + if (this.getEngineConfigString("bmtOnly")!=null) { + this.isBMTonly= Boolean.valueOf(this.getEngineConfigString("bmtOnly")); + } + System.out.println("["+this.getName()+"] align BMT only set to " + this.isBMTonly); + + if (this.getEngineConfigString("skipBMTC")!=null) { + this.skipBMTC= Boolean.valueOf(this.getEngineConfigString("skipBMTC")); + } + System.out.println("["+this.getName()+"] skip BMTC set to " + this.skipBMTC); + + if (this.getEngineConfigString("svtAlignTopBottomSeparately")!=null) { + this.svtTopBottomSep= Boolean.valueOf(this.getEngineConfigString("svtAlignTopBottomSeparately")); + } + System.out.println("["+this.getName()+"] align SVT top-bottom separately set to " + this.svtTopBottomSep); + + if (this.getEngineConfigString("alignVariables")!=null) { + this.alignVars = this.getEngineConfigString("alignVariables"); + } + System.out.println("["+this.getName()+"] align variables set to " + this.alignVars); + + if (this.getEngineConfigString("cosmics")!=null) { + this.isCosmics= Boolean.valueOf(this.getEngineConfigString("cosmics")); + } + System.out.println("["+this.getName()+"] use cosmics bank set to " + this.isCosmics); + + if(this.getEngineConfigString("maxDocaCut") != null) { + this.maxDocaCut = Double.parseDouble(this.getEngineConfigString("maxDocaCut")); + } + if(isCosmics) this.maxDocaCut = Double.MAX_VALUE; + System.out.println("["+this.getName()+"] track DOCA cut set to " + this.maxDocaCut + " mm"); + + + if (this.getEngineConfigString("maxResidual")!=null) { + this.maxResidualCutBMTZ = this.maxResidualCutBMTC = this.maxResidualCutSVT = Double.valueOf(this.getEngineConfigString("maxResidual")); + } + if (this.getEngineConfigString("maxResidualSVT")!=null) { + this.maxResidualCutSVT = Double.valueOf(this.getEngineConfigString("maxResidualSVT")); + } + if (this.getEngineConfigString("maxResidualBMTZ")!=null) { + this.maxResidualCutBMTZ = Double.valueOf(this.getEngineConfigString("maxResidualBMTZ")); + } + if (this.getEngineConfigString("maxResidualBMTC")!=null) { + this.maxResidualCutBMTC = Double.valueOf(this.getEngineConfigString("maxResidualBMTC")); + } + System.out.println("["+this.getName()+"] max residual for SVT clusters set to " + this.maxResidualCutSVT + " mm"); + System.out.println("["+this.getName()+"] max residual for BMTC clusters set to " + this.maxResidualCutBMTC + " mm"); + System.out.println("["+this.getName()+"] max residual for BMTZ clusters set to " + this.maxResidualCutBMTZ + " mm"); + + if (this.getEngineConfigString("minClustersSVT")!=null) { + this.minClustersSVT = Integer.valueOf(this.getEngineConfigString("minClustersSVT")); + } + if (this.getEngineConfigString("minClustersBMTZ")!=null) { + this.minClustersBMTZ = Integer.valueOf(this.getEngineConfigString("minClustersBMTZ")); + } + if (this.getEngineConfigString("minClustersBMTC")!=null) { + this.minClustersBMTC = Integer.valueOf(this.getEngineConfigString("minClustersBMTC")); + } + System.out.println("["+this.getName()+"] min number of SVT clusters set to " + this.minClustersSVT); + System.out.println("["+this.getName()+"] min number of BMTC clusters set to " + this.minClustersBMTC); + System.out.println("["+this.getName()+"] min number of BMTZ clusters set to " + this.minClustersBMTZ); + + if (this.getEngineConfigString("useBeamspot")!=null) { + this.includeBeamspot = Boolean.valueOf(this.getEngineConfigString("useBeamspot")); + } + System.out.println("["+this.getName()+"] treat beamspot as an additional measurement set to " + this.includeBeamspot); + + if (this.getEngineConfigString("debug")!=null) { + this.debug = Boolean.parseBoolean(this.getEngineConfigString("debug")); + } + + this.nAlignables = ((this.svtTopBottomSep ? NSVTSENSORS : NSVTSENSORS/2) + (this.isSVTonly ? 0: 18) + (includeBeamspot? 1 : 0)); + this.setAlignVars(alignVars); + + return true; + } + + + + private void setAlignVars(String alignVars) { + orderTx = -1; + orderTy = -1; + orderTz = -1; + orderRx = -1; + orderRy = -1; + orderRz = -1; + if (alignVars.length() > 2 && !alignVars.contains(" ")) { + for (int i = 0; i < alignVars.length() / 2; i++) { + String s = alignVars.substring(2 * i, 2 * i + 2); + if (s.equals("Tx")) { + orderTx = i; + } else if (s.equals("Ty")) { + orderTy = i; + } else if (s.equals("Tz")) { + orderTz = i; + } else if (s.equals("Rx")) { + orderRx = i; + } else if (s.equals("Ry")) { + orderRy = i; + } else if (s.equals("Rz")) { + orderRz = i; + } + nAlignVars = i + 1; + } + if (debug) { + System.out.println(nAlignVars + " alignment variables requested"); + } + return; + } + //old version + String split[] = alignVars.split("[ \t]+"); + int i = 0; + for (String s : split) { + if (s.equals("Tx")) { + orderTx = i; + i++; + } else if (s.equals("Ty")) { + orderTy = i; + i++; + } else if (s.equals("Tz")) { + orderTz = i; + i++; + } else if (s.equals("Rx")) { + orderRx = i; + i++; + } else if (s.equals("Ry")) { + orderRy = i; + i++; + } else if (s.equals("Rz")) { + orderRz = i; + i++; + } + } + nAlignVars = i; + if(debug) System.out.println(nAlignVars + " alignment variables requested"); + } + + +} diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/banks/RecoBankReader.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/banks/RecoBankReader.java index 2b9eff4114..24cb47c9c3 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/banks/RecoBankReader.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/banks/RecoBankReader.java @@ -21,7 +21,10 @@ import org.jlab.rec.cvt.hit.Strip; import org.jlab.rec.cvt.svt.SVTGeometry; import org.jlab.rec.cvt.track.Seed; +import org.jlab.rec.cvt.track.StraightTrack; +import org.jlab.rec.cvt.track.Track; import org.jlab.rec.cvt.trajectory.Helix; +import org.jlab.rec.cvt.trajectory.Ray; /** * @@ -30,14 +33,14 @@ public class RecoBankReader { - public static List readBSTHitBank(DataEvent event) { + public static List readBSTHitBank(DataEvent event, String bankname) { - if(!event.hasBank("BST::Hits")) + if(!event.hasBank(bankname)) return null; else { List hits = new ArrayList<>(); - DataBank bank = event.getBank("BST::Hits"); + DataBank bank = event.getBank(bankname); for(int i = 0; i < bank.rows(); i++) { int id = bank.getShort("ID", i); int sector = bank.getByte("sector", i); @@ -67,14 +70,14 @@ public static List readBSTHitBank(DataEvent event) { } } - public static List readBMTHitBank(DataEvent event) { + public static List readBMTHitBank(DataEvent event, String bankname) { - if(!event.hasBank("BMT::Hits")) + if(!event.hasBank(bankname)) return null; else { List hits = new ArrayList<>(); - DataBank bank = event.getBank("BMT::Hits"); + DataBank bank = event.getBank(bankname); for(int i = 0; i < bank.rows(); i++) { int id = bank.getShort("ID", i); int sector = bank.getByte("sector", i); @@ -100,14 +103,14 @@ public static List readBMTHitBank(DataEvent event) { } } - public static Map readBSTClusterBank(DataEvent event, List svthits) { + public static Map readBSTClusterBank(DataEvent event, List svthits, String bankname) { - if(!event.hasBank("BST::Clusters")) + if(!event.hasBank(bankname)) return null; else { Map clusters = new HashMap<>(); - DataBank bank = event.getBank("BST::Clusters"); + DataBank bank = event.getBank(bankname); for(int i = 0; i < bank.rows(); i++) { int id = bank.getShort("ID", i); int tid = bank.getShort("trkID", i); @@ -168,14 +171,14 @@ public static Map readBSTClusterBank(DataEvent event, List readBMTClusterBank(DataEvent event, List bmthits) { + public static Map readBMTClusterBank(DataEvent event, List bmthits, String bankname) { - if(!event.hasBank("BMT::Clusters")) + if(!event.hasBank(bankname)) return null; else { Map clusters = new HashMap<>(); - DataBank bank = event.getBank("BMT::Clusters"); + DataBank bank = event.getBank(bankname); for(int i = 0; i < bank.rows(); i++) { int id = bank.getShort("ID", i); int tid = bank.getShort("trkID", i); @@ -218,7 +221,7 @@ public static Map readBMTClusterBank(DataEvent event, List readBMTClusterBank(DataEvent event, List readBSTCrossBank(DataEvent event, Map bstClusters) { + public static Map readBSTCrossBank(DataEvent event, Map bstClusters, String bankname) { - if(!event.hasBank("BST::Crosses")) + if(!event.hasBank(bankname)) return null; else { Map crosses = new HashMap<>(); - DataBank bank = event.getBank("BST::Crosses"); + DataBank bank = event.getBank(bankname); for(int i = 0; i < bank.rows(); i++) { int id = bank.getShort("ID", i); int tid = bank.getShort("trkID", i); @@ -306,14 +309,14 @@ public static Map readBSTCrossBank(DataEvent event, Map readBMTCrossBank(DataEvent event, Map bmtClusters) { + public static Map readBMTCrossBank(DataEvent event, Map bmtClusters, String bankname) { - if(!event.hasBank("BMT::Crosses")) + if(!event.hasBank(bankname)) return null; else { Map crosses = new HashMap<>(); - DataBank bank = event.getBank("BMT::Crosses"); + DataBank bank = event.getBank(bankname); for(int i = 0; i < bank.rows(); i++) { int id = bank.getShort("ID", i); int tid = bank.getShort("trkID", i); @@ -353,14 +356,14 @@ public static Map readBMTCrossBank(DataEvent event, Map readCVTSeedsBank(DataEvent event, double xb, double yb, Map svtCrosses, Map bmtCrosses) { + public static Map readCVTSeedsBank(DataEvent event, double xb, double yb, Map svtCrosses, Map bmtCrosses, String bankname) { - if(!event.hasBank("CVT::Seeds") || svtCrosses==null) + if(!event.hasBank(bankname) || svtCrosses==null) return null; else { Map seeds = new HashMap<>(); - DataBank bank = event.getBank("CVT::Seeds"); + DataBank bank = event.getBank(bankname); for(int i = 0; i < bank.rows(); i++) { int tid = bank.getShort("ID", i); double pt = bank.getFloat("pt", i); @@ -429,14 +432,14 @@ public static Map readCVTSeedsBank(DataEvent event, double xb, do } public static Map readCVTTracksBank(DataEvent event, double xb, double yb, Map cvtSeeds, - Map svtCrosses, Map bmtCrosses) { + Map svtCrosses, Map bmtCrosses, String bankname) { - if(!event.hasBank("CVT::Tracks")) + if(!event.hasBank(bankname)) return null; else { Map seeds = new HashMap<>(); - DataBank bank = event.getBank("CVT::Tracks"); + DataBank bank = event.getBank(bankname); for(int i = 0; i < bank.rows(); i++) { int tid = bank.getShort("ID", i); double pt = bank.getFloat("pt", i); @@ -577,4 +580,37 @@ public static Map readCVTTracksBank(DataEvent event, double xb, d } */ + public static List readCVTCosmicsBank(DataEvent event, String bankname) { + + if(!event.hasBank(bankname)) + return null; + else { + List tracks = new ArrayList<>(); + + DataBank bank = event.getBank(bankname); + for(int i = 0; i < bank.rows(); i++) { + int tid = bank.getShort("ID", i); + double chi2 = bank.getFloat("chi2", i); + int ndf = bank.getShort("ndf", i); + double yxSlope = bank.getFloat("trkline_yx_slope", i); + double yxInterc = bank.getFloat("trkline_yx_interc", i)*10; + double yzSlope = bank.getFloat("trkline_yz_slope", i); + double yzInterc = bank.getFloat("trkline_yz_interc", i)*10; + + Ray ray = new Ray(yxSlope, yxInterc, yzSlope, yzInterc); + StraightTrack track = new StraightTrack(ray); + track.setId(tid); + track.setChi2(chi2); + track.setNDF(ndf); + + for (int j = 0; j < 18; j++) { + int crossId = bank.getShort("Cross"+(j+1)+"_ID", i); + DetectorType det = crossId>1000 ? DetectorType.BMT : DetectorType.BST; + if(crossId>0) track.add(new Cross(det, BMTType.UNDEFINED, 0, 0, crossId)); + } + tracks.add(track); + } + return tracks; + } + } } diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/banks/RecoBankWriter.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/banks/RecoBankWriter.java index 581fb57def..156728653f 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/banks/RecoBankWriter.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/banks/RecoBankWriter.java @@ -995,7 +995,7 @@ public static DataBank fillStraightTrackKFTrajectoryBank(DataEvent event, List sVThits, List bMThits, List sVTclusters, List bMTclusters, diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/bmt/BMTGeometry.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/bmt/BMTGeometry.java index 671c3b26ae..6f53f00663 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/bmt/BMTGeometry.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/bmt/BMTGeometry.java @@ -300,7 +300,7 @@ public double getDPhi(int layer, int sector) { public double getThickness() { return BMTConstants.HDRIFT; - } + } /** * Return offset of the selected tile, identified by layer and sector @@ -805,7 +805,7 @@ public double getPhi(int layer, int sector, Point3D traj) { Point3D local = this.toLocal(traj, layer, sector); return local.toVector3D().phi(); } - + /** * Calculate Theta Lorentz based on solenoid scale and drift settings * @param layer diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/bmt/CCDBConstantsLoader.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/bmt/CCDBConstantsLoader.java index 3cfad3e2dd..3e8cbb96d1 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/bmt/CCDBConstantsLoader.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/bmt/CCDBConstantsLoader.java @@ -59,7 +59,7 @@ public static final synchronized void Load(DatabaseConstantProvider dbprovider) double[][] CRCEDGE1 = new double[NREGIONS][NSECTORS]; // the angle of the first edge of each PCB detector A, B, C double[][] CRCEDGE2 = new double[NREGIONS][NSECTORS]; // the angle of the second edge of each PCB detector A, B, C double[] CRCXPOS = new double[NREGIONS]; // Distance on the PCB between the PCB first edge and the edge of the first strip in mm - + int GRID_SIZE=405; double[] THETA_L_grid = new double [GRID_SIZE]; double[] ELEC_grid = new double [GRID_SIZE]; @@ -211,7 +211,7 @@ public static final synchronized void Load(DatabaseConstantProvider dbprovider) ELEC_grid[i]=dbprovider.getDouble("/calibration/mvt/lorentz/Edrift",i); MAG_grid[i]=dbprovider.getDouble("/calibration/mvt/lorentz/Bfield",i); } - + // alignment and offsets double xpos = dbprovider.getDouble("/geometry/cvt/mvt/position/x", 0 ); double ypos = dbprovider.getDouble("/geometry/cvt/mvt/position/y", 0 ); @@ -254,7 +254,7 @@ public static final synchronized void Load(DatabaseConstantProvider dbprovider) BMTConstants.TOLOCAL[layer-1][sector-1] = transform.inverse(); } - + BMTConstants.setCRCRADIUS(CRCRADIUS); BMTConstants.setCRZRADIUS(CRZRADIUS); BMTConstants.setCRZNSTRIPS(CRZNSTRIPS); diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/cross/Cross.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/cross/Cross.java index fcc5cc01e9..0b1ce4fea8 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/cross/Cross.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/cross/Cross.java @@ -380,7 +380,7 @@ public void update(Point3D trackPos, Vector3D trackDir) { * Sets the cross parameters: the position and direction unit vector */ public void updateBMTCross(Point3D trackPos, Vector3D trackDir) { - + if(this.getDetector()!=DetectorType.BMT) return; Point3D crossPoint = this.getBMTCrossPoint(trackPos); @@ -590,7 +590,7 @@ public String toString() { public String printInfo() { String s = " cross: " + this.getDetector() + " " + this.getType() + " ID " + this.getId() + " Sector " + this.getSector() + " Region " + this.getRegion() + " sort "+this.getOrderedRegion()+" cosmic region "+this.getSVTCosmicsRegion(); - if(this.getPoint0()!=null) s += " Point " + this.getPoint().toString(); + if(this.getPoint()!=null) s += " Point " + this.getPoint().toString(); if(this.getPoint0()!=null) s += " Point0 " + this.getPoint0().toString(); if(this.getDir()!=null) s += " Direction "+ this.getDir().toString(); return s; diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/hit/Hit.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/hit/Hit.java index 3c7116bb38..b24add63dd 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/hit/Hit.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/hit/Hit.java @@ -279,7 +279,7 @@ public void setAssociatedClusterID(int _AssociatedClusterID) { public int getAssociatedTrackID() { return AssociatedTrackID; } - + public void setAssociatedTrackID(int associatedTrackID) { AssociatedTrackID = associatedTrackID; diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/CVTEngine.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/CVTEngine.java index e278b0c9c7..a362de3291 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/CVTEngine.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/CVTEngine.java @@ -276,9 +276,9 @@ public void setBmtzmaxclussize(int bmtzmaxclussize) { @Override public boolean processDataEvent(DataEvent event) { - + Swim swimmer = new Swim(); - + int run = this.getRun(event); IndexedTable svtStatus = this.getConstantsManager().getConstants(run, "/calibration/svt/status"); diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/CosmicTracksRec.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/CosmicTracksRec.java index ade5c0d319..5a53ef846e 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/CosmicTracksRec.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/CosmicTracksRec.java @@ -208,4 +208,4 @@ public List getBMTcrosses() { else return CVTcrosses.get(1); } -} \ No newline at end of file +} diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/RecUtilities.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/RecUtilities.java index 0204c84d20..140467feae 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/RecUtilities.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/RecUtilities.java @@ -941,4 +941,4 @@ Track recovTrkMisClusSearch(Seed seed, org.jlab.clas.tracking.trackrep.Helix hlx } return fittedTrack; } -} \ No newline at end of file +} diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/TracksFromTargetRec.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/TracksFromTargetRec.java index d9ffecf6c5..4097b27fa3 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/TracksFromTargetRec.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/services/TracksFromTargetRec.java @@ -37,7 +37,7 @@ * @author ziegler */ public class TracksFromTargetRec { - + private final RecUtilities recUtil = new RecUtilities(); private List SVThits; @@ -368,6 +368,7 @@ public List getTracks(DataEvent event, boolean initFromMc, boolean kfFilt for(int it = 0; it < tracks.size(); it++) { int id = it + 1; tracks.get(it).setId(id); +// System.out.println("Fit " + tracks.get(it).toString()); tracks.get(it).findTrajectory(swimmer, Geometry.getInstance().geOuterSurfaces()); tracks.get(it).update_Crosses(id, xb, yb); tracks.get(it).update_Clusters(id); @@ -417,8 +418,8 @@ public List getSeedsFromBanks(DataEvent event) { this.init(); - SVThits = RecoBankReader.readBSTHitBank(event); - BMThits = RecoBankReader.readBMTHitBank(event); + SVThits = RecoBankReader.readBSTHitBank(event, "BST::Hits"); + BMThits = RecoBankReader.readBMTHitBank(event, "BMT::Hits"); if(SVThits!= null) { Collections.sort(SVThits); } @@ -430,19 +431,19 @@ public List getSeedsFromBanks(DataEvent event) { } - SVTclustersHM = RecoBankReader.readBSTClusterBank(event, SVThits); - BMTclustersHM = RecoBankReader.readBMTClusterBank(event, BMThits); + SVTclustersHM = RecoBankReader.readBSTClusterBank(event, SVThits, "BST::Clusters"); + BMTclustersHM = RecoBankReader.readBMTClusterBank(event, BMThits, "BMT::Clusters"); - SVTcrossesHM = RecoBankReader.readBSTCrossBank(event, SVTclustersHM); - BMTcrossesHM = RecoBankReader.readBMTCrossBank(event, BMTclustersHM); + SVTcrossesHM = RecoBankReader.readBSTCrossBank(event, SVTclustersHM, "BST::Crosses"); + BMTcrossesHM = RecoBankReader.readBMTCrossBank(event, BMTclustersHM, "BMT::Crosses"); - CVTseedsHM = RecoBankReader.readCVTSeedsBank(event, xb, yb, SVTcrossesHM, BMTcrossesHM); + CVTseedsHM = RecoBankReader.readCVTSeedsBank(event, xb, yb, SVTcrossesHM, BMTcrossesHM, "CVT::Seeds"); if(CVTseedsHM == null) { return null; } else { - CVTseedsHMOK = RecoBankReader.readCVTTracksBank(event, xb, yb, CVTseedsHM, SVTcrossesHM, BMTcrossesHM); + CVTseedsHMOK = RecoBankReader.readCVTTracksBank(event, xb, yb, CVTseedsHM, SVTcrossesHM, BMTcrossesHM, "CVT::Tracks"); if(CVTseedsHMOK == null) { return null; } else { diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/svt/SVTGeometry.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/svt/SVTGeometry.java index 9f2b529f5c..43eab9ca95 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/svt/SVTGeometry.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/svt/SVTGeometry.java @@ -87,7 +87,7 @@ public static double getActiveSensorLength() { public static double getModuleLength() { return SVTConstants.MODULELEN; } - + public Line3D getStrip0(int layer, int sector, int strip) { Line3d line = this._svtStripFactory.getShiftedStrip(layer-1, sector-1, strip-1); return new Line3D(line.origin().x,line.origin().y,line.origin().z, @@ -246,7 +246,7 @@ public Surface getSurface(int layer, int sector, Strip strip) { surface.hemisphere = 0; surface.setLayer(layer); surface.setSector(sector); - surface.setError(0); + surface.setError(0); for(String key : SVTConstants.MATERIALPROPERTIES.keySet()) { double[] p = SVTConstants.MATERIALPROPERTIES.get(key); surface.addMaterial(key, p[0]/2, p[1], p[2], p[3], p[4], Units.MM); diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/track/Track.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/track/Track.java index eea7abe5e6..25f36eaa9d 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/track/Track.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/track/Track.java @@ -515,6 +515,8 @@ public void findTrajectory(Swim swimmer, List outer) { stVec.setEdge(Geometry.getInstance().getSVT().distanceToEdge(traj.layer, traj.sector, pos)); } else if(MLayer.getDetectorType(index) == DetectorType.BMT) { +// System.out.println(pos.toString() + " " + dir.toString()); +// if(Double.isNaN(pos.x())) System.exit(1); Vector3D localDir = Geometry.getInstance().getBMT().getLocalTrack(traj.layer, traj.sector, pos, dir); stVec.setTrkPhiAtSurface(localDir.phi()); stVec.setTrkThetaAtSurface(localDir.theta()); @@ -658,6 +660,9 @@ public String toString() { this.getNDF(), this.getChi2(), this.getSeed().getStatus(), this.getKFIterations(), this.getStatus()); for(Cross c: this) str = str + c.toString() + "\n"; for(Cluster c: this.getSeed().getClusters()) str = str + c.toString() + "\n"; + if(this.trajs!=null) { + for(HitOnTrack h : trajs.values()) str = str + h.toString() + "\n"; + } return str; } diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/track/TrackSeeder.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/track/TrackSeeder.java index e342c15d6a..7020f5e517 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/track/TrackSeeder.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/track/TrackSeeder.java @@ -446,7 +446,7 @@ public List findCandUsingMicroMegas(Seed trkCand, List bmt_crosses) BMTTrkSeed.setHelix(trkCand.getHelix()); BMTTrkSeed.setCrosses(matches); AllSeeds.put(st,BMTTrkSeed); - + //if (AllSeeds.size() > 200) { // AllSeeds.clear(); // return AllSeeds; diff --git a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/TrajectoryFinder.java b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/TrajectoryFinder.java index 0e4fa6d690..156ce51cda 100644 --- a/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/TrajectoryFinder.java +++ b/reconstruction/cvt/src/main/java/org/jlab/rec/cvt/trajectory/TrajectoryFinder.java @@ -31,7 +31,7 @@ public class TrajectoryFinder { public TrajectoryFinder() { } - + public Trajectory findTrajectory(int id, Ray ray, ArrayList candCrossList) { Trajectory traj = new Trajectory(ray); traj.setId(id);