lcsim/src/org/lcsim/util/step
diff -u -r1.4 -r1.5
--- CamKalmanStepper.java 17 Jul 2007 14:41:22 -0000 1.4
+++ CamKalmanStepper.java 20 Jul 2007 14:08:43 -0000 1.5
@@ -98,6 +98,7 @@
kMass = parm[7] ;
kRungeKutta = true;
set(parm,covar); // create external track parameters from given arguments
+ System.out.println("First value of kRp[0]"+ kRp[0]+" kRp[1]"+kRp[1]+" kRp[2]"+ kRp[2]);
}
/** set rp and covariant Matrix
@@ -135,9 +136,20 @@
ncount1++;
setBField1Dim(b);
double dxyz = xyz-kRp[ixyz];
- int nsteps = (int)(Math.abs(dxyz/0.01)+0.5);
+ if(ixyz==0)
+ System.out.println("x="+xyz+" kRp[0]="+kRp[ixyz]+" ixyz="+ixyz);
+ if(ixyz==0)
+ System.out.println(" new x="+(kRp[0]+dxyz));
+ // int nsteps = (int)(Math.abs(dxyz/0.01)+0.5); debug
+ int nsteps = (int)(Math.abs(dxyz/0.1)+0.5);
if(ncount1<4)System.out.println("In propagateTo dxyz="+dxyz+" nsteps="+nsteps);
- for (int i=0;i<nsteps;i++) stepByDx(dxyz/nsteps);
+ for (int i=0;i<nsteps;i++){
+ if(getStopTkELow()){ System.out.println("not enough energy left");
+ break;}
+ stepByDx(ixyz,dxyz/nsteps);
+ if(ixyz==0) System.out.println("!!!! In Propagate:after stepBy new x ="+kRp[0]+" step number ="+i);
+ }
+
return (ktransport=true);
}
/** stepby ds on the particle pathLength
@@ -149,7 +161,8 @@
double E = Math.sqrt(P*P+kMass*kMass) ;
double k=kB2C*kQ*kB/P ;
double dpdx = getDeDx(P, kMass,"Iron")*E/P*ds ;
-
+
+ if(ncount1<6) System.out.println("I'm in StepByDs");
if(dpdx>P) return stopTkELow=true ;
double rdL= getRadLength("Iron"); //im mm
@@ -160,10 +173,16 @@
{0.,0.,0.,wk[0],0.,0.},{0.,0.,0.,0.,wk[1],0.},{0.,0.,0.,0.,0.,wk[2]}};
Matrix Qk = new Matrix(valsq);
double[][]fifi=new double[6][6];
+ if(ncount1<6) System.out.println("1st value kRp[0]" + kRp[0]+" kRp[3]= " + kRp[3]+" kRp[4]=" + kRp[4]);
+// if(ncount1<6) System.out.println(" 1st value P="+ P+ " ds= "+ds);
kRp[0]+=kRp[3]/P*ds ;
+ // System.out.println("finding out certain variables kRp[3]="+ kRp[3]+ " P="+P+" ds="+ds);
+
kRp[1]+=kRp[4]/P*ds ;
kRp[2]+=kRp[5]/P*ds ;
k*=ds;
+ if(ncount1<6) System.out.print("StepByDs ncount: "+ncount1+ " kRp[0]="+kRp[0]+" kRp[3]= " + kRp[3]+" kRp[4]=" + kRp[4]);
+
kRp[3]+=(k/(1.+0.25*k*k)*kRp[4]-0.5*k*k*kRp[3])-dpdx*kRp[3]/P ;
kRp[4]+=(-k/(1.+0.25*k*k)*kRp[3]-0.5*k*k*kRp[4])-dpdx*kRp[4]/P ;
kRp[5]+=-dpdx*kRp[5]/P;
@@ -184,33 +203,21 @@
fifi[4][4] =-k*kRp[3]*kRp[4]/(P*P) ;
fifi[4][5] =-k*kRp[3]*kRp[5]/(P*P) ;
Matrix fi = new Matrix(fifi,6,6) ;
- if (ncount1<2)System.out.println("Print Matrix fi- ncount1="+ncount1);
- if(ncount1<2)fi.print(6,6);
+ if (ncount1==0)System.out.println("Print Matrix fi- ncount1="+ncount1);
+ if(ncount1==0)fi.print(6,6);
// fi.constructWithCopy(fifi);
fi=fi.times(ds);
Matrix fiMtx = fi.copy();
fi.plus(Matrix.identity(6,6));
//update covariant matrix
Matrix PkMinus= new Matrix(kCovariance,6,6); //
- if(ncount1<2)System.out.println("Print PkMinus- ncount1="+ncount1);
- if(ncount1<2)PkMinus.print(6,6);
+ if(ncount1==0)System.out.println("Print PkMinus- ncount1="+ncount1);
+ if(ncount1==0)PkMinus.print(6,6);
PkMinus =((fiMtx.times(PkMinus)).times(fiMtx.transpose())).plus(Qk);
- if(ncount1<2)System.out.println("Print PkMinus after MS inclusion");
- if(ncount1<2)PkMinus.print(6,6);
+ if(ncount1==0)System.out.println("Print PkMinus after MS inclusion");
+ if(ncount1==0)PkMinus.print(6,6);
return (ktransport=true);
}
-/** stepby dx on the x direction
- *@param dx one step in x
- */
- private boolean stepByDx(double dx){ // one step by dx
- double P = Math.sqrt(kRp[3]*kRp[3]+ kRp[4]*kRp[4]+kRp[5]*kRp[5]);
- if(Math.abs(kRp[3])<1.E-010){
- System.out.println("Zero Px, can't propagate by dx");
- return (stopTkELow=true);
- }
- return stepByDs(dx*P/kRp[3]);
- }
-
/** stepby dxyz on the xyz direction
*@param dxyz one step in xyz
*/
@@ -220,7 +227,7 @@
System.out.println("Zero Pxyz, can't propagate in direction"+ixyz);
return (stopTkELow=true);
}
- return stepByDs(dxyz*P/kRp[3]);
+ return stepByDs(dxyz*P/kRp[(ixyz+3)]);
}
/** stepby dr along the radius direction