lcsim/sandbox/RobKutschke/TRFTests/v1
diff -u -r1.3 -r1.4
--- FullFitKalman.java 18 Sep 2007 16:50:17 -0000 1.3
+++ FullFitKalman.java 11 Oct 2007 14:55:13 -0000 1.4
@@ -16,6 +16,9 @@
import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
import org.lcsim.recon.tracking.trfcyl.ThinCylMs;
+import org.lcsim.recon.tracking.trfzp.SurfZPlane;
+import org.lcsim.recon.tracking.trfzp.ThinZPlaneMs;
+
import org.lcsim.util.aida.AIDA;
@@ -32,9 +35,9 @@
* a single track.
*
*@author $Author: kutschke $
- *@version $Id: FullFitKalman.java,v 1.3 2007/09/18 16:50:17 kutschke Exp $
+ *@version $Id: FullFitKalman.java,v 1.4 2007/10/11 14:55:13 kutschke Exp $
*
- * Date $Date: 2007/09/18 16:50:17 $
+ * Date $Date: 2007/10/11 14:55:13 $
*
*/
@@ -42,6 +45,11 @@
{
private AIDA aida = AIDA.defaultInstance();
+
+ // Flags to control: multiple scattering, energy loss and adding the hit.
+ private boolean doMs = true;
+ private boolean doEloss = true;
+ private boolean doMeas = true;
// static methods
@@ -108,6 +116,14 @@
{ return _pprop; }
//
+
+ public void setDoMs( boolean b){
+ doMs = b;
+ }
+
+ public void setDoEloss( boolean b){
+ doEloss = b;
+ }
/**
*Fit the specified track.
@@ -191,6 +207,11 @@
System.out.println ("Params: " + e_save.vector() );
return icount+1;
}
+
+ //if ( icount != 0 ) {
+ //int istat = interact( trh, hit, dir );
+ //}
+
// Add the hit.
int fstat = _addfit.addHit(trh,hit);
@@ -204,50 +225,7 @@
}
if ( fstat>0 ) return 10000 + 1000*fstat + icount;
- // This fails since the interactor is lost in copying.
- // if ( hit.surface().getInteractor() != null ){
- //System.out.println ("Have interactor: " + hit.surface() );
- //}
-
- if ( hit.surface().pureType().equals(SurfCylinder.staticType()) ){
-
- SurfCylinder s = (SurfCylinder) hit.surface();
- double r = s.radius();
- if ( r > 1.3 && RKDebug.Instance().getDoMs() && r < 100. ){
- TrackError eold = trh.newTrack().error();
- //System.out.println( "Error before: " + eold );
-
- aida.histogram1D("Fit scat radius forward:",300,0.,150.).fill(r);
-
- double l_over_radl = ( r< 10) ? 0.000916 : 0.013747;
- ThinCylMs scat = new ThinCylMs(l_over_radl*RKDebug.Instance().getMsFac());
- ETrack et = trh.newTrack();
- ETrack ets = new ETrack(et);
- double chnew = trh.chisquared();
-
- scat.interact(et);
- hit.update(et);
- trh.setFit(et,chnew);
-
- /*
- for ( int i=0; i<5; ++i ){
- double ex1 = et.error().get(i,i);
- double ex2 = ets.error().get(i,i);
- double sigsq = ex1-ex2;
- //double pull = -9.;
- double sig=-1.;
- if ( sigsq > 0. ){
- sig = Math.sqrt(sigsq);
- //pull = (et.vector(i)-ets.vector(i))/sig;
- }
- aida.cloud1D("Forward Fit Delta param:"+i).fill(et.vector(i)-ets.vector(i));
- if ( sig > 0. ) aida.cloud1D("Forward Fit Delta error:"+i).fill(sig);
- }
- //System.out.println( "Error after: " + trh.newTrack().error().minus(eold) );
- */
- }
-
- }
+ int istat = interact( trh, hit, dir );
++icount;
}
@@ -280,17 +258,8 @@
// Extract the next hit pointer.
Hit hit = (Hit)ihit.next();
- // Propagate to the next surface.
- //System.out.println ("Prop to: " + hit.surface() );
- //System.out.println ("Before prop: \n" );
- // psm.Print( trh.newTrack().error() );
-
PropStat pstat = trh.propagate(_pprop,hit.surface(),dir);
- //System.out.println ("Pstat: " + pstat );
- //System.out.println ("Pstat: " + pstat.success() );
- //System.out.println ("After prop: \n" );
- // psm.Print( trh.newTrack().error() );
if ( ! pstat.success() ) {
System.out.println ("Error: "
+ RKDebug.Instance().getTrack() + " "
@@ -302,12 +271,13 @@
return icount+1;
}
+ //if ( icount != 0 ) {
+ //int istat = interact( trh, hit, dir );
+ //}
+
// Add the hit.
int fstat = _addfit.addHit(trh,hit);
- //System.out.println ("fstat: " + fstat );
- //System.out.println ("After addhit: \n" );
- //psm.Print( trh.newTrack().error() );
if ( fstat>0 ){
System.out.println ("Error: "
@@ -328,61 +298,14 @@
}
chold = chnew;
- if ( hit.surface().pureType().equals(SurfCylinder.staticType()) ){
-
- SurfCylinder s = (SurfCylinder) hit.surface();
- double r = s.radius();
- if ( r > 1.3 && RKDebug.Instance().getDoMs() ){
- TrackError eold = trh.newTrack().error();
-
- aida.histogram1D("Fit scat radius backward:",300, 0., 150.).fill(r);
-
- double l_over_radl = ( r< 10) ? 0.000916 : 0.013747;
- ThinCylMs scat = new ThinCylMs(l_over_radl*RKDebug.Instance().getMsFac());
- ETrack et = trh.newTrack();
- ETrack ets = new ETrack(et);
- double chnew1 = trh.chisquared();
-
- scat.interact(et);
- hit.update(et);
- trh.setFit(et,chnew1);
-
- /*
- for ( int i=0; i<5; ++i ){
- double ex1 = et.error().get(i,i);
- double ex2 = ets.error().get(i,i);
- double sigsq = ex1-ex2;
- //double pull = -9.;
- double sig=-1.;
- if ( sigsq > 0. ){
- sig = Math.sqrt(sigsq);
- //pull = (et.vector(i)-ets.vector(i))/sig;
- }
- //System.out.println ("Delta: "
- // + i + " "
- // + r + " "
- // + (et.vector(i)-ets.vector(i)) + " "
- // + (ex1-ex2)
- // );
- aida.cloud1D("Backward Fit Delta param:"+i).fill(et.vector(i)-ets.vector(i));
- aida.cloud1D("Backward Fit Delta error:"+i).fill(sig);
-
- }
- //System.out.println( "Error after: " + trh.newTrack().error().minus(eold) );
- */
- }
- }
- /*
- ETrack et = trh.newTrack();
- hit.surface().interact(et);
- hit.update(et);
- trh.setFit(et,chnew);
- */
+ // Apply multiple scattering and energy loss
+ int istat = interact( trh, hit, dir );
++icount;
}
+ // Propagate to the PCAO.
SurfDCA sdca = new SurfDCA( 0., 0. );
PropStat pstat = trh.propagate( _pprop, sdca,dir );
if ( ! pstat.success() ) return -1;
@@ -390,7 +313,75 @@
return 0;
}
-
+
+ private int interact ( HTrack trh, Hit hit, PropDir dir ){
+ if ( hit.surface().pureType().equals(SurfCylinder.staticType()) ){
+
+ SurfCylinder s = (SurfCylinder) hit.surface();
+ double r = s.radius();
+ if ( r > 1.3 && doMs && r < 100. ){
+ TrackError eold = trh.newTrack().error();
+ //System.out.println( "Error before: " + eold );
+
+ aida.histogram1D("Fit scat radius forward:",300,0.,150.).fill(r);
+
+ double l_over_radl = ( r< 10) ? 0.000916 : 0.013747;
+ ThinCylMs scat = new ThinCylMs(l_over_radl*RKDebug.Instance().getMsFac());
+ ETrack et = trh.newTrack();
+ ETrack ets = new ETrack(et);
+ double chnew = trh.chisquared();
+
+ scat.interact(et);
+ hit.update(et);
+ trh.setFit(et,chnew);
+
+ /*
+ for ( int i=0; i<5; ++i ){
+ double ex1 = et.error().get(i,i);
+ double ex2 = ets.error().get(i,i);
+ double sigsq = ex1-ex2;
+ //double pull = -9.;
+ double sig=-1.;
+ if ( sigsq > 0. ){
+ sig = Math.sqrt(sigsq);
+ //pull = (et.vector(i)-ets.vector(i))/sig;
+ }
+ aida.cloud1D("Forward Fit Delta param:"+i).fill(et.vector(i)-ets.vector(i));
+ if ( sig > 0. ) aida.cloud1D("Forward Fit Delta error:"+i).fill(sig);
+ }
+ //System.out.println( "Error after: " + trh.newTrack().error().minus(eold) );
+ */
+ } // end if MS enabled
+
+ } // end CYlinder MS
+
+ if ( hit.surface().pureType().equals(SurfZPlane.staticType()) ){
+
+ SurfZPlane s = (SurfZPlane) hit.surface();
+ double z = s.z();
+ if ( doMs && Math.abs(z) < 165.2 ){
+ TrackError eold = trh.newTrack().error();
+
+ aida.histogram1D("Fit scat z forward:",300,-150,150.).fill(z);
+
+ double l_over_radl = ( Math.abs(z)< 25) ? 0.000916 : 0.013747;
+ ThinZPlaneMs scat = new ThinZPlaneMs(l_over_radl*RKDebug.Instance().getMsFac());
+ ETrack et = trh.newTrack();
+ ETrack ets = new ETrack(et);
+ double chnew = trh.chisquared();
+
+ scat.interact(et);
+ hit.update(et);
+ trh.setFit(et,chnew);
+
+ } // end if MS enabled.
+
+ } // end ZPlane MS
+
+ // Successful return;
+ return 0;
+ }
+
/**
*output stream
*
@@ -402,3 +393,4 @@
}
}
+
