2 added + 8 modified, total 10 files
hps-java/src/main/java/org/lcsim/HPSUsers/mgraham
diff -N PropXYXY_Test.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ PropXYXY_Test.java 24 May 2011 18:07:42 -0000 1.1
@@ -0,0 +1,524 @@
+/*
+ * PropXYXY_Test.java
+ *
+ * Created on July 24, 2007, 10:15 PM
+ *
+ * $Id: PropXYXY_Test.java,v 1.1 2011/05/24 18:07:42 mgraham Exp $
+ */
+
+package org.lcsim.HPSUsers.mgraham;
+
+//import junit.framework.TestCase;
+import org.lcsim.HPSKalmanFilter.PropXYXY;
+import org.lcsim.recon.tracking.spacegeom.SpacePath;
+import org.lcsim.recon.tracking.spacegeom.SpacePoint;
+import org.lcsim.recon.tracking.trfbase.ETrack;
+import org.lcsim.recon.tracking.trfbase.PropDir;
+import org.lcsim.recon.tracking.trfbase.PropStat;
+import org.lcsim.recon.tracking.trfbase.Propagator;
+import org.lcsim.recon.tracking.trfbase.Surface;
+import org.lcsim.recon.tracking.trfbase.TrackDerivative;
+import org.lcsim.recon.tracking.trfbase.TrackError;
+import org.lcsim.recon.tracking.trfbase.TrackVector;
+import org.lcsim.recon.tracking.trfbase.VTrack;
+import org.lcsim.recon.tracking.trfutil.Assert;
+import org.lcsim.recon.tracking.trfxyp.SurfXYPlane;
+import org.lcsim.util.Driver;
+
+/**
+ *
+ * @author Norman Graf
+ */
+public class PropXYXY_Test extends Driver
+{
+ private boolean debug=true;
+ /** Creates a new instance of PropXYXY_Test */
+ public void testPropXYXY()
+ {
+ String ok_prefix = "PropXYXY (I): ";
+ String error_prefix = "PropXYXY test (E): ";
+
+ if(debug) System.out.println( ok_prefix
+ + "-------- Testing component PropXYXY. --------" );
+
+
+ //********************************************************************
+
+ if(debug) System.out.println( ok_prefix + "Test constructor." );
+ PropXYXY prop = new PropXYXY(0.5);
+ if(debug) System.out.println( prop );
+
+ //********************************************************************
+
+ // Here we propagate some tracks both forward and backward and then
+ // each back to the original track. We check that the returned
+ // track parameters match those of the original track.
+ if(debug) System.out.println( ok_prefix + "Check reversibility." );
+
+ double u1[] ={10., 10., 10., 10., 10., 10.,10.,10.};
+ double u2[] ={10., 10., 10., 10., 10., 10.,11.,11.};
+ double phi1[]={ Math.PI/2., Math.PI/2., Math.PI/2., Math.PI/2., Math.PI/2., Math.PI/2.,0,0};
+ double phi2[]={5*Math.PI/6., 5*Math.PI/6., 5*Math.PI/6., 5*Math.PI/6., 5*Math.PI/6., 5*Math.PI/6.,0,0};
+ int sign_du[]={ 1, -1, -1, 1, 1, 1 , 1 ,1};
+ double v[] ={ -2., 2., 40., 40., -2., -2., 2. ,2 };
+ double z[] ={ 3., 3., 3., 3., 3., 3., 3. ,3};
+ double dvdu[]={-1.5, -1.5, 1.5, 1.5, -1.5, 0., 1.5 ,1.5 };
+ double dzdu[]={ 2.3, -2.3, -2.3, 2.3, 0., 2.3, 2.3 ,-2.3};
+ double qp[] ={ 0.05, -0.05, 0.05, -0.05, 0.05, 0.05, 0.05 ,0.5};
+
+ double maxdiff = 1.e-7;
+ int ntrk = 8;
+ int i;
+ for ( i=0; i<ntrk; ++i )
+ {
+ if(debug) System.out.println( "********** Propagate track " + i + ". **********" );
+ PropStat pstat = new PropStat();
+ SurfXYPlane sxyp1 = new SurfXYPlane(u1[i],phi1[i]);
+ SurfXYPlane sxyp2= new SurfXYPlane(u2[i],phi2[i]);
+ TrackVector vec1 = new TrackVector();
+ vec1.set(SurfXYPlane.IV ,v[i]); // v
+ vec1.set(SurfXYPlane.IZ ,z[i]); // z
+ vec1.set(SurfXYPlane.IDVDU ,dvdu[i]); // dv/du
+ vec1.set(SurfXYPlane.IDZDU ,dzdu[i]); // dz/du
+ vec1.set(SurfXYPlane.IQP ,qp[i]); // q/p
+ VTrack trv1 = new VTrack(sxyp1.newPureSurface(),vec1);
+ if (sign_du[i]==1) trv1.setForward();
+ else trv1.setBackward();
+ if(debug) System.out.println( " starting: " + trv1 );
+
+ VTrack trv2f = new VTrack(trv1);
+ pstat = prop.vecDirProp(trv2f,sxyp2,PropDir.FORWARD_MOVE);
+ if(debug) System.out.println( pstat );
+ if(debug) System.out.println("pstat= "+pstat);
+ Assert.assertTrue( pstat.forward() );
+ if(debug) System.out.println( " forward: " + trv2f );
+ Assert.assertTrue(check_dz(trv1,trv2f)>=0.);
+
+ VTrack trv2b = new VTrack(trv1);
+ pstat = prop.vecDirProp(trv2b,sxyp2,PropDir.BACKWARD);
+ if(debug) System.out.println( pstat );
+ Assert.assertTrue( pstat.backward() );
+ if(debug) System.out.println( " backward: " + trv2b );
+ Assert.assertTrue(check_dz(trv1,trv2b)<=0.);
+
+ VTrack trv2fb = new VTrack(trv2f);
+ pstat = prop.vecDirProp(trv2fb,sxyp1,PropDir.BACKWARD_MOVE);
+ if(debug) System.out.println( pstat );
+ Assert.assertTrue( pstat.backward() );
+ if(debug) System.out.println( " f return: " + trv2fb );
+ Assert.assertTrue(check_dz(trv2f,trv2fb)<=0.);
+
+ VTrack trv2bf = new VTrack(trv2b);
+ pstat = prop.vecDirProp(trv2bf,sxyp1,PropDir.FORWARD);
+ if(debug) System.out.println( pstat );
+ Assert.assertTrue( pstat.forward() );
+ if(debug) System.out.println( " b return: " + trv2bf );
+ Assert.assertTrue(check_dz(trv2b,trv2bf)>=0.);
+
+ double difff =
+ sxyp1.vecDiff(trv2fb.vector(),trv1.vector()).amax();
+ double diffb =
+ sxyp1.vecDiff(trv2bf.vector(),trv1.vector()).amax();
+ if(debug) System.out.println( "diffs: " + difff + ' ' + diffb );
+ Assert.assertTrue( difff < maxdiff );
+ Assert.assertTrue( diffb < maxdiff );
+
+ }
+
+ //********************************************************************
+
+ // Repeat the above with errors.
+ if(debug) System.out.println( ok_prefix + "Check reversibility with errors." );
+ double evv[] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
+ double evz[] = { 0.01, -0.01, 0.01, -0.01, 0.01, -0.01 };
+ double ezz[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, };
+ double evdv[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004 };
+ double ezdv[] = { 0.04, -0.04, 0.04, -0.04, 0.04, -0.04, };
+ double edvdv[] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
+ double evdz[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004 };
+ double edvdz[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004 };
+ double ezdz[] = { 0.04, -0.04, 0.04, -0.04, 0.04, -0.04 };
+ double edzdz[] = { 0.02, 0.02, 0.02, 0.02, 0.02, 0.02 };
+ double evqp[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004 };
+ double ezqp[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004 };
+ double edvqp[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004 };
+ double edzqp[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004 };
+ double eqpqp[] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
+
+ maxdiff = 1.e-6;
+
+ for ( i=0; i<ntrk; ++i )
+ {
+ if(debug) System.out.println( "********** Propagate track " + i + ". **********" );
+ PropStat pstat = new PropStat();
+ SurfXYPlane sxyp1 = new SurfXYPlane(u1[i],phi1[i]);
+ SurfXYPlane sxyp2 = new SurfXYPlane(u2[i],phi2[i]);
+ TrackVector vec1 = new TrackVector();
+ vec1.set(SurfXYPlane.IV, v[i]); // v
+ vec1.set(SurfXYPlane.IZ, z[i]); // z
+ vec1.set(SurfXYPlane.IDVDU, dvdu[i]); // dv/du
+ vec1.set(SurfXYPlane.IDZDU, dzdu[i]); // dz/du
+ vec1.set(SurfXYPlane.IQP, qp[i]); // q/p
+
+ TrackError err1 = new TrackError();
+ err1.set(SurfXYPlane.IV,SurfXYPlane.IV, evv[i]);
+ err1.set(SurfXYPlane.IV,SurfXYPlane.IZ, evz[i]);
+ err1.set(SurfXYPlane.IZ,SurfXYPlane.IZ, ezz[i]);
+ err1.set(SurfXYPlane.IV,SurfXYPlane.IDVDU, evdv[i]);
+ err1.set(SurfXYPlane.IZ,SurfXYPlane.IDVDU, ezdv[i]);
+ err1.set(SurfXYPlane.IDVDU,SurfXYPlane.IDVDU, edvdv[i]);
+ err1.set(SurfXYPlane.IV,SurfXYPlane.IDZDU, evdz[i]);
+ err1.set(SurfXYPlane.IZ,SurfXYPlane.IDZDU, ezdz[i]);
+ err1.set(SurfXYPlane.IDVDU,SurfXYPlane.IDZDU, edvdz[i]);
+ err1.set(SurfXYPlane.IDZDU,SurfXYPlane.IDZDU, edzdz[i]);
+ err1.set(SurfXYPlane.IV,SurfXYPlane.IQP, evqp[i]);
+ err1.set(SurfXYPlane.IZ,SurfXYPlane.IQP, ezqp[i]);
+ err1.set(SurfXYPlane.IDVDU,SurfXYPlane.IQP, edvqp[i]);
+ err1.set(SurfXYPlane.IDZDU,SurfXYPlane.IQP, edzqp[i]);
+ err1.set(SurfXYPlane.IQP,SurfXYPlane.IQP, eqpqp[i]);
+ ETrack trv1 = new ETrack(sxyp1.newPureSurface(),vec1,err1);
+ if(sign_du[i]==1) trv1.setForward();
+ else trv1.setBackward();
+ if(debug) System.out.println( " starting: " + trv1 );
+
+ ETrack trv2f = new ETrack(trv1);
+ pstat = prop.errDirProp(trv2f,sxyp2,PropDir.FORWARD);
+ Assert.assertTrue( pstat.forward() );
+ if(debug) System.out.println( " forward: " + trv2f );
+ ETrack trv2b = new ETrack(trv1);
+ pstat = prop.errDirProp(trv2b,sxyp2,PropDir.BACKWARD);
+ Assert.assertTrue( pstat.backward() );
+ if(debug) System.out.println( " backward: " + trv2b );
+ ETrack trv2fb = new ETrack(trv2f);
+ pstat = prop.errDirProp(trv2fb,sxyp1,PropDir.BACKWARD);
+ Assert.assertTrue( pstat.backward() );
+ if(debug) System.out.println( " f return: " + trv2fb );
+ ETrack trv2bf = new ETrack(trv2b);
+ pstat = prop.errDirProp(trv2bf,sxyp1,PropDir.FORWARD);
+ Assert.assertTrue( pstat.forward() );
+ if(debug) System.out.println( " b return: " + trv2bf );
+ double difff =
+ sxyp1.vecDiff(trv2fb.vector(),trv1.vector()).amax();
+ double diffb =
+ sxyp1.vecDiff(trv2bf.vector(),trv1.vector()).amax();
+ if(debug) System.out.println( "vec diffs: " + difff + ' ' + diffb );
+ Assert.assertTrue( difff < maxdiff );
+ Assert.assertTrue( diffb < maxdiff );
+ TrackError dfb = trv2fb.error().minus(trv1.error());
+ TrackError dbf = trv2bf.error().minus(trv1.error());
+ double edifff = dfb.amax();
+ double ediffb = dbf.amax();
+ if(debug) System.out.println( "err diffs: " + edifff + ' ' + ediffb );
+ Assert.assertTrue( edifff < maxdiff );
+ Assert.assertTrue( ediffb < maxdiff );
+ }
+
+ //********************************************************************
+
+ if(debug) System.out.println( ok_prefix + "Test Nearest Propagation" );
+
+ PropStat pstat = new PropStat();
+ SurfXYPlane sxyp1 = new SurfXYPlane(2.,Math.PI/3.);
+ SurfXYPlane sxyp2 = new SurfXYPlane(3.,Math.PI/3.);
+
+ TrackVector vec1 = new TrackVector();
+
+ vec1.set(SurfXYPlane.IV ,1.); // v
+ vec1.set(SurfXYPlane.IZ ,1.); // z
+ vec1.set(SurfXYPlane.IDVDU ,1.); // dv/du
+ vec1.set(SurfXYPlane.IDZDU ,1.); // dz/du
+ vec1.set(SurfXYPlane.IQP ,0.01); // q/p
+
+ VTrack trv1 = new VTrack(sxyp1.newPureSurface(),vec1);
+ trv1.setForward();
+
+ if(debug) System.out.println( " starting: " + trv1 );
+ VTrack trv2n = new VTrack(trv1);
+ pstat = prop.vecDirProp(trv2n,sxyp2,PropDir.NEAREST);
+ Assert.assertTrue( pstat.forward() );
+ if(debug) System.out.println( " nearest: " + trv2n );
+
+ trv1.setBackward();
+
+ if(debug) System.out.println( " starting: " + trv1 );
+ trv2n = new VTrack(trv1);
+ pstat = prop.vecDirProp(trv2n,sxyp2,PropDir.NEAREST);
+ Assert.assertTrue( pstat.backward() );
+ if(debug) System.out.println( " nearest: " + trv2n );
+
+ //********************************************************************
+
+ if(debug) System.out.println( ok_prefix + "Test XXX_MOVE and Same Surface Propagation." );
+
+ VTrack trvt0 = new VTrack(sxyp1.newPureSurface(),vec1);
+ trvt0.setForward();
+ VTrack trvt1 = new VTrack(trvt0);
+ PropStat tst = prop.vecDirProp(trvt1,sxyp1,PropDir.NEAREST);
+ Assert.assertTrue( tst.success() && trvt1.equals(trvt0) );
+ tst = prop.vecDirProp(trvt1,sxyp1,PropDir.FORWARD);
+ Assert.assertTrue( tst.success() && trvt1.equals(trvt0) );
+ tst = prop.vecDirProp(trvt1,sxyp1,PropDir.BACKWARD);
+ Assert.assertTrue( tst.success() && trvt1.equals(trvt0) );
+
+ trvt1 = new VTrack(trvt0);
+ tst = prop.vecDirProp(trvt1,sxyp1,PropDir.NEAREST_MOVE);
+ Assert.assertTrue( tst.success() && !trvt1.equals(trvt0) );
+ trvt1 = new VTrack(trvt0);
+ tst = prop.vecDirProp(trvt1,sxyp1,PropDir.FORWARD_MOVE);
+ Assert.assertTrue( tst.success() && !trvt1.equals(trvt0) );
+ trvt1 = new VTrack(trvt0);
+ tst = prop.vecDirProp(trvt1,sxyp1,PropDir.BACKWARD_MOVE);
+ Assert.assertTrue( tst.success() && !trvt1.equals(trvt0) );
+
+ //********************************************************************
+
+ if(debug) System.out.println( ok_prefix + "Test Zero B field propagation" );
+ {
+
+ PropXYXY prop0 = new PropXYXY(0.);
+ if(debug) System.out.println( prop0 );
+ Assert.assertTrue( prop0.bField() == 0. );
+
+ double u=10.,phi=0.;
+ Surface srf = new SurfXYPlane(u,phi);
+ VTrack trv0 = new VTrack(srf);
+ TrackVector vec = new TrackVector();
+ vec.set(SurfXYPlane.IV, 2.);
+ vec.set(SurfXYPlane.IZ, 10.);
+ vec.set(SurfXYPlane.IDVDU, 4.);
+ vec.set(SurfXYPlane.IDZDU, 2.);
+ trv0.setVector(vec);
+ trv0.setForward();
+ u=4.;
+
+ Surface srf_to = new SurfXYPlane(u,phi);
+
+ VTrack trv = new VTrack(trv0);
+ pstat = prop0.vecDirProp(trv,srf_to,PropDir.FORWARD);
+ Assert.assertTrue( !pstat.success() );
+
+ trv = new VTrack(trv0);
+ pstat = prop0.vecDirProp(trv,srf_to,PropDir.BACKWARD);
+ Assert.assertTrue( pstat.success() );
+
+ trv = new VTrack(trv0);
+ trv.setBackward();
+ pstat = prop0.vecDirProp(trv,srf_to,PropDir.BACKWARD);
+ Assert.assertTrue( !pstat.success() );
+
+ trv = new VTrack(trv0);
+ trv.setBackward();
+ pstat = prop0.vecDirProp(trv,srf_to,PropDir.FORWARD);
+ Assert.assertTrue( pstat.success() );
+
+ trv = new VTrack(trv0);
+ trv.setForward();
+ pstat = prop0.vecDirProp(trv,srf_to,PropDir.NEAREST);
+ Assert.assertTrue( pstat.success() );
+
+ Assert.assertTrue( pstat.backward() );
+ Assert.assertTrue( trv.vector(SurfXYPlane.IDVDU) == trv0.vector(SurfXYPlane.IDVDU) );
+ Assert.assertTrue( trv.vector(SurfXYPlane.IDZDU) == trv0.vector(SurfXYPlane.IDZDU) );
+ Assert.assertTrue(trv.surface().pureEqual(srf_to));
+
+ check_zero_propagation(trv0,trv,pstat);
+
+ srf_to = new SurfXYPlane(4.,Math.PI/16.);
+ trv = new VTrack(trv0);
+ trv.setForward();
+ pstat = prop0.vecDirProp(trv,srf_to,PropDir.NEAREST);
+ Assert.assertTrue( pstat.success() );
+ check_zero_propagation(trv0,trv,pstat);
+
+ srf_to = new SurfXYPlane(14.,Math.PI/4.);
+ trv = new VTrack(trv0);
+ trv.setForward();
+ pstat = prop0.vecDirProp(trv,srf_to,PropDir.NEAREST);
+ Assert.assertTrue( pstat.success() );
+ check_zero_propagation(trv0,trv,pstat);
+
+ srf_to = new SurfXYPlane(14.,Math.PI/2.);
+ trv = new VTrack(trv0);
+ trv.setForward();
+ pstat = prop0.vecDirProp(trv,srf_to,PropDir.NEAREST);
+ Assert.assertTrue( pstat.success() );
+ check_zero_propagation(trv0,trv,pstat);
+
+ srf_to = new SurfXYPlane(14.,Math.PI);
+ trv = new VTrack(trv0);
+ trv.setForward();
+ pstat = prop0.vecDirProp(trv,srf_to,PropDir.NEAREST);
+ Assert.assertTrue( pstat.success() );
+ check_zero_propagation(trv0,trv,pstat);
+
+ srf_to = new SurfXYPlane(14.,Math.PI*5./4.);
+ trv = new VTrack(trv0);
+ trv.setSurface(new SurfXYPlane(14.,Math.PI*7./4.));
+ trv.setForward();
+ VTrack tmp = new VTrack(trv);
+ VTrack der = new VTrack(trv);
+ pstat = prop0.vecDirProp(trv,srf_to,PropDir.NEAREST);
+ Assert.assertTrue( pstat.success() );
+ check_zero_propagation(tmp,trv,pstat);
+ check_derivatives(prop0,der,srf_to);
+
+ }
+
+
+ //********************************************************************
+
+ if(debug) System.out.println( ok_prefix + "Test cloning." );
+ Assert.assertTrue( prop.newPropagator() != null);
+
+ //********************************************************************
+
+ if(debug) System.out.println( ok_prefix + "Test the field." );
+ Assert.assertTrue( prop.bField() == 2.0 );
+
+ //********************************************************************
+
+ {
+ if(debug) System.out.println("===========================================\n");
+ PropXYXY p1 = new PropXYXY(2.);
+ PropXYXY p2 = new PropXYXY(-2.);
+ SurfXYPlane xy1 = new SurfXYPlane(5.,0.);
+ SurfXYPlane xy2 = new SurfXYPlane(3.,0.);
+ ETrack tre0 = new ETrack( xy2.newSurface() );
+ TrackVector vec = new TrackVector();
+ TrackError err = new TrackError();
+
+ vec.set(0, 0.);
+ vec.set(1, 0.);
+ vec.set(2, 0.0);
+ vec.set(3, 0.1);
+ vec.set(4, -0.1);
+ err.set(0,0, 0.1);
+ err.set(1,1, 0.1);
+ err.set(2,2, 0.1);
+ err.set(3,3, 0.1);
+ err.set(4,4, 0.1);
+ tre0.setVector(vec);
+ tre0.setError(err);
+ tre0.setBackward();
+
+ ETrack tre1 = new ETrack(tre0);
+ ETrack tre2 = new ETrack(tre0);
+ PropStat pstat2=p1.errDirProp(tre1,xy1,PropDir.BACKWARD);
+ Assert.assertTrue(pstat2.success());
+ pstat2=p2.errDirProp(tre2,xy1,PropDir.BACKWARD);
+ Assert.assertTrue(pstat2.success());
+
+ if(debug) System.out.println("tre1= \n"+tre1+'\n'+"tre2= \n"+tre2+'\n');
+ if(debug) System.out.println("===========================================\n");
+
+
+
+
+
+
+ //********************************************************************
+
+ if(debug) System.out.println( ok_prefix
+ + "------------- All tests passed. -------------" );
+ }
+ }
+
+ static void check_zero_propagation( VTrack trv0, VTrack trv, PropStat pstat)
+ {
+
+ SpacePoint sp = trv.spacePoint();
+ SpacePoint sp0 = trv0.spacePoint();
+
+ SpacePath sv = trv.spacePath();
+ SpacePath sv0 = trv0.spacePath();
+
+ Assert.assertTrue( Math.abs(sv0.dx() - sv.dx())<1e-7 );
+ Assert.assertTrue( Math.abs(sv0.dy() - sv.dy())<1e-7 );
+ Assert.assertTrue( Math.abs(sv0.dz() - sv.dz())<1e-7 );
+
+ double x0 = sp0.x();
+ double y0 = sp0.y();
+ double z0 = sp0.z();
+ double x1 = sp.x();
+ double y1 = sp.y();
+ double z1 = sp.z();
+
+ double dx = sv.dx();
+ double dy = sv.dy();
+ double dz = sv.dz();
+
+ double prod = dx*(x1-x0)+dy*(y1-y0)+dz*(z1-z0);
+ double moda = Math.sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0) + (z1-z0)*(z1-z0));
+ double modb = Math.sqrt(dx*dx+dy*dy+dz*dz);
+ double st = pstat.pathDistance();
+ Assert.assertTrue( Math.abs(prod-st) < 1.e-7 );
+ Assert.assertTrue( Math.abs(Math.abs(prod) - moda*modb) < 1.e-7 );
+ }
+
+ static void check_derivatives( Propagator prop, VTrack trv0, Surface srf)
+ {
+ for(int i=0;i<4;++i)
+ for(int j=0;j<4;++j)
+ check_derivative(prop,trv0,srf,i,j);
+ }
+
+ static void check_derivative( Propagator prop, VTrack trv0, Surface srf,int i,int j)
+ {
+
+ double dx = 1.e-3;
+ VTrack trv = new VTrack(trv0);
+ TrackVector vec = trv.vector();
+ boolean forward = trv0.isForward();
+
+ VTrack trv_0 = new VTrack(trv0);
+ TrackDerivative der = new TrackDerivative();
+ PropStat pstat = prop.vecProp(trv_0,srf,der);
+ Assert.assertTrue(pstat.success());
+
+ TrackVector tmp= new TrackVector(vec);
+ tmp.set(j, tmp.get(j)+dx);
+ trv.setVector(tmp);
+ if(forward) trv.setForward();
+ else trv.setBackward();
+
+ VTrack trv_pl = new VTrack(trv);
+ pstat = prop.vecProp(trv_pl,srf);
+ Assert.assertTrue(pstat.success());
+
+ TrackVector vecpl = trv_pl.vector();
+
+ tmp= new TrackVector(vec);
+ tmp.set(j,tmp.get(j)-dx);
+ trv.setVector(tmp);
+ if(forward) trv.setForward();
+ else trv.setBackward();
+
+ VTrack trv_mn = new VTrack(trv);
+ pstat = prop.vecProp(trv_mn,srf);
+ Assert.assertTrue(pstat.success());
+
+ TrackVector vecmn = trv_mn.vector();
+
+ double dy = (vecpl.get(i)-vecmn.get(i))/2.;
+
+ double dydx = dy/dx;
+
+ double didj = der.get(i,j);
+
+ if( Math.abs(didj) > 1e-10 )
+ Assert.assertTrue( Math.abs((dydx - didj)/didj) < 1e-4 );
+ else
+ Assert.assertTrue( Math.abs(dydx) < 1e-4 );
+ }
+ //**********************************************************************
+ //Checker of correct z propagation
+
+ static double check_dz(VTrack trv1, VTrack trv2)
+ {
+ double z1 = trv1.vector().get(SurfXYPlane.IZ);
+ double z2 = trv2.vector().get(SurfXYPlane.IZ);
+ double dzdu = trv1.vector().get(SurfXYPlane.IDZDU);
+ int sign_du = (trv1.isForward())?1:-1;
+ return (z2-z1)*dzdu*sign_du;
+ }
+}
hps-java/src/main/java/org/lcsim/HPSUsers/mgraham
diff -u -r1.1 -r1.2
--- KalmanGeom.java 16 May 2011 17:45:09 -0000 1.1
+++ KalmanGeom.java 24 May 2011 18:07:42 -0000 1.2
@@ -11,6 +11,7 @@
import java.util.Set;
import org.lcsim.HPSKalmanFilter.KalmanSurface;
import org.lcsim.HPSKalmanFilter.PropDCAXY;
+import org.lcsim.HPSKalmanFilter.PropXYXY;
//import org.lcsim.HPSKalmanFilter.ShapeDispatcher;
import org.lcsim.HPSKalmanFilter.ShapeDispatcher;
@@ -60,7 +61,7 @@
import org.lcsim.recon.tracking.trfxyp.ClusXYPlane1;
import org.lcsim.recon.tracking.trfxyp.ClusXYPlane2;
import org.lcsim.recon.tracking.trfxyp.HitXYPlane2;
-import org.lcsim.recon.tracking.trfxyp.PropXYXY;
+
import org.lcsim.recon.tracking.trfxyp.SurfXYPlane;
/**
@@ -70,9 +71,9 @@
*
*
*@author $Author: mgraham $
- *@version $Id: KalmanGeom.java,v 1.1 2011/05/16 17:45:09 mgraham Exp $
+ *@version $Id: KalmanGeom.java,v 1.2 2011/05/24 18:07:42 mgraham Exp $
*
- * Date $Date: 2011/05/16 17:45:09 $
+ * Date $Date: 2011/05/24 18:07:42 $
*
*/
/* To make the Kalman filter work for any detector, the type of hit that is
@@ -84,7 +85,8 @@
* looping through detector elements, etc. that may be useful for finding
* intersections at some point. (ecfine) */
public class KalmanGeom {
- boolean _DEBUG=false;
+// boolean _DEBUG=true;
+ boolean _DEBUG=true;
// TRF wants distances in cm, not mm.
private double mmTocm = 0.1;
// this flag is a temporary fix to avoid making a surface behind the origin in x.
@@ -403,8 +405,8 @@
double dist1 = dotProduct(strip1.origin(), strip1.w());
double dist2 = dotProduct(strip2.origin(), strip2.w());
-
- System.out.println("TrackerHit Position = [" + thit.getPosition()[0] + "," + thit.getPosition()[1] + "," + thit.getPosition()[2]);
+ if(track.slope()<0)
+ System.out.println("TrackerHit Position = [" + thit.getPosition()[0] + "," + thit.getPosition()[1] + "," + thit.getPosition()[2]);
@@ -510,8 +512,8 @@
double dzv1 = tmpcov1.e(1, 2);
double dzv2 = tmpcov2.e(1, 2);
if(_DEBUG){
- System.out.println("z1 = " + z1 + "+/-" + Math.sqrt(dz1) + ", v1 = " + v1 + "+/-" + Math.sqrt(dv1) + ", cov(z,v) = " + dzv1);
- System.out.println("z2 = " + z2 + "+/-" + Math.sqrt(dz2) + ", v2 = " + v2 + "+/-" + Math.sqrt(dv2) + ", cov(z,v) = " + dzv2);
+ if(track.slope()<0) System.out.println("z1 = " + z1 + "+/-" + Math.sqrt(dz1) + ", v1 = " + v1 + "+/-" + Math.sqrt(dv1) + ", cov(z,v) = " + dzv1);
+ if(track.slope()<0) System.out.println("z2 = " + z2 + "+/-" + Math.sqrt(dz2) + ", v2 = " + v2 + "+/-" + Math.sqrt(dv2) + ", cov(z,v) = " + dzv2);
}
// ClusXYPlane2 cluster1 = new ClusXYPlane2(dist1, phi1, v1, z1, dv1, dz1,dzv1);
// ClusXYPlane2 cluster2 = new ClusXYPlane2(dist2, phi2, v2, z2, dv2, dz2,dzv2);
@@ -527,9 +529,9 @@
hit1.setParentPointer(cluster1);
hit2.setParentPointer(cluster2);
if(_DEBUG){
- System.out.println("hit1 predicted vector: " + hit1.toString());
- System.out.println("hit2 predicted vector: " + hit2.toString());
- System.out.println("hit position from trackerhit: [" + thit.getPosition()[0] + ", " +
+ if(track.slope()<0) System.out.println("hit1 predicted vector: " + hit1.toString());
+ if(track.slope()<0) System.out.println("hit2 predicted vector: " + hit2.toString());
+ if(track.slope()<0) System.out.println("hit position from trackerhit: [" + thit.getPosition()[0] + ", " +
thit.getPosition()[1] + ", " + thit.getPosition()[2] + "]");
}
hps-java/src/main/java/org/lcsim/HPSUsers/mgraham
diff -u -r1.1 -r1.2
--- KalmanFilterDriver.java 16 May 2011 17:45:09 -0000 1.1
+++ KalmanFilterDriver.java 24 May 2011 18:07:42 -0000 1.2
@@ -64,7 +64,9 @@
System.out.println("geom field = " + geom.bz + ", trackUtils field = " + trackUtils.bz);
fitk = new FullFitKalman(prop);
-//
+ PropXYXY_Test foobar=new PropXYXY_Test();
+ foobar.testPropXYXY();
+ //
}
public void process(EventHeader event) {
@@ -143,7 +145,7 @@
if (kfchi2 < 100) {
aida.histogram1D(kaldirPos + "Kalman Chi2", 50, 0, 100).fill(kfchi2);
- aida.histogram1D(kaldirPos + "Kalman Fit DCA", 50, -15, 15).fill(kfDCA);
+ aida.histogram1D(kaldirPos + "Kalman Fit DCA", 50, -0.5, 0.5).fill(kfDCA);
aida.histogram1D(kaldirPos + "Kalman Fit z0", 50, -3, 3).fill(kfz0);
aida.histogram1D(kaldirPos + "Kalman Fit sin(phi0)", 50, -0.5, 0.5).fill(kfphi0);
aida.histogram1D(kaldirPos + "Kalman Fit tanlambda", 50, -0.1, 0.1).fill(kfslope);
@@ -153,7 +155,7 @@
aida.histogram1D(kaldirPos + "Kalman-Nominal sin(phi0)", 50, -0.1, 0.1).fill(kfphi0 - nomphi0);
aida.histogram1D(kaldirPos + "Kalman-Nominal slope", 50, -0.2, 0.2).fill(kfslope - nomslope);
aida.histogram1D(kaldirPos + "Kalman-Nominal qOverpt", 50, -0.2, 0.2).fill(kfqOverpt - nomqOverpt);
- aida.histogram1D(kaldirPos + "Nominal Fit DCA", 50, -15, 15).fill(nomDCA);
+ aida.histogram1D(kaldirPos + "Nominal Fit DCA", 50, -0.5, 0.5).fill(nomDCA);
aida.histogram1D(kaldirPos + "Nominal Fit z0", 50, -3, 3).fill(nomz0);
aida.histogram1D(kaldirPos + "Nominal Fit sin(phi0)", 50, -0.5, 0.5).fill(nomphi0);
aida.histogram1D(kaldirPos + "Nominal Fit tanlambda", 50, -0.1, 0.1).fill(nomslope);
hps-java/src/main/java/org/lcsim/HPSKalmanFilter
diff -N PropXYXY.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ PropXYXY.java 24 May 2011 18:07:42 -0000 1.1
@@ -0,0 +1,839 @@
+package org.lcsim.HPSKalmanFilter;
+
+import org.lcsim.recon.tracking.trfutil.TRFMath;
+import org.lcsim.recon.tracking.trfutil.Assert;
+
+import org.lcsim.recon.tracking.trfbase.Propagator;
+import org.lcsim.recon.tracking.trfbase.PropDirected;
+import org.lcsim.recon.tracking.trfbase.PropDir;
+import org.lcsim.recon.tracking.trfbase.PropStat;
+import org.lcsim.recon.tracking.trfbase.VTrack;
+import org.lcsim.recon.tracking.trfbase.Surface;
+import org.lcsim.recon.tracking.trfbase.TrackVector;
+import org.lcsim.recon.tracking.trfbase.TrackDerivative;
+import org.lcsim.recon.tracking.trfxyp.SurfXYPlane;
+
+/**
+ * Propagates tracks from one XYPlane to another in a constant field.
+ *<p>
+ * Propagation will fail if either the origin or destination is
+ * not a XYPlane.
+ * Propagator works incorrectly for tracks with very small curvatures
+ *
+ *
+ *@author Norman A. Graf
+ *@version 1.0
+ *
+ */
+public class PropXYXY extends PropDirected
+{
+
+ // attributes
+
+ private double _bfac;
+
+ private static final int IV = SurfXYPlane.IV;
+ private static final int IZ = SurfXYPlane.IZ;
+ private static final int IDVDU = SurfXYPlane.IDVDU;
+ private static final int IDZDU = SurfXYPlane.IDZDU;
+ private static final int IQP = SurfXYPlane.IQP;
+
+ // static methods
+
+ /**
+ *Return a String representation of the class' type name.
+ *Included for completeness with the C++ version.
+ *
+ * @return A String representation of the class' type name.
+ */
+ public static String typeName()
+ { return "PropXYXY"; }
+
+ /**
+ *Return a String representation of the class' type name.
+ *Included for completeness with the C++ version.
+ *
+ * @return A String representation of the class' type name.
+ */
+ public static String staticType()
+ { return typeName(); }
+
+
+ /**
+ *Construct an instance from a constant solenoidal magnetic field in Tesla.
+ *
+ * @param bfield The magnetic field strength in Tesla.
+ */
+ public PropXYXY(double bfield)
+ {
+ _bfac = TRFMath.BFAC*bfield;
+ }
+
+ /**
+ *Clone an instance.
+ *
+ * @return A Clone of this instance.
+ */
+ public Propagator newPropagator( )
+ {
+ return new PropXYXY( bField() );
+ }
+
+ /**
+ *Propagate a track without error in the specified direction
+ *and return the derivative matrix in deriv.
+ *
+ * @param trv The VTrack to propagate.
+ * @param srf The Surface to which to propagate.
+ * @param dir The direction in which to propagate.
+ * @param deriv The track derivatives to update at the surface srf.
+ * @return The propagation status.
+ **/
+ public PropStat vecDirProp( VTrack trv, Surface srf,
+ PropDir dir,TrackDerivative deriv )
+ {
+ PropStat pstat = vec_propagatexyxy_( _bfac, trv, srf, dir, deriv );
+ return pstat;
+ }
+
+ /**
+ *Propagate a track without error in the specified direction.
+ *
+ * @param trv The VTrack to propagate.
+ * @param srf The Surface to which to propagate.
+ * @param dir The direction in which to propagate.
+ * @return The propagation status.
+ */
+ @Override
+ public PropStat vecDirProp( VTrack trv, Surface srf,
+ PropDir dir)
+ {
+ TrackDerivative deriv = null;
+ return vecDirProp(trv, srf, dir, deriv);
+ }
+
+ /**
+ *Return the strength of the magnetic field in Tesla.
+ *
+ * @return The strength of the magnetic field in Tesla.
+ */
+ public double bField()
+ {
+ return _bfac/TRFMath.BFAC;
+ }
+
+
+ /**
+ *Return a String representation of the class' type name.
+ *Included for completeness with the C++ version.
+ *
+ * @return A String representation of the class' type name.
+ */
+ public String type()
+ { return staticType(); }
+
+
+ /**
+ *output stream
+ *
+ * @return A String representation of this instance.
+ */
+ @Override
+ public String toString()
+ {
+ return "XYPlane-XYPlane propagation with constant "
+ + bField() + " Tesla field";
+ }
+
+
+ // Private function to determine dphi of the propagation
+ double direction(int flag_forward, int sign_dphi,
+ double du,
+ double norm, double cat,
+ double sinphi, double cosphi)
+ {
+
+ int sign_cat = 0;
+ if(du*flag_forward>0.) sign_cat = 1;
+ if(du*flag_forward<0.) sign_cat = -1;
+ if(du == 0.)
+ {
+ if(sinphi >=0. ) sign_cat = 1;
+ if(sinphi < 0. ) sign_cat = -1;
+ }
+
+ double sin_dphi = norm*sinphi + cat*cosphi*sign_cat;
+ double cos_dphi = -norm*cosphi + cat*sinphi*sign_cat;
+
+ int sign_sindphi = 0;
+ if(sin_dphi> 0.) sign_sindphi = 1;
+ if(sin_dphi< 0.) sign_sindphi = -1;
+ if(sin_dphi == 0.) sign_sindphi = sign_dphi;
+
+ if( Math.abs(cos_dphi)>1. )
+ {
+ double sn= -1.;
+ if ( cos_dphi > 0 ) sn=1.;
+ cos_dphi= sn*Math.sqrt(1.-sin_dphi*sin_dphi);
+ }
+ if( Math.abs(sin_dphi)>1. )
+ {
+ double sn= -1.;
+ if ( sin_dphi > 0 ) sn=1.;
+ sin_dphi= sn*Math.sqrt(1.-cos_dphi*cos_dphi);
+ }
+
+ double dphi = Math.PI*(sign_dphi - sign_sindphi) + sign_sindphi*Math.acos(cos_dphi);
+ return dphi;
+
+ }
+ //**********************************************************************
+
+ // Private function to propagate a track without error
+ // The corresponding track parameters are:
+ // u (cm) is fixed
+ // 0 - v (cm)
+ // 1 - z (cm)
+ // 2 - dv/du
+ // 3 - dz/du
+ // 4 - q/p p is momentum of a track, q is its charge
+ // If deriv is nonzero, return the derivative matrix there.
+
+ PropStat
+ vec_propagatexyxy_( double B, VTrack trv, Surface srf,
+ PropDir dir1,
+ TrackDerivative deriv )
+ {
+
+ // construct return status
+ PropStat pstat = new PropStat();
+
+ PropDir dir = dir1; //need to check constness of this
+ boolean move = Propagator.reduceDirection(dir);
+ if(move) dir = reduce(dir);
+ // fetch the originating surface and vector
+ Surface srf1 = trv.surface();
+ // TrackVector vec1 = trv.get_vector();
+
+ // Check origin is a XYPlane.
+ Assert.assertTrue( srf1.pureType().equals(SurfXYPlane.staticType()) );
+ if ( !srf1.pureType( ).equals(SurfXYPlane.staticType()) )
+ return pstat;
+ SurfXYPlane sxyp1 = (SurfXYPlane) srf1;
+
+ // Check destination is an XYPlane.
+ Assert.assertTrue( srf.pureType().equals(SurfXYPlane.staticType()) );
+ if ( !srf.pureType( ).equals(SurfXYPlane.staticType()) )
+ return pstat;
+ SurfXYPlane sxyp2 = (SurfXYPlane) srf;
+
+ // If surfaces are the same and no XXX_MOVE is requested we can return now.
+ boolean same = srf.pureEqual(srf1);
+ if ( same && ! move )
+ {
+ if(deriv!=null)
+ {
+ deriv.setIdentity();
+ }
+ pstat.setSame();
+ return pstat;
+ }
+
+
+ if ( same && dir1.equals(PropDir.NEAREST_MOVE) ) dir = PropDir.FORWARD;
+
+ //v00.63.04 cng 01/12/01
+ if( Math.abs(B) < 1e-7 ) return zeroBField(trv,sxyp1,sxyp2,dir1,deriv);
+
+ // Fetch the u's and phi's of the planes and the starting track vector.
+ TrackVector vec = trv.vector();
+ int iphi = SurfXYPlane.NORMPHI;
+ int idist = SurfXYPlane.DISTNORM;
+ double phi_0 = sxyp1.parameter(iphi);
+ double u_0 = sxyp1.parameter(idist);
+ double phi_n = sxyp2.parameter(iphi);
+ double u_n = sxyp2.parameter(idist);
+ if ( same && move )
+ {
+ int forw= trv.isForward()? 1 : -1 ;
+ int forw_dir = (dir == PropDir.FORWARD) ? 1 : -1 ;
+ int sn = vec.get(IDVDU) > 0. ? 1 : -1 ;
+ u_0 += 1.e-7 * (double)(forw*forw_dir*sn);
+ }
+
+ double b1_0 = vec.get(IV); // v
+ double b2_0 = vec.get(IZ); // z
+ double b3_0 = vec.get(IDVDU); // dv/du
+ double b4_0 = vec.get(IDZDU); // dz/du
+ double b5_0 = vec.get(IQP); // q/p
+
+ double phi_u = phi_0 - phi_n;
+
+ double cosphi_u = Math.cos(phi_u);
+ double sinphi_u = Math.sin(phi_u);
+// System.out.println( "PropXYXY._vec_propagate: phi_u = "+phi_u);
+ // check if du == 0 ( that is track moves parallel to the destination plane )
+ double du_du_0 = cosphi_u - b3_0*sinphi_u;
+ if(du_du_0==0.) return pstat;
+
+ double a_hat_u = 1./du_du_0;
+ double a_hat_u2 = a_hat_u*a_hat_u;
+
+ double u = u_0*cosphi_u - b1_0*sinphi_u;
+ double b1 = b1_0*cosphi_u + u_0*sinphi_u;
+ double b2 = b2_0;
+ double b3 = (b3_0*cosphi_u + sinphi_u)*a_hat_u;
+ double b4 = b4_0*a_hat_u;
+ double b5 = b5_0;
+
+ int sign_du0 = 0;
+ if(trv.isForward()) sign_du0 = 1;
+ if(trv.isBackward()) sign_du0 = -1;
+ if(sign_du0 == 0)
+ {
+ System.out.println( "PropXYXY._vec_propagate: Unknown direction of a track ");
+ System.exit(1);
+ }
+ int sign_du = 0;
+ if(du_du_0*sign_du0 > 0) sign_du = 1;
+ if(du_du_0*sign_du0 < 0) sign_du = -1;
+
+ // check that q/p != 0
+ Assert.assertTrue(b5 !=0. );
+
+ // 1/curv of the track is r
+ double r = 1/(b5*B)*Math.sqrt(1 + b3_0*b3_0)/Math.sqrt(1 + b3_0*b3_0 + b4_0*b4_0);
+ double b3_hat = Math.sqrt(1 + b3*b3);
+ double b34_hat = Math.sqrt(1 + b3*b3 + b4*b4);
+ double b3_hat2 = b3_hat*b3_hat;
+ double b34_hat2 = b34_hat*b34_hat;
+ double cosphi = -b3*sign_du/b3_hat;
+ double sinphi = sign_du/b3_hat;
+ double rsinphi = 1./(b5*B)*sign_du/b34_hat;
+ double rcosphi = -b3/(b5*B)*sign_du/b34_hat;
+
+ double du = u_n - u;
+ double norm = du/r - cosphi;
+ System.out.println( "PropXYXY._vec_propagate: u_n = "+u_n+"; u = "+u+"; r = "+r);
+
+ // xy-xy propagation failed : noway to the new plane
+ if(Math.abs(norm)>1.) return pstat;
+
+ double cat = Math.sqrt(1.-norm*norm);
+ int flag_forward = 0;
+ int sign_dphi = 0;
+
+ if( dir.equals(PropDir.NEAREST))
+ {
+ // try forward propagation
+ flag_forward = 1;
+ if(b5*B>0.) sign_dphi = 1;
+ if(b5*B<0.) sign_dphi = -1;
+ double dphi1=
+ direction(flag_forward,sign_dphi,du,norm, cat, sinphi, cosphi);
+ // try backward propagation
+ flag_forward = -1;
+ if(b5*B>0.) sign_dphi = -1;
+ if(b5*B<0.) sign_dphi = 1;
+ double dphi2=
+ direction(flag_forward,sign_dphi,du,norm, cat, sinphi, cosphi);
+ if( Math.abs(dphi2)>Math.abs(dphi1) )
+ {
+ flag_forward = -flag_forward;
+ sign_dphi = - sign_dphi;
+ }
+ }
+
+ else if( dir.equals(PropDir.FORWARD))
+ {
+ flag_forward = 1;
+ if(b5*B>0.) sign_dphi = 1;
+ if(b5*B<0.) sign_dphi = -1;
+ }
+ else if( dir.equals(PropDir.BACKWARD))
+ {
+ flag_forward = -1;
+ if(b5*B>0.) sign_dphi = -1;
+ if(b5*B<0.) sign_dphi = 1;
+ }
+ else
+ {
+ System.out.println("PropXYXY._vec_propagate: Unknown direction.");
+ System.exit(1);
+ }
+
+ int sign_cat = 0;
+ if(du*sign_dphi*b5*B>0.) sign_cat = 1;
+ if(du*sign_dphi*b5*B<0.) sign_cat = -1;
+ if(du == 0.)
+ {
+ if(sinphi >=0. ) sign_cat = 1;
+ if(sinphi < 0. ) sign_cat = -1;
+ }
+
+ double sin_dphi = norm*sinphi + cat*cosphi*sign_cat;
+ double cos_dphi = -norm*cosphi + cat*sinphi*sign_cat;
+ if(cos_dphi>1.) cos_dphi=1.;
+ if(cos_dphi<-1.) cos_dphi= -1.;
+
+ int sign_sindphi = 0;
+ if(sin_dphi> 0.) sign_sindphi = 1;
+ if(sin_dphi< 0.) sign_sindphi = -1;
+ if(sin_dphi == 0.) sign_sindphi = sign_dphi;
+
+ double dphi = Math.PI*(sign_dphi - sign_sindphi) + sign_sindphi*Math.acos(cos_dphi);
+
+ // check that I didn't make any mistakes
+
+ Assert.assertTrue(Math.abs(Math.sin(dphi) - sin_dphi)<1.e-5);
+
+ // check if dun == 0 ( that is track moves parallel to the destination plane)
+ double du_n_du = cos_dphi - b3*sin_dphi;
+ if(du_n_du==0.) return pstat;
+
+ double a_hat_dphi = 1./du_n_du;
+ double a_hat_dphi2 = a_hat_dphi*a_hat_dphi;
+ double c_hat_dphi = sin_dphi + b3*cos_dphi ;
+
+ double b1_n = b1 + rsinphi*(1-cos_dphi) - rcosphi*sin_dphi;
+ double b2_n = b2 + b4/(b5*B)*sign_du/b34_hat*dphi;
+ double b3_n = c_hat_dphi*a_hat_dphi;
+ double b4_n = b4*a_hat_dphi;
+ double b5_n = b5;
+ System.out.println( "PropXYXY._vec_propagate: sign_cat="+sign_cat+"; cosphi="+cosphi
+ +"; sinphi="+sinphi+"; cat="+cat);
+ System.out.println( "PropXYXY._vec_propagate: a_hat_dphi="+a_hat_dphi+"; cos_dphi="+cos_dphi
+ +"; norm="+norm);
+ System.out.println( "PropXYXY._vec_propagate: b4_0="+b4_0+"; b4="+b4+"; b4_n="+b4_n);
+
+ // double u_n_0 = u_n*cosphi_u + b1_n*sinphi_u;
+
+ // check if track crossed original plane during the propagation
+ // switch (dir) {
+ // if( dir.equals(PropDir.FORWARD:
+ // if((u_n_0 - u_0)*sign_du0<0) return pstat;
+ // break;
+ // if( dir.equals(PropDir.BACKWARD:
+ // if((u_n_0 - u_0)*sign_du0>0) return pstat;
+ // break;
+ // }
+
+ int sign_dun = 0;
+ if(du_n_du*sign_du > 0) sign_dun = 1;
+ if(du_n_du*sign_du < 0) sign_dun = -1;
+
+ vec.set(IV ,b1_n);
+ vec.set(IZ ,b2_n);
+ vec.set(IDVDU ,b3_n);
+ vec.set(IDZDU ,b4_n);
+ vec.set(IQP ,b5_n);
+
+ // Update trv
+ trv.setSurface(srf.newPureSurface());
+ trv.setVector(vec);
+
+ // set new direction of the track
+ if(sign_dun == 1) trv.setForward();
+ if(sign_dun == -1) trv.setBackward();
+
+ // Calculate sT.
+ double crv = B*b5;
+ double dv = b1_n - b1;
+ double dxy = Math.sqrt( du*du + dv*dv );
+ double arg = 0.5*crv*dxy;
+ double dst = dxy*TRFMath.asinrat(arg);
+
+ // Calculate s.
+ double dz = b2_n - b2;
+ Assert.assertTrue( flag_forward==1 || flag_forward==-1 );
+ double ds = ((double)(flag_forward))*Math.sqrt(dst*dst+dz*dz);
+
+ // Set the return status.
+ pstat.setPathDistance(ds);
+ //(flag_forward==1)?pstat.set_forward():pstat.set_backward();
+
+ // exit now if user did not ask for error matrix.
+ if ( deriv == null ) return pstat;
+
+ // du_d0
+
+ double du_db1_0 = - sinphi_u;
+
+ // db1_d0
+
+ double db1_db1_0 = cosphi_u;
+
+ // db3_d0
+
+ double db3_db3_0 = a_hat_u2;
+
+ // db4_d0
+
+ double db4_db3_0 = b4_0*sinphi_u*a_hat_u2;
+ double db4_db4_0 = a_hat_u;
+
+ // dr_d
+
+ double dr_db3 = r*b3*b4*b4/(b3_hat2*b34_hat2);
+ double dr_db4 = -r*b4/b34_hat2;
+ double dr_db5 = -r/b5;
+
+ // dcosphi_d
+
+ double dcosphi_db3 = - sign_du/b3_hat - cosphi*b3/b3_hat2;
+
+ // dsinphi_d
+
+ double dsinphi_db3 = - sinphi*b3/b3_hat2;
+
+ // dcat_d
+
+ double dcat_db3 = norm/cat*(du/(r*r)*dr_db3 + dcosphi_db3 );
+ double dcat_db4 = norm/cat* du/(r*r)*dr_db4;
+ double dcat_db5 = norm/cat* du/(r*r)*dr_db5;
+ double dcat_du = norm/(cat*r);
+
+ // dnorm_d
+
+ double dnorm_db3 = - du/(r*r)*dr_db3 - dcosphi_db3;
+ double dnorm_db4 = - du/(r*r)*dr_db4;
+ double dnorm_db5 = - du/(r*r)*dr_db5;
+ double dnorm_du = - 1./r;
+
+ // dcos_dphi_d
+
+ double dcos_dphi_db3 = - cosphi*dnorm_db3 - norm*dcosphi_db3 +
+ sign_cat*(sinphi*dcat_db3 + cat*dsinphi_db3);
+ double dcos_dphi_db4 = - cosphi*dnorm_db4 + sign_cat*sinphi*dcat_db4;
+ double dcos_dphi_db5 = - cosphi*dnorm_db5 + sign_cat*sinphi*dcat_db5;
+ double dcos_dphi_du = - cosphi*dnorm_du + sign_cat*sinphi*dcat_du;
+
+ // dsin_dphi_d
+
+ double dsin_dphi_db3 = sinphi*dnorm_db3 + norm*dsinphi_db3 +
+ sign_cat*(cosphi*dcat_db3 + cat*dcosphi_db3);
+ double dsin_dphi_db4 = sinphi*dnorm_db4 + sign_cat*cosphi*dcat_db4;
+ double dsin_dphi_db5 = sinphi*dnorm_db5 + sign_cat*cosphi*dcat_db5;
+ double dsin_dphi_du = sinphi*dnorm_du + sign_cat*cosphi*dcat_du;
+
+ // ddphi_d
+
+ double ddphi_db3;
+ double ddphi_db4;
+ double ddphi_db5;
+ double ddphi_du;
+ if(Math.abs(sin_dphi)>0.5)
+ {
+ ddphi_db3 = - dcos_dphi_db3/sin_dphi;
+ ddphi_db4 = - dcos_dphi_db4/sin_dphi;
+ ddphi_db5 = - dcos_dphi_db5/sin_dphi;
+ ddphi_du = - dcos_dphi_du /sin_dphi;
+ }
+ else
+ {
+ ddphi_db3 = dsin_dphi_db3/cos_dphi;
+ ddphi_db4 = dsin_dphi_db4/cos_dphi;
+ ddphi_db5 = dsin_dphi_db5/cos_dphi;
+ ddphi_du = dsin_dphi_du /cos_dphi;
+ }
+
+ // da_hat_dphi_d
+
+ double da_hat_dphi_db3 = - a_hat_dphi2*
+ (dcos_dphi_db3 - sin_dphi - b3*dsin_dphi_db3);
+ double da_hat_dphi_db4 = - a_hat_dphi2*(dcos_dphi_db4 - b3*dsin_dphi_db4);
+ double da_hat_dphi_db5 = - a_hat_dphi2*(dcos_dphi_db5 - b3*dsin_dphi_db5);
+ double da_hat_dphi_du = - a_hat_dphi2*(dcos_dphi_du - b3*dsin_dphi_du );
+
+ // dc_hat_dphi_d
+
+ double dc_hat_dphi_db3 = b3*dcos_dphi_db3 + dsin_dphi_db3 + cos_dphi;
+ double dc_hat_dphi_db4 = b3*dcos_dphi_db4 + dsin_dphi_db4;
+ double dc_hat_dphi_db5 = b3*dcos_dphi_db5 + dsin_dphi_db5;
+ double dc_hat_dphi_du = b3*dcos_dphi_du + dsin_dphi_du ;
+
+ // db1_n_d
+
+ double db1_n_db1 = 1;
+ double db1_n_db3 = (dr_db3*sinphi+r*dsinphi_db3)*(1-cos_dphi)
+ - rsinphi*dcos_dphi_db3
+ - dr_db3*cosphi*sin_dphi - r*dcosphi_db3*sin_dphi
+ - rcosphi*dsin_dphi_db3;
+ double db1_n_db4 = dr_db4*sinphi*(1-cos_dphi) - rsinphi*dcos_dphi_db4
+ - dr_db4*cosphi*sin_dphi - rcosphi*dsin_dphi_db4;
+ double db1_n_db5 = dr_db5*sinphi*(1-cos_dphi) - rsinphi*dcos_dphi_db5
+ - dr_db5*cosphi*sin_dphi - rcosphi*dsin_dphi_db5;
+ double db1_n_du = - rsinphi*dcos_dphi_du - rcosphi*dsin_dphi_du;
+
+ // db2_n_d
+
+ double db2_n_db2 = 1.;
+ double db2_n_db3 = 1./(b5*B)*b4*sign_du/b34_hat*
+ ( - dphi*b3/b34_hat2 + ddphi_db3);
+ double db2_n_db4 = 1./(b5*B)*sign_du/b34_hat*
+ ( - dphi*b4*b4/b34_hat2 + b4*ddphi_db4 + dphi );
+ double db2_n_db5 = 1./(b5*B)*b4*sign_du/b34_hat*( ddphi_db5 - dphi/b5);
+ double db2_n_du = 1./(b5*B)*b4*sign_du/b34_hat*ddphi_du;
+
+ // db3_n_d
+
+ double db3_n_db3 = a_hat_dphi*dc_hat_dphi_db3 + da_hat_dphi_db3*c_hat_dphi;
+ double db3_n_db4 = a_hat_dphi*dc_hat_dphi_db4 + da_hat_dphi_db4*c_hat_dphi;
+ double db3_n_db5 = a_hat_dphi*dc_hat_dphi_db5 + da_hat_dphi_db5*c_hat_dphi;
+ double db3_n_du = a_hat_dphi*dc_hat_dphi_du + da_hat_dphi_du *c_hat_dphi;
+
+ // db4_n_d
+
+ double db4_n_db3 = b4*da_hat_dphi_db3;
+ double db4_n_db4 = b4*da_hat_dphi_db4 + a_hat_dphi;
+ double db4_n_db5 = b4*da_hat_dphi_db5;
+ double db4_n_du = b4*da_hat_dphi_du;
+
+ // db5_n_d
+
+ // db1_n_d0
+
+ double db1_n_db1_0 = db1_n_du * du_db1_0 + db1_n_db1*db1_db1_0;
+ double db1_n_db2_0 = 0.;
+ double db1_n_db3_0 = db1_n_db3*db3_db3_0 + db1_n_db4*db4_db3_0;
+ double db1_n_db4_0 = db1_n_db4*db4_db4_0;
+ double db1_n_db5_0 = db1_n_db5;
+
+ // db2_n_d0
+
+ double db2_n_db1_0 = db2_n_du *du_db1_0;
+ double db2_n_db2_0 = db2_n_db2;
+ double db2_n_db3_0 = db2_n_db3*db3_db3_0 + db2_n_db4*db4_db3_0;
+ double db2_n_db4_0 = db2_n_db4*db4_db4_0;
+ double db2_n_db5_0 = db2_n_db5;
+
+ // db3_n_d0
+
+ double db3_n_db1_0 = db3_n_du*du_db1_0;
+ double db3_n_db2_0 = 0.;
+ double db3_n_db3_0 = db3_n_db3*db3_db3_0 + db3_n_db4*db4_db3_0;
+ double db3_n_db4_0 = db3_n_db4*db4_db4_0;
+ double db3_n_db5_0 = db3_n_db5;
+
+ // db4_n_d0
+
+ double db4_n_db1_0 = db4_n_du*du_db1_0;
+ double db4_n_db2_0 = 0.;
+ double db4_n_db3_0 = db4_n_db3*db3_db3_0 + db4_n_db4*db4_db3_0;
+ double db4_n_db4_0 = db4_n_db4*db4_db4_0;
+ double db4_n_db5_0 = db4_n_db5;
+
+ // db5_n_d0
+
+ double db5_n_db1_0 = 0.;
+ double db5_n_db2_0 = 0.;
+ double db5_n_db3_0 = 0.;
+ double db5_n_db4_0 = 0.;
+ double db5_n_db5_0 = 1.;
+
+
+ deriv.set(IV,IV ,db1_n_db1_0);
+ deriv.set(IV,IZ ,db1_n_db2_0);
+ deriv.set(IV,IDVDU ,db1_n_db3_0);
+ deriv.set(IV,IDZDU ,db1_n_db4_0);
+ deriv.set(IV,IQP ,db1_n_db5_0);
+ deriv.set(IZ,IV ,db2_n_db1_0);
+ deriv.set(IZ,IZ ,db2_n_db2_0);
+ deriv.set(IZ,IDVDU ,db2_n_db3_0);
+ deriv.set(IZ,IDZDU ,db2_n_db4_0);
+ deriv.set(IZ,IQP ,db2_n_db5_0);
+ deriv.set(IDVDU,IV ,db3_n_db1_0);
+ deriv.set(IDVDU,IZ ,db3_n_db2_0);
+ deriv.set(IDVDU,IDVDU ,db3_n_db3_0);
+ deriv.set(IDVDU,IDZDU ,db3_n_db4_0);
+ deriv.set(IDVDU,IQP ,db3_n_db5_0);
+ deriv.set(IDZDU,IV ,db4_n_db1_0);
+ deriv.set(IDZDU,IZ ,db4_n_db2_0);
+ deriv.set(IDZDU,IDVDU ,db4_n_db3_0);
+ deriv.set(IDZDU,IDZDU ,db4_n_db4_0);
+ deriv.set(IDZDU,IQP ,db4_n_db5_0);
+ deriv.set(IQP,IV ,db5_n_db1_0);
+ deriv.set(IQP,IZ ,db5_n_db2_0);
+ deriv.set(IQP,IDVDU ,db5_n_db3_0);
+ deriv.set(IQP,IDZDU ,db5_n_db4_0);
+ deriv.set(IQP,IQP ,db5_n_db5_0);
+
+ return pstat;
+
+ }
+
+ //cng
+ PropStat zeroBField( VTrack trv, SurfXYPlane sxyp1,
+ SurfXYPlane sxyp2, PropDir dir1,
+ TrackDerivative deriv)
+ {
+ // construct return status
+ PropStat pstat = new PropStat();
+
+ PropDir dir = dir1; //need to check constness of this
+ boolean move = Propagator.reduceDirection(dir);
+ boolean same = sxyp2.pureEqual(sxyp1);
+
+ // There is only one solution. Can't XXX_MOVE
+ if ( same && move ) return pstat;
+
+ if ( same )
+ {
+ if(deriv!=null)
+ {
+ deriv.setIdentity();
+ }
+ pstat.setSame();
+ return pstat;
+ }
+
+ TrackVector vec = trv.vector();
+ double v0 = vec.get(IV);
+ double z0 = vec.get(IZ);
+ double dvdu0 = vec.get(IDVDU);
+ double dzdu0 = vec.get(IDZDU);
+
+ double du0 =1.;
+ if( trv.isBackward() ) du0 = -1.;
+
+ double phi0 = sxyp1.parameter(SurfXYPlane.NORMPHI);
+ double cphi0 = Math.cos(phi0);
+ double sphi0 = Math.sin(phi0);
+ double u0 = sxyp1.parameter(SurfXYPlane.DISTNORM);
+
+ double phi1 = sxyp2.parameter(SurfXYPlane.NORMPHI);
+ double cphi1 = Math.cos(phi1);
+ double sphi1 = Math.sin(phi1);
+ double u1 = sxyp2.parameter(SurfXYPlane.DISTNORM);
+
+ double a = (cphi0 - dvdu0*sphi0)*du0;
+ double b = (sphi0 + dvdu0*cphi0)*du0;
+ double c = dzdu0*du0;
+
+ double x0 = u0*cphi0 - v0*sphi0;
+ double y0 = u0*sphi0 + v0*cphi0;
+
+ double ap = u1*cphi1;
+ double bp = u1*sphi1;
+ double cp = 0;
+
+ double xp = ap;
+ double yp = bp;
+ double zp = 0.0;
+
+ double denom = a*ap + b*bp + c*cp;
+
+ if( denom == 0. ) return pstat;
+
+ double S = ( (xp-x0)*ap + (yp-y0)*bp + (zp-z0)*cp )/denom;
+
+ double x1 = x0 + S*a;
+ double y1 = y0 + S*b;
+ double z1 = z0 + S*c;
+
+ boolean forward = S > 0. ? true : false;
+
+ if( dir == PropDir.FORWARD && !forward ) return pstat;
+ if( dir == PropDir.BACKWARD && forward ) return pstat;
+
+ double v1 = y1*cphi1 - x1*sphi1;
+
+ double v01 = y0*cphi1 - x0*sphi1;
+ double u01 = y0*sphi1 + x0*cphi1;
+ double z01 = z0;
+
+ if( u01 == u1 ) return pstat;
+
+ double dvdu1 = (v1-v01)/(u1-u01);
+ double dzdu1 = (z1-z01)/(u1-u01);
+
+ vec.set(IV, v1);
+ vec.set(IZ, z1);
+ vec.set(IDVDU, dvdu1);
+ vec.set(IDZDU, dzdu1);
+
+ // Update trv
+ trv.setSurface(sxyp2.newPureSurface());
+ trv.setVector(vec);
+ // set new direction of the track
+ if( b*sphi1 + a*cphi1 >0 ) trv.setForward();
+ else trv.setBackward();
+
+ // Calculate s.
+ double ds = S*Math.sqrt(a*a+b*b+c*c);
+
+ // Set the return status.
+ pstat.setPathDistance(ds);
+
+ if( deriv == null ) return pstat;
+
+ double dx0dv0 = -sphi0;
+ double dy0dv0 = cphi0;
+
+ double dadv_du = -sphi0*du0;
+ double dbdv_du = cphi0*du0;
+ double dcdz_du = du0;
+
+ double ddenomdv_du = dadv_du*ap + dbdv_du*bp;
+ double ddenomdz_du = dcdz_du*cp;
+
+
+ double dSdv0 = -(dx0dv0*ap + dy0dv0*bp)/denom;
+ double dSdz0 = -cp/denom;
+ double dSdv_du = -S/denom*ddenomdv_du;
+ double dSdz_du = -S/denom*ddenomdz_du;
+
+ double dy1dv0 = dy0dv0 + dSdv0*b;
+ double dx1dv0 = dx0dv0 + dSdv0*a;
+ double dx1dv_du = dSdv_du*a + S*dadv_du;
+ double dy1dv_du = dSdv_du*b + S*dbdv_du;
+
+ double du01dv0 = dy0dv0*sphi1 + dx0dv0*cphi1;
+ double dv01dv0 = dy0dv0*cphi1 - dx0dv0*sphi1;
+
+ double dv1dv0 = dy1dv0*cphi1 - dx1dv0*sphi1;
+ double dv1dv_du = dy1dv_du*cphi1 - dx1dv_du*sphi1;
+
+ double dz1dz0 = 1 + dSdz0*c;
+ double dz1dv0 = dSdv0*c;
+ double dz1dv_du = dSdv_du*c;
+ double dz1dz_du = dSdz_du*c + S*dcdz_du;
+
+ double dv_du1dv0 = ((dv1dv0-dv01dv0)*(u1-u01) + du01dv0*(v1-v01))/(u1-u01)/(u1-u01);
+ double dv_du1dv_du = dv1dv_du/(u1-u01);
+
+ double dz_du1dv0 = (dz1dv0*(u1-u01) + du01dv0*(z1-z01))/(u1-u01)/(u1-u01);
+ double dz_du1dz0 = (dz1dz0 - 1.)/(u1-u01);
+ double dz_du1dv_du = dz1dv_du/(u1-u01);
+ double dz_du1dz_du = dz1dz_du/(u1-u01);
+
+ //Set the track derivatives...
+ deriv.setIdentity();
+
+ deriv.set(IV,IV, dv1dv0);
+ deriv.set(IV,IDVDU, dv1dv_du);
+
+ deriv.set(IZ,IV, dz1dv0);
+ deriv.set(IZ,IZ, dz1dz0);
+ deriv.set(IZ,IDVDU, dz1dv_du);
+ deriv.set(IZ,IDZDU, dz1dz_du);
+
+ deriv.set(IDVDU,IV, dv_du1dv0);
+ deriv.set(IDVDU,IDVDU, dv_du1dv_du);
+
+ deriv.set(IDZDU,IV, dz_du1dv0);
+ deriv.set(IDZDU,IZ, dz_du1dz0);
+ deriv.set(IDZDU,IDVDU, dz_du1dv_du);
+ deriv.set(IDZDU,IDZDU, dz_du1dz_du);
+
+ deriv.set(IQP,IQP, 1.0);
+
+ return pstat;
+
+ }
+ //cng
+
+
+
+
+}
hps-java/src/main/java/org/lcsim/HPSKalmanFilter
diff -u -r1.2 -r1.3
--- FullFitKalman.java 16 May 2011 17:45:09 -0000 1.2
+++ FullFitKalman.java 24 May 2011 18:07:42 -0000 1.3
@@ -47,14 +47,14 @@
* a single track.
*
*@author $Author: mgraham $
- *@version $Id: FullFitKalman.java,v 1.2 2011/05/16 17:45:09 mgraham Exp $
+ *@version $Id: FullFitKalman.java,v 1.3 2011/05/24 18:07:42 mgraham Exp $
*
- * Date $Date: 2011/05/16 17:45:09 $
+ * Date $Date: 2011/05/24 18:07:42 $
*
*/
public class FullFitKalman extends FullFitter {
- boolean _DEBUG = false;
+ boolean _DEBUG = true;
private AIDA aida = AIDA.defaultInstance();
// Flags to control: multiple scattering, energy loss and adding the hit.
private boolean doMs = true;
@@ -190,7 +190,7 @@
System.out.println("pstat = " + pstat);
System.out.println("After prop: " + trh.newTrack());
}
- if(!pstat.success())
+ if (!pstat.success())
return -666;
// System.out.println("trh= \n"+trh+", hit= \n"+hit);
// System.out.println("_addfit= "+_addfit);
@@ -213,11 +213,14 @@
//
}
// Multiple scattering--this is a silly way to do it..
- ThinXYPlaneMs interactor = new ThinXYPlaneMs(.01);
+// ThinXYPlaneMs interactor = new ThinXYPlaneMs(.01);
+
+ System.out.println("trh error pre MS: " + trh.newTrack().error().toString());
+ ThinXYPlaneMs interactor = new ThinXYPlaneMs(.005);
ETrack tre2 = trh.newTrack();
interactor.interact(tre2);
trh.setFit(tre2, trh.chisquared());
-
+ System.out.println("trh error pre MS: " + trh.newTrack().error().toString());
if (fstat > 0)
return 10000 + 1000 * fstat + icount;
hps-java/src/main/java/org/lcsim/HPSKalmanFilter
diff -u -r1.2 -r1.3
--- TrackUtils.java 16 May 2011 17:45:09 -0000 1.2
+++ TrackUtils.java 24 May 2011 18:07:42 -0000 1.3
@@ -6,6 +6,8 @@
import Jama.Matrix;
import hep.physics.matrix.SymmetricMatrix;
+import org.lcsim.contrib.Pelham.analysis.HelixParamCalculator;
+import org.lcsim.event.MCParticle;
import org.lcsim.fit.helicaltrack.HelicalTrackFit;
import org.lcsim.recon.tracking.trfbase.TrackError;
import org.lcsim.recon.tracking.trfbase.TrackVector;
@@ -33,7 +35,9 @@
TrackVector tv = new TrackVector();
tv.set(0, r_signed * mmTocm);
tv.set(1, t.z0() * mmTocm);
- tv.set(2, t.phi0());
+// tv.set(2, t.phi0());
+ //I think this should by tanphi
+ tv.set(2, Math.tan(t.phi0()));
// tv.set(3, t.cth() );
tv.set(3, t.slope());
// tv.set(4, t.curvature()/bz);
@@ -48,13 +52,13 @@
// SurfDCA s = new SurfDCA((t.dca() * Math.sin(t.phi0())),
// (-t.dca() * Math.cos(t.phi0())));
- SurfDCA s = new SurfDCA(0., 0.);
+ SurfDCA s = new SurfDCA(r_signed * mmTocm, 0.);
VTrack vt = new VTrack(s, tv);
if (_DEBUG) {
System.out.println("making VTrack with: ");
System.out.println(" r_signed = " + r_signed * mmTocm);
System.out.println(" z0 = " + t.z0() * mmTocm);
- System.out.println(" phi0 = " + t.phi0());
+ System.out.println(" tanphi0 = " + Math.tan(t.phi0()));
System.out.println(" tanlamda = " + t.slope());
System.out.println(" q/pt = " + -t.curvature() / mmTocm / (TRFMath.BFAC * bz));
// System.out.println(" q/pt = " + t.curvature() * 0.299999 * bz);
@@ -69,6 +73,21 @@
}
return vt;
}
+ //make a VTrack from an MC particle
+ public VTrack makeVTrack(MCParticle mcp) {
+ TrackVector tv = new TrackVector();
+ HelixParamCalculator helix = new HelixParamCalculator(mcp, bz);
+ double r_signed = -(helix.getDCA());
+ double curvemc = (-1) / helix.getRadius();
+ tv.set(0, r_signed * mmTocm);
+ tv.set(1, helix.getZ0() * mmTocm);
+ tv.set(2, helix.getPhi0());
+ tv.set(3, helix.getSlopeSZPlane());
+ tv.set(4, curvemc / mmTocm / (TRFMath.BFAC * bz)); // curvature to q/pt
+ SurfDCA s = new SurfDCA(r_signed * mmTocm, 0.);
+ VTrack vt = new VTrack(s, tv);
+ return vt;
+ }
public void setBZ(double b) {
bz = b;
hps-java/src/main/java/org/lcsim/HPSUsers/Example
diff -u -r1.6 -r1.7
--- HeavyPhotonLLDriver.java 21 Apr 2011 17:07:38 -0000 1.6
+++ HeavyPhotonLLDriver.java 24 May 2011 18:07:42 -0000 1.7
@@ -19,8 +19,8 @@
//********************
// set pname to the detector you want to run...HPS1pt8 is the default test detector (uses strategy file 1.3)
// see "setPName" method for list
-// public String pName = "HPS1pt9";
- public String pName = "HPS3pt4";
+ public String pName = "HPS1pt8";
+// public String pName = "HPS3pt4";
JasAnalysisDriver jad;
DetailedAnalysisDriver dad;
FastTrackAnalysisDriver ftad;
@@ -55,9 +55,9 @@
MultiTrackReco trd = new MultiTrackReco(strategyPrefix, detType, bfield, axialStrategy, finalStrategy, llStrategy, pairs, passLayers);
add(trd);
- jad = new JasAnalysisDriver(nlayers);
+ jad = new JasAnalysisDriver(nlayers,fitlayers,detType);
ftad = new FastTrackAnalysisDriver(nlayers,fitlayers,detType);
- dad = new DetailedAnalysisDriver(nlayers);
+ dad = new DetailedAnalysisDriver(nlayers,fitlayers,detType);
// add(ftad); //uncomment to run the fasttrackanalysisdriver (text file output)
// add(jad);
add(dad);
hps-java/src/main/java/org/lcsim/HPSUsers/Example
diff -u -r1.2 -r1.3
--- DetailedAnalysisDriver.java 19 May 2011 17:21:48 -0000 1.2
+++ DetailedAnalysisDriver.java 24 May 2011 18:07:42 -0000 1.3
@@ -121,14 +121,22 @@
FileWriter fw;
PrintWriter pw;
double[] beamsize = {0.001, 0.02, 0.02};
-
- public DetailedAnalysisDriver(int layers) {
- nlayers[0] = layers;
-
+ Integer _minLayers=8;
+ // flipSign is a kludge...
+// HelicalTrackFitter doesn't deal with B-fields in -ive Z correctly
+// so we set the B-field in +iveZ and flip signs of fitted tracks
+// note: this should be -1 for Test configurations and +1 for Full (v3.X and lower) configurations
+// this is set by the _config variable (detType in HeavyPhotonDriver)
+ int flipSign=1;
+ String _config="3pt4";
+ public DetailedAnalysisDriver(int trackerLayers,int mintrkLayers, String config) {
+ nlayers[0] = trackerLayers;
+ _minLayers=mintrkLayers;
+ _config=config;
+ if(_config.contains("Test")) flipSign=-1;
// Define the efficiency histograms
IHistogramFactory hf = aida.histogramFactory();
-
peffFindable = hf.createProfile1D("Findable Efficiency vs p", "", 20, 0., beamP);
thetaeffFindable = hf.createProfile1D("Findable Efficiency vs theta", "", 20, 80, 100);
phieffFindable = hf.createProfile1D("Findable Efficiency vs phi", "", 25, -0.25, 0.25);
@@ -414,9 +422,9 @@
double slopeErr = Math.sqrt(track.getErrorMatrix().e(HelicalTrackFit.slopeIndex, HelicalTrackFit.slopeIndex));
double curveErr = Math.sqrt(track.getErrorMatrix().e(HelicalTrackFit.curvatureIndex, HelicalTrackFit.curvatureIndex));
_nchRec++;
- if (track.getCharge() < 0)
+ if (track.getCharge()*flipSign < 0)
_neleRec++;
- if (track.getCharge() > 0)
+ if (track.getCharge()*flipSign > 0)
_nposRec++;
SeedTrack stEle = (SeedTrack) track;
@@ -480,9 +488,9 @@
// Make plots for fake, non-fake, and all tracks
if (purity < 0.5) {
- if (track.getCharge() < 0)
+ if (track.getCharge()*flipSign < 0)
_neleFake++;
- if (track.getCharge() > 0)
+ if (track.getCharge()*flipSign > 0)
_nposFake++;
cthfake.fill(cth, 1.0);
phifake.fill(phi, 1.0);
@@ -491,9 +499,9 @@
fillTrackInfo(trackdir, "fake tracks", track.getChi2(), nhits, p, pperp, px, py, pz, phi, cth, d0, xoca, yoca, z0);
} else {
- if (track.getCharge() < 0)
+ if (track.getCharge()*flipSign < 0)
_neleTru++;
- if (track.getCharge() > 0)
+ if (track.getCharge()*flipSign > 0)
_nposTru++;
cthfake.fill(cth, 0.0);
phifake.fill(phi, 0.0);
@@ -593,7 +601,7 @@
if (htlayer == 1)
l1DeltaZ.put(track, z - zTr);
- if (purity == 1 && track.getCharge() > 0 && nhits == 10) {
+ if (purity == 1 && track.getCharge()*flipSign > 0 && nhits == 10) {
if (clusterlist.get(0).rawhits().size() == 1 && clusterlist.get(1).rawhits().size() == 1) {
aida.cloud1D(hitdir + tkresid + "SingleStrip--Track delta y: Layer " + htlayer).fill(y - yTr);
aida.cloud1D(hitdir + tkresid + "SingleStrip--Track delta z: Layer " + htlayer).fill(z - zTr);
@@ -855,7 +863,7 @@
for (Track track1 : tracklist) {
Track ele = null;
Track pos = null;
- int ch1 = track1.getCharge();
+ int ch1 = track1.getCharge()*flipSign;
int index = tracklist.indexOf(track1);
List<Track> subtracklist = tracklist.subList(index, tracklist.size());
for (Track track2 : subtracklist) {
@@ -969,7 +977,7 @@
// Find the number of layers hit by this mc particle
// System.out.println("MC pt=" + pt);
int nhits = findable.LayersHit(mcp);
- boolean isFindable = findable.InnerTrackerIsFindable(mcp, nlayers[0] - 2);
+ boolean isFindable = findable.InnerTrackerIsFindable(mcp, _minLayers);
Set<SimTrackerHit> mchitlist = mcHittomcP.allTo(mcp);
Set<HelicalTrackCross> hitlist = hittomc.allTo(mcp);
@@ -998,7 +1006,7 @@
if (trktomc.allTo(d).isEmpty())
bothreco = false;
// if (findable.LayersHit(d) != nlayers[0])
- if (!findable.InnerTrackerIsFindable(d, nlayers[0] - 2))
+ if (!findable.InnerTrackerIsFindable(d, _minLayers))
bothfindable = false;
}
double vtxWgt = 0;
hps-java/src/main/java/org/lcsim/HPSUsers/Example
diff -u -r1.4 -r1.5
--- JasAnalysisDriver.java 20 Apr 2011 20:07:39 -0000 1.4
+++ JasAnalysisDriver.java 24 May 2011 18:07:42 -0000 1.5
@@ -118,8 +118,20 @@
PrintWriter pw;
boolean isBeamConstrain = false;
double[] beamsize = {0.001, 0.02, 0.02};
-
- public JasAnalysisDriver(int layers) {
+ Integer _minLayers=8;
+ // flipSign is a kludge...
+// HelicalTrackFitter doesn't deal with B-fields in -ive Z correctly
+// so we set the B-field in +iveZ and flip signs of fitted tracks
+// note: this should be -1 for Test configurations and +1 for Full (v3.X and lower) configurations
+// this is set by the _config variable (detType in HeavyPhotonDriver)
+ int flipSign=1;
+ String _config="3pt4";
+ public JasAnalysisDriver(int trackerLayers,int mintrkLayers, String config) {
+ nlayers[0] = trackerLayers;
+ _minLayers=mintrkLayers;
+ _config=config;
+ if(_config.contains("Test")) flipSign=-1;
+
// Define the efficiency histograms
IHistogramFactory hf = aida.histogramFactory();
@@ -151,8 +163,7 @@
VxEffFindable = hf.createProfile1D("Aprime Efficiency vs Vx: Findable", "", 25, 0., 50.);
VyEffFindable = hf.createProfile1D("Aprime Efficiency vs Vy: Findable", "", 40, -0.2, 0.2);
- VzEffFindable = hf.createProfile1D("Aprime Efficiency vs Vz: Findable", "", 40, -0.2, 0.2);
- nlayers[0] = layers;
+ VzEffFindable = hf.createProfile1D("Aprime Efficiency vs Vz: Findable", "", 40, -0.2, 0.2);
}
public void process(
@@ -416,7 +427,7 @@
List<BilliorTrack> btlist = new ArrayList<BilliorTrack>();
for (Track track1 : tracklist) {
for (Track track2 : tracklist) {
- if (track1 != track2 && track1.getCharge() > 0 && track2.getCharge() < 0) {
+ if (track1 != track2 && track1.getCharge()*flipSign > 0 && track2.getCharge()*flipSign < 0) {
/*
vlist.clear();
vlist.add(track1);
hps-java/src/main/java/org/lcsim/HPSBase
diff -u -r1.3 -r1.4
--- MultiTrackReco.java 19 Apr 2011 17:14:26 -0000 1.3
+++ MultiTrackReco.java 24 May 2011 18:07:42 -0000 1.4
@@ -67,6 +67,7 @@
hthdriver.setStereoPair("Tracker", pair[0], pair[1]);
if (detType.contains("Test")) {
+ System.out.println("Setting separation for a Test detector");
hthdriver.setMaxSeperation(10.01);
hthdriver.setTolerance(0.01);
} else {
CVSspam 0.2.8
