Print

Print


Commit in lcsim/sandbox/RobKutschke/TRFTests/v1 on MAIN
FitTest_sid00Driver.java+95-371.1 -> 1.2
New diagnostics.  Fix starting covariance matrix.  Test RKCovSmear.

lcsim/sandbox/RobKutschke/TRFTests/v1
FitTest_sid00Driver.java 1.1 -> 1.2
diff -u -r1.1 -r1.2
--- FitTest_sid00Driver.java	30 Aug 2007 14:16:33 -0000	1.1
+++ FitTest_sid00Driver.java	7 Sep 2007 19:04:15 -0000	1.2
@@ -29,9 +29,9 @@
  * no multiple scattering or energy loss.
  *
  *@author $Author: kutschke $
- *@version $Id: FitTest_sid00Driver.java,v 1.1 2007/08/30 14:16:33 kutschke Exp $
+ *@version $Id: FitTest_sid00Driver.java,v 1.2 2007/09/07 19:04:15 kutschke Exp $
  *
- * Date $Date: 2007/08/30 14:16:33 $
+ * Date $Date: 2007/09/07 19:04:15 $
  *
  */
 
@@ -49,12 +49,15 @@
     RKTrackGenParams param = null;
     RKTrackGen gen = null;
 
-    private TrackError vstartc = new TrackError();
-    private TrackError vstartz = new TrackError();
+    private TrackError vstartc_in  = new TrackError();
+    private TrackError vstartz_in  = new TrackError();
+    private TrackError vstart_out  = new TrackError();
 
     private RKTrackFitDiag DiagForward  = null;
     private RKTrackFitDiag DiagBackward = null;
 
+    private RKCovSmear covsmear = null;
+
     public FitTest_sid00Driver(){
     }
 
@@ -66,7 +69,7 @@
 	// Limits for track generation parameters.
 	param = new RKTrackGenParams();
 	param.setbz(rkgeom.getBz());
-	param.setptmin(  10.);
+	param.setptmin(  10.0);
 	//param.setczmin(  0.93);
 	//param.setczmax(  0.993);
 	param.setczmin(  -0.993);
@@ -88,19 +91,34 @@
 	seed += 876633217;
 	RKHit.setSeed(seed);
 
-	// Intial covariance matrix.
-	double sig=10.0;
-	vstartc.set(0,0,sig);
-	vstartc.set(1,1,500.);
-	vstartc.set(2,2,sig);
-	vstartc.set(3,3,sig);
-	vstartc.set(4,4,sig);
-
-	vstartz.set(0,0,sig);
-	vstartz.set(1,1,sig);
-	vstartz.set(2,2,sig);
-	vstartz.set(3,3,sig);
-	vstartz.set(4,4,sig);
+	// Seed for smearing covariance matrices.
+	seed += 876633217;
+	covsmear = new RKCovSmear( seed );
+
+	// Intial covariance matrices.
+	double sig=0.0625;
+	double fac1=10.;
+	double fac2=100.;
+	vstartc_in.set(0,0,sig/fac1);
+	vstartc_in.set(1,1,sig/fac1);
+	vstartc_in.set(2,2,sig/fac2);
+	vstartc_in.set(3,3,sig/fac2);
+	vstartc_in.set(4,4,sig/fac2);
+
+	vstartz_in.set(0,0,sig);
+	vstartz_in.set(1,1,sig);
+	vstartz_in.set(2,2,sig);
+	vstartz_in.set(3,3,sig);
+	vstartz_in.set(4,4,sig);
+
+	sig=1.0;
+	vstart_out.set(0,0,sig);
+	vstart_out.set(1,1,sig);
+	vstart_out.set(2,2,sig);
+	vstart_out.set(3,3,sig);
+	vstart_out.set(4,4,sig);
+
+
 
     }
 
@@ -112,8 +130,9 @@
 			   + event.getDetector().getName() );
 
 	//int ndo = 1000;
-	int ndo = 1000;
+	int ndo = 10000;
 	for ( int i=0; i<ndo; ++i ){
+	    RKDebug.Instance().setTrack(i);
 	    DoIt(i);
 	}
 
@@ -122,7 +141,7 @@
 
     private void DoIt( int itrk){
 
-	System.out.println ("New track: " + itrk );
+	//	System.out.println ("New track: " + itrk );
 
 	// Generate a track.
 	RKTrack rktrk = gen.newRKTrack();
@@ -130,6 +149,9 @@
 	// Monitor generated track info.
 	plotGenTrk(rktrk);
 
+	// Load track info into the debug-info container.
+	RKDebug.Instance().setRKTrack(rktrk);
+
 	// Generate hits on this track in my own format.
 	RKHitList hlist = GenerateOneArcHits( rktrk );
 
@@ -147,7 +169,8 @@
 	toTRF trftrk = new toTRF(rktrk);
 	VTrack vtdca = trftrk.atDcaZaxis();
 
-	TrackError vstart = vstartz;
+	TrackError vstart = vstart_out;
+	RKDebug.Instance().setStartType(0);
 
 	// Prep for the outward fit.
 	ETrack et = new ETrack( vtdca.surface().newPureSurface(), vtdca.vector(), vstart );
@@ -160,6 +183,14 @@
 	// Do the outward fit.
 	//System.out.println ("Start forward fit.");
 	int fstatf = fitk.fitForward(ht);
+	if ( fstatf  != 0 ){
+	    System.out.println("Forwards status: " + 
+			       fstatf + " " + 
+			       hlist.nHits() + " " + 
+			       hlist.nCyl() + " " + 
+			       hlist.nZp() + " " + 
+			       rktrk.cz() );
+	}
 
 	// Final state reference parameters for the outward fit.
 	VTrack vtend = hlist.outerMostHit().getVTrack();
@@ -170,14 +201,32 @@
 	    double[] limits = { 0.00005, 0.15, 0.0003, 0.001, 0.0003};
 	    DiagForward = new RKTrackFitDiag ( "Forward", "Cyl", nbins, limits );
 	}
-	DiagForward.fill( fstatf, ht, vtend, hlist.nDof(), vstart);
+	DiagForward.fill( fstatf, ht, vtend, hlist.nDof(), vstart, rktrk);
 
+	// Check code for smearing tracks with covariance matrix.
+	TrackError te = ht.newTrack().error();
+	TrackVector vv = covsmear.Smear( ht.newTrack().error(), ht.newTrack().vector() );
+	for ( int i=0; i<5; ++i ){
+	    double xx = te.get(i,i);
+	    double sig = Math.sqrt(Math.abs(xx));
+	    double pull = (vv.get(i)-ht.newTrack().vector(i))/sig;
+	    aida.histogram1D("Smear Check",100, -5., 5.).fill(pull);
+	}
+
+	
 	if ( vtend.surface().pureType().equals(SurfCylinder.staticType()) ){
-	    vstart = vstartc;
+	    vstart = vstartc_in;
+	    RKDebug.Instance().setStartType(1);
+	} else{
+	    vstart = vstartz_in;
+	    RKDebug.Instance().setStartType(2);
 	}
 
 	// Prep for inward fit.
-	ETrack etin = new ETrack( vtend.surface().newPureSurface(), vtend.vector(), vstart );
+	ETrack etin = new ETrack( vtend.surface().newPureSurface(), 
+				  ht.newTrack().vector(),
+				  vstart );
+				  //vtend.vector(), vstart );
 
 	// Not sure why this is needed.
 	if ( vtend.surface().type() == SurfZPlane.staticType() ){
@@ -204,6 +253,11 @@
 			       hlist.nCyl() + " " + 
 			       hlist.nZp() + " " + 
 			       rktrk.cz() );
+	    if ( fstatb < 100 ){
+		aida.histogram1D("Backward: cz for failed propagation",50, -1., 1. ).fill(rktrk.cz());
+	    } else{
+		aida.histogram1D("Backward: cz for failed addhit",50, -1., 1. ).fill(rktrk.cz());
+	    }
 	}
 	if ( fstatb != 0 ){
 	    System.out.println ("Bad fit: "
@@ -222,7 +276,11 @@
 				    h.getPsi()
 				    );
 	    }
-
+	    aida.cloud2D("Bad fit nZ vs nCyl").fill( hlist.nCyl(), hlist.nZp() );
+	    aida.histogram1D( "Bad fit: Start type",5,0.,5.).fill(RKDebug.Instance().getStartType());
+	} else{
+	    aida.cloud2D("Good fit nZ vs nCyl").fill( hlist.nCyl(), hlist.nZp() );
+	    aida.histogram1D( "Good fit: Start type",5,0.,5.).fill(RKDebug.Instance().getStartType());
 	}
 
 	// Diagnostics for inward fit.
@@ -231,7 +289,7 @@
 	    double[] limits = { 0.002, 0.002, 0.0001, 0.0001, 0.0005};
 	    DiagBackward = new RKTrackFitDiag ( "Backward", "DCA", nbins, limits );
 	}
-	DiagBackward.fill( fstatb, htin, vtdca, hlist.nDof(), vstart);
+	DiagBackward.fill( fstatb, htin, vtdca, hlist.nDof(), vstart, rktrk);
 
     }
 
@@ -260,8 +318,8 @@
 		if ( sum > pathlimit ) break;
 		if ( s.inBounds(vtdca.vector(1)) ){
 		    rkl.addHit( new RKHit( rktrk, s, vtdca, ps, sum ) );
-		    aida.cloud2D( "C r vs z").fill( vtdca.vector(1), s.radius*mmTocm );
-		    aida.cloud2D( "Both r vs z").fill( vtdca.vector(1), s.radius*mmTocm );
+		    aida.cloud2D( "C r vs z",-1).fill( vtdca.vector(1), s.radius*mmTocm );
+		    aida.cloud2D( "Both r vs z",-1).fill( vtdca.vector(1), s.radius*mmTocm );
 		    /*
 		    double r = s.radius*mmTocm;
 		    System.out.printf ("Adding Cyl: %10.4f %10.4f %10.4f %10.4f %10.4f %10.4f\n",
@@ -288,9 +346,9 @@
 	} else{
 	    System.out.println ("Can't get here ...");
 	}
-	aida.cloud2D("Path length vs cz").fill(rktrk.cz(),xsum);
+	aida.cloud2D("Path length vs cz",-1).fill(rktrk.cz(),xsum);
 	aida.cloud1D("Sum").fill(xsum);
-	aida.cloud2D("Hits vs cz").fill( rktrk.cz(), ncmeas );
+	aida.cloud2D("Hits vs cz",-1).fill( rktrk.cz(), ncmeas );
 
 	sum = 0;
 	for ( RKSurf s: rkgeom.getZ( rktrk.z0(), rktrk.cz() ) ){
@@ -301,8 +359,8 @@
 		double r = Math.sqrt(vtz.vector(0)*vtz.vector(0) + vtz.vector(1)*vtz.vector(1) );
 		if ( s.inBounds(r) ){
 		    rkl.addHit( new RKHit( rktrk, s, vtz, ps, sum ) );
-		    aida.cloud2D( "Z r vs z").fill( s.zc*mmTocm, r );
-		    aida.cloud2D( "Both r vs z").fill( s.zc*mmTocm, r );
+		    aida.cloud2D( "Z r vs z",-1).fill( s.zc*mmTocm, r );
+		    aida.cloud2D( "Both r vs z",-1).fill( s.zc*mmTocm, r );
 		    nmeas += s.rodim;
 		    nzmeas += s.rodim;
 		    /*
@@ -316,9 +374,9 @@
 		}
 	    }	    
 	}
-	aida.cloud2D("Z Hits vs cz").fill( rktrk.cz(), nzmeas );
-	aida.cloud2D("Nmeas vs cz").fill( rktrk.cz(), nmeas );
-	aida.cloud2D("Nz vs Nc").fill( ncmeas, nzmeas );
+	aida.cloud2D("Z Hits vs cz",-1).fill( rktrk.cz(), nzmeas );
+	aida.cloud2D("Nmeas vs cz",-1).fill( rktrk.cz(), nmeas );
+	aida.cloud2D("Nz vs Nc",-1).fill( ncmeas, nzmeas );
 	return rkl;
     }
 
@@ -339,8 +397,8 @@
 	aida.histogram1D("d0", nbins, d0min, d0max).fill(t.d0());
 	aida.histogram1D("z0", nbins, z0min, z0max).fill(t.z0());
 	aida.histogram1D("p", nbins, 0., 50.).fill(t.p());
-	aida.cloud2D("PCA0").fill(t.x0(),t.y0());
-	aida.cloud2D("D0 vs Z0").fill(t.z0(),t.d0());
+	aida.cloud2D("PCA0",-1).fill(t.x0(),t.y0());
+	aida.cloud2D("D0 vs Z0",-1).fill(t.z0(),t.d0());
     }
 
 }
CVSspam 0.2.8