lcsim/sandbox/RobKutschke/TRFTests/v1

FitTest_sid00Driver.java 1.3 -> 1.4

diff -u -r1.3 -r1.4
--- FitTest_sid00Driver.java	13 Sep 2007 22:36:02 -0000	1.3
+++ FitTest_sid00Driver.java	18 Sep 2007 16:50:17 -0000	1.4
@@ -21,17 +21,22 @@

 import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
 import org.lcsim.recon.tracking.trfzp.SurfZPlane;

+import org.lcsim.recon.tracking.trfcyl.ThinCylMs;
+import org.lcsim.recon.tracking.trfcyl.ThinCylMsSim;
+

 import org.lcsim.recon.tracking.trffit.HTrack;

+import hep.physics.vec.Hep3Vector;
+

 /**
  * Driver for exercises learning how to use the TRF code.
  * This version exercises the fitter with only measurement errors,
  * no multiple scattering or energy loss.
  *
  *@author $Author: kutschke $

- *@version $Id: FitTest_sid00Driver.java,v 1.3 2007/09/13 22:36:02 kutschke Exp $

+ *@version $Id: FitTest_sid00Driver.java,v 1.4 2007/09/18 16:50:17 kutschke Exp $

- * Date $Date: 2007/09/13 22:36:02 $

+ * Date $Date: 2007/09/18 16:50:17 $

  *
  */

@@ -49,15 +54,24 @@

     RKTrackGenParams param = null;
     RKTrackGen gen = null;

-    private TrackError vstartc_in  = new TrackError();
-    private TrackError vstartz_in  = new TrackError();
-    private TrackError vstart_out  = 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;

+    private ThinCylMsSim scatSim0 = null;
+    private ThinCylMsSim scatSim1 = null;
+    private ThinCylMsSim scatSim2 = null;
+    private double scatThick0 = 0.006136;
+    private double scatThick1 = 0.000916;
+    private double scatThick2 = 0.013747;
+
+    private boolean doms = true;
+

     public FitTest_sid00Driver(){
     }

@@ -69,18 +83,28 @@

 	// Limits for track generation parameters.
 	param = new RKTrackGenParams();
 	param.setbz(rkgeom.getBz());

-	param.setptmin(  0.5);

+	param.setptmin(  1.0 );
+	//param.setptmin(  20.0);
+	//param.setptmax(  20.0);

 	//param.setczmin(  0.93);
 	//param.setczmax(  0.993);

-	param.setczmin(  -0.993);
-	param.setczmax(   0.993);

+	//param.setczmin(  -0.993);
+	//param.setczmax(   0.993);
+	//param.setczmin(  0.);
+	//param.setczmax(   0.);

 	//param.setczmax(  -0.93);
 	//param.setczmin(  -0.97);
 	//param.setczmax(  -0.93);

-	//param.setczmin(  -0.4);
-	//param.setczmax(   0.4);

+	param.setczmin(  -0.6);
+	param.setczmax(   0.6);

 	//param.setczmin( -0.995);
 	//param.setczmax(  0.995);

+
+	//param.setd0min(  0.0);
+	//param.setd0max(  0.0);
+	//param.setz0min(  0.0);
+	//param.setz0max(  0.0);
+

 	System.out.println (param);
 
 	// Seed for track random number generator

@@ -118,6 +142,32 @@

 	vstart_out.set(3,3,sig);
 	vstart_out.set(4,4,sig);

+	RKDebug.Instance().setDoMs(doms);
+	RKDebug.Instance().setMsFac(1.0);
+
+	scatSim0 = new ThinCylMsSim(scatThick0*RKDebug.Instance().getMsFac());
+	scatSim1 = new ThinCylMsSim(scatThick1*RKDebug.Instance().getMsFac());
+	scatSim2 = new ThinCylMsSim(scatThick2*RKDebug.Instance().getMsFac());
+
+	// Cylinder track: Pt= 1 GeV, phat=(1,0.,0.) as it goes through some cylinder.
+	SurfCylinder stmp = new SurfCylinder(1.0);
+	TrackVector vtmp = new TrackVector();
+	vtmp.set(0,0.0);
+	vtmp.set(1,0.0);
+	vtmp.set(2,0.0);
+	vtmp.set(3,0.0);
+	vtmp.set(3,1.0);
+	VTrack ttmp = new VTrack( stmp, vtmp );
+	
+	// Desynchronize the random number generators.
+	//   Want moderate size prime numbers to avoid any of the
+	//   periodicities that might occur in the detector.
+	for ( int i=0; i<101; ++i ){
+	    scatSim0.interact(ttmp);
+	}
+	for ( int i=0; i<211; ++i ){
+	    scatSim2.interact(ttmp);
+	}

@@ -156,7 +206,8 @@

 	RKHitList hlist  = GenerateMixedOneArcHits( rktrk );
 
 	// Diagnostics for the generated hits.

-	CheckHitList( hlist, rktrk );

+	// Only useful when there is not scattering or eloss.
+	//CheckHitList( hlist, rktrk );

 
 	// Check for sufficient hits.
 	if ( ! hlist.fitable() ){

@@ -173,6 +224,7 @@

 	VTrack vtdca = trftrk.atDcaZaxis();
 
 	TrackError vstart = vstart_out;

 	RKDebug.Instance().setStartType(0);
 
 	// Prep for the outward fit.

@@ -206,6 +258,7 @@

 	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 ){

@@ -214,6 +267,7 @@

 	    double pull = (vv.get(i)-ht.newTrack().vector(i))/sig;
 	    aida.histogram1D("Smear Check",100, -5., 5.).fill(pull);
 	}

+	*/

 
 	// Prep for the inward fit.
 	if ( vtend.surface().pureType().equals(SurfCylinder.staticType()) ){

@@ -380,6 +434,7 @@

 	toTRF trftrk = new toTRF(rktrk);
 	VTrack vtdca = trftrk.atDcaZaxis();

 	// This code only looks at the first arc of a track.
 	double pathlimit = trftrk.halfarc();

@@ -438,7 +493,6 @@

 		rks_cyl = cyls.get(next_cyl);
 		ps_cyl  = prop.vecDirProp( vt_cyl, rks_cyl.getCylinder(), PropDir.FORWARD);

 		// If we do not get to this cylinder, try the next one.
 		// This can happen if starting point is outside this one but inside a later one.
 		if ( !ps_cyl.success() ) {

@@ -567,7 +621,7 @@

 		if ( sum > pathlimit ) break;
 
 		// Track parameters at this step.

-		vtdca = vt_cyl;

+		vtdca = new VTrack(vt_cyl);

 
 		// Add the hit.
 		rkl.addHit( new RKHit( rktrk, rks_cyl, vtdca, ps_cyl, sum ) );

@@ -580,6 +634,84 @@

 		nmeas  += rks_cyl.rodim;
 		ncmeas += rks_cyl.rodim;

+		if ( r_cyl > 1.3 && doms ){
+
+		    aida.histogram1D("Gen scat radius", 300, 0., 150.).fill(r_cyl);
+
+		    // Interact.
+		    VTrack test = new VTrack(vtdca);
+		    //rks_cyl.getCylinder().simInteract(vtdca);
+		    if ( r_cyl < 10. ){
+			scatSim1.interact(vtdca);
+		    }else{
+			scatSim2.interact(vtdca);
+		    }
+
+		    /*
+		    // Diagnostics.
+		    for ( int i=0;i<5; ++i){
+			double d = test.vector(i) - vtdca.vector(i);
+			aida.cloud1D("Thin MS for param " +i ).fill(d);
+		    }
+		    VTUtil vnew = new VTUtil(vtdca);
+		    VTUtil vold = new VTUtil(test);
+		
+		    double d = vnew.momentum()-vold.momentum();
+		    aida.cloud1D("Thin MS for p").fill(d);
+	
+		    Hep3Vector v1 = vnew.asHep3Vector();
+		    Hep3Vector v2 = vold.asHep3Vector();
+		
+		    double dot = v1.x()*v2.x() + v1.y()*v2.y() + v1.z()*v2.z();
+		    dot = dot/v1.magnitude()/v2.magnitude()-1.;
+		    aida.cloud1D("Thin MS for dot").fill(dot);
+		    */
+
+		    //if ( r_cyl <1.3 ){
+		    if ( r_cyl > 10. ){
+			double scatThick = (r_cyl<10.) ? scatThick1 : scatThick2;
+			ThinCylMs scat1 = new ThinCylMs(scatThick);
+			TrackError vtmp = new TrackError();
+			ETrack et = new ETrack( test.surface().newPureSurface(), test.vector(), vtmp);
+			ETrack ets = new ETrack(et);
+			scat1.interact(et);
+
+			/*
+			for ( int i=0;i<5; ++i){
+			    if ( et.error(i,i) >0. ){
+				double sig = Math.sqrt(et.error(i,i));
+				double pull = (test.vector(i) - vtdca.vector(i))/sig;
+				if ( rktrk.q() > 0. ){
+				    aida.histogram1D("Gen Pull param: q=+1 " + i,200,-10.,10.).fill(pull);
+				}else {
+				    aida.histogram1D("Gen Pull param: q=-1 " + i,200,-10.,10.).fill(pull);
+				}
+				if ( i == 4 ){
+				    if ( rktrk.q() > 0. ){
+					aida.histogram1D("Gen Pull curve: q=+1 r= " + r_cyl,200,-10.,10.).fill(pull);
+					aida.cloud1D("Gen Sig(4) for q=+1 layer " + r_cyl).fill(sig);
+				    }else {
+					aida.histogram1D("Gen Pull curve: q=-1 r= " + r_cyl,200,-10.,10.).fill(pull);
+					aida.cloud1D("Gen Sig(4) for q=-1 layer " + r_cyl).fill(sig);
+				    }
+				}
+			    }
+			
+			}
+			*/
+			
+			// This occurs for d0=0, z0=0, tanlambada=0 and any phi or curvature.
+			// I am not sure if this is a problem or not but I think it is.
+			if ( et.error(4,4) == 0. && et.error(4,3)==0. ){
+			    System.out.println ("Scatter matrix at r: " + r_cyl + "\n" + et.error() );
+			}
+
+		    }
+
+		    
+
+		}
+

 	    } else if ( addzp ) {
 
 		// Advance counter to next surface.

@@ -629,9 +761,17 @@

+    private String GenDir = null;

     // Monitoring histograms for the generated track.
     private void plotGenTrk( RKTrack t ){

+	String pwd = aida.tree().pwd();
+	if ( GenDir == null ){
+	    GenDir = "GeneratorInfo";
+	    aida.tree().mkdir(GenDir);
+	}
+	aida.tree().cd(GenDir);
+

 	int nbins = 50;
 
 	double d0min = -5.;

@@ -648,6 +788,8 @@

 	aida.histogram1D("p", nbins, 0., 50.).fill(t.p());
 	aida.cloud2D("PCA0",-1).fill(t.x0(),t.y0());
 	aida.cloud2D("D0 vs Z0",-1).fill(t.z0(),t.d0());

+
+	aida.tree().cd(pwd);

Commit in `lcsim/sandbox/RobKutschke/TRFTests/v1` on MAIN
`FitTest_sid00Driver.java`	+155	-13	1.3 -> 1.4
`FullFitKalman.java`	+118	-5	1.2 -> 1.3
`RKGeom.java`	+3	-2	1.3 -> 1.4
	+276	-20