3 removed files
lcsim/src/org/lcsim/contrib/FermiLab/Cam/CamKalman
diff -N CamKalmanStepper.java
--- CamKalmanStepper.java 23 Aug 2007 21:09:45 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,460 +0,0 @@
-package org.lcsim.contrib.FermiLab.Cam.CamKalman;
-import java.util.Random;
-
-import org.lcsim.geometry.Detector;
-import org.lcsim.material.MaterialCalculator;
-import org.lcsim.material.MaterialManager;
-import org.lcsim.material.MaterialNotFoundException;
-import org.lcsim.recon.cluster.util.CalHitMapMgr;
-import org.lcsim.util.aida.AIDA;
-
-import Jama.Matrix;
- /**
- * CamKalmanStepper- A Kalman Stepper
- * trasnsported from JAS2
- *JINST - published October 2006
- *arXiv: physics/0604197
- *@Author-C. Milstene <[log in to unmask]>
- * Code developped consulting Hans Wenzel
- * More comments and references to come
- * Hard-coded material Silicon for the dE/dx
- * (Testing program-TestKalm.java contains
- * a main- and is the (Driver/MainProject))
- */
-public class CamKalmanStepper{
- 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 CamKalmanStepper(){
- 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 CamKalmanStepper(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 CamKalmanStepper(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);
- }
- }
- }
- /** 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");
- return stopTkELow=true ;
- }
- double rdL= getRadLength("Iron"); //im mm
- if(Debug)
- System.out.println(" number 1: dpdx="+ dpdx);
- wk=vwk((getRadLength("Iron"))*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],0.,0.},{0.,0.,0.,0.,wk[1],0.},{0.,0.,0.,0.,0.,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]);
- k*=ds;
- 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;
- 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);
- fi=fi.times(ds);
- 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)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);
- 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("!!!!!!!!!!!!!!!!!!!!!!!!!!!!rrpb[j]="+ rrpb[j]+ " mState[j]="+mState[j]+ " j="+j);
- rrpb[j]=(mState[j]-rrpb[j]); //[3x1]
- }
- 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];
- 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[4]*nnx[0]-kRp[3]*nnx[1];
- 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;}
-}
lcsim/src/org/lcsim/contrib/FermiLab/Cam/CamKalman
diff -N Test1.java
--- Test1.java 29 Nov 2007 02:25:57 -0000 1.4
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,260 +0,0 @@
-package org.lcsim.contrib.FermiLab.Cam.CamKalman;
-import java.util.Random;
-
-import org.lcsim.spacegeom.CartesianPoint;
-import org.lcsim.spacegeom.CartesianVector;
-import org.lcsim.spacegeom.SpacePoint;
-import org.lcsim.spacegeom.SpaceVector;
-import org.lcsim.util.aida.AIDA;
-import org.lcsim.util.swim.HelixSwimmer;
-
-/*new modifications as of 06/29/07
- * by Caroline Milstene and Kudzanayi Munetsi-Mugomba
- * Test 1 will Test the CamKalmanStepper program
- * xe, ye, ze,
- * ntrylim represent the number of tracks
- *
-*/
-
-public class Test1{
-
- private static double [] sigma={5e-1,5e-1,5e-1,1e-4,1e-4,1e-4};
- private static int npoints;
- private static double ncount=0;
- private static double ntry;
- private static double [] tracksm = new double [8];
- private static int npointslim=10;// hadron calorimeter
- // private static int npointslim=5; // Silicon Tracker
- private static int klim=4;
- private static int ntrylim =10;
- private static double[][] cov = new double[6][6];
- private static double[][] mcov = new double[3][3];
- protected static double [] rp0=new double[8];
- public static AIDA aida = AIDA.defaultInstance();
- private static double []mes = new double[3];
- private static double B=5;
- private static double sign=1.;
- private static double kQ=+1;
-
- public static void main(String args[]){
- double track [] = new double [6];
- track[0]=1.5; track[1]=0.5; track[2]=0.02; track[3]=10.; track[4]=0.; track[5]=0.1;
- double mMass=0.106;
- rp0[6]=1.; rp0[7]=mMass;
- for (int i=0; i<6; i++)
- {
- rp0[i]= track[i];
- }
- Test1(rp0,ntrylim, npointslim);
- }
-
- public static void Test1(double [] rppp,int ntrylim,int npointslim){
-
- double k1B2C= 0.299792458*1E-03;
- double xx, yy, zz;
- double[] rpp1=new double[6];
- boolean randnum=true;
- double []se= new double[1000];
- double []xe = new double[1000]; double []ye = new double[1000]; double []ze = new double[1000];
- double[]x=new double[1000]; double[]y=new double[1000]; double[]z=new double[1000];
-
- double[]x_tr=new double[1000]; double[]y_tr=new double[1000]; double[]z_tr=new double[1000];
- double[]px_tr=new double[1000];double[]py_tr=new double[1000];double[]pz_tr=new double[1000];
- for(int i=0;i<6;i++) cov[i][i]=50;
- double pxhlx=rppp[3], pyhlx=rppp[4],pzhlx=rppp[5];
- double xhlx=rppp[0],yhlx=rppp[1],zhlx=rppp[2];
- SpaceVector p = new CartesianVector(pxhlx,pyhlx,pzhlx);
- SpacePoint r0 = new CartesianPoint(xhlx, yhlx, zhlx);
- // measurement point created by an Helix swimmer
- HelixSwimmer hswim=new HelixSwimmer(B);
- int q=(int)kQ;
- hswim.setTrack( p, r0, q);
- se[0]=Math.sqrt(pxhlx*pxhlx+pyhlx*pyhlx+pzhlx*pzhlx);
- r0= hswim.getPointAtLength(se[0]);
- // 1st point on track
- xe[0]=r0.x();ye[0]=r0.y();ze[0]=r0.z();
- // measurement point created by the Kalman stepper
- x_tr[0]=rppp[0]; y_tr[0]=rppp[1]; z_tr[0]=rppp[2];
- px_tr[0]=rppp[3]; py_tr[0]=rppp[4]; pz_tr[0]=rppp[5];
- CamKalmanStepper kalmStpr = new CamKalmanStepper(rppp, cov);
- kalmStpr.setSign(sign);
- aida.cloud2D("Helix x versus y").fill(r0.x(), r0.y());
- aida.cloud2D("Helix x versus z").fill(r0.x(), r0.z());
- aida.cloud2D("Kalman Transport x versus y").fill(x_tr[0], y_tr[0]);
- aida.cloud2D("Kalman Transport x versus z").fill(x_tr[0], z_tr[0]);
-
- aida.cloud2D("Kalman Transport all the steps x versus y").fill(x_tr[0], y_tr[1]);
- aida.cloud2D("Kalman Transport all the steps x versus z").fill(x_tr[0], z_tr[2]);
- System.out.println(" Kalman/Helix 1st point :x="+rppp[0]+" y="+rppp[1]+" z="+rppp[2]);
-
- aida.cloud2D("Px vs. Py - Kalman -Fe").fill(rppp[3], rppp[4]);
- aida.cloud2D("tracksm[3] vs. tracksm[4]").fill(tracksm[3], tracksm[4]);
-
- double ds=20.;
-//
-// create measured points on track
- for(int j=1; j<npointslim; j++){
- // swimmer
- se[j]=se[(j-1)]+ds;
- //spacePoints
- r0= hswim.getPointAtLength(se[j]);
- // 2nd to npoints on track
- xe[j]=r0.x();ye[j]=r0.y();ze[j]=r0.z();
- aida.cloud2D("Helix x versus y").fill(r0.x(), r0.y());
- aida.cloud2D("Helix x versus z").fill(r0.x(), r0.z());
- }
- // Kalman transport
- ds=20.;
- for(int j=1; j<npointslim; j++){
- //stepper -make klim steps to cover ds=100mm
- for(int k=1;k<klim;k++){ // 10 steps
- double ds2=ds/klim;
- kalmStpr.stepByDs(ds2);
- rpp1=kalmStpr.getParameter(); // get x,y,z,px,py,pz
- x_tr[k]=rpp1[0];y_tr[k]=rpp1[1];z_tr[k]=rpp1[2];
- px_tr[k]=rpp1[3]; py_tr[k]=rpp1[4]; pz_tr[k]=rpp1[5];
- if(kalmStpr.getStopTkELow()){
- break;
- }
- aida.cloud2D("Kalman Transport all the steps x versus y").fill(rpp1[0], rpp1[1]);
- aida.cloud2D("Kalman Transport all the steps x versus z").fill(rpp1[0], rpp1[2]);
-
- aida.cloud2D("Kalman Transport all the steps px versus py").fill(rpp1[3], rpp1[4]);
- aida.cloud2D("Kalman Transport all the steps px versus pz").fill(rpp1[3], rpp1[5]);
- }
- if(kalmStpr.getStopTkELow())break;
- aida.cloud2D("Kalman Transport x versus y").fill(rpp1[0], rpp1[1]);
- aida.cloud2D("Kalman Transport x versus z").fill(rpp1[0], rpp1[2]);
- }
- // Kalman transport + update
- // re-do Kalman level 1 before updating starting at point 0
- // with starting momenta- then update
- for(int i=0;i<6;i++) cov[i][i]=50;
- double [ ]rpp2= new double[6];
- double [ ]rpp3=new double[8];
- rpp3[0]=x_tr[0]; rpp3[1]=y_tr[0]; rpp3[2]=z_tr[0];
- rpp3[3]=px_tr[0]; rpp3[4]=py_tr[0];rpp3[5]=pz_tr[0];
- rpp3[6]=1.;rpp3[7]=0.105;
- kalmStpr.setSign(1);
- kalmStpr = new CamKalmanStepper(rpp3, cov);
- // CamKalmanStepper kalmStpr2 = new CamKalmanStepper(rpp3, cov);
- kalmStpr.resetStopTkELow();
- //update with measured points
- System.out.println("PRINT updating");
- for(int k=0;k<3;k++) mcov[k][k]=1.e-03;
- mes[0]= xe[0]; mes[1]= ye[0]; mes[2]= ze[0];
- kalmStpr.upDate(mes, mcov);
- rpp2=kalmStpr.getParameter(); // get x,y,z,px,py,pz
- aida.cloud2D("kalmStpr after update all the steps x vs. y").fill(rpp2[0], rpp2[1]);
- aida.cloud2D("kalmStpr after update all the steps x vs. z").fill(rpp2[0], rpp2[2]);
- aida.cloud2D("Kalman x versus y after update").fill(rpp2[0], rpp2[1]);
- aida.cloud2D("Kalman x versus after update z").fill(rpp2[0], rpp2[2]);
-
-//
- for(int j=1; j<npointslim; j++){
- //stepper -make 4 steps to cover ds=20mm
- for(int k=1;k<klim;k++){ // 4 steps
- double ds2=ds/klim;
- kalmStpr.stepByDs(ds2);
- rpp2=kalmStpr.getParameter();
- if(kalmStpr.getStopTkELow())break;
- aida.cloud2D("kalmStpr after update all the steps x vs. y").fill(rpp2[0], rpp2[1]);
- aida.cloud2D("kalmStpr after update all the steps x vs. z").fill(rpp2[0], rpp2[2]);
-//
- }
- if(kalmStpr.getStopTkELow())break;
- //update with measured points
- System.out.println("PRINT updating");
- mes[0]= xe[j]; mes[1]= ye[j]; mes[2]= ze[j];
- kalmStpr.upDate(mes, mcov);
- rpp2=kalmStpr.getParameter(); // get x,y,z,px,py,pz
- aida.cloud2D("Kalman x versus y after update").fill(rpp2[0], rpp2[1]);
- aida.cloud2D("Kalman x versus after update z").fill(rpp2[0], rpp2[2]);
- }
-//========================================================================================================
- // create ntrylim tracks out of it by smearing
-
- if(ntrylim==-1)return;
- for(int ntry=0; ntry<ntrylim; ntry++){
- kalmStpr.resetStopTkELow(); // at start of track
- npoints=0; // reset the number of points for each track
- Random rndm =new Random();
- randnum=true;
- // smear the IP-(xyz,px,py,pz) within the precision
- if(randnum == true){
- tracksm[0]=x_tr[0]+(2*rndm.nextDouble()-1.)*0.0001; tracksm[1]=y_tr[0]+(2*rndm.nextDouble()-1.)*0.0001;
- tracksm[2]=z_tr[0]+(2*rndm.nextDouble()-1.)*0.0001;
- tracksm[3]=px_tr[0]+(2*rndm.nextDouble()-1.)*0.0001; tracksm[4]=py_tr[0]+(2*rndm.nextDouble()-1.)*0.0001;
- tracksm[5]=pz_tr[0]+(2*rndm.nextDouble()-1.)*0.0001;
- }else{
- tracksm[0]=x_tr[0]+rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[1])); tracksm[1]=y_tr[0]+rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[1]));
- tracksm[2]=z_tr[0]+rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[2]));
- tracksm[3]=px_tr[0]+rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[3])); tracksm[4]=py_tr[0]+rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[4]));
- tracksm[5]=pz_tr[0]+rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[5]));
- }
- tracksm[6]=1.; tracksm[7]=0.106;
-
- // measured covariance next track
- for(int j=0;j<3;j++) mcov[j][j]=1e-3;
- for(int j=0; j<6;j++) cov[j][j]=50.;
- // smeared Helix
- for(int k=0;k<npointslim;k++){
- xx=xe[k]; yy=ye[k]; zz=ze[k];
- x[k]=xx; y[k]=yy; z[k]=zz;
- x[k]+=rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[0]));
- y[k]+=rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[1]));
- z[k]+=rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[2]));
- aida.cloud2D("Smeared Helix x versus y").fill(x[k], y[k]);
- aida.cloud2D("Smeared Helix x versus z").fill(x[k], z[k]);
- }
- // Transport Matrix
- kalmStpr = new CamKalmanStepper(tracksm, cov);
- kalmStpr.setSign(sign);
- kalmStpr.resetCovariance();
- for(int i=0;i<6;i++) cov[i][i]=50;
- ds=20.;
- System.out.println("npointslim="+ npointslim);
- for( int i=0; i<npointslim; i++){ //propagate to xe[i], ye[i], ze[i]
- // if (ncount=0)
- npoints = i;
- for(int k=0;k<klim;k++){
- if(kalmStpr.getStopTkELow()==true) break;
- double ds2=ds/klim;
- kalmStpr.stepByDs(ds2);
- double [] ttt=new double[6];
- ttt=kalmStpr.getParameter();
-
- for(int m =0;m<6;m++){
- tracksm[m]=ttt[m];
- //+rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[m]));
- tracksm[6]=1.; tracksm[7]=0.106;
-
- aida.cloud2D(" Smeared Kalman all steps x versus y").fill(tracksm[0], tracksm[1]);
- aida.cloud2D(" Smeared Kalman all steps x versus z").fill(tracksm[0], tracksm[2]);
- }
- aida.cloud2D(" Smeared Kalman x versus y").fill(tracksm[0], tracksm[1]);
- aida.cloud2D(" Smeared Kalman x versus z").fill(tracksm[0], tracksm[2]);
-
- if(kalmStpr.getStopTkELow()==true) break;
- //update with measured points
- kalmStpr.resetStopTkELow(); // at start of track
- kalmStpr.resetCovariance();
- mes[0]= x[i]; mes[1]= y[i]; mes[2]= z[i];
- kalmStpr.upDate(mes, mcov);
- ttt=kalmStpr.getParameter();
- for(int m =0;m<6;m++){
- tracksm[m]=ttt[m];
- tracksm[6]=1.; tracksm[7]=0.106;
-
- aida.cloud2D("kalmStpr smeared tracks after update x vs. y").fill(tracksm[0], tracksm[1]);
- aida.cloud2D("kalmStpr smeared tracks after update x vs. z").fill(tracksm[0],tracksm[2]);
- }
- }
- if(kalmStpr.getStopTkELow()==true)break;}
- kalmStpr.resetStopTkELow(); // at start of track
- kalmStpr.resetCovariance();
- }
- }
- }
-
-
lcsim/src/org/lcsim/contrib/FermiLab/Cam/CamKalman
diff -N Test2.java
--- Test2.java 16 May 2008 05:09:36 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,276 +0,0 @@
-package org.lcsim.contrib.FermiLab.Cam.CamKalman;
-
-import java.util.*;
-import org.lcsim.util.aida.*;
-import Jama.*;
-import Jama.util.*;
-import org.lcsim.geometry.Detector;
-import org.lcsim.util.step.*;
-import org.lcsim.util.swim.HelixSwimmer;
-import org.lcsim.spacegeom.CartesianPoint;
-import org.lcsim.spacegeom.CartesianVector;
-import org.lcsim.spacegeom.SpaceVector;
-import org.lcsim.spacegeom.SpacePoint;
-
-/*new modifications as of 06/29/07
- * by Caroline Milstene and Kudzanayi Munetsi-Mugomba
- * Test 1 will Test the CamKalmanStepper program
- * xe, ye, ze,
- * ntrylim represent the number of tracks
- *
-*/
-
-public class Test2{
-
- private static double [] sigma={5e-1,5e-1,5e-1,1e-4,1e-4,1e-4};
- private static int npoints;
- private static double ncount=0;
- private static double ntry;
- private static double [] tracksm = new double [8];
- private static int npointslim=60;// hadron calorimeter
- // private static int npointslim=5; // Silicon Tracker
- private static int klim=4;
- private static int ntrylim =20;
- private static double[][] cov = new double[6][6];
- private static double[][] mcov = new double[3][3];
- protected static double [] rp0=new double[8];
- public static AIDA aida = AIDA.defaultInstance();
- private static double []mes = new double[3];
- private static double B=5;
- private static double sign=1.;
- private static double kQ=+1;
-
- public static void main(String args[]){
- double track [] = new double [6];
- track[0]=1.5; track[1]=0.5; track[2]=0.02; track[3]=0.5; track[4]=0.; track[5]=0.1;
- double mMass=0.106;
- rp0[6]=1.; rp0[7]=mMass;
- for (int i=0; i<6; i++)
- {
- rp0[i]= track[i];
- }
- Test2(rp0,ntrylim, npointslim);
- }
-
- public static void Test2(double [] vect8,int ntrylim,int npointslim){
-
- double k1B2C= 0.299792458*1E-03;
- double xx, yy, zz;
- double[] vectSix=new double[6];
- boolean randnum=true;
- double []se= new double[1000];
- double []xe = new double[1000]; double []ye = new double[1000]; double []ze = new double[1000];
- double[]x=new double[1000]; double[]y=new double[1000]; double[]z=new double[1000];
-
- double[]x_tr=new double[1000]; double[]y_tr=new double[1000]; double[]z_tr=new double[1000];
- double[]px_tr=new double[1000];double[]py_tr=new double[1000];double[]pz_tr=new double[1000];
- for(int i=0;i<6;i++) cov[i][i]=50;
- double pxhlx=vect8[3], pyhlx=vect8[4],pzhlx=vect8[5];
- double xhlx=vect8[0],yhlx=vect8[1],zhlx=vect8[2];
- double pxKalm=vect8[3], pyKalm=vect8[4],pzKalm=vect8[5];
- double xKalm=vect8[0],yKalm=vect8[1],zKalm=vect8[2];
- SpaceVector p = new CartesianVector(pxhlx,pyhlx,pzhlx);
- SpacePoint r0 = new CartesianPoint(xhlx, yhlx, zhlx);
- // measurement point created by an Helix swimmer
- HelixSwimmer hswim=new HelixSwimmer(B);
- int q=(int)kQ;
- hswim.setTrack( p, r0, q);
- se[0]=Math.sqrt(pxhlx*pxhlx+pyhlx*pyhlx+pzhlx*pzhlx);
- r0= hswim.getPointAtLength(se[0]);
- // 1st point on track
- xe[0]=r0.x();ye[0]=r0.y();ze[0]=r0.z();
- // measurement point created by the Kalman stepper
- x_tr[0]=xKalm; y_tr[0]=yKalm; z_tr[0]=zKalm;
- px_tr[0]=pxKalm; py_tr[0]=pyKalm; pz_tr[0]=pzKalm;
- double []vectEight=new double[8];
- vectEight[0]=xKalm; vectEight[1]=yKalm;vectEight[2]=zKalm;
- vectEight[3]=pxKalm; vectEight[4]=pyKalm;vectEight[5]=pzKalm;
- vectEight[6]=kQ; vectEight[7]=0.105;
- CamKalmanStepper kalmStpr = new CamKalmanStepper(vectEight, cov);
- kalmStpr.setSign(sign);
- aida.cloud2D("Helix x versus y").fill(r0.x(), r0.y());
- aida.cloud2D("Helix x versus z").fill(r0.x(), r0.z());
- aida.cloud2D("Kalman Transport x versus y").fill(x_tr[0], y_tr[0]);
- aida.cloud2D("Kalman Transport x versus z").fill(x_tr[0], z_tr[0]);
-
- aida.cloud2D("Kalman Transport all the steps x versus y").fill(x_tr[0], y_tr[1]);
- aida.cloud2D("Kalman Transport all the steps x versus z").fill(x_tr[0], z_tr[2]);
- System.out.println(" Kalman/Helix 1st point :x="+vect8[0]+" y="+vect8[1]+" z="+vect8[2]);
-
- aida.cloud2D("Px vs. Py - Kalman -Fe").fill(vect8[3], vect8[4]);
- aida.cloud2D("tracksm[3] vs. tracksm[4]").fill(tracksm[3], tracksm[4]);
-
- double ds=20.;
-//
-// Helix- Create measured points on track
-//
- for(int j=1; j<npointslim; j++){
- // swimmer
- se[j]=se[(j-1)]+ds;
- //spacePoints
- r0= hswim.getPointAtLength(se[j]);
- // 2nd to npoints on track
- xe[j]=r0.x();ye[j]=r0.y();ze[j]=r0.z();
- aida.cloud2D("Helix x versus y").fill(r0.x(), r0.y());
- aida.cloud2D("Helix x versus z").fill(r0.x(), r0.z());
- }
- //
- // Kalman transport
- //
- ds=20.;
- for(int j=1; j<npointslim; j++){
- //stepper -make klim steps to cover ds=100mm
- for(int k=1;k<klim;k++){ // 10 steps
- double ds2=ds/klim;
- kalmStpr.stepByDs(ds2);
- vectEight=kalmStpr.getNewRp(); // get x,y,z,px,py,pz, q,mass
- x_tr[k+(j-1)*klim]=vectEight[0];y_tr[k+(j-1)*klim]=vectEight[1];z_tr[k+(j-1)*klim]=vectEight[2];
- px_tr[k+(j-1)*klim]=vectEight[3]; py_tr[k+(j-1)*klim]=vectEight[4]; pz_tr[k+(j-1)*klim]=vectEight[5];
- if(kalmStpr.getStopTkELow()){
- break;
- }
- aida.cloud2D("Kalman Transport all the steps x versus y").fill(vectEight[0], vectEight[1]);
- aida.cloud2D("Kalman Transport all the steps x versus z").fill(vectEight[0], vectEight[2]);
-
- aida.cloud2D("Kalman Transport all the steps px versus py").fill(vectEight[3], vectEight[4]);
- aida.cloud2D("Kalman Transport all the steps px versus pz").fill(vectEight[3], vectEight[5]);
- }
- if(kalmStpr.getStopTkELow())break;
- aida.cloud2D("Kalman Transport x versus y").fill(vectEight[0], vectEight[1]);
- aida.cloud2D("Kalman Transport x versus z").fill(vectEight[0], vectEight[2]);
- }
- // Kalman transport + update
- // re-do Kalman level 1 before updating starting at point 0
- // with starting momenta- then update
- for(int i=0;i<6;i++) cov[i][i]=50;
- vectEight=new double[8];
- // First -IP-
- vectEight[0]=xKalm; vectEight[1]=yKalm; vectEight[2]=zKalm;
- vectEight[3]=pxKalm; vectEight[4]=pyKalm;vectEight[5]=pzKalm;
- vectEight[6]=1.;vectEight[7]=0.105;
- kalmStpr.setSign(1);
- kalmStpr = new CamKalmanStepper(vectEight, cov);
- // CamKalmanStepper kalmStpr2 = new CamKalmanStepper(vectEight, cov);
- kalmStpr.resetStopTkELow();
- //update with measured points
- System.out.println("PRINT updating");
- for(int k=0;k<3;k++) mcov[k][k]=1.e-03;
- mes[0]= xe[0]; mes[1]= ye[0]; mes[2]= ze[0];
- kalmStpr.upDate(mes, mcov);
- vectEight=kalmStpr.getNewRp(); // get x,y,z,px,py,pz,q,mass
- aida.cloud2D("kalmStpr after update all the steps x vs. y").fill(vectEight[0], vectEight[1]);
- aida.cloud2D("kalmStpr after update all the steps x vs. z").fill(vectEight[0], vectEight[2]);
- aida.cloud2D("Kalman x versus y after update").fill(vectEight[0], vectEight[1]);
- aida.cloud2D("Kalman x versus after update z").fill(vectEight[0], vectEight[2]);
-//
- for(int j=1; j<npointslim; j++){
- //stepper -make 4 steps to cover ds=20mm
- for(int k=1;k<klim;k++){ // 4 steps
- double ds2=ds/klim;
- kalmStpr.stepByDs(ds2);
- vectEight=kalmStpr.getNewRp();
- if(kalmStpr.getStopTkELow()) break;
- aida.cloud2D("kalmStpr after update all the steps x vs. y").fill(vectEight[0], vectEight[1]);
- aida.cloud2D("kalmStpr after update all the steps x vs. z").fill(vectEight[0], vectEight[2]);
-//
- }
- if(kalmStpr.getStopTkELow())break;
- //update with measured points
- System.out.println("PRINT updating");
- mes[0]= xe[j]; mes[1]= ye[j]; mes[2]= ze[j];
- kalmStpr.upDate(mes, mcov);
- vectEight=kalmStpr.getNewRp(); // get x,y,z,px,py,pz,q,mass
- aida.cloud2D("Kalman x versus y after update").fill(vectEight[0], vectEight[1]);
- aida.cloud2D("Kalman x versus after update z").fill(vectEight[0], vectEight[2]);
- }
-//========================================================================================================
- // create ntrylim tracks out of it by smearing
-
- if(ntrylim==-1)return;
- for(int ntry=0; ntry<ntrylim; ntry++){
- // at start of each track (number ntry)
- kalmStpr.resetStopTkELow();
- kalmStpr.resetCovariance();
- for(int i=0;i<6;i++) cov[i][i]=50;
- npoints=-1; // reset the number of points for each track
- Random rndm =new Random();
- randnum=true;
- // smear the IP-(xyz,px,py,pz) within the precision
-
- if(randnum == true){
- tracksm[0]=xKalm+Math.abs(2*rndm.nextDouble()-1.)*0.0001;
- tracksm[1]=yKalm+2*(rndm.nextDouble()-1.)*0.0001;
- tracksm[2]=zKalm+(2*rndm.nextDouble()-1.)*0.0001;
- if(tracksm[0]<0.)System.out.println
- ( "==========smaller than 0. xKalm="+xKalm+"yKalm="+yKalm+" tracksm[0]="+tracksm[0]);
- tracksm[3]=pxKalm+(2*rndm.nextDouble()-1.)*0.0001; tracksm[4]=pyKalm+(2*rndm.nextDouble()-1.)*0.0001;
- tracksm[5]=pzKalm+(2*rndm.nextDouble()-1.)*0.0001;
- }else{
- tracksm[0] = xKalm+rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[1]));
- tracksm[1] = yKalm+rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[1]));
- tracksm[2] = zKalm+rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[2]));
- tracksm[3] = pxKalm+rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[3]));
- tracksm[4] = pyKalm+rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[4]));
- tracksm[5] = pzKalm+rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[5]));
- }
- tracksm[6]=1.; tracksm[7]=0.106;
- //
- // measured covariance next track
- //
- for(int j=0;j<3;j++) mcov[j][j]=1e-3;
-
- // smeared Helix
- for(int k=0;k<npointslim;k++){
- xx=xe[k]; yy=ye[k]; zz=ze[k];
- x[k]=xx; y[k]=yy; z[k]=zz;
- x[k]+= 2*rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[0]));
- y[k]+= 2*rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[1]));
- z[k]+= 2*rndm.nextGaussian()*Math.abs(Math.sqrt(sigma[2]));
- aida.cloud2D("Smeared Helix x versus y").fill(x[k], y[k]);
- aida.cloud2D("Smeared Helix x versus z").fill(x[k], z[k]);
- }// end of npointslim
- //
- // Transport Matrix
- //
- if((ntry==8)|| (ntry==9)){
- System.out.print(" tracksm[0]="+tracksm[0]+" tracksm[1]="+tracksm[1]+" tracksm[2]="+tracksm[2]);
- System.out.print(" tracksm[3]="+tracksm[3]+" tracksm[4]="+tracksm[4]+" tracksm[5]="+tracksm[5]);
- System.out.println(" ntry="+ntry);
- }
- // New Instance for that track
- kalmStpr = new CamKalmanStepper(tracksm, cov);
- kalmStpr.setSign(sign);
- ds=20.;
- int kpoints=-1;
- System.out.println("npointslim="+ npointslim);
- for( int i=0; i<npointslim; i++){ //propagate to xe[i], ye[i], ze[i]
- // if (ncount=0)
- npoints = i;
- for(int k=0;k<klim;k++){
- double ds2=ds/klim;
- kalmStpr.stepByDs(ds2);
- tracksm = kalmStpr.getNewRp();
- if(tracksm[0]<0.)System.out.println(" smaller than 0 tracksm[0]="+tracksm[0]+
- " i="+i+" k="+k+" ntry="+ntry);
- aida.cloud2D(" Smeared Kalman all steps x versus y").fill(tracksm[0], tracksm[1]);
- aida.cloud2D(" Smeared Kalman all steps x versus z").fill(tracksm[0], tracksm[2]);
- if(kalmStpr.getStopTkELow()==true) {
- break;
- }
- }
- if(kalmStpr.getStopTkELow()==true) break;
- aida.cloud2D(" Smeared Kalman x versus y").fill(tracksm[0], tracksm[1]);
- aida.cloud2D(" Smeared Kalman x versus z").fill(tracksm[0], tracksm[2]);
-
- //update with measured points
- kalmStpr.resetCovariance();
- mes[0]= x[i]; mes[1]= y[i]; mes[2]= z[i];
- kalmStpr.upDate(mes, mcov);
- tracksm=kalmStpr.getNewRp();
- aida.cloud2D("kalmStpr smeared tracks after update x vs. y").fill(tracksm[0], tracksm[1]);
- aida.cloud2D("kalmStpr smeared tracks after update x vs. z").fill(tracksm[0],tracksm[2]);
- }// End of npointslim
- } // End of each track- Each ntry
- }// End of test2 method
- }// End of the java/class software
-
\ No newline at end of file
CVSspam 0.2.8
