lcsim/sandbox/RobKutschke/TRFTests/v1
diff -u -r1.4 -r1.5
--- FitTest_sid00Driver.java 18 Sep 2007 16:50:17 -0000 1.4
+++ FitTest_sid00Driver.java 11 Oct 2007 16:54:53 -0000 1.5
@@ -1,5 +1,8 @@
import java.util.ArrayList;
+import java.util.Iterator;
import java.util.List;
+import java.util.Set;
+import Jama.Matrix;
import org.lcsim.event.EventHeader;
@@ -23,10 +26,23 @@
import org.lcsim.recon.tracking.trfcyl.ThinCylMs;
import org.lcsim.recon.tracking.trfcyl.ThinCylMsSim;
+import org.lcsim.recon.tracking.trfcyl.CylEloss;
+import org.lcsim.recon.tracking.trfcyl.CylElossSim;
+
+import org.lcsim.recon.tracking.trfzp.ThinZPlaneMs;
+import org.lcsim.recon.tracking.trfzp.ThinZPlaneMsSim;
+
+import org.lcsim.recon.tracking.trfeloss.DeDxBethe;
import org.lcsim.recon.tracking.trffit.HTrack;
import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.VecOp;
+
+import org.lcsim.event.EventHeader.LCMetaData;
+import org.lcsim.event.SimTrackerHit;
+
/**
* Driver for exercises learning how to use the TRF code.
@@ -34,9 +50,9 @@
* no multiple scattering or energy loss.
*
*@author $Author: kutschke $
- *@version $Id: FitTest_sid00Driver.java,v 1.4 2007/09/18 16:50:17 kutschke Exp $
+ *@version $Id: FitTest_sid00Driver.java,v 1.5 2007/10/11 16:54:53 kutschke Exp $
*
- * Date $Date: 2007/09/18 16:50:17 $
+ * Date $Date: 2007/10/11 16:54:53 $
*
*/
@@ -70,7 +86,17 @@
private double scatThick1 = 0.000916;
private double scatThick2 = 0.013747;
- private boolean doms = true;
+ private ThinZPlaneMsSim zscatSim1 = null;
+ private ThinZPlaneMsSim zscatSim2 = null;
+
+ private CylElossSim celoss0 = null;
+ private CylElossSim celoss1 = null;
+ private CylElossSim celoss2 = null;
+ private CylElossSim zeloss1 = null;
+ private CylElossSim zeloss2 = null;
+
+ private boolean doms = true;
+ private boolean doeloss = false;
public FitTest_sid00Driver(){
}
@@ -83,23 +109,23 @@
// Limits for track generation parameters.
param = new RKTrackGenParams();
param.setbz(rkgeom.getBz());
- param.setptmin( 1.0 );
+ //param.setptmin( 1.0 );
+ //param.setptmax( 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.65);
+ //param.setczmax( 0.93);
//param.setczmin( 0.);
//param.setczmax( 0.);
+ //param.setczmin( 0.7);
//param.setczmax( -0.93);
//param.setczmin( -0.97);
//param.setczmax( -0.93);
- param.setczmin( -0.6);
- param.setczmax( 0.6);
- //param.setczmin( -0.995);
- //param.setczmax( 0.995);
-
+ //param.setczmin( -0.6);
+ //param.setczmax( 0.6);
+ param.setczmin( -0.993);
+ param.setczmax( 0.993);
+
//param.setd0min( 0.0);
//param.setd0max( 0.0);
//param.setz0min( 0.0);
@@ -107,16 +133,21 @@
System.out.println (param);
- // Seed for track random number generator
- long seed = -987835122;
+ // Seed for track random number generator.
+ long seed = -987835122;
+
+ // An increment used to generate independent streams of random
+ // numbers - I have no proof that this really works.
+ long seedIncrement = 876633217;
+
gen = new RKTrackGen(param,seed);
// Seed for hit generator.
- seed += 876633217;
+ seed += seedIncrement;
RKHit.setSeed(seed);
// Seed for smearing covariance matrices.
- seed += 876633217;
+ seed += seedIncrement;
covsmear = new RKCovSmear( seed );
// Intial covariance matrices.
@@ -143,12 +174,18 @@
vstart_out.set(4,4,sig);
RKDebug.Instance().setDoMs(doms);
+ RKDebug.Instance().setDoEloss(doeloss);
RKDebug.Instance().setMsFac(1.0);
+ RKDebug.Instance().setPrintOutliers(true);
+
scatSim0 = new ThinCylMsSim(scatThick0*RKDebug.Instance().getMsFac());
scatSim1 = new ThinCylMsSim(scatThick1*RKDebug.Instance().getMsFac());
scatSim2 = new ThinCylMsSim(scatThick2*RKDebug.Instance().getMsFac());
+ zscatSim1 = new ThinZPlaneMsSim(scatThick1*RKDebug.Instance().getMsFac());
+ zscatSim2 = new ThinZPlaneMsSim(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();
@@ -156,18 +193,48 @@
vtmp.set(1,0.0);
vtmp.set(2,0.0);
vtmp.set(3,0.0);
- vtmp.set(3,1.0);
+ vtmp.set(4,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);
- }
+ // Seed the random number generators in the Siminteractors.
+ seed += seedIncrement;
+ scatSim0.setSeed(seed);
+ seed += seedIncrement;
+ scatSim1.setSeed(seed);
+ seed += seedIncrement;
+ scatSim2.setSeed(seed);
+ seed += seedIncrement;
+ zscatSim1.setSeed(seed);
+ seed += seedIncrement;
+ zscatSim2.setSeed(seed);
+
+ double densityBe = rkgeom.getCylinders().get(0).density;
+ double densitySi = rkgeom.getCylinders().get(1).density;
+ DeDxBethe dedxBe = new DeDxBethe(densityBe);
+ DeDxBethe dedxSi = new DeDxBethe(densitySi);
+ System.out.println ("Densities: " + densityBe + " " + densitySi );
+
+ double thickbeampipe = rkgeom.getCylinders().get(0).thick;
+ double thickpixel = rkgeom.getCylinders().get(1).thick;
+ double thickstrip = rkgeom.getCylinders().get(6).thick;
+
+ celoss0 = new CylElossSim( thickbeampipe, dedxBe );
+ celoss1 = new CylElossSim( thickpixel, dedxSi );
+ celoss2 = new CylElossSim( thickstrip, dedxSi );
+
+ zeloss1 = new CylElossSim( thickpixel, dedxSi );
+ zeloss2 = new CylElossSim( thickstrip, dedxSi );
+
+ seed += seedIncrement;
+ celoss0.setSeed(seed);
+ seed += seedIncrement;
+ celoss1.setSeed(seed);
+ seed += seedIncrement;
+ celoss2.setSeed(seed);
+ seed += seedIncrement;
+ zeloss1.setSeed(seed);
+ seed += seedIncrement;
+ zeloss2.setSeed(seed);
}
@@ -185,13 +252,66 @@
RKDebug.Instance().setTrack(i);
DoIt(i);
}
+ RKDebug.Instance().setTrack(-99);
+
+ // This is the code to debug the zplane Ms bug.
+ //HackTest();
+
+ /*
+ List<List<SimTrackerHit>> simTrackerHitCollections = event.get(SimTrackerHit.class);
+ for ( List<SimTrackerHit> simTrackerHits : simTrackerHitCollections )
+ {
+ LCMetaData meta = event.getMetaData(simTrackerHits);
+ System.out.println ("Hits: " + meta.getName() + " " + simTrackerHits.size() );
+ }
+
+ Set keys = event.keys();
+ System.out.println ("Keys: \n" + keys );
+ System.out.println ("Event: \n" + event.getClass());
+ System.out.println ("Event: \n" + event);
+
+
+ for ( Iterator i=keys.iterator(); i.hasNext(); ){
+ String s = (String) i.next();
+ System.out.println (s);
+ }
+
+
+ double densityBe = rkgeom.getCylinders().get(0).density;
+
+ DeDxBethe dedxBe = new DeDxBethe(densityBe);
+ double xl = 0.2;
+ double xh = 10.2;
+ int nn = 501;
+ double dx = (xh-xl)/nn;
+ double len=0.1;
+
+ for ( int i=0; i<nn; ++i){
+ double p0 = xl + (i+0.5)*dx;
+ double p = dedxBe.loseEnergy(p0,len);
+ //double dp = dedxBe.dEdX( p);
+ double dp = p0-p;
+
+ aida.cloud2D("Bethe-Block test:",-1).fill(p0,dp);
+
+ double p2 = dedxBe.loseEnergy(p,-len);
+ double dp2 = p2-p0;
+ aida.cloud2D("Bethe-Block Round Trip Error",-1).fill(p0,dp2);
+ double sig = dedxBe.sigmaEnergy(p0,len);
+ aida.cloud2D("Bethe-Block Sigma",-1).fill(p0,sig);
+ aida.cloud2D("Bethe-Block test:",-1).fill(p0,dp+sig);
+ aida.cloud2D("Bethe-Block test:",-1).fill(p0,dp-sig);
+
+
+ }
+ */
}
private void DoIt( int itrk){
- System.out.println ("New track: " + itrk );
+ //System.out.println ("New track: " + itrk );
// Generate a track.
RKTrack rktrk = gen.newRKTrack();
@@ -218,6 +338,8 @@
// Instantiate fitter.
Propagator prop = rkgeom.newPropagator();
FullFitKalman fitk = new FullFitKalman(prop);
+ fitk.setDoMs(doms);
+ fitk.setDoEloss(doeloss);
// Get starting track parameters for outward fit ( at PCAO ).
toTRF trftrk = new toTRF(rktrk);
@@ -334,8 +456,8 @@
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());
+ //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.
@@ -365,7 +487,6 @@
RKHitList rkl = new RKHitList();
Propagator prop = rkgeom.newPropagator();
double sum = 0.;
- //System.out.println("Start dca: " + vtdca +"\n" );
for ( RKSurf s: rkgeom.getCylinders() ){
PropStat ps = prop.vecDirProp( vtdca, s.getCylinder(), PropDir.FORWARD);
if ( ps.success() ){
@@ -434,7 +555,6 @@
toTRF trftrk = new toTRF(rktrk);
VTrack vtdca = trftrk.atDcaZaxis();
-
// This code only looks at the first arc of a track.
double pathlimit = trftrk.halfarc();
@@ -481,6 +601,7 @@
int next_zp_sav = next_zp;
// Check cylindrical surfaces until the next good one is found.
+ //while ( next_cyl < -1 ){
while ( next_cyl < max_cyl ){
// Default status is failure.
@@ -596,6 +717,7 @@
} else {
// This should only happen when both lists are exhausted.
+ // Or when one of the hit types is turned off.
if ( next_cyl<max_cyl || next_zp < max_zp ){
System.out.println ("Don't know how I got here ... "
+ next_cyl + " "
@@ -636,80 +758,84 @@
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. ){
+ // For now, skip the beam pipe since we have not finished
+ // the bookkeeping for scattering surfaces with no hits.
+ if ( r_cyl > 1.3 ){
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);
+ 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() );
- }
+ aida.histogram1D("Gen scat radius", 300, 0., 150.).fill(r_cyl);
}
+ }
+
+ if ( r_cyl > 1.3 && doeloss ){
+ // Save the state before eloss
+ VTrack before = new VTrack(vtdca);
+ VTUtil vtu_before = new VTUtil ( before);
+ // Perform eloss
+ if ( r_cyl < 10. ){
+ celoss1.interact(vtdca);
+ }else{
+ celoss2.interact(vtdca);
+ }
+ VTUtil vtu_after = new VTUtil ( vtdca);
+
+ // Diagnostics
+ //if ( r_cyl < 10. ){
+ if ( r_cyl < 2. ){
+ CylEloss ce = new CylEloss( celoss1.thickness(), celoss1.dEdX() );
+
+ TrackError vtmp = new TrackError();
+ ETrack et = new ETrack( before.surface().newPureSurface(), before.vector(), vtmp);
+ ETrack ets = new ETrack(et);
+ ce.interact_dir(et,PropDir.FORWARD);
+ VTUtil vtu2 = new VTUtil ( et );
+
+ if ( vtdca.vector(4) > 0. ){
+ aida.cloud1D("Eloss wrt orig q=+1").fill( vtdca.vector(4)-before.vector(4));
+ aida.cloud1D("Eloss wrt new q=+1" ).fill( vtdca.vector(4)-et.vector(4));
+ } else{
+ aida.cloud1D("Eloss wrt orig q=-1").fill( vtdca.vector(4)-before.vector(4));
+ aida.cloud1D("Eloss wrt new q=-1" ).fill( vtdca.vector(4)-et.vector(4));
+ }
+ aida.cloud1D("Delta E base").fill(vtu2.p() - vtu_before.p());
+ aida.cloud1D("Delta E straggle").fill(vtu_after.p() - vtu2.p());
+ System.out.println (" Energy: "
+ + (vtu2.p()-vtu_before.p())
+ );
+ //System.out.println ( et.error() );
+
+ for ( int i=0; i<5; ++i ){
+ if ( vtdca.vector(4) > 0. ){
+ aida.cloud1D("Cyl Eloss q=+1 "+i).fill(vtdca.vector(i)-before.vector(i));
+ } else{
+ aida.cloud1D("Cyl Eloss q=-1 "+i).fill(vtdca.vector(i)-before.vector(i));
+ }
+ if ( et.error(i,i) > 0. ){
+ double sig = Math.sqrt(et.error(i,i));
+ double pull = (vtdca.vector(i)-before.vector(i))/sig;
+ aida.cloud1D("Cyl Pull " + i).fill(pull);
+ aida.cloud1D("Cyl Sigma " +i).fill(sig*1.e6);
+ }
+ }
+ }
}
} else if ( addzp ) {
@@ -720,7 +846,6 @@
// 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;
@@ -739,6 +864,20 @@
nmeas += rks_zp.rodim;
nzmeas += rks_zp.rodim;
+ if ( doms ){
+
+ // Interact.
+ VTrack test = new VTrack(vtdca);
+ if ( Math.abs(z_zp) < 25. ){
+ zscatSim1.interact(vtdca);
+ }else{
+ zscatSim2.interact(vtdca);
+ }
+ aida.histogram1D("Gen scat z", 340, -170., 170.).fill(z_zp);
+
+ }
+
+
} else {
// Neither intersected but try next surfaces for both lists.
@@ -850,6 +989,78 @@
}
-}
+ private void HackTest(){
+ double pp = 1.0;
+ TrackVector v0 = new TrackVector();
+ v0.set(0,0.0);
+ v0.set(1,0.0);
+ v0.set(2,0.2);
+ v0.set(3,0.0);
+ v0.set(4,1./pp);
+
+ SurfZPlane s = new SurfZPlane(7.6);
+ VTrack vt0 = new VTrack( s, v0 );
+ vt0.setForward();
+
+
+ VTUtil vu0 = new VTUtil(vt0);
+ Hep3Vector p = vu0.asHep3Vector();
+ Hep3Vector phat0 = VecOp.unit(p);
+
+ Hep3Vector n = new BasicHep3Vector(0.,0.,1.);
+ Hep3Vector u = new BasicHep3Vector(phat0.v());
+ Hep3Vector v = VecOp.unit(VecOp.cross(n,u));
+ Hep3Vector w = VecOp.unit(VecOp.cross(p,v));
+
+ double radl = scatThick2*10.;
+ ThinZPlaneMsSim sim = new ThinZPlaneMsSim(radl);
+ ThinZPlaneMs scat = new ThinZPlaneMs(radl);
+
+ double ltrue = radl/vu0.costh();
+
+ double rms = (0.0136)*Math.sqrt(ltrue)*(1.+0.038*Math.log(ltrue))/pp;
+
+ for ( int j=0; j<1000; ++j){
+ VTrack vt = new VTrack( vt0 );
+ vt.setForward();
+ sim.interact(vt);
+
+ double ccz = 1./Math.sqrt( 1. + vt.vector(2)*vt.vector(2) + vt.vector(2)*vt.vector(3) );
+ double ccx = vt.vector(2)*ccz;
+ double ccy = vt.vector(3)*ccz;
+ double nnn = Math.sqrt( ccx*ccx + ccy*ccy + ccz*ccz );
+
+ VTUtil vu = new VTUtil(vt);
+ Hep3Vector p1 = vu.asHep3Vector();
+ Hep3Vector phat = VecOp.unit(p1);
+
+ double cv = VecOp.dot(phat,v);
+ double cw = VecOp.dot(phat,w);
+
+ double tv = Math.asin(cv);
+ double tw = Math.asin(cw);
+ aida.cloud1D("cv").fill(cv);
+ aida.cloud1D("cw").fill(cw);
+ aida.cloud1D("tv").fill(tv);
+ aida.cloud1D("tw").fill(tw);
+ aida.cloud1D("tv pull").fill(tv/rms);
+ aida.cloud1D("tw pull").fill(tw/rms);
+
+ TrackError vtmp = new TrackError();
+ ETrack et = new ETrack( vt0.surface().newPureSurface(), vt0.vector(), vtmp);
+ scat.interact(et);
+
+ for ( int i=0;i<5; ++i){
+ double diff = (vt.vector(i) - vt0.vector(i));
+ aida.cloud1D("Plane scat test param: " + i,-1).fill(diff);
+ if ( et.error(i,i) >0. ){
+ double sig = Math.sqrt(et.error(i,i));
+ double pull = diff/sig;
+ aida.histogram1D("Plane scat pull test param: " + i,200,-10.,10.).fill(pull);
+ }
+ }
+ }
+ }
+}