lcsim/sandbox/RobKutschke/TRFTests/v1
diff -u -r1.2 -r1.3
--- FitTest_sid00Driver.java 7 Sep 2007 19:04:15 -0000 1.2
+++ FitTest_sid00Driver.java 13 Sep 2007 22:36:02 -0000 1.3
@@ -29,9 +29,9 @@
* no multiple scattering or energy loss.
*
*@author $Author: kutschke $
- *@version $Id: FitTest_sid00Driver.java,v 1.2 2007/09/07 19:04:15 kutschke Exp $
+ *@version $Id: FitTest_sid00Driver.java,v 1.3 2007/09/13 22:36:02 kutschke Exp $
*
- * Date $Date: 2007/09/07 19:04:15 $
+ * Date $Date: 2007/09/13 22:36:02 $
*
*/
@@ -69,7 +69,7 @@
// Limits for track generation parameters.
param = new RKTrackGenParams();
param.setbz(rkgeom.getBz());
- param.setptmin( 10.0);
+ param.setptmin( 0.5);
//param.setczmin( 0.93);
//param.setczmax( 0.993);
param.setczmin( -0.993);
@@ -141,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();
@@ -153,7 +153,10 @@
RKDebug.Instance().setRKTrack(rktrk);
// Generate hits on this track in my own format.
- RKHitList hlist = GenerateOneArcHits( rktrk );
+ RKHitList hlist = GenerateMixedOneArcHits( rktrk );
+
+ // Diagnostics for the generated hits.
+ CheckHitList( hlist, rktrk );
// Check for sufficient hits.
if ( ! hlist.fitable() ){
@@ -181,7 +184,6 @@
}
// Do the outward fit.
- //System.out.println ("Start forward fit.");
int fstatf = fitk.fitForward(ht);
if ( fstatf != 0 ){
System.out.println("Forwards status: " +
@@ -213,7 +215,7 @@
aida.histogram1D("Smear Check",100, -5., 5.).fill(pull);
}
-
+ // Prep for the inward fit.
if ( vtend.surface().pureType().equals(SurfCylinder.staticType()) ){
vstart = vstartc_in;
RKDebug.Instance().setStartType(1);
@@ -222,13 +224,13 @@
RKDebug.Instance().setStartType(2);
}
- // Prep for inward fit.
+ // Seed inward fit with the output of the outward fit.
ETrack etin = new ETrack( vtend.surface().newPureSurface(),
ht.newTrack().vector(),
vstart );
- //vtend.vector(), vstart );
- // Not sure why this is needed.
+ // When the starting surface is a zplane, we need to tell the ETrack
+ // which direction it is forward.
if ( vtend.surface().type() == SurfZPlane.staticType() ){
if ( rktrk.cz() > 0 ){
etin.setForward();
@@ -244,7 +246,6 @@
}
// Do inward fit.
- // System.out.println ("Start backward fit.");
int fstatb = fitk.fitBackward(htin);
if ( fstatb != 0 ){
System.out.println("Backwards status: " +
@@ -320,8 +321,8 @@
rkl.addHit( new RKHit( rktrk, s, vtdca, ps, sum ) );
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",
r*Math.cos(vtdca.vector(0)),
r*Math.sin(vtdca.vector(0)),
@@ -336,18 +337,8 @@
}
}
}
- double xsum = 0.;
-
- // NaN checker?
- if ( sum > 0 ){
- xsum = sum;
- } else if ( sum < 0 ){
- xsum = sum;
- } else{
- System.out.println ("Can't get here ...");
- }
- aida.cloud2D("Path length vs cz",-1).fill(rktrk.cz(),xsum);
- aida.cloud1D("Sum").fill(xsum);
+ aida.cloud2D("Path length vs cz",-1).fill(rktrk.cz(),sum);
+ aida.cloud1D("Sum").fill(sum);
aida.cloud2D("Hits vs cz",-1).fill( rktrk.cz(), ncmeas );
sum = 0;
@@ -377,9 +368,267 @@
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;
+ }
+
+ // Generate perfect hits on the outward going arc only.
+ // Do not apply energy loss or multiple scattering.
+ private RKHitList GenerateMixedOneArcHits( RKTrack rktrk ){
+
+ // Convert to TRF style track.
+ toTRF trftrk = new toTRF(rktrk);
+ VTrack vtdca = trftrk.atDcaZaxis();
+
+ // This code only looks at the first arc of a track.
+ double pathlimit = trftrk.halfarc();
+
+ int nmeas = 0;
+ int nzmeas = 0;
+ int ncmeas = 0;
+ double sum = 0.;
+ RKHitList rkl = new RKHitList();
+ Propagator prop = rkgeom.newPropagator();
+
+ List<RKSurf> cyls = rkgeom.getCylinders();
+ List<RKSurf> zps = rkgeom.getZ( rktrk.z0(), rktrk.cz() );
+
+ int next_cyl = 0;
+ int next_zp = 0;
+ int max_cyl = cyls.size();
+ int max_zp = zps.size();
+
+ // At top of the loop, vtdca contains the track parameters at the
+ // departure point for this step.
+ while ( (next_cyl<max_cyl) || ( next_zp < max_zp ) ){
+
+ VTrack vt_cyl = null;
+ VTrack vt_zp = null;
+
+ double path_cyl = 0.;
+ double path_zp = 0.;
+
+ int status_cyl = 0;
+ int status_zp = 0;
+
+ RKSurf rks_cyl = null;
+ RKSurf rks_zp = null;
+
+ PropStat ps_cyl = null;
+ PropStat ps_zp = null;
+
+ double r_cyl = 0.;
+ double z_cyl = 0.;
+ double r_zp = 0.;
+ double z_zp = 0.;
+
+ int next_cyl_sav = next_cyl;
+ int next_zp_sav = next_zp;
+
+ // Check cylindrical surfaces until the next good one is found.
+ while ( next_cyl < max_cyl ){
+
+ // Default status is failure.
+ status_cyl = 0;
+
+ // Parameters at the starting surface for this step.
+ vt_cyl = new VTrack(vtdca);
+
+ // Move to the next surface.
+ 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() ) {
+ ++next_cyl;
+ continue;
+ }
+
+ // If we moved backwards, try the next cylinder outward.
+ path_cyl = ps_cyl.pathDistance();
+ if ( path_cyl < 0. ) {
+ //System.out.println ("Negative path at cyl..." );
+ status_cyl = 3;
+ ++next_cyl;
+ continue;
+ }
+ status_cyl = 1;
+
+ // Needed for downstream diagnostics.
+ r_cyl = rks_cyl.radius*mmTocm;
+ z_cyl = vt_cyl.vector(1);
+
+ if ( rks_cyl.inBounds(z_cyl) ){
+
+ // We have the next good surface from this list.
+ status_cyl = 2;
+ break;
+
+ } else{
+
+ // If not in bounds, try the next surface.
+ ++next_cyl;
+ continue;
+ }
+
+ }
+
+ // Check the next z surface.
+ while ( next_zp < max_zp ){
+
+ // Default status is failure.
+ status_zp = 0;
+
+ // Parameters at the start of this step.
+ vt_zp = new VTrack(vtdca);
+
+ // Move to the next surface.
+ rks_zp = zps.get(next_zp);
+ ps_zp = prop.vecDirProp( vt_zp, rks_zp.getZPlane(), PropDir.FORWARD );
+
+
+ // This happens when the z surface being tested is behind the starting point.
+ // Try the next surface.
+ if ( !ps_zp.success() ) {
+ ++next_zp;
+ continue;
+ }
+
+ // Normally this will not happen. If it does, try the next surface.
+ path_zp = ps_zp.pathDistance();
+ if ( path_zp < 0. ) {
+ status_zp = 3;
+ ++next_zp;
+ continue;
+ }
+ status_zp = 1;
+
+ // Needed for downstream diagnostics
+ r_zp = Math.sqrt(vt_zp.vector(0)*vt_zp.vector(0) +
+ vt_zp.vector(1)*vt_zp.vector(1) );
+ z_zp = rks_zp.zc*mmTocm;
+
+ if ( rks_zp.inBounds(r_zp) ){
+
+ // We have the next good surface from this list.
+ status_zp = 2;
+ break;
+ }else{
+
+ // Out of bounds so try the next surface.
+ ++next_zp;
+ continue;
+ }
+
+ }
+
+ // At this point we have zero, one or two good solutions.
+ // If two, choose the one with the shortest step length.
+ boolean addcyl = false;
+ boolean addzp = false;
+ if ( status_cyl == 2 && status_zp == 2 ){
+ if ( path_cyl < path_zp ) {
+ addcyl = true;
+ } else {
+ addzp = true;
+ }
+ } else if ( status_cyl == 2 ){
+ addcyl = true;
+
+ } else if ( status_zp == 2 ){
+ addzp = true;
+ } else {
+
+ // This should only happen when both lists are exhausted.
+ if ( next_cyl<max_cyl || next_zp < max_zp ){
+ System.out.println ("Don't know how I got here ... "
+ + next_cyl + " "
+ + max_cyl + " "
+ + next_zp + " "
+ + max_zp
+ );
+ }
+ }
+
+
+ // Add the chosen hit.
+ if ( addcyl ){
+
+ // Advance counter to next surface.
+ ++next_cyl;
+
+ // Restore the z index to it's position at the start of this step.
+ next_zp = next_zp_sav;
+
+ // Did we exceed the half-arc limit?
+ sum += path_cyl;
+ if ( sum > pathlimit ) break;
+
+ // Track parameters at this step.
+ vtdca = vt_cyl;
+
+ // Add the hit.
+ rkl.addHit( new RKHit( rktrk, rks_cyl, vtdca, ps_cyl, sum ) );
+
+ // Some diagnostics
+ aida.cloud2D( "C r vs z",-1).fill( z_cyl, r_cyl );
+ aida.cloud2D( "Both r vs z",-1).fill( z_cyl, r_cyl);
+
+ // Bookkeeping.
+ nmeas += rks_cyl.rodim;
+ ncmeas += rks_cyl.rodim;
+
+ } else if ( addzp ) {
+
+ // Advance counter to next surface.
+ ++next_zp;
+
+ // Restore the cylinder index to it's position at the start of this step.
+ next_cyl = next_cyl_sav;
+
+
+ // Did we exceed the half-arc limit?
+ sum += path_zp;
+ if ( sum > pathlimit ) break;
+
+ // Track parameters at this step.
+ vtdca = vt_zp;
+
+ // Add the hit.
+ rkl.addHit( new RKHit( rktrk, rks_zp, vtdca, ps_zp, sum ) );
+
+ // Some diagnostics
+ aida.cloud2D( "Z r vs z",-1).fill( z_zp, r_zp );
+ aida.cloud2D( "Both r vs z",-1).fill( z_zp, r_zp );
+
+ // Bookkeeping.
+ nmeas += rks_zp.rodim;
+ nzmeas += rks_zp.rodim;
+
+ } else {
+
+ // Neither intersected but try next surfaces for both lists.
+ ++next_cyl;
+ ++next_zp;
+ }
+
+
+ }
+
+ // More diagnostics.
+ aida.cloud2D("Path length vs cz",-1).fill(rktrk.cz(),sum);
+ aida.cloud2D("C Hits vs cz",-1).fill( rktrk.cz(), ncmeas );
+ 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;
}
+
// Monitoring histograms for the generated track.
private void plotGenTrk( RKTrack t ){
@@ -401,4 +650,64 @@
aida.cloud2D("D0 vs Z0",-1).fill(t.z0(),t.d0());
}
+
+ // Check that the new hit generator matches the old one.
+ // Only valid when scattering and eloss are turned off.
+ void CheckHitList( RKHitList hlist, RKTrack rktrk ){
+
+ RKHitList hlist2 = GenerateOneArcHits( rktrk );
+ int d1 = hlist.nHits() - hlist2.nHits();
+ int d2 = hlist.nCyl() - hlist2.nCyl();
+ int d3 = hlist.nZp() - hlist2.nZp();
+ int d4 = hlist.nDof() - hlist2.nDof();
+
+ if ( (d1==0) && (d2==0) && (d3==0) && (d4 ==0) ){
+ List<RKHit> l1 = hlist.getForward();
+ List<RKHit> l2 = hlist2.getForward();
+ for ( int i=0;i<l1.size(); ++i ){
+ RKHit h1 = l1.get(i);
+ RKHit h2 = l2.get(i);
+ double diff = h2.getPath()-h1.getPath();
+ aida.cloud1D("Difference in Path length: ").fill(diff);
+
+ }
+ }else{
+ System.out.println ( "Diffs: " + d1 + " " + d2 + " " + d3 + " " + d4 + " " + rktrk.cz() );
+ System.out.println ( "Hits1: "
+ + hlist2.nHits() + " "
+ + hlist2.nCyl() + " "
+ + hlist2.nZp() + " "
+ + hlist2.nDof() + " "
+ );
+ System.out.println ( "Hits2: "
+ + hlist.nHits() + " "
+ + hlist.nCyl() + " "
+ + hlist.nZp() + " "
+ + hlist.nDof() + " "
+ );
+ int i=0;
+ for ( RKHit h: hlist.getForward() ){
+ System.out.printf ( "Hits1: %3d %-20s %10.4f %10.4f\n",
+ i++,
+ h.getSurface().getTypeAsString(),
+ h.getSurface().radius,
+ h.getSurface().zc
+ );
+ }
+ i=0;
+ for ( RKHit h: hlist2.getForward() ){
+ System.out.printf ( "Hits2: %3d %-20s %10.4f %10.4f\n",
+ i++,
+ h.getSurface().getTypeAsString(),
+ h.getSurface().radius,
+ h.getSurface().zc
+ );
+ }
+
+ }
+
+ }
+
}
+
+