lcsim/src/org/lcsim/util/step
diff -N CamKalmanStepper2.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ CamKalmanStepper2.java 28 Nov 2007 21:50:19 -0000 1.1
@@ -0,0 +1,497 @@
+package org.lcsim.util.step;
+import org.lcsim.geometry.compact.converter.lcdd.util.*;
+import org.lcsim.units.clhep.PhysicalConstants.*;
+import java.util.*;
+import java.text.*;
+import org.lcsim.util.aida.*;
+import Jama.*;
+import Jama.util.*;
+import org.lcsim.recon.muon.CoilSubLayer;
+import org.lcsim.recon.cluster.util.CalHitMapMgr;
+import org.lcsim.geometry.subdetector.CylindricalBarrelCalorimeter;
+import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.subdetector.*;
+import org.lcsim.geometry.field.*;
+import org.lcsim.geometry.layer.*;
+import org.lcsim.geometry.Layered;
+import org.lcsim.geometry.LayeredSubdetector;
+import org.lcsim.material.*;
+ /**
+ * CamKalmanStepper2- A Kalman Stepper
+ * trasnsported from JAS2
+ *JINST - published October 2006
+ *arXiv: physics/0604197
+ *@Author-C. Milstene <[log in to unmask]>
+ * Code developped with Fedor Ignatov
+ * and Hans Wenzel consulting
+ * More comments and references to come
+ * Hard-coded material Iron for the dE/dx
+ * (Testing program-TestKalm.java contains
+ * a main- and is the (Driver/MainProject))
+ */
+public class CamKalmanStepper2{
+ private boolean Debug=false;
+ private final double kAlmost1=0.999;
+ private final double kAlmost0=1E-33;
+ private final double kdistsmall=0.001;
+ private final double kVeryBig =(1./kAlmost0);
+ private final double kB2C=0.299979245*1.E-03;
+ private final double kAlmost0Field =1.E-33;
+ private final double kMostProbablePt=0.35;
+ private final double cmTomm=10.;
+ private final double mevTogev=1./1000.;
+ private final double gevTomev=1000.;
+ private double [] zRhoZoverA;
+ private double kQ;
+ private double kB=5.; // default- reset in propagateTo
+ private double kMass;
+ private String MaterialName;
+ private boolean kRungeKutta;
+ public boolean kupDate = false;
+ private static boolean ktransport=true;
+ private static boolean stopTkELow=false;
+ private double [] kRp = new double [6];
+ private double [] kRPlus ;
+ private double [] measState;
+ private double [][]measCov=new double[3][3];
+ private double [] nwRp = new double[8];
+ private Matrix PkMinus = new Matrix(6,6);
+ private Matrix PkPlus = new Matrix(6,6) ;
+ private Matrix Hk ;
+ private double[][] kCovariance = new double [6][6];// construct a 6 by 6 array of zeroes
+ private Matrix fi = new Matrix(6,6); // construct
+ private AIDA aida = AIDA.defaultInstance();
+ private CalHitMapMgr evtptr = null;
+ private boolean StopTkElow=false;
+ private Random rndcm = new Random();
+ private Detector det ;
+ private static int ncount1=0;
+ /** default constructor*/
+ public CamKalmanStepper2(){
+ evtptr = CalHitMapMgr.getInstance();
+ kQ = -1. ;
+ kMass = 0.105;
+ kRungeKutta = true;
+ for(int i=0; i<6; i++)kRp[i]=0.;
+ for(int i=0; i<6; i++){
+ for(int j=0; j<6; j++)kCovariance[i][j]=0.; // or any starting value
+ }
+ }
+ /** Copy Constructor-
+ *@param rp an Array with the particle position and Momnetum
+ * (x,y,z,px,py,pz) - a phase-space point
+ */
+ public CamKalmanStepper2(double [] rp){
+ evtptr = CalHitMapMgr.getInstance();
+ kQ = rp[6] ;
+ kMass = rp[7] ;
+ kRungeKutta = true;
+ for(int i=0; i<6; i++)this.kRp[i]=rp[i];
+ for(int i=0; i<6; i++){
+ for(int j=0; j<6; j++)this.kCovariance[i][j]=kCovariance[i][j];
+ }
+ }
+ /** Copy Constructor-
+ *@param rp an array with the particle position and momentum(x,y,z,px,py,pz)
+ *@param cova a covariant matrix
+ */
+ public CamKalmanStepper2(double [] parm, double [][]covar){ // Copy constructor
+ evtptr = CalHitMapMgr.getInstance();
+ kQ = parm[6] ;
+ kMass = parm[7] ;
+ kRungeKutta = true;
+ set(parm,covar); // create external track parameters from given arguments
+ if(Debug)
+ System.out.println("First value of kRp[0]"+ kRp[0]+" kRp[1]"+kRp[1]+" kRp[2]"+ kRp[2]);
+ }
+/** set rp and covariant Matrix
+ *@param rp an array with the particle position and momentum
+ *@param covar the covariant matrix
+ */
+ public void set(double [] rp, double [] [] covar){
+ for(int i=0; i<6; i++)this.kRp[i]=rp[i];
+ for(int i=0; i<6; i++){
+ for(int j=0;j<6;j++){
+ kCovariance[i][j] = covar[i][j];
+ }
+ }
+ PkMinus= new Matrix(kCovariance,6,6);
+ System.out.println(" PkMinus "); // NOW
+ PkMinus.print(6,6); // NOW
+ }
+ /** Propagate the particle from a position x
+ * to s new position in a magnetic field b
+ *@param x the position
+ *@param b the magnetic field
+ */
+ public boolean propagateTo(double r, double b){
+ setBField1Dim(b);
+ double dr = r-Math.sqrt(kRp[0]*kRp[0]+kRp[1]*kRp[1]);
+ int nsteps = (int)(Math.abs(dr/0.005)+0.5);
+ for (int i=0;i<nsteps;i++) stepByDr(dr/nsteps);
+ return (ktransport=true);
+}
+ /** Propagate the particle from a position
+ * xyz defined either by x, y, or z
+ * to s new position in a magnetic field b
+ *@param xyz the position
+ *@param b the magnetic field
+ */
+
+public boolean propagateTo(int ixyz,double xyz, double b){
+ ncount1++;
+ setBField1Dim(b);
+ double dxyz = xyz-kRp[ixyz];
+ if(ixyz==0&& Debug)
+ System.out.println("x="+xyz+" kRp[0]="+kRp[ixyz]+" ixyz="+ixyz);
+ if((ixyz==0)&& (ncount1<5)&&Debug)
+ 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.2)+0.5);
+ if(Debug)
+ System.out.println("In propagateTo dxyz="+dxyz+" nsteps="+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&&Debug) System.out.println("!!!! In Propagate:after stepBy new x ="+kRp[0]+" step number ="+i);
+ }
+ return (ktransport=true);
+}
+ /** stepby ds on the particle pathLength
+ *@param ds one step length
+ */
+
+
+ public boolean stepByDs(double ds){ // one step by ds
+ double[]wk=new double[3];
+ double P = Math.sqrt(kRp[3]*kRp[3]+ kRp[4]*kRp[4]+kRp[5]*kRp[5]) ;
+ double E = Math.sqrt(P*P+kMass*kMass) ;
+ double k=kB2C*kQ*kB/P ;
+ double dpdx = getDeDx(P, kMass,"Iron")*E/P*ds ;
+ // double dpdx=0.;
+ if(Debug)
+ System.out.println("I'm in StepByDs"+" k="+k);
+ if(dpdx>P){
+ System.out.println(" dpdx> p -track stopped not enough energy");
+ System.out.println(" dpdx="+dpdx+" P="+P);
+ return stopTkELow=true ;
+ }
+ if(ncount1<6) System.out.println("I'm in StepByDs");
+ if(dpdx>P) return stopTkELow=true ;
+ double rdL= getRadLength("Iron"); //im mm
+ if(Debug)
+ System.out.println(" number 1: dpdx="+ dpdx);
+ wk=vwk((getRadLength("Iron"))*Math.abs(ds));
+ if(Debug)
+ System.out.println("wk[0]="+wk[0]+" wk[1]="+wk[1]+" wk[2]="+wk[2]);
+ //Next to be checked again
+
+ double[] [] valsq={{0.,0.,0.,0.,0.,0.},{0.,0.,0.,0.,0.,0.},{0.,0.,0.,0.,0.,0.},
+ {0.,0.,0.,(wk[0]*wk[0]),0.,0.},{0.,0.,0.,0.,(wk[1]*wk[1]),0.},{0.,0.,0.,0.,0.,(wk[2]*wk[2])}};
+ // MSk is the multiple scattering
+ Matrix MSk = new Matrix(valsq);
+ double[][]fifi=new double[6][6];
+ if(Debug)
+ System.out.println("1st value kRp[0]" + kRp[0]+" kRp[3]= " + kRp[3]+" kRp[4]=" + kRp[4]+ " ds="+ds);
+ if(Debug)
+ System.out.println(" number 2: dpdx="+ dpdx);
+ kRp[0]+=(kRp[3]/P)*ds ;
+ kRp[1]+=(kRp[4]/P)*ds ;
+ kRp[2]+=(kRp[5]/P)*ds ;
+ P=Math.sqrt(kRp[3]*kRp[3]+ kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ // 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);
+ // System.out.println("finding out certain variables kRp[3]="+ kRp[3]+ " P="+P+" ds="+ds);
+ System.out.println(" !!!!!!! HELLO FEDOR !!!!");
+ k*=ds;
+ double pxx = kRp[3];
+ double pyy = kRp[4];
+ kRp[3]+=(k/(1.+0.25*k*k/2))*pyy-0.5*k*k*pxx-dpdx*pxx/P ;
+ kRp[4]+=(-k/(1.+0.25*k*k/2))*pxx-0.5*k*k*pyy-dpdx*pyy/P ;
+ double p2=Math.sqrt(kRp[3]*kRp[3]+ kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ if(ncount1<6) System.out.print("StepByDs ncount: "+ncount1+ " kRp[0]="+kRp[0]+" kRp[3]= " + kRp[3]+" kRp[4]=" + kRp[4]);
+ kRp[5]+=-dpdx*kRp[5]/P;
+ // System.out.println(" p2-P="+(p2-P)+" dpdx="+dpdx+" p2="+p2);
+ fifi[0][3] =-kRp[3]*kRp[3]/(P*P*P)+1./P ;
+ fifi[0][4] =-kRp[3]*kRp[4]/(P*P*P) ;
+ fifi[0][5] =-kRp[3]*kRp[5]/(P*P*P) ;
+ fifi[1][3] =-kRp[4]*kRp[3]/(P*P*P) ;
+ fifi[1][4] =-kRp[4]*kRp[4]/(P*P*P)+1./P ;
+ fifi[1][5] =-kRp[4]*kRp[5]/(P*P*P) ;
+ fifi[2][3] =-kRp[5]*kRp[3]/(P*P*P) ;
+ fifi[2][4] =-kRp[5]*kRp[4]/(P*P*P) ;
+ fifi[2][5] =-kRp[5]*kRp[3]/(P*P*P)+1./P ;
+ fifi[3][3] =-k*kRp[4]*kRp[3]/(P*P) ;
+ fifi[3][4] =-k*kRp[4]*kRp[4]/(P*P)+k ;
+ fifi[3][5] =-k*kRp[4]*kRp[5]/(P*P) ;
+ fifi[4][3] =k*kRp[3]*kRp[3]/(P*P)-k ;
+ 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) ;
+ ncount1++;
+ if ((ncount1==1)&& Debug)System.out.println("Print Matrix fi- ncount1="+ncount1);
+ if(ncount1==0)fi.print(6,6);
+ // fi.constructWithCopy(fifi);
+// fi=fi.times(ds); Just Checking
+ Matrix fiMtx = fi.copy();
+ Matrix II = Matrix.identity(6,6);
+ fi=fi.plus(II);
+ if(Debug){
+ System.out.println("Matrix (fi*ds+I(6,6))");
+ fi.print(6,6);
+ }
+ // update PkMinus
+ if(Debug){
+ System.out.println("1/Print PkMinus reenter");
+ PkMinus.print(6,6);
+ }
+ //PkMinus =((fiMtx.times(PkMinus)).times(fiMtx.transpose())).plus(MSk);
+ if(Debug)System.out.println("2/ fi*PkMinus");
+ PkMinus =((fi.times(PkMinus)));
+ if(Debug)System.out.println(" fi*Pkminus");
+ if(Debug)PkMinus.print(6,6);
+ if(Debug)System.out.println( "3/((fi.times(PkMinus)).times(fi.transpose())) ");
+ PkMinus = PkMinus.times(fi.transpose());
+ if(Debug)PkMinus.print(6,6);
+ if(Debug)System.out.println(" MSk");
+ if(Debug)MSk.print(6,6);
+ if(Debug)System.out.println("PkMinus+MSk");
+ // PkMinus = PkMinus.plus(MSk);
+ if(Debug)PkMinus.print(6,6);
+ return (ktransport=true);
+}
+ /** stepby dxyz on the xyz direction
+ *@param dxyz one step in xyz
+ */
+ private boolean stepByDx(int ixyz,double dxyz){ // one step by dx
+ double P = Math.sqrt(kRp[3]*kRp[3]+ kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ double dstoo=Math.abs(dxyz*P/kRp[(ixyz+3)]);
+ if(Math.abs(kRp[ixyz+3])<1.E-010){
+ System.out.println("Zero Pxyz, can't propagate in direction"+ixyz);
+ return (stopTkELow=true);
+ }
+ if(dstoo<=kdistsmall)return(ktransport=true);
+ return stepByDs(dstoo);
+ }
+
+/** stepby dr along the radius direction
+ *@param dr one step in r
+ */
+ private boolean stepByDr(double dr){ // one step by dr
+ double P = Math.sqrt(kRp[3]*kRp[3]+ kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ double a = kRp[3]*kRp[3]+kRp[4]*kRp[4]/P/P;
+ double b = kRp[0]*kRp[3]+kRp[1]*kRp[4]/P ;
+ double c =-(dr*dr+2*dr*Math.sqrt(kRp[0]*kRp[0]+kRp[1]*kRp[1]));
+ double ds=(-b+Math.sqrt(b*b-a*c))/a ;
+ return stepByDs(ds);
+ }
+/** update the phase space point and covariant Matrix
+ *Using the measured position and errors
+ *@param mstate a 3dim array with the measure position
+ *@param mCov the measurement error Matrix
+ */
+ public boolean upDate(double []mState,double[][]mCov){
+ // update vector with current extrapolation
+ double [] cState= new double[6];
+ for(int i=0;i<6;i++) {cState[i]=kRp[i];}
+ // measured position x,y,z
+ measState = new double [3];
+ double [][]val1 = {{1.,0.,0.,0.,0.,0.},{0.,1.,0.,0.,0.,0.},
+ {0.,0.,1.,0.,0.,0.}};
+ Hk=new Matrix(val1);
+ if((ncount1==1)&& Debug)System.out.println( "print HK Matrix");
+ if((ncount1==1)&&Debug)Hk.print(3,6);
+ // measurement covariante Matrix
+ for(int i=0;i<3;i++) measCov[i][i]=mCov[i][i];
+ /** update error covariance Matrix- Kalman Gain Matrix
+ *Problem here !!!!!!
+ */
+ //Matrix KkPart = (Hk.times(PkMinus).times(Hk.transpose()));
+ Matrix KkPart= Hk.times(PkMinus);
+ if(Debug){
+ System.out.println(" Hk.times(PkMinus)");
+ KkPart.print(3,6);
+ }
+ KkPart=KkPart.times(Hk.transpose());
+ if(Debug){
+ System.out.println(" Hk.times(PkMinus).times(Hk.transpose()");
+ KkPart.print(3,3);}
+ Matrix Meas=new Matrix(measCov,3,3);
+ if((ncount1<4)&&Debug)System.out.println( "print Meas Matrix");
+ KkPart.plusEquals(Meas) ;
+ if(Debug){System.out.println("KKPart+Meas");
+ KkPart.print(3,3);}
+ Matrix Kk = PkMinus.times(Hk.transpose()).times(KkPart.inverse());
+ if(Debug){
+ System.out.println(" Kk ");
+ Kk.print(6,3);
+ }
+ Matrix PkPlus =(Matrix.identity(6,6).minus(Kk.times(Hk))).times(PkMinus);
+ if(Debug){
+ System.out.println(" PkPlus ");
+ PkPlus.print(6,6);
+ }
+ double []rrpa=new double[3]; double []rrpb = new double[3];
+ if(Debug){
+ System.out.print("!!!cState="+cState[0]+" cState[1]="+cState[1]+ " cState[2]="+cState[2]);
+ System.out.println("cState[3]="+cState[3]+ " cState[4]="+cState[4]+" cState[5]"+cState[5]);
+ if(Debug){System.out.println("Hk Matrix");
+ Hk.print(3,6);}
+ }
+ rrpb=Hk.times(cState);
+ for(int j=0;j<3;j++){
+ if(getStopTkELow()==true)break;
+ if(Debug)
+ System.out.println("!!!!!!!!!!!!!!bef.subtr!!!!!!!!rrpb[j]="+ rrpb[j]+ " mState[j]="+mState[j]+ " j="+j);
+ rrpb[j]=(mState[j]-rrpb[j]); //[3x1]
+ if(Debug)
+ System.out.println("!!!!!!!!!!!!!!!!!!NEW rrpb[j]="+ rrpb[j]+ " mState[j]="+mState[j]+ " j="+j);
+ }
+ double[]gainFactor=new double[6];
+ gainFactor=Kk.times(rrpb); // [6]
+ for(int i=0;i<3;i++){
+ if(getStopTkELow()==true)break;
+ if(Debug)
+ System.out.println("!!!!!!!!!!!!!!!!!!!!!!!!!!!gainFactor[i]="+ gainFactor[i]+ " cState[i]="+cState[i]+ " i="+i);
+ kRp[i] =cState[i]+gainFactor[i]; // scaled factor 0.25 taken out
+ // kRp[i] +=Kk.times(mState[i]-Hk.times(cState[i]));
+ }
+ return (kupDate = true);
+ }
+/** F=Get the dedx for a given material calculated from
+ *the Bethe-Block-(Sternheimer) calculation
+ *@param pmom (double)absolute particle momentum in GeV/c
+ *@param mass (double)particle mass in GeV
+ *@param materialName a String with the material gone through
+ *@return dedx (double)the energy loss in GeV/mm
+ */
+ public double getDeDx(double pmom, double mass, String materialName){
+ if(Math.abs(pmom/mass )< 1.e-3){ return 0.;}
+ double dedx = 1e-10;
+ double [] p={gevTomev*kRp[3],gevTomev*kRp[4],gevTomev*kRp[5]};
+ if(Debug)System.out.println("MeV px="+p[0]+" py="+p[1]+" pz="+p[2]);
+ MaterialManager manager = MaterialManager.instance();
+ try{
+ org.lcsim.material.Material mat = manager.findMaterial(materialName);
+ double oneCm=1.;
+ dedx = MaterialCalculator.computeBetheBloch(mat, p,(1000*mass), kQ, oneCm);// dedx MeV/cm
+ if(Debug)System.out.println("getDeDx !!! dedx(MeV/cm)="+dedx);
+ dedx*=mevTogev/cmTomm;
+ if(Debug)System.out.println("getDeDx !!! dedx(GeV/mm)="+dedx);
+ } catch (MaterialNotFoundException e){}
+ if(dedx<=0.) return 1e-10;
+ return dedx;
+ }
+
+ public double getRadLength(String materialName){
+ MaterialManager manager = MaterialManager.instance();
+ double radLength=1.E-14;
+ double density=1.E-13;
+ try{
+ org.lcsim.material.Material mat = manager.findMaterial(materialName);
+ double zeff = mat.getZeff();
+ double aeff = mat.getAeff();
+ density= mat.getDensity();
+ radLength=MaterialCalculator.computeRadiationLengthTsai(aeff,zeff);
+ } catch (MaterialNotFoundException e){}
+ if(radLength<=0.) return 1.e-14;
+ return (radLength/density); // in cm
+ }
+ public void setMatNam(String mtn){MaterialName=mtn;}
+ public String getMatNam(){return MaterialName;}
+ public double getX(){return kRp[0];}
+ public double getY(){return kRp[1];}
+ public double getZ(){return kRp[2];}
+
+ public double [] getParameter(){return kRp ;}
+ public double[] getNewRp(){
+ for(int i=0;i<6;i++) nwRp[i]=kRp[i];
+ System.out.println("in CamKalm:xyz"+ nwRp[0]+" "+nwRp[1]+" "+nwRp[2] );
+ System.out.println("in CamKalm:pxyz"+ nwRp[3]+" "+nwRp[4]+" "+nwRp[5] );
+ nwRp[6]=kQ; nwRp[7]=kMass;
+ return nwRp;
+ }
+ public double [][] getCovariance(){return measCov;}
+ public void resetCovariance(){
+ for(int i=0; i<6; i++){
+ for(int j=0;j<6;j++){
+ kCovariance[i][j] = 0.;
+ }
+ }
+ }
+ public void setBField1Dim(double bZ){ kB=bZ;}
+
+ public void setRungeKutta(){kRungeKutta=true;}
+
+ public double getBField1Dim(){return kB;}
+ public void setSign(double xQ){kQ=xQ;}
+ public double[] getBfield(){
+ double []bField=new double[3];
+ try{
+ bField = det.getFieldMap().getField(kRp);
+ }catch(Exception e) {}
+ return bField;
+ }
+
+/**@return vectorial product p x n_unit(rndm:-1:1) -CM- for mult.scat.-March05*/
+ private double [] pCrossN()
+ {
+ double []alfa = new double[3];
+ double [] nnx = new double[3];
+ double alfNorm;
+ /** new direction of vector momentum for ms-C.M.-March 05*/
+ nnx[0]=-1.+2*rndcm.nextDouble();
+ nnx[1]=-1.+2*rndcm.nextDouble();
+ nnx[2]=-1.+2*rndcm.nextDouble();
+ alfa[0] = kRp[4]*nnx[2]-kRp[5]*nnx[1];
+ alfa[1] = kRp[5]*nnx[0]-kRp[3]*nnx[2];
+ alfa[2] = kRp[3]*nnx[1]-kRp[4]*nnx[0];
+ alfNorm = Math.sqrt(alfa[0]*alfa[0]+alfa[1]*alfa[1]+alfa[2]*alfa[2]);
+ alfa[0] = alfa[0]/alfNorm; // without units
+ alfa[1] = alfa[1]/alfNorm;
+ alfa[2] = alfa[2]/alfNorm;
+ if(Debug)System.out.println(" setpxn="+alfa[0]+alfa[1]+alfa[2]);
+ return alfa;
+ }
+ /**@param numRadL(double) the number of radiation Length
+ * @return Multiple Scattering contributions to Momentum componants-C.M.-March05*/
+ private double [] vwk(double numRadL)
+ {
+ double [] alf = new double[3];
+ double []thetak = new double[3];
+ double [] wk = new double[3];
+ double theta0;
+ double pabs = Math.sqrt(kRp[3]*kRp[3]+kRp[4]*kRp[4]+kRp[5]*kRp[5]);
+ if(Debug)System.out.println("numRadL ="+numRadL);
+ if(numRadL==0.) theta0=0.;
+ else{theta0 = 13.6*Math.sqrt(numRadL)*(1.+0.038*Math.log(numRadL));}
+ if(Debug)System.out.println("theta0 first="+theta0);
+ if(pabs==0.)theta0=0.;
+ else{theta0=theta0/(1000.*pabs);}
+ if(Debug)System.out.println("Gauss Random"+rndcm.nextGaussian());
+ alf = pCrossN();
+ if(Debug)System.out.println("alfa="+alf[0]+" "+alf[1]+" "+alf[2]);
+ thetak[0]=rndcm.nextGaussian()*theta0*alf[0];
+ thetak[1]=rndcm.nextGaussian()*theta0*alf[1];
+ thetak[2]=rndcm.nextGaussian()*theta0*alf[2];
+ if(Debug)System.out.println("thetak = "+thetak[0]+" "+thetak[1]+" "+thetak[2]);
+ wk[0]=2*pabs*thetak[0];
+ wk[1]=2*pabs*thetak[1];
+ wk[2]=2*pabs*thetak[2];
+ if(Debug)System.out.println("wk ="+wk[0]+" "+wk[1]+" "+wk[2]);
+ return wk;
+ }
+ /** Set Boolean Indicator If Energy Left
+ */
+ public void setStopTkELow(boolean b)
+ { stopTkELow=b;}
+ /** Reset Boolean Indicator Energy Left for the step
+ */
+ public void resetStopTkELow()
+ {setStopTkELow(false);}
+ /** Return Boolean Indicator Energy Left for the step
+ */
+ public boolean getStopTkELow()
+ {return stopTkELow;}
+}