ExampleProject/src
diff -N TestCamKalmFe.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TestCamKalmFe.java 14 Nov 2008 23:32:55 -0000 1.1
@@ -0,0 +1,265 @@
+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.util.swim.Helix;
+//import org.lcsim.event.LCIOParameters;
+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 12 will Test the CamKalmanStepper program
+ * described in 2006_JINST_1_p1003.pdf.
+ * http://www.iop.org/EJ/abstract/1748-0221/1/10/P10003
+ * and In arXiv: physics/0604197; and physics/0605015
+ * ntrylim represent the number of tracks
+ */
+//
+//
+public class TestCamKalmFe{
+ private static double [] sigma = {1.e-3,1.e-3,1.e-3,1.e-3,1.e-3,1.e-3};
+ private static int _layerN,_nStepsDone;
+ private static double ncount = 0;
+ private static double ntry;
+ private static HelixSwimmer hswim;
+ private static CamKalmanStepper2 CamStpr;
+ private static int _nLayersLim =20;// hadron calorimeter
+ private static double ds = 50.; //path length 50cm.
+ private static int _nStepsToDo = 6; //(px,py) fixed for 5mm in BField=(30mm/5)=6steps
+ private static int _nTracksLim =100;
+ 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 xVx , yVx , zVx ;
+ private static double pxVx , pyVx , pzVx;
+ private static double B = 5.;
+ private static double sign = 1.;
+ private static double kQ = 1.;
+ public static void main(String args[]){
+ System.out.println (" Num. of layers="+_nLayersLim+" ds="+ds+" Num of Tracks="+_nTracksLim);
+ double track [] = new double [6];
+ track[0] = 1.5+0.15; track[1] = 0.5+0.3; track[2] = 0.02+0.2;
+ // track[3],track [4],track[5] below are the Momentum componants
+ // Remark: Here dedx is in GeV/mm, P>dedx otherwise the particle
+ // has not enough energy to move in the material
+ track[3] =3.0 ; track[4] =0.5 ; track[5] = 0.01 ;
+ double muonMass = 0.107;
+ rp0[6] =1.; rp0[7] = muonMass;
+ for (int i=0; i<6; i++){ double xytmp=track[i]; rp0[i]= xytmp;}
+ TestCamKalmFe(rp0,_nTracksLim, _nLayersLim);
+ }
+ // TESTCAMKALMFE
+ //
+ public static void TestCamKalmFe(double [] vect8,int nTracksLim,int nLayersLim){
+ boolean smearedGaussHelix = true;
+ boolean smearedVertex =true;
+ boolean randnum = true;
+ int klb=1; // nStepsToLayback
+ int q;
+ double [] tracksm = new double [8];
+ double[]x = new double[10000];
+ double[]y = new double[10000];
+ double[]z = new double[10000];
+ int [] subPlaneIndex = new int[nLayersLim];
+ double ds2;
+ double k1B2C = 0.299792458*1E-03;
+ double xx, yy, zz;
+ String matterName = "Iron";
+ Random rndm = new Random();
+ for(int i=0;i<6;i++) cov[i][i]=50;
+ pxVx = vect8[3] ; pyVx = vect8[4] ; pzVx = vect8[5];
+ xVx = vect8[0] ; yVx = vect8[1] ; zVx = vect8[2];
+ tracksm[0] = xVx;tracksm[1]=yVx;tracksm[2]=zVx;
+ tracksm[3] = pxVx; tracksm[4]=pyVx;tracksm[5]=pzVx;
+ tracksm[6] = kQ; tracksm[7]=0.105;
+ CamStpr = new CamKalmanStepper2(tracksm, cov);
+ SpaceVector p = new CartesianVector(pxVx,pyVx,pzVx);
+ SpacePoint r0 = new CartesianPoint(xVx, yVx, zVx);
+ q =(int)kQ;
+ // measurement point created by an Helix swimmer
+ hswim=new HelixSwimmer(B);
+ // 1st Source point on track-before smearing
+ System.out.println(" ((((Starting Point Source before smearing))))");
+ System.out.println(" BBBB xVx="+xVx+" yVx="+ yVx+" zVx="+zVx);
+ System.out.println(" BBBB pxVx="+pxVx+" pyVx="+pyVx+" pzVx="+pzVx);
+//
+//==========================================================================
+//
+// create _nTracksLim tracks out of it by smearing
+
+ for(int ntry=0; ntry<_nTracksLim; ntry++){
+ // at start of each track (number ntry)
+ CamStpr.resetStopTkELow(); CamStpr.resetCovariance(50.);
+ randnum=false;
+ // SMEAR VERTEX- NOT NEEDED IF THERE ARE DATA
+ if(randnum == true){ // randum smearing
+ if(smearedVertex){
+ xVx = xVx+(2*rndm.nextDouble()-1.)*sigma[0];
+ yVx = yVx+(2*rndm.nextDouble()-1.)*sigma[1];
+ zVx = zVx+(2*rndm.nextDouble()-1.)*sigma[2];
+ pxVx = pxVx+(2*rndm.nextDouble()-1.)*sigma[3];
+ pyVx = pyVx+(2*rndm.nextDouble()-1.)*sigma[4];
+ pzVx = pzVx+(2*rndm.nextDouble()-1.)*0.0001;
+ if(tracksm[0]<0.)System.out.println
+ ( "====== smaller than 0. xVx="+xVx+"yVx="+yVx+"tracksm[0]="+tracksm[0]);
+ }// End smeared Vx
+ }else{ // Gaussian smearing
+ if(smearedVertex){
+ xVx = xVx + rndm.nextGaussian()*Math.abs(sigma[0]);
+ yVx = yVx + rndm.nextGaussian()*Math.abs(sigma[1]);
+ zVx = zVx + rndm.nextGaussian()*Math.abs(sigma[2]);
+ pxVx = pxVx+ rndm.nextGaussian()*Math.abs(sigma[3]);
+ pyVx = pyVx+ rndm.nextGaussian()*Math.abs(sigma[4]);
+ pzVx = pzVx+ rndm.nextGaussian()*Math.abs(sigma[5]);
+ } // End of Smeared Vertex (Gaussian)
+ }// End smeared Vx
+ tracksm[0] = xVx ; tracksm[1] = yVx ; tracksm[2] = zVx ;
+ tracksm[3] = pxVx; tracksm[4] = pyVx; tracksm[5] = pzVx;
+ aida.cloud2D(" vertex Smeared x versus y").fill(tracksm[0],tracksm[1]);
+ aida.cloud2D(" vertex Smeared x versus z").fill(tracksm[0],tracksm[2]);
+ tracksm[6] = 1.; tracksm[7] = 0.107;
+// TRACK NEW INSTANCE
+ CamStpr.setRpCov(tracksm,cov);
+ CamStpr.setSign(kQ);
+ //***update with measured points
+ tracksm = CamStpr.getNewRp();
+ //*** measured covariance next track
+ for(int kk = 0;kk<3;kk++){
+ for(int k = 0;k<3;k++){
+ if(kk != k) mcov[kk][k]=0.;
+ if(kk == k){ double covs=sigma[kk]*sigma[k]; mcov[kk][k]=covs;}
+ }
+ }
+// DATA: USE THE SMEARED HELIX FROM THE SAME RANDOMIZED VERTEX
+// NOT NEEDED IF THERE ARE REAL DATA
+ p = new CartesianVector(pxVx,pyVx,pzVx);
+ r0 = new CartesianPoint(xVx, yVx, zVx);
+ q = (int)kQ;
+ //*** set hswim for the new randomized atarting points
+ hswim.setTrack( p, r0, q);
+ //*** measurement point role played by an Helix swimmer
+ // Helix- Create measured points on track
+ for(int k = 0;k<_nLayersLim;k++){
+ //spacePoints
+ r0 = hswim.getPointAtLength((k+1)*ds);
+ //*** Layer0 to LayerN on track randomize the data points
+ xx = r0.x();yy = r0.y();zz = r0.z();
+ p = hswim.getMomentumAtLength((k+1)*ds);
+ x[k] = xx ; y[k] = yy ; z[k] = zz;
+ if(smearedGaussHelix){
+ double xxxx = rndm.nextGaussian()*Math.abs(sigma[0]);
+ double yyyy = rndm.nextGaussian()*Math.abs(sigma[1]);
+ double zzzz = rndm.nextGaussian()*Math.abs(sigma[2]);
+ x[k] += xxxx;
+ y[k] += yyyy;
+ z[k] += zzzz;
+ }// End Smearing Helix
+ 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 _nLayersLim
+ //
+// FORWARD
+//
+//*** System.out.println(" AAAAAA START Propagate + UPDATE AAAA for tracks "+ntry);
+ int [ ][ ] numSteps = new int[_nLayersLim][2];
+ int inStepsToDo =10;
+
+ for( int i=0; i<_nLayersLim; i++){//propagate to next layer
+ double pNorm = Math.sqrt(tracksm[3]*tracksm[3]+tracksm[4]*tracksm[4]+tracksm[5]*tracksm[5]);
+ inStepsToDo = (pNorm>1.)?10:50;
+// TRANSPORT
+ numSteps[i] = CamStpr.propagateInBarrelByDs(ds,B,inStepsToDo,matterName);
+ tracksm = CamStpr.getNewRp();
+ _layerN = i; // Layer number 0,1,..5
+ if(CamStpr.getStopTkELow()==true)break;
+// UPDATING
+ double mes0 = x[i], mes1 = y[i], mes2 = z[i];
+ mes[0] = mes0; mes[1] = mes1; mes[2]= mes2;
+ CamStpr.upDate(mes, mcov);
+ tracksm = CamStpr.getNewRp();
+ cov = CamStpr.getCovUpdate();
+ aida.cloud2D("CamStpr smeared tracks update x vs. y").fill(tracksm[0], tracksm[1]);
+ aida.cloud2D("CamStpr smeared tracks update x vs. z").fill(tracksm[0],tracksm[2]);
+ }// End of _nLayersLim
+//
+// BACWARD
+//
+ // update before starting to go backward (using the mesurements of the last point
+ //
+ // reset covariance matrices ====
+ CamStpr.resetCovariance(50.);
+ for(int k=0;k<3;k++){ mcov[k][k]=sigma[k]*sigma[k];}
+ CamStpr.setRpCov(tracksm, cov); // keep position change cov mtx
+ double dsBack =-ds; // (.mm)
+ if(CamStpr.getStopTkELow()== false) {
+ CamStpr.upDate(mes, mcov);
+ tracksm = CamStpr.getNewRp();
+ aida.cloud2D("CamStpr Back smeared x versus y update").fill(tracksm[0], tracksm[1]);
+ aida.cloud2D("CamStpr Back smeared x versus z update ").fill(tracksm[0], tracksm[2]);
+ }
+
+ cov = CamStpr.getCovUpdate();
+//
+// TRANSPORT BACK
+//
+ for( int j=(_layerN-1);j>=0; j--){ //propagate to x[i], y[i], z[i]
+ if(j==_layerN-1&&CamStpr.getStopTkELow() ) CamStpr.setStopTkELow(false);
+ CamStpr.propagateInBarrelBackByDs(dsBack,B,matterName,numSteps[j+1]);
+ tracksm = CamStpr.getNewRp();
+ xx = x[j]; yy = y[j]; zz=z[j] ;
+ mes[0] = xx ; mes[1] = yy; mes[2]= zz ;
+// UPDATE BACK
+ CamStpr.upDate(mes, mcov);
+ tracksm = CamStpr.getNewRp(); // get x,y,z,px,py,pz,q,mass
+ cov = CamStpr.getCovUpdate();
+ aida.cloud2D("CamStpr Back smeared x versus y update").fill(tracksm[0], tracksm[1]);
+ aida.cloud2D("CamStpr Back smeared x versus z update ").fill(tracksm[0], tracksm[2]);
+ aida.cloud1D("CamStprBack smeared X ( measure-updated) ").fill((mes[0]-tracksm[0]));
+ aida.cloud1D("CamStprBack smeared Y ( measure-updated) ").fill((mes[1]-tracksm[1]));
+ aida.cloud1D("CamStprBack smeared Z ( measure-updated) ").fill((mes[2]-tracksm[2]));
+ }
+ //*** System.out.println("AAAAAAA Back to the Vertex AAAAAAA ")
+ int [ ] nStepsToVx = new int[2];
+ nStepsToVx[0] = numSteps[0][0]; nStepsToVx[1] = numSteps[0][1];
+ double [] par = {xVx, yVx,zVx,pxVx,pyVx,pzVx};
+ double pLast = Math.sqrt(tracksm[3]*tracksm[3]+tracksm[4]*tracksm[4]+tracksm[5]*tracksm[5]);
+ //double dsToVx = (xVx-tracksm[0])*(pLast/tracksm[3]);
+ double dsToVx = dsBack;
+ /*
+ System.out.print(" ####dsToVx="+dsToVx+" nStepsDone="+nStepsToVx[1]);
+ System.out.println(" track number="+ntry+"dsBack="+dsBack);
+ System.out.println(" nStepsToVx[0] ="+nStepsToVx[0]+" nStepsToVx[1] ="+nStepsToVx[1]);
+ System.out.println(" dsBack ="+dsBack+" B="+B);
+ */
+ CamStpr.propagateInBarrelBackByDs(dsToVx,B,matterName,nStepsToVx);
+ tracksm = CamStpr.getNewRp();
+ cov=CamStpr.getCovUpdate();
+ if((ntry<=2) ||(ntry>(_nTracksLim-2))){
+ System.out.println("FINAL ### tracksm[0] = "+tracksm[0]+" par[0]="+par[0]);
+ System.out.println(" tracksm[1] = "+tracksm[1]+" par[1]="+par[1]);
+ System.out.println(" tracksm[2] = "+tracksm[2]+" par[2]="+par[2]);
+ System.out.println(" tracksm[3] = "+tracksm[3]+" par[3]="+par[3]);
+ System.out.println(" xVx="+xVx+"yVx="+yVx+" zVx="+zVx);
+ System.out.println(" pxVx="+pxVx+"pyVx="+pyVx+" pzVx="+pzVx);
+ }
+ if((ntry<=2) ||(ntry>(_nTracksLim-2))){
+ for(int n=0; n<3;n++)
+ System.out.println( "####n="+n+" par[n]- tracksm[n]="+(par[n]-tracksm[n]));
+ }
+ aida.cloud1D("delta_x ").fill((par[0]-tracksm[0]));
+ aida.cloud1D("delta_y").fill(par[1]-tracksm[1]);
+ aida.cloud1D("delta_z ").fill((par[2]-tracksm[2]));
+ aida.cloud1D("delta_pz ").fill((par[5]-tracksm[5]));
+ double ptBackVx = Math.sqrt(tracksm[3]*tracksm[3]+tracksm[4]*tracksm[4]);
+ double ptVx = Math.sqrt(par[3]*par[3]+par[4]*par[4]);
+ aida.cloud1D("delta_pt").fill(ptVx-ptBackVx);
+ } // End of each track- Each ntry
+ }// End of test8 method
+ }// End of the java/class software
\ No newline at end of file