Print

Print


Commit in lcsim/src/org/lcsim/util/step on MAIN
CamKalmanStepper.java+29-221.4 -> 1.5
Kudzanayi Munetsi-Mugomba 07-20-07 corrections in PropagateTo

lcsim/src/org/lcsim/util/step
CamKalmanStepper.java 1.4 -> 1.5
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
CVSspam 0.2.8