10 removed files
lcsim/src/org/lcsim/contrib/KFFitter
diff -N History
--- History 13 Apr 2007 20:32:06 -0000 1.4
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,84 +0,0 @@
-## History file
-## Kalman Filter Fitter (KFFitter)
-## Please summarize changes to the KFFitter code here.
-## Most recent first please.
-
-13 April 2007 Fred Blanc V01-01-02
-
- - Added KFVertex class to handle a vertex constraint.
- The constraint is applied to d0=0 and z0=0, with an uncertainty
- of +-2 microns
- The constraint it applied at the end of the outside-in fit, in the KFTrack class.
-
- - Account for the uncertainty on the dE/dx correction, in the 'filter' function
- of KFSite
-
- - removed from package the unused KFMatrix class
-
- - cleaned up several obsolete debugging statements and a few debugging functions
-
-
-09 March 2007 Fred Blanc V01-01-01
-
- (this version of the code was used to produce the results presented at the
- SiD tracking meeting on the day of this commit)
-
- Many new feature:
-
- - KFFitterDriver.java
- add fit results to event via the BaseTrack class
- - KFTrack.java
- use again the track finder seed track as intial parameters
- refined the initial covariance matrix
- - KFSite.java
- refined the multiple scattering and energy loss corrections
- - KFPoint.java
- added the gaussian smearing of SimTrackerHits:
- 5micron in rphi and z for the Vertex Detector
- 7micron in rphi and "large" for z in the Tracker
- - KFTrackParameters.java
- fixed bug in conversion between "fit" and "point" parametrizations (in _alpha and _z0)
- turned off the calls to the helix parametrization, as it is not needed (but the code is still there)
- - test/KFTest.java
- created several diagnostic plots of the fit variables, and their pulls
-
- - in all files:
- turn on/off all the log messages with a flag in the KFConstant class
- made "muon" the default hypothesis
-
-01 February 2007 Fred Blanc V01-01-00
- - using the track fit parameters (d0,phi0,omega,z0,s) at all stages of the
- propagation and filtering.
- - added effects of:
- 1) multiple scattering: add MS error matrix to track covariance matrix before filtering
- 2) energy loss: subtract or add dE/dx at the end of the filtering step
- - calculate the seed track using the first 3 hits on a track
-
-
-13 December 2006 Fred Blanc V01-00-00
- - Create KFFitter package. Included files:
- KFConstants.java stores a list of constants
- KFFitterDriver.java main driver
- KFMain.java used to run on the command line
- java org.lcsim.contrib.KFFitter.KFMain <data file> <number of events to run>
- KFMatrix.java unused (using Jama instead)
- KFPoint.java defines a point in space and its error matrix
- KFSite.java defines a Kalman site and handles the filtering of propagates track
- KFSites can be a measurement and/or a material site
- KFTrack.java contains a list of KFSites.
- handles the swimming trough list of KFSites
- KFTrackParameters.java defines the track parametrization and its error matrix
- Supports three parameterizations:
- 1) fit parameters: d0, phi0, omega, z0, tanl (+ path length)
- 2) point parameters: x, y, z, px, py, pz
- 3) helix parameters: xc, yc, rc, z0, tanl
-
- - KFFitterDriver.java is the main driver
- - Track parameterization is transformed to the r,phi,z,px,py,pz for
- the filtering stage (this will changes in a coming version)
- - Example code is given in test/TestKFF.java
- - The default functionality is: 1) loop over tracks reconstructed
- with the SODTrackFinder package, 2) apply the Kalman Filter
- at sites defined as measurement points, or material points
- (material points are not yet used in the first version of the
- code).
lcsim/src/org/lcsim/contrib/KFFitter
diff -N KFConstants.java
--- KFConstants.java 13 Apr 2007 20:32:06 -0000 1.4
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,40 +0,0 @@
-package org.lcsim.contrib.KFFitter;
-
-/**
- * Definition of constants in KFFitter package
- * @author Fred Blanc and Steve Wagner
- * @version $Id: KFConstants.java,v 1.4 2007/04/13 20:32:06 fblanc Exp $
- */
-
-public class KFConstants
-{
-
- public KFConstants()
- {
- }
-
- public double Pi() {return pi;}
- public double twoPi() {return twopi;}
- public double vtxResolution() {return vtxres;}
- public double sodResolution() {return sodres;}
- public double pointThreeTimesB() {return ptconst;}
- public double pionMass() {return piMass;}
- public double muonMass() {return muMass;}
- public double electronMass() {return eMass;}
-
- public boolean debug() {return DeBug;}
-
- private static double pi = 3.14159265359;
- private static double twopi = 6.283185307;
- private static double vtxres = 0.0005;
- private static double sodres = 0.0007;
- private static double ptconst = 0.0149896229; // = 0.299792458*B
-
- private static double piMass = 0.13957018;
- private static double muMass = 0.105658369;
- private static double eMass = 0.000510998;
-
- private static boolean DeBug = false;
-
-}
-
lcsim/src/org/lcsim/contrib/KFFitter
diff -N KFFitterDriver.java
--- KFFitterDriver.java 8 Dec 2007 02:39:10 -0000 1.6
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,200 +0,0 @@
-package org.lcsim.contrib.KFFitter;
-
-import java.util.*;
-import java.io.*;
-import Jama.Matrix;
-import Jama.util.Maths;
-import hep.physics.matrix.SymmetricMatrix;
-import org.lcsim.event.EventHeader;
-import org.lcsim.util.Driver;
-import org.lcsim.event.EventHeader.LCMetaData;
-import org.lcsim.event.TrackerHit;
-import org.lcsim.event.SimTrackerHit;
-import org.lcsim.event.MCParticle;
-import org.lcsim.event.base.BaseTrack;
-import org.lcsim.util.aida.AIDA;
-import hep.aida.*;
-
-import org.lcsim.contrib.SODTracker.SODTrackFinderDriver;
-import org.lcsim.contrib.SODTracker.SODTrack;
-
-import org.lcsim.contrib.KFFitter.KFTrack;
-import org.lcsim.contrib.KFFitter.KFConstants;
-
-/**
- * Driver for KFFitter package
- * @author Fred Blanc and Steve Wagner
- * @version $Id: KFFitterDriver.java,v 1.6 2007/12/08 02:39:10 jeremy Exp $
- */
-
-public class KFFitterDriver extends Driver
-{
- private AIDA aida = AIDA.defaultInstance();
- private java.util.Random gtor;
-
- public KFFitterDriver()
- {
- add(new SODTrackFinderDriver());
- gtor = new java.util.Random(10101);
-
- }
-
- public void process(EventHeader event)
- {
- super.process(event); // this takes care that the child Drivers are loaded and processed.
-
- List<BaseTrack> listTrk = new ArrayList<BaseTrack>();
-
- // 1. Extract list of SODTracks from event
- List<SODTrack> SODTrackList = (List<SODTrack>) event.get("SODTracks");
-
- // 2. Loop over SODTracks and fit them
- for ( SODTrack SODtrk : SODTrackList ) {
-
- double trkPT = SODtrk.getPT();
- aida.cloud1D("SODPt").fill(trkPT);
-
- // create new KFTrack and initialize it to the SODTrack fit parameters
- KFTrack track = new KFTrack();
- track.defineSeedTrack(SODtrk);
-
- List<SimTrackerHit> hits = (List<SimTrackerHit>) SODtrk.getHits();
- int nbHits = hits.size();
-
- double chi2 = SODtrk.getChi2();
-
- double ndf = nbHits;
-
- double genPtot = 0.;
- double genPt = 0.;
- double genD0 = 0.;
- double genZ0 = 0.;
- double genTanl = 0.;
-
- if(chi2<100.) { // run fitter only if good track seed
- for ( SimTrackerHit hit : hits ) {
- // assign hits to KFTrack
- track.addHit(hit,gtor.nextGaussian(),gtor.nextGaussian());
-
- genPtot = hit.getMCParticle().getMomentum().magnitude();
- double genPx = hit.getMCParticle().getPX();
- double genPy = hit.getMCParticle().getPY();
- double genPz = hit.getMCParticle().getPZ();
- genPt = Math.sqrt(genPx*genPx+genPy*genPy);
- if(genPt>0) {genTanl = genPz/genPt;}
- double genX0 = hit.getMCParticle().getOriginX();
- double genY0 = hit.getMCParticle().getOriginY();
- genD0 = Math.sqrt(genX0*genX0+genY0+genY0);
- genZ0 = hit.getMCParticle().getOriginZ();
- }
-
- // *************************
- // run Kalman Filter Fit
- // *************************
- track.fit();
-
- // *****************************
- // add track in event record
- // *****************************
- BaseTrack fittedTrack = new BaseTrack();
-
- KFTrackParameters kalmanTrack = track.getKFTrackParameters();
- kalmanTrack.setPointToPOCA();
-
- double[] pnt = {0.,0.,0.};
- pnt[0] = kalmanTrack.getPointParameters()[0];
- pnt[1] = kalmanTrack.getPointParameters()[1];
- pnt[2] = kalmanTrack.getPointParameters()[2];
-
- fittedTrack.setReferencePoint(pnt);
- fittedTrack.setRefPointIsDCA(true);
- fittedTrack.setFitSuccess(true);
-
- double[] fitPar = kalmanTrack.getFitParameters();
- // unit convertion from cm -> mm
- fitPar[0] *= 10.;
- fitPar[2] *= -0.1;
- fitPar[3] *= 10.;
-
- double magfield = 5.0;
- fittedTrack.setTrackParameters(fitPar, magfield);
-
- Matrix covMat = kalmanTrack.getErrorMatrix();
-
- // setup unit-conversion matrix (cm -> mm)
- double[][] diag = new double[5][5];
- diag[0][0] = 10.; // d0
- diag[1][1] = 1.; // phi0
- diag[2][2] = -0.1; // omega
- diag[3][3] = 10.; // z0
- diag[4][4] = 1.; // tanl
- Matrix convMat = new Matrix(diag);
-
- Matrix covMatmm = convMat.times(covMat).times(convMat);
-
- SymmetricMatrix errMat = new SymmetricMatrix(Maths.fromJamaMatrix(covMatmm));
- fittedTrack.setCovarianceMatrix(errMat);
-
- double KFchi2 = 1.; // fix me
- int KFndf = 1; // fix me
- fittedTrack.setChisq(KFchi2);
- fittedTrack.setNDF(KFndf);
-
- if(Math.abs(fitPar[0])<100.) {
- boolean fillCould = false;
- int nBins = 25;
- double min = -5.;
- double max = 5.;
- listTrk.add(fittedTrack);
-
- // fill pull plots for d0, omega, z0, tanl, and pt
- KFConstants _cst = new KFConstants();
-
- double pullD0 = (fitPar[0]-genD0)/Math.sqrt(errMat.e(0,0));
- if(fillCould) aida.cloud1D("KFD0Pull").fill(pullD0);
- aida.histogram1D("KFD0PullH",nBins,min,max).fill(pullD0);
-
- double genOmega = 0.1*_cst.pointThreeTimesB()/genPt;
- double pullOmega = (Math.abs(fitPar[2])-genOmega)/Math.sqrt(errMat.e(2,2));
- if(fillCould) aida.cloud1D("KFOmegaPull").fill(pullOmega);
- aida.histogram1D("KFOmegaPullH",nBins,min,max).fill(pullOmega);
-
- double fitPt = 0.1*_cst.pointThreeTimesB()/Math.abs(fitPar[2]);
- double errPt = Math.sqrt(errMat.e(2,2))*0.1*_cst.pointThreeTimesB()/(fitPar[2]*fitPar[2]);
- double sSquared = fitPar[4]*fitPar[4];
- double fitPtot = fitPt*Math.sqrt(1.+sSquared);
- double errPtot = fitPt*
- Math.sqrt(((1.+sSquared)/(fitPar[2]*fitPar[2]))*errMat.e(2,2)+
- (sSquared/(1.+sSquared))*errMat.e(4,4));
-
- aida.cloud1D("KFerrPt").fill(errPt);
- aida.cloud1D("KFerrPtot").fill(errPtot);
-
- double pullPt = (fitPt-genPt)/errPt;
- if(fillCould) aida.cloud1D("KFPtPull").fill(pullPt);
- aida.histogram1D("KFPtPullH",nBins,min,max).fill(pullPt);
- double pullPtot = (fitPtot-genPtot)/errPtot;
- if(fillCould) aida.cloud1D("KFPtotPull").fill(pullPtot);
- aida.histogram1D("KFPtotPullH",nBins,min,max).fill(pullPtot);
- //System.out.println("Driver: pt and ptot = "+genPt+" "+genPtot);
-
- double pullZ0 = (fitPar[3]-genZ0)/Math.sqrt(errMat.e(3,3));
- if(fillCould) aida.cloud1D("KFZ0Pull").fill(pullZ0);
- aida.histogram1D("KFZ0PullH",nBins,min,max).fill(pullZ0);
- double pullTanl = (fitPar[4]-genTanl)/Math.sqrt(errMat.e(4,4));
- //double pullTanl = (Math.abs(fitPar[4])-genTanl)/Math.sqrt(errMat.e(4,4));
- if(fillCould) aida.cloud1D("KFTanlPull").fill(pullTanl);
- aida.histogram1D("KFTanlPullH",nBins,min,max).fill(pullTanl);
- }
-
- } else {
- System.out.println("Driver: skip this track");
- System.out.println();
- } // if(chi2<100.) {
- } // for ( SODTrack SODtrk : SODTrackList ) {
-
- // write list of Kalman-fitted tracks in the event
- event.put("KFTracks", listTrk, BaseTrack.class, 0);
- }
-
-}
lcsim/src/org/lcsim/contrib/KFFitter
diff -N KFMain.java
--- KFMain.java 11 Sep 2007 00:21:03 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,102 +0,0 @@
-package org.lcsim.contrib.KFFitter;
-
-import java.util.*;
-import org.lcsim.util.Driver;
-import org.lcsim.event.*;
-
-import org.lcsim.util.loop.*;
-import java.io.File;
-
-/**
- * Main function used in command-line running of KFFitter package
- * @author Fred Blanc and Steve Wagner
- * @version $Id: KFMain.java,v 1.3 2007/09/11 00:21:03 tonyj Exp $
- */
-
-public class KFMain extends Driver
-{
- private int nEvents = 0;
-
- public KFMain() throws Exception
- {
- this(new String[0]);
- }
-
- public KFMain(String[] args) throws Exception
- {
- add(new KFFitterDriver());
- }
-
- public void process(EventHeader event)
- {
- super.process(event);
- if (nEvents%10==0) {
- System.out.println("done processing event "+nEvents);
- }
- ++nEvents;
- }
-
- public void mainCopy(String[] args) throws Exception
- {
- if (args.length >= 1 &&
- args[0].matches("^\\s*[/\\.\\w-_\\+]+\\.slcio")) {
- KFMain driver = this;
- convert(args[0],driver);
- }
- else {
- System.out.println("usage: java KFMain <input.slcio> <nb events>");
- return;
- }
- }
-
- public static void main(String[] args) throws Exception
- {
- if (args.length >= 1 &&
- args[0].matches("^\\s*[/\\.\\w-_\\+]+\\.slcio")) {
- KFMain driver = new KFMain(args);
- int nbEventsToProcess = -1;
- if (args.length == 2 && args[1].matches("[0-9]+")) { // matches a number at least 1 time
- nbEventsToProcess = Integer.parseInt(args[1]);
- System.out.print("Processing "); System.out.print(nbEventsToProcess); System.out.println(" events");
- }
- convert(args[0],driver,nbEventsToProcess);
- }
- else {
- System.out.println("usage: java KFMain <input.slcio> <nb events>");
- return;
- }
- }
-
- private static void convert(String iFile) throws Exception
- {
- convert(iFile, new KFMain());
- }
-
- private static void convert(String iFile, KFMain driver) throws Exception
- {
- int nLoop = -1;
- convert(iFile, new KFMain(), nLoop);
- }
-
- private static void convert(String iFile, KFMain driver, int nbLoops) throws Exception
- {
- // Setup LCSim event loop
- LCSimLoop loop = new LCSimLoop();
-
- File input = new File(iFile);
-
- loop.setLCIORecordSource(input);
- loop.add(driver);
- // loop.add(new LCIODriver(output));
-
- loop.loop(nbLoops, null);
- loop.dispose();
-
- System.out.println("Processed "+driver.getNEvents()+" events.");
- }
-
- public int getNEvents()
- {
- return nEvents;
- }
-}
lcsim/src/org/lcsim/contrib/KFFitter
diff -N KFPoint.java
--- KFPoint.java 13 Apr 2007 20:32:06 -0000 1.5
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,98 +0,0 @@
-package org.lcsim.contrib.KFFitter;
-
-import Jama.Matrix;
-
-/**
- * Define a point in the KFFitter package
- * Used for holding 3D coordinates and error matrix
- * @author Fred Blanc and Steve Wagner
- * @version $Id: KFPoint.java,v 1.5 2007/04/13 20:32:06 fblanc Exp $
- */
-
-public class KFPoint
-{
-
- public KFPoint()
- {
- _x = 0.;
- _y = 0.;
- _z = 0.;
- _r = 0.;
- _phi = 0.;
-
- _rphiErr = 0.;
- _zErr = 0.;
-
- }
-
- public void makePoint(double X, double Y, double Z)
- {
- _x = X;
- _y = Y;
- _z = Z;
-
- _r = Math.sqrt(_x*_x+_y*_y);
- _phi = Math.atan2(_y,_x);
-
- _rphiErr = 0.0007;
- if(_r<10.){_rphiErr = 0.0005;};
-
- _zErr = 0.0005;
- if(_r>10.){_zErr = 1000.;}
-
- }
-
- public void smearRPhiPoint(double rdm1, double rdm2)
- {
- _phi = _phi+(_rphiErr/_r)*rdm1;
- _x = _r*Math.cos(_phi);
- _y = _r*Math.sin(_phi);
-
- if(_zErr<0.01){
- _z = _z + _zErr * rdm2;
- }
- }
-
- public double getX() { return _x; }
- public double getY() { return _y; }
- public double getZ() { return _z; }
- public double getR() { return _r; }
- public double getPhi() { return _phi; }
- public double getDistanceFromOrigin() {return Math.sqrt(_r*_r+_z*_z);}
-
- public double getRPhiErr() { return _rphiErr;}
- public double getZErr() { return _zErr;}
-
- public Matrix getPhiZErrorMatrix()
- {
- double[][] vals = {{_rphiErr*_rphiErr/(_r*_r),0.},{0.,_zErr*_zErr}};
- Matrix errMat = new Matrix(vals);
- return errMat;
- }
-
- public Matrix getXYZErrorMatrix()
- {
- double s2 = _rphiErr*_rphiErr;
-
- double xx = Math.sin(_phi)*Math.sin(_phi)*s2;
- double xy = -Math.sin(_phi)*Math.cos(_phi)*s2;
- double yx = xy;
- double yy = Math.cos(_phi)*Math.cos(_phi)*s2;
-
- double zz = _zErr*_zErr;
-
- double[][] vals = {{xx,xy,0.},{yx,yy,0.},{0.,0.,zz}};
- Matrix errMat = new Matrix(vals);
-
- return errMat;
- }
-
- private double _x;
- private double _y;
- private double _z;
-
- private double _r;
- private double _phi;
- private double _rphiErr;
- private double _zErr;
-}
lcsim/src/org/lcsim/contrib/KFFitter
diff -N KFSite.java
--- KFSite.java 13 Apr 2007 20:32:06 -0000 1.5
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,460 +0,0 @@
-package org.lcsim.contrib.KFFitter;
-
-import org.lcsim.event.SimTrackerHit;
-import org.lcsim.event.TrackerHit;
-
-import java.util.ArrayList;
-import java.util.List;
-
-import Jama.Matrix;
-
-/**
- * Define a site used in the Kalman filter fit
- * A site can be any combination of a measurement, material location
- * or a B-field discontinuity. The class handles the filtering step
- * of the Kalman filter.
- * @author Fred Blanc and Steve Wagner
- * @version $Id: KFSite.java,v 1.5 2007/04/13 20:32:06 fblanc Exp $
- */
-
-public class KFSite
-{
-
- public KFSite()
- {
- }
-
- public void setLocation(SimTrackerHit hit, double rdm1, double rdm2)
- {
- double[] pos = hit.getPoint();
- // convert to cm
- for(int i=0; i<3; i++) {pos[i]*=0.1;}
- _siteLocation.makePoint(pos[0],pos[1],pos[2]);
- _siteLocation.smearRPhiPoint(rdm1, rdm2);
-
- if(_cst.debug()) System.out.println("KFSite: Add Pt at coord = "+pos[0]+" "+pos[1]+" "+pos[2]);
-
- _PhiZErrorMatrix = _siteLocation.getPhiZErrorMatrix();
-
- _PhiZErrorMatrixInverse = _PhiZErrorMatrix.inverse();
-
- if(_cst.debug()) {
- System.out.println("KFSite: PhiZ cov matrix:");
- _PhiZErrorMatrix.print(12,10);
- }
- }
-
- public KFTrackParameters filter(KFTrackParameters trk)
- {
-
- double[] tpp = trk.getFitParameters();
- Matrix stateVect = new Matrix(tpp,5);
- if(_cst.debug()) {
- System.out.println("KFSite: initial stateVect:");
- stateVect.print(15,12);
- System.out.println("KFSite: initial stateVect in point parameter space:");
- trk.printPointParams();
- }
-
- double d0 = tpp[0];
- double phi0 = tpp[1];
- double omega = tpp[2];
- double z0 = tpp[3];
- double tanl = tpp[4];
- double l = 0.;
-
- Matrix trkCovNoNoise = trk.getErrorMatrix();
- if(_cst.debug()) {
- System.out.println("KFSite: trkCovNoNoise:");
- trkCovNoNoise.print(15,12);
- }
-
- Matrix Q = new Matrix(5,5,0.);
-
- if(_materialSite && _applyMSCorrection) {
- //add here process noise due to dE/dx and scattering effects
-
- // Consider 3 unit vectors:
- // e1 = (px,py,pz)/p
- // e2 = (-py,px,0)/pt
- // e3 = (-px*pz,-py*pz,px*px+py*py)/(p*pt)
- // Original direction = e1
- // Scattered direction =
- // sqrt(1-sin^2(t2)*sin^2(t3))*e1 + sin(t2)*e2 + sin(t3)*e3
- //
- // 1. Determine d(d0,phi0,omega,z0,s)/d(theta1) and d(d0,phi0,omega,z0,s)/d(theta2)
- // Create tracks with momentum rotated by small amounts about e2 and e3
- // Measure difference in fit parameters due to these variations
- // --> get 5x2 matrix SCT = dp/dtheta
- // 2. Determine average scattering angle theta -> sigmaSCT
- // --> 2x2 covariance matrix sigSCT = ( sigmaSCT^2 0 )
- // ( 0 sigmaSCT^2 )
- // 3. Q_scattering = SCT x sigSCT x SCT^T
-
- // 1. Determine derivatives
- double sctVar = 0.000001; // variation of scattering angle
- double[] ptPar = trk.getPointParameters();
- double x = ptPar[0];
- double y = ptPar[1];
- double z = ptPar[2];
- double px = ptPar[3];
- double py = ptPar[4];
- double pz = ptPar[5];
- int chg = trk.getCharge();
-
- double pt = Math.sqrt(px*px+py*py);
- double p = Math.sqrt(px*px+py*py+pz*pz);
-
- double pxt1p = Math.cos(sctVar)*px - Math.sin(sctVar)*py*p/pt;
- double pyt1p = Math.cos(sctVar)*py + Math.sin(sctVar)*px*p/pt;
- double pzt1p = Math.cos(sctVar)*pz;
-
- double pxt1m = Math.cos(sctVar)*px + Math.sin(sctVar)*py*p/pt;
- double pyt1m = Math.cos(sctVar)*py - Math.sin(sctVar)*px*p/pt;
- double pzt1m = Math.cos(sctVar)*pz;
-
- double pxt2p = Math.cos(sctVar)*px - Math.sin(sctVar)*px*pz/pt;
- double pyt2p = Math.cos(sctVar)*py - Math.sin(sctVar)*py*pz/pt;
- double pzt2p = Math.cos(sctVar)*pz + Math.sin(sctVar)*pt;
-
- double pxt2m = Math.cos(sctVar)*px + Math.sin(sctVar)*px*pz/pt;
- double pyt2m = Math.cos(sctVar)*py + Math.sin(sctVar)*py*pz/pt;
- double pzt2m = Math.cos(sctVar)*pz - Math.sin(sctVar)*pt;
-
-
- KFTrackParameters trkVar1P = new KFTrackParameters();
- trkVar1P.setTrackPointParameters(chg,x,y,z,pxt1p,pyt1p,pzt1p);
- double[] fit1P = trkVar1P.getFitParameters();
- KFTrackParameters trkVar1M = new KFTrackParameters();
- trkVar1M.setTrackPointParameters(chg,x,y,z,pxt1m,pyt1m,pzt1m);
- double[] fit1M = trkVar1M.getFitParameters();
- double[] v1 = new double [5];
- for(int i1=0;i1<5;i1++) {v1[i1]=(fit1P[i1]-fit1M[i1])/(2.*sctVar);}
-
- KFTrackParameters trkVar2P = new KFTrackParameters();
- trkVar2P.setTrackPointParameters(chg,x,y,z,pxt2p,pyt2p,pzt2p);
- double[] fit2P = trkVar2P.getFitParameters();
- KFTrackParameters trkVar2M = new KFTrackParameters();
- trkVar2M.setTrackPointParameters(chg,x,y,z,pxt2m,pyt2m,pzt2m);
- double[] fit2M = trkVar2M.getFitParameters();
- double[] v2 = new double [5];
- for(int i2=0;i2<5;i2++) {v2[i2]=(fit2P[i2]-fit2M[i2])/(2.*sctVar);}
-
- double[][] vals = {{v1[0],v1[1],v1[2],v1[3],v1[4]},{v2[0],v2[1],v2[2],v2[3],v2[4]}};
- Matrix SCTderiv = new Matrix(vals);
- Matrix SCTderivT = SCTderiv.transpose();
- if(_cst.debug()) {
- System.out.println("KFSite: SCTderiv:");
- SCTderiv.print(15,12);
- }
-
- // 2. Average scattering angle
- double mass = _cst.muonMass(); // assume muon hypothesis
- double beta = p/Math.sqrt(p*p+mass*mass); // use particle momentum
- double cosEntranceAngle = Math.abs((x*px+y*py+z*pz)/(Math.sqrt(x*x+y*y+z*z)*p));
- if(_cst.debug()) System.out.println("KFSite: MSct cos(entrance angle) = "+cosEntranceAngle);
- double pathLength = _radThickness/cosEntranceAngle;
- double sctRMS = (0.0136/(beta*p))*Math.sqrt(pathLength)*(1+0.038*Math.log(pathLength));
-
- double[][] valsSCT = {{sctRMS*sctRMS,0.},{0.,sctRMS*sctRMS}};
- Matrix sigSCT = new Matrix(valsSCT);
- if(_cst.debug()) {
- System.out.println("KFSite: sigSCT:");
- sigSCT.print(15,12);
- }
-
- // 3. determine Q matrix for scattering
- Q = SCTderivT.times(sigSCT).times(SCTderiv);
- if(_cst.debug()) {
- System.out.println("KFSite: Q:");
- Q.print(15,12);
- }
- }
-
- Matrix trkCov = trkCovNoNoise.plus(Q);
- if(_cst.debug()) {
- System.out.println("KFSite: trkCov:");
- trkCov.print(15,12);
- }
-
- Matrix filteredCov = new Matrix(5,5);
-
- if(_measurementSite) {
- // define matrix relating the measured parameters (phi,z) at radius R and the
- // state vector (fit parameters). In the current setup, the state vector is
- // given by the fit parameters (d0,phi0,omega,z0,s)
-
- double[] vP = new double [5];
- double[] vZ = new double [5];
-
- double[] variations = new double[] {1e-05,1e-05,1e-10,1e-05,1e-06};
-
- double[] fitPars = trk.getFitParameters();
- double pthlength = 0.;
-
- for(int j=0; j<5; j++) {
- double p1p = fitPars[0]; double p1m = fitPars[0];
- if(j==0) { p1p += 0.5*variations[0]; p1m -= 0.5*variations[0]; }
- double p2p = fitPars[1]; double p2m = fitPars[1];
- if(j==1) { p2p += 0.5*variations[1]; p2m -= 0.5*variations[1]; }
- double p3p = fitPars[2]; double p3m = fitPars[2];
- if(j==2) { p3p += 0.5*variations[2]; p3m -= 0.5*variations[2]; }
- double p4p = fitPars[3]; double p4m = fitPars[3];
- if(j==3) { p4p += 0.5*variations[3]; p4m -= 0.5*variations[3]; }
- double p5p = fitPars[4]; double p5m = fitPars[4];
- if(j==4) { p5p += 0.5*variations[4]; p5m -= 0.5*variations[4]; }
-
- KFTrackParameters trkVarPlus = new KFTrackParameters();
- trkVarPlus.setTrackFitParameters(p1p, p2p, p3p, p4p, p5p,pthlength);
- KFTrackParameters trkVarPlusPropagated =
- trkVarPlus.makePredictionAtRadius(_siteLocation.getR());
- KFTrackParameters trkVarMinus = new KFTrackParameters();
- trkVarMinus.setTrackFitParameters(p1m, p2m, p3m, p4m, p5m,pthlength);
- KFTrackParameters trkVarMinusPropagated =
- trkVarMinus.makePredictionAtRadius(_siteLocation.getR());
-
- vP[j] = (trkVarPlusPropagated.getPhi()-trkVarMinusPropagated.getPhi())/variations[j];
- vZ[j] = (trkVarPlusPropagated.getZ()-trkVarMinusPropagated.getZ())/variations[j];
-
- if(_cst.debug()) System.out.println("KFSite: phi+ = "+trkVarPlusPropagated.getPhi()+" phi- = "+trkVarMinusPropagated.getPhi()+" z+ = "+trkVarPlusPropagated.getZ()+" z- = "+trkVarMinusPropagated.getZ());
- }
-
- double[][] vals = {{vP[0],vP[1],vP[2],vP[3],vP[4]},{vZ[0],vZ[1],vZ[2],vZ[3],vZ[4]}};
- Matrix H = new Matrix(vals);
-
- Matrix HT = H.transpose();
- if(_cst.debug()) {
- System.out.println("KFSite: H:");
- H.print(12,10);
- System.out.println("KFSite: HT:");
- HT.print(12,10);
- }
-
- Matrix trkCovInverse = trkCov.inverse();
-
- if(_cst.debug()) {
- System.out.println("KFSite: trkCovInverse:");
- trkCovInverse.print(15,12);
- }
-
- double[] ptParam = trk.getPointParameters();
- double refPhi = Math.atan2(ptParam[1],ptParam[0]);
- double refZ = ptParam[2];
- double[] mp = {_siteLocation.getPhi()-refPhi,_siteLocation.getZ()-refZ};
- //double[] mp = {_siteLocation.getPhi(),_siteLocation.getZ()};
- Matrix measVect = new Matrix(mp,2);
-
- if(_cst.debug()) {
- System.out.println("KFSite: measVect:");
- measVect.print(12,10);
- System.out.println("KFSite: PhiZ cov matrix inverted:");
- _PhiZErrorMatrixInverse.print(12,10);
- }
-
- Matrix cv1 = _PhiZErrorMatrixInverse.times(H);
- Matrix cv2 = HT.times(cv1);
- if(_cst.debug()) {
- System.out.println("KFSite: measCovInverse:");
- cv2.print(15,12);
- }
- Matrix filteredCovInv = trkCovInverse.plus(cv2);
- filteredCov = filteredCovInv.inverse();
-
- if(_cst.debug()) {
- System.out.println("KFSite: Filtered Cov:");
- filteredCov.print(15,12);
- }
-
- Matrix tmp1 = filteredCovInv.times(stateVect);
- Matrix tmp2 = HT.times(_PhiZErrorMatrixInverse).times(measVect);
- Matrix tmp3 = tmp1.plus(tmp2);
- Matrix filteredTrkPar = filteredCov.times(tmp3);
-
- if(_cst.debug()) {
- System.out.println("KFSite: Filtered State Vector:");
- filteredTrkPar.print(15,12);
- }
-
- double[][] par = filteredTrkPar.getArrayCopy();
- d0 = par[0][0];
- phi0 = par[1][0];
- omega = par[2][0];
- z0 = par[3][0];
- tanl = par[4][0];
-
- if(_cst.debug()) {
- System.out.println("KFSite: Fitted Pt = "+_cst.pointThreeTimesB()/Math.abs(omega));
- System.out.println("KFSite: pars = "+d0+" "+phi0+" "+omega+" "+z0+" "+tanl);
- }
-
- } //if(_measurementSite) {
-
- KFTrackParameters tmpTrk = new KFTrackParameters();
- tmpTrk.setTrackFitParameters(d0,phi0,omega,z0,tanl,l);
- tmpTrk.setPointToRadius(_siteLocation.getR());
-
- double[] tmpPar = tmpTrk.getPointParameters();
- int tmpCh = tmpTrk.getCharge();
-
- double tmpPtot = tmpTrk.getPtot();
- double ratioCorr = 1.;
- Matrix QE = new Matrix(5,5,0.);
- if(_materialSite && _applyDEDXCorrection) {
- double dEdx = dEdx(tmpPtot);
- if(_cst.debug()) {
- System.out.println("KFSite: Ptot = "+tmpPtot);
- System.out.println("KFSite: dEdx [GeV] = "+_thickness*dEdx);
- }
- double tmpEtot = Math.sqrt(_cst.muonMass()*_cst.muonMass()+tmpPtot*tmpPtot);
- double tmpEtotCorr = tmpEtot;
-
- double[] pP = tmpTrk.getPointParameters();
- double cosEntranceAngle = Math.abs((pP[0]*pP[3]+pP[1]*pP[4]+pP[2]*pP[5])/
- (Math.sqrt(pP[0]*pP[0]+pP[1]*pP[1]+pP[2]*pP[2])*tmpPtot));
- if(_cst.debug()) System.out.println("KFSite: dEdx cos(entrance angle) = "+cosEntranceAngle);
- if(_fitOutward) {
- tmpEtotCorr -= dEdx*_thickness/cosEntranceAngle; // apply energy _loss_ when fitting inside-out
- } else {
- tmpEtotCorr += dEdx*_thickness/cosEntranceAngle; // apply energy _gain_ when fitting outside-in
- }
- ratioCorr = tmpEtotCorr/tmpEtot;
-
- // Determine uncertainty on dEdx. Use value from
- // H. Bicshel, Rev.Mod.Phys. 60 (1988) 663, table VI, p.682:
- // w = 31.54keV for 320um Si => 0.000985625GeV/cm
- double dESig = 0.000985625*_thickness/cosEntranceAngle;
- if(_cst.debug()) System.out.println("KFSite: dEdx sigma = "+dESig);
-
- Matrix sigDEDX = new Matrix(1,1,dESig*dESig);
- double deltadEdx = dEdx/1000.;
- double ratioCorrPlus = (tmpEtot+(dEdx+deltadEdx)*_thickness/cosEntranceAngle)/tmpEtot;
- double ratioCorrMinus = (tmpEtot+(dEdx-deltadEdx)*_thickness/cosEntranceAngle)/tmpEtot;
-
- KFTrackParameters trkVardEdxPlus = new KFTrackParameters();
- trkVardEdxPlus.setTrackPointParameters(tmpCh,
- tmpPar[0],
- tmpPar[1],
- tmpPar[2],
- ratioCorrPlus*tmpPar[3],
- ratioCorrPlus*tmpPar[4],
- ratioCorrPlus*tmpPar[5]);
- double[] fitPardEdxPlus = trkVardEdxPlus.getFitParameters();
- KFTrackParameters trkVardEdxMinus = new KFTrackParameters();
- trkVardEdxMinus.setTrackPointParameters(tmpCh,
- tmpPar[0],
- tmpPar[1],
- tmpPar[2],
- ratioCorrMinus*tmpPar[3],
- ratioCorrMinus*tmpPar[4],
- ratioCorrMinus*tmpPar[5]);
- double[] fitPardEdxMinus = trkVardEdxMinus.getFitParameters();
- double[] jacdEdx = new double[5];
- for(int i=0;i<5;i++) {jacdEdx[i]=(fitPardEdxPlus[i]-fitPardEdxMinus[i])/(2.*deltadEdx);}
- Matrix DEDXderiv = new Matrix(jacdEdx,1);
- if(_cst.debug()) {
- System.out.println("KFSite: DEDXderiv:");
- DEDXderiv.print(15,12);
- }
- Matrix DEDXderivT = DEDXderiv.transpose();
- if(_cst.debug()) {
- System.out.println("KFSite: DEDXderivT:");
- DEDXderivT.print(15,12);
- }
- QE = DEDXderivT.times(sigDEDX).times(DEDXderiv);
- }
-
- KFTrackParameters filtTrk = new KFTrackParameters();
-
- filtTrk.setTrackPointParameters(tmpCh,
- tmpPar[0],
- tmpPar[1],
- tmpPar[2],
- ratioCorr*tmpPar[3],
- ratioCorr*tmpPar[4],
- ratioCorr*tmpPar[5]);
-
- if(_cst.debug()) {
- System.out.println("KFSite: QE:");
- QE.print(15,12);
- }
-
- if(_measurementSite) {
- filtTrk.setErrorMatrix(filteredCov.plus(QE));
- } else {
- filtTrk.setErrorMatrix(trkCov.plus(QE));
- }
- filtTrk.setPointToRadius(_siteLocation.getR());
-
- return filtTrk;
- }
-
- public void setTrackParOut(KFTrackParameters trk) {_trackParamOutward = trk;}
- public void setTrackParIn(KFTrackParameters trk) {_trackParamInward = trk;}
-
- public KFPoint getSiteLocation() {return _siteLocation;}
-
- public void setFitOutward() {_fitOutward = true;}
- public void setFitInward() {_fitOutward = false;}
-
- public void isMeasurementSite(boolean bool) {_measurementSite=bool;}
- public void isMaterialSite(boolean bool,double rT, double t) {
- _materialSite=bool;
- _radThickness = rT;
- _thickness = t;
- }
- public void isBFieldDiscontinuitySite(boolean bool) {_BFieldDiscontinuitySite=bool;}
-
- public void applyMSCorrection(boolean bool) {_applyMSCorrection=bool;}
- public void applyDEDXCorrection(boolean bool) {
- if(_cst.debug()) System.out.print("KFSite: Setting dedx bool from: "+_applyDEDXCorrection+" to ");
- _applyDEDXCorrection=bool;
- if(_cst.debug()) System.out.println(_applyDEDXCorrection);
- }
-
- private double dEdx(double p)
- {
- double dedx = 0.0;
- double rho = 2.3296; // g/cm^3 (Si)
- double Z = 14.; // (Si)
- double A = 28.0855; // g/mol (Si)
- double beta = p/Math.sqrt(p*p+_cst.muonMass()*_cst.muonMass());
- double gamma = 1./Math.sqrt(1-beta*beta);
- double eta = beta*gamma;
- double s = _cst.electronMass()/_cst.muonMass();
- double Wmax = 2.*_cst.electronMass()*eta*eta/(1.+2.*s*Math.sqrt(1+eta*eta)+s*s);
- double I = 0.000000172;
-
- double frac = 2.*_cst.electronMass()*eta*eta*Wmax/(I*I);
-
- double delta = 0.;
- double C = 0.;
-
- dedx = 0.0001535*rho*(Z/A)*(Math.log(frac)-2*beta*beta-delta-2*C/Z)/(beta*beta);
-
- return dedx;
- }
-
- private KFTrackParameters _trackParamOutward;
- private KFTrackParameters _trackParamOutwardCombined;
- private KFTrackParameters _trackParamInward;
- private KFTrackParameters _trackParamInwardCombined;
-
- private KFPoint _siteLocation = new KFPoint();
-
- private Matrix _PhiZErrorMatrix = new Matrix(2,2);
- private Matrix _PhiZErrorMatrixInverse = new Matrix(2,2);
-
- private boolean _measurementSite = false;
- private boolean _materialSite = false;
- private boolean _BFieldDiscontinuitySite = false;
-
- private boolean _applyMSCorrection = true;
- private boolean _applyDEDXCorrection = true;
-
- private boolean _fitOutward = true;
-
- private double _radThickness = 0.;
- private double _thickness = 0.;
-
- private KFConstants _cst = new KFConstants();
-}
lcsim/src/org/lcsim/contrib/KFFitter
diff -N KFTrack.java
--- KFTrack.java 13 Apr 2007 20:32:06 -0000 1.5
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,331 +0,0 @@
-package org.lcsim.contrib.KFFitter;
-
-import org.lcsim.event.Track;
-import org.lcsim.event.SimTrackerHit;
-import org.lcsim.event.TrackerHit;
-import org.lcsim.event.MCParticle;
-import org.lcsim.contrib.SODTracker.SODTrack;
-
-import java.util.ArrayList;
-import java.util.*;
-
-import Jama.Matrix;
-
-/**
- * Define a track used in the Kalman filter fit
- * A track is a collection of KFSites.
- * The KFTrack class handles the loop over sites to apply the Kalman fit
- * @author Fred Blanc and Steve Wagner
- * @version $Id: KFTrack.java,v 1.5 2007/04/13 20:32:06 fblanc Exp $
- */
-
-public class KFTrack
-{
- public KFTrack()
- {
- _siteList = new java.util.LinkedList();
- _orderedSiteList = new java.util.LinkedList();
- }
-
- public void addHit(SimTrackerHit hit, double rdm1, double rdm2)
- {
- // define site for this hit
- KFSite kfhit = new KFSite();
- kfhit.setLocation(hit,rdm1,rdm2);
-
- kfhit.isMeasurementSite(true);
- double radius = kfhit.getSiteLocation().getR();
- if(_cst.debug()) System.out.println("Track addHit at Radius "+radius);
- double distance = kfhit.getSiteLocation().getDistanceFromOrigin();
- if(_cst.debug()) System.out.println("Track addHit at Distance "+distance);
- double radThickness = 0.00121;
- if (radius<2.0) {radThickness = 0.00301;} // add beam pipe 0.18%
- if (radius>10.) {radThickness = 0.00532;} // tracker
- if (radius>10.&&radius<30.) {radThickness = 0.01064;} // + 2 layers of support tube 0.266%
- //double thickness = 0.0113;
- //if (radius>10.) {thickness = 0.031;}
- double thickness = 0.0170;
- if (radius>10.) {thickness = 0.047;}
- kfhit.isMaterialSite(true,radThickness,thickness);
-
- kfhit.applyMSCorrection(true);
- kfhit.applyDEDXCorrection(true);
-
- _siteList.add(kfhit);
-
- }
-
- private void sortSiteList()
- {
- if(_cst.debug()) System.out.println("KFTrack sortSiteList");
- // loop over _siteList
- ListIterator hits = _siteList.listIterator();
- int hitCounter = 0;
- double[] dist = new double[100];
- if(_cst.debug()) System.out.print(" init list radii = ");
- for (; hits.hasNext();) {
- KFSite hit = (KFSite)hits.next();
- // record ordering parameter (e.g. distance from origin)
- dist[hitCounter] = hit.getSiteLocation().getDistanceFromOrigin();
- if(_cst.debug()) System.out.print(dist[hitCounter]+"; ");
- hitCounter++;
- }
- if(_cst.debug()) System.out.println();
- // loop over _siteList and fill _orderedSiteList
- double currentDist = 0.;
- int index = 0;
- if(_cst.debug()) System.out.print(" ordered list radii = ");
- for (int j=0;j<hitCounter;j++) {
- double nextSmallestDist = 100000.;
- for (int i=0;i<hitCounter;i++) {
- if(dist[i]<nextSmallestDist&&dist[i]>currentDist) {
- index=i;
- nextSmallestDist=dist[i];
- }
- }
- KFSite nextSite = (KFSite)_siteList.get(index);
- _orderedSiteList.add(nextSite);
- if(_cst.debug()) System.out.print(nextSite.getSiteLocation().getDistanceFromOrigin()+"; ");
- currentDist = nextSmallestDist;
- }
- if(_cst.debug()) System.out.println();
-
- // check ordering of hits
- ListIterator hts = _orderedSiteList.listIterator();
- if(_cst.debug()) System.out.print(" check order : ");
- for (; hts.hasNext();){
- KFSite ht = (KFSite)hts.next();
- double R = ht.getSiteLocation().getR();
- if(_cst.debug()) System.out.print(R+"; ");
- }
- if(_cst.debug()) System.out.println();
- }
-
- public void defineSeedTrack(SODTrack seedTrack)
- {
- // load track parameters d0,omega,...
- double d = seedTrack.getTrackParameter(0)/10.; // divided by 10 to get cm
- double phi = seedTrack.getTrackParameter(1);
- double omega = seedTrack.getTrackParameter(2)*(-10.); // mult. by -10 to get cm-1
- double z = seedTrack.getTrackParameter(3)/10.; // divided by 10 to get cm
- double s = seedTrack.getTrackParameter(4);
-
- double l = 0.0; // path length set to zero (=POCA)
-
- if(_cst.debug()) System.out.println("KFTrack: define track par = "+d+" "+phi+" "+omega+" "+z+" "+s);
-
- _seedTrack.setTrackFitParameters(d,phi,omega,z,s,l);
-
- }
-
- public void fit()
- {
- // sort list of hits by distance from origin
- sortSiteList();
-
- // initialize track param covariance matrix
- double[][] diag = new double[5][5];
- diag[0][0] = 1.; // d0
- diag[1][1] = 0.25; // phi0
- diag[2][2] = 0.1; // omega
- diag[3][3] = 10.; // z0
- diag[4][4] = 1.; // tanl
- Matrix initCov = new Matrix(diag);
-
- if(_cst.debug()) {
- System.out.println("KFTrack: init cov matrix");
- initCov.print(6,2);
- System.out.print("KFTrack: ");
- _seedTrack.printFitParams();
- }
-
- ListIterator hits = _orderedSiteList.listIterator();
-
- //***************
- // Fit inside-out
- //***************
- if(_cst.debug()) System.out.println("KFTrack: Fitting Track outward");
-
- // initialize track parameters for outward fit
- KFTrackParameters currentTrackOut = new KFTrackParameters(_seedTrack);
-
- /*
- // Determine seed track parameters from 3 innermost hits
- KFTrackParameters inOutSeedTrk =
- seedTrkFromHits(_seedTrack.getCharge(),
- (KFSite)_orderedSiteList.get(0),
- (KFSite)_orderedSiteList.get(1),
- (KFSite)_orderedSiteList.get(2));
- currentTrackOut = inOutSeedTrk;
- */
-
- // determine track params at fixed small radius
- //_seedTrack.setPointToPOCA();
- double refRadius = 1.0;
- currentTrackOut.setPointToRadius(refRadius);
- if(_cst.debug()) {
- System.out.print("KFTrack: set radius to 1 : ");
- currentTrackOut.printFitParams();
- }
- // initialize error matrix
- currentTrackOut.setErrorMatrix(initCov);
-
- for (; hits.hasNext();){
-
- KFSite hit = (KFSite)hits.next();
- double R = hit.getSiteLocation().getR();
- if(_cst.debug()) System.out.println("KFTrack: hit radius = "+R);
-
- if(_cst.debug()) {
- System.out.print("KFTrack: current pars = ");
- currentTrackOut.printFitParams();
- }
-
- // propagate track to radius of hit
- KFTrackParameters predictedTrack = currentTrackOut.makePredictionAtRadius(R);
-
- if(_cst.debug()) {
- System.out.print("KFTrack: predicted pars = ");
- predictedTrack.printFitParams();
- predictedTrack.printPointParams();
- }
-
- // filter predicted track parameters with hit at radius R
- hit.setFitOutward();
- currentTrackOut = hit.filter(predictedTrack);
-
- // record track parameters in KFSite
- hit.setTrackParOut(currentTrackOut);
- }
-
- //***************
- // Fit outside-in
- //***************
- if(_cst.debug()) System.out.println("KFTrack: Fitting Track inward");
-
- // initialize track parameters for inward fit
- KFTrackParameters currentTrackIn = new KFTrackParameters(currentTrackOut);
-
- // determine track params at fixed small radius
- refRadius = 130.0;
- currentTrackIn.setPointToRadius(refRadius);
- if(_cst.debug()) {
- System.out.print("KFTrack: set radius to 130 : ");
- currentTrackIn.printFitParams();
- }
-
- // initialize error matrix
- currentTrackIn.setErrorMatrix(initCov);
-
- for (; hits.hasPrevious();){
-
- KFSite hit = (KFSite)hits.previous();
- double R = hit.getSiteLocation().getR();
- if(_cst.debug()) {
- System.out.println("KFTrack: hit radius = "+R);
- System.out.print("KFTrack: current pars = ");
- currentTrackIn.printFitParams();
- }
-
- // propagate track to radius of hit
- KFTrackParameters predictedTrack = currentTrackIn.makePredictionAtRadius(R);
-
- if(_cst.debug()) {
- System.out.print("KFTrack: predicted pars = ");
- predictedTrack.printFitParams();
- }
-
- // filter predicted track parameters with hit at radius R
- hit.setFitInward();
- currentTrackIn = hit.filter(predictedTrack);
-
- // record track parameters in KFSite
- hit.setTrackParIn(currentTrackIn);
-
- }
-
- if(_applyVertexConstraint) {
- KFVertex origin = new KFVertex();
- currentTrackIn = origin.filter(currentTrackIn);
- }
-
- _KFTrackAtPOCA = currentTrackIn;
-
- _fittedPtAtPOCA = currentTrackIn.getPt();
-
- if(_cst.debug()) {
- System.out.println("KFTrack: Final fitted pt = "+_fittedPtAtPOCA);
- System.out.print("KFTrack: Final fitted par = ");
- currentTrackIn.printFitParams();
- System.out.print("KFTrack: Final cov matrix = ");
- currentTrackIn.getErrorMatrix().print(15,12);
- Matrix cov = currentTrackIn.getErrorMatrix();
- double[][] covArray = cov.getArrayCopy();
- double[] error = new double[5];
- for(int c=0; c<5; c++) {error[c] = Math.sqrt(covArray[c][c]);}
- for(int a=0; a<5; a++) {
- for(int b=0; b<5; b++) {
- covArray[a][b] /= error[a]*error[b];
- }
- }
- Matrix corMat = new Matrix(covArray);
- System.out.print("KFTrack: Final correlation matrix = ");
- corMat.print(15,12);
- }
- }
-
- private KFTrackParameters seedTrkFromHits(int q, KFSite p1, KFSite p2, KFSite p3)
- {
- double x1 = p1.getSiteLocation().getX();
- double y1 = p1.getSiteLocation().getY();
- double z1 = p1.getSiteLocation().getZ();
- double x2 = p2.getSiteLocation().getX();
- double y2 = p2.getSiteLocation().getY();
- double z2 = p2.getSiteLocation().getZ();
- double x3 = p3.getSiteLocation().getX();
- double y3 = p3.getSiteLocation().getY();
-
- double r1 = Math.sqrt(x1*x1+y1*y1);
- double r2 = Math.sqrt(x2*x2+y2*y2);
-
- double denom = (x1-x2)*(y3-y2)-(x2-x3)*(y2-y1);
- double alpha = 0.5*((x1-x3)*(x2-x3)+(y1-y3)*(y2-y3))/denom;
- double beta = 0.5*((x1-x3)*(x1-x2)+(y1-y3)*(y1-y2))/denom;
-
- double xc = 0.5*(x1+x2)+alpha*(y2-y1);
- double yc = 0.5*(y1+y2)+beta*(x1-x2);
-
- double r0 = Math.sqrt(xc*xc+yc*yc);
- double rc = Math.sqrt((xc-x1)*(xc-x1)+(yc-y1)*(yc-y1));
-
- double d0 = r0-rc;
- double phi0 = Math.atan2(q*xc,-q*yc);
- double omega = -q/rc;
-
- double s = (z2-z1)/(r2-r1);
- double z0 = z1+(Math.abs(d0)-r1)*s;
-
- double l = 0.;
-
- KFTrackParameters seedTrk = new KFTrackParameters();
- seedTrk.setTrackFitParameters(d0,phi0,omega,z0,s,l);
-
- return seedTrk;
- }
-
- public KFTrackParameters getSeedTrackParameters() { return _seedTrack; }
- public KFTrackParameters getKFTrackParameters() { return _KFTrackAtPOCA; }
- public double getFittedPtAtPOCA() { return _fittedPtAtPOCA; }
-
- private java.util.LinkedList _siteList;
- private java.util.LinkedList _orderedSiteList;
- private KFTrackParameters _seedTrack = new KFTrackParameters();
- private KFTrackParameters _KFTrackAtPOCA = new KFTrackParameters();
-
- private double _fittedPtAtPOCA;
-
- private boolean _applyVertexConstraint = true;
-
- private KFConstants _cst = new KFConstants();
-
-}
lcsim/src/org/lcsim/contrib/KFFitter
diff -N KFTrackParameters.java
--- KFTrackParameters.java 13 Apr 2007 20:32:06 -0000 1.5
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,276 +0,0 @@
-package org.lcsim.contrib.KFFitter;
-
-import Jama.Matrix;
-
-/**
- * Define track parameters
- * Allows multiple representations, and the transformation from one to the other.
- * Holds the covariance matrix in the standrd (d0,phi0,omega,z0,tanl) representation
- * @author Fred Blanc and Steve Wagner
- * @version $Id: KFTrackParameters.java,v 1.5 2007/04/13 20:32:06 fblanc Exp $
- */
-public class KFTrackParameters
-{
-
- KFTrackParameters()
- {
-
- }
-
- KFTrackParameters(KFTrackParameters trk)
- {
- double[] par = trk.getFitParameters();
- double l = trk.getPathLength();
- this.setTrackFitParameters(par[0], par[1], par[2], par[3], par[4], l);
-
- this._FitErrorMatrix = trk.getErrorMatrix();
- }
-
- private int _charge = 0;
-
- private double _d0 = 0., _phi0 = 0., _omega = 0., _z0 = 0., _s = 0.;
- private double _l = 0., _alpha = 0.;
- private double[] _fitParams = new double[] {0.,0.,0.,0.,0.};
-
- private double _xc = 0., _yc = 0., _r0 = 0., _rc = 0, _pt = 0.;
- private double[] _helixParams = new double[] {0.,0.,0.,0.,0.};
-
- private double _px = 0., _py = 0., _pz = 0., _p = 0.;
- private double _x = 0., _y = 0., _z = 0.;
- private double[] _pointParams = new double[] {0.,0.,0.,0.,0.,0.};
-
- public void setTrackFitParameters(double dzero, double phizero, double omeg,
- double zzero, double spar, double pathLength)
- {
- _charge = (int) (-1*Math.signum(omeg));
- //System.out.println("charge = "+_charge);
- _d0 = dzero;
- _phi0 = phizero;
- _omega = omeg; if(_omega==0.) {System.out.println("Null curvature!");}
- _z0 = zzero;
- _s = spar;
-
- _fitParams[0] = _d0;
- _fitParams[1] = _phi0;
- _fitParams[2] = _omega;
- _fitParams[3] = _z0;
- _fitParams[4] = _s;
-
- _l = pathLength;
- _alpha = _l*_omega;
-
- makePointFromFitRepresentation();
- //makeHelixFromFitRepresentation();
-
- }
-
- public void setTrackPointParameters(int ch, double x, double y, double z, double px, double py, double pz)
- {
- _charge = ch;
- _x = x;
- _y = y;
- _z = z;
- _px = px;
- _py = py;
- _pz = pz;
-
- _p = Math.sqrt(_px*_px + _py*_py + _pz*_pz);
-
- _pointParams[0] = _x;
- _pointParams[1] = _y;
- _pointParams[2] = _z;
- _pointParams[3] = _px;
- _pointParams[4] = _py;
- _pointParams[5] = _pz;
-
- makeFitFromPointRepresentation();
- //makeHelixFromPointRepresentation();
-
- }
-
- private void makePointFromFitRepresentation()
- {
- _pt = _cst.pointThreeTimesB()/Math.abs(_omega);
- _px = _pt*(Math.cos(_alpha)*Math.cos(_phi0)-Math.sin(_alpha)*Math.sin(_phi0));
- _py = _pt*(Math.sin(_alpha)*Math.cos(_phi0)+Math.cos(_alpha)*Math.sin(_phi0));
- _pz = _pt*_s;
- _p = Math.sqrt(_px*_px + _py*_py + _pz*_pz);
-
- _x = _d0*_charge*Math.sin(_phi0)+(_charge/Math.abs(_omega))*
- (Math.sin(_phi0)-Math.cos(_alpha)*Math.sin(_phi0)-Math.sin(_alpha)*Math.cos(_phi0));
- _y = -_d0*_charge*Math.cos(_phi0)-(_charge/Math.abs(_omega))*
- (Math.cos(_phi0)+Math.sin(_alpha)*Math.sin(_phi0)-Math.cos(_alpha)*Math.cos(_phi0));
- //_z = _z0+_s*Math.abs(_alpha/_omega);
- _z = _z0+_s*_l;
-
- _pointParams[0] = _x;
- _pointParams[1] = _y;
- _pointParams[2] = _z;
- _pointParams[3] = _px;
- _pointParams[4] = _py;
- _pointParams[5] = _pz;
- }
-
- private void makeHelixFromFitRepresentation()
- {
- _xc = _charge*Math.sin(_phi0)*(_d0+1./Math.abs(_omega));
- _yc = -_charge*Math.cos(_phi0)*(_d0+1./Math.abs(_omega));
- _r0 = Math.sqrt(_xc*_xc+_yc*_yc);
- _rc = 1./Math.abs(_omega);
- _pt = _cst.pointThreeTimesB()/Math.abs(_omega);
- //_s = _s;
- //_z0 = _z0;
-
- _helixParams[0] = _xc;
- _helixParams[1] = _yc;
- _helixParams[2] = _rc;
- _helixParams[3] = _z0;
- _helixParams[4] = _s;
-
- }
-
- private void makeHelixFromPointRepresentation()
- {
- _xc = _x+_charge*_py/_cst.pointThreeTimesB();
- _yc = _y-_charge*_px/_cst.pointThreeTimesB();
- _r0 = Math.sqrt(_xc*_xc+_yc*_yc);
- _rc = Math.sqrt((_x-_xc)*(_x-_xc)+(_y-_yc)*(_y-_yc));
- _pt = Math.sqrt(_px*_px+_py*_py); // = _rc*_cst.pointThreeTimesB();
- _s = _pz/_pt;
- _z0 = _z+Math.atan2(_y*_xc-_x*_yc,_r0*_r0-_x*_xc-_y*_yc)*_pz/_cst.pointThreeTimesB();
-
- _helixParams[0] = _xc;
- _helixParams[1] = _yc;
- _helixParams[2] = _rc;
- _helixParams[3] = _z0;
- _helixParams[4] = _s;
-
- }
-
- private void makeFitFromPointRepresentation()
- {
- double tmpXc = _x+_charge*_py/_cst.pointThreeTimesB();
- double tmpYc = _y-_charge*_px/_cst.pointThreeTimesB();
- double tmpPt = Math.sqrt(_px*_px+_py*_py);
- _d0 = Math.sqrt(tmpXc*tmpXc+tmpYc*tmpYc)-tmpPt/_cst.pointThreeTimesB();
- _phi0 = Math.atan2(_charge*tmpXc,-_charge*tmpYc);
- _omega = -_charge*_cst.pointThreeTimesB()/tmpPt;
- _alpha = Math.atan2(-_charge*(_x*_px+_y*_py),
- tmpPt*tmpPt/_cst.pointThreeTimesB()+_charge*(_x*_py-_y*_px));
- _z0 = _z+_alpha*_pz/(_charge*_cst.pointThreeTimesB());
- _s = _pz/tmpPt;
- _l = _alpha/_omega;
-
- _fitParams[0] = _d0;
- _fitParams[1] = _phi0;
- _fitParams[2] = _omega;
- _fitParams[3] = _z0;
- _fitParams[4] = _s;
- }
-
- private void makeFitFromHelixRepresentation()
- {
- _d0 = Math.sqrt(_xc*_xc+_yc*_yc)-_rc;
- _phi0 = Math.atan2(_charge*_xc,-_charge*_yc);
- _omega = -_charge/_rc;
- //_z0 = _z0;
- //_s = _s;
-
- _fitParams[0] = _d0;
- _fitParams[1] = _phi0;
- _fitParams[2] = _omega;
- _fitParams[3] = _z0;
- _fitParams[4] = _s;
- }
-
- public KFTrackParameters makePredictionAtRadius(double newRadius) {
-
- // in fit representation, only the path length changes when propagating
- // to a specific radius
-
- double newPathLength = propagate(newRadius);
-
- KFTrackParameters newTrack = new KFTrackParameters();
- newTrack.setTrackFitParameters(_d0,_phi0,_omega,_z0,_s,newPathLength);
- newTrack.setErrorMatrix(_FitErrorMatrix);
-
- return newTrack;
- }
-
- public double propagate(double r)
- {
- double par = 1.-((r*r-_d0*_d0)*_omega*_omega)/(2.*(1.+_d0*Math.abs(_omega)));
- double newl = Math.abs(Math.acos(par)/_omega);
- return newl;
- }
-
- public void setPointToRadius(double r) {
- double newPathLength = propagate(r);
- _l = newPathLength;
- _alpha = _l*_omega;
- makePointFromFitRepresentation();
- //makeHelixFromFitRepresentation();
- }
-
- public void setPointToPOCA() {
- _l = 0.;
- makePointFromFitRepresentation();
- //makeHelixFromFitRepresentation();
- }
-
- public double[] getFitParameters() {return _fitParams;}
- public double[] getPointParameters() {return _pointParams;}
- public double[] getHelixParameters() {return _helixParams;}
-
- public int getCharge() {return _charge;}
- public double getPathLength() {return _l;}
- public double getPt() {return _pt;}
- public double getPtot() {return _p;}
- public double getPhi() {return Math.atan2(_y,_x);}
- public double getZ() {return _z;}
-
- public void setErrorMatrix(Matrix cov) { _FitErrorMatrix = cov; }
-
- public Matrix getErrorMatrix() {return _FitErrorMatrix;}
-
- private Matrix _FitErrorMatrix = new Matrix(5,5);
-
- private KFConstants _cst = new KFConstants();
-
- public void resetTrackParameters() {
- int i;
- _charge = 0;
- _d0 = 0.; _phi0 = 0.; _omega = 0.; _z0 = 0.; _s = 0.;
- _l = 0; _alpha = 0.;
- for(i=0;i<5;i++){_fitParams[i]=0.;};
-
- _p = 0.; _pt = 0.; _px = 0.; _py = 0.; _pz = 0.;
- _x = 0.; _y = 0.; _z = 0.;
- for(i=0;i<6;i++){_pointParams[i]=0.;};
-
- _xc = 0.; _yc = 0.; _rc = 0.; _r0 = 0.;
- for(i=0;i<5;i++){_helixParams[i]=0.;};
-
- }
-
- public void printFitParams() {
- System.out.print("Fit params = ");
- System.out.print(_fitParams[0]+" ");
- System.out.print(_fitParams[1]+" ");
- System.out.print(_fitParams[2]+" ");
- System.out.print(_fitParams[3]+" ");
- System.out.print(_fitParams[4]+" ");
- System.out.println(_l);
- }
-
- public void printPointParams() {
- System.out.print("Point params = ");
- System.out.print(_pointParams[0]+" ");
- System.out.print(_pointParams[1]+" ");
- System.out.print(_pointParams[2]+" ");
- System.out.print(_pointParams[3]+" ");
- System.out.print(_pointParams[4]+" ");
- System.out.print(_pointParams[5]+" ");
- System.out.println(getPhi());
- }
-}
lcsim/src/org/lcsim/contrib/KFFitter
diff -N KFVertex.java
--- KFVertex.java 13 Apr 2007 20:32:06 -0000 1.1
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,139 +0,0 @@
-package org.lcsim.contrib.KFFitter;
-
-import java.util.ArrayList;
-import java.util.List;
-
-import Jama.Matrix;
-
-/**
- * Define a vertex used in the Kalman filter fit
- * A vertex is a point, with a given resolution
- * The class handles the filtering step of the Kalman filter.
- * @author Fred Blanc and Steve Wagner
- * @version $Id: KFVertex.java,v 1.1 2007/04/13 20:32:06 fblanc Exp $
- */
-
-public class KFVertex
-{
-
- public KFVertex()
- {
- }
-
- public KFTrackParameters filter(KFTrackParameters trk)
- {
-
- double[][] vals = {{4e-08,0.},{0.,4e-08}};
- Matrix errMat = new Matrix(vals);
- _DZErrorMatrix = errMat;
- _DZErrorMatrixInverse = _DZErrorMatrix.inverse();
-
-
- double[] tpp = trk.getFitParameters();
- Matrix stateVect = new Matrix(tpp,5);
- if(_cst.debug()) {
- System.out.println("KFVertex: initial stateVect:");
- stateVect.print(15,12);
- System.out.println("KFVertex: initial stateVect in point parameter space:");
- trk.printPointParams();
- }
-
- double d0 = tpp[0];
- double phi0 = tpp[1];
- double omega = tpp[2];
- double z0 = tpp[3];
- double tanl = tpp[4];
- double l = 0.;
-
- Matrix trkCov = trk.getErrorMatrix();
- if(_cst.debug()) {
- System.out.println("KFVertex: trkCov:");
- trkCov.print(15,12);
- }
-
- // define matrix relating the measured parameters (d,z) and the
- // state vector (fit parameters). In the current setup, the state vector is
- // given by the fit parameters (d0,phi0,omega,z0,s)
-
- double[][] valsErr = {{1.,0.,0.,0.,0.},{0.,0.,0.,1.,0.}};
- Matrix H = new Matrix(valsErr);
-
- Matrix HT = H.transpose();
- if(_cst.debug()) {
- System.out.println("KFVertex: H:");
- H.print(12,10);
- System.out.println("KFVertex: HT:");
- HT.print(12,10);
- }
-
- Matrix trkCovInverse = trkCov.inverse();
-
- if(_cst.debug()) {
- System.out.println("KFVertex: trkCovInverse:");
- trkCovInverse.print(15,12);
- }
-
- double refD = trk.getFitParameters()[0];
- double refZ = trk.getFitParameters()[3];
- double[] mp = {-refD,-refZ};
- Matrix measVect = new Matrix(mp,2);
-
- if(_cst.debug()) {
- System.out.println("KFVertex: measVect:");
- measVect.print(12,10);
- System.out.println("KFVertex: PhiZ cov matrix inverted:");
- _DZErrorMatrixInverse.print(12,10);
- }
-
- Matrix cv1 = _DZErrorMatrixInverse.times(H);
- Matrix cv2 = HT.times(cv1);
- if(_cst.debug()) {
- System.out.println("KFVertex: measCovInverse:");
- cv2.print(15,12);
- }
- Matrix filteredCovInv = trkCovInverse.plus(cv2);
- Matrix filteredCov = filteredCovInv.inverse();
-
- if(_cst.debug()) {
- System.out.println("KFVertex: Filtered Cov:");
- filteredCov.print(15,12);
- }
-
- Matrix tmp1 = filteredCovInv.times(stateVect);
- Matrix tmp2 = HT.times(_DZErrorMatrixInverse).times(measVect);
- Matrix tmp3 = tmp1.plus(tmp2);
- Matrix filteredTrkPar = filteredCov.times(tmp3);
-
- if(_cst.debug()) {
- System.out.println("KFVertex: Filtered State Vector:");
- filteredTrkPar.print(15,12);
- }
-
- double[][] par = filteredTrkPar.getArrayCopy();
- d0 = par[0][0];
- phi0 = par[1][0];
- omega = par[2][0];
- z0 = par[3][0];
- tanl = par[4][0];
-
- if(_cst.debug()) {
- System.out.println("KFVertex: Fitted Pt = "+_cst.pointThreeTimesB()/Math.abs(omega));
- System.out.println("KFVertex: pars = "+d0+" "+phi0+" "+omega+" "+z0+" "+tanl);
- }
-
- KFTrackParameters filtTrk = new KFTrackParameters();
- filtTrk.setTrackFitParameters(d0,phi0,omega,z0,tanl,l);
- filtTrk.setErrorMatrix(filteredCov);
-
- return filtTrk;
- }
-
- public KFPoint getVertexLocation() {return _vtxLocation;}
-
- private KFPoint _vtxLocation = new KFPoint();
-
- private Matrix _DZErrorMatrix = new Matrix(2,2);
- private Matrix _DZErrorMatrixInverse = new Matrix(2,2);
-
- private KFConstants _cst = new KFConstants();
-}
lcsim/src/org/lcsim/contrib/KFFitter/test
diff -N TestKFF.java
--- TestKFF.java 24 Jul 2007 22:23:25 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,64 +0,0 @@
-package org.lcsim.contrib.KFFitter.test;
-import java.util.List;
-
-import java.util.*;
-import java.io.*;
-import hep.physics.matrix.SymmetricMatrix;
-import org.lcsim.event.EventHeader;
-import org.lcsim.event.MCParticle;
-import org.lcsim.event.TrackerHit;
-import org.lcsim.event.base.BaseTrack;
-import org.lcsim.util.Driver;
-import org.lcsim.util.aida.AIDA;
-
-import org.lcsim.contrib.KFFitter.KFFitterDriver;
-
-
-public class TestKFF extends Driver {
- private AIDA aida = AIDA.defaultInstance();
- public TestKFF(){
- add(new KFFitterDriver());
- // try{
- // FileOutputStream fos = new FileOutputStream("output.dat");
- // PrintStream ps = new PrintStream(fos);
- // System.setOut(ps);
- // }
- // catch(FileNotFoundException x){
- // System.out.println("blow me!! output.dat doesn't exist");
- // x.printStackTrace();
- // }
- };
- public void process(EventHeader event) {
- super.process(event); // this takes care that the child Drivers are loaded and processed.
-
- List<BaseTrack> KFTrackList = (List<BaseTrack>) event.get("KFTracks");
- if(KFTrackList.size()<100) {
- aida.cloud1D("nKFTracks").fill(KFTrackList.size());
- for(BaseTrack trk : KFTrackList) {
- double[] p = trk.getMomentum();
- double pt = Math.sqrt(p[0]*p[0]+p[1]*p[1]);
- double ptot = Math.sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
- double[] trkPar = trk.getTrackParameters();
- SymmetricMatrix cov = trk.getErrorMatrix();
-
- aida.cloud1D("KFptot").fill(ptot);
- aida.cloud1D("KFpt").fill(pt);
- aida.cloud1D("KFd0").fill(trkPar[0]);
- aida.cloud1D("KFphi0").fill(trkPar[1]);
- aida.cloud1D("KFomega").fill(Math.abs(trkPar[2]));
- aida.cloud1D("KFz0").fill(trkPar[3]);
- aida.cloud1D("KFtanl").fill(trkPar[4]);
- //aida.cloud1D("KFtanl").fill(Math.abs(trkPar[4]));
-
- aida.cloud1D("KFd0Err").fill(Math.sqrt(cov.e(0,0)));
- aida.cloud1D("KFphiErr").fill(Math.sqrt(cov.e(1,1)));
- aida.cloud1D("KFomegaErr").fill(Math.sqrt(cov.e(2,2)));
- aida.cloud1D("KFz0Err").fill(Math.sqrt(cov.e(3,3)));
- aida.cloud1D("KFtanlErr").fill(Math.sqrt(cov.e(4,4)));
-
- // aida.cloud1D("KFd0").fill();
- // aida.cloud1D("KFd0").fill();
- }
- }
- }
-}
CVSspam 0.2.8