10 modified files
lcsim/src/org/lcsim/contrib/KFFitter
diff -u -r1.2 -r1.3
--- History 1 Feb 2007 22:29:34 -0000 1.2
+++ History 10 Mar 2007 00:27:25 -0000 1.3
@@ -3,6 +3,34 @@
## Please summarize changes to the KFFitter code here.
## Most recent first please.
+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.
lcsim/src/org/lcsim/contrib/KFFitter
diff -u -r1.2 -r1.3
--- KFConstants.java 30 Jan 2007 17:01:58 -0000 1.2
+++ KFConstants.java 10 Mar 2007 00:27:25 -0000 1.3
@@ -1,5 +1,11 @@
package org.lcsim.contrib.KFFitter;
+/**
+ * Definition of constants in KFFitter package
+ * @author Fred Blanc and Steve Wagner
+ * @version $Id: KFConstants.java,v 1.3 2007/03/10 00:27:25 fblanc Exp $
+ */
+
public class KFConstants
{
@@ -12,13 +18,21 @@
public double vtxResolution() {return vtxres;}
public double sodResolution() {return sodres;}
public double pointThreeTimesB() {return ptconst;}
- public double pionMass() {return pimass;}
+ public double pionMass() {return piMass;}
+ public double muonMass() {return muMass;}
+
+ 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.015; // = 0.3*B
+ private static double ptconst = 0.0149896229; // = 0.299792458*B
+
+ private static double piMass = 0.13957018;
+ private static double muMass = 0.105658369;
+
+ private static boolean DeBug = false;
- private static double pimass = 0.13957018;
}
+
lcsim/src/org/lcsim/contrib/KFFitter
diff -u -r1.3 -r1.4
--- KFFitterDriver.java 1 Feb 2007 22:29:34 -0000 1.3
+++ KFFitterDriver.java 10 Mar 2007 00:27:25 -0000 1.4
@@ -2,32 +2,55 @@
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.geometry.TrackerIDDecoder;
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.4 2007/03/10 00:27:25 fblanc 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);
+
+ //IAnalysisFactory af = IAnalysisFactory.create();
+ //ITree tree = af.createTreeFactory().create();
+ //IHistogramFactory hf = af.createHistogramFactory(tree);
+ //IHistogram1D histo1 = hf.createHistogram1D("histo1", 25, -5.0, 5.0);
+
}
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");
@@ -35,58 +58,159 @@
for ( SODTrack SODtrk : SODTrackList ) {
double trkPT = SODtrk.getPT();
- System.out.print("Driver: trk pT = "); System.out.println(trkPT);
+ aida.cloud1D("SODPt").fill(trkPT);
+ /*
+ System.out.println("Driver: trk pT = "+trkPT);
System.out.print("Driver: SOD track = ");
- System.out.print(SODtrk.getTrackParameter(0)); System.out.print(" ");
- System.out.print(SODtrk.getTrackParameter(1)); System.out.print(" ");
- System.out.print(SODtrk.getTrackParameter(2)); System.out.print(" ");
- System.out.print(SODtrk.getTrackParameter(3)); System.out.print(" ");
- System.out.print(SODtrk.getTrackParameter(4)); System.out.print(" ");
+ System.out.print(SODtrk.getTrackParameter(0)+" ");
+ System.out.print(SODtrk.getTrackParameter(1)+" ");
+ System.out.print(SODtrk.getTrackParameter(2)+" ");
+ System.out.print(SODtrk.getTrackParameter(3)+" ");
+ System.out.print(SODtrk.getTrackParameter(4)+" ");
System.out.println(SODtrk.getCharge());
+ */
// create new KFTrack and initialize it to the SODTrack fit parameters
KFTrack track = new KFTrack();
- track.defineSeedTrack(SODtrk);
- //track.defineSeedTrack(); // For testing purpose
+ track.defineSeedTrack(SODtrk);
+ //track.defineSeedTrackDebug(); // For testing purpose
List<SimTrackerHit> hits = (List<SimTrackerHit>) SODtrk.getHits();
int nbHits = hits.size();
- System.out.print("Driver: SOD hit list size = "); System.out.println(nbHits);
-
- double doca = 0.1*SODtrk.getTrackParameter(0);
- System.out.print("Driver: SOD doca = "); System.out.print(doca); System.out.println(" cm");
+ //System.out.println("Driver: SOD hit list size = "+nbHits);
double chi2 = SODtrk.getChi2();
- System.out.print("Driver: SOD chi2 = "); System.out.println(chi2);
- aida.cloud1D("chi2").fill(chi2);
+ //System.out.println("Driver: SOD chi2 = "+chi2);
+
double ndf = nbHits;
- if(ndf>0.) { double r = chi2/ndf; aida.cloud1D("chi2 per ndf").fill(r); }
+ //if(ndf>0.) { double r = chi2/ndf; aida.cloud1D("chi2 per ndf").fill(r); }
+
+ 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 ) {
- track.addHit(hit); // assign hits to KFTrack
- //track.addHitDebug(hit); // assign hits to KFTrack
+ // assign hits to KFTrack
+ //track.addHit(hit,0.,0.);
+ track.addHit(hit,gtor.nextGaussian(),gtor.nextGaussian());
+ //track.addHitDebug(hit);
+
+ 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;}
+ //if(genPt>0) {genTanl = Math.abs(genPz/genPt);}
+ double genX0 = hit.getMCParticle().getOriginX();
+ double genY0 = hit.getMCParticle().getOriginY();
+ genD0 = Math.sqrt(genX0*genX0+genY0+genY0);
+ genZ0 = hit.getMCParticle().getOriginZ();
}
- System.out.print("Driver: track params before fit = ");
- track.getTrackParameters().printFitParams();
- track.fit(); // run Kalman Filter Fit
- System.out.println();
-
- double KFpt = track.getFittedPtAtPOCA();
-
- // fill histograms
- aida.cloud1D("KFpt all").fill(KFpt);
- aida.cloud1D("SODpt all").fill(trkPT);
- aida.cloud2D("pt (helix vs KF) all").fill(KFpt,trkPT);
-
- if(Math.abs(KFpt-20.0)<0.1) {
- aida.cloud1D("KFpt").fill(track.getFittedPtAtPOCA());
- }
- if(Math.abs(trkPT-20.0)<0.1) {
- aida.cloud1D("SODpt").fill(trkPT);
- }
- if(Math.abs(KFpt-20.0)<0.1&&Math.abs(trkPT-20.0)<0.1) {
- aida.cloud2D("pt (helix vs KF)").fill(KFpt,trkPT);
+ //System.out.print("Driver: track params before fit = ");
+ //track.getSeedTrackParameters().printFitParams();
+ // *************************
+ // run Kalman Filter Fit
+ // *************************
+ track.fit();
+ //System.out.println();
+
+ //System.out.print("Driver: track params after fit = ");
+ //track.getKFTrackParameters().printFitParams();
+
+ // *****************************
+ // 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();
+ 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 {
@@ -94,5 +218,9 @@
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 -u -r1.1 -r1.2
--- KFMain.java 13 Dec 2006 22:45:12 -0000 1.1
+++ KFMain.java 10 Mar 2007 00:27:25 -0000 1.2
@@ -7,6 +7,12 @@
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.2 2007/03/10 00:27:25 fblanc Exp $
+ */
+
public class KFMain extends Driver
{
private int nEvents = 0;
lcsim/src/org/lcsim/contrib/KFFitter
diff -u -r1.1 -r1.2
--- KFMatrix.java 13 Dec 2006 22:45:12 -0000 1.1
+++ KFMatrix.java 10 Mar 2007 00:27:25 -0000 1.2
@@ -3,6 +3,13 @@
//import org.lcsim.spacegeom.Matrix;
import Jama.Matrix;
+/**
+ * Matrix class in KFFitter package
+ * (unused)
+ * @author Fred Blanc and Steve Wagner
+ * @version $Id: KFMatrix.java,v 1.2 2007/03/10 00:27:25 fblanc Exp $
+ */
+
public class KFMatrix
{
lcsim/src/org/lcsim/contrib/KFFitter
diff -u -r1.3 -r1.4
--- KFPoint.java 1 Feb 2007 22:29:34 -0000 1.3
+++ KFPoint.java 10 Mar 2007 00:27:25 -0000 1.4
@@ -2,6 +2,13 @@
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.4 2007/03/10 00:27:25 fblanc Exp $
+ */
+
public class KFPoint
{
@@ -38,6 +45,18 @@
}
+ public void smearRPhiPoint(double rdm1, double rdm2)
+ {
+ //System.out.println("KFPoint: Smearing "+rdm1+" "+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; }
lcsim/src/org/lcsim/contrib/KFFitter
diff -u -r1.3 -r1.4
--- KFSite.java 1 Feb 2007 22:29:34 -0000 1.3
+++ KFSite.java 10 Mar 2007 00:27:25 -0000 1.4
@@ -8,6 +8,15 @@
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.4 2007/03/10 00:27:25 fblanc Exp $
+ */
+
public class KFSite
{
@@ -15,34 +24,24 @@
{
}
- public void setLocation(SimTrackerHit hit)
+ 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]);
- System.out.print("KFSite: Add Pt at coord = ");
- System.out.print(pos[0]); System.out.print(" ");
- System.out.print(pos[1]); System.out.print(" ");
- System.out.println(pos[2]);
-
- //System.out.print("KFSite: read back radius = ");
- //System.out.println(_siteLocation.getR());
-
- //double rphierr = _siteLocation.getRPhiErr();
- //double zerr = _siteLocation.getZErr();
- //double[][] vals = {{rphierr,0.},{0.,zerr}};
_PhiZErrorMatrix = _siteLocation.getPhiZErrorMatrix();
_PhiZErrorMatrixInverse = _PhiZErrorMatrix.inverse();
- System.out.println("KFSite: PhiZ cov matrix:");
- _PhiZErrorMatrix.print(12,10);
-
- //System.out.println("KFSite: PhiZ cov matrix inverted:");
- //_PhiZErrorMatrixInverse.print(12,10);
-
+ if(_cst.debug()) {
+ System.out.println("KFSite: PhiZ cov matrix:");
+ _PhiZErrorMatrix.print(12,10);
+ }
}
public void setLocation(double x, double y, double z)
@@ -57,9 +56,7 @@
_siteLocation.makePoint(x,y,z);
System.out.print("KFSite: Add Pt at coord = ");
- System.out.print(x); System.out.print(" ");
- System.out.print(y); System.out.print(" ");
- System.out.println(z);
+ System.out.println(x+" "+y+" "+z);
_PhiZErrorMatrix = _siteLocation.getPhiZErrorMatrix();
@@ -77,8 +74,12 @@
double[] tpp = trk.getFitParameters();
Matrix stateVect = new Matrix(tpp,5);
- System.out.println("KFSite: initial stateVect:");
- stateVect.print(15,12);
+ 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];
@@ -88,8 +89,10 @@
double l = 0.;
Matrix trkCovNoNoise = trk.getErrorMatrix();
- System.out.println("KFSite: trkCovNoNoise:");
- trkCovNoNoise.print(15,12);
+ if(_cst.debug()) {
+ System.out.println("KFSite: trkCovNoNoise:");
+ trkCovNoNoise.print(15,12);
+ }
Matrix Q = new Matrix(5,5,0.);
@@ -165,28 +168,39 @@
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();
- System.out.println("KFSite: SCTderiv:");
- SCTderiv.print(15,12);
+ if(_cst.debug()) {
+ System.out.println("KFSite: SCTderiv:");
+ SCTderiv.print(15,12);
+ }
// 2. Average scattering angle
- double mass = _cst.pionMass(); // assume pion hypothesis
+ double mass = _cst.muonMass(); // assume muon hypothesis
double beta = p/Math.sqrt(p*p+mass*mass); // use particle momentum
- double sctRMS = (0.0136/(beta*p))*Math.sqrt(_radThickness)*(1+0.038*Math.log(_radThickness));
+ 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);
- System.out.println("KFSite: sigSCT:");
- sigSCT.print(15,12);
+ if(_cst.debug()) {
+ System.out.println("KFSite: sigSCT:");
+ sigSCT.print(15,12);
+ }
// 3. determine Q matrix for scattering
Q = SCTderivT.times(sigSCT).times(SCTderiv);
- System.out.println("KFSite: Q:");
- Q.print(15,12);
+ if(_cst.debug()) {
+ System.out.println("KFSite: Q:");
+ Q.print(15,12);
+ }
}
Matrix trkCov = trkCovNoNoise.plus(Q);
- System.out.println("KFSite: trkCov:");
- trkCov.print(15,12);
+ if(_cst.debug()) {
+ System.out.println("KFSite: trkCov:");
+ trkCov.print(15,12);
+ }
Matrix filteredCov = new Matrix(5,5);
@@ -227,25 +241,26 @@
vP[j] = (trkVarPlusPropagated.getPhi()-trkVarMinusPropagated.getPhi())/variations[j];
vZ[j] = (trkVarPlusPropagated.getZ()-trkVarMinusPropagated.getZ())/variations[j];
- //System.out.print("KFSite: vP = ");
- //System.out.print(vP[j]);
- //System.out.print(" ; vZ = ");
- //System.out.println(vZ[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();
- System.out.println("KFSite: H:");
- H.print(12,10);
- System.out.println("KFSite: HT:");
- HT.print(12,10);
+ 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();
- System.out.println("KFSite: trkCovInverse:");
- trkCovInverse.print(15,12);
+ if(_cst.debug()) {
+ System.out.println("KFSite: trkCovInverse:");
+ trkCovInverse.print(15,12);
+ }
double[] ptParam = trk.getPointParameters();
double refPhi = Math.atan2(ptParam[1],ptParam[0]);
@@ -254,21 +269,26 @@
//double[] mp = {_siteLocation.getPhi(),_siteLocation.getZ()};
Matrix measVect = new Matrix(mp,2);
- System.out.println("KFSite: measVect:");
- measVect.print(12,10);
-
- System.out.println("KFSite: PhiZ cov matrix inverted:");
- _PhiZErrorMatrixInverse.print(12,10);
+ 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);
- System.out.println("KFSite: measCovInverse:");
- cv2.print(15,12);
+ if(_cst.debug()) {
+ System.out.println("KFSite: measCovInverse:");
+ cv2.print(15,12);
+ }
Matrix filteredCovInv = trkCovInverse.plus(cv2);
filteredCov = filteredCovInv.inverse();
- System.out.println("KFSite: Filtered Cov:");
- filteredCov.print(15,12);
+ if(_cst.debug()) {
+ System.out.println("KFSite: Filtered Cov:");
+ filteredCov.print(15,12);
+ }
//Matrix tmp1 = trkCovInverse.times(stateVect);
//Matrix tmp2 = HT.times(_PhiZErrorMatrixInverse).times(measVect);
@@ -278,8 +298,10 @@
Matrix tmp3 = tmp1.plus(tmp2);
Matrix filteredTrkPar = filteredCov.times(tmp3);
- System.out.println("KFSite: Filtered State Vector:");
- filteredTrkPar.print(15,12);
+ if(_cst.debug()) {
+ System.out.println("KFSite: Filtered State Vector:");
+ filteredTrkPar.print(15,12);
+ }
double[][] par = filteredTrkPar.getArrayCopy();
d0 = par[0][0];
@@ -288,26 +310,40 @@
z0 = par[3][0];
tanl = par[4][0];
- System.out.print("KFSite: Fitted Pt = "); System.out.println(0.015/Math.abs(omega));
+ 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);
+ //System.out.print("KFSite: fitPar 1:"); tmpTrk.printFitParams();
tmpTrk.setPointToRadius(_siteLocation.getR());
+ //double tmpRad = _siteLocation.getR();
+ //System.out.println("KFSite: site radius ="+tmpRad);
+ //System.out.print("KFSite: fitPar 2:"); tmpTrk.printFitParams();
+ //System.out.print("KFSite: pointPar 2:"); tmpTrk.printPointParams();
double tmpPtot = tmpTrk.getPtot();
double ratioCorr = 1.;
- if(_applyDEDXCorrection) {
+ if(_materialSite && _applyDEDXCorrection) {
double dEdx = dEdx(tmpPtot);
- System.out.print("KFSite: Ptot = "); System.out.println(tmpPtot);
- System.out.print("KFSite: dEdx = "); System.out.println(dEdx);
-
- double tmpEtot = Math.sqrt(_cst.pionMass()*_cst.pionMass()+tmpPtot*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 -= _thickness*dEdx; // apply energy _loss_ when fitting inside-out
+ tmpEtotCorr -= dEdx*_thickness/cosEntranceAngle; // apply energy _loss_ when fitting inside-out
} else {
- tmpEtotCorr += _thickness*dEdx; // apply energy _gain_ when fitting outside-in
+ tmpEtotCorr += dEdx*_thickness/cosEntranceAngle; // apply energy _gain_ when fitting outside-in
}
ratioCorr = tmpEtotCorr/tmpEtot;
}
@@ -325,9 +361,6 @@
ratioCorr*tmpPar[4],
ratioCorr*tmpPar[5]);
- //KFTrackParameters filtTrk = new KFTrackParameters();
- //filtTrk.setTrackFitParameters(d0,phi0,omega,z0,tanl,l);
-
if(_measurementSite) {
filtTrk.setErrorMatrix(filteredCov);
} else {
@@ -356,11 +389,9 @@
public void applyMSCorrection(boolean bool) {_applyMSCorrection=bool;}
public void applyDEDXCorrection(boolean bool) {
- System.out.print("KFSite: Setting dedx bool from: ");
- System.out.print(_applyDEDXCorrection);
- System.out.print(" to ");
+ if(_cst.debug()) System.out.print("KFSite: Setting dedx bool from: "+_applyDEDXCorrection+" to ");
_applyDEDXCorrection=bool;
- System.out.println(_applyDEDXCorrection);
+ if(_cst.debug()) System.out.println(_applyDEDXCorrection);
}
private double dEdx(double p)
@@ -369,10 +400,10 @@
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.pionMass()*_cst.pionMass());
+ 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 = 0.000510998/_cst.pionMass();
+ double s = 0.000510998/_cst.muonMass();
double Wmax = 2.*0.000510998*eta*eta/(1.+2.*s*Math.sqrt(1+eta*eta)+s*s);
double I = 0.000000172;
lcsim/src/org/lcsim/contrib/KFFitter
diff -u -r1.3 -r1.4
--- KFTrack.java 1 Feb 2007 22:29:34 -0000 1.3
+++ KFTrack.java 10 Mar 2007 00:27:25 -0000 1.4
@@ -11,6 +11,14 @@
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.4 2007/03/10 00:27:25 fblanc Exp $
+ */
+
public class KFTrack
{
public KFTrack()
@@ -19,25 +27,31 @@
_orderedSiteList = new java.util.LinkedList();
}
- public void addHit(SimTrackerHit hit)
+ public void addHit(SimTrackerHit hit, double rdm1, double rdm2)
{
// define site for this hit
KFSite kfhit = new KFSite();
- kfhit.setLocation(hit);
+ kfhit.setLocation(hit,rdm1,rdm2);
+
kfhit.isMeasurementSite(true);
double radius = kfhit.getSiteLocation().getR();
- double radThickness = 0.00107;
- if (radius>10.) {radThickness = 0.00532;}
- double thickness = 0.01;
- if (radius>10.) {thickness = 0.03;}
+ 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);
- //System.out.print("KFTrack: read back radius = ");
- //System.out.println(kfhit.getSiteLocation().getR());
}
@@ -79,19 +93,24 @@
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++) {
@@ -100,9 +119,22 @@
nextSmallestDist=dist[i];
}
}
- _orderedSiteList.add(_siteList.get(index));
+ 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)
@@ -116,18 +148,13 @@
double l = 0.0; // path length set to zero (=POCA)
- System.out.print("KFTrack: define track par = ");
- System.out.print(d); System.out.print(" ");
- System.out.print(phi); System.out.print(" ");
- System.out.print(omega); System.out.print(" ");
- System.out.print(z); System.out.print(" ");
- System.out.println(s);
+ 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 defineSeedTrack()
+ public void defineSeedTrackDebug()
{
// ************************
// For testing purpose only
@@ -141,19 +168,14 @@
double d = 0.;
double phi = 0;
//double omega = 0.00300; // 5GeV
- //double omega = 0.00075; // 20GeV
- double omega = 0.00015; // 100GeV
+ double omega = 0.00075; // 20GeV
+ //double omega = 0.00015; // 100GeV
double z = 0.;
double s = 0.;
double l = 0.0; // path length set to zero (=POCA)
- System.out.print("KFTrack: define track par = ");
- System.out.print(d); System.out.print(" ");
- System.out.print(phi); System.out.print(" ");
- System.out.print(omega); System.out.print(" ");
- System.out.print(z); System.out.print(" ");
- System.out.println(s);
+ System.out.println("KFTrack: define track par = "+d+" "+phi+" "+omega+" "+z+" "+s);
_seedTrack.setTrackFitParameters(d,phi,omega,z,s,l);
@@ -164,60 +186,51 @@
// sort list of hits by distance from origin
sortSiteList();
- int initChg = _seedTrack.getCharge();
- // initialize track param correlation matrix
- double[][] cor = new double[5][5];
- cor[0][0] = 1.; // d0
- cor[0][1] = cor[1][0] = initChg*0.6; // d0 - phi0
- cor[0][2] = cor[2][0] = -initChg*0.3; // d0 - omega
- cor[1][1] = 1.; // phi0
- cor[1][2] = cor[2][1] = -0.8; // phi0 - omega
- cor[2][2] = 1.; // omega
- cor[3][3] = 1.; // z0
- cor[3][4] = cor[4][3] = -initChg*0.9; // z0 - tanl
- cor[4][4] = 1.; // tanl
- Matrix initCor = new Matrix(cor);
-
// initialize track param covariance matrix
double[][] diag = new double[5][5];
- diag[0][0] = 10.; // d0
- diag[1][1] = 1.; // phi0
- diag[2][2] = 100.; // omega
- diag[3][3] = 100.; // z0
- diag[4][4] = 100.; // tanl
- Matrix initDiag = new Matrix(diag);
-
- //Matrix initCov = new Matrix();
- Matrix initCov = initDiag.times(initCor).times(initDiag);
- System.out.println("KFTrack: init cov matrix");
- initCov.print(6,2);
-
- System.out.print("KFTrack: ");
- _seedTrack.printFitParams();
+ 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);
+
+ //Matrix initCov = initDiag.times(initCor).times(initDiag);
+ 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
//***************
- System.out.println("KFTrack: Fitting Track outward");
+ if(_cst.debug()) System.out.println("KFTrack: Fitting Track outward");
// 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));
+ KFTrackParameters currentTrackOut = new KFTrackParameters(inOutSeedTrk);
+ */
// initialize track parameters for outward fit
- //KFTrackParameters currentTrackOut = new KFTrackParameters(_seedTrack);
- KFTrackParameters currentTrackOut = new KFTrackParameters(inOutSeedTrk);
+ KFTrackParameters currentTrackOut = new KFTrackParameters(_seedTrack);
+
// determine track params at fixed small radius
//_seedTrack.setPointToPOCA();
double refRadius = 1.0;
currentTrackOut.setPointToRadius(refRadius);
- System.out.print("KFTrack: set radius to 1 : ");
- currentTrackOut.printFitParams();
+ if(_cst.debug()) {
+ System.out.print("KFTrack: set radius to 1 : ");
+ currentTrackOut.printFitParams();
+ }
// initialize error matrix
currentTrackOut.setErrorMatrix(initCov);
@@ -228,17 +241,21 @@
//siteLocation = hit.getSiteLocation();
//double R = siteLocation.getR();
double R = hit.getSiteLocation().getR();
- System.out.print("KFTrack: hit radius = "); System.out.println(R);
+ if(_cst.debug()) System.out.println("KFTrack: hit radius = "+R);
- System.out.print("KFTrack: current pars = ");
- currentTrackOut.printFitParams();
+ if(_cst.debug()) {
+ System.out.print("KFTrack: current pars = ");
+ currentTrackOut.printFitParams();
+ }
// propagate track to radius of hit
KFTrackParameters predictedTrack = currentTrackOut.makePredictionAtRadius(R);
- System.out.print("KFTrack: predicted pars = ");
- predictedTrack.printFitParams();
- predictedTrack.printPointParams();
+ if(_cst.debug()) {
+ System.out.print("KFTrack: predicted pars = ");
+ predictedTrack.printFitParams();
+ predictedTrack.printPointParams();
+ }
// filter predicted track parameters with hit at radius R
hit.setFitOutward();
@@ -251,24 +268,18 @@
//***************
// Fit outside-in
//***************
- System.out.println("KFTrack: Fitting Track inward");
-
- // Determine seed track parameters from 3 outtermost hits
- //int hitListSize = _orderedSiteList.size();
- //KFTrackParameters outInSeedTrk =
- // seedTrkFromHits(_seedTrack.getCharge(),
- // (KFSite)_orderedSiteList.get(hitListSize-1),
- // (KFSite)_orderedSiteList.get(hitListSize-2),
- // (KFSite)_orderedSiteList.get(hitListSize-3));
- //KFTrackParameters currentTrackIn = new KFTrackParameters(outInSeedTrk);
+ 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);
- System.out.print("KFTrack: set radius to 130 : ");
- currentTrackIn.printFitParams();
+ if(_cst.debug()) {
+ System.out.print("KFTrack: set radius to 130 : ");
+ currentTrackIn.printFitParams();
+ }
+
// initialize error matrix
currentTrackIn.setErrorMatrix(initCov);
@@ -279,16 +290,19 @@
//siteLocation = hit.getSiteLocation();
//double R = siteLocation.getR();
double R = hit.getSiteLocation().getR();
- System.out.print("KFTrack: hit radius = "); System.out.println(R);
-
- System.out.print("KFTrack: current pars = ");
- currentTrackIn.printFitParams();
+ 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);
- System.out.print("KFTrack: predicted pars = ");
- predictedTrack.printFitParams();
+ if(_cst.debug()) {
+ System.out.print("KFTrack: predicted pars = ");
+ predictedTrack.printFitParams();
+ }
// filter predicted track parameters with hit at radius R
hit.setFitInward();
@@ -299,25 +313,29 @@
}
+ _KFTrackAtPOCA = currentTrackIn;
+
_fittedPtAtPOCA = currentTrackIn.getPt();
- System.out.print("KFTrack: Final fitted pt = ");
- System.out.println(_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];
+
+ 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);
}
- Matrix corMat = new Matrix(covArray);
- System.out.print("KFTrack: Final corrlation matrix = ");
- corMat.print(15,12);
}
private KFTrackParameters seedTrkFromHits(int q, KFSite p1, KFSite p2, KFSite p3)
@@ -334,19 +352,6 @@
double r1 = Math.sqrt(x1*x1+y1*y1);
double r2 = Math.sqrt(x2*x2+y2*y2);
- System.out.println("KFTrack: hit coordinates = ");
- System.out.print(x1); System.out.print(" ");
- System.out.print(y1); System.out.print(" ");
- System.out.print(z1); System.out.print(" ");
- System.out.println(r1);
- System.out.print(x2); System.out.print(" ");
- System.out.print(y2); System.out.print(" ");
- System.out.print(z2); System.out.print(" ");
- System.out.println(r2);
- System.out.print(x3); System.out.print(" ");
- System.out.print(y3); System.out.print(" ");
- System.out.println(Math.sqrt(x3*x3+y3*y3));
-
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;
@@ -372,13 +377,17 @@
return seedTrk;
}
- public KFTrackParameters getTrackParameters() { return _seedTrack; }
+ 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 KFConstants _cst = new KFConstants();
+
}
lcsim/src/org/lcsim/contrib/KFFitter
diff -u -r1.3 -r1.4
--- KFTrackParameters.java 1 Feb 2007 22:29:34 -0000 1.3
+++ KFTrackParameters.java 10 Mar 2007 00:27:25 -0000 1.4
@@ -2,6 +2,13 @@
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.4 2007/03/10 00:27:25 fblanc Exp $
+ */
public class KFTrackParameters
{
@@ -36,7 +43,7 @@
double zzero, double spar, double pathLength)
{
_charge = (int) (-1*Math.signum(omeg));
- //System.out.print("charge = "); System.out.println(_charge);
+ //System.out.println("charge = "+_charge);
_d0 = dzero;
_phi0 = phizero;
_omega = omeg; if(_omega==0.) {System.out.println("Null curvature!");}
@@ -56,8 +63,9 @@
//this.printFitParams();
makePointFromFitRepresentation();
- makeHelixFromFitRepresentation();
+ //makeHelixFromFitRepresentation();
+ //makeFitFromPointRepresentation();
//System.out.print("KFTrackParameters: control fit pars = ");
//this.printFitParams();
@@ -83,7 +91,7 @@
_pointParams[5] = _pz;
makeFitFromPointRepresentation();
- makeHelixFromPointRepresentation();
+ //makeHelixFromPointRepresentation();
//System.out.print("KFTrackParameters: control fit pars = ");
//this.printFitParams();
@@ -102,7 +110,8 @@
(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+_alpha*_s/Math.abs(_omega);
+ //_z = _z0+_s*Math.abs(_alpha/_omega);
+ _z = _z0+_s*_l;
_pointParams[0] = _x;
_pointParams[1] = _y;
@@ -129,10 +138,10 @@
_helixParams[4] = _s;
//System.out.print("KFTrackParameters: helix par = ");
- //System.out.print(_xc); System.out.print(" ");
- //System.out.print(_yc); System.out.print(" ");
- //System.out.print(_rc); System.out.print(" ");
- //System.out.print(_z0); System.out.print(" ");
+ //System.out.print(_xc+" ");
+ //System.out.print(_yc+" ");
+ //System.out.print(_rc+" ");
+ //System.out.print(_z0+" ");
//System.out.println(_s);
}
@@ -153,10 +162,10 @@
_helixParams[4] = _s;
//System.out.print("KFTrackParameters: helix par = ");
- //System.out.print(_xc); System.out.print(" ");
- //System.out.print(_yc); System.out.print(" ");
- //System.out.print(_rc); System.out.print(" ");
- //System.out.print(_z0); System.out.print(" ");
+ //System.out.print(_xc+" ");
+ //System.out.print(_yc+" ");
+ //System.out.print(_rc+" ");
+ //System.out.print(_z0+" ");
//System.out.println(_s);
}
@@ -168,8 +177,15 @@
_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(_y*_py+_x*_px,_charge*tmpPt*tmpPt/_cst.pointThreeTimesB()+_x*_py-_y*_px);
- _z0 = _z-_alpha*_pz/_cst.pointThreeTimesB();
+ _alpha = Math.atan2(-_charge*(_x*_px+_y*_py),
+ tmpPt*tmpPt/_cst.pointThreeTimesB()+_charge*(_x*_py-_y*_px));
+ _z0 = _z+_alpha*_pz/(_charge*_cst.pointThreeTimesB());
+ //double t1=_y*_py+_x*_px;
+ //double t2=_charge*tmpPt*tmpPt/_cst.pointThreeTimesB()+_x*_py-_y*_px;
+ //System.out.println("Trk: atan("+t1+" ; "+t2);
+ //System.out.println("Trk: x= "+_x+" "+_y+" "+_z+" "+_px+" "+_py+" "+_pz);
+ //System.out.println("Trk: "+tmpPt+" "+_charge+" "+_cst.pointThreeTimesB());
+ //System.out.println("Trk: alpha = "+_alpha+" z0 = "+_z0);
_s = _pz/tmpPt;
_l = _alpha/_omega;
@@ -221,14 +237,15 @@
public void setPointToRadius(double r) {
double newPathLength = propagate(r);
_l = newPathLength;
+ _alpha = _l*_omega;
makePointFromFitRepresentation();
- makeHelixFromFitRepresentation();
+ //makeHelixFromFitRepresentation();
}
public void setPointToPOCA() {
_l = 0.;
makePointFromFitRepresentation();
- makeHelixFromFitRepresentation();
+ //makeHelixFromFitRepresentation();
}
public double[] getFitParameters() {return _fitParams;}
@@ -268,24 +285,24 @@
public void printFitParams() {
System.out.print("Fit params = ");
- System.out.print(_fitParams[0]); System.out.print(" ");
- System.out.print(_fitParams[1]); System.out.print(" ");
- System.out.print(_fitParams[2]); System.out.print(" ");
- System.out.print(_fitParams[3]); System.out.print(" ");
- System.out.print(_fitParams[4]); System.out.print(" ");
- //System.out.print(_alpha); System.out.print(" ");
+ 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.print(_alpha);
System.out.println(_l);
}
public void printPointParams() {
System.out.print("Point params = ");
- System.out.print(_pointParams[0]); System.out.print(" ");
- System.out.print(_pointParams[1]); System.out.print(" ");
- System.out.print(_pointParams[2]); System.out.print(" ");
- System.out.print(_pointParams[3]); System.out.print(" ");
- System.out.print(_pointParams[4]); System.out.print(" ");
- System.out.print(_pointParams[5]); System.out.print(" ");
- //System.out.print(_alpha); System.out.print(" ");
+ 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.print(_alpha);
System.out.println(getPhi());
}
}
lcsim/src/org/lcsim/contrib/KFFitter/test
diff -u -r1.1 -r1.2
--- TestKFF.java 13 Dec 2006 22:45:12 -0000 1.1
+++ TestKFF.java 10 Mar 2007 00:27:26 -0000 1.2
@@ -2,11 +2,14 @@
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;
@@ -26,5 +29,35 @@
};
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