24 removed files
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/AddFitKalman.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/AddFitKalman.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,340 +0,0 @@
-package org.hps.recon.tracking.kalman;
-
-import java.util.ArrayList;
-import java.util.List;
-
-import org.lcsim.recon.tracking.trfbase.ETrack;
-import org.lcsim.recon.tracking.trfbase.Hit;
-import org.lcsim.recon.tracking.trfbase.TrackError;
-import org.lcsim.recon.tracking.trfbase.TrackVector;
-import org.lcsim.recon.tracking.trfutil.Assert;
-
-import Jama.Matrix;
-// Fit tracks using Kalman filter.
-
-/**
- *AddFitKalman uses a Kalman filter to update the track fit.
- *
- *
- *@author Norman A. Graf
- *@version 1.0
- *
- */
-public class AddFitKalman extends AddFitter {
-
- boolean _DEBUG = false;
- // Maximum allowed hit dimension.
- private static final int MAXDIM = 3;
- // Maximum number of track vectors.
- private static final int MAXVECTOR = 1;
- // Maximum number of track errors.
- private static final int MAXERROR = 1;
-
- //private: // nested classes
- // The nested class Box holds the vectors, matrices and symmetric
- // matrices needed for adding a hit in the main class.
- class Box {
- /*
- private: // typedefs
- typedef Ptr<TrfVector,AutoPolicy> VectorPtr;
- typedef Ptr<TrfSMatrix,AutoPolicy> SMatrixPtr;
- typedef Ptr<TrfMatrix,AutoPolicy> MatrixPtr;
- typedef vector<VectorPtr> VectorList;
- typedef vector<SMatrixPtr> SMatrixList;
- typedef vector<MatrixPtr> MatrixList;
- */
-
- // enums
- // number of vectors
- private static final int NVECTOR = 2;
- // number of errors
- private static final int NERROR = 3;
- // number of vectors
- private static final int NDERIV = 2;
- // number of gains
- private static final int NGAIN = 2;
- // attributes
- // dimension of the vector, matrix, etc
- private int _size;
- // array of vectors
- List _vectors;
- // array of error matrices
- List _errors;
- // array of derivatives (Nx5 matrices)
- List _derivs;
- // array of gains (5xN matrices)
- List _gains;
-
- // methods
- // constructor
- public Box(int size) {
- _size = size;
- _vectors = new ArrayList();
- _errors = new ArrayList();
- _derivs = new ArrayList();
- _gains = new ArrayList();
-
- int icnt;
- // Track vectors
- for (icnt = 0; icnt < NVECTOR; ++icnt)
- _vectors.add(new Matrix(size, 1));
- // Track errors
- for (icnt = 0; icnt < NERROR; ++icnt)
- _errors.add(new Matrix(size, size));
- // Track derivatives
- for (icnt = 0; icnt < NDERIV; ++icnt)
- _derivs.add(new Matrix(size, 5));
- // Gains
- for (icnt = 0; icnt < NDERIV; ++icnt)
- _gains.add(new Matrix(5, size));
-
- }
-
- // return the dimension
- public int get_size() {
- return _size;
- }
- // fetch a vector
-
- public Matrix get_vector(int ivec) {
- return (Matrix) _vectors.get(ivec);
- }
- // fetch an error
-
- public Matrix get_error(int ierr) {
- return (Matrix) _errors.get(ierr);
- }
- // fetch a derivative
-
- public Matrix get_deriv(int ider) {
- return (Matrix) _derivs.get(ider);
- }
- // fetch a gain
-
- public Matrix get_gain(int igai) {
- return (Matrix) _gains.get(igai);
- }
- } // end of Box inner class
-
- // static methods
- //
- /**
- *Return a String representation of the class' the type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' the type name.
- */
- public static String typeName() {
- return "AddFitKalman";
- }
-
- //
- /**
- *Return a String representation of the class' the type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' the type name.
- */
- public static String staticType() {
- return typeName();
- }
- // attributes
- // Array of boxes for each supported size.
- private List _boxes;
- // Track vectors.
- private List _tvectors;
- // Track errors.
- private List _terrors;
-
- //methods
- //
- /**
- *Construct a default instance.
- * Allocate space needed for hits of dimension from 1 to MAXDIM.
- *
- */
- public AddFitKalman() {
- _boxes = new ArrayList();
- _tvectors = new ArrayList();
- _terrors = new ArrayList();
-
- // Create boxes for hit containers.
- for (int dim = 1; dim < MAXDIM; ++dim)
- _boxes.add(new Box(dim));
-
- int icnt;
-
- for (icnt = 0; icnt < MAXVECTOR; ++icnt)
- _tvectors.add(new Matrix(5, 1));
-
- for (icnt = 0; icnt < MAXERROR; ++icnt)
- _terrors.add(new Matrix(5, 5));
-
- }
-
- //
- /**
- *Return a String representation of the class' the type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' the type name.
- */
- public String type() {
- return staticType();
- }
-
- //
- /**
- *Add a hit and fit with the new hit.
- * Use a Kalman filter to add a hit to a track.
- * The hit is updated with the input track.
- * Note: We make direct use of the underlying vector and matrix
- * classes here. It will probably be neccessary to modify
- * this routine if these are changed.
- *
- * @param tre The ETrack to update.
- * @param chsq The chi-square for the fit.
- * @param hit The Hit to add to the track.
- * @return 0 if successful.
- */
- public int addHitFit(ETrack tre, double chsq, Hit hit) {
- // Update the hit with the input track.
- hit.update(tre);
-
- // Fetch hit size.
- int dim = hit.size();
- Assert.assertTrue(dim <= MAXDIM);
-
- // Fetch the box holding the needed hit containers.
- // The chice of boxes depends on the size of the hit.
- Box box = (Box) _boxes.get(dim - 1);
- Assert.assertTrue(box.get_size() == dim);
-
- // Fetch the hit containers.
- Matrix diff = box.get_vector(0);
- Matrix hit_res = box.get_vector(1);
- Matrix hit_err = box.get_error(0);
- Matrix hit_err_tot = box.get_error(1);
- Matrix hit_res_err = box.get_error(2);
- Matrix dhit_dtrk = box.get_deriv(0);
- Matrix new_dhit_dtrk = box.get_deriv(1);
- Matrix trk_err_dhit_dtrk = box.get_gain(0);
- Matrix gain = box.get_gain(1);
-
- // Fetch the track containers.
- Matrix new_vec = (Matrix) _tvectors.get(0);
- Matrix new_err = (Matrix) _terrors.get(0);
-
- hit_err = hit.measuredError().matrix();
- //System.out.println("hit_err= \n"+hit_err);
-
- // Fetch track prediction of hit.
- diff = hit.differenceVector().matrix();
- //System.out.println("diff= \n"+diff);
- dhit_dtrk = hit.dHitdTrack().matrix();
- //System.out.println("dhit_dtrk= \n"+dhit_dtrk);
- // Fetch track info.
- Matrix trk_vec = tre.vector().matrix();
- Matrix trk_err = tre.error().getMatrix(); //need to fix this!
-
- //System.out.println("trk_vec= \n"+trk_vec);
- //System.out.println("trk_err= \n"+trk_err);
-
- // Build gain matrix.
- hit_err_tot = hit.predictedError().matrix().plus(
- hit_err);
- //System.out.println("hit_err_tot= \n"+hit_err_tot);
- //if ( invert(hit_err_tot)!=0 ) return 3;
- hit_err_tot = hit_err_tot.inverse();
- //System.out.println("hit_err_tot inverse= \n"+hit_err_tot);
- trk_err_dhit_dtrk = trk_err.times(dhit_dtrk.transpose());
- gain = trk_err_dhit_dtrk.times(hit_err_tot);
- //System.out.println("trk_err_dhit_dtrk= \n"+trk_err_dhit_dtrk);
- //System.out.println("gain= \n"+gain);
-
- // if ( get_debug() ) {
- //System.out.println("\n");
- //System.out.println(" trk_vec: " + "\n"+ trk_vec + "\n");
- //System.out.println(" trk_err: " + "\n"+ trk_err + "\n");
- //System.out.println(" dhit_dtrk: " + "\n"+ dhit_dtrk + "\n");
- // }
-
- // We need to return dhit_dtrk to its original state for the
- // next call.
- // dhit_dtrk = dhit_dtrk.transpose(); //need to check this!
- dhit_dtrk.transpose();
- // Build new track vector.
- new_vec = trk_vec.minus(gain.times(diff));
- //System.out.println("new_vec= \n"+new_vec);
-
- // Build new error;
- new_err = trk_err.minus(trk_err_dhit_dtrk.times(hit_err_tot.times(trk_err_dhit_dtrk.transpose())));
- //System.out.println("new_err= \n"+new_err);
- // Check the error.
- {
- int nbad = 0;
- for (int i = 0; i < 5; ++i) {
- if (new_err.get(i, i) < 0.0)
- ++nbad;
- double eii = new_err.get(i, i);
- for (int j = 0; j < i; ++j) {
- double ejj = new_err.get(j, j);
- double eij = new_err.get(j, i);
- if (Math.abs(eij * eij) >= eii * ejj)
- ++nbad;
- }
- }
- if (nbad > 0)
- return 5;
- }
-
- // Create track vector with new values.
- tre.setVectorAndKeepDirection(new TrackVector(new_vec));
- tre.setError(new TrackError(new_err));
-
- // Calculate residual vector.
-
- // Update the hit with the new track.
- //System.out.println("update the hit");
- hit.update(tre);
-
- hit_res = hit.differenceVector().matrix();
- new_dhit_dtrk = hit.dHitdTrack().matrix();
- //System.out.println("new_dhit_dtrk= \n"+new_dhit_dtrk);
-
- // Calculate residual covariance and invert.
- hit_res_err = hit_err.minus(dhit_dtrk.times(new_err.times(dhit_dtrk.transpose())));
- // System.out.println("hit_res_err= \n"+hit_res_err);
- hit_res_err = hit_res_err.inverse();
- //System.out.println("hit_res_err inverse= \n"+hit_res_err);
- // Update chi-square.
- // result should be 1x1 matrix, so should be able to do the following
- //System.out.println( " hr*hre*hrT=\n"+(hit_res.transpose()).times(hit_res_err.times(hit_res)) );
- double dchsq = (hit_res.transpose()).times(hit_res_err.times(hit_res)).get(0, 0);
- //System.out.println("chsq= "+chsq+", dchsq= "+dchsq);
- chsq = chsq + dchsq;
- setChisquared(chsq);
-
- if (_DEBUG) {
- System.out.println(" gain: " + "\n" + gain + "\n");
- System.out.println(" new_vec: " + "\n" + new_vec + "\n");
- System.out.println(" new_err: " + "\n" + new_err + "\n");
- System.out.println("new_dhit_dtrk: " + "\n" + new_dhit_dtrk + "\n");
- System.out.println(" hit_res: " + "\n" + hit_res + "\n");
- System.out.println(" hit_res_err: " + "\n" + hit_res_err + "\n");
- System.out.println(" dchsq & chsq: " + "\n" + dchsq + " " + chsq + "\n");
- }
-
- return 0;
-
- }
-
- /**
- *output stream
- *
- * @return The String representation of this instance.
- */
- public String toString() {
- return getClass().getName();
- }
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/AddFitter.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/AddFitter.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,178 +0,0 @@
-package org.hps.recon.tracking.kalman;
-import org.lcsim.recon.tracking.trfbase.ETrack;
-import org.lcsim.recon.tracking.trfbase.Hit;
-import org.lcsim.recon.tracking.trfbase.Surface;
-import org.lcsim.recon.tracking.trffit.HTrack;
-import org.lcsim.recon.tracking.trfutil.Assert;
-/**
- * Adds a hit prediction to a track, refits the track and uses the
- * new track parameters to update the hit prediction.
- *<p>
- * Two methods are provided: one updates an ETrack and its chi-square
- * and the other updates an HTrack.
- *<p>
- * This base class will always return an error. Subclasses will typically
- * implement the ETrack method which is then invoked by the HTrack method
- * defined here.
- *
- *@author Norman A. Graf
- *@version 1.0
- *
- */
-
-public class AddFitter
-{
-
- // Static methods.
-
- //
-
- /**
- *Return String representation of this class' type name.
- *Included for completeness with the C++ versin.
- *
- * @return String representation of this class' type.
- */
- public static String typeName()
- { return "AddFitter"; }
-
- //
-
- /**
- *Return String representation of this class' type name.
- *Included for completeness with the C++ versin.
- *
- * @return String representation of this class' type.
- */
- public static String staticType()
- { return typeName(); }
-
- // workaround
- private double _chsq;
-
- // methods
-
- //
-
- /**
- *Construct a default instance.
- *
- */
- public AddFitter()
- {
- _chsq = 0.;
- }
-
- //
-
- /**
- *Return the generic type.
- * This is only needed at this level.
- *
- * @return String representation of this class' type.
- */
- public String genericType()
- { return staticType(); }
-
- //
-
- /**
- *Add a hit and fit with the new hit.
- * Return status 0 if fit is successful, negative value for a local
- * error and positive for an error in add_hit_fit.
- * The default method calls add_hit_fit and return its status.
- * If the fit is successful, then the track fit is updated and the hit
- * is added to the end of its list.
- *
- * @param trh The HTrack to which the hit will be added.
- * @param hit The Hit to add.
- * @return 0 if hit update and fit are successful.
- */
- public int addHit(HTrack trh, Hit hit)
- {
- // Fetch the starting fit and chi-square.
- ETrack tre = trh.newTrack();
- double chsq = trh.chisquared();
-
- // check the track and hit are at the same surface
- Surface tsrf = tre.surface();
- Surface hsrf = hit.surface();
- Assert.assertTrue( tsrf.pureEqual(hsrf) );
- if ( ! tsrf.pureEqual(hsrf) ) return -1;
-
- // Check the track is fully fit before adding hit.
- // Unless this is the first hit.
- if ( trh.hits().size() !=0 )
- {
- Assert.assertTrue( trh.isFit() );
- if ( ! trh.isFit() ) return -2;
- }
- // Fit with the new point; exit if error occurs.
- int stat = addHitFit(tre,chsq,hit); //chsq is return argument in c++
- // need to fix this
- if ( stat != 0 ) return stat;
-
-
- // Update the track with the new fit and hit.
- trh.addHit(hit);
- trh.setFit(tre,_chsq);
-
- return 0;
-
- }
-
-
- /**
- *Set the chi-squared for the fit.
- *
- * @param chsq The value of chi-square to set for this fit.
- */
- public void setChisquared(double chsq)
- {
- _chsq = chsq;
- }
-
-
- /**
- *Return the chi-squared for the fit.
- *
- * @return The value of chi-square for this fit.
- */
- public double chisquared()
- {
- return _chsq;
- }
- //
-
-
- /**
- *Refit a track and update its chi-square by adding the specified hit.
- * Return status 0 if fit is successful, positive value for error.
- * This is the method implemented by subclasses.
- * If the fit fails, the track and chi-square may return any value
- * and the hit may be updated with any track.
- * If the fit is successful, the fit track and chi-square are
- * returned and the hit is updated with the fit track.
- * Normally this is not invoked directly but is called by add_hit.
- * The default method here returns an error and throws and AssertException.
- *
- * @param tre The ETrack to which the hit will be added.
- * @param chisq The value of chi-square to set for this fit.
- * @param phit The Hit to add.
- * @return 1.
- */
- public int addHitFit(ETrack tre, double chisq, Hit phit)
- { Assert.assertTrue( false );
- return 1;}
-
-
- /**
- *output stream
- *
- * @return AString representation of this object.
- */
- public String toString()
- {
- return getClass().getName();
- }
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/BoxHelper.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/BoxHelper.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,45 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-
-package org.hps.recon.tracking.kalman;
-
-import org.lcsim.detector.solids.Box;
-import org.lcsim.detector.solids.ISolid;
-import org.lcsim.event.Track;
-
-/**
- *
- * @author ecfine
- */
-public class BoxHelper implements ShapeHelper{
- public void printShape(ISolid solid) {
- checkBox(solid);
- System.out.println("Shape: Box");
- }
-
- public void printLocalCoords(ISolid solid) {
-
- }
-
- public void printGlobalCoords(ISolid solid) {
-
- }
-
- private void checkBox(ISolid solid){
- if (solid instanceof Box) {
- } else {
- System.out.println("Error! This is not a box! ");
- return;
- }
- }
-
- public void findIntersection(ISolid solid, Track track) {
- throw new UnsupportedOperationException("Not supported yet.");
- }
-
- public KalmanSurface getKalmanSurf(ISolid solid) {
- throw new UnsupportedOperationException("Not supported yet.");
- }
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/FullFitKalman.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/FullFitKalman.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,573 +0,0 @@
-package org.hps.recon.tracking.kalman;
-
-import java.util.Iterator;
-import java.util.List;
-import java.util.ListIterator;
-
-import org.hps.recon.tracking.kalman.util.RKDebug;
-import org.hps.recon.tracking.kalman.util.SurfaceCode;
-import org.lcsim.recon.tracking.trfbase.ETrack;
-import org.lcsim.recon.tracking.trfbase.Hit;
-import org.lcsim.recon.tracking.trfbase.PropDir;
-import org.lcsim.recon.tracking.trfbase.PropStat;
-import org.lcsim.recon.tracking.trfbase.Propagator;
-import org.lcsim.recon.tracking.trfbase.Surface;
-import org.lcsim.recon.tracking.trfbase.TrackError;
-import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
-import org.lcsim.recon.tracking.trfcyl.ThinCylMs;
-import org.lcsim.recon.tracking.trffit.FullFitter;
-import org.lcsim.recon.tracking.trffit.HTrack;
-import org.lcsim.recon.tracking.trfzp.SurfZPlane;
-import org.lcsim.recon.tracking.trfzp.ThinZPlaneMs;
-import org.lcsim.util.aida.AIDA;
-
-/**
- *
- * Copied from org.lcsim.contrib.RobKutschke.TRF.trffit. Added some very naive
- * multiple scattering, which will need to be updated.
- *
- *
- * Full track fit using Kalman filter. The propagator is specified
- * when the fitter is constructed. The starting surface, vector and
- * error matrix are taken from the input track. Errors should be
- * increased appropriately if the fitter is applied repeatedly to
- * a single track.
- *
- *@author $Author: mgraham $
- *@version $Id: FullFitKalman.java,v 1.6 2011/11/16 18:00:03 mgraham Exp $
- *
- * Date $Date: 2011/11/16 18:00:03 $
- *
- */
-public class FullFitKalman extends FullFitter {
-
- boolean _DEBUG = true;
- private AIDA aida = AIDA.defaultInstance();
- // Flags to control: multiple scattering, energy loss and adding the hit.
- private boolean doMs = true;
- private boolean doEloss = true;
- private boolean doMeas = true;
- private double dedxscale = 1.;
- private double dedxsigma = 0.0;
-
- // static methods
- //
- /**
- *Return a String representation of the class' the type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' the type name.
- */
- public static String typeName() {
- return "FullFitKalman";
- }
-
- //
- /**
- *Return a String representation of the class' the type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' the type name.
- */
- public static String staticType() {
- return typeName();
- }
- // The propagator.
- private Propagator _pprop;
- // The add fitter.
- private AddFitKalman _addfit;
- int AddFitKalmanDebugLevel = 0;
- //
-
- /**
- *Construct an instance specifying a propagator.
- *
- * @param prop The Propagator to be used during the fit.
- */
- public FullFitKalman(Propagator prop) {
- _pprop = prop;
- _addfit = new AddFitKalman();
-
-// try{
-// ToyConfig config = ToyConfig.getInstance();
-// AddFitKalmanDebugLevel = config.getInt( "AddFitKalmanDebugLevel",
-// AddFitKalmanDebugLevel );
-// dedxscale = config.getDouble("dEdXScale");
-// dedxsigma = config.getDouble("dEdXSigma");
-//
-//
-// } catch (ToyConfigException e){
-// System.out.println (e.getMessage() );
-// System.out.println ("Stopping now." );
-// System.exit(-1);
-// }
-// System.out.println ("FullfitKalman dedxscale: " + dedxscale );
-
- }
-
- //
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' the type name.
- */
- public String type() {
- return staticType();
- }
-
- //
- /**
- *Return the propagator.
- *
- * @return The Propagator used in the fit.
- */
- public Propagator propagator() {
- return _pprop;
- }
-
- //
- public void setDoMs(boolean b) {
- doMs = b;
- }
-
- public void setDoEloss(boolean b) {
- doEloss = b;
- }
-
- /**
- *Fit the specified track.
- *
- * @param trh The HTrack to fit.
- * @return 0 if successful.
- */
- public int fit(HTrack trh) {
- // Copy the hits from the track.
- List hits = trh.hits();
- if (_DEBUG)
- System.out.println("FullFitKalman::Hits has " + hits.size() + " elements");
- // Delete the list of hits from the track.
- while (trh.hits().size() > 0)
- trh.dropHit();
-// System.out.println("Hits has "+hits.size()+" elements");
-
- // Set direction to be nearest.
- PropDir dir = PropDir.NEAREST;
-// PropDir dir = PropDir.FORWARD;
-// RKDebug.Instance().setPropDir(dir);
-
- // Loop over hits and fit.
- int icount = 0;
- for (Iterator ihit = hits.iterator(); ihit.hasNext();) {
-
- // Extract the next hit pointer.
- Hit hit = (Hit) ihit.next();
- if (_DEBUG) {
- System.out.println("******* ADDING NEW HIT *********** ");
- System.out.println("Hit " + icount + " is: \n");
- System.out.println("Before prop: " + trh.newTrack());
- System.out.println("Propogating to surface : " + hit.surface());
- }
- // System.out.println("Before prop spacepoint: " + hit.surface().spacePoint(trh.newTrack().vector()));
-
- // propagate to the surface
- PropStat pstat = trh.propagate(_pprop, hit.surface(), dir);
-
- if (_DEBUG) {
- System.out.println("pstat = " + pstat);
- System.out.println("After prop: " + trh.newTrack());
- }
- if (!pstat.success())
- return -666;
-// System.out.println("trh= \n"+trh+", hit= \n"+hit);
-// System.out.println("_addfit= "+_addfit);
-// System.out.println("After prop spacepoint: " + hit.cluster().surface().spacePoint(trh.newTrack().vector()));
-// System.out.println("trh surface: " + trh.newTrack().surface().toString());
-// System.out.println("hit surface" + hit.surface().toString());
-
-// fit track
- if (_DEBUG)
- System.out.println("Predicted Hit Vector before addHit: " + hit.predictedVector());
- int fstat = _addfit.addHit(trh, hit);
- if (_DEBUG)
- System.out.println("Predicted Hit Vector after addHit: " + hit.predictedVector());
-//
- if (_DEBUG) {
- System.out.println("Measured Hit Vector after addHit: " + hit.cluster());
- System.out.println("After addhit: " + fstat + " chi^2 = "
- + trh.chisquared() + "\n" + trh.newTrack());
-// System.out.println("After addhit spacepoint: " + hit.cluster().surface().spacePoint(trh.newTrack().vector()));
-//
- }
-
- // Multiple scattering--this is a silly way to do it..
- ThinXYPlaneMs interactor = new ThinXYPlaneMs(.02);
-
- if (_DEBUG)
- System.out.println("trh error pre MS: " + trh.newTrack().error().toString());
-// ThinXYPlaneMs interactor = new ThinXYPlaneMs(.005);
- ETrack tre2 = trh.newTrack();
- interactor.interact(tre2);
- trh.setFit(tre2, trh.chisquared());
- if (_DEBUG)
- System.out.println("trh error pre MS: " + trh.newTrack().error().toString());
- if (fstat > 0)
- return 10000 + 1000 * fstat + icount;
-
-
- icount++;
-
- }
- if (_DEBUG)
- System.out.println("fit completed");
- return 0;
- }
-
- public int fitForward(HTrack trh) {
- if (_DEBUG)
- System.out.println("fitting forward...");
-
- PropDir dir = PropDir.FORWARD;
-// RKDebug.Instance().setPropDir(dir);
-
-
- // Copy the hits from the track.
- List hits = trh.hits();
-
- // Delete the list of hits from the track.
- while (trh.hits().size() > 0)
- trh.dropHit();
-
- double sumde = 0.;
-
- // Loop over hits and fit.
- int icount = 0;
- for (Iterator ihit = hits.iterator(); ihit.hasNext();) {
- Surface s_save = trh.newTrack().surface().newPureSurface();
- ETrack e_save = trh.newTrack();
-
- // Extract the next hit pointer.
- Hit hit = (Hit) ihit.next();
-
- int from = (new SurfaceCode(s_save)).getCode();
- int to = (new SurfaceCode(hit.surface())).getCode();
-
- // Propagate to the next surface.
- PropStat pstat = trh.propagate(_pprop, hit.surface(), dir);
- if (!pstat.success()) {
- if (AddFitKalmanDebugLevel > 0) {
- System.out.println("Error: "
- // + RKDebug.Instance().getTrack() + " "
- // + RKDebug.Instance().getPropDir() + " "
- + icount);
- System.out.println("From surface 5: " + s_save);
- System.out.println("To surface 5: " + hit.surface());
- System.out.println("Params: " + e_save.vector());
- }
-// aida.histogram1D("/Bugs/Fit/Failed Fwd prop from Surface",5,0,5).fill( from );
-// aida.histogram1D("/Bugs/Fit/Failed Fwd prop to Surface",5,0,5).fill( to );
-// aida.cloud2D("/Bugs/Fit/Failed Fwd prop to vs from Surface").fill( from, to );
- if (_DEBUG)
- System.out.println("pstat not success :(");
- return icount + 1;
- }
-//
- if (icount != 0) {
- int istat = interact(trh, hit, dir);
- }
-
-
- // Add the hit.
- int fstat = _addfit.addHit(trh, hit);
- if (_DEBUG)
- System.out.println("Hit added!");
- if (fstat > 0) {
- if (AddFitKalmanDebugLevel > 0) {
-
-// System.out.println("Error: "
-// + RKDebug.Instance().getTrack() + " "
-// + RKDebug.Instance().getPropDir() + " ");
- System.out.println("From surface 4: " + s_save);
- System.out.println("To surface 4: " + hit.surface());
- }
-// aida.histogram1D("/Bugs/Fit/Failed Fwd addhit from Surface",5,5,5).fill( from );
-// aida.histogram1D("/Bugs/Fit/Failed Fwd addhit to Surface",5,0,5).fill( to );
-// aida.cloud2D("/Bugs/Fit/Failed Fwd addhit to vs from Surface").fill( from, to );
-//
- }
- /*
- if (fstat > 0)
- return 10000 + 1000 * fstat + icount;
-
- VTUtil before = new VTUtil(trh.newTrack());
- int istat = interact(trh, hit, dir);
- VTUtil after = new VTUtil(trh.newTrack());
-
- double de = before.e() - after.e();
- sumde += de;
-
- SurfCylinder ss = (SurfCylinder) trh.newTrack().surface();
- if (_DEBUG)System.out.printf("Forw dedx: %10.4f %12.8f %12.8f\n", ss.radius(), de, sumde);
- */
-
- // Multiple scattering--this is a silly way to do it..
- ThinXYPlaneMs interactor = new ThinXYPlaneMs(.02);
-
- if (_DEBUG)
- System.out.println("trh error pre MS: " + trh.newTrack().error().toString());
-// ThinXYPlaneMs interactor = new ThinXYPlaneMs(.005);
- ETrack tre2 = trh.newTrack();
- interactor.interact(tre2);
- trh.setFit(tre2, trh.chisquared());
- if (_DEBUG)
- System.out.println("trh error pre MS: " + trh.newTrack().error().toString());
- if (fstat > 0)
- return 10000 + 1000 * fstat + icount;
-
-
- ++icount;
- }
-
-
- //System.out.println ("Forward fit sumde: " + sumde );
-// aida.cloud1D("Forward dedx check:").fill(sumde-RKDebug.Instance().getDeGen());
- if (_DEBUG)
- System.out.println("Completed forward fit");
- return 0;
-
- }
-
- public int fitBackward(HTrack trh) {
- PropDir dir = PropDir.BACKWARD;
-// RKDebug.Instance().setPropDir(dir);
-
- //RKPrintSymMatrix psm = new RKPrintSymMatrix();
-
- // Copy the hits from the track.
- List hits = trh.hits();
-
- // Delete the list of hits from the track.
- while (trh.hits().size() > 0)
- trh.dropHit();
-
- double chold = 0.;
-
- double sumde = 0.;
-
-// RKZot zot = new RKZot(trh);
- boolean zottable = true;
-
- int nc = 0;
- int nz = 0;
- int nu = 0;
- String thishit;
-
-
- ListIterator bkglist = hits.listIterator(hits.size());
- // Loop over hits and fit.
- int icount = 0;
- _DEBUG = false;
- while (bkglist.hasPrevious()) {
- Hit hit = (Hit) bkglist.previous();
-
- if (_DEBUG) {
- System.out.println("******* Fit Backward: ADDING NEW HIT *********** ");
- System.out.println("Fit Backward: Hit " + icount + " is: \n");
- System.out.println("Fit Backward: Before prop: " + trh.newTrack());
- System.out.println("Fit Backward: Propogating to surface : " + hit.surface());
- }
- // System.out.println("Before prop spacepoint: " + hit.surface().spacePoint(trh.newTrack().vector()));
- dir = PropDir.BACKWARD;
- if (icount == 0)
- dir = PropDir.FORWARD;//if it's the first hit, propagate forward.
- // propagate to the surface
- PropStat pstat = trh.propagate(_pprop, hit.surface(), dir);
-
- if (_DEBUG) {
- System.out.println("Fit Backward: pstat = " + pstat);
- System.out.println("Fit Backward: After prop: " + trh.newTrack());
- }
-
- if (!pstat.success()) {
- if (_DEBUG) {
- System.out.println("*****************************************************************************");
- System.out.println("FullFitKalman::fitBackward pstat =" + pstat.toString() + " for hit# " + icount);
- System.out.println("FullFitKalman::fitBackward hit = " + hit.toString());
- System.out.println("FullFitKalman::fitBackward track = " + trh.toString());
- System.out.println("FullFitKalman::fitBackward propagator = " + _pprop.toString());
- }
- return -666;
- }
-// System.out.println("trh= \n"+trh+", hit= \n"+hit);
-// System.out.println("_addfit= "+_addfit);
-// System.out.println("After prop spacepoint: " + hit.cluster().surface().spacePoint(trh.newTrack().vector()));
-// System.out.println("trh surface: " + trh.newTrack().surface().toString());
-// System.out.println("hit surface" + hit.surface().toString());
-
-// fit track
- if (_DEBUG)
- System.out.println("Fit Backward: Predicted Hit Vector before addHit: " + hit.predictedVector());
- int fstat = _addfit.addHit(trh, hit);
- if (icount > 2 && trh.chisquared() < 0.1)
- _DEBUG = true;
- if (_DEBUG)
- System.out.println("Fit Backward: Predicted Hit Vector after addHit: " + hit.predictedVector());
-//
- if (_DEBUG) {
- System.out.println("Fit Backward: Measured Hit Vector after addHit: " + hit.cluster());
- System.out.println("Fit Backward: After addhit: " + fstat + " chi^2 = "
- + trh.chisquared() + "\n" + trh.newTrack());
-// System.out.println("After addhit spacepoint: " + hit.cluster().surface().spacePoint(trh.newTrack().vector()));
-//
- }
-
- // Multiple scattering--this is a silly way to do it..
- if (isEven(icount)) {
- ThinXYPlaneMs interactor = new ThinXYPlaneMs(0.01);
-
-
- if (_DEBUG)
- System.out.println("Fit Backward: trh error pre MS: " + trh.newTrack().error().toString());
- ETrack tre2 = trh.newTrack();
- interactor.interact(tre2);
- trh.setFit(tre2, trh.chisquared());
- if (_DEBUG)
- System.out.println("Fit Backward: trh error post MS: " + trh.newTrack().error().toString());
- if (fstat > 0) {
- System.out.println("FullFitKalman::fitBackward fstat= = " + fstat);
- return 10000 + 1000 * fstat + icount;
- }
- }
- icount++;
- }
-
- return 0;
- }/// }
-
- private int interact(HTrack trh, Hit hit, PropDir dir) {
-
- if (hit.surface().pureType().equals(SurfCylinder.staticType())) {
-
- SurfCylinder s = (SurfCylinder) hit.surface();
- double r = s.radius();
- if (doMs) {
- TrackError eold = trh.newTrack().error();
-
-// aida.histogram1D("/Bugs/Fit/Fit scat radius:",300,0.,150.).fill(r);
-
- double l_over_radl = 0.;
- if (r < 1.3) {
- l_over_radl = 0.006136;
- } else if (r < 10.) {
- l_over_radl = 0.000916;
- } else {
- l_over_radl = 0.013747;
- }
-
- ThinCylMs scat = new ThinCylMs(l_over_radl * RKDebug.Instance().getMsFac());
- ETrack et = trh.newTrack();
- ETrack ets = new ETrack(et);
- double chnew = trh.chisquared();
-
- scat.interact(et);
- hit.update(et);
- trh.setFit(et, chnew);
-
-
- for (int i = 0; i < 5; ++i) {
- double ex1 = et.error().get(i, i);
- double ex2 = ets.error().get(i, i);
- double sigsq = ex1 - ex2;
- double pull = -9.;
- double sig = -1.;
- if (sigsq > 0.) {
- sig = Math.sqrt(sigsq);
- pull = (et.vector(i) - ets.vector(i)) / sig;
- }
-// aida.cloud1D("/Bugs/Fit/Forward Fit Delta param:"+i).fill(et.vector(i)-ets.vector(i));
-// if ( sig > 0. ) aida.cloud1D("/Bugs/Fit/Forward Fit Delta error:"+i).fill(sig);
- }
-// System.out.println( "Error after: " + trh.newTrack().error().minus(eold) );
-// */
- } // end if MS enabled
-//
-//
- } // end CYlinder MS
-//
- if (hit.surface().pureType().equals(SurfZPlane.staticType())) {
-
- SurfZPlane s = (SurfZPlane) hit.surface();
- double z = s.z();
- if (doMs) {
- TrackError eold = trh.newTrack().error();
-
-//// aida.histogram1D("/Bugs/Fit/Fit scat z forward:",300,-150,150.).fill(z);
-
- double l_over_radl = (Math.abs(z) < 25) ? 0.000916 : 0.013747;
- ThinZPlaneMs scat = new ThinZPlaneMs(l_over_radl * RKDebug.Instance().getMsFac());
- ETrack et = trh.newTrack();
- ETrack ets = new ETrack(et);
- double chnew = trh.chisquared();
-
- scat.interact(et);
- hit.update(et);
- trh.setFit(et, chnew);
-
- } // end if MS enabled.
-
- } // end ZPlane MS
-
-// if ( hit.surface().pureType().equals(SurfXYPlane.staticType()) ){
-//
-// SurfXYPlane s = (SurfXYPlane) hit.surface();
-// double z = s.z();
-// if ( doMs ){
-// TrackError eold = trh.newTrack().error();
-//
-////// aida.histogram1D("/Bugs/Fit/Fit scat z forward:",300,-150,150.).fill(z);
-//
-// double l_over_radl = ( Math.abs(z)< 25) ? 0.000916 : 0.013747;
-// ThinZPlaneMs scat = new ThinZPlaneMs(l_over_radl*RKDebug.Instance().getMsFac());
-// ETrack et = trh.newTrack();
-// ETrack ets = new ETrack(et);
-// double chnew = trh.chisquared();
-//
-// scat.interact(et);
-// hit.update(et);
-// trh.setFit(et,chnew);
-//
-// } // end if MS enabled.
-//
-// } // end ZPlane MS
-
- // Successful return;
- return 0;
- }
-
- // There is no hit at this site so we only need to do the interaction, not update hits.
- private int interactonly(HTrack trh, double r, double l_over_radl) {
-
- ThinCylMs scat = new ThinCylMs(l_over_radl * RKDebug.Instance().getMsFac());
- ETrack et = trh.newTrack();
- double chnew = trh.chisquared();
-
- scat.interact(et);
- trh.setFit(et, chnew);
-
- // Successful return;
- return 0;
- }
-
-// /**
-// *output stream
-// *
-// * @return The String representation of this instance.
-// */
- public String toString() {
- return getClass().getName();
- }
-
- private boolean isEven(int num) {
- return num % 2 == 0;
- }
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixParamCalculator.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/HelixParamCalculator.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,206 +0,0 @@
-/*
- * HelixParameterCalculator.java
- *
- * Created on July 7th, 2008, 11:09 AM
- *
- *
- */
-package org.hps.recon.tracking.kalman;
-
-import hep.physics.vec.BasicHep3Vector;
-
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.MCParticle;
-
-/**
- * Class used for calculating MC particle track paramters
- * @author Pelham Keahey
- *
- */
-public class HelixParamCalculator {
-
- MCParticle mcp;
- private double R,BField,theta,arclength;
- //Some varibles are usesd in other calculations, thus they are global
- /*
- xc, yc --coordinates
- mcdca -- MC distance of closest approach
- mcphi0 --azimuthal angle
- tanL -- Slope SZ plane ds/dz
- x0,y0 are the position of the particle at the dca
- */
- private double xc,yc,mcdca,mcphi0,tanL;
- private double x0,y0;
- /**
- * Constructor that is fed a magnetic field and MCPARTICLE
- * @param mcpc
- * @param cBField
- */
- public HelixParamCalculator(MCParticle mcpc,double cBField)
- {
- //mc and event global varibles used for calculation
- mcp=mcpc;
-
- //Returns the MagneticField at point (0,0,0), assumes constant magfield
- BField = cBField;
-
- //Calculate theta, the of the helix projected into an SZ plane, from the z axis
- double px = mcp.getPX();
- double py = mcp.getPY();
- double pz = mcp.getPZ();
- double pt = Math.sqrt(px*px + py*py);
- double p = Math.sqrt(pt*pt + pz*pz);
- double cth = pz / p;
- theta = Math.acos(cth);
-
- //Calculate Radius of the Helix
- R = ((mcp.getCharge())*(mcp.getMomentum().magnitude()*Math.sin(theta))/(.0003*BField));
-
- //Slope in the Dz/Ds sense, tanL Calculation
- tanL = mcp.getPZ()/(Math.sqrt(mcp.getPX()*mcp.getPX()+mcp.getPY()*mcp.getPY()));
-
- //Distance of closest approach Calculation
- xc = mcp.getOriginX() + R * Math.sin(Math.atan2(mcp.getPY(),mcp.getPX()));
- yc = mcp.getOriginY() - R * Math.cos(Math.atan2(mcp.getPY(),mcp.getPX()));
- double xcyc = Math.sqrt(xc*xc + yc*yc);
- if(mcp.getCharge()>0)
- {
- mcdca = R - xcyc;
- }
- else
- {
- mcdca = R + xcyc;
- }
-
-
- //azimuthal calculation of the momentum at the DCA, phi0, Calculation
- mcphi0 = Math.atan2(xc/(R-mcdca), -yc/(R-mcdca));
- if(mcphi0<0)
- {
- mcphi0 += 2*Math.PI;
- }
- //z0 Calculation, z position of the particle at dca
- x0 = -mcdca*Math.sin(mcphi0);
- y0 = mcdca*Math.sin(mcphi0);
- arclength = (((mcp.getOriginX()-x0)*Math.cos(mcphi0))+((mcp.getOriginY()-y0)*Math.sin(mcphi0)));
-
- }
- /**
- * Calculates the B-Field from event
- * @param mcpc
- * @param eventc
- */
- public HelixParamCalculator(MCParticle mcpc,EventHeader eventc)
- {
- this(mcpc,eventc.getDetector().getFieldMap().getField(new BasicHep3Vector(0.,0.,0.)).z());
- }
- /**
- * Return the magneticfield at point 0,0,0
- * @return double BField
- */
- public double getMagField()
- {
- return BField;
- }
- /**
- * Return the radius of the Helix track
- * @return double R
- */
- public double getRadius()
- {
- return R;
- }
- /**
- * Return the theta angle for the projection of the helix in the SZ plane
- * from the z axis
- * @return double theta
- */
- public double getTheta()
- {
- return theta;
- }
- /**
- * Return the particle's momentum
- * @return double mcp momentum
- */
- public double getMCMomentum()
- {
- return mcp.getMomentum().magnitude();
- }
- /**
- * Return the curvature (omega)
- * @return double omega
- */
- public double getMCOmega()
- {
- return mcp.getCharge()/((mcp.getMomentum().magnitude()*Math.sin(theta))/(.0003*BField));
- }
- /**
- * Return the transvers momentum of the MC particle, Pt
- * @return double Pt
- */
- public double getMCTransverseMomentum()
- {
- return (mcp.getMomentum().magnitude())*Math.sin(theta);
- }
- /**
- * Return the slope of the helix in the SZ plane, tan(lambda)
- * @return double tanL
- */
- public double getSlopeSZPlane()
- {
- return tanL;
- }
- /**
- * Return the distance of closest approach
- * @return double mcdca
- */
- public double getDCA()
- {
- return mcdca;
- }
- /**
- * Return the azimuthal angle of the momentum when at the position of closest approach
- * @return double mcphi0
- */
- public double getPhi0()
- {
- return mcphi0;
- }
- /**
- * Return the z position at the distance of closest approach
- * @return double z0 position
- */
- public double getZ0()
- {
- double x0 = -mcdca*Math.sin(mcphi0);
- double y0 = mcdca*Math.sin(mcphi0);
- double s = (((mcp.getOriginX()-x0)*Math.cos(mcphi0))+((mcp.getOriginY()-y0)*Math.sin(mcphi0)));
- return mcp.getOriginZ()-(s*tanL);
- }
- /**
- * Return the arclength of the helix from the ORIGIN TO THE DCA
- * @return double arclength
- */
- public double getArcLength()
- {
- return arclength;
- }
- /**
- * Return the x position of the particle when at the dca
- * @return double arclength
- */
- public double getX0()
- {
- return x0;
- }
- /**
- * Return the y position of the particle at the dca
- * @return double arclength
- */
- public double getY0()
- {
- return y0;
- }
-
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanFilterDriver.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanFilterDriver.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,132 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-package org.hps.recon.tracking.kalman;
-
-import java.util.List;
-
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.Track;
-import org.lcsim.event.TrackerHit;
-import org.lcsim.fit.helicaltrack.HelicalTrackFit;
-import org.lcsim.geometry.Detector;
-import org.lcsim.recon.tracking.trfbase.ETrack;
-import org.lcsim.recon.tracking.trfbase.Propagator;
-import org.lcsim.recon.tracking.trfbase.TrackError;
-import org.lcsim.recon.tracking.trfbase.VTrack;
-import org.lcsim.recon.tracking.trffit.HTrack;
-import org.lcsim.util.Driver;
-
-/**
- *
- * @author ecfine ([log in to unmask])
- */
-
-
-/* Takes Matched Tracks from the Track Reconstruction Driver, and makes TRF
- * HTracks. Adds hits from the original tracks to the new HTrack, and runs a
- * Kalman filter over the HTrack.
- *
- * This only does a forward fit, and multiple scattering is only included at
- * interacting planes, and just by assuming that each plane is .01 radiation
- * lengths. Energy loss is not accounted for. Additionally, the method for
- * constructing hits assumes that every hit occurs at an XY plane, while ideally,
- * there would be a method which checks the lcsim detector geometry and then
- * decides what kind of surface to model the shape as. There may be some methods
- * in the ShapeHelper interface (particularly in TrdHelper) that would be useful
- * for this, but nothing completed. Additionally, to run realistic multiple
- * scattering with non-interacting detector elements, there would need to be a
- * way to find the intercepts between a specific track and the detector elements.
- *
- * Also, magnetic field is just set at 1.0 in each class. It should be taken from
- * the detector geometry. */
-
-public class KalmanFilterDriver extends Driver{
- ShapeDispatcher shapeDis = new ShapeDispatcher();
- TrackUtils trackUtils = new TrackUtils();
- Propagator prop = null;
- FullFitKalman fitk = null;
- KalmanGeom geom = null;
- Detector detector = null;
- HTrack ht = null;
- double bz = 0.5;
-
-
-
- public void detectorChanged(Detector det){
- detector = det;
- geom = new KalmanGeom(detector); // new geometry containing detector info
- prop = geom.newPropagator();
- System.out.println("geom field = " + geom.bz + ", trackUtils field = " + trackUtils.bz);
-// trackUtils.setBZ(geom.bz);
- fitk = new FullFitKalman(prop);
-//
- }
-
-
- public void process(EventHeader event){
- /* Get the tracklist for each event, and then for each track
- * get the starting track parameters and covariance matrix for an
- * outward fit from the seedtracker. */
- if (event.hasItem("MatchedTracks")) {
- List<Track> trklist = (List<Track>) event.get("MatchedTracks");
- System.out.println("number of tracks: " + trklist.size());
- for (int i = 0; i < trklist.size(); i++) {
- /* Start with a HelicalTrackFit, turn it into a VTrack,
- * turn that into an ETrack, and turn that into an HTrack.
- * Then add detector hits from the original track. */
- if(trklist.get(i).getTrackerHits().size()<4) {
- System.out.println("Continue, this track has only " + trklist.get(i).getTrackerHits().size());
- continue;
- }
- System.out.println("Making tracks...");
- Track track = trklist.get(i);
- HelicalTrackFit helicalTrack = shapeDis.trackToHelix(track);
- VTrack vt = trackUtils.makeVTrack(helicalTrack);
- TrackError initialError = trackUtils.getInitalError(helicalTrack);
- ETrack et = new ETrack (vt, initialError);
- ht = new HTrack(et);
-
- /* Add hits from original track */
- for (int k = 0; k < track.getTrackerHits().size(); k++) {
- TrackerHit thit = track.getTrackerHits().get(k);
- System.out.println("Adding hit...");
- //ht = geom.addTrackerHit(thit, ht, helicalTrack, vt);
- // phansson
- // removing unused arguments to avoid confusion (see function for details)
- ht = geom.addTrackerHit(thit, ht);
- }
-
-
-
- /* Once we have an HTrack with the ordered list of hits, we pass
- * this to the Kalman fitter. */
- System.out.println("Running Kalman fit...");
- int fstarf = fitk.fit(ht);
-// System.out.println("geom field = " + geom.bz + ", trackUtils field = " + trackUtils.bz);
-// List<RawTrackerHit> rawHits = event.get(RawTrackerHit.class, "RawTrackerHitMaker_RawTrackerHits");
-// System.out.println("SimTrackerHits info: ");
-// for(int j = 0; j < rawHits.size(); j++){
-// List<SimTrackerHit> simHits = rawHits.get(j).getSimTrackerHits();
-// for(int l = 0; l < simHits.size(); l++){
-// System.out.println("point = [" + simHits.get(l).getPoint()[0]
-// + ", " + simHits.get(l).getPoint()[1] +
-// ", " + simHits.get(l).getPoint()[2] +
-// "], momentum = [" + simHits.get(l).getMomentum()[0]
-// + ", " + simHits.get(l).getMomentum()[1]
-// +", " + simHits.get(l).getMomentum()[2] + "]");
-// }
-// }
- }
-
-
- }
-// } else {
-// System.out.println("No tracks!");
-// }
-
- }
-
-
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanGeom.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanGeom.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,743 +0,0 @@
-package org.hps.recon.tracking.kalman;
-
-import hep.physics.vec.BasicHep3Matrix;
-import hep.physics.vec.Hep3Matrix;
-import hep.physics.vec.Hep3Vector;
-import hep.physics.vec.VecOp;
-
-import java.util.ArrayList;
-import java.util.List;
-
-import org.hps.recon.tracking.HPSTransformations;
-import org.hps.recon.tracking.kalman.util.PropDCAZ;
-import org.lcsim.detector.IDetectorElement;
-import org.lcsim.detector.ILogicalVolume;
-import org.lcsim.detector.IPhysicalVolume;
-import org.lcsim.detector.ITransform3D;
-import org.lcsim.detector.solids.ISolid;
-import org.lcsim.detector.solids.Point3D;
-import org.lcsim.detector.tracker.silicon.ChargeCarrier;
-import org.lcsim.detector.tracker.silicon.SiSensor;
-import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
-import org.lcsim.event.RawTrackerHit;
-import org.lcsim.event.TrackerHit;
-import org.lcsim.fit.helicaltrack.HelicalTrackCross;
-import org.lcsim.fit.helicaltrack.HelicalTrackFit;
-import org.lcsim.fit.helicaltrack.HelicalTrackStrip;
-import org.lcsim.geometry.Detector;
-import org.lcsim.recon.tracking.trfbase.ETrack;
-import org.lcsim.recon.tracking.trfbase.Hit;
-import org.lcsim.recon.tracking.trfbase.PropDispatch;
-import org.lcsim.recon.tracking.trfbase.Propagator;
-import org.lcsim.recon.tracking.trfbase.VTrack;
-import org.lcsim.recon.tracking.trfcyl.PropCyl;
-import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
-import org.lcsim.recon.tracking.trfcylplane.PropCylZ;
-import org.lcsim.recon.tracking.trfcylplane.PropZCyl;
-import org.lcsim.recon.tracking.trfdca.PropCylDCA;
-import org.lcsim.recon.tracking.trfdca.PropDCACyl;
-import org.lcsim.recon.tracking.trfdca.SurfDCA;
-import org.lcsim.recon.tracking.trffit.HTrack;
-import org.lcsim.recon.tracking.trfxyp.ClusXYPlane1;
-import org.lcsim.recon.tracking.trfxyp.PropXYXY;
-import org.lcsim.recon.tracking.trfxyp.SurfXYPlane;
-import org.lcsim.recon.tracking.trfzp.PropZZ;
-import org.lcsim.recon.tracking.trfzp.SurfZPlane;
-
-/**
- *
- * Extract needed information from the geometry system and serve it
- * to the callers in a convenient form.
- *
- *
- *@author $Author: phansson $
- *@version $Id: KalmanGeom.java,v 1.6 2013/10/15 00:33:54 phansson Exp $
- *
- * Date $Date: 2013/10/15 00:33:54 $
- *
- */
-/* To make the Kalman filter work for any detector, the type of hit that is
- * added to the HTrack has to be determined based on the lcsim detector
- * geometry. Additionally, this class has to be able to find the intersections
- * between a track and the detector elements to be able to implement more
- * realistic multiple scattering. Currently this class does not support
- * any of this functionality, and just adds hits. However there are methods for
- * looping through detector elements, etc. that may be useful for finding
- * intersections at some point. (ecfine) */
-public class KalmanGeom {
-
- // TRF wants distances in cm, not mm.
- private double mmTocm = 0.1;
- // this flag is a temporary fix to avoid making a surface behind the origin in x.
- private double flag = 0;
-// // private Array radius = null;
-// private Subdetector sd_tbar = null;
-// private Subdetector sd_vbar = null;
-// private Subdetector sd_tec = null;
-// private Subdetector sd_vec = null;
-// private Subdetector sd_tfor = null;
-// private Subdetector sd_bpipe = null;
-//
-// private IDetectorElement de_tbar = null;
-// private IDetectorElement de_vbar = null;
-// private IDetectorElement de_tec = null;
-// private IDetectorElement de_vec = null;
-// private IDetectorElement de_tfor = null;
- // Nominal magnetic field.
- public double bz = 0.5;
- // Lists of interesting surfaces.
- public List<KalmanSurface> Surf = new ArrayList<KalmanSurface>();
- KalmanSurface surf = null;
- // Specific surface type A to surface type B propagators.
- private PropCyl propcyl = null;
- private PropZZ propzz = null;
- private PropDCACyl propdcacyl = null;
- private PropCylDCA propcyldca = null;
- private PropZCyl propzcyl = null;
- private PropCylZ propcylz = null;
- private PropXYXY propxyxy = null;
- private PropDCAZ propdcaz = null;
- private PropDCAXY propdcaxy = null;
- private PropXYDCA propxydca = null;
-// private PropZDCA propzdca = null;
- // The master propagator that can go between any pair of surfaces.
- private PropDispatch pDispatch = null;
- // The run time configuration system.
- //ToyConfig config;
- // Information from the run time configuration system.
-// double respixel;
-// double resstrip;
-// boolean trackerbarrel2d;
-// double zres2dTrackerBarrel;
-// boolean vtxFwdEquivStrips;
- // Does this detector have forward tracking.
- boolean hasforward = false;
- // Stuff for adding surfaces
- static ArrayList physicalVolumes = new ArrayList();
- ILogicalVolume logical;
- ShapeDispatcher shapeDispatcher = new ShapeDispatcher();
- Detector detector = null;
-
- public KalmanGeom(Detector det) {
- detector = det;
- System.out.println("New detector: " + detector.getName());
- logical = detector.getTrackingVolume().getLogicalVolume();
-
-
- // Extract information from the run time configuration system.
- // try{
- // ToyConfig config = ToyConfig.getInstance();
- // respixel = config.getDouble("respixel");
- // resstrip = config.getDouble("resstrip");
- // trackerbarrel2d = config.getBoolean("trackerbarrel2d");
- // zres2dTrackerBarrel = config.getDouble("zres2dTrackerBarrel");
- // vtxFwdEquivStrips = config.getBoolean("vtxFwdEquivStrips");
- //
- // } catch (ToyConfigException e){
- // System.out.println (e.getMessage() );
- // System.out.println ("Stopping now." );
- // System.exit(-1);
- // }
-
-
-
- // Map<String, Subdetector> subDetMap = detector.getSubdetectors();
- //
- // Subdetector sd_tbar = subDetMap.get("TrackerBarrel");
- // Subdetector sd_vbar = subDetMap.get("VertexBarrel");
- // Subdetector sd_tec = subDetMap.get("TrackerEndcap");
- // Subdetector sd_vec = subDetMap.get("VertexEndcap");
- // Subdetector sd_tfor = subDetMap.get("TrackerForward");
- // Subdetector sd_bpipe = subDetMap.get("BeamPipe");
- // System.out.println ("Checking .... " + sd_tbar + " | " + sd_tfor);
- // Don't use subdetectors anymore...
-
- // Check for forward tracking system.
- // if ( sd_tfor == null ) {
- // System.out.println ("Checking 1 .... " );
- // if ( detector.getName().compareTo("sid00") != 0 ){
- // System.out.println("Expected to find a TrackerForward Subdetector but did not!");
- // System.exit(-1);
- // }
- // }else{
- // System.out.println ("Checking 2 .... " );
- //
- // hasforward = true;
- // } // I may want to implement something like this at some point but not now.
-
- // if ( hasforward ){
- // de_tfor = sd_tfor.getDetectorElement();
- // }
-
- /* Cycle through detector and add all surfaces to the surface list. */
- // addAllSurfaces();
-
- System.out.println("Number of surfaces: " + Surf.size());
- for (KalmanSurface surf : Surf) {
- surf.Print();
- }
-
-
- double[] origin = {0., 0., 0.};
-
-// Map<String,Field> fields = detector.getFields();
-// Set<String> keys = fields.keySet();
-// for ( String key : keys ){
-// Field field = fields.get(key);
-// String classname = field.getClass().getName();
-// String shortname = classname.replaceAll( "org.lcsim.geometry.field.", "");
-// if ( shortname.compareTo("Solenoid") != 0 ){
-// System.out.println("Expected, but did not find, a solenoid: " + shortname );
-//// System.exit(-1);
-// }
-// Solenoid s = (Solenoid)field;
-// bz = s.getInnerField()[2];
-// if ( bz == 0. ){
-// System.out.println("This code will not work with a magnetic field of 0: " + shortname );
-//// System.exit(-1);
-// }
-// break;
-// }
-
- // Instantiate surface-pair specific propagators.
- propcyl = new PropCyl(bz);
- propzz = new PropZZ(bz);
- propdcacyl = new PropDCACyl(bz);
- propcyldca = new PropCylDCA(bz);
- propzcyl = new PropZCyl(bz);
- propcylz = new PropCylZ(bz);
- propxyxy = new PropXYXY(bz);
- propdcaxy = new PropDCAXY(bz);
- propxydca = new PropXYDCA(bz);
- // propdcaz = new PropDCAZ(bz);
- // propzdca = new PropZDCA(bz);
- // Again, need xy plane propagators!
-
- // Instantiate and configure the general purpose propagator.
- pDispatch = new PropDispatch();
- pDispatch.addPropagator(SurfZPlane.staticType(), SurfZPlane.staticType(), propzz);
- pDispatch.addPropagator(SurfCylinder.staticType(), SurfCylinder.staticType(), propcyl);
- pDispatch.addPropagator(SurfDCA.staticType(), SurfCylinder.staticType(), propdcacyl);
- pDispatch.addPropagator(SurfCylinder.staticType(), SurfDCA.staticType(), propcyldca);
- pDispatch.addPropagator(SurfZPlane.staticType(), SurfCylinder.staticType(), propzcyl);
- pDispatch.addPropagator(SurfCylinder.staticType(), SurfZPlane.staticType(), propcylz);
- pDispatch.addPropagator(SurfXYPlane.staticType(), SurfXYPlane.staticType(), propxyxy);
- pDispatch.addPropagator(SurfDCA.staticType(), SurfXYPlane.staticType(), propdcaxy);
- pDispatch.addPropagator(SurfXYPlane.staticType(), SurfDCA.staticType(), propxydca);
- // pDispatch.addPropagator( SurfDCA.staticType(), SurfZPlane.staticType(), propdcaz);
- // pDispatch.addPropagator( SurfZPlane.staticType(), SurfDCA.staticType(), propzdca);
-
- }
-
- public Propagator newPropagator() {
- // Clone not supported for pDispatch.
- //return pDispatch.newPropagator();
- return (Propagator) pDispatch;
- }
-
-// public List<RKSurf> getCylinders(){
-// return Surf;
-// }
-// public List<RKSurf> getZplus(){
-// return ZSurfplus;
-// }
-//
-// public List<RKSurf> getZMinus(){
-// return ZSurfminus;
-// }
- public double getBz() {
- return bz;
- }
-
- public void setBz(double bfield) {
- bz = bfield;
- }
-
-// Return a list of z surfaces, going foward along the track.
-// public List<RKSurf> getZ( double z0, double cz ){
-// List<RKSurf> zlist = new ArrayList<RKSurf>();
-// if ( cz > 0 ){
-// for ( Iterator ihit=ZSurfAll.iterator(); ihit.hasNext(); ){
-// RKSurf s = (RKSurf) ihit.next();
-// if ( s.zc >= z0 ){
-// zlist.add(s);
-// }
-// }
-// } else if ( cz < 0 ){
-// for ( ListIterator ihit=ZSurfAll.listIterator(ZSurfAll.size());
-// ihit.hasPrevious(); ){
-// RKSurf s = (RKSurf) ihit.previous();
-// if ( s.zc <= z0 ){
-// zlist.add(s);
-// }
-// }
-// }
-//
-// return zlist;
-// }
- // Adds all surfaces to the surface list
- private void addAllSurfaces() {
- addDaughterSurfaces(logical);
- }
-
- private void addDaughterSurfaces(ILogicalVolume logical) {
- if (logical.getDaughters().size() == 0) {
- // need to avoid making surface for target. for now, flag..
- flag++;
- if (flag > 5) {
- surf = shapeDispatcher.getKalmanSurf(logical.getSolid());
- Surf.add(surf);
- }
- }
- for (int n = 0; n < logical.getNumberOfDaughters(); n++) {
- IPhysicalVolume physical = logical.getDaughter(n);
- addDaughterSurfaces(physicalToLogical(physical));
- physicalVolumes.remove(physicalVolumes.size() - 1);
- }
- }
-
- // Given a point, find the solid that contains it. Doesn't work.
- private ISolid findSolidFromPoint(Point3D hitPoint) {
- ISolid solid = checkDaughterSurfaces(logical, hitPoint);
- return solid;
- }
-
- private ISolid checkDaughterSurfaces(ILogicalVolume logical, Point3D hitPoint) {
- if (logical.getDaughters().size() == 0) {
- if (pointIsOnSolid(logical.getSolid(), hitPoint)) {
- return logical.getSolid();
- }
- }
- for (int n = 0; n < logical.getNumberOfDaughters(); n++) {
- IPhysicalVolume physical = logical.getDaughter(n);
- checkDaughterSurfaces(physicalToLogical(physical), hitPoint);
- physicalVolumes.remove(physicalVolumes.size() - 1);
- }
- System.out.print("This hit isn't on a solid!");
- return null;
- }
-
- // Given a TrackerHit, return the surface the hit should be placed on. Doesn't work.
- private KalmanSurface findTrackerHitSurface(TrackerHit thit) {
- double[] position = thit.getPosition();
- Point3D hitPoint = new Point3D(position[0], position[1], position[2]);
- ISolid hitSolid = findSolidFromPoint(hitPoint);
- KalmanSurface hitSurf = getKSurfFromSolid(hitSolid);
- return hitSurf;
- }
-
- private KalmanSurface getKSurfFromSolid(ISolid hitSolid) {
- throw new UnsupportedOperationException("Not yet implemented");
- }
-
- // Adds the physical volume to the physicalVolumes ArrayList.
- private ILogicalVolume physicalToLogical(IPhysicalVolume physical) {
- physicalVolumes.add(physical.getTransform());
- return physical.getLogicalVolume();
- }
-
- // Adds intercepts between a HTrack and the detector. Doesn't work.
- public void addIntercepts(HTrack ht, Detector detector) {
- logical = detector.getTrackingVolume().getLogicalVolume();
- addDaughterIntercepts(logical, ht);
- }
-
- private void addDaughterIntercepts(ILogicalVolume logical, HTrack ht) {
- if (logical.getDaughters().size() == 0) {
- ht = shapeDispatcher.addIntercept(logical.getSolid(), ht);
- }
- for (int n = 0; n < logical.getNumberOfDaughters(); n++) {
- IPhysicalVolume physical = logical.getDaughter(n);
- addDaughterIntercepts(physicalToLogical(physical), ht);
- physicalVolumes.remove(physicalVolumes.size() - 1);
- }
- }
-
- // Add a hit from an lcsim Track to a trf HTrack.
- public HTrack addTrackerHit(TrackerHit thit, HTrack htrack, HelicalTrackFit track,
- VTrack vtrack) {
- /* This should check the kind of solid that got hit and model it as the
- * correct surface based on the solid dimensions. For now it just
- * assumes XY planes. */
-
- double[] position = thit.getPosition();
- if (position.length != 3) {
- System.out.println("Position has more than 3 coordinates?!");
- }
-
-
-// // For checking surfaces
-// KalmanSurface surface= findTrackerHitSurface(thit);
-// PropStat prop = new PropStat();
-// prop.setForward();
-// double path = (position[2] - track.z0())/track.slope();
-// KalmanHit khit = new KalmanHit(track, surface, vtrack, prop, path);
-// Hit hit = khit.MakeHit();
-// Point3D hitPoint = new Point3D(position[0], position[1], position[2]);
-// ISolid hitSolid = findSolidFromPoint(hitPoint);
-
-
- //PELLE:
- // The HelicalTrackFit track and VTrack vtrack are not used
-
-
-
- // Need an ETrack to add hits from clusters.
- ETrack tre = htrack.newTrack();
-
- // This also assumes the hit in question is a Helical Track Cross.
- HelicalTrackCross htrackcross = (HelicalTrackCross) thit;
- HelicalTrackStrip strip1 = (HelicalTrackStrip) htrackcross.getStrips().get(0);
- HelicalTrackStrip strip2 = (HelicalTrackStrip) htrackcross.getStrips().get(1);
-
- double dist1 = dotProduct(strip1.origin(), strip1.w());
- double dist2 = dotProduct(strip2.origin(), strip2.w());
-
-
- double phi1 = Math.atan2(strip1.w().y() * dist1 * strip1.origin().y(), strip1.w().x() * dist1 * strip1.origin().x());
- double phi2 = Math.atan2(strip2.w().y() * dist2 * strip2.origin().y(), strip2.w().x() * dist2 * strip2.origin().x());
- // phi needs to be between 0 and 2PI. Additionally, it seemed like
- // rounding errors caused problems for very small negative values of phi,
- // so these are approximated as 0.
- if (phi1 < 0) {
- if (phi1 > -0.00000001) {
- phi1 = 0;
- } else {
- phi1 = phi1 + 2 * Math.PI;
- }
- }
- if (phi2 < 0) {
- if (phi2 > -0.00000001) {
- phi2 = 0;
- } else {
- phi2 = phi2 + 2 * Math.PI;
- }
- }
-
- System.out.println("Origin of strip1: " + strip1.origin());
- System.out.println("u,v,w of strip1 : " + strip1.u().toString() + "," + strip1.v().toString() + "," + strip1.w().toString());
- System.out.println("Origin of strip2: " + strip2.origin());
- System.out.println("u,v,w of strip2 : " + strip2.u().toString() + "," + strip2.v().toString() + "," + strip2.w().toString());
- System.out.println("dist1: " + dist1 + ", phi1: " + phi1);
- System.out.println("dist2: " + dist2 + ", phi2: " + phi2);
- // ClusXYPlane needs wz, wv, avz, phi and dist to create a cluster.
- double wz1 = strip1.u().z();
- double wz2 = strip2.u().z();
- double wv1 = strip1.u().x() * (-Math.sin(phi1)) + strip1.u().y() * Math.cos(phi1);
- double wv2 = strip2.u().x() * (-Math.sin(phi2)) + strip2.u().y() * Math.cos(phi2);
-
-
- double avz1 = strip1.umeas()
- - (dist1 * (Math.cos(phi1) * strip1.u().x() + Math.sin(phi1) * strip1.u().y())
- - (strip1.origin().x() * strip1.u().x() + strip1.origin().y() * strip1.u().y()
- + strip1.origin().z() * strip1.u().z()));
- double avz2 = strip2.umeas()
- - (dist2 * (Math.cos(phi2) * strip2.u().x() + Math.sin(phi2) * strip2.u().y())
- - (strip2.origin().x() * strip2.u().x() + strip2.origin().y() * strip2.u().y()
- + strip2.origin().z() * strip2.u().z()));
- double davz1 = strip1.du();
- double davz2 = strip2.du();
-
- System.out.println("avz1 = " + avz1 + " = " + strip1.umeas() + " - " + (dist1 * (Math.cos(phi1) * strip1.u().x() + Math.sin(phi1) * strip1.u().y())
- - (strip1.origin().x() * strip1.u().x() + strip1.origin().y() * strip1.u().y()
- + strip1.origin().z() * strip1.u().z())));
- System.out.println("avz2 = " + avz2 + " = " + strip2.umeas() + " - " + (dist2 * (Math.cos(phi2) * strip2.u().x() + Math.sin(phi2) * strip2.u().y())
- - (strip2.origin().x() * strip2.u().x() + strip2.origin().y() * strip2.u().y()
- + strip2.origin().z() * strip2.u().z())));
-
-
-
-// // Just to check and make sure Vtrf makes sense:
-// if((dist1*Math.cos(phi1) - strip1.origin().x())*Math.cos(phi1) ==
-// (dist1*Math.sin(phi1) - strip1.origin().y())*Math.sin(phi1)){
-// System.out.println("Vtrf1 coherent!");
-// } else{
-// System.out.println("Vtrf1 doesn't make sense :(");
-// }
-// if((dist2*Math.cos(phi2) - strip2.origin().x())*Math.cos(phi2) ==
-// (dist2*Math.sin(phi2) - strip2.origin().y())*Math.sin(phi2)){
-// System.out.println("Vtrf2 coherent!");
-// } else{
-// System.out.println("Vtrf2 doesn't make sense :(");
-// }
-
- System.out.println("wz1 = " + wz1 + ", wv1 = " + wv1 + ", avz1 = " + avz1 + ", davz1 =" + davz1);
- System.out.println("wz2 = " + wz2 + ", wv2 = " + wv2 + ", avz2 = " + avz2 + ", davz2 =" + davz2);
-
- // Create new clusters and get hit predictions.
- ClusXYPlane1 cluster1 = new ClusXYPlane1(dist1, phi1, wv1, wz1, avz1, davz1);
- ClusXYPlane1 cluster2 = new ClusXYPlane1(dist2, phi2, wv2, wz2, avz2, davz2);
- List hits1 = cluster1.predict(tre);
- List hits2 = cluster2.predict(tre);
- Hit hit1 = (Hit) hits1.get(0);
- Hit hit2 = (Hit) hits2.get(0);
-
- System.out.println("hit1 predicted vector: " + hit1.predictedVector());
- System.out.println("hit2 predicted vector: " + hit2.predictedVector());
- System.out.println("hit position from trackerhit: [" + thit.getPosition()[0] + ", "
- + thit.getPosition()[1] + ", " + thit.getPosition()[2] + "]");
-
- hit1.setParentPointer(cluster1);
- hit2.setParentPointer(cluster2);
- htrack.addHit(hit1);
- htrack.addHit(hit2);
-
- return htrack;
-
- }
-
-
- // Add a hit from an lcsim Track to a trf HTrack.
- public HTrack addTrackerHit(TrackerHit thit, HTrack htrack) {
- /* This should check the kind of solid that got hit and model it as the
- * correct surface based on the solid dimensions. For now it just
- * assumes XY planes. */
-
- System.out.println("\nAdd Tracker Hit at position " + thit.getPosition()[0] + "," + thit.getPosition()[1] + "," + thit.getPosition()[2]);
-
- //PELLE:
- // The HelicalTrackFit track and VTrack vtrack are not used so was removed
- // to avoid confusion.
- // They might be needed to be able to figure out what solid to
- // to extrapolate between if this gets implemented.
-
-
-
- double[] position = thit.getPosition();
- if (position.length != 3) {
- System.out.println("Position has more than 3 coordinates?!");
- }
-
-
-// // For checking surfaces
-// KalmanSurface surface= findTrackerHitSurface(thit);
-// PropStat prop = new PropStat();
-// prop.setForward();
-// double path = (position[2] - track.z0())/track.slope();
-// KalmanHit khit = new KalmanHit(track, surface, vtrack, prop, path);
-// Hit hit = khit.MakeHit();
-// Point3D hitPoint = new Point3D(position[0], position[1], position[2]);
-// ISolid hitSolid = findSolidFromPoint(hitPoint);
-
-
-
-
- // Need an ETrack to add hits from clusters.
- ETrack tre = htrack.newTrack();
-
- // This also assumes the hit in question is a Helical Track Cross.
- HelicalTrackCross htrackcross = (HelicalTrackCross) thit;
- HelicalTrackStrip strip1 = (HelicalTrackStrip) htrackcross.getStrips().get(0);
- HelicalTrackStrip strip2 = (HelicalTrackStrip) htrackcross.getStrips().get(1);
-
- Hep3Vector u1 = strip1.u();
- Hep3Vector u2 = strip2.u();
- Hep3Vector v1 = strip1.v();
- Hep3Vector v2 = strip2.v();
- Hep3Vector w1 = strip1.w();
- Hep3Vector w2 = strip2.w();
-
- Hep3Vector origin1 = strip1.origin();
- Hep3Vector origin2 = strip2.origin();
-
- double u1meas = strip1.umeas();
- double u2meas = strip2.umeas();
-
-
- System.out.println("strip1 origin: " + origin1);
- System.out.println("strip2 origin: " + origin2);
- System.out.printf("strip1 u=%s\tv=%s\tw=%s\n",u1.toString(),v1.toString(),w1.toString());
- System.out.printf("strip2 u=%s\tv=%s\tw=%s\n",u2.toString(),v2.toString(),w2.toString());
- System.out.println("strip1 umeas: " + u1meas);
- System.out.println("strip2 umeas: " + u2meas);
-
- System.out.println("Rotate strip2 u,v,w into strip1 frame");
-
- Hep3Matrix strip2ToTrk = getStripToTrackRotation(strip2);
- System.out.println("Strip2ToTrk matrix " + strip2ToTrk.toString());
- System.out.println("Get the rotation matrix for going from track frame to strip1 frame");
- Hep3Matrix trackToStrip1 = getTrackToStripRotation(strip1);
- Hep3Matrix strip2ToStrip1 = VecOp.mult(trackToStrip1,strip2ToTrk);
-
-
- u2 = VecOp.mult(strip2ToStrip1, u2);
- v2 = VecOp.mult(strip2ToStrip1, v2);
- w2 = VecOp.mult(strip2ToStrip1, w2);
-
- System.out.printf("strip2 u=%s\tv=%s\tw=%s\n",u2.toString(),v2.toString(),w2.toString());
-
-
-
- double dist1 = dotProduct(origin1, w1);
- double dist2 = dotProduct(origin2, w2);
-
-
- double phi1 = Math.atan2(w1.y() * dist1 * origin1.y(), w1.x() * dist1 * origin1.x());
- double phi2 = Math.atan2(w2.y() * dist2 * origin2.y(), w2.x() * dist2 * origin2.x());
- // phi needs to be between 0 and 2PI. Additionally, it seemed like
- // rounding errors caused problems for very small negative values of phi,
- // so these are approximated as 0.
- if (phi1 < 0) {
- if (phi1 > -0.00000001) {
- phi1 = 0;
- } else {
- phi1 = phi1 + 2 * Math.PI;
- }
- }
- if (phi2 < 0) {
- if (phi2 > -0.00000001) {
- phi2 = 0;
- } else {
- phi2 = phi2 + 2 * Math.PI;
- }
- }
- System.out.println("dist1: " + dist1 + ", phi1: " + phi1);
- System.out.println("dist2: " + dist2 + ", phi2: " + phi2);
- // ClusXYPlane needs wz, wv, avz, phi and dist to create a cluster.
- double wz1 = u1.z();
- double wz2 = u2.z();
- double wv1 = u1.x() * (-Math.sin(phi1)) + u1.y() * Math.cos(phi1);
- double wv2 = u2.x() * (-Math.sin(phi2)) + u2.y() * Math.cos(phi2);
-
-
- double avz1 = u1meas
- - (dist1 * (Math.cos(phi1) * u1.x() + Math.sin(phi1) * u1.y())
- - (origin1.x() * u1.x() + origin1.y() * u1.y()
- + origin1.z() * u1.z()));
- double avz2 = u2meas
- - (dist2 * (Math.cos(phi2) * u2.x() + Math.sin(phi2) * u2.y())
- - (origin2.x() * u2.x() + origin2.y() * u2.y()
- + origin2.z() * u2.z()));
- double davz1 = strip1.du();
- double davz2 = strip2.du();
-
- System.out.println("avz1 = " + avz1 + " = " + u1meas + " - " + (dist1 * (Math.cos(phi1) * u1.x() + Math.sin(phi1) * u1.y())
- - (origin1.x() * u1.x() + origin1.y() * u1.y()
- + origin1.z() * u1.z())));
- System.out.println("avz2 = " + avz2 + " = " + u2meas + " - " + (dist2 * (Math.cos(phi2) * u2.x() + Math.sin(phi2) * u2.y())
- - (origin2.x() * u2.x() + origin2.y() * u2.y()
- + origin2.z() * u2.z())));
-
-
-
-// // Just to check and make sure Vtrf makes sense:
-// if((dist1*Math.cos(phi1) - origin1.x())*Math.cos(phi1) ==
-// (dist1*Math.sin(phi1) - origin1.y())*Math.sin(phi1)){
-// System.out.println("Vtrf1 coherent!");
-// } else{
-// System.out.println("Vtrf1 doesn't make sense :(");
-// }
-// if((dist2*Math.cos(phi2) - origin2.x())*Math.cos(phi2) ==
-// (dist2*Math.sin(phi2) - origin2.y())*Math.sin(phi2)){
-// System.out.println("Vtrf2 coherent!");
-// } else{
-// System.out.println("Vtrf2 doesn't make sense :(");
-// }
-
- System.out.println("wz1 = " + wz1 + ", wv1 = " + wv1 + ", avz1 = " + avz1 + ", davz1 =" + davz1);
- System.out.println("wz2 = " + wz2 + ", wv2 = " + wv2 + ", avz2 = " + avz2 + ", davz2 =" + davz2);
-
- // Create new clusters and get hit predictions.
- ClusXYPlane1 cluster1 = new ClusXYPlane1(dist1, phi1, wv1, wz1, avz1, davz1);
- ClusXYPlane1 cluster2 = new ClusXYPlane1(dist2, phi2, wv2, wz2, avz2, davz2);
- List hits1 = cluster1.predict(tre);
- List hits2 = cluster2.predict(tre);
- Hit hit1 = (Hit) hits1.get(0);
- Hit hit2 = (Hit) hits2.get(0);
-
- System.out.println("hit1 predicted vector: " + hit1.predictedVector());
- System.out.println("hit2 predicted vector: " + hit2.predictedVector());
- System.out.println("hit position from trackerhit: [" + thit.getPosition()[0] + ", "
- + thit.getPosition()[1] + ", " + thit.getPosition()[2] + "]");
-
- hit1.setParentPointer(cluster1);
- hit2.setParentPointer(cluster2);
- htrack.addHit(hit1);
- htrack.addHit(hit2);
-
- return htrack;
-
- }
-
-
-
-
-
- private boolean pointIsOnSolid(ISolid solid, Point3D hitPoint) {
- return shapeDispatcher.pointIsOnSolid(solid, hitPoint);
- }
-
- private double dotProduct(Hep3Vector v1, Hep3Vector v2) {
- double dotProduct = v1.x() * v2.x() + v1.y() * v2.y() + v1.z() * v2.z();
- return dotProduct;
- }
-
-
- private Hep3Matrix getStripToTrackRotation(HelicalTrackStrip strip) {
- //This function transforms the vec to the track coordinates that the supplied strip has
-
- //Transform from strip frame (u,v,w) to JLab frame (done through the RawTrackerHit)
- ITransform3D stripToDet = GetLocalToGlobal(strip);
- //Get rotation matrix
- Hep3Matrix stripToDetMatrix = (BasicHep3Matrix) stripToDet.getRotation().getRotationMatrix();
- //Transformation between JLab and tracking coordinates
- Hep3Matrix detToTrackMatrix = (BasicHep3Matrix) HPSTransformations.getMatrix();
-
- if (true) {
- System.out.println("Getting the rotation to go from strip (u,v,w) to track coordinates");
- System.out.println("stripToDet (JLab) translation:");
- System.out.println(stripToDet.getTranslation().toString());
- System.out.println("stripToDet Rotation:");
- System.out.println(stripToDet.getRotation().toString());
- System.out.println("detToTrack Rotation:");
- System.out.println(detToTrackMatrix.toString());
-
-
- }
-
- return (Hep3Matrix) VecOp.mult(detToTrackMatrix,stripToDetMatrix);
- }
-
- private ITransform3D GetLocalToGlobal(HelicalTrackStrip strip) {
- //Transform from sensor frame (u,v,w) to tracking frame
- RawTrackerHit rth = (RawTrackerHit) strip.rawhits().get(0);
- IDetectorElement ide = rth.getDetectorElement();
- SiSensor sensor = ide.findDescendants(SiSensor.class).get(0);
- SiSensorElectrodes electrodes = sensor.getReadoutElectrodes(ChargeCarrier.HOLE);
- return electrodes.getLocalToGlobal();
- }
-
- private Hep3Matrix getTrackToStripRotation(HelicalTrackStrip strip) {
- //This function transforms the hit to the sensor coordinates
-
- //Transform from JLab frame to sensor frame (done through the RawTrackerHit)
- ITransform3D detToStrip = GetGlobalToLocal(strip);
- //Get rotation matrix
- Hep3Matrix detToStripMatrix = (BasicHep3Matrix) detToStrip.getRotation().getRotationMatrix();
- //Transformation between the JLAB and tracking coordinate systems
- Hep3Matrix detToTrackMatrix = (BasicHep3Matrix) HPSTransformations.getMatrix();
-
- if (true) {
- System.out.println("Getting the rotation to go from track to strip (u,v,w)");
- System.out.println("gblToLoc translation:");
- System.out.println(detToStrip.getTranslation().toString());
- System.out.println("gblToLoc Rotation:");
- System.out.println(detToStrip.getRotation().toString());
- System.out.println("detToTrack Rotation:");
- System.out.println(detToTrackMatrix.toString());
- }
-
- return (Hep3Matrix) VecOp.mult(detToStripMatrix, VecOp.inverse(detToTrackMatrix));
- }
-
- private ITransform3D GetGlobalToLocal(HelicalTrackStrip strip) {
- //Transform from JLab frame (RawTrackerHit) to sensor frame (i.e. u,v,w)
- RawTrackerHit rth = (RawTrackerHit) strip.rawhits().get(0);
- IDetectorElement ide = rth.getDetectorElement();
- SiSensor sensor = ide.findDescendants(SiSensor.class).get(0);
- SiSensorElectrodes electrodes = sensor.getReadoutElectrodes(ChargeCarrier.HOLE);
- return electrodes.getGlobalToLocal();
- }
-
-
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanSurface.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/KalmanSurface.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,431 +0,0 @@
-package org.hps.recon.tracking.kalman;
-
-import hep.physics.vec.Hep3Vector;
-
-import java.util.List;
-
-import org.lcsim.detector.IDetectorElement;
-import org.lcsim.detector.IGeometryInfo;
-import org.lcsim.detector.ILogicalVolume;
-import org.lcsim.detector.material.IMaterial;
-import org.lcsim.detector.solids.ISolid;
-import org.lcsim.detector.solids.Tube;
-import org.lcsim.geometry.compact.Subdetector;
-import org.lcsim.geometry.subdetector.PolyconeSupport;
-import org.lcsim.geometry.subdetector.PolyconeSupport.ZPlane;
-import org.lcsim.material.Material;
-import org.lcsim.recon.tracking.trfbase.Surface;
-import org.lcsim.recon.tracking.trfcyl.BSurfCylinder;
-import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
-import org.lcsim.recon.tracking.trfcyl.ThinCylMs;
-import org.lcsim.recon.tracking.trfcyl.ThinCylMsSim;
-import org.lcsim.recon.tracking.trfxyp.SurfXYPlane;
-import org.lcsim.recon.tracking.trfzp.SurfZPlane;
-
-/**
- *
- * Holds summary information about one sensitive surface. The information
- * is extracted from the org.lcsim geometry system. This is a routine
- * that needs to know that org.lcsim uses mm for distances but TRF uses cm.
- *
- *
- *@author $Author: jeremy $
- *@version $Id: KalmanSurface.java,v 1.1 2011/06/01 17:08:04 jeremy Exp $
- *
- * Date $Date: 2011/06/01 17:08:04 $
- *
- */
-
-public class KalmanSurface{
-
- // TRF wants distances in cm, not mm.
- private double mmTocm = 0.1;
-
-// // Name of this Subdetector.
-// public String name = "";
-
- // Readout dimension:
- // 0 = inactive
- // 1 = strips
- // 2 = pixels
- public int rodim;
-
- // References to subdetector and detector element info.
-// public Subdetector sd = null;
-// public IDetectorElement de = null;
-
- public double radius = 0.;
- public double thick = 0.;
- public double zmin = 0.;
- public double zmax = 0.;
- public double zc = 0.;
- public double rmin = 0.;
- public double rmax = 0.;
-
- // Does a 1D measurement measure the first or second coordinate?
- public int ixy = -1;
- static public int ixy_Undef = -1;
- static public int ixy_x = 0;
- static public int ixy_y = 1;
- static public int ixy_phi = 0;
- static public int ixy_z = 1;
-
- private double[] resolutions;
-
- // Different inputs have different places to find materials.
- public Material material = null;
- public IMaterial imaterial = null;
-
- // Name of the material.
- public String matname = "Undefined";
-
- // Radiation length of the material;
- public double radl = 0;
-
- // Thickness in units of radiation length.
- public double thick_radl = 0.;
-
- // Density.
- public double density = 0.0;
-
- // Tolerance for comparing floating point numbers.
- private double tolerance = 1.e-4;
-
- // Types
- private int type;
- // Types
- static public int type_tube=0;
- static public int type_zdisc=1;
- static public int type_xyplane=2;
- Surface surf = null;
-
- KalmanSurface(String string, double distnorm, double normphi) {
- if(string.equals("xyplane")){
- System.out.println("dist norm = " + distnorm + ", normphi = " + normphi);
- surf = new SurfXYPlane(distnorm, normphi);
- type = type_xyplane;
- } else{
- throw new UnsupportedOperationException("Not yet implemented");
- }
- }
-
- public SurfCylinder getCylinder(){
- if ( type != type_tube ){
-// System.out.println("Cannot return a cylinder from a non-tube surface: " + name );
- System.exit(-1);
- }
- SurfCylinder s = new SurfCylinder( radius*mmTocm );
- ThinCylMs scat = new ThinCylMs ( thick_radl );
- ThinCylMsSim sim = new ThinCylMsSim( thick_radl );
- s.setInteractor(scat);
- s.setSimInteractor(sim);
- return s;
- }
-
- public BSurfCylinder getBCylinder(){
- if ( type != type_tube ){
-// System.out.println("Cannot return a bounded cylinder from a non-tube surface: " + name );
- System.exit(-1);
- }
- return new BSurfCylinder( radius*mmTocm, zmin*mmTocm, zmax*mmTocm );
- }
-
- public boolean inBounds ( double trfarg ){
- if ( type == type_tube ){
- double z = trfarg/mmTocm;
- return (z>=zmin)&&(z<=zmax);
- } else if ( type == type_zdisc ){
- double r = trfarg/mmTocm;
- return (r>=rmin)&&(r<=rmax);
- }
- return false;
- }
-
- public int getType(){
- return type;
- }
-
- public String getTypeAsString(){
- if ( type == type_tube ){
- return "Tube ";
- } else if ( type == type_zdisc ){
- return "Zdisc";
- }
- return "Unknown";
- }
-
- public SurfZPlane getZPlane(){
- if ( type != type_zdisc ){
-// System.out.println("Cannot return a ZPlane from a non-zplane surface: " + name );
- System.exit(-1);
- }
- return new SurfZPlane( zc*mmTocm );
- }
-
- public double[] getResolutions(){
- return resolutions;
- }
-
- public double getResolution0(){
- boolean ok = true;
- if ( resolutions == null ){
- ok = false;
- } else {
- if ( resolutions.length < 1 ){
- ok = false;
- }
- }
- if ( !ok ){
-// System.out.println ("KalmanSurface: Illegal request for getResolution0 " + name );
- return 0.;
- }
- return resolutions[0];
- }
-
- public double getResolution1(){
- boolean ok = true;
- if ( resolutions == null ){
- ok = false;
- } else {
- if ( resolutions.length < 1 ){
- ok = false;
- }
- }
- if ( !ok ){
-// System.out.println ("KalmanSurface: Illegal request for getResolution1 " + name );
- return 0.;
- }
- return resolutions[1];
- }
-
-
- public KalmanSurface( IDetectorElement de, int rodim, int ixy, double[] res ){
- if ( de == null ){
-// System.out.println("KalmanSurface: Cannot instantiate RKSurf with null DetectorElement." );
- System.exit(-1);
- }
-
-// this.name = de.getName();
- this.rodim = rodim;
- this.ixy = ixy;
- resolutions = new double[res.length];
- for ( int i=0; i<res.length; ++i){
- resolutions[i] = res[i]*mmTocm;
- }
- Build( de);
- }
-
- public KalmanSurface( IDetectorElement de, int rodim, int ixy, double res ){
- if ( de == null ){
- System.out.println("KalmanSurface: Cannot instantiate KalmanSurface with null DetectorElement." );
- System.exit(-1);
- }
-
-// this.name = de.getName();
- this.rodim = rodim;
- this.ixy = ixy;
- resolutions = new double[1];
- resolutions[0] = res*mmTocm;
- Build( de);
- }
-
- public KalmanSurface( IDetectorElement de, int rodim, int ixy, double res0, double res1 ){
- if ( de == null ){
- System.out.println("KalmanSurface: Cannot instantiate KalmanSurface with null DetectorElement." );
- System.exit(-1);
- }
-
-// this.name = de.getName();
- this.rodim = rodim;
- this.ixy = ixy;
- resolutions = new double[2];
- resolutions[0] = res0*mmTocm;
- resolutions[1] = res1*mmTocm;
- Build( de);
- }
-
-
-
- public KalmanSurface( Subdetector sd ){
- if ( sd == null ){
- System.out.println("KalmanSurface: Cannot instantiate RKSurf with null Subdetector." );
- System.exit(-1);
- }
-
-// this.sd = sd;
- rodim = 0;
-// name = sd.getName();
- Build( sd );
- }
-
- // Adjust z position. Used in the equivstrips model.
- public void hackZc( double dz){
- zc+=dz;
- }
-
- public void Print(){
- if ( type == type_tube ){
-// System.out.printf ( "%-40s %1d %3d %10.2f %10.2f %10.2f %10.2f %-20s %10.2f %10.6f %10.6f\n",
-// name,
-// rodim,
-// ixy,
-// radius,
-// thick,
-// zmin,
-// zmax,
-// matname,
-// radl,
-// thick_radl,
-// density
-// );
- } else {
-// System.out.printf ( "%-40s %1d %3d %10.2f %10.2f %10.2f %10.2f %-20s %10.2f %10.6f %10.6f\n",
-// name,
-// rodim,
-// ixy,
-// zc,
-// thick,
-// rmin,
-// rmax,
-// matname,
-// radl,
-// thick_radl,
-// density
-// );
- }
- }
-
-
- private void Build( Subdetector sd ){
-
- String classname = sd.getClass().getName();
- String shortname = classname.replaceAll("org.lcsim.geometry.subdetector.","");
-
- if ( shortname.compareTo("PolyconeSupport") == 0 ) {
- PolyconeSupport pc = (PolyconeSupport) sd;
- AddPolycone(pc);
-
- } else {
-// System.out.println("KalmanSurface: Do not know how to add this subdetector: "
-// + name + " " + classname);
- System.exit(-1);
- }
-
- matname = material.getName();
- radl = material.getRadiationLength();
-
- thick_radl = thick/radl;
-
- }
-
- private void Build( IDetectorElement de ){
-
- //System.out.println ("Starting: " + name );
- IGeometryInfo g = de.getGeometry();
- if ( g == null ){
-// System.out.println("Missing geometry for detector element: " + name );
- System.exit(-1);
- }
- ILogicalVolume lv = g.getLogicalVolume();
- if ( lv == null ){
-// System.out.println("Missing logical volume for detector element: " + name );
- System.exit(-1);
- }
- IMaterial mat = lv.getMaterial();
- ISolid solid = lv.getSolid();
- Hep3Vector center = g.getPosition();
-
- String solidname = getSolidTypeName(solid);
-
- if ( solidname.compareTo("Tube") == 0 ){
- Tube tube = (Tube) solid;
- if ( tube.getZHalfLength() > 1.0 ){
- AddTube( tube, center );
- }else{
- AddZDisc( tube, center);
- }
-
- } else {
-// System.out.println ( "KalmanSurface: Do not recognize this shape for this DetectorElement: "
-// + solidname + " " + name );
- System.exit(-1);
- }
- matname = mat.getName();
- radl = mat.getRadiationLength();
- density = mat.getDensity();
-
- thick_radl = thick/radl;
-
- }
-
- private void AddPolycone ( PolyconeSupport pc ){
- type = type_tube;
- List<ZPlane> zplanes = pc.getZPlanes();
- boolean ok = false;
- for ( int i=1; i<zplanes.size(); ++i ){
- ZPlane zp1 = zplanes.get(i-1);
- ZPlane zp2 = zplanes.get(i);
- if ( zp1.getZ() < 0. && zp2.getZ() > 0. ){
- if ( Math.abs(zp1.getRMax()-zp1.getRMax()) > tolerance ||
- Math.abs(zp1.getRMin()-zp1.getRMin()) > tolerance ){
-// System.out.println ("RKSurf: Mismatched radii in beampipe: "
-// + name );
- }
- radius = 0.5*(zp1.getRMax() + zp1.getRMin());
- thick = zp1.getRMax() - zp1.getRMin();
- zmin = zp1.getZ();
- zmax = zp2.getZ();
- ok = true;
- }
- }
-
- if ( !ok ){
-// System.out.println( "RKSurf: Could not parse PolyconeSupport for subdetector"
-// + name );
- System.exit(-1);
- }
-
- // Depracated but replacement appears not to be in place.
- material = pc.getMaterial();
-
- radl = material.getRadiationLength();
- density = material.getDensity();
-
- }
-
- private void AddTube ( Tube tube, Hep3Vector center ){
- type = type_tube;
- if ( Math.abs( center.x()) > tolerance ||
- Math.abs( center.y()) > tolerance ||
- Math.abs( center.z()) > tolerance ){
- System.out.println("RKSurf: Only do tubes centered on (0.0,0.0,0.) " );
- System.exit(-1);
- }
-
- zmin = -tube.getZHalfLength();
- zmax = tube.getZHalfLength();
- radius = 0.5*(tube.getOuterRadius()+tube.getInnerRadius());
- thick = tube.getOuterRadius()-tube.getInnerRadius();
-
- }
-
- private void AddZDisc ( Tube tube, Hep3Vector center ){
- type = type_zdisc;
- zmin = center.z()-tube.getZHalfLength();
- zmax = center.z()+tube.getZHalfLength();
- radius = 0.5*(tube.getOuterRadius()+tube.getInnerRadius());
- thick = 2.*tube.getZHalfLength();
-
- zc = center.z();
- rmin = tube.getInnerRadius();
- rmax = tube.getOuterRadius();
- }
-
- // Utility function to extract a short version of the class name.
- private String getSolidTypeName( ISolid s ){
- String fullname = s.getClass().getName();
- String name = fullname.replaceAll( "org.lcsim.detector.solids.", "");
- return name;
- }
-
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PlayingWithTracksDriver.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PlayingWithTracksDriver.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,64 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-
-package org.hps.recon.tracking.kalman;
-
-import java.util.List;
-
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.MCParticle;
-import org.lcsim.event.Track;
-import org.lcsim.event.TrackerHit;
-import org.lcsim.fit.helicaltrack.HelicalTrack2DHit;
-import org.lcsim.fit.helicaltrack.HelicalTrack3DHit;
-import org.lcsim.fit.helicaltrack.HelicalTrackCross;
-import org.lcsim.fit.helicaltrack.HelicalTrackHit;
-import org.lcsim.util.Driver;
-
-/**
- *
- * @author ecfine
- */
-public class PlayingWithTracksDriver extends Driver{
-
- public PlayingWithTracksDriver() {
- }
-
-
- public void process(EventHeader event){
-
- System.out.println("event number " + event.getEventNumber());
- System.out.println("event detector name: " +event.getDetectorName());
- List<MCParticle> mcparticles = event.getMCParticles();
- System.out.println("number of mcparticles: " + mcparticles.size());
- MCParticle firstParticle = mcparticles.get(0);
- System.out.println("energy of first particle: " + firstParticle.getEnergy());
- if (event.hasItem("MatchedTracks")) {
- List<Track> trklist = (List<Track>) event.get("MatchedTracks");
- System.out.println("number of tracks: " + trklist.size());
- for (int i = 0; i < trklist.size(); i ++){
- Track trk = trklist.get(i);
- for (int k = 0; k < trk.getTrackerHits().size(); k++){
- TrackerHit hit = (TrackerHit) trk.getTrackerHits().get(k);
- if (hit instanceof HelicalTrack2DHit){
- System.out.println("HelicalTrack2DHit");
- } else
- if (hit instanceof HelicalTrack3DHit){
- System.out.println("HelicalTrack3DHit");
- } else if(hit instanceof HelicalTrackCross){
- System.out.println("HelicalTrackCross");
- } else
- if (hit instanceof HelicalTrackHit){
- System.out.println("HelicalTrackHit");
- } else {
- System.out.println("Don't know this kind of track...");
- }
- }
- }
- } else {
- System.out.println("No tracks!");
- }
- }
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PrintDetectorElements.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PrintDetectorElements.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,134 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-
-package org.hps.recon.tracking.kalman;
-
-import java.util.ArrayList;
-import java.util.List;
-
-import org.lcsim.detector.ILogicalVolume;
-import org.lcsim.detector.IPhysicalVolume;
-import org.lcsim.detector.ITransform3D;
-import org.lcsim.detector.solids.Point3D;
-import org.lcsim.geometry.Detector;
-
-/**
- *
- * @author ecfine
- */
-
-public class PrintDetectorElements{
- ShapeDispatcher shapeDispatcher = new ShapeDispatcher();
- static int indent = 0;
- List<Point3D> coords = null;
- static ArrayList physicalVolumes = new ArrayList();
- Boolean printing = true;
-
-
- public void run(Detector detector, Boolean print) {
- printing = print;
- if (printing) {
- System.out.println();
- System.out.println();
- System.out.println("Detector name is " + detector.getDetectorName());
- }
- ILogicalVolume logical = getLogicalVolume(detector);
- loopThroughDaughters(logical);
- }
-
- private ILogicalVolume getLogicalVolume(Detector det) {
- ILogicalVolume logical = det.getTrackingVolume().getLogicalVolume();
- return logical;
- }
-
- private void loopThroughDaughters(ILogicalVolume logical) {
- if (printing) {
- printIndent();
- System.out.println("Logical Volume: " + logical.getName());
- indent ++;
- printIndent();
- System.out.println("Solid: " + logical.getSolid().getName());
- printIndent();
- System.out.println("Material: " + logical.getMaterial().getName());
- printIndent();
- System.out.println("Volume: " + logical.getSolid().getCubicVolume());
- printIndent();
- shapeDispatcher.printShape(logical.getSolid());
- printIndent();
- shapeDispatcher.printLocalCoords(logical.getSolid());
- shapeDispatcher.printGlobalCoords(logical.getSolid());
- System.out.println();
- indent ++;
- }
- for (int n = 0; n < logical.getNumberOfDaughters(); n++) {
- IPhysicalVolume physical = logical.getDaughter(n);
- loopThroughDaughters(physicalToLogical(physical));
- indent = (indent - 2);
- physicalVolumes.remove(physicalVolumes.size() - 1);
- }
- }
-
-
- // Prints the information in the physical volume and returns the logical
- // volume associated with the physical volume. Adds the physical volume to
- // the physicalVolumes ArrayList. Obtains the local rotation and caluculates
- // global rotation.
- private ILogicalVolume physicalToLogical(IPhysicalVolume physical) {
- physicalVolumes.add(physical.getTransform());
- if (printing){
- printIndent();
- System.out.println("Physical Volume: " + physical.getName());
-
- indent ++;
- printIndent();
- System.out.println("Sensitive? " + physical.isSensitive());
- printIndent();
- System.out.println("Local Translation: " + physical.getTranslation());
- printIndent();
- System.out.println("Local Rotation: [" + physical.getRotation().getComponent(0,0) +
- " " + physical.getRotation().getComponent(0, 1) +
- " " + physical.getRotation().getComponent(0, 2));
- printIndent();
- System.out.println(" " + physical.getRotation().getComponent(1,0) +
- " " + physical.getRotation().getComponent(1, 1) +
- " " + physical.getRotation().getComponent(1, 2));
- printIndent();
- System.out.println(" " + physical.getRotation().getComponent(2,0) +
- " " + physical.getRotation().getComponent(2, 1) +
- " " + physical.getRotation().getComponent(2, 2) + "]");
- }
- getGlobalTransform(physical);
- indent --;
- return physical.getLogicalVolume();
- }
-
- private void getGlobalTransform(IPhysicalVolume physical) {
- ITransform3D transform = (ITransform3D) physicalVolumes.get(0);
- for (int i = 1; i < physicalVolumes.size(); i++){
- ITransform3D lastTransform = (ITransform3D) physicalVolumes.get(i);
- transform.multiplyBy(lastTransform);
- }
- printIndent();
- System.out.println("Global Transform: " + transform.getTranslation());
- printIndent();
- System.out.println(" [" + transform.getRotation().getComponent(0,0) +
- " " + transform.getRotation().getComponent(0, 1) +
- " " + transform.getRotation().getComponent(0, 2));
- printIndent();
- System.out.println(" " + transform.getRotation().getComponent(1,0) +
- " " + transform.getRotation().getComponent(1, 1) +
- " " + transform.getRotation().getComponent(1, 2));
- printIndent();
- System.out.println(" " + transform.getRotation().getComponent(2,0) +
- " " + transform.getRotation().getComponent(2, 1) +
- " " + transform.getRotation().getComponent(2, 2) + "]");
- }
-
- public void printIndent () {
- for (int k = indent; k > 0; k --) {
- System.out.print(" ");
- }
- }
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PrintDetectorElementsDriver.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PrintDetectorElementsDriver.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,24 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-
-package org.hps.recon.tracking.kalman;
-
-import org.lcsim.geometry.Detector;
-import org.lcsim.util.Driver;
-
-
-/**
- *
- * @author ecfine
- */
-public class PrintDetectorElementsDriver extends Driver{
- final boolean printing = true; // Determines whether results are printed
-
-
- public void detectorChanged(Detector detector) {
- PrintDetectorElements loop = new PrintDetectorElements();
- loop.run(detector, printing);
- }
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PrintIntersections.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PrintIntersections.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,130 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-
-package org.hps.recon.tracking.kalman;
-
-import java.util.ArrayList;
-
-import org.lcsim.detector.ILogicalVolume;
-import org.lcsim.detector.IPhysicalVolume;
-import org.lcsim.detector.ITransform3D;
-import org.lcsim.event.Track;
-import org.lcsim.geometry.Detector;
-
-/**
- *
- * @author ecfine
- */
-
-public class PrintIntersections{
- ShapeDispatcher shapeDispatcher = new ShapeDispatcher();
- static int indent = 0;
- static ArrayList physicalVolumes = new ArrayList();
- Track track;
-
- public void run(Detector detector, Track trk) {
- System.out.println();
- System.out.println();
- System.out.println("Detector name is " + detector.getDetectorName());
- track = trk;
- ILogicalVolume logical = getLogicalVolume(detector);
- loopThroughDaughters(logical);
-
- }
-
- private ILogicalVolume getLogicalVolume(Detector det) {
- ILogicalVolume logical = det.getTrackingVolume().getLogicalVolume();
- return logical;
- }
-
- private void loopThroughDaughters(ILogicalVolume logical) {
- printIndent();
- System.out.println("Logical Volume: " + logical.getName());
- indent ++;
-// printIndent();
-// System.out.println("Solid: " + logical.getSolid().getName());
-// printIndent();
-// System.out.println("Material: " + logical.getMaterial().getName());
-// printIndent();
-// System.out.println("Volume: " + logical.getSolid().getCubicVolume());
-// printIndent();
-// shapeDispatcher.printShape(logical.getSolid());
-// printIndent();
-// shapeDispatcher.printLocalCoords(logical.getSolid());
-// shapeDispatcher.printGlobalCoords(logical.getSolid());
-// System.out.println();
-// indent ++;
- if (logical.getDaughters().size() == 0){
- System.out.println("finding intersections");
- shapeDispatcher.findIntersection(logical.getSolid(), track);
- }
- for (int n = 0; n < logical.getNumberOfDaughters(); n++) {
- IPhysicalVolume physical = logical.getDaughter(n);
- loopThroughDaughters(physicalToLogical(physical));
- indent = (indent - 1);
- physicalVolumes.remove(physicalVolumes.size() - 1);
- }
- }
-
-
- // Prints the information in the physical volume and returns the logical
- // volume associated with the physical volume. Adds the physical volume to
- // the physicalVolumes ArrayList. Obtains the local rotation and caluculates
- // global rotation.
- private ILogicalVolume physicalToLogical(IPhysicalVolume physical) {
- physicalVolumes.add(physical.getTransform());
- printIndent();
- System.out.println("Physical Volume: " + physical.getName());
-// indent ++;
-// printIndent();
-// System.out.println("Sensitive? " + physical.isSensitive());
-// printIndent();
-// System.out.println("Local Translation: " + physical.getTranslation());
-// printIndent();
-// System.out.println("Local Rotation: [" + physical.getRotation().getComponent(0,0) +
-// " " + physical.getRotation().getComponent(0, 1) +
-// " " + physical.getRotation().getComponent(0, 2));
-// printIndent();
-// System.out.println(" " + physical.getRotation().getComponent(1,0) +
-// " " + physical.getRotation().getComponent(1, 1) +
-// " " + physical.getRotation().getComponent(1, 2));
-// printIndent();
-// System.out.println(" " + physical.getRotation().getComponent(2,0) +
-// " " + physical.getRotation().getComponent(2, 1) +
-// " " + physical.getRotation().getComponent(2, 2) + "]");
-
-// getGlobalTransform(physical);
-// indent --;
- return physical.getLogicalVolume();
- }
-
- private void getGlobalTransform(IPhysicalVolume physical) {
- ITransform3D transform = (ITransform3D) physicalVolumes.get(0);
- for (int i = 1; i < physicalVolumes.size(); i++){
- ITransform3D lastTransform = (ITransform3D) physicalVolumes.get(i);
- transform.multiplyBy(lastTransform);
- }
- printIndent();
- System.out.println("Global Transform: " + transform.getTranslation());
- printIndent();
- System.out.println(" [" + transform.getRotation().getComponent(0,0) +
- " " + transform.getRotation().getComponent(0, 1) +
- " " + transform.getRotation().getComponent(0, 2));
- printIndent();
- System.out.println(" " + transform.getRotation().getComponent(1,0) +
- " " + transform.getRotation().getComponent(1, 1) +
- " " + transform.getRotation().getComponent(1, 2));
- printIndent();
- System.out.println(" " + transform.getRotation().getComponent(2,0) +
- " " + transform.getRotation().getComponent(2, 1) +
- " " + transform.getRotation().getComponent(2, 2) + "]");
- }
-
- public void printIndent() {
- for (int k = indent; k > 0; k --) {
- System.out.print(" ");
- }
- }
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PropDCACyl.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PropDCACyl.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,732 +0,0 @@
-package org.hps.recon.tracking.kalman;
-
-
-import org.lcsim.recon.tracking.trfbase.PropDir;
-import org.lcsim.recon.tracking.trfbase.PropDirected;
-import org.lcsim.recon.tracking.trfbase.PropStat;
-import org.lcsim.recon.tracking.trfbase.Propagator;
-import org.lcsim.recon.tracking.trfbase.Surface;
-import org.lcsim.recon.tracking.trfbase.TrackDerivative;
-import org.lcsim.recon.tracking.trfbase.TrackVector;
-import org.lcsim.recon.tracking.trfbase.VTrack;
-import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
-import org.lcsim.recon.tracking.trfdca.SurfDCA;
-import org.lcsim.recon.tracking.trfutil.Assert;
-import org.lcsim.recon.tracking.trfutil.TRFMath;
-/**
- * Propagates tracks from a DCA surface to a Cylinder.
- *<p>
- * Propagation will fail if either the origin is not a DCA surface,
- * or the destination is not a Cylinder.
- *<p>
- * The default direction for propagation is forward.
- *<p>
- * Propagation to a cylinder at the radius of the DCA will succeed
- * and the track parameters are valid but the errors are not valid.
- *
- *
- *@author Norman A. Graf
- *@version 1.0
- *
- */
-public class PropDCACyl extends PropDirected
-{
-
- // static variables
- // Assign track parameter indices
-
-
- public static final int IRSIGNED = SurfDCA.IRSIGNED;
- public static final int IZ_DCA = SurfDCA.IZ;
- public static final int IPHID = SurfDCA.IPHID;
- public static final int ITLM_DCA = SurfDCA.ITLM;
- public static final int IQPT_DCA = SurfDCA.IQPT;
-
- public static final int IPHI = SurfCylinder.IPHI;
- public static final int IZ_CYL = SurfCylinder.IZ;
- public static final int IALF = SurfCylinder.IALF;
- public static final int ITLM_CYL = SurfCylinder.ITLM;
- public static final int IQPT_CYL = SurfCylinder.IQPT;
-
- // Attributes ***************************************************
-
-
-
- // bfield * BFAC
- private double _bfac;
-
-
- // Methods ******************************************************
-
- // static methods
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public static String typeName()
- { return "PropDCACyl";
- }
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public static String staticType()
- { return typeName();
- }
-
-
-
- /**
- *Construct an instance from a constant solenoidal magnetic field in Tesla.
- *
- * @param bfield The magnetic field strength in Tesla.
- */
- public PropDCACyl(double bfield)
- {
- super(PropDir.FORWARD);
- _bfac=bfield * TRFMath.BFAC;
- }
-
- /**
- *Clone an instance.
- *
- * @return A Clone of this instance.
- */
- public Propagator newPropagator()
- {
- return new PropDCACyl( bField() );
- }
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public String type()
- { return staticType();
- }
-
- /**
- *Propagate a track without error in the specified direction.
- *and return the derivative matrix in deriv.
- *
- * The track parameters for a cylinder are:
- * phi z alpha tan(lambda) curvature
- *
- * @param trv The VTrack to propagate.
- * @param srf The Surface to which to propagate.
- * @param dir The direction in which to propagate.
- * @param deriv The track derivatives to update at the surface srf.
- * @return The propagation status.
- */
- public PropStat vecDirProp( VTrack trv, Surface srf,
- PropDir dir, TrackDerivative deriv )
- {
- return dcaCylPropagate( _bfac, trv, srf, dir, deriv );
- }
-
- /**
- *Propagate a track without error in the specified direction.
- *
- * The track parameters for a cylinder are:
- * phi z alpha tan(lambda) curvature
- *
- * @param trv The VTrack to propagate.
- * @param srf The Surface to which to propagate.
- * @param dir The direction in which to propagate.
- * @return The propagation status.
- */
- public PropStat vecDirProp( VTrack trv, Surface srf,
- PropDir dir )
- {
- TrackDerivative deriv =null;
- return vecDirProp(trv, srf, dir, deriv);
-
- }
-
- // propagate a track with error in the specified direction
- // PropStat err_dir_prop( ETrack& tre, Surface& srf,
- // PropDir dir ) ;
-
- //
-
- /**
- *Return the strength of the magnetic field in Tesla.
- *
- * @return The strength of the magnetic field in Tesla.
- */
- public double bField()
- {
- return _bfac/TRFMath.BFAC;
- }
-
-
- /**
- *output stream
- * @return A String representation of this instance.
- */
- public String toString()
- {
- return "DCA propagation to a Cylinder with constant "
- + bField() + " Tesla field";
- }
-
-
-
- /**
- * Propagate from dca to cylinder.
- * Why is this public?
- *
- * @param _bfac The numerical factor (including the field)
- * @param trv The VTrack to propagate.
- * @param srf The cylindrical surface to which to propagate.
- * @param dir The direction in which to propagate.
- * @param deriv The track derivatives to update at the surface srf.
- * @return The propagation status.
- */
- public PropStat dcaCylPropagate( double _bfac,
- VTrack trv,
- Surface srf,
- PropDir dir,
- TrackDerivative deriv )
- {
-
- // construct return status
- PropStat pstat = new PropStat();
-
- // fetch the originating Surface and check it is a DCA surface
- Surface srf0 = trv.surface();
- Assert.assertTrue( srf0.pureType().equals(SurfDCA.staticType() ));
- if ( ! srf0.pureType( ).equals(SurfDCA.staticType()) ){
- System.out.println("Failing PropDCACyl because it's not a DCA surface");
- return pstat;
- }
- SurfDCA srf_dca = ( SurfDCA ) srf0;
-
- // Check that dca surface has zero tilt.
-
- boolean tilted = srf_dca.dXdZ() != 0 || srf_dca.dYdZ() != 0;
- Assert.assertTrue(!tilted);
- if(tilted) {
- System.out.println("Failing PropDCACyl DCA surface is tilted");
- return pstat;
- }
-
- // check that the destination surface is a Cylinder
- Assert.assertTrue( srf.pureType().equals(SurfCylinder.staticType()) );
- if ( !srf.pureType( ).equals(SurfCylinder.staticType()) ){
- System.out.println("Failing PropDCACyl because it's not a Cyl surface");
- return pstat;
- }
- SurfCylinder srf_cyl = ( SurfCylinder ) srf;
- //SurfCylinder srf_cyl = new SurfCylinder( (SurfCylinder) srf);
-
- // fetch the originating TrackVector
- TrackVector vec_dca = trv.vector();
- double r1p = Math.abs(vec_dca.get(IRSIGNED)); // r
- double z1 = vec_dca.get(IZ_DCA); // z
- double phid1p= vec_dca.get(IPHID); // phi_direction
- double tlam1 = vec_dca.get(ITLM_DCA); // tan(lambda)
- double qpt1 = vec_dca.get(IQPT_DCA); // q/pT
-
- double xv = srf_dca.x();
- double yv = srf_dca.y();
-
- // calculate alf1 and phi1
- double sign_alf1p;
- double sign_r1p=0.;
-
- if ( !TRFMath.isZero( vec_dca.get(IRSIGNED) ) )
- {
- sign_alf1p = vec_dca.get(IRSIGNED)/Math.abs(vec_dca.get(IRSIGNED));
- sign_r1p = sign_alf1p;
- }
- else
- {
- sign_alf1p = 0.0;
- }
-
- if( (!TRFMath.isZero(xv) || !TRFMath.isZero(yv)) && TRFMath.isZero(sign_alf1p) ) sign_alf1p=1;
-
- double alf1p = sign_alf1p * TRFMath.PI2; // alpha
- double phi1p = phid1p - alf1p; // phi_position
- phi1p = TRFMath.fmod2( phi1p, TRFMath.TWOPI );
-
- // calculate salf1, calf1 and crv1
- double salf1p=0.;
- salf1p = alf1p>0. ? 1.:( alf1p < 0. ? -1. : 0.) ;
-
- double cosphi1p=Math.cos(phi1p);
- double sinphi1p=Math.sin(phi1p);
-
- double x1=r1p*cosphi1p+xv;
- double y1=r1p*sinphi1p+yv;
-
- double phi1 = Math.atan2(y1,x1); // phi position in (xv,yv) coord system
- double r1 = Math.sqrt(x1*x1+y1*y1); // r1 in (xv,yv) coord system
- double alf1 = alf1p + phi1p - phi1; // alf1 in (xv,yv) coord system
- if( TRFMath.isZero(x1) && TRFMath.isZero(y1) )
- {
- phi1 = phid1p;
- alf1 = 0.;
- }
- if( sign_r1p == 0 && !TRFMath.isZero(r1) )
- sign_r1p=1;
-
- // calculate salf1, calf1 and crv1
- double salf1 = Math.sin( alf1 );
- double calf1 = Math.cos( alf1 );
-
- if( TRFMath.isZero(alf1) )
- {
- salf1=0.;
- calf1=1.;
- }
- if( TRFMath.isEqual(Math.abs(alf1),TRFMath.PI2) )
- {
- salf1= alf1 > 0 ? 1. : -1. ;
- calf1=0.;
- }
-
- double crv1 = _bfac * qpt1;
- double sign_crv = 1.;
-
- // fetch r2 of the destination Cylinder
- double r2 = srf_cyl.parameter(SurfCylinder.RADIUS);
-
- // calculate tlam2, qpt2 and crv2 of the destination Cylinder
- double tlam2 = tlam1;
- double qpt2 = qpt1;
- double crv2 = crv1;
-
- // calculate salf2 of the destination Cylinder
- double salf2 = r1/r2*salf1 + 0.5*crv1/r2*(r2*r2-r1*r1);
-
- // If salf2 is close to 1 or -1, set it to that value.
- double diff = Math.abs( Math.abs(salf2) - 1.0 );
- if ( diff < 1.e-10 ) salf2 = salf2>0 ? 1.0 : -1.0;
- // if salf2 > 1, track does not cross cylinder
- if ( Math.abs(salf2) > 1.0 ) {
- System.out.println("Failing PropDCACyl track does not cross cylinder: salf2 = "+salf2);
- System.out.println("r1 = "+r1+"; r2 = "+r2+"; salf1="+salf1+"; crv1 = "+crv1);
- System.out.println("r1p = "+r1p+"; x1 = "+x1+"; y1 ="+y1+"; phi1p = "+phi1p);
- return pstat;
- }
-
- // there are two possibilities for alf2
- double alf21 = Math.asin( salf2 );
- double alf22 = alf21>0 ? Math.PI-alf21 : -Math.PI-alf21;
- double calf21 = Math.cos( alf21 );
- double calf22 = Math.cos( alf22 );
- //double phi21 = phi1 + atan2( -sign_crv*calf21, r2*crv2-salf2 );
- //double phi22 = phi1 + atan2( -sign_crv*calf22, r2*crv2-salf2 );
- // double cnst = salf1-r1*crv1 > 0 ? TRFMath.PI2 : -TRFMath.PI2;
- // cnst = (r1 == 0.) ? 0. : cnst;
-
- double cnst = Math.atan2(salf1-r1*crv1,calf1);
- if( TRFMath.isEqual(Math.abs(alf1),TRFMath.PI2) )
- {
- cnst = salf1-r1*crv1 > 0 ? TRFMath.PI2 : -TRFMath.PI2;
- cnst = (r1==0.) ? 0. : cnst;
- }
-
- double phi21 = phi1 + cnst - Math.atan2( salf2-r2*crv2, calf21 );
- double phi22 = phi1 + cnst - Math.atan2( salf2-r2*crv2, calf22 );
-
- if( TRFMath.isZero(crv1) )
- {
- phi21 = phi1 + cnst - alf21;
- phi22 = phi1 + cnst - alf22;
- }
-
- if ( TRFMath.isZero(calf21) )
- {
- phi21 = phi1;
- phi22 = phi1;
- }
-
- // construct an sT object for each solution
- ST_DCACyl sto1 = new ST_DCACyl(r1,phi1,alf1,crv1,r2,phi21,alf21);
- ST_DCACyl sto2 = new ST_DCACyl(r1,phi1,alf1,crv1,r2,phi22,alf22);
- // check the two solutions are nonzero and have opposite sign
- // or at least one is nonzero
-
- // choose the correct solution
- boolean use_first_solution;
-
- if( dir.equals(PropDir.NEAREST))
- use_first_solution = Math.abs(sto2.st()) > Math.abs(sto1.st());
-
- else if( dir.equals(PropDir.FORWARD))
- use_first_solution = sto1.st() > 0.0;
-
- else if( dir.equals(PropDir.BACKWARD))
- use_first_solution = sto1.st() < 0.0;
-
- else
- {
- use_first_solution = false;
- System.out.println( "PropCyl._vec_propagate: Unknown direction.");
- System.exit(1);
- }
-
- // assign phi2, alf2 and sto2 for the chosen solution
- double phi2, alf2;
- ST_DCACyl sto = new ST_DCACyl();
- double calf2;
- if ( use_first_solution )
- {
- sto = sto1;
- phi2 = phi21;
- alf2 = alf21;
- calf2 = calf21;
- }
- else
- {
- sto = sto2;
- phi2 = phi22;
- alf2 = alf22;
- calf2 = calf22;
- }
-
- // check alpha range
- Assert.assertTrue( Math.abs(alf2) <= Math.PI );
-
- // fetch sT
- double st = sto.st();
- double s = st*Math.sqrt(1+tlam1*tlam1);
-
- // calculate z2 of the destination Cylinder
- double z2 = z1 + st * tlam1;
-
- // construct the destination TrackVector
- TrackVector vec_cyl = new TrackVector();
- /*
-#ifdef TRF_PHIRANGE
- vec_cyl(IPHI) = fmod1(phi2, TWOPI);
-#else
- vec_cyl(IPHI) = phi2;
-#endif
- */
- vec_cyl.set(IPHI , phi2);
- vec_cyl.set(IZ_CYL , z2);
- vec_cyl.set(IALF , alf2);
- vec_cyl.set(ITLM_CYL , tlam2);
- vec_cyl.set(IQPT_CYL , qpt2);
-/*
- // For axial tracks, zero z and tan(lambda).
-
- if(trv.is_axial()) {
- vec_cyl(SurfDCA::IZ) = 0.;
- vec_cyl(SurfDCA::ITLM) = 0.;
- }
- */
-
- // set the surface of trv to the destination Cylinder
- trv.setSurface( srf_cyl.newPureSurface() );
-
- // set the vector of trv to the destination TrackVector (Cyl. coord.)
- trv.setVector(vec_cyl);
-
- // set the direction of trv
- if ( Math.abs(alf2) <= TRFMath.PI2 ) trv.setForward();
- else trv.setBackward();
-
- // set the return status
- pstat.setPathDistance(s);
-
- // exit now if user did not ask for error matrix.
- if ( deriv == null ) {
- System.out.println("PropDCACyl did not ask for error matrix");
- return pstat;
- }
-
-
- // calculate derivatives
-
- //double dalf2_dr1 = (salf1-r1*crv1)/(r2*calf2);
- //double dalf2_dr1 = (1.0-r1*crv1*salf1)/(r2*calf2);
- // commented out by SK ( DCA to DCA(xv,yv) change)
- double salf12 = sign_r1p*salf1;
- if( sign_r1p == 0 && salf1 == 0. )
- salf12 = 1;
- double dalf2_dr1_sign = (salf12-crv1*r1*sign_r1p)/r2/calf2;
- double dalf2_dr1 = (salf1-crv1*r1)/r2/calf2;
- double dalf2_dalf1_or = 1/r2*calf1/calf2; // over r1
- double dalf2_dalf1 = r1*dalf2_dalf1_or;
- double dalf2_dcrv1 = (r2*r2-r1*r1)/(2.0*r2*calf2);
-
- //double dphi2_dr1 = -sign_crv*
- // (r2*crv2*salf2-1.0)*(r1*crv1-salf1)/
- // (r2*calf2*(1.0 + r2*r2*crv2*crv2 - 2.0*r2*crv2*salf2));
- // commented out on 10.16.2001 by SK ( DCA to DCA(xv,yv) change)
- // double dphi2_dr1 = -sign_crv *
- // (r2*crv2*salf2-1.0)*(sign_alf1*r1*crv1-1.)/
- // (r2*calf2*(1.0 + r2*r2*crv2*crv2 - 2.0*r2*crv2*salf2));
- //
- // double dphi2_dcrv1 = -sign_crv*
- // ((1.0-r2*crv2*salf2)*(r2*r2-r1*r1)-2.0*r2*r2*calf2*calf2)/
- // (2.0*r2*calf2*(1.0 + r2*r2*crv2*crv2 - 2.0*r2*crv2*salf2));
-
- double dphi2_dphi1=1.;
- double dphi2_dr1_sign = -crv1*calf1*sign_r1p/(1.-2.*r1*crv1*salf1+r1*r1*crv1*crv1)
- - dalf2_dr1_sign*(1-salf2*crv2*r2)/(1.-2.*r2*crv2*salf2+r2*r2*crv2*crv2);
-
- double dphi2_dr1 = -crv1*calf1/(1.-2.*r1*crv1*salf1+r1*r1*crv1*crv1)
- - dalf2_dr1*(1-salf2*crv2*r2)/(1.-2.*r2*crv2*salf2+r2*r2*crv2*crv2);
- double dphi2_dalf1_m1_or = (-salf1*crv1-r1*crv1*crv1+2.*crv1*salf1)/(1.-2.*r1*crv1*salf1+r1*r1*crv1*crv1)
- - dalf2_dalf1_or*(1-salf2*crv2*r2)/(1.-2.*r2*crv2*salf2+r2*r2*crv2*crv2); // minus 1
- double dphi2_dalf1 = dphi2_dalf1_m1_or*r1 + 1;
- double dphi2_dcrv1 = -r1*calf1/(1.-2.*r1*crv1*salf1+r1*r1*crv1*crv1)
- -(dalf2_dcrv1*(1.-r2*crv2*salf2) - r2*calf2)/(1.-2.*r2*crv2*salf2+r2*r2*crv2*crv2);
-
-
- //double dst_dr1 = sto.d_st_dr1(dphi2_dr1, dalf2_dr1 );
- // double dst_dr1 = sto.d_st_dr1(sign_alf1*dphi2_dr1,sign_alf1*dalf2_dr1 );
- // double dst_crv1 = sto.d_st_dcrv1(dphi2_dcrv1, dalf2_dcrv1);
- double dst_dr1_sign = sto.d_st_dr1_sign(sign_r1p,dphi2_dr1_sign,dalf2_dr1_sign);
- double dst_dr1 = sto.d_st_dr1(dphi2_dr1,dalf2_dr1);
- double dst_dcrv1 = sto.d_st_dcrv1(dphi2_dcrv1, dalf2_dcrv1);
- double dst_dalf1_or = sto.d_st_dalf1_or(r1,dphi2_dalf1_m1_or, dalf2_dalf1_or);
-
-
- // build the derivative matrix
- // derivatives to prime coordinates
-
- double dphi1_dphi1p_tr = r1p*Math.cos(phi1p-phi1); // times r1
- double dphi1_dr1p_tr = sign_r1p*Math.sin(phi1p-phi1); // times r1
-
- double dalf1_dphi1p_tr = r1 - dphi1_dphi1p_tr;
- double dalf1_dr1p_tr = -dphi1_dr1p_tr;
-
- double dr1_dr1p_sign = Math.cos(phi1p-phi1);
- double dr1_dphi1p = r1p*Math.sin(phi1-phi1p);
-
-
- double dphi2_dphi1p = dphi2_dr1*dr1_dphi1p +dphi2_dalf1 - dphi1_dphi1p_tr*dphi2_dalf1_m1_or;
- double dphi2_dr1p = dphi2_dr1_sign*dr1_dr1p_sign - dphi2_dalf1_m1_or*dphi1_dr1p_tr;
-
- double dalf2_dphi1p = dalf2_dr1*dr1_dphi1p + dalf2_dalf1_or*dalf1_dphi1p_tr;
- double dalf2_dr1p = dalf2_dr1_sign*dr1_dr1p_sign + dalf2_dalf1_or*dalf1_dr1p_tr;
-
- double dst_dphi1p = dst_dr1*dr1_dphi1p + dst_dalf1_or*dalf1_dphi1p_tr;
- double dst_dr1p = dst_dr1_sign*dr1_dr1p_sign + dst_dalf1_or*dalf1_dr1p_tr;
-
- //deriv(IPHI, IRSIGNED, sign_alf1 * dphi2_dr1);
- deriv.set(IPHI, IRSIGNED, dphi2_dr1p);
- deriv.set(IPHI, IZ_DCA, 0.0);
- deriv.set(IPHI, IPHID, dphi2_dphi1p);
- deriv.set(IPHI, ITLM_DCA, 0.0);
- deriv.set(IPHI, IQPT_DCA, _bfac * dphi2_dcrv1);
-
- deriv.set(IZ_CYL, IRSIGNED, tlam1 * dst_dr1p);
- deriv.set(IZ_CYL, IZ_DCA, 1.0);
- deriv.set(IZ_CYL, IPHID, tlam1 * dst_dphi1p);
- deriv.set(IZ_CYL, ITLM_DCA, st);
- deriv.set(IZ_CYL, IQPT_DCA, tlam1 * _bfac * dst_dcrv1);
-
- //deriv.set(IALF, IRSIGNED, sign_alf1 * dalf2_dr1);
- deriv.set(IALF, IRSIGNED, dalf2_dr1p);
- deriv.set(IALF, IZ_DCA, 0.0);
- deriv.set(IALF, IPHID, dalf2_dphi1p);
- deriv.set(IALF, ITLM_DCA, 0.0);
- deriv.set(IALF, IQPT_DCA, _bfac * dalf2_dcrv1);
-
- deriv.set(ITLM_CYL, IRSIGNED, 0.0);
- deriv.set(ITLM_CYL, IZ_DCA, 0.0);
- deriv.set(ITLM_CYL, IPHID, 0.0);
- deriv.set(ITLM_CYL, ITLM_DCA, 1.0);
- deriv.set(ITLM_CYL, IQPT_DCA, 0.0);
-
- deriv.set(IQPT_CYL, IRSIGNED, 0.0);
- deriv.set(IQPT_CYL, IZ_DCA, 0.0);
- deriv.set(IQPT_CYL, IPHID, 0.0);
- deriv.set(IQPT_CYL, ITLM_DCA, 0.0);
- deriv.set(IQPT_CYL, IQPT_DCA, 1.0);
-
-
-/*
- // For axial tracks, zero all derivatives of or with respect to z or
- // tan(lambda), that are not already zero. This will force the error
- // matrix to have zero errors for z and tan(lambda).
-
- if(trv.is_axial()) {
- deriv.set(IZ_CYL, IRSIGNED, 0.);
- deriv.set(IZ_CYL, IZ_DCA, 0.);
- deriv.set(IZ_CYL, IPHID, 0.);
- deriv.set(IZ_CYL, ITLM_DCA, 0.);
- deriv.set(IZ_CYL, IQPT_DCA, 0.);
-
- deriv.set(ITLM_CYL, ITLM_DCA, 0.);
- }
- */
- return pstat;
-
- }
-
-
- // Private class STCalc.
- //
- // An STCalc_ object calculates sT (the signed transverse path length)
- // and its derivatives w.r.t. alf1 and crv1. It is constructed from
- // the starting (r1, phi1, alf1, crv1) and final track parameters
- // (r2, phi2, alf2) assuming these are consistent. Methods are
- // provided to retrieve sT and the two derivatives.
-
- class ST_DCACyl
- {
-
- private boolean _big_crv;
- private double _st;
- // double _dst_dalf2;
- private double _dst_dr1;
- private double _dst_dcrv1;
- private double _dst_dphi2;
- double _dst_dphi2_or;
- private double _crv1; // was public? need to look into this
- // constructor
- public ST_DCACyl()
- {
- }
- public ST_DCACyl(double r1, double phi1, double alf1, double crv1,
- double r2, double phi2, double alf2)
- {
-
- _crv1 = crv1;
- Assert.assertTrue( r1 >= 0.0 );
- Assert.assertTrue( r2 >= 0.0 );
- double rmax = r1+r2;
-
- // Calculate the change in xy direction
- double phi_dir_diff = TRFMath.fmod2(phi2+alf2-phi1-alf1,TRFMath.TWOPI);
- Assert.assertTrue( Math.abs(phi_dir_diff) <= Math.PI );
-
- // Evaluate whether |C| is "big"
- _big_crv = rmax*Math.abs(crv1) > 0.001;
-
- // If the curvature is big we can use
- // sT = (phi_dir2 - phi_dir1)/crv1
- if ( _big_crv )
- {
- Assert.assertTrue( crv1 != 0.0 );
- _st = phi_dir_diff/crv1;
- }
-
- // Otherwise, we calculate the straight-line distance
- // between the points and use an approximate correction
- // for the (small) curvature.
- else
- {
-
- // evaluate the distance
- double d = Math.sqrt( r1*r1 + r2*r2 - 2.0*r1*r2*Math.cos(phi2-phi1) );
- double arg = 0.5*d*crv1;
- double arg2 = arg*arg;
- double st_minus_d = d*arg2*( 1.0/6.0 + 3.0/40.0*arg2 );
- _st = d + st_minus_d;
-
- // evaluate the sign
- // We define a metric xsign = abs( (dphid-d*C)/(d*C) ).
- // Because sT*C = dphid and d = abs(sT):
- // xsign = 0 for sT > 0
- // xsign = 2 for sT < 0
- // Numerical roundoff will smear these predictions.
- double sign = 0.0;
- if ( crv1*_st != 0. )
- {
- double xsign = Math.abs( (phi_dir_diff - _st*crv1) / (_st*crv1) );
- if ( xsign < 0.5 ) sign = 1.0;
- if ( xsign > 1.5 && xsign < 3.0 ) sign = -1.0;
- }
- // If the above is indeterminate, assume zero curvature.
- // In this case abs(alpha) decreases monotonically
- // with sT. Track passing through origin has alpha = 0 on one
- // side and alpha = +/-pi on the other. If both points are on
- // the same side, we use dr/ds > 0 for |alpha|<pi/2.
- if ( sign == 0. )
- {
- sign = 1.0;
- if ( Math.abs(alf2) > Math.abs(alf1) ) sign = -1.0;
- if ( Math.abs(alf2) == Math.abs(alf1) )
- {
- if ( Math.abs(alf2) < TRFMath.PI2 )
- {
- if ( r2 < r1 ) sign = -1.0;
- }
- else
- {
- if ( r2 > r1 ) sign = -1.0;
- }
- }
- }
-
- // Correct _st using the above sign.
- Assert.assertTrue( Math.abs(sign) == 1.0 );
- _st = sign*_st;
-
- // save derivatives
- // _dst_dalf2 = (_st/d)*2.0*(r1-r2*cos(phi2-phi1));
- //_dst_dphi2 = (_st/d)*2.0*r1*r2*sin(phi2-phi1);
- // _dst_dcrv1 = sign*d*d*arg*( 1.0/6.0 + 3.0/20.0*arg2);
- // double root = (1.0 + 0.5*arg*arg + 3.0/8.0*arg*arg*arg*arg );
- // _dst_dphi2 = sign*(r1*r2*Math.sin(phi2-phi1))*root/d;
- // _dst_dr1 = sign*(r1-r2*Math.cos(phi2-phi1))*root/d;
-
- if ( TRFMath.isZero(d) )
- {
- _dst_dcrv1 = 0.0;
- _dst_dphi2 = sign*r1;
- _dst_dphi2_or = sign;
- _dst_dr1 = 0.0;
- }
- else
- {
- _dst_dcrv1 = sign*d*d*arg*( 1.0/6.0 + 3.0/20.0*arg2);
- double root = (1.0 + 0.5*arg*arg + 3.0/8.0*arg*arg*arg*arg );
- _dst_dphi2_or = sign*(r2*Math.sin(phi2-phi1))*root/d;
- _dst_dphi2 = _dst_dphi2_or*r1;
- _dst_dr1 = sign*(r1-r2*Math.cos(phi2-phi1))*root/d;
- }
- }
-
- }
- public double st()
- {
- return _st;
- }
-
- public double d_st_dr1_sign(double sign_alf1,double dphi2_dr1_sign,double dalf2_dr1_sign)
- {
- if ( _big_crv ) return ( dphi2_dr1_sign + dalf2_dr1_sign ) / _crv1;
- else return _dst_dphi2 * dphi2_dr1_sign + _dst_dr1*sign_alf1;
- }
-
- public double d_st_dr1( double d_phi2_dr1, double d_alf2_dr1 )
- {
- if ( _big_crv ) return ( d_phi2_dr1 + d_alf2_dr1 ) / _crv1;
- else return _dst_dphi2 * d_phi2_dr1 + _dst_dr1;
- }
- public double d_st_dcrv1( double d_phi2_dcrv1, double d_alf2_dcrv1 )
- {
- if ( _big_crv ) return ( d_phi2_dcrv1 + d_alf2_dcrv1 - _st ) / _crv1;
- else return _dst_dphi2 * d_phi2_dcrv1 + _dst_dcrv1;
- }
-
- public double d_st_dalf1_or(double r1,double dphi2_dalf1_m1_or, double dalf2_dalf1_or)
- {
- if ( _big_crv ) return ( dphi2_dalf1_m1_or + dalf2_dalf1_or) / _crv1;
- else return _dst_dphi2_or * (dphi2_dalf1_m1_or*r1+1);
- }
-
- }
-
-}
-
-
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PropDCAXY.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PropDCAXY.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,296 +0,0 @@
-package org.hps.recon.tracking.kalman;
-
-// May need to remove next line when moved to proper package.
-import org.lcsim.recon.tracking.trfbase.PropDir;
-import org.lcsim.recon.tracking.trfbase.PropDirected;
-import org.lcsim.recon.tracking.trfbase.PropStat;
-import org.lcsim.recon.tracking.trfbase.Propagator;
-import org.lcsim.recon.tracking.trfbase.Surface;
-import org.lcsim.recon.tracking.trfbase.TrackDerivative;
-import org.lcsim.recon.tracking.trfbase.VTrack;
-import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
-//import org.lcsim.recon.tracking.trfdca.PropDCACyl;
-import org.lcsim.recon.tracking.trfcylplane.PropCylXY;
-import org.lcsim.recon.tracking.trfdca.SurfDCA;
-import org.lcsim.recon.tracking.trfutil.Assert;
-import org.lcsim.recon.tracking.trfutil.TRFMath;
-import org.lcsim.recon.tracking.trfxyp.SurfXYPlane;
-
-/**
- * Propagates tracks from a DCA surface to an XYPlane.
- * It does so by creating an intermediate Cylinder surface,
- * and using existing code to perform the transformation in two steps.
- *<p>
- * Propagation will fail if either the origin is not a DCA surface,
- * or the destination is not a XYPlane.
- *<p>
- * The default direction for propagation is forward.
- *<p>
- * PropDCACyl and PropCylXY do not work in all circumstances. See those
- * codes for details.
- *
- *
- *@author $Author: mgraham $
- *@version $Id: PropDCAXY.java,v 1.3 2011/11/16 18:00:03 mgraham Exp $
- *
- * Date $Date: 2011/11/16 18:00:03 $
- *
- */
-public class PropDCAXY extends PropDirected {
-
- // static variables
- // Assign track parameter indices
- public static final int IRSIGNED = SurfDCA.IRSIGNED;
- public static final int IZ_DCA = SurfDCA.IZ;
- public static final int IPHID = SurfDCA.IPHID;
- public static final int ITLM_DCA = SurfDCA.ITLM;
- public static final int IQPT_DCA = SurfDCA.IQPT;
- public static final int IPHI = SurfCylinder.IPHI;
- public static final int IZ_CYL = SurfCylinder.IZ;
- public static final int IALF = SurfCylinder.IALF;
- public static final int ITLM_CYL = SurfCylinder.ITLM;
- public static final int IQPT_CYL = SurfCylinder.IQPT;
- private static final int IV = SurfXYPlane.IV;
- private static final int IZC = SurfXYPlane.IZ;
- private static final int IDVDU = SurfXYPlane.IDVDU;
- private static final int IDZDU = SurfXYPlane.IDZDU;
- private static final int IQP_XY = SurfXYPlane.IQP;
- // Attributes ***************************************************
- // bfield * BFAC
- private double _bfac;
-
- // Methods ******************************************************
- // static methods
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public static String typeName() {
- return "PropDCAXY";
- }
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public static String staticType() {
- return typeName();
- }
-
- /**
- *Construct an instance from a constant solenoidal magnetic field in Tesla.
- *
- * @param bfield The magnetic field strength in Tesla.
- */
- public PropDCAXY(double bfield) {
- super(PropDir.FORWARD);
- _bfac = bfield * TRFMath.BFAC;
- }
-
- /**
- *Clone an instance.
- *
- * @return A Clone of this instance.
- */
- public Propagator newPropagator() {
- return new PropDCAXY(bField());
- }
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public String type() {
- return staticType();
- }
-
- /**
- *Propagate a track without error in the specified direction.
- *and return the derivative matrix in deriv.
- *
- * The track parameters for a cylinder are:
- * phi z alpha tan(lambda) curvature
- *
- * @param trv The VTrack to propagate.
- * @param srf The Surface to which to propagate.
- * @param dir The direction in which to propagate.
- * @param deriv The track derivatives to update at the surface srf.
- * @return The propagation status.
- */
- public PropStat vecDirProp(VTrack trv, Surface srf,
- PropDir dir, TrackDerivative deriv) {
- return dcaXYPropagate(_bfac, trv, srf, dir, deriv);
- }
-
- /**
- *Propagate a track without error in the specified direction.
- *
- * The track parameters for a cylinder are:
- * phi z alpha tan(lambda) curvature
- *
- * @param trv The VTrack to propagate.
- * @param srf The Surface to which to propagate.
- * @param dir The direction in which to propagate.
- * @return The propagation status.
- */
- public PropStat vecDirProp(VTrack trv, Surface srf,
- PropDir dir) {
- TrackDerivative deriv = null;
- return vecDirProp(trv, srf, dir, deriv);
-
- }
-
- // propagate a track with error in the specified direction
- // PropStat err_dir_prop( ETrack& tre, Surface& srf,
- // PropDir dir ) ;
- //
- /**
- *Return the strength of the magnetic field in Tesla.
- *
- * @return The strength of the magnetic field in Tesla.
- */
- public double bField() {
- return _bfac / TRFMath.BFAC;
- }
-
- /**
- *output stream
- * @return A String representation of this instance.
- */
- public String toString() {
- return "DCA propagation to an XYPlane with constant "
- + bField() + " Tesla field";
- }
-
- /**
- * Propagate from dca to cylinder.
- * Why is this public?
- *
- * @param _bfac The numerical factor (including the field)
- * @param trv The VTrack to propagate.
- * @param srf The cylindrical surface to which to propagate.
- * @param dir The direction in which to propagate.
- * @param deriv The track derivatives to update at the surface srf.
- * @return The propagation status.
- */
- public PropStat dcaXYPropagate(double _bfac,
- VTrack trv,
- Surface srf,
- PropDir dir,
- TrackDerivative deriv) {
-// System.out.println("running PropDCAXY");
-
- // Construct return status
- PropStat pstat = new PropStat();
-
- // Fetch the originating Surface and check it is a DCA surface
- Surface srf0 = trv.surface();
- Assert.assertTrue(srf0.pureType().equals(SurfDCA.staticType()));
- if (!srf0.pureType().equals(SurfDCA.staticType())) {
- System.out.println("Failing PropDCAXY because it's not a DCA");
- return pstat;
- }
- SurfDCA srf_dca = (SurfDCA) srf0;
-
- // Check that dca surface has zero tilt.
- boolean tilted = srf_dca.dXdZ() != 0 || srf_dca.dYdZ() != 0;
- Assert.assertTrue(!tilted);
- if (tilted) {
- System.out.println("Failing PropDCAXY because it's tilted");
- return pstat;
- }
- // Check that the destination surface is an XYPlane
- Assert.assertTrue(srf.pureType().equals(SurfXYPlane.staticType()));
- if (!srf.pureType().equals(SurfXYPlane.staticType())) {
- System.out.println("Failing PropDCAXY because it's not a XY plane");
- return pstat;
- }
- SurfXYPlane srf_xyp = (SurfXYPlane) srf;
-
- // Instantiate the intermediate propagators.
- PropDCACyl prop1 = new PropDCACyl(bField());
- PropCylXY prop2 = new PropCylXY(bField());
-
- // Radius of a cylinder going through the PCA ...
- // But:
- // PropDCACyl requires that the cylinder be at a different radius than the DCA
- // so put it 100 microns outside the dca.
- double rhack = Math.abs(trv.vector(IRSIGNED)) + 0.01;
-
- // Instantiate the surface at the intermediate step.
- SurfCylinder srf_cyl = new SurfCylinder(rhack);
-
- // Setup for receiving the derivatives.
- TrackDerivative d1 = null;
- TrackDerivative d2 = null;
- if (deriv != null) {
- d1 = new TrackDerivative();
- d2 = new TrackDerivative();
- }
-
- // Save the track z direction for later use.
- boolean forward = trv.vector(ITLM_DCA) >= 0.;
-
-
-
- // Do the propagation in two steps.
- PropStat p1 = prop1.vecDirProp(trv, srf_cyl, dir, d1);
- if (!p1.success()) {
- System.out.println("Failing PropDCAXY propagation PropDCACyl failed");
- System.out.println("trv = "+trv.toString());
- System.out.println("cylinder = "+srf_cyl.toString());
-// trv.setTrackBackward();
-// p1=prop1.vecDirProp(trv, srf_cyl, dir, d1);
- if(!p1.success())return p1;
-// System.out.println("SUCCESS!!!!!!");
- }
- // System.out.println("propagated to fake cylinder, track is: " + trv.toString());
- PropStat p2 = prop2.vecDirProp(trv, srf_xyp, dir, d2);
-// System.out.println("propagated to plane, track is: " + trv.toString());
- if (!p2.success()) {
- System.out.println("Failing PropDCAXY propagation PropCylXY failed");
- return p2;
- }
-
- // Forward/backward is defined wrt positive z axis.
-
- //mg...I'm not sure if this is what we want...seems to screw up PropXYXY
- //mg (months later) ... how can this screw up PropXYXY???
-/*
- if ( forward ){
- trv.setForward();
- } else{
- trv.setBackward();
- }
- */
- trv.setForward(); //mg set to setForward by default
-
- // Set the transformation matrix to the product of the those from
- // the two intermediate steps.
- if (deriv != null) {
- deriv.set(d2.times(d1));
- }
-
- // Set properties of the return value.
- if (p2.forward()) {
- pstat.setForward();
- } else if (p2.backward()) {
- pstat.setBackward();
- } else if (p2.same()) {
- pstat.setSame();
- }
- double s = p1.pathDistance() + p2.pathDistance();
- pstat.setPathDistance(s);
-
-
-
- return pstat;
-
- }
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PropXYCyl.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PropXYCyl.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,627 +0,0 @@
-package org.hps.recon.tracking.kalman;
-
-import org.lcsim.recon.tracking.trfbase.PropDir;
-import org.lcsim.recon.tracking.trfbase.PropDirected;
-import org.lcsim.recon.tracking.trfbase.PropStat;
-import org.lcsim.recon.tracking.trfbase.Propagator;
-import org.lcsim.recon.tracking.trfbase.Surface;
-import org.lcsim.recon.tracking.trfbase.TrackDerivative;
-import org.lcsim.recon.tracking.trfbase.TrackVector;
-import org.lcsim.recon.tracking.trfbase.VTrack;
-import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
-import org.lcsim.recon.tracking.trfutil.Assert;
-import org.lcsim.recon.tracking.trfutil.TRFMath;
-import org.lcsim.recon.tracking.trfxyp.SurfXYPlane;
-
-/**
- * Propagates tracks from an XYPlane to a Cylinder in a constant field.
- *<p>
- * Propagation will fail if either the origin is not an XYPlane
- * or destination is not a Cylinder.
- * Propagator works incorrectly for tracks with very small curvatures.
- *<p>
- *@author Norman A. Graf
- *@version 1.0
- *
- */
-public class PropXYCyl extends PropDirected
-{
-
- // Assign track parameter indices.
-
- private static final int IV = SurfXYPlane.IV;
- private static final int IZC = SurfXYPlane.IZ;
- private static final int IDVDU = SurfXYPlane.IDVDU;
- private static final int IDZDU = SurfXYPlane.IDZDU;
- private static final int IQP_XY = SurfXYPlane.IQP;
-
- private static final int IPHI = SurfCylinder.IPHI;
- private static final int IZ = SurfCylinder.IZ;
- private static final int IALF = SurfCylinder.IALF;
- private static final int ITLM = SurfCylinder.ITLM;
- private static final int IQPT = SurfCylinder.IQPT;
-
-
-
- // attributes
-
- // BFAC * bfield
- private double _bfac;
-
- // static methods
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public static String typeName()
- { return "PropXYCyl";
- }
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public static String staticType()
- { return typeName();
- }
-
-
- /**
- *Construct an instance from a constant solenoidal magnetic field in Tesla.
- *
- * @param bfield The magnetic field strength in Tesla.
- */
- public PropXYCyl(double bfield)
- {
- _bfac = TRFMath.BFAC*bfield;
- }
-
- /**
- *Clone an instance.
- *
- * @return A Clone of this instance.
- */
- public Propagator newPropagator( )
- {
- return new PropXYCyl( bField() );
- }
-
- /**
- *Propagate a track without error in the specified direction.
- *
- * The track parameters for a cylinder are:
- * phi z alpha tan(lambda) curvature
- *
- * @param trv The VTrack to propagate.
- * @param srf The Surface to which to propagate.
- * @param dir The direction in which to propagate.
- * @return The propagation status.
- */
- public PropStat vecDirProp(VTrack trv, Surface srf,
- PropDir dir )
- {
- TrackDerivative deriv = null;
- return vecDirProp(trv, srf, dir, deriv);
- }
-
- /**
- *Propagate a track without error in the specified direction
- *and return the derivative matrix in deriv.
- *
- * The track parameters for a cylinder are:
- * phi z alpha tan(lambda) curvature
- *
- * @param trv The VTrack to propagate.
- * @param srf The Surface to which to propagate.
- * @param dir The direction in which to propagate.
- * @param deriv The track derivatives to update at the surface srf.
- * @return The propagation status.
- */
- public PropStat vecDirProp(VTrack trv, Surface srf,
- PropDir dir, TrackDerivative deriv )
- {
- return vecPropagateXYCyl( _bfac, trv, srf, dir, deriv);
-
- }
-
- /**
- *Return the strength of the magnetic field in Tesla.
- *
- * @return The strength of the magnetic field in Tesla.
- */
- public double bField()
- {
- return _bfac/TRFMath.BFAC;
- }
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public String type()
- { return staticType();
- }
-
- /**
- *output stream
- * @return A String representation of this instance.
- */
- public String toString()
- {
- return "XYPlane-Cylinder propagation with constant "
- + bField() +" Tesla field";
-
- }
-
-
- //**********************************************************************
-
- // Private function to propagate a track without error
- //
- // The track parameters for a cylinder are:
- // phi z alpha tan(lambda) curvature
- // (NOT [. . . lambda .] as in TRF and earlier version of TRF++.)
- //
- // If pderiv is nonzero, return the derivative matrix there.
- // On Cylinder:
- // r (cm) is fixed
- // 0 - phi
- // 1 - z (cm)
- // 2 - alp
- // 3 - tlam
- // 4 - q/pt pt is transverse momentum of a track, q is its charge
- // On XYPlane:
- // u (cm) is fixed
- // 0 - v (cm)
- // 1 - z (cm)
- // 2 - dv/du
- // 3 - dz/du
- // 4 - q/p p is momentum of a track, q is its charge
-
- PropStat
- vecPropagateXYCyl( double bfac, VTrack trv, Surface srf,
- PropDir dir,
- TrackDerivative deriv )
- {
-// System.out.println("In vecPropagateXYCyl");
- // construct return status
- PropStat pstat = new PropStat();
-
- // fetch the originating surface and vector
- Surface srf1 = trv.surface();
- TrackVector vec1 = trv.vector();
-
- // Check origin is a XYPlane.
- Assert.assertTrue( srf1.pureType().equals(SurfXYPlane.staticType()) );
- if ( !srf1.pureType( ).equals(SurfXYPlane.staticType()) )
- return pstat;
- SurfXYPlane sxyp1 = ( SurfXYPlane) srf1;
-// System.out.println("In vecPropagateXYCyl: XYPlane is ok");
-
- // Check destination is a cylinder.
- Assert.assertTrue( srf.pureType().equals(SurfCylinder.staticType()) );
- if ( !srf.pureType( ).equals(SurfCylinder.staticType()) )
- return pstat;
- SurfCylinder scy2 = ( SurfCylinder) srf;
-// System.out.println("In vecPropagateXYCyl: Cylinder Surface is ok");
-
- // Fetch the dist and phi of the plane and the starting track vector.
- int iphi = SurfXYPlane.NORMPHI;
- int idist = SurfXYPlane.DISTNORM;
- double phi = sxyp1.parameter(iphi);
- double r = sxyp1.parameter(idist);
-
- TrackVector vec = trv.vector();
- double v = vec.get(IV); // v
- double z = vec.get(IZC); // z
- double b = vec.get(IDVDU); // dv/du
- double a = vec.get(IDZDU); // dz/du
- double e = vec.get(IQP_XY); // q/p
-
- // Fetch the radii and the starting track vector.
- int irad = SurfCylinder.RADIUS;
- double r2 = scy2.parameter(irad);
-
- int sign_du = 0;
- if(trv.isForward()) sign_du = 1;
- if(trv.isBackward()) sign_du = -1;
- if(sign_du == 0)
- {
- System.out.println("PropXYCyl._vec_propagate: Unknown direction of a track ");
- return pstat;
- }
-
- // Calculate cylindrical coordinates
-
- double cnst=sign_du>0?0.: Math.PI;
- double atn=Math.atan(v/r);
- Assert.assertTrue(r!=0.);
-
- double phi1= TRFMath.fmod2(phi+atn,TRFMath.TWOPI);
- double r1=Math.sqrt(r*r+v*v);
- double z1=z;
- double alp1= TRFMath.fmod2(Math.atan(b)-atn+cnst,TRFMath.TWOPI);
- double tlm1= a*sign_du/Math.sqrt(1+b*b);
- double qpt1=e*Math.sqrt((1+a*a+b*b)/(1+b*b));
-
-
- // Check alpha range.
- alp1 = TRFMath.fmod2( alp1, TRFMath.TWOPI );
- Assert.assertTrue( Math.abs(alp1) <= Math.PI );
- //if ( trv.is_forward() ) Assert.assertTrue( Math.abs(alp1) <= TRFMath.PI2 );
- //else Assert.assertTrue( Math.abs(alp1) > TRFMath.PI2 );
-
- // Calculate the cosine of lambda.
- //double clam1 = 1.0/Math.sqrt(1+tlm1*tlm1);
-
- // Calculate curvature: C = _bfac*(q/p)/Math.cos(lambda)
- // and its derivatives
- // Assert.assertTrue( clam1 != 0.0 );
- // double dcrv1_dqp1 = bfac/clam1;
- // double crv1 = dcrv1_dqp1*qp1;
- // double dcrv1_dtlm1 = crv1*clam1*clam1*tlm1;
-
- // Calculate the curvature = _bfac*(q/pt)
- double dcrv1_dqpt1 = bfac;
- double crv1 = dcrv1_dqpt1*qpt1;
- //double dcrv1_dtlm1 = 0.0;
-
- // Evaluate the new track vector.
- // See dla log I-044
-
- // lambda and curvature do not change
- double tlm2 = tlm1;
- double crv2 = crv1;
- double qpt2 = qpt1;
-
- // We can evaluate Math.sin(alp2), leaving two possibilities for alp2
- // 1st solution: alp21, phi21, phid21, tht21
- // 2nd solution: alp22, phi22, phid22, tht22
- // evaluate phi2 to choose
- double salp1 = Math.sin( alp1 );
- double calp1 = Math.cos( alp1 );
- double salp2 = r1/r2*salp1 + 0.5*crv1/r2*(r2*r2-r1*r1);
- // if salp2 > 1, track does not cross cylinder
-// System.out.println("In vecPropagateXYCyl: Math.abs(salp2) = "+Math.abs(salp2));
- if ( Math.abs(salp2) > 1.0 ) return pstat;
- double alp21 = Math.asin( salp2 );
- double alp22 = alp21>0 ? Math.PI-alp21 : -Math.PI-alp21;
- double calp21 = Math.cos( alp21 );
- double calp22 = Math.cos( alp22 );
- double phi20 = phi1 + Math.atan2( salp1-r1*crv1, calp1 );
- double phi21 = phi20 - Math.atan2( salp2-r2*crv2, calp21 ); // phi position
- double phi22 = phi20 - Math.atan2( salp2-r2*crv2, calp22 );
-
- // Construct an sT object for each solution.
- STCalcXY sto1 = new STCalcXY(r1,phi1,alp1,crv1,r2,phi21,alp21);
- STCalcXY sto2 = new STCalcXY(r1,phi1,alp1,crv1,r2,phi22,alp22);
- // Check the two solutions are nonzero and have opposite sign
- // or at least one is nonzero.
-
- // Choose the correct solution
- boolean use_first_solution = false;
-
- if (dir.equals(PropDir.NEAREST))
- {
- use_first_solution = Math.abs(sto2.st()) > Math.abs(sto1.st());
- }
- else if (dir.equals(PropDir.FORWARD))
- {
- use_first_solution = sto1.st() > 0.0;
- }
- else if (dir.equals(PropDir.BACKWARD))
- {
- use_first_solution = sto1.st() < 0.0;
- }
- else
- {
- System.out.println("PropCyl._vec_propagate: Unknown direction.");
- System.exit(1);
- }
-
- // Assign phi2, alp2 and sto2 for the chosen solution.
- double phi2, alp2;
- STCalcXY sto;
- double calp2;
- if ( use_first_solution )
- {
- sto = sto1;
- phi2 = phi21;
- alp2 = alp21;
- calp2 = calp21;
- }
- else
- {
- sto = sto2;
- phi2 = phi22;
- alp2 = alp22;
- calp2 = calp22;
- }
-
- // fetch sT.
- double st = sto.st();
-
- // use sT to evaluate z2
- double z2 = z1 + tlm1*st;
-
- // Check alpha range.
- Assert.assertTrue( Math.abs(alp2) <= Math.PI );
-
- // put new values in vec
- vec.set(IPHI , phi2);
- vec.set(IZ , z2);
- vec.set(IALF , alp2);
- vec.set(ITLM , tlm2);
- vec.set(IQPT , qpt2);
-
- // Update trv
- trv.setSurface(srf.newPureSurface());
- trv.setVector(vec);
- if ( Math.abs(alp2) <= TRFMath.PI2 ) trv.setForward();
- else trv.setBackward();
-
- // Set the return status.
- double s = st*Math.sqrt(1.0+tlm2*tlm2);
- pstat.setPathDistance(s);
- //st > 0 ? pstat.set_forward() : pstat.set_backward();
-
- // exit now if user did not ask for error matrix.
- if ( deriv == null ) return pstat;
-
- // Calculate derivatives.
- // dphi1
-
- double dphi1_dv= r/(r*r+v*v);
-
- // dz1
-
- double dz1_dz=1.;
-
- // dalf1
-
- double dalp1_db= 1./(1.+b*b);
- double dalp1_dv= -r/(r*r+v*v);
-
- // dr1
-
- double dr1_dv= v/Math.sqrt(v*v+r*r);
-
- // dtlm1
-
- double dtlm1_da= sign_du/Math.sqrt(1+b*b);
- double dtlm1_db= -a*sign_du*b/(1+b*b)/Math.sqrt(1+b*b);
-
- // dcrv1
-
- double dqpt1_de= Math.sqrt((1+a*a+b*b)/(1+b*b));
- double dqpt1_da= a*e/Math.sqrt((1+b*b)*(1+a*a+b*b));
- double dqpt1_db= -e*b*a*a/Math.sqrt(1+a*a+b*b)/Math.sqrt(1+b*b)/(1+b*b);
-
- double dcrv1_de= dqpt1_de*bfac;
- double dcrv1_da= dqpt1_da*bfac;
- double dcrv1_db= dqpt1_db*bfac;
-
- // alpha_2
- double da2da1 = r1*calp1/r2/calp2;
- double da2dc1 = (r2*r2-r1*r1)*0.5/r2/calp2;
- double da2dr1 = (salp1-crv2*r1)/r2/calp2;
-
- // phi2
- double rcsal1 = r1*crv1*salp1;
- double rcsal2 = r2*crv2*salp2;
- double den1 = 1.0 + r1*r1*crv1*crv1 - 2.0*rcsal1;
- double den2 = 1.0 + r2*r2*crv2*crv2 - 2.0*rcsal2;
- double dp2dp1 = 1.0;
- double dp2da1 = (1.0-rcsal1)/den1 - (1.0-rcsal2)/den2*da2da1;
- double dp2dc1 = -r1*calp1/den1 + r2*calp2/den2
- - (1.0-rcsal2)/den2*da2dc1;
- double dp2dr1= -crv1*calp1/den1-(1.0-rcsal2)*da2dr1/den2;
-
- // z2
- double dz2dz1 = 1.0;
- double dz2dl1 = st;
- double dz2da1 = tlm1*sto.d_st_dalp1(dp2da1, da2da1);
- double dz2dc1 = tlm1*sto.d_st_dcrv1(dp2dc1, da2dc1);
- double dz2dr1 = tlm1*sto.d_st_dr1( dp2dr1, da2dr1);
-
- // final derivatives
-
- // phi2
- double dphi2_dv=dp2dp1*dphi1_dv+dp2da1*dalp1_dv+dp2dr1*dr1_dv;
- double dphi2_db=dp2da1*dalp1_db+dp2dc1*dcrv1_db;
- double dphi2_da=dp2dc1*dcrv1_da;
- double dphi2_de=dp2dc1*dcrv1_de;
-
- // alp2
- double dalp2_dv= da2da1*dalp1_dv+da2dr1*dr1_dv;
- double dalp2_db= da2da1*dalp1_db+da2dc1*dcrv1_db;
- double dalp2_da= da2dc1*dcrv1_da;
- double dalp2_de= da2dc1*dcrv1_de;
-
- // crv2
- double dqpt2_da=dqpt1_da;
- double dqpt2_db=dqpt1_db;
- double dqpt2_de=dqpt1_de;
-
-
- // tlm2
- double dtlm2_da= dtlm1_da;
- double dtlm2_db= dtlm1_db;
-
- // z2
- double dz2_dz= dz2dz1*dz1_dz;
- double dz2_dv= dz2dr1*dr1_dv+dz2da1*dalp1_dv;
- double dz2_db= dz2da1*dalp1_db+dz2dl1*dtlm1_db+dz2dc1*dcrv1_db;
- double dz2_da= dz2dl1*dtlm1_da+dz2dc1*dcrv1_da;
- double dz2_de= dz2dc1*dcrv1_de;
-
-
- // Build derivative matrix.
-
- deriv.set(IPHI,IV , dphi2_dv);
- deriv.set(IPHI,IDVDU, dphi2_db);
- deriv.set(IPHI,IDZDU, dphi2_da);
- deriv.set(IPHI,IQP_XY, dphi2_de);
- deriv.set(IZ,IV , dz2_dv);
- deriv.set(IZ,IZC , dz2_dz);
- deriv.set(IZ,IDVDU , dz2_db);
- deriv.set(IZ,IDZDU , dz2_da);
- deriv.set(IZ,IQP_XY, dz2_de);
- deriv.set(IALF,IV , dalp2_dv);
- deriv.set(IALF,IDVDU , dalp2_db);
- deriv.set(IALF,IDZDU , dalp2_da);
- deriv.set(IALF,IQP_XY , dalp2_de);
- deriv.set(ITLM,IDVDU , dtlm2_db);
- deriv.set(ITLM,IDZDU , dtlm2_da);
- deriv.set(IQPT,IDVDU , dqpt2_db);
- deriv.set(IQPT,IDZDU , dqpt2_da);
- deriv.set(IQPT,IQP_XY , dqpt2_de);
-
- return pstat;
-
- }
-
- //**********************************************************************
- // helpers
- //**********************************************************************
-
- // Private class STCalcXY.
- //
- // An STCalcXY_ object calculates sT (the signed transverse path length)
- // and its derivatives w.r.t. alf1 and crv1. It is constructed from
- // the starting (r1, phi1, alf1, crv1) and final track parameters
- // (r2, phi2, alf2) assuming these are consistent. Methods are
- // provided to retrieve sT and the two derivatives.
-
- class STCalcXY
- {
-
- private boolean _big_crv;
- private double _st;
- private double _dst_dphi21;
- private double _dst_dcrv1;
- private double _dst_dr1;
- private double _cnst1,_cnst2;
- public double _crv1;
-
- // constructor
- public STCalcXY()
- {
- }
- public STCalcXY(double r1, double phi1, double alf1, double crv1,
- double r2, double phi2, double alf2)
- {
-
- _crv1 = crv1;
- Assert.assertTrue( r1 > 0.0 );
- Assert.assertTrue( r2 > 0.0 );
- double rmax = r1+r2;
-
- // Calculate the change in xy direction
- double phi_dir_diff = TRFMath.fmod2(phi2+alf2-phi1-alf1,TRFMath.TWOPI);
- Assert.assertTrue( Math.abs(phi_dir_diff) <= Math.PI );
-
- // Evaluate whether |C| is" big"
- _big_crv = rmax*Math.abs(crv1) > 0.001 || Math.abs(r1-r2)<1.e-9;
- if( Math.abs(crv1) < 1.e-10 ) _big_crv=false;
-
- // If the curvature is big we can use
- // sT = (phi_dir2 - phi_dir1)/crv1
- if ( _big_crv )
- {
- Assert.assertTrue( crv1 != 0.0 );
- _st = phi_dir_diff/crv1;
- }
-
- // Otherwise, we calculate the straight-line distance
- // between the points and use an approximate correction
- // for the (small) curvature.
- else
- {
-
- // evaluate the distance
- double d = Math.sqrt( r1*r1 + r2*r2 - 2.0*r1*r2*Math.cos(phi2-phi1) );
- double arg = 0.5*d*crv1;
- double arg2 = arg*arg;
- double st_minus_d = d*arg2*( 1.0/6.0 + 3.0/40.0*arg2 );
- _st = d + st_minus_d;
-
- // evaluate the sign
- // We define a metric xsign = abs( (dphid-d*C)/(d*C) ).
- // Because sT*C = dphid and d = abs(sT):
- // xsign = 0 for sT > 0
- // xsign = 2 for sT < 0
- // Numerical roundoff will smear these predictions.
- double xsign = (crv1 == 0. ? 0.: Math.abs( (phi_dir_diff - _st*crv1) / (_st*crv1)) );
- double sign = 0.0;
- if ( crv1 != 0. )
- {
- if ( xsign < 0.5 ) sign = 1.0;
- if ( xsign > 1.5 && xsign < 3.0 ) sign = -1.0;
- }
- // If the above is indeterminate, assume zero curvature.
- // In this case abs(alpha) decreases monotonically
- // with sT. Track passing through origin has alpha = 0 on one
- // side and alpha = +/-pi on the other. If both points are on
- // the same side, we use dr/ds > 0 for |alpha|<pi/2.
- if ( sign == 0. )
- {
- sign = 1.0;
- if ( Math.abs(alf2) > Math.abs(alf1) ) sign = -1.0;
- if ( Math.abs(alf2) == Math.abs(alf1) )
- {
- if ( Math.abs(alf2) < TRFMath.PI2 )
- {
- if ( r2 < r1 ) sign = -1.0;
- }
- else
- {
- if ( r2 > r1 ) sign = -1.0;
- }
- }
- }
-
- // Correct _st using the above sign.
- Assert.assertTrue( Math.abs(sign) == 1.0 );
- _st = sign*_st;
-
- // save derivatives
- _dst_dcrv1 = sign*d*d*arg*( 1.0/6.0 + 3.0/20.0*arg2);
- double root = (1.0 + 0.5*arg*arg + 3.0/8.0*arg*arg*arg*arg );
- _dst_dphi21 = sign*(r1*r2*Math.sin(phi2-phi1))*root/d;
- _dst_dr1= (1.+arg2/2.*(1+3./4.*arg2))/d*sign;
- _cnst1=r1-r2*Math.cos(phi2-phi1);
- _cnst2=r1*r2*Math.sin(phi2-phi1);
- }
- }
-
- public double st()
- { return _st;
- }
-
- public double d_st_dalp1(double d_phi2_dalf1, double d_alf2_dalf1 )
- {
- if ( _big_crv ) return ( d_phi2_dalf1 + d_alf2_dalf1 - 1.0 ) / _crv1;
- else return _dst_dphi21 * d_phi2_dalf1;
-
- }
- public double d_st_dcrv1(double d_phi2_dcrv1, double d_alf2_dcrv1 )
- {
- if ( _big_crv ) return ( d_phi2_dcrv1 + d_alf2_dcrv1 - _st ) / _crv1;
- else return _dst_dcrv1 + _dst_dphi21*d_phi2_dcrv1;
-
- }
- public double d_st_dr1( double d_phi2_dr1, double d_alf2_dr1 )
- {
- if ( _big_crv ) return ( d_phi2_dr1 + d_alf2_dr1 ) / _crv1;
- else return _dst_dr1*(_cnst1+_cnst2*d_phi2_dr1);
- }
- }
-}
\ No newline at end of file
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PropXYDCA.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PropXYDCA.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,306 +0,0 @@
-package org.hps.recon.tracking.kalman;
-
-// May need to remove next line when moved to proper package.
-import org.lcsim.recon.tracking.trfbase.PropDir;
-import org.lcsim.recon.tracking.trfbase.PropDirected;
-import org.lcsim.recon.tracking.trfbase.PropStat;
-import org.lcsim.recon.tracking.trfbase.Propagator;
-import org.lcsim.recon.tracking.trfbase.Surface;
-import org.lcsim.recon.tracking.trfbase.TrackDerivative;
-import org.lcsim.recon.tracking.trfbase.VTrack;
-import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
-//import org.lcsim.recon.tracking.trfcylplane.PropXYCyl;
-import org.lcsim.recon.tracking.trfdca.PropCylDCA;
-import org.lcsim.recon.tracking.trfdca.SurfDCA;
-import org.lcsim.recon.tracking.trfutil.Assert;
-import org.lcsim.recon.tracking.trfutil.TRFMath;
-import org.lcsim.recon.tracking.trfxyp.SurfXYPlane;
-
-
-/**
- * Propagates tracks from an XY surface to an DCA surface.
- * It does so by creating an intermediate Cylinder surface,
- * and using existing code to perform the transformation in two steps.
- *<p>
- * Propagation will fail if either the origin is not an XY surface,
- * or the destination is not a DCA surface.
- *<p>
- * The default direction for propagation is forward.
- *<p>
- * PropDCACyl and PropCylXY do not work in all circumstances. See those
- * codes for details.
- *
- *
- *@author $Author: mgraham $
- *@version $Id: PropXYDCA.java,v 1.1 2011/07/07 20:57:38 mgraham Exp $
- *
- * Date $Date: 2011/07/07 20:57:38 $
- *
- */
-public class PropXYDCA extends PropDirected
-{
-
- // static variables
- // Assign track parameter indices
-
- public static final int IRSIGNED = SurfDCA.IRSIGNED;
- public static final int IZ_DCA = SurfDCA.IZ;
- public static final int IPHID = SurfDCA.IPHID;
- public static final int ITLM_DCA = SurfDCA.ITLM;
- public static final int IQPT_DCA = SurfDCA.IQPT;
-
- public static final int IPHI = SurfCylinder.IPHI;
- public static final int IZ_CYL = SurfCylinder.IZ;
- public static final int IALF = SurfCylinder.IALF;
- public static final int ITLM_CYL = SurfCylinder.ITLM;
- public static final int IQPT_CYL = SurfCylinder.IQPT;
-
- private static final int IV = SurfXYPlane.IV;
- private static final int IZC = SurfXYPlane.IZ;
- private static final int IDVDU = SurfXYPlane.IDVDU;
- private static final int IDZDU = SurfXYPlane.IDZDU;
- private static final int IQP_XY = SurfXYPlane.IQP;
-
-
- // Attributes ***************************************************
-
- // bfield * BFAC
- private double _bfac;
-
-
- // Methods ******************************************************
-
- // static methods
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public static String typeName()
- { return "PropXYDCA";
- }
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public static String staticType()
- { return typeName();
- }
-
-
-
- /**
- *Construct an instance from a constant solenoidal magnetic field in Tesla.
- *
- * @param bfield The magnetic field strength in Tesla.
- */
- public PropXYDCA(double bfield)
- {
- super(PropDir.FORWARD);
- _bfac=bfield * TRFMath.BFAC;
- }
-
- /**
- *Clone an instance.
- *
- * @return A Clone of this instance.
- */
- public Propagator newPropagator()
- {
- return new PropXYDCA( bField() );
- }
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public String type()
- { return staticType();
- }
-
- /**
- *Propagate a track without error in the specified direction.
- *and return the derivative matrix in deriv.
- *
- * The track parameters for a cylinder are:
- * phi z alpha tan(lambda) curvature
- *
- * @param trv The VTrack to propagate.
- * @param srf The Surface to which to propagate.
- * @param dir The direction in which to propagate.
- * @param deriv The track derivatives to update at the surface srf.
- * @return The propagation status.
- */
- public PropStat vecDirProp( VTrack trv, Surface srf,
- PropDir dir, TrackDerivative deriv )
- {
- return dcaXYPropagate( _bfac, trv, srf, dir, deriv );
- }
-
- /**
- *Propagate a track without error in the specified direction.
- *
- * The track parameters for a cylinder are:
- * phi z alpha tan(lambda) curvature
- *
- * @param trv The VTrack to propagate.
- * @param srf The Surface to which to propagate.
- * @param dir The direction in which to propagate.
- * @return The propagation status.
- */
- public PropStat vecDirProp( VTrack trv, Surface srf,
- PropDir dir )
- {
- TrackDerivative deriv =null;
- return vecDirProp(trv, srf, dir, deriv);
-
- }
-
- // propagate a track with error in the specified direction
- // PropStat err_dir_prop( ETrack& tre, Surface& srf,
- // PropDir dir ) ;
-
- //
-
- /**
- *Return the strength of the magnetic field in Tesla.
- *
- * @return The strength of the magnetic field in Tesla.
- */
- public double bField()
- {
- return _bfac/TRFMath.BFAC;
- }
-
-
- /**
- *output stream
- * @return A String representation of this instance.
- */
- public String toString()
- {
- return "XYPlane propagation to a DCA surface with constant "
- + bField() + " Tesla field";
- }
-
-
-
- /**
- * Propagate from XY to dca.
- * Why is this public?
- *
- * @param _bfac The numerical factor (including the field)
- * @param trv The VTrack to propagate.
- * @param srf The cylindrical surface to which to propagate.
- * @param dir The direction in which to propagate.
- * @param deriv The track derivatives to update at the surface srf.
- * @return The propagation status.
- */
- public PropStat dcaXYPropagate( double _bfac,
- VTrack trv,
- Surface srf,
- PropDir dir,
- TrackDerivative deriv )
- {
-// System.out.println("running PropXYDCA");
-
- // Construct return status
- PropStat pstat = new PropStat();
-
- // Fetch the originating Surface and check it is a XYPlane surface
- Surface srf0 = trv.surface();
- Assert.assertTrue( srf0.pureType().equals(SurfXYPlane.staticType() ));
- if ( ! srf0.pureType( ).equals(SurfXYPlane.staticType()) )
- return pstat;
- SurfXYPlane srf_xy = ( SurfXYPlane ) srf0;
- //System.out.println("XY Plane surface is OK...");
- // check the final surface is a DCA surface
- Assert.assertTrue( srf.pureType().equals(SurfDCA.staticType() ));
- if ( ! srf.pureType( ).equals(SurfDCA.staticType()) )
- return pstat;
- SurfDCA srf_dca = ( SurfDCA ) srf;
- //System.out.println("DCA surface is a DCA surface...");
- // Check that dca surface has zero tilt.
- boolean tilted = srf_dca.dXdZ() != 0 || srf_dca.dYdZ() != 0;
- Assert.assertTrue(!tilted);
- if(tilted) return pstat;
- //System.out.println("...and is not tilted");
-
- // Instantiate the intermediate propagators.
- PropXYCyl prop1 = new PropXYCyl( bField() );
- PropCylDCA prop2 = new PropCylDCA( bField() );
-
-
- // Radius of a cylinder going through the PCA ...
- // But:
- // PropDCACyl requires that the cylinder be at a different radius than the DCA
- // so put it 100 microns outside the dca.
-// double rhack = Math.abs(trv.vector(IRSIGNED))+0.01;
-// double rhack = Math.sqrt(srf_dca.x()*srf_dca.x()+srf_dca.y()*srf_dca.y())+0.01;
- //take this back to 10cm (layer 1)...if this is too small, track misses surface???
- double rhack = Math.sqrt(10.0);
- // Instantiate the surface at the intermediate step.
- SurfCylinder srf_cyl = new SurfCylinder(rhack);
-
- // Setup for receiving the derivatives.
- TrackDerivative d1 = null;
- TrackDerivative d2 = null;
- if ( deriv != null ){
- d1 = new TrackDerivative();
- d2 = new TrackDerivative();
- }
-
- // Save the track z direction for later use.
- boolean forward = trv.vector(ITLM_DCA)>=0.;
-
-
-//System.out.println("Propagating!");
- // Do the propagation in two steps.
- PropStat p1 = prop1.vecDirProp(trv, srf_cyl, dir, d1);
- if ( !p1.success() ) return p1;
-// System.out.println("propagated to fake cylinder, track is: " + trv.toString());
- PropStat p2 = prop2.vecDirProp(trv,srf_dca, dir, d2);
-// System.out.println("propagated to DCA, track is: " + trv.toString());
- if ( !p2.success() ) return p2;
-
- // Forward/backward is defined wrt positive z axis.
-/* mg...I'm not sure if this is what we want...seems to screw up PropXYXY
- if ( forward ){
- trv.setForward();
- } else{
- trv.setBackward();
- }
- */
-
- trv.setForward();
- // Set the transformation matrix to the product of the those from
- // the two intermediate steps.
- if ( deriv != null ){
- deriv.set(d2.times(d1));
- }
-
- // Set properties of the return value.
- if ( p2.forward() ){
- pstat.setForward();
- }else if ( p2.backward() ){
- pstat.setBackward();
- } else if ( p2.same() ){
- pstat.setSame();
- }
- double s = p1.pathDistance() + p2.pathDistance();
- pstat.setPathDistance(s);
-
-
-
- return pstat;
-
- }
-
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PropXYXY.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/PropXYXY.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,855 +0,0 @@
-package org.hps.recon.tracking.kalman;
-
-import org.lcsim.recon.tracking.trfbase.PropDir;
-import org.lcsim.recon.tracking.trfbase.PropDirected;
-import org.lcsim.recon.tracking.trfbase.PropStat;
-import org.lcsim.recon.tracking.trfbase.Propagator;
-import org.lcsim.recon.tracking.trfbase.Surface;
-import org.lcsim.recon.tracking.trfbase.TrackDerivative;
-import org.lcsim.recon.tracking.trfbase.TrackVector;
-import org.lcsim.recon.tracking.trfbase.VTrack;
-import org.lcsim.recon.tracking.trfutil.Assert;
-import org.lcsim.recon.tracking.trfutil.TRFMath;
-import org.lcsim.recon.tracking.trfxyp.SurfXYPlane;
-
-/**
- * Propagates tracks from one XYPlane to another in a constant field.
- *<p>
- * Propagation will fail if either the origin or destination is
- * not a XYPlane.
- * Propagator works incorrectly for tracks with very small curvatures
- *
- *
- *@author Norman A. Graf
- *@version 1.0
- *
- */
-public class PropXYXY extends PropDirected {
-
- // attributes
- private boolean _debug=false;
- private double _bfac;
- private static final int IV = SurfXYPlane.IV;
- private static final int IZ = SurfXYPlane.IZ;
- private static final int IDVDU = SurfXYPlane.IDVDU;
- private static final int IDZDU = SurfXYPlane.IDZDU;
- private static final int IQP = SurfXYPlane.IQP;
-
- // static methods
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public static String typeName() {
- return "PropXYXY";
- }
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public static String staticType() {
- return typeName();
- }
-
- /**
- *Construct an instance from a constant solenoidal magnetic field in Tesla.
- *
- * @param bfield The magnetic field strength in Tesla.
- */
- public PropXYXY(double bfield) {
- _bfac = TRFMath.BFAC * bfield;
- }
-
- /**
- *Clone an instance.
- *
- * @return A Clone of this instance.
- */
- public Propagator newPropagator() {
- return new PropXYXY(bField());
- }
-
- /**
- *Propagate a track without error in the specified direction
- *and return the derivative matrix in deriv.
- *
- * @param trv The VTrack to propagate.
- * @param srf The Surface to which to propagate.
- * @param dir The direction in which to propagate.
- * @param deriv The track derivatives to update at the surface srf.
- * @return The propagation status.
- **/
- public PropStat vecDirProp(VTrack trv, Surface srf,
- PropDir dir, TrackDerivative deriv) {
- PropStat pstat = vec_propagatexyxy_(_bfac, trv, srf, dir, deriv);
- return pstat;
- }
-
- /**
- *Propagate a track without error in the specified direction.
- *
- * @param trv The VTrack to propagate.
- * @param srf The Surface to which to propagate.
- * @param dir The direction in which to propagate.
- * @return The propagation status.
- */
- @Override
- public PropStat vecDirProp(VTrack trv, Surface srf,
- PropDir dir) {
- TrackDerivative deriv = null;
- return vecDirProp(trv, srf, dir, deriv);
- }
-
- /**
- *Return the strength of the magnetic field in Tesla.
- *
- * @return The strength of the magnetic field in Tesla.
- */
- public double bField() {
- return _bfac / TRFMath.BFAC;
- }
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public String type() {
- return staticType();
- }
-
- /**
- *output stream
- *
- * @return A String representation of this instance.
- */
- @Override
- public String toString() {
- return "XYPlane-XYPlane propagation with constant "
- + bField() + " Tesla field";
- }
-
- // Private function to determine dphi of the propagation
- double direction(int flag_forward, int sign_dphi,
- double du,
- double norm, double cat,
- double sinphi, double cosphi) {
-
- int sign_cat = 0;
- if (du * flag_forward > 0.)
- sign_cat = 1;
- if (du * flag_forward < 0.)
- sign_cat = -1;
- if (du == 0.) {
- if (sinphi >= 0.)
- sign_cat = 1;
- if (sinphi < 0.)
- sign_cat = -1;
- }
-
- double sin_dphi = norm * sinphi + cat * cosphi * sign_cat;
- double cos_dphi = -norm * cosphi + cat * sinphi * sign_cat;
-
- int sign_sindphi = 0;
- if (sin_dphi > 0.)
- sign_sindphi = 1;
- if (sin_dphi < 0.)
- sign_sindphi = -1;
- if (sin_dphi == 0.)
- sign_sindphi = sign_dphi;
-
- if (Math.abs(cos_dphi) > 1.) {
- double sn = -1.;
- if (cos_dphi > 0)
- sn = 1.;
- cos_dphi = sn * Math.sqrt(1. - sin_dphi * sin_dphi);
- }
- if (Math.abs(sin_dphi) > 1.) {
- double sn = -1.;
- if (sin_dphi > 0)
- sn = 1.;
- sin_dphi = sn * Math.sqrt(1. - cos_dphi * cos_dphi);
- }
-
- double dphi = Math.PI * (sign_dphi - sign_sindphi) + sign_sindphi * Math.acos(cos_dphi);
- return dphi;
-
- }
- //**********************************************************************
-
- // Private function to propagate a track without error
- // The corresponding track parameters are:
- // u (cm) is fixed
- // 0 - v (cm)
- // 1 - z (cm)
- // 2 - dv/du
- // 3 - dz/du
- // 4 - q/p p is momentum of a track, q is its charge
- // If deriv is nonzero, return the derivative matrix there.
- PropStat vec_propagatexyxy_(double B, VTrack trv, Surface srf,
- PropDir dir1,
- TrackDerivative deriv) {
-
- // construct return status
- PropStat pstat = new PropStat();
-
- PropDir dir = dir1; //need to check constness of this
- boolean move = Propagator.reduceDirection(dir);
- if (move)
- dir = reduce(dir);
- // fetch the originating surface and vector
- Surface srf1 = trv.surface();
- // TrackVector vec1 = trv.get_vector();
-
- // Check origin is a XYPlane.
- Assert.assertTrue(srf1.pureType().equals(SurfXYPlane.staticType()));
- if (!srf1.pureType().equals(SurfXYPlane.staticType()))
- return pstat;
- SurfXYPlane sxyp1 = (SurfXYPlane) srf1;
-
- // Check destination is an XYPlane.
- Assert.assertTrue(srf.pureType().equals(SurfXYPlane.staticType()));
- if (!srf.pureType().equals(SurfXYPlane.staticType()))
- return pstat;
- SurfXYPlane sxyp2 = (SurfXYPlane) srf;
-
-
- // If surfaces are the same and no XXX_MOVE is requested we can return now.
- boolean same = srf.pureEqual(srf1);
- if(_debug)System.out.println("same=" + same + "; move=" + move + "; " + dir1.toString() + " " + dir.toString());
-
- if (same && !move) {
- if (deriv != null) {
- deriv.setIdentity();
- }
- pstat.setSame();
- return pstat;
- }
-
-
- if (same && dir1.equals(PropDir.NEAREST_MOVE))
- dir = PropDir.FORWARD;
-
- //v00.63.04 cng 01/12/01
- if (Math.abs(B) < 1e-7)
- return zeroBField(trv, sxyp1, sxyp2, dir1, deriv);
-
- // Fetch the u's and phi's of the planes and the starting track vector.
- TrackVector vec = trv.vector();
- int iphi = SurfXYPlane.NORMPHI;
- int idist = SurfXYPlane.DISTNORM;
- double phi_0 = sxyp1.parameter(iphi);
- double u_0 = sxyp1.parameter(idist);
- double phi_n = sxyp2.parameter(iphi);
- double u_n = sxyp2.parameter(idist);
- if (same && move) {
- int forw = trv.isForward() ? 1 : -1;
- int forw_dir = (dir == PropDir.FORWARD) ? 1 : -1;
- int sn = vec.get(IDVDU) > 0. ? 1 : -1;
- u_0 += 1.e-7 * (double) (forw * forw_dir * sn);
- }
-
- double b1_0 = vec.get(IV); // v
- double b2_0 = vec.get(IZ); // z
- double b3_0 = vec.get(IDVDU); // dv/du
- double b4_0 = vec.get(IDZDU); // dz/du
- double b5_0 = vec.get(IQP); // q/p
-
- double phi_u = phi_0 - phi_n;
-
- double cosphi_u = Math.cos(phi_u);
- double sinphi_u = Math.sin(phi_u);
-// System.out.println( "PropXYXY._vec_propagate: phi_u = "+phi_u);
- // check if du == 0 ( that is track moves parallel to the destination plane )
- double du_du_0 = cosphi_u - b3_0 * sinphi_u;
- if (du_du_0 == 0.)
- return pstat;
-
- double a_hat_u = 1. / du_du_0;
- double a_hat_u2 = a_hat_u * a_hat_u;
-
- double u = u_0 * cosphi_u - b1_0 * sinphi_u;
- double b1 = b1_0 * cosphi_u + u_0 * sinphi_u;
- double b2 = b2_0;
- double b3 = (b3_0 * cosphi_u + sinphi_u) * a_hat_u;
- double b4 = b4_0 * a_hat_u;
- double b5 = b5_0;
-
- int sign_du0 = 0;
- if (trv.isForward())
- sign_du0 = 1;
- if (trv.isBackward())
- sign_du0 = -1;
- if (sign_du0 == 0) {
- System.out.println("PropXYXY._vec_propagate: Unknown direction of a track ");
- System.exit(1);
- }
- int sign_du = 0;
- if (du_du_0 * sign_du0 > 0)
- sign_du = 1;
- if (du_du_0 * sign_du0 < 0)
- sign_du = -1;
-
- // check that q/p != 0
- Assert.assertTrue(b5 != 0.);
-
- // 1/curv of the track is r
- double r = 1 / (b5 * B) * Math.sqrt(1 + b3_0 * b3_0) / Math.sqrt(1 + b3_0 * b3_0 + b4_0 * b4_0);
- double b3_hat = Math.sqrt(1 + b3 * b3);
- double b34_hat = Math.sqrt(1 + b3 * b3 + b4 * b4);
- double b3_hat2 = b3_hat * b3_hat;
- double b34_hat2 = b34_hat * b34_hat;
- double cosphi = -b3 * sign_du / b3_hat;
- double sinphi = sign_du / b3_hat;
- double rsinphi = 1. / (b5 * B) * sign_du / b34_hat;
- double rcosphi = -b3 / (b5 * B) * sign_du / b34_hat;
-
- double du = u_n - u;
- double norm = du / r - cosphi;
- if(_debug)System.out.println("PropXYXY._vec_propagate: u_n = " + u_n + "; u = " + u + "; r = " + r);
-
- // xy-xy propagation failed : noway to the new plane
- if (Math.abs(norm) > 1.)
- return pstat;
-
- double cat = Math.sqrt(1. - norm * norm);
- int flag_forward = 0;
- int sign_dphi = 0;
-
- if (dir.equals(PropDir.NEAREST)) {
- // try forward propagation
- flag_forward = 1;
- if (b5 * B > 0.)
- sign_dphi = 1;
- if (b5 * B < 0.)
- sign_dphi = -1;
- double dphi1 =
- direction(flag_forward, sign_dphi, du, norm, cat, sinphi, cosphi);
- // try backward propagation
- flag_forward = -1;
- if (b5 * B > 0.)
- sign_dphi = -1;
- if (b5 * B < 0.)
- sign_dphi = 1;
- double dphi2 =
- direction(flag_forward, sign_dphi, du, norm, cat, sinphi, cosphi);
- if (Math.abs(dphi2) > Math.abs(dphi1)) {
- flag_forward = -flag_forward;
- sign_dphi = -sign_dphi;
- }
- } else if (dir.equals(PropDir.FORWARD)) {
- flag_forward = 1;
- if (b5 * B > 0.)
- sign_dphi = 1;
- if (b5 * B < 0.)
- sign_dphi = -1;
- } else if (dir.equals(PropDir.BACKWARD)) {
- flag_forward = -1;
- if (b5 * B > 0.)
- sign_dphi = -1;
- if (b5 * B < 0.)
- sign_dphi = 1;
- } else {
- System.out.println("PropXYXY._vec_propagate: Unknown direction.");
- System.exit(1);
- }
-
- int sign_cat = 0;
- if (du * sign_dphi * b5 * B > 0.)
- sign_cat = 1;
- if (du * sign_dphi * b5 * B < 0.)
- sign_cat = -1;
- if (du == 0.) {
- if (sinphi >= 0.)
- sign_cat = 1;
- if (sinphi < 0.)
- sign_cat = -1;
- }
-
- double sin_dphi = norm * sinphi + cat * cosphi * sign_cat;
- double cos_dphi = -norm * cosphi + cat * sinphi * sign_cat;
- if (cos_dphi > 1.)
- cos_dphi = 1.;
- if (cos_dphi < -1.)
- cos_dphi = -1.;
-
- int sign_sindphi = 0;
- if (sin_dphi > 0.)
- sign_sindphi = 1;
- if (sin_dphi < 0.)
- sign_sindphi = -1;
- if (sin_dphi == 0.)
- sign_sindphi = sign_dphi;
-
- double dphi = Math.PI * (sign_dphi - sign_sindphi) + sign_sindphi * Math.acos(cos_dphi);
-
- // check that I didn't make any mistakes
-
- Assert.assertTrue(Math.abs(Math.sin(dphi) - sin_dphi) < 1.e-5);
-
- // check if dun == 0 ( that is track moves parallel to the destination plane)
- double du_n_du = cos_dphi - b3 * sin_dphi;
- if (du_n_du == 0.)
- return pstat;
-
- double a_hat_dphi = 1. / du_n_du;
- double a_hat_dphi2 = a_hat_dphi * a_hat_dphi;
- double c_hat_dphi = sin_dphi + b3 * cos_dphi;
-
- double b1_n = b1 + rsinphi * (1 - cos_dphi) - rcosphi * sin_dphi;
- double b2_n = b2 + b4 / (b5 * B) * sign_du / b34_hat * dphi;
- double b3_n = c_hat_dphi * a_hat_dphi;
- double b4_n = b4 * a_hat_dphi;
- double b5_n = b5;
- if(_debug)System.out.println(trv.toString());
- if(_debug) System.out.println("PropXYXY._vec_propagate: sign_du0=" + sign_du0 + "; du_du_0=" + du_du_0);
- if(_debug)System.out.println("PropXYXY._vec_propagate: sign_du=" + sign_du + "; b3_hat=" + b3_hat);
- if(_debug)System.out.println("PropXYXY._vec_propagate: sign_cat=" + sign_cat + "; cosphi=" + cosphi
- + "; sinphi=" + sinphi + "; cat=" + cat);
- if(_debug)System.out.println("PropXYXY._vec_propagate: a_hat_dphi=" + a_hat_dphi + "; cos_dphi=" + cos_dphi
- + "; norm=" + norm);
- if(_debug)System.out.println("PropXYXY._vec_propagate: b3=" + b3 + "; sin_dphi=" + sin_dphi);
- if(_debug)System.out.println("PropXYXY._vec_propagate: b4_0=" + b4_0 + "; b4=" + b4 + "; b4_n=" + b4_n);
-
- // double u_n_0 = u_n*cosphi_u + b1_n*sinphi_u;
-
- // check if track crossed original plane during the propagation
- // switch (dir) {
- // if( dir.equals(PropDir.FORWARD:
- // if((u_n_0 - u_0)*sign_du0<0) return pstat;
- // break;
- // if( dir.equals(PropDir.BACKWARD:
- // if((u_n_0 - u_0)*sign_du0>0) return pstat;
- // break;
- // }
-
- int sign_dun = 0;
- if (du_n_du * sign_du > 0)
- sign_dun = 1;
- if (du_n_du * sign_du < 0)
- sign_dun = -1;
-
- vec.set(IV, b1_n);
- vec.set(IZ, b2_n);
- vec.set(IDVDU, b3_n);
- vec.set(IDZDU, b4_n);
- vec.set(IQP, b5_n);
- if(_debug)System.out.println("PropXYXY._vec_propagate: sign_dun = "+sign_dun+"; du_n_du = "+du_n_du+"; sign_du="+sign_du);
-
- // Update trv
- trv.setSurface(srf.newPureSurface());
- trv.setVector(vec);
-
- // set new direction of the track
- if (sign_dun == 1)
- trv.setForward();
- if (sign_dun == -1)
- trv.setBackward();
-
- // Calculate sT.
- double crv = B * b5;
- double dv = b1_n - b1;
- double dxy = Math.sqrt(du * du + dv * dv);
- double arg = 0.5 * crv * dxy;
- double dst = dxy * TRFMath.asinrat(arg);
-
- // Calculate s.
- double dz = b2_n - b2;
- Assert.assertTrue(flag_forward == 1 || flag_forward == -1);
- double ds = ((double) (flag_forward)) * Math.sqrt(dst * dst + dz * dz);
-
- // Set the return status.
- pstat.setPathDistance(ds);
- //(flag_forward==1)?pstat.set_forward():pstat.set_backward();
-
- // exit now if user did not ask for error matrix.
- if (deriv == null)
- return pstat;
-
- // du_d0
-
- double du_db1_0 = -sinphi_u;
-
- // db1_d0
-
- double db1_db1_0 = cosphi_u;
-
- // db3_d0
-
- double db3_db3_0 = a_hat_u2;
-
- // db4_d0
-
- double db4_db3_0 = b4_0 * sinphi_u * a_hat_u2;
- double db4_db4_0 = a_hat_u;
-
- // dr_d
-
- double dr_db3 = r * b3 * b4 * b4 / (b3_hat2 * b34_hat2);
- double dr_db4 = -r * b4 / b34_hat2;
- double dr_db5 = -r / b5;
-
- // dcosphi_d
-
- double dcosphi_db3 = -sign_du / b3_hat - cosphi * b3 / b3_hat2;
-
- // dsinphi_d
-
- double dsinphi_db3 = -sinphi * b3 / b3_hat2;
-
- // dcat_d
-
- double dcat_db3 = norm / cat * (du / (r * r) * dr_db3 + dcosphi_db3);
- double dcat_db4 = norm / cat * du / (r * r) * dr_db4;
- double dcat_db5 = norm / cat * du / (r * r) * dr_db5;
- double dcat_du = norm / (cat * r);
-
- // dnorm_d
-
- double dnorm_db3 = -du / (r * r) * dr_db3 - dcosphi_db3;
- double dnorm_db4 = -du / (r * r) * dr_db4;
- double dnorm_db5 = -du / (r * r) * dr_db5;
- double dnorm_du = -1. / r;
-
- // dcos_dphi_d
-
- double dcos_dphi_db3 = -cosphi * dnorm_db3 - norm * dcosphi_db3
- + sign_cat * (sinphi * dcat_db3 + cat * dsinphi_db3);
- double dcos_dphi_db4 = -cosphi * dnorm_db4 + sign_cat * sinphi * dcat_db4;
- double dcos_dphi_db5 = -cosphi * dnorm_db5 + sign_cat * sinphi * dcat_db5;
- double dcos_dphi_du = -cosphi * dnorm_du + sign_cat * sinphi * dcat_du;
-
- // dsin_dphi_d
-
- double dsin_dphi_db3 = sinphi * dnorm_db3 + norm * dsinphi_db3
- + sign_cat * (cosphi * dcat_db3 + cat * dcosphi_db3);
- double dsin_dphi_db4 = sinphi * dnorm_db4 + sign_cat * cosphi * dcat_db4;
- double dsin_dphi_db5 = sinphi * dnorm_db5 + sign_cat * cosphi * dcat_db5;
- double dsin_dphi_du = sinphi * dnorm_du + sign_cat * cosphi * dcat_du;
-
- // ddphi_d
-
- double ddphi_db3;
- double ddphi_db4;
- double ddphi_db5;
- double ddphi_du;
- if (Math.abs(sin_dphi) > 0.5) {
- ddphi_db3 = -dcos_dphi_db3 / sin_dphi;
- ddphi_db4 = -dcos_dphi_db4 / sin_dphi;
- ddphi_db5 = -dcos_dphi_db5 / sin_dphi;
- ddphi_du = -dcos_dphi_du / sin_dphi;
- } else {
- ddphi_db3 = dsin_dphi_db3 / cos_dphi;
- ddphi_db4 = dsin_dphi_db4 / cos_dphi;
- ddphi_db5 = dsin_dphi_db5 / cos_dphi;
- ddphi_du = dsin_dphi_du / cos_dphi;
- }
-
- // da_hat_dphi_d
-
- double da_hat_dphi_db3 = -a_hat_dphi2
- * (dcos_dphi_db3 - sin_dphi - b3 * dsin_dphi_db3);
- double da_hat_dphi_db4 = -a_hat_dphi2 * (dcos_dphi_db4 - b3 * dsin_dphi_db4);
- double da_hat_dphi_db5 = -a_hat_dphi2 * (dcos_dphi_db5 - b3 * dsin_dphi_db5);
- double da_hat_dphi_du = -a_hat_dphi2 * (dcos_dphi_du - b3 * dsin_dphi_du);
-
- // dc_hat_dphi_d
-
- double dc_hat_dphi_db3 = b3 * dcos_dphi_db3 + dsin_dphi_db3 + cos_dphi;
- double dc_hat_dphi_db4 = b3 * dcos_dphi_db4 + dsin_dphi_db4;
- double dc_hat_dphi_db5 = b3 * dcos_dphi_db5 + dsin_dphi_db5;
- double dc_hat_dphi_du = b3 * dcos_dphi_du + dsin_dphi_du;
-
- // db1_n_d
-
- double db1_n_db1 = 1;
- double db1_n_db3 = (dr_db3 * sinphi + r * dsinphi_db3) * (1 - cos_dphi)
- - rsinphi * dcos_dphi_db3
- - dr_db3 * cosphi * sin_dphi - r * dcosphi_db3 * sin_dphi
- - rcosphi * dsin_dphi_db3;
- double db1_n_db4 = dr_db4 * sinphi * (1 - cos_dphi) - rsinphi * dcos_dphi_db4
- - dr_db4 * cosphi * sin_dphi - rcosphi * dsin_dphi_db4;
- double db1_n_db5 = dr_db5 * sinphi * (1 - cos_dphi) - rsinphi * dcos_dphi_db5
- - dr_db5 * cosphi * sin_dphi - rcosphi * dsin_dphi_db5;
- double db1_n_du = -rsinphi * dcos_dphi_du - rcosphi * dsin_dphi_du;
-
- // db2_n_d
-
- double db2_n_db2 = 1.;
- double db2_n_db3 = 1. / (b5 * B) * b4 * sign_du / b34_hat
- * (-dphi * b3 / b34_hat2 + ddphi_db3);
- double db2_n_db4 = 1. / (b5 * B) * sign_du / b34_hat
- * (-dphi * b4 * b4 / b34_hat2 + b4 * ddphi_db4 + dphi);
- double db2_n_db5 = 1. / (b5 * B) * b4 * sign_du / b34_hat * (ddphi_db5 - dphi / b5);
- double db2_n_du = 1. / (b5 * B) * b4 * sign_du / b34_hat * ddphi_du;
-
- // db3_n_d
-
- double db3_n_db3 = a_hat_dphi * dc_hat_dphi_db3 + da_hat_dphi_db3 * c_hat_dphi;
- double db3_n_db4 = a_hat_dphi * dc_hat_dphi_db4 + da_hat_dphi_db4 * c_hat_dphi;
- double db3_n_db5 = a_hat_dphi * dc_hat_dphi_db5 + da_hat_dphi_db5 * c_hat_dphi;
- double db3_n_du = a_hat_dphi * dc_hat_dphi_du + da_hat_dphi_du * c_hat_dphi;
-
- // db4_n_d
-
- double db4_n_db3 = b4 * da_hat_dphi_db3;
- double db4_n_db4 = b4 * da_hat_dphi_db4 + a_hat_dphi;
- double db4_n_db5 = b4 * da_hat_dphi_db5;
- double db4_n_du = b4 * da_hat_dphi_du;
-
- // db5_n_d
-
- // db1_n_d0
-
- double db1_n_db1_0 = db1_n_du * du_db1_0 + db1_n_db1 * db1_db1_0;
- double db1_n_db2_0 = 0.;
- double db1_n_db3_0 = db1_n_db3 * db3_db3_0 + db1_n_db4 * db4_db3_0;
- double db1_n_db4_0 = db1_n_db4 * db4_db4_0;
- double db1_n_db5_0 = db1_n_db5;
-
- // db2_n_d0
-
- double db2_n_db1_0 = db2_n_du * du_db1_0;
- double db2_n_db2_0 = db2_n_db2;
- double db2_n_db3_0 = db2_n_db3 * db3_db3_0 + db2_n_db4 * db4_db3_0;
- double db2_n_db4_0 = db2_n_db4 * db4_db4_0;
- double db2_n_db5_0 = db2_n_db5;
-
- // db3_n_d0
-
- double db3_n_db1_0 = db3_n_du * du_db1_0;
- double db3_n_db2_0 = 0.;
- double db3_n_db3_0 = db3_n_db3 * db3_db3_0 + db3_n_db4 * db4_db3_0;
- double db3_n_db4_0 = db3_n_db4 * db4_db4_0;
- double db3_n_db5_0 = db3_n_db5;
-
- // db4_n_d0
-
- double db4_n_db1_0 = db4_n_du * du_db1_0;
- double db4_n_db2_0 = 0.;
- double db4_n_db3_0 = db4_n_db3 * db3_db3_0 + db4_n_db4 * db4_db3_0;
- double db4_n_db4_0 = db4_n_db4 * db4_db4_0;
- double db4_n_db5_0 = db4_n_db5;
-
- // db5_n_d0
-
- double db5_n_db1_0 = 0.;
- double db5_n_db2_0 = 0.;
- double db5_n_db3_0 = 0.;
- double db5_n_db4_0 = 0.;
- double db5_n_db5_0 = 1.;
-
-
- deriv.set(IV, IV, db1_n_db1_0);
- deriv.set(IV, IZ, db1_n_db2_0);
- deriv.set(IV, IDVDU, db1_n_db3_0);
- deriv.set(IV, IDZDU, db1_n_db4_0);
- deriv.set(IV, IQP, db1_n_db5_0);
- deriv.set(IZ, IV, db2_n_db1_0);
- deriv.set(IZ, IZ, db2_n_db2_0);
- deriv.set(IZ, IDVDU, db2_n_db3_0);
- deriv.set(IZ, IDZDU, db2_n_db4_0);
- deriv.set(IZ, IQP, db2_n_db5_0);
- deriv.set(IDVDU, IV, db3_n_db1_0);
- deriv.set(IDVDU, IZ, db3_n_db2_0);
- deriv.set(IDVDU, IDVDU, db3_n_db3_0);
- deriv.set(IDVDU, IDZDU, db3_n_db4_0);
- deriv.set(IDVDU, IQP, db3_n_db5_0);
- deriv.set(IDZDU, IV, db4_n_db1_0);
- deriv.set(IDZDU, IZ, db4_n_db2_0);
- deriv.set(IDZDU, IDVDU, db4_n_db3_0);
- deriv.set(IDZDU, IDZDU, db4_n_db4_0);
- deriv.set(IDZDU, IQP, db4_n_db5_0);
- deriv.set(IQP, IV, db5_n_db1_0);
- deriv.set(IQP, IZ, db5_n_db2_0);
- deriv.set(IQP, IDVDU, db5_n_db3_0);
- deriv.set(IQP, IDZDU, db5_n_db4_0);
- deriv.set(IQP, IQP, db5_n_db5_0);
-
- return pstat;
-
- }
-
- //cng
- PropStat zeroBField(VTrack trv, SurfXYPlane sxyp1,
- SurfXYPlane sxyp2, PropDir dir1,
- TrackDerivative deriv) {
- // construct return status
- PropStat pstat = new PropStat();
-
- PropDir dir = dir1; //need to check constness of this
- boolean move = Propagator.reduceDirection(dir);
- boolean same = sxyp2.pureEqual(sxyp1);
-
- // There is only one solution. Can't XXX_MOVE
- if (same && move)
- return pstat;
-
- if (same) {
- if (deriv != null) {
- deriv.setIdentity();
- }
- pstat.setSame();
- return pstat;
- }
-
- TrackVector vec = trv.vector();
- double v0 = vec.get(IV);
- double z0 = vec.get(IZ);
- double dvdu0 = vec.get(IDVDU);
- double dzdu0 = vec.get(IDZDU);
-
- double du0 = 1.;
- if (trv.isBackward())
- du0 = -1.;
-
- double phi0 = sxyp1.parameter(SurfXYPlane.NORMPHI);
- double cphi0 = Math.cos(phi0);
- double sphi0 = Math.sin(phi0);
- double u0 = sxyp1.parameter(SurfXYPlane.DISTNORM);
-
- double phi1 = sxyp2.parameter(SurfXYPlane.NORMPHI);
- double cphi1 = Math.cos(phi1);
- double sphi1 = Math.sin(phi1);
- double u1 = sxyp2.parameter(SurfXYPlane.DISTNORM);
-
- double a = (cphi0 - dvdu0 * sphi0) * du0;
- double b = (sphi0 + dvdu0 * cphi0) * du0;
- double c = dzdu0 * du0;
-
- double x0 = u0 * cphi0 - v0 * sphi0;
- double y0 = u0 * sphi0 + v0 * cphi0;
-
- double ap = u1 * cphi1;
- double bp = u1 * sphi1;
- double cp = 0;
-
- double xp = ap;
- double yp = bp;
- double zp = 0.0;
-
- double denom = a * ap + b * bp + c * cp;
-
- if (denom == 0.)
- return pstat;
-
- double S = ((xp - x0) * ap + (yp - y0) * bp + (zp - z0) * cp) / denom;
-
- double x1 = x0 + S * a;
- double y1 = y0 + S * b;
- double z1 = z0 + S * c;
-
- boolean forward = S > 0. ? true : false;
-
- if (dir == PropDir.FORWARD && !forward)
- return pstat;
- if (dir == PropDir.BACKWARD && forward)
- return pstat;
-
- double v1 = y1 * cphi1 - x1 * sphi1;
-
- double v01 = y0 * cphi1 - x0 * sphi1;
- double u01 = y0 * sphi1 + x0 * cphi1;
- double z01 = z0;
-
- if (u01 == u1)
- return pstat;
-
- double dvdu1 = (v1 - v01) / (u1 - u01);
- double dzdu1 = (z1 - z01) / (u1 - u01);
-
- vec.set(IV, v1);
- vec.set(IZ, z1);
- vec.set(IDVDU, dvdu1);
- vec.set(IDZDU, dzdu1);
-
- // Update trv
- trv.setSurface(sxyp2.newPureSurface());
- trv.setVector(vec);
- // set new direction of the track
- if (b * sphi1 + a * cphi1 > 0)
- trv.setForward();
- else
- trv.setBackward();
-
- // Calculate s.
- double ds = S * Math.sqrt(a * a + b * b + c * c);
-
- // Set the return status.
- pstat.setPathDistance(ds);
-
- if (deriv == null)
- return pstat;
-
- double dx0dv0 = -sphi0;
- double dy0dv0 = cphi0;
-
- double dadv_du = -sphi0 * du0;
- double dbdv_du = cphi0 * du0;
- double dcdz_du = du0;
-
- double ddenomdv_du = dadv_du * ap + dbdv_du * bp;
- double ddenomdz_du = dcdz_du * cp;
-
-
- double dSdv0 = -(dx0dv0 * ap + dy0dv0 * bp) / denom;
- double dSdz0 = -cp / denom;
- double dSdv_du = -S / denom * ddenomdv_du;
- double dSdz_du = -S / denom * ddenomdz_du;
-
- double dy1dv0 = dy0dv0 + dSdv0 * b;
- double dx1dv0 = dx0dv0 + dSdv0 * a;
- double dx1dv_du = dSdv_du * a + S * dadv_du;
- double dy1dv_du = dSdv_du * b + S * dbdv_du;
-
- double du01dv0 = dy0dv0 * sphi1 + dx0dv0 * cphi1;
- double dv01dv0 = dy0dv0 * cphi1 - dx0dv0 * sphi1;
-
- double dv1dv0 = dy1dv0 * cphi1 - dx1dv0 * sphi1;
- double dv1dv_du = dy1dv_du * cphi1 - dx1dv_du * sphi1;
-
- double dz1dz0 = 1 + dSdz0 * c;
- double dz1dv0 = dSdv0 * c;
- double dz1dv_du = dSdv_du * c;
- double dz1dz_du = dSdz_du * c + S * dcdz_du;
-
- double dv_du1dv0 = ((dv1dv0 - dv01dv0) * (u1 - u01) + du01dv0 * (v1 - v01)) / (u1 - u01) / (u1 - u01);
- double dv_du1dv_du = dv1dv_du / (u1 - u01);
-
- double dz_du1dv0 = (dz1dv0 * (u1 - u01) + du01dv0 * (z1 - z01)) / (u1 - u01) / (u1 - u01);
- double dz_du1dz0 = (dz1dz0 - 1.) / (u1 - u01);
- double dz_du1dv_du = dz1dv_du / (u1 - u01);
- double dz_du1dz_du = dz1dz_du / (u1 - u01);
-
- //Set the track derivatives...
- deriv.setIdentity();
-
- deriv.set(IV, IV, dv1dv0);
- deriv.set(IV, IDVDU, dv1dv_du);
-
- deriv.set(IZ, IV, dz1dv0);
- deriv.set(IZ, IZ, dz1dz0);
- deriv.set(IZ, IDVDU, dz1dv_du);
- deriv.set(IZ, IDZDU, dz1dz_du);
-
- deriv.set(IDVDU, IV, dv_du1dv0);
- deriv.set(IDVDU, IDVDU, dv_du1dv_du);
-
- deriv.set(IDZDU, IV, dz_du1dv0);
- deriv.set(IDZDU, IZ, dz_du1dz0);
- deriv.set(IDZDU, IDVDU, dz_du1dv_du);
- deriv.set(IDZDU, IDZDU, dz_du1dz_du);
-
- deriv.set(IQP, IQP, 1.0);
-
- return pstat;
-
- }
- //cng
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/ShapeDispatcher.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/ShapeDispatcher.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,117 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-
-package org.hps.recon.tracking.kalman;
-
-import org.lcsim.detector.solids.Box;
-import org.lcsim.detector.solids.ISolid;
-import org.lcsim.detector.solids.Point3D;
-import org.lcsim.detector.solids.Trd;
-import org.lcsim.detector.solids.Tube;
-import org.lcsim.event.Track;
-import org.lcsim.fit.helicaltrack.HelicalTrackFit;
-import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
-import org.lcsim.recon.tracking.seedtracker.SeedTrack;
-import org.lcsim.recon.tracking.trffit.HTrack;
-
-/**
- * This class takes a solid and dispatches it to its Helper class,
- * for any of the possible Helper methods.
- *
- * @author ecfine
- */
-public class ShapeDispatcher implements ShapeHelper{
-
- public void printShape(ISolid solid) {
- if (solid instanceof Trd) {
- TrdHelper trapHelper = new TrdHelper();
- trapHelper.printShape(solid);
- } else if (solid instanceof Tube) {
-// TubeHelper tubeHelper = new TubeHelper();
-// tubeHelper.printShape(solid);
- } else if (solid instanceof Box) {
- BoxHelper boxHelper = new BoxHelper();
- boxHelper.printShape(solid);
- } else {
- System.out.println("Error! I don't recognize this shape!");
- }
- }
-
- public void printLocalCoords(ISolid solid) {
- if (solid instanceof Trd) {
- TrdHelper trapHelper = new TrdHelper();
- trapHelper.printLocalCoords(solid);
- } else if (solid instanceof Tube) {
-// TubeHelper tubeHelper = new TubeHelper();
-// tubeHelper.printLocalCoords(solid);
- } else if (solid instanceof Box) {
- BoxHelper boxHelper = new BoxHelper();
- boxHelper.printLocalCoords(solid);
- } else {
- System.out.println("Error! I don't recognize this shape!");
- }
- }
-
- public void printGlobalCoords(ISolid solid) {
- if (solid instanceof Trd) {
- TrdHelper trapHelper = new TrdHelper();
- trapHelper.printGlobalCoords(solid);
- } else if (solid instanceof Tube) {
-// TubeHelper tubeHelper = new TubeHelper();
-// tubeHelper.printGlobalCoords(solid);
- } else if (solid instanceof Box) {
- BoxHelper boxHelper = new BoxHelper();
- boxHelper.printGlobalCoords(solid);
- } else {
- System.out.println("Error! I don't recognize this shape!");
- }
- }
-
- public HTrack addIntercept(ISolid solid, HTrack ht) {
- HTrack track = null;
- if (solid instanceof Trd) {
- TrdHelper trapHelper = new TrdHelper();
-// track = trapHelper.addIntercept(solid, ht);
- } else {
- System.out.println("This shape is not supported yet");
-
- }
- return track;
- }
-
- public HelicalTrackFit trackToHelix(Track track) {
- SeedTrack stEle = (SeedTrack) track;
- SeedCandidate seedEle = stEle.getSeedCandidate();
- HelicalTrackFit ht = seedEle.getHelix();
- return ht;
- }
-
- public KalmanSurface getKalmanSurf(ISolid solid) {
- KalmanSurface surf = null;
- if (solid instanceof Trd) {
- TrdHelper trapHelper = new TrdHelper();
- surf = trapHelper.getKalmanSurf(solid);
- } else {
- System.out.println("This solid is not supported yet");
- }
- return surf;
- }
-
- public void findIntersection(ISolid solid, Track track) {
- throw new UnsupportedOperationException("Not supported yet.");
- }
-
- public boolean pointIsOnSolid(ISolid solid, Point3D hitPoint) {
- boolean pointOnSolid = false;
- if (solid instanceof Trd){
- TrdHelper trapHelper = new TrdHelper();
- pointOnSolid = trapHelper.pointIsOnSolid(solid, hitPoint);
- } else {
- System.out.print("This solid is not supported yet");
- }
- return pointOnSolid;
-
- }
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/ShapeHelper.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/ShapeHelper.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,35 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-
-package org.hps.recon.tracking.kalman;
-
-import org.lcsim.detector.solids.ISolid;
-import org.lcsim.event.Track;
-
-/**
- * Interface for dealing with lcsim solids and converting them to trf surfaces.
- * Most of the classes which extend this are basically empty, with the exception
- * of the TrdHelper, which has some geometry things, and ShapeDispatcher, which
- * dispatches solids to the correct helper.
- *
- * @author ecfine
- */
-public interface ShapeHelper {
-
- /* Prints the shape of the solid */
- public void printShape(ISolid solid);
-
- /* Prints the local coordinates of the solid */
- public void printLocalCoords(ISolid solid);
-
- /* Prints the global coordinates of the solid */
- public void printGlobalCoords(ISolid solid);
-
- /* Finds the intersection of a solid and a track */
- public void findIntersection(ISolid solid, Track track);
-
- public KalmanSurface getKalmanSurf(ISolid solid);
-
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/ThinXYPlaneMs.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/ThinXYPlaneMs.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,212 +0,0 @@
-package org.hps.recon.tracking.kalman;
-// ThinXYPlaneMs
-
-import org.lcsim.recon.tracking.trfbase.ETrack;
-import org.lcsim.recon.tracking.trfbase.Interactor;
-import org.lcsim.recon.tracking.trfbase.PropDir;
-import org.lcsim.recon.tracking.trfbase.Surface;
-import org.lcsim.recon.tracking.trfbase.TrackError;
-import org.lcsim.recon.tracking.trfbase.TrackVector;
-import org.lcsim.recon.tracking.trfutil.Assert;
-import org.lcsim.recon.tracking.trfutil.TRFMath;
-import org.lcsim.recon.tracking.trfxyp.SurfXYPlane;
-/**
- *
- * Class for modifying the covariance matrix of a track to account
- * for multiple scattering at a thin XY-plane whose material is
- *represented by the number of radiation lengths.
- *
- *@author Norman A. Graf
- *@version 1.0
- *
- */
-public class ThinXYPlaneMs extends Interactor
-{
-
- // static methods
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public static String typeName()
- { return "ThinXYPlaneMs"; }
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public static String staticType()
- { return typeName(); }
-
- // static attributes
- // Assign track parameter indices.
-
- private static final int IV = SurfXYPlane.IV;
- private static final int IZ = SurfXYPlane.IZ;
- private static final int IDVDU = SurfXYPlane.IDVDU;
- private static final int IDZDU = SurfXYPlane.IDZDU;
- private static final int IQP = SurfXYPlane.IQP;
-
- //attributes
- private double _radLength;
-
- //non-static methods
-
- // Constructor. The Interactor is constructed with the
- // appropriate number of radiation lengths.
-
- /**
- * Construct an instance from the number of radiation
- * lengths of the thin xy plane material.
- * The Interactor is constructed with the
- * appropriate number of radiation lengths.
- *
- * @param radLength The thickness of the material in radiation lengths.
- */
- public ThinXYPlaneMs(double radLength)
- {
- _radLength = radLength;
- }
-
- /**
- *Interact the given track in this xy plane,
- *using the thin material approximation for multiple scattering.
- *Note that the track parameters are not modified. Only the
- *covariance matrix is updated to reflect the uncertainty caused
- *by traversing the thin xy plane of material.
- *
- * @param tre The ETrack to scatter.
- */
- public void interact(ETrack tre)
- {
-
- // This can only be used with XYplanes... check that we have one..
-
- SurfXYPlane XYpl = new SurfXYPlane( 10.0, 0.0 );
- Surface srf = tre.surface();
- Assert.assertTrue( srf.pureType().equals(XYpl.pureType()) );
-
- TrackError cleanError = tre.error();
- TrackError newError = new TrackError(cleanError);
-
- // set the rms scattering appropriate to this momentum
-
- // Theta = (0.0136 GeV)*(z/p)*(sqrt(radLength))*(1+0.038*log(radLength))
-
-
-
- TrackVector theVec = tre.vector();
- double trackMomentum = theVec.get(IQP);
-
- double f = theVec.get(IDVDU);
- double g = theVec.get(IDZDU);
-
- double theta = Math.atan(Math.sqrt(f*f + g*g));
- double phi = f!=0. ? Math.atan(Math.sqrt((g*g)/(f*f))) : TRFMath.PI2;
- if ( f==0 && g<0 ) phi=3*TRFMath.PI2;
- if((f<0)&&(g>0))
- phi = Math.PI - phi;
- if((f<0)&&(g<0))
- phi = Math.PI + phi;
- if((f>0)&&(g<0))
- phi = 2*Math.PI - phi;
-
- double trueLength = _radLength/Math.cos(theta);
-
- double stdTheta = (0.0136)*trackMomentum*Math.sqrt(trueLength)*
- (1 + 0.038*Math.log(trueLength));
-
- double uhat = Math.sqrt(1-Math.sin(theta)*Math.sin(theta));
- double vhat = Math.sin(theta)*Math.cos(phi);
- double zhat = Math.sin(theta)*Math.sin(phi);
-
- double Qvu,Qzu,Qvz;
-
- // The MS covariance matrix can now be set.
-
- // **************** code for matrix ***********************//
-
- // Insert values for upper triangle... use symmetry to set lower.
-
- double norm = Math.sqrt(vhat*vhat + zhat*zhat);
-
- Qvu = (zhat/(norm*uhat))*(zhat/(norm*uhat));
- Qvu += Math.pow((vhat/norm)*(1 + (norm*norm)/(uhat*uhat)),2.0);
- Qvu *= stdTheta;
- Qvu *= stdTheta;
-
- Qzu = (vhat/(norm*uhat))*(vhat/(norm*uhat));
- Qzu += Math.pow((zhat/norm)*(1 + (norm*norm)/(uhat*uhat)),2.0);
- Qzu *= stdTheta;
- Qzu *= stdTheta;
-
- Qvz = - vhat*zhat/(uhat*uhat);
- Qvz += vhat*zhat/(norm*norm)*Math.pow((1 + (norm*norm)/(uhat*uhat)),2.0);
- Qvz *= stdTheta;
- Qvz *= stdTheta;
-
- newError.set(IDVDU,IDVDU, newError.get(IDVDU,IDVDU) + Qvu);
- newError.set(IDZDU,IDZDU, newError.get(IDZDU,IDZDU) + Qzu);
- newError.set(IDVDU,IDZDU, newError.get(IDVDU,IDZDU) + Qvz);
-
-
- tre.setError( newError );
-
- }
-
- // method for adding the interaction with direction
- public void interact_dir(ETrack tre, PropDir dir)
- {
- interact(tre);
- }
-
- /**
- *Return the number of radiation lengths.
- *
- * @return The thickness of the scattering material in radiation lengths.
- */
- public double radLength()
- {
- return _radLength;
- }
-
- /**
- *Make a clone of this object.
- *
- * @return A Clone of this instance.
- */
- public Interactor newCopy()
- {
- return new ThinXYPlaneMs(_radLength);
- }
-
-
- /**
- *Return a String representation of the class' type name.
- *Included for completeness with the C++ version.
- *
- * @return A String representation of the class' type name.
- */
- public String type()
- { return staticType(); }
-
- /**
- *output stream
- *
- * @return A String representation of this instance.
- */
- public String toString()
- {
- return "ThinXYPlaneMs with "+_radLength+" radiation lengths";
- }
-
-
-
-
-}
-
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/TrackUtils.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/TrackUtils.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,150 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-package org.hps.recon.tracking.kalman;
-
-import hep.physics.matrix.SymmetricMatrix;
-import hep.physics.vec.BasicHep3Vector;
-import hep.physics.vec.Hep3Vector;
-
-import org.lcsim.event.MCParticle;
-import org.lcsim.fit.helicaltrack.HelicalTrackFit;
-import org.lcsim.fit.helicaltrack.HelixParamCalculator;
-import org.lcsim.fit.helicaltrack.HelixUtils;
-import org.lcsim.recon.tracking.trfbase.TrackError;
-import org.lcsim.recon.tracking.trfbase.TrackSurfaceDirection;
-import org.lcsim.recon.tracking.trfbase.TrackVector;
-import org.lcsim.recon.tracking.trfbase.VTrack;
-import org.lcsim.recon.tracking.trfdca.SurfDCA;
-import org.lcsim.recon.tracking.trfutil.TRFMath;
-
-import Jama.Matrix;
-
-/**
- * Class for converting lcsim tracks to trf tracks.
- * @author ecfine
- */
-public class TrackUtils {
-
- private double mmTocm = 0.1; // TRF wants distances in cm, not mm.
- public double bz = 0.5;
- private double _epsilon = 1e-4;
- boolean _DEBUG = false;
-
- // Converts a HelicalTrackFit to a VTrack, prints track params.
- public VTrack makeVTrack(HelicalTrackFit t) {
-
- // Opposite sign convention.
- double r_signed = -(t.dca());
-
- TrackVector tv = new TrackVector();
- tv.set(0, r_signed * mmTocm);
- tv.set(1, t.z0() * mmTocm);
- tv.set(2, Math.tan(t.phi0()));
- tv.set(3, t.slope());
- tv.set(4, -t.curvature() / mmTocm / (TRFMath.BFAC * bz)); // curvature to q/pt
- double q = t.pT(bz) / (t.R() * mmTocm * TRFMath.BFAC * bz);
- if ((q / t.pT(bz) + tv.get(4)) > _epsilon) {
- System.out.println("something wrong with curvature? q/pt = "
- + (q / t.pT(bz)) + ", -tv.get(4) = " + -tv.get(4));
- }
-// SurfDCA s = new SurfDCA((t.dca() * Math.sin(t.phi0())),
-// (-t.dca() * Math.cos(t.phi0())));
- double dcaX = mmTocm * HelixUtils.PointOnHelix(t, 0).x();
- double dcaY = mmTocm * HelixUtils.PointOnHelix(t, 0).y();
- if(_DEBUG)System.out.println("DCA X="+dcaX+"; Y="+dcaY);
-// SurfDCA s = new SurfDCA(dcaX, dcaY);
-// I think this should be 0,0
- SurfDCA s = new SurfDCA(0, 0);
-
- //mg this is just wrong...
-// SurfDCA s = new SurfDCA(r_signed * mmTocm, 0.);
-
-// VTrack vt = new VTrack(s, tv);
- VTrack vt = new VTrack(s, tv, TrackSurfaceDirection.TSD_FORWARD);
- if (_DEBUG) {
- System.out.println("making VTrack with: ");
- System.out.println(" r_signed = " + r_signed * mmTocm);
- System.out.println(" z0 = " + t.z0() * mmTocm);
- System.out.println(" tanphi0 = " + Math.tan(t.phi0()));
- System.out.println(" tanlamda = " + t.slope());
- System.out.println(" q/pt = " + -t.curvature() / mmTocm / (TRFMath.BFAC * bz));
-// System.out.println(" q/pt = " + t.curvature() * 0.299999 * bz);
- System.out.println("from HelicalTrackFit with: ");
- System.out.println(" dca = " + t.dca() + "+/-" + Math.sqrt(t.covariance().e(HelicalTrackFit.dcaIndex, HelicalTrackFit.dcaIndex)));
- System.out.println(" z0 = " + t.z0() + "+/-" + Math.sqrt(t.covariance().e(HelicalTrackFit.z0Index, HelicalTrackFit.z0Index)));
- System.out.println(" phi0 = " + t.phi0() + "+/-" + Math.sqrt(t.covariance().e(HelicalTrackFit.phi0Index, HelicalTrackFit.phi0Index)));
- System.out.println(" slope = " + t.slope() + "+/-" + Math.sqrt(t.covariance().e(HelicalTrackFit.slopeIndex, HelicalTrackFit.slopeIndex)));
- System.out.println(" curvature = " + t.curvature() + "+/-" + Math.sqrt(t.covariance().e(HelicalTrackFit.curvatureIndex, HelicalTrackFit.curvatureIndex)));
-// System.out.println("SurfDCA at " + t.dca() * Math.sin(t.phi0()) + ", " +
-// -t.dca() * Math.cos(t.phi0()));
- }
- return vt;
- }
- //make a VTrack from an MC particle
-
- public VTrack makeVTrack(MCParticle mcp) {
- TrackVector tv = new TrackVector();
- HelixParamCalculator helix = new HelixParamCalculator(mcp, bz);
- double r_signed = -(helix.getDCA());
- double curvemc = (-1) / helix.getRadius();
- tv.set(0, r_signed * mmTocm);
- tv.set(1, helix.getZ0() * mmTocm);
- tv.set(2, Math.tan(helix.getPhi0()));
- tv.set(3, helix.getSlopeSZPlane());
- tv.set(4, curvemc / mmTocm / (TRFMath.BFAC * bz)); // curvature to q/pt
-// SurfDCA s = new SurfDCA(r_signed * mmTocm, 0.);
- SurfDCA s = new SurfDCA(0, 0.);
- VTrack vt = new VTrack(s, tv);
- return vt;
-
- }
-
- public void setBZ(double b) {
- bz = b;
- }
-
- public Hep3Vector getMomentum(VTrack vt){
- double qoverpt=vt.qOverP();
- double phi0=Math.atan(vt.vector().get(2));
- double slope=vt.vector().get(3);
- double Pt = Math.abs(1./qoverpt);
- double px = Pt * Math.cos(phi0);
- double py = Pt * Math.sin(phi0);
- double pz = Pt * slope;
- Hep3Vector p=new BasicHep3Vector(px,py,pz);
-
- return p;
- }
-
- // Returns a TrackError from a HelicalTrackFit.
- public TrackError getInitalError(HelicalTrackFit t) {
- SymmetricMatrix oldError = t.covariance();
- SymmetricMatrix newError = new SymmetricMatrix(5);
- double edca = oldError.e(HelicalTrackFit.dcaIndex, HelicalTrackFit.dcaIndex);
- newError.setElement(0, 0, edca * mmTocm * mmTocm * 1e6);
- double ez0 = oldError.e(HelicalTrackFit.z0Index, HelicalTrackFit.z0Index);
- newError.setElement(1, 1, ez0 * mmTocm * mmTocm * 1e6);
- double ephi0 = oldError.e(HelicalTrackFit.phi0Index, HelicalTrackFit.phi0Index);
- newError.setElement(2, 2, ephi0 * 1e6);
- double eslope = oldError.e(HelicalTrackFit.slopeIndex, HelicalTrackFit.slopeIndex);
- newError.setElement(3, 3, eslope * 100000);
- double ecurve = oldError.e(HelicalTrackFit.curvatureIndex, HelicalTrackFit.curvatureIndex);
- newError.setElement(4, 4, (ecurve / Math.pow(TRFMath.BFAC * bz * mmTocm, 2)) * 1e6);
-
-// for (int i = 0; i < error.getNColumns(); i++){
-// double newElement = error.e(i, i) * 10000;
-// error.setElement(i, i, newElement);
-// }
- Matrix errorMatrix = new Matrix(5, 5);
- for (int k = 0; k < 5; k++) {
- for (int j = 0; j < 5; j++) {
- errorMatrix.set(k, j, newError.e(k, j));
- }
- }
- if(_DEBUG) System.out.println("Setting initial error:\n"+errorMatrix.toString());
- TrackError trackerror = new TrackError(errorMatrix);
- return trackerror;
- }
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/TrdHelper.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/TrdHelper.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,281 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-
-package org.hps.recon.tracking.kalman;
-
-import hep.physics.vec.Hep3Vector;
-
-import java.util.ArrayList;
-import java.util.List;
-
-import org.lcsim.detector.IRotation3D;
-import org.lcsim.detector.ITransform3D;
-import org.lcsim.detector.solids.ISolid;
-import org.lcsim.detector.solids.Point3D;
-import org.lcsim.detector.solids.Polygon3D;
-import org.lcsim.detector.solids.Trd;
-import org.lcsim.event.Track;
-import org.lcsim.fit.helicaltrack.HelicalTrackFit;
-
-/**
- *
- * @author ecfine
- */
-public class TrdHelper implements ShapeHelper{
- PrintDetectorElements printDetectorElements = new PrintDetectorElements();
-
- /* Prints that solid is trapezoid */
- public void printShape(ISolid solid) {
- checkTrapezoid(solid);
- System.out.println("Shape: Trapezoid");
- }
-
- /* Prints the local coordinates of the trapezoid */
- public void printLocalCoords(ISolid solid){
- checkTrapezoid(solid);
- Trd trap = (Trd) solid;
- System.out.println("Local Coordinates: " + trap.getVertices().get(0));
- for (int n = 1; n < trap.getVertices().size(); n ++) {
- printDetectorElements.printIndent();
- System.out.println(" " + trap.getVertices().get(n));
- }
- }
-
- /* Prints the global coordinates of the trapezoid */
- public void printGlobalCoords(ISolid solid){
- checkTrapezoid(solid);
- Trd trap = (Trd) solid;
- for(int i = 0; i < trap.getVertices().size(); i++) {
- Hep3Vector currentPoint = (Hep3Vector) trap.getVertices().get(i);
- Hep3Vector transformedPoint = (Hep3Vector) currentPoint;
- for(int k = PrintDetectorElements.physicalVolumes.size(); k > 0; k--) {
- ITransform3D lastTransform = (ITransform3D) PrintDetectorElements.physicalVolumes.get(k-1);
- currentPoint = lastTransform.transformed(currentPoint);
- }
- transformedPoint = currentPoint;
- printDetectorElements.printIndent();
- if (i == 0){
- System.out.println("Global Coordinates: [ " + transformedPoint.x() +
- ", " + transformedPoint.y() + ", " + transformedPoint.z() +"]");
- } else{
- System.out.println(" [ " + transformedPoint.x() +
- ", " + transformedPoint.y() + ", " + transformedPoint.z() +"]");
- }
- }
- }
-
- public void findIntersection(ISolid solid, Track track) {
- checkTrapezoid(solid);
- Trd trap = (Trd) solid;
- ShapeDispatcher shapedis = new ShapeDispatcher();
- HelicalTrackFit helix = shapedis.trackToHelix(track);
- List<Polygon3D> faces = trap.getFaces();
-
- // For each face, find the vertices
- for (int i = 0; i < faces.size(); i++) {
- Polygon3D face = faces.get(i);
- Hep3Vector localNormal = face.getNormal();
- Hep3Vector normal = getGlobalRot(localNormal); // get the normal in global coords
- List<Hep3Vector> vertices = new ArrayList<Hep3Vector>();
-
- // For each vertex, find and print the Global coordinates
- for(int k = 0; k < face.getVertices().size(); k++){
- Hep3Vector localVertex = (Hep3Vector) face.getVertices().get(k);
-// System.out.println("local vertex: " + localVertex);
- Hep3Vector vertex = getGlobalCoords(localVertex);
- System.out.println("vertex: " + vertex);
- vertices.add(vertex);
- }
-// System.out.println("local normal: " + localNormal);
- System.out.println("normal: " + normal);
- if (0.01 > normal.z() && normal.z() > -0.01){
- if (org.lcsim.fit.helicaltrack.HelixUtils.isInterceptingBoundedXYPlane(helix, normal, vertices)) {
- System.out.println("INTERSECTING BOUNDED PLANE!");
- } else {
- System.out.println("No intersection");
- }
- if (org.lcsim.fit.helicaltrack.HelixUtils.isInterceptingXYPlane(helix, normal, vertices.get(0))) {
- System.out.println("INTERSECTING PLANE!");
- } else {
- System.out.println("No intersection");
- }
- } else {
- System.out.println("wrong orientation");
- }
- }
- }
-
-
- // This is (clearly) a work in progress...
- public KalmanSurface getKalmanSurf(ISolid solid) {
- KalmanSurface surf = null;
-// checkTrapezoid(solid);
-// Trd trap = (Trd) solid;
-// List<Polygon3D> faces = trap.getFaces();
-// /* This needs to decide what surface the trapezoid should be modeled as.
-// * I suspect xy plane and z plane are the main choices. For now, it's
-// * just going to default to an xy plane, but that should be revisited. */
-//// List vertices = trap.getVertices();
-//
-// // From a list of 8 vertices, need to match the pairs with the smallest difference
-// // between them. Or alternatively, need to find the two biggest squares. That
-// // may actually be a better way to do things... could do it using faces that way?
-//
-// ArrayList areas = new ArrayList();
-// for (int i = 0; i < faces.size(); i++) {
-// Polygon3D face = faces.get(i);
-// double area = findArea(face);
-// areas.add(area);
-// }
-// for (int k )
-//
-//
-// Hep3Vector localNormal = face.getNormal();
-// Hep3Vector normal = getGlobalRot(localNormal); // get the normal in global coords
-// List<Hep3Vector> vertices = new ArrayList<Hep3Vector>();
-// Hep3Vector planeNormal = null;
-// // Only want to make surface for xy plane faces
-// if (0.01 > normal.z() && normal.z() > -0.01){
-// // For each vertex, find the Global coordinates & add to vertices
-// for(int k = 0; k < face.getVertices().size(); k++){
-// Hep3Vector localVertex = (Hep3Vector) face.getVertices().get(k);
-// Hep3Vector vertex = getGlobalCoords(localVertex);
-// vertices.add(vertex);
-// planeNormal = normal;
-// }
-// /* From the 8 vertices, we want one origin point, since
-// * the origin and a vector define the plane. Need a single x & y value.
-// * Here, we find the smallest difference in x and y values between two
-// * points. Once we find this pair, we can average the x and y values
-// * for our origin point. */
-// double minxDiff = (vertices.get(0).x() - vertices.get(1).x());
-// double minyDiff = (vertices.get(0).y() - vertices.get(1).y());
-// int secondPoint = 1;
-// for (int j = 2; j < vertices.size(); j++){
-// double xDiff = vertices.get(0).x() - vertices.get(j).x();
-// double yDiff = vertices.get(0).y() - vertices.get(j).y();
-// if (xDiff < minxDiff && yDiff < minyDiff) {
-// minxDiff = xDiff;
-// minyDiff = yDiff;
-// secondPoint = j;
-// }
-// }
-// // Average the x, y, and z values of the two points to find origin.
-// Hep3Vector doubleMidPoint = VecOp.add(vertices.get(0), vertices.get(secondPoint));
-// Hep3Vector origin = VecOp.mult(.5, doubleMidPoint);
-//
-// /* Make a new surface from the origin and plane normal as found
-// * above. */
-//
-// double distnorm = (planeNormal.x() * origin.x() + planeNormal.y() * origin.y()
-// + planeNormal.z() * origin.z())/planeNormal.magnitude();
-// double normphi = Math.atan2(planeNormal.y(), planeNormal.x());
-// if (normphi < 0){
-// normphi = normphi + TRFMath.TWOPI;
-// }
-// System.out.println("Plane normal = " + planeNormal + ", origin = " + origin);
-// surf = new KalmanSurface("xyplane", distnorm, normphi);
-// }
-//
-//
-//
-
- return surf;
- }
-
- boolean pointIsOnSolid(ISolid solid, Point3D hitPoint) {
- boolean onSolid = true;
- checkTrapezoid(solid);
- Trd trap = (Trd) solid;
- List faces = trap.getFaces();
- // Want to find all projections for each face, and see if projections
- // of the hit point are inside them. This is almost certainly not the
- // most efficient way to see whether the point is in the solid, but
- // I'm not really sure how else to deal with it.
- for(int i = 0; i < faces.size(); i++){
- Polygon3D face = (Polygon3D) faces.get(i);
- List vertices = face.getVertices();
- ArrayList xProjVertices = new ArrayList();
- ArrayList yProjVertices = new ArrayList();
- ArrayList zProjVertices = new ArrayList();
- for (int k = 0; k < vertices.size(); k ++){
- Point3D currentPoint = (Point3D) vertices.get(k);
- double[] xProj = new double[2];
-
-
- }
- }
- return onSolid;
- }
-
- // Finds the area of an arbitrary quadrilateral. See http://softsurfer.com/
- // Archive/algorithm_0101/algorithm_0101.htm#Quandrilaterals for more info.
- private double findArea(Polygon3D face) {
- double[] v2_v0 = new double[3];
- v2_v0[0] = face.getVertices().get(2).x() - face.getVertices().get(0).x();
- v2_v0[1] = face.getVertices().get(2).y() - face.getVertices().get(0).y();
- v2_v0[2] = face.getVertices().get(2).z() - face.getVertices().get(0).z();
- double[] v3_v1 = new double[3];
- v3_v1[0] = face.getVertices().get(3).x() - face.getVertices().get(1).x();
- v3_v1[1] = face.getVertices().get(3).y() - face.getVertices().get(1).y();
- v3_v1[2] = face.getVertices().get(3).z() - face.getVertices().get(1).z();
- double area = .5 * magnitude(crossProduct(v2_v0, v3_v1));
- return area;
- }
-
- // Finds the cross product of two 3d vectors, represented as double[]s.
- private double[] crossProduct(double[] v1, double[] v2){
- double[] cross = new double[3];
- cross[0] = v1[1]*v2[2] - v2[1]*v1[2];
- cross[1] = v1[2]*v2[0] - v2[2]*v1[0];
- cross[2] = v1[0]*v2[1] - v2[0]*v1[1];
- return cross;
- }
-
- // Finds the magnitude of a 3d vector, represented as a double[]
- private double magnitude(double[] v){
- double mag = Math.sqrt(Math.pow(v[0], 2) + Math.pow(v[1], 2) + Math.pow(v[2], 2));
- return mag;
- }
-//
-// public HTrack addIntercept(ISolid solid, HTrack ht) {
-//
-// }
-
-
-
-
- /* Finds the global coordinates of a local point */
- private Hep3Vector getGlobalCoords(Hep3Vector local){
- Hep3Vector currentPoint = local;
- for(int k = PrintIntersections.physicalVolumes.size(); k > 0; k--) {
- ITransform3D lastTransform = (ITransform3D) PrintIntersections.physicalVolumes.get(k-1);
- currentPoint = lastTransform.transformed(currentPoint);
- }
- return currentPoint;
- }
-
- /* Finds the global rotation of a local point. Use this for normals, where
- the translation is irrelevant. */
- private Hep3Vector getGlobalRot(Hep3Vector local){
- Hep3Vector currentPoint = local;
- for(int k = PrintIntersections.physicalVolumes.size(); k > 0; k--) {
- ITransform3D lastTransform = (ITransform3D) PrintIntersections.physicalVolumes.get(k-1);
- IRotation3D lastRotation = lastTransform.getRotation();
- currentPoint = lastRotation.rotated(currentPoint);
- }
- return currentPoint;
- }
-
- /* Checks that solid is really trapezoid */
- public void checkTrapezoid(ISolid solid){
- if (solid instanceof Trd) {
- } else {
- System.out.println("Error! This shape is not a trapezoid!");
- return;
- }
- }
-
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/TrfDetector.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/TrfDetector.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,44 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-
-package org.hps.recon.tracking.kalman;
-
-import java.util.List;
-import java.util.Map;
-
-import org.lcsim.recon.tracking.trflayer.Detector;
-import org.lcsim.recon.tracking.trflayer.Layer;
-
-/**
- *
- * @author ecfine
- */
-public class TrfDetector extends Detector {
- // List of layer names in original order.
- private List _names;
-
- // Map layer names to layers.
- private Map _layermap;
-
-
-
- public int addLayer(String name, Layer lyr)
- {
- // Set unchecked.
-// if ( _check == CHECKED_OK ) _check = UNCHECKED;
-
- // Check name does not already appear in map.
- if ( isAssigned(name) ){
- System.out.println("Name already assigned, layer not added.");
- return 0;
- }
- // Add layer.
- _names.add(name);
- _layermap.put(name, lyr);
- System.out.println("Layer named " + name + " added");
- return 1;
- }
-
-}
java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/TubeHelper.java 2014-03-27 03:22:56 UTC (rev 397)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/kalman/TubeHelper.java 2014-03-27 18:50:25 UTC (rev 398)
@@ -1,44 +0,0 @@
-/*
- * To change this template, choose Tools | Templates
- * and open the template in the editor.
- */
-package org.hps.recon.tracking.kalman;
-
-import org.lcsim.detector.solids.ISolid;
-import org.lcsim.detector.solids.Tube;
-import org.lcsim.event.Track;
-
-/**
- *
- * @author ecfine
- */
-public class TubeHelper implements ShapeHelper {
-
- public void printShape(ISolid solid) {
- checkTube(solid);
- System.out.println("Shape: Tube");
-
- }
-
- public void printLocalCoords(ISolid solid) {
- }
-
- public void printGlobalCoords(ISolid solid) {
- }
-
- public void checkTube(ISolid solid) {
- if (solid instanceof Tube) {
- } else {
- System.out.print("Error! This isn't a tube!");
- return;
- }
- }
-
- public void findIntersection(ISolid solid, Track track) {
- }
-
- public KalmanSurface getKalmanSurf(ISolid solid) {
- KalmanSurface surf = null;
- return surf;
- }
-}
SVNspam 0.1