13 added files
lcsim/src/org/lcsim/contrib/SODTracker
diff -N History
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ History 10 Mar 2006 21:21:34 -0000 1.1
@@ -0,0 +1,26 @@
+## History file
+## Silicon Outer Detector (SOD) Tracker
+## Please summarize changes to the SODTracker code here.
+## Most recent first please.
+
+10 March 2006 Fred Blanc V01-00-00
+ - Create SODTracker package. Included files:
+ SODFittedCir.java
+ SODHel.java
+ SODHit.java
+ SODHitAdder.java
+ SODL1Path.java
+ SODNoiseMaker.java
+ SODReadGeom.java
+ SODTrack.java
+ SODTrackFinder.java
+ SODTrackFinderDriver.java
+ SODmatinv.java
+ - SODTrackFinderDriver.java is the main driver
+ - Example code is given in test/TestSOD.java
+ - The default functionality is: 1) find track hits in the
+ vertex barrel from MCParticle list, 2) add hits from the outer
+ detector are added, 3) the track is fitted to a circle, and
+ 4) the track is added to the event as a list of 'SODTrack',
+ which implements 'org.lcsim.event.Track'.
+
lcsim/src/org/lcsim/contrib/SODTracker
diff -N SODFittedCir.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SODFittedCir.java 10 Mar 2006 21:21:34 -0000 1.1
@@ -0,0 +1,383 @@
+package org.lcsim.contrib.SODTracker;
+
+import java.util.*;
+
+public class SODFittedCir extends SODHel
+ {
+ public SODFittedCir()
+ {
+// new
+ current_hel=new SODHel();
+ next_hel=new SODHel();
+ _hitlist=new java.util.LinkedList();
+ nhits=0;
+ }
+ public void make_fitted_cir(java.util.LinkedList l,SODHel h)
+ {
+ make_hel_with_ref(h.D0(),h.Phi0(),h.Omega(),h.Z0(),h.Tanl(),h.Xref(),h.Yref());
+ current_hel=h;
+ next_hel=h;
+// current_hel.make_hel(h.D0(),h.Phi0(),h.Omega(),h.Z0(),h.Tanl());
+// current_hel.change_ref(xnew,ynew);
+ chisq=1000000.0;
+ for (ListIterator ii = l.listIterator();ii.hasNext();){
+ SODHit shit=(SODHit)ii.next(); _hitlist.add(shit); nhits++;
+ }
+ fail=1300;
+ fail=IterateFit();
+ }
+ public void clone_fitted_cir(SODFittedCir fc)
+ {
+ make_hel_with_ref(fc.D0(),fc.Phi0(),fc.Omega(),fc.Z0(),fc.Tanl(),fc.Xref(),fc.Yref());
+ current_hel.make_hel(fc.D0(),fc.Phi0(),fc.Omega(),fc.Z0(),fc.Tanl());
+ next_hel.make_hel(fc.D0(),fc.Phi0(),fc.Omega(),fc.Z0(),fc.Tanl());
+ chisq=fc.Chisq();
+ fail=fc.Fail();
+ java.util.LinkedList l=fc.hitlist();
+ for (ListIterator ii = l.listIterator();ii.hasNext();){
+ SODHit shit=(SODHit)ii.next(); _hitlist.add(shit); nhits++;
+ }
+ }
+ public void add_hit_to_fit(SODHit h)
+ {
+ current_hel.make_hel_with_ref(D0(),Phi0(),Omega(),Z0(),Tanl(),Xref(),Yref());
+ next_hel=current_hel;
+ _hitlist.add(h); nhits++;
+ fail=IterateFit();
+// if (0!=fail){System.out.print("IterateFit FAILED ");System.out.println(fail);}
+ }
+
+public int Fail(){return fail;}
+
+public int IterateFit(){
+ int ftemp=1301; // not enough hits
+ int prtflg=0;
+ if(nfree>=nhits) {return ftemp;}
+ int ndof=nhits-nfree;
+// ftemp=0;
+// if(niter>=1) {
+// float prevchisq=0.0;
+// for (int i=0; i< niter; i++) {
+// itofit=i+1;
+// ftemp=DoFit();
+// if(ftemp!=0) {break;}
+// if(fabs(chisq-prevchisq)<0.01*chisq) {break;}
+// prevchisq=chisq;
+// }//endof iter loop
+// }else{
+// float prevchisq=0.0;
+// chisq=1000000.0;
+// int iter=0;
+// while(fabs(chisq-prevchisq)>0.01) {
+// iter++;
+// prevchisq=chisq;
+// ftemp=DoFit();
+// if(ftemp!=0) break;
+// if(iter>=1000) break;
+// }//endof (fabs(chisq-oldchisq).gt.0.01)
+// }//endof (niter>=1)
+// int ndof=nhits-nfree;
+// prob=Dcxprobab(ndof,chisq);
+// rcs=chisq/ndof;
+ int max_converge=15;
+ int cnt_converge=50;
+ double[] chisq_mem=new double[max_converge];
+ for (int iloop=0; iloop<max_converge; ++iloop){chisq_mem[iloop]=1000000.0;}
+ for (int iloop=0; iloop<max_converge; ++iloop){
+ ftemp=DoFit(); chisq_mem[iloop]=chisq;
+ if (0 != ftemp)break;
+// if ((iloop>0)&&(chisq_mem[iloop]>chisq_mem[iloop-1])){System.out.println("FIT DIVERGES");}
+ if ((iloop>0)&&(chisq_mem[iloop]<chisq_mem[iloop-1])
+ &&((chisq_mem[iloop-1]-chisq_mem[iloop])<0.01)){
+// System.out.println("FIT CONVERGES");
+ cnt_converge=iloop;
+ break;
+ }
+ if ((iloop>0)&&(chisq_mem[iloop]==chisq_mem[iloop-1])){
+// System.out.println("FIT CONVERGES");
+ cnt_converge=iloop;
+ break;
+ }
+ current_hel=next_hel;
+ }
+ prtflg=0;
+ if (prtflg!=0){
+ System.out.print("my_three_chisq");
+ for (int iloop=0; iloop<max_converge; ++iloop){
+ System.out.print(" ");System.out.print(chisq_mem[iloop]);
+ }
+ System.out.println(" ");
+ }
+ prtflg=0;
+ if (prtflg!=0){
+ System.out.print("cnt_converge ");System.out.print(cnt_converge);
+ System.out.print(" ");System.out.print(chisq);
+ System.out.print(" ");System.out.print(ndof);
+// System.out.println(" ");
+ System.out.print(" "); current_hel.printHelix();
+ }
+ if ((0==ftemp)&&(cnt_converge<max_converge)){
+ double xdof=ndof; double rcs=chisq/xdof;
+ if (rcs<100.0){
+ make_hel_with_ref(current_hel.D0(),current_hel.Phi0(),current_hel.Omega(),
+ current_hel.Z0(),current_hel.Tanl(),current_hel.Xref(),current_hel.Yref());
+ }else{
+ fail=1304; // rcs too big
+ chisq=chisq_mem[0]; //set chisq to that of input helix
+ }
+ }else{
+ if (0==ftemp)fail=1303; //fit didn't converge in max_converge steps
+ chisq=chisq_mem[0]; //set chisq to that of input helix
+ }
+ return ftemp;
+}//endof IterateFit
+
+public int DoFit(){
+ int ftemp=1302; //matrix didn't invert
+ int prtflg=0;
+ chisq=0.0;
+ norder=3;
+// double[][] A = new double[3][3];
+ double[][] A = new double[10][10];
+ for (int ii=0; ii<10; ii++){for (int jj=0; jj<10; jj++){A[ii][jj]=0.0;}}
+ double[][] Z = new double[3][3];
+ double[] B = new double[3];
+ double[] D = new double[3];
+ for (int ii=0; ii<norder; ii++){
+ D[ii]=0.0; B[ii]=0.0; for (int jj=0; jj<norder; jj++){A[ii][jj]=0.0;Z[ii][jj]=0.0;}
+ }
+// System.out.println(" ");
+ for (ListIterator ii = _hitlist.listIterator(); ii.hasNext();){
+ SODHit shit=(SODHit)ii.next();
+ double [] derivs=new double[4];
+ if ( shit.calc_circ_derivatives(current_hel,derivs) ){
+ prtflg=0;
+ if (prtflg!=0){
+ System.out.print("derivs ");
+ for (int jj=0; jj<4; ++jj){System.out.print(" ");System.out.print(derivs[jj]);}
+ System.out.println(" ");
+ }//end prtflg
+ double [] numb_derivs=new double[4];
+ boolean whaa=shit.calc_numb_derivatives(current_hel,numb_derivs);
+ prtflg=0;
+ if (prtflg!=0){
+ System.out.print("nderiv ");
+ for (int jj=0; jj<4; ++jj){System.out.print(" ");System.out.print(numb_derivs[jj]);}
+ System.out.println(" ");
+ }
+ chisq+=derivs[0]*derivs[0];
+ for(int ipar=0; ipar<norder; ipar++){
+// D[ipar]+=derivs[0]*derivs[ipar+1];
+ D[ipar]+=numb_derivs[0]*numb_derivs[ipar+1];
+ //inner parameter loop
+ for(int jpar=0; jpar<norder; jpar++){
+// A[ipar][jpar]+=derivs[ipar+1]*derivs[jpar+1];
+ A[ipar][jpar]+=numb_derivs[ipar+1]*numb_derivs[jpar+1];
+ }//endof inner parameter loop
+ }//endof outer parameter loop
+ }// calc_derivatives succeeded
+ }//endof pointloop
+ prtflg=0;
+ if (prtflg!=0){
+ System.out.print("chisq before fit ");System.out.print(chisq);
+// System.out.println(" ");
+ System.out.print(" "); current_hel.printHelix();
+ }
+// A[0][0]=1.;A[0][1]=1.;A[1][1]=100.;A[0][2]=-1.;A[1][2]=2.;A[2][2]=10000.;
+// A[1][0]=A[0][1];A[2][0]=A[0][2];A[2][1]=A[1][2];
+ for (int ii=0; ii<norder; ii++){for (int jj=0; jj<norder; jj++){Z[ii][jj]=A[ii][jj];}}
+ SODmatinv iflg= new SODmatinv(); double det=0.0; iflg.matinv(A,norder,det);
+// change impedence match int iflg=matinv(A,Z); ftemp=iflg;
+// double det=0.0; int iflg=matinv(A,norder,det); ftemp=iflg;
+ prtflg=0;
+ if (prtflg!=0){
+ for (int ii=0; ii<norder; ii++){
+ System.out.print("A ");System.out.print(A[ii][0]);
+ System.out.print(" ");System.out.print(A[ii][1]);
+ System.out.print(" ");System.out.println(A[ii][2]);
+ }
+ for (int ii=0; ii<norder; ii++){
+ System.out.print("Z ");System.out.print(Z[ii][0]);
+ System.out.print(" ");System.out.print(Z[ii][1]);
+ System.out.print(" ");System.out.println(Z[ii][2]);
+ }
+ }//end prtflg
+//if matrix inverted, calc new chisq
+ if ( check_inversion(A,Z) ){
+ ftemp=0;
+ for(int ii=0;ii<norder;ii++){for(int jj=0;jj<norder;jj++){B[ii]+=A[ii][jj]*D[jj];}}
+ SODHel new_hel=new SODHel();
+ new_hel.make_hel(current_hel.D0()-0.5*B[0],current_hel.Phi0()-0.5*B[1],
+ current_hel.Omega()-0.5*B[2],current_hel.Z0(),current_hel.Tanl());
+// new_hel.printHelix();
+ next_hel=new_hel;
+ double new_chisq=0.0;
+// for (int ii=0; ii<norder; ii++){
+// D[ii]=0.0; B[ii]=0.0; for (int jj=0; jj<norder; jj++){A[ii][jj]=0.0;Z[ii][jj]=0.0;}
+// }
+// System.out.println(" ");
+// for (ListIterator ii = _hitlist.listIterator(); ii.hasNext();){
+// SODHit shit=(SODHit)ii.next();
+// double [] derivs=new double[4];
+// if ( shit.calc_circ_derivatives(new_hel,derivs) ){
+// prtflg=1;
+// if (prtflg!=0){
+// System.out.print("derivs ");System.out.print(derivs[0]);
+// System.out.print(" ");System.out.print(derivs[1]);
+// System.out.print(" ");System.out.print(derivs[2]);
+// System.out.print(" ");System.out.println(derivs[3]);
+// }//end prtflg
+// new_chisq+=derivs[0]*derivs[0];
+// for(int ipar=0; ipar<norder; ipar++){
+// D[ipar]+=derivs[0]*derivs[ipar+1];
+// //inner parameter loop
+// for(int jpar=0; jpar<norder; jpar++){
+// A[ipar][jpar]+=derivs[ipar+1]*derivs[jpar+1];
+// }//endof inner parameter loop
+// }//endof outer parameter loop
+// }// calc_derivatives succeeded
+// }//endof pointloop
+// System.out.print("chisq before fit ");System.out.print(new_chisq);
+// System.out.println(" ");
+// System.out.print(" "); new_hel.printHelix();
+// if (new_chisq<chisq)System.out.println("FITTING!!! ");
+//try again!
+// for (int ii=0; ii<norder; ii++){for (int jj=0; jj<norder; jj++){Z[ii][jj]=A[ii][jj];}}
+// SODmatinv jflg= new SODmatinv(); det=0.0; jflg.matinv(A,norder,det);
+// if ( check_inversion(A,Z) ){
+// for(int ii=0;ii<norder;ii++){for(int jj=0;jj<norder;jj++){B[ii]+=A[ii][jj]*D[jj];}}
+// SODHel newest_hel=new SODHel();
+// newest_hel.make_hel(new_hel.D0()-0.5*B[0],new_hel.Phi0()-0.5*B[1],new_hel.Omega()-0.5*B[2],
+// new_hel.Z0(),new_hel.Tanl());
+// newest_hel.printHelix();
+// double next_chisq=0.0;
+// System.out.println(" ");
+// for (ListIterator ii = _hitlist.listIterator(); ii.hasNext();){
+// SODHit shit=(SODHit)ii.next();
+// double [] derivs=new double[4];
+// if ( shit.calc_circ_derivatives(newest_hel,derivs) ){
+// prtflg=1;
+// if (prtflg!=0){
+// System.out.print("derivs ");System.out.print(derivs[0]);
+// System.out.print(" ");System.out.print(derivs[1]);
+// System.out.print(" ");System.out.print(derivs[2]);
+// System.out.print(" ");System.out.println(derivs[3]);
+// }//end prtflg
+// next_chisq+=derivs[0]*derivs[0];
+// }// calc_derivatives succeeded
+// }//endof pointloop
+// System.out.print("my_three_chisq");
+// System.out.print(" ");System.out.print(current_hel.D0() );
+// System.out.print(" ");System.out.print(new_hel.D0() );
+// System.out.print(" ");System.out.print(newest_hel.D0() );
+// System.out.print(" ");System.out.print(current_hel.Phi0() );
+// System.out.print(" ");System.out.print(current_hel.Omega() );
+// System.out.print(" ");System.out.print(chisq);
+// System.out.print(" ");System.out.print(new_chisq);
+// System.out.print(" ");System.out.println(next_chisq);
+// if (next_chisq<new_chisq)System.out.println("STILL FITTING!!! ");
+// }//end next matrix inverted
+ }//end matrix inverted
+ return ftemp;
+}
+
+ int matinv(double[][] A, int snorder, double sdet){
+ int iflg=0;
+ double[] amat = new double[6];
+ amat[0]=A[0][0];amat[1]=A[0][1];amat[2]=A[1][1];
+ amat[3]=A[0][2];amat[4]=A[1][2];amat[5]=A[2][2];
+ double[] ainv = new double[6];
+ sdet=dchThreeInv(amat,ainv);
+ if (0.0!=sdet)iflg=1;
+ A[0][0]=ainv[0];A[0][1]=ainv[1];A[1][1]=ainv[2];
+ A[0][2]=ainv[3];A[1][2]=ainv[4];A[2][2]=ainv[5];
+ A[1][0]=ainv[1];A[2][0]=ainv[3];A[2][1]=ainv[4];
+ return iflg;
+ }
+
+// int dchThreeInv( double matrix[6], double invmat[6]) {
+ double dchThreeInv( double[] matrix, double[] invmat) {
+
+ /* Declare variables. */
+ double det, detinv;
+// double cofactor[6];
+ double[] cofactor = new double[6];
+
+ /**************************************************************************/
+
+// System.out.print("matrix ");System.out.print(matrix[0]);
+// System.out.print(" ");System.out.print(matrix[1]);
+// System.out.print(" ");System.out.print(matrix[2]);
+// System.out.print(" ");System.out.print(matrix[3]);
+// System.out.print(" ");System.out.print(matrix[4]);
+// System.out.print(" ");System.out.println(matrix[5]);
+ cofactor[0] = (matrix[2]*matrix[5] - matrix[4]*matrix[4]);
+ cofactor[1] = -(matrix[1]*matrix[5] - matrix[3]*matrix[4]);
+ cofactor[3] = (matrix[1]*matrix[4] - matrix[2]*matrix[3]);
+ det = matrix[0] * cofactor[0] +
+ matrix[1] * cofactor[1] +
+ matrix[3] * cofactor[3];
+ if (det == 0.0) return det;
+ detinv = 1./det;
+ cofactor[2] = (matrix[0]*matrix[5] - matrix[3]*matrix[3]);
+ cofactor[4] = -(matrix[0]*matrix[4] - matrix[1]*matrix[3]);
+ cofactor[5] = (matrix[0]*matrix[2] - matrix[1]*matrix[1]);
+
+ invmat[0] = cofactor[0] * detinv;
+ invmat[1] = cofactor[1] * detinv;
+ invmat[2] = cofactor[2] * detinv;
+ invmat[3] = cofactor[3] * detinv;
+ invmat[4] = cofactor[4] * detinv;
+ invmat[5] = cofactor[5] * detinv;
+// System.out.print("invmat ");System.out.print(invmat[0]);
+// System.out.print(" ");System.out.print(invmat[1]);
+// System.out.print(" ");System.out.print(invmat[2]);
+// System.out.print(" ");System.out.print(invmat[3]);
+// System.out.print(" ");System.out.print(invmat[4]);
+// System.out.print(" ");System.out.println(invmat[5]);
+
+ return det;
+ }
+ boolean check_inversion( double[][] A, double[][] Z) {
+ int prtflg=0;
+ boolean rtflg=false; int norder=3;
+ double[][] I = new double[3][3];
+ for(int ipar=0; ipar<norder; ipar++){
+ for(int jpar=0; jpar<norder; jpar++){
+ I[ipar][jpar]=A[ipar][0]*Z[0][jpar]+A[ipar][1]*Z[1][jpar]+A[ipar][2]*Z[2][jpar];
+ }
+ }
+ if ((Math.abs(I[0][0]-1.0)<epsilon)&&(Math.abs(I[1][1]-1.0)<epsilon)&&(Math.abs(I[2][2]-1.0)<epsilon)&&
+ (Math.abs(I[0][1] )<epsilon)&&(Math.abs(I[0][2] )<epsilon)&&(Math.abs(I[1][2] )<epsilon))rtflg=true;
+// if (rtflg)System.out.println("MATRIX INVERTED!!!!!");
+ prtflg=0;
+ if (prtflg!=0){
+ for (int ii=0; ii<norder; ii++){
+ System.out.print("I ");System.out.print(I[ii][0]);
+ System.out.print(" ");System.out.print(I[ii][1]);
+ System.out.print(" ");System.out.println(I[ii][2]);
+ }
+ }//end prtflg
+ return rtflg;
+ }
+ public void printFit()
+ {
+ System.out.print(" ");System.out.print(chisq);
+ System.out.print(" ");System.out.print(nhits);
+ System.out.print(" "); current_hel.printHelix();
+ }
+ public double Chisq(){return chisq;}
+ public int Nhits(){return nhits;}
+ public java.util.LinkedList hitlist(){return _hitlist;}
+ private double chisq;
+ private int norder;
+ private int fail;
+ private int nhits;
+ private static int nfree=3;
+ private SODHel current_hel;
+ private SODHel next_hel;
+ private java.util.LinkedList _hitlist;
+ private static double pi=3.14159265359;
+ private static double twopi=6.283185307;
+ private static double epsilon=0.000000001;
+}
lcsim/src/org/lcsim/contrib/SODTracker
diff -N SODHel.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SODHel.java 10 Mar 2006 21:21:34 -0000 1.1
@@ -0,0 +1,206 @@
+package org.lcsim.contrib.SODTracker;
+
+import org.lcsim.util.step.TrkParams;
+import java.util.*;
+
+public class SODHel
+{
+ public SODHel()
+ {
+ }
+ public void make_hel(TrkParams t)
+ {
+ d0=t.getParam(0);phi0=t.getParam(1);omega=t.getParam(2);
+ z0=t.getParam(3);tanl=t.getParam(4);q=1;if(omega<0.0)q=-1;
+ aomega=Math.abs(omega);sphi0=Math.sin(phi0);cphi0=Math.cos(phi0);
+ ominfl=true; if(omin>aomega){ominfl=false;}
+// xc=-sphi0*(d0+1.0/omega);yc=cphi0*(d0+1.0/omega);
+ x0=-sphi0*d0; y0=cphi0*d0; xc=x0-sphi0/omega; yc=y0+cphi0/omega;
+ }
+ public void make_hel(double d,double p,double o,double z,double t)
+ {
+ d0=d;phi0=p;omega=o;z0=z;tanl=t;q=1;if(omega<0.0)q=-1;
+ if (phi0>pi)phi0-=twopi; if (phi0<(-pi))phi0+=twopi;
+ aomega=Math.abs(omega);sphi0=Math.sin(phi0);cphi0=Math.cos(phi0);
+ ominfl=true; if(omin>aomega){ominfl=false;}
+// xc=-sphi0*(d0+1.0/omega);yc=cphi0*(d0+1.0/omega);
+// x0=-sphi0*d0; y0=cphi0*d0; xc=x0-sphi0/omega; yc=y0+cphi0/omega;
+ xref=0.0; yref=0.0; x0=X0(); y0=Y0(); xc=Xc(); yc=Yc();
+ }
+ public void make_hel_with_ref(double d,double p,double o,double z,double t, double xnew, double ynew)
+ {
+ d0=d;phi0=p;omega=o;z0=z;tanl=t;q=1;if(omega<0.0)q=-1;
+ aomega=Math.abs(omega);sphi0=Math.sin(phi0);cphi0=Math.cos(phi0);
+ ominfl=true; if(omin>aomega){ominfl=false;}
+// xc=-sphi0*(d0+1.0/omega);yc=cphi0*(d0+1.0/omega);
+// x0=-sphi0*d0; y0=cphi0*d0; xc=x0-sphi0/omega; yc=y0+cphi0/omega;
+ xref=xnew; yref=xnew; x0=X0(); y0=Y0(); xc=Xc(); yc=Yc();
+ }
+ public void change_ref(double xnew,double ynew)
+ {
+// the helix is always constructed with (xref,yref)=(0,0) but may
+// need to change it to improve matrix inversion problems
+ double rnew=Math.sqrt(xnew*xnew+ynew*ynew);
+ double phinew=Math.atan2(ynew,xnew);
+ double phi1 = phi1_at_radius(rnew);
+ double dphi=phinew-phi1;
+ if(dphi<-pi)dphi+=twopi;if(dphi>pi)dphi-=twopi;
+ _resid=dphi*rnew; d0=_resid;
+ phinew=phi0+omega*len; phi0=phinew;
+ xref=xnew; yref=ynew;
+ x0=X0(); y0=Y0(); xc=Xc(); yc=Yc();
+ printHelix();
+// System.out.print("change xref yref ");System.out.print(xnew);
+// System.out.print(" ");System.out.println(ynew);
+ }
+ public double Xref(){return xref;}
+ public double Yref(){return yref;}
+ public void printHelix()
+ {
+// System.out.print("SODHel ");
+ System.out.print(d0);
+ System.out.print(" ");System.out.print(phi0);
+ System.out.print(" ");System.out.print(omega);
+ System.out.print(" ");System.out.print(z0);
+ System.out.print(" ");System.out.print(tanl);
+// System.out.print(" ");System.out.print(xref);
+// System.out.print(" ");System.out.print(yref);
+ System.out.println(" ");
+ }
+ public double D0(){return d0;}
+ public double Phi0(){return phi0;}
+ public double Omega(){return omega;}
+ public double Z0(){return z0;}
+ public double Tanl(){return tanl;}
+ public int Q(){return q;}
+// public double Xc(){return xc;}
+// public double Yc(){return yc;}
+ public double Xc(){return (X0()-sphi0/omega);}
+ public double Yc(){return (Y0()+cphi0/omega);}
+ public double X0(){return (xref-sphi0*d0);}
+ public double Y0(){return (yref+cphi0*d0);}
+ public double Xh(double l){
+ if(ominfl){double phit=phi0+omega*l;return (xc+Math.sin(phit)/omega);}
+ else{return (-sphi0*d0+cphi0*l-0.5*l*l*omega*sphi0);}//ominfl
+ }//endof Xh
+ public double Yh(double l){
+ if(ominfl){double phit=phi0+omega*l;return (yc-Math.cos(phit)/omega);
+ }else{return (cphi0*d0+sphi0*l+0.5*l*l*omega*cphi0);}//ominfl
+ }//endof Yh
+ public double Zh(double l){return (z0+tanl*l);}//endof Zh
+ public double Px(double l){
+ if(ominfl){double phit=phi0+omega*l; return ptcon*Math.cos(phit)/aomega;}
+ else{return ptmax*cphi0;}//ominfl
+ }//endof Px
+ public double Py(double l){
+ if(ominfl){double phit=phi0+omega*l; return ptcon*Math.sin(phit)/aomega;}
+ else{return ptmax*sphi0;}//ominfl
+ }//endof Py
+ public double Pz(){
+ if(ominfl){return ptcon*tanl/aomega;}
+ else{return ptmax*tanl;}//ominfl
+ }//endof Pz
+ public double Ptot(){
+ if(ominfl){return ptcon*Math.sqrt(1.0+tanl*tanl)/aomega;}
+ else{return ptmax*Math.sqrt(1.0+tanl*tanl);}//ominfl
+ }//endof Ptot
+ public double Pt(){
+ if(ominfl){return ptcon/aomega;}
+ else{return ptmax;}//ominfl
+ }//endof Pt
+ public double phi1_at_radius(double rl)
+ {
+// System.out.print("xref yref ");System.out.print(xref);
+// System.out.print(" ");System.out.println(yref);
+ if ((0.0==xref)&&(0.0==yref)){
+ double orcsq=(1.0+d0*omega)*(1.0+d0*omega);
+ double orlsq=(rl*omega)*(rl*omega);
+ double cphil=(1.0+orcsq-orlsq)/(2.0*(1.0+d0*omega));
+ double phil1=Math.acos(cphil);
+ double l1=phil1/omega; if (l1<0.0)l1=-1.0*l1;
+ double phit=phi0+phil1; if (omega<0.0)phit=phi0-phil1;
+// xh=xc+Math.sin(phit)/omega;
+// yh=yc-Math.cos(phit)/omega;
+ xh=Xh(l1); yh=Yh(l1);
+ double rh=Math.sqrt(xh*xh+yh*yh);
+ double phih=Math.atan2(yh,xh);
+ len=l1;
+// System.out.println("or here?");
+ return phih;
+ }else{
+// System.out.println("am i here?");
+ xh=0.0; yh=0.0; len=0.0; double phih=0.0; return phih;
+ }
+ }
+ public boolean calc_circ_derivatives(double rl,double phihit,double [] derivs)
+ {
+ double phi1 = phi1_at_radius(rl);
+ double dphi=phihit-phi1;
+ if(dphi<-pi)dphi+=twopi;if(dphi>pi)dphi-=twopi;
+ _resid=dphi*rl;
+ derivs[0]= _resid;
+// derivs[1]= 0.9*.0140;derivs[2]= 1.0*.0140;derivs[3]= 1.2*.0140;
+// impedence matching (lots of circle and phi strip approxes)
+ double phi=phi1;
+ double fac=1.0;
+ if(_resid<0.0) {fac=-fac;}
+ if(d0<0.0) {fac=-fac;}
+ double vx=-Math.sin(phi), vy=Math.cos(phi);
+ double tx= Math.cos(phi), ty=Math.sin(phi); //am I sure???
+ double f0=0.0;
+//
+ double dddd0=-vx*sphi0+vy*cphi0; derivs[1]=dddd0*fac;
+ double dddp0=-(yh-y0+f0*ty)*vx+(xh-x0+f0*tx)*vy;
+ dddp0=dddp0*(1.0+d0*omega);
+ derivs[2]=dddp0*fac;
+ double dddom;
+//? if (ominfl){
+//? dddom=((len*Math.cos(phi)-xh+x0)*vx+(len*Math.sin(phi)-yh+y0)*vy)/omega;
+//? dddom+=f0*len*cosl*(-Math.sin(phi)*vx+Math.cos(phi)*vy);
+//? }else{
+ dddom=0.5*len*len*(-vx*sphi0+vy*cphi0);
+//??? dddom=0.5*len*len*dddd0;
+// dddom=0.5*len*len*(-Math.sin(phi)*vx+Math.cos(phi)*vy);
+//ght double chsq=(xh-x0)*(xh-x0)+(yh-y0)*(yh-y0);
+//ght dddom=0.5*chsq*(ty*vx-tx*vy);
+// double xdddom1=((len*Math.cos(phi)-xh+x0)*vx+(len*Math.sin(phi)-yh+y0)*vy)/omega;
+// double xdddom2=len*(-Math.sin(phi)*vx+Math.cos(phi)*vy);
+// f0=(dddom-xdddom1)/xdddom2;
+// System.out.print("dddom ");System.out.print(dddom);
+// System.out.print(" ");System.out.print(xdddom1);
+// System.out.print(" ");System.out.print(xdddom2);
+// System.out.print(" ");System.out.println(f0);
+// double dang=Math.cos(phi)*Math.sin(phihit)-Math.sin(phi)*Math.cos(phihit);
+// System.out.print("dang omega ");System.out.print(dang);
+// System.out.print(" ");System.out.println(omega);
+//? }
+ derivs[3]=dddom*fac;
+// derivs[1]=derivs[1]/100.;derivs[2]=derivs[2]/10000.;derivs[3]=derivs[3]/10000.;
+ return true;
+ }
+ private double d0;
+ private double phi0;
+ private double omega;
+ private double z0;
+ private double tanl;
+ private double aomega;
+ private double sphi0;
+ private double cphi0;
+ private int q;
+ private boolean ominfl;
+ private static double omin=0.000005;
+ private static double pi=3.14159265359;
+ private static double twopi=6.283185307;
+ private static double ptmax=3000.0;
+ private static double ptcon=0.015;
+ private double xref;
+ private double yref;
+ private double x0;
+ private double y0;
+ private double xc;
+ private double yc;
+ private double _resid;
+ private double len;
+ private double xh;
+ private double yh;
+}
lcsim/src/org/lcsim/contrib/SODTracker
diff -N SODHit.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SODHit.java 10 Mar 2006 21:21:34 -0000 1.1
@@ -0,0 +1,228 @@
+package org.lcsim.contrib.SODTracker;
+
+import org.lcsim.event.TrackerHit;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.MCParticle;
+import org.lcsim.util.step.TrkParams;
+import java.util.*;
+
+public class SODHit
+{
+ public SODHit()
+ {
+ _used_on_hel=-1;
+ _resid=-10000.0;
+ _is_foreground_hit=false;
+ }
+// public SODHit(TrackerHit h)
+ public void make_hit(TrackerHit h)
+ {
+ _used_on_hel=0;
+ _is_foreground_hit=true;
+//jas2 double[] hitpos = h.getPoint();
+ double[] hitpos = h.getPosition();
+ _hitx = hitpos[0];
+ _hity = hitpos[1];
+ _hitrad = Math.sqrt(_hitx*_hitx+_hity*_hity);
+ _hitlay=0;
+// if _hitlay=0 it's not a SODHIT and other data may be invalid,
+// SO CHECK ON IT!!!
+ if ((_hitrad>20.00)&&(_hitrad<20.10))_hitlay=1;
+ if ((_hitrad>46.25)&&(_hitrad<46.35))_hitlay=2;
+ if ((_hitrad>72.50)&&(_hitrad<72.60))_hitlay=3;
+ if ((_hitrad>98.75)&&(_hitrad<98.85))_hitlay=4;
+ if ((_hitrad>125.0)&&(_hitrad<125.1))_hitlay=5;
+// not really true right now (let sit hits be SODHit's also and
+// protect "elsewhere")
+ if ((_hitrad>1.15)&&(_hitrad<1.25))_hitlay=1;
+ if ((_hitrad>2.35)&&(_hitrad<2.45))_hitlay=2;
+ if ((_hitrad>3.55)&&(_hitrad<3.65))_hitlay=3;
+ if ((_hitrad>4.75)&&(_hitrad<4.85))_hitlay=4;
+ if ((_hitrad>5.95)&&(_hitrad<6.05))_hitlay=5;
+ _hitphi = Math.atan2(_hity,_hitx);
+// System.out.print("drawmetrk ");System.out.print(_hitx);
+// System.out.print(" ");System.out.println(_hity);
+// TrackerHit _savehit = h;
+ }
+ public void make_hit(int lay, double xin, double yin, double rgauss, SimTrackerHit trackerHit, MCParticle mcp)
+ {
+ _used_on_hel=0;
+ _is_foreground_hit=true;
+ _hitx = xin*mmtocm;
+ _hity = yin*mmtocm;
+ _hitrad = Math.sqrt(_hitx*_hitx+_hity*_hity);
+ _hitlay=lay+1;
+ _hitphi = Math.atan2(_hity,_hitx);
+ double myres = sodres; if (_hitrad<10.0){myres = vtxres;}
+ _hitx-=rgauss*myres*Math.sin(_hitphi);
+ _hity+=rgauss*myres*Math.cos(_hitphi);
+ _hitphi = Math.atan2(_hity,_hitx);
+ _trackerHit = trackerHit;
+ _mcParticle = mcp;
+ }
+ public void make_hit(int lay, double xin, double yin, double rgauss, SimTrackerHit trackerHit)
+ {
+ _used_on_hel=0;
+ _is_foreground_hit=true;
+ _hitx = xin*mmtocm;
+ _hity = yin*mmtocm;
+ _hitrad = Math.sqrt(_hitx*_hitx+_hity*_hity);
+ _hitlay=lay+1;
+ _hitphi = Math.atan2(_hity,_hitx);
+ double myres = sodres; if (_hitrad<10.0){myres = vtxres;}
+ _hitx-=rgauss*myres*Math.sin(_hitphi);
+ _hity+=rgauss*myres*Math.cos(_hitphi);
+ _hitphi = Math.atan2(_hity,_hitx);
+ _trackerHit = trackerHit;
+ }
+ public void make_hit(double xrad, double xphi)
+ {
+ _used_on_hel=0;
+ _is_foreground_hit=false;
+ _hitx = xrad*Math.cos(xphi);
+ _hity = xrad*Math.sin(xphi);
+ _hitrad = xrad;
+ _hitlay=0;
+// if _hitlay=0 it's not a SODHIT and other data may be invalid,
+// SO CHECK ON IT!!!
+ if ((_hitrad>20.00)&&(_hitrad<20.10))_hitlay=1;
+ if ((_hitrad>46.25)&&(_hitrad<46.35))_hitlay=2;
+ if ((_hitrad>72.50)&&(_hitrad<72.60))_hitlay=3;
+ if ((_hitrad>98.75)&&(_hitrad<98.85))_hitlay=4;
+ if ((_hitrad>125.0)&&(_hitrad<125.1))_hitlay=5;
+ _hitphi = xphi;
+// System.out.print("drawmebg ");System.out.print(_hitx);
+// System.out.print(" ");System.out.println(_hity);
+// TrackerHit _savehit = h;
+ }
+ public int GetUsedOnHel(){return _used_on_hel;}
+ public void SetUsedOnHel(int i){_used_on_hel=i;}
+ public int layer(){return _hitlay;}
+ public double x(){return _hitx;}
+ public double y(){return _hity;}
+ public double r(){return _hitrad;}
+ public double phi(){return _hitphi;}
+// the following only make sense after calc_resid (test on true)
+ public double get_resid(){return _resid;}
+ public boolean calc_circ_derivatives(SODHel t, double [] derivs)
+ {
+// calc_resid(t); derivs[0]= get_resid()/e();
+ t.calc_circ_derivatives(_hitrad,_hitphi,derivs);
+ double estrip=e();
+// c++ for(int i=0; i<derivs.length(); i++) {derivs[i]/=estrip;}
+ for(int i=0; i<4; i++) {derivs[i]/=estrip;}
+// derivs[0]= derivs[0]/e();
+ return true;
+ }
+ public boolean calc_numb_derivatives(SODHel t, double [] derivs)
+ {
+ calc_resid(t); double resid_base=get_resid(); double error=e();
+ derivs[0]=resid_base/error; double delta=0.0001;
+ SODHel new_hel=new SODHel();
+ new_hel.make_hel(t.D0()+delta,t.Phi0(),t.Omega(),t.Z0(),t.Tanl());
+ calc_resid(new_hel);derivs[1]=(get_resid()-resid_base)/(delta*error);
+ new_hel.make_hel(t.D0(),t.Phi0()+delta,t.Omega(),t.Z0(),t.Tanl());
+ calc_resid(new_hel);derivs[2]=(get_resid()-resid_base)/(delta*error);
+ delta=0.01*t.Omega();
+ new_hel.make_hel(t.D0(),t.Phi0(),t.Omega()+delta,t.Z0(),t.Tanl());
+ calc_resid(new_hel);derivs[3]=(get_resid()-resid_base)/(delta*error);
+ return true;
+ }
+ public boolean calc_resid(SODHel t)
+ {
+ double phi1 = t.phi1_at_radius(_hitrad);
+ double dphi=_hitphi-phi1;
+ if(dphi<-pi)dphi+=twopi;if(dphi>pi)dphi-=twopi;
+ _resid=dphi*_hitrad;
+ int prtflg=0;
+ if (prtflg!=0){
+ System.out.print("hittup ");System.out.print(_hitlay);
+ System.out.print(" ");System.out.print(_hitrad);
+ System.out.print(" ");System.out.print(_resid);
+ System.out.print(" ");System.out.print(t.D0());
+ System.out.print(" ");System.out.print(t.Phi0());
+ System.out.print(" ");System.out.println(t.Omega());
+ }//end printing
+ return true;
+ }
+ public boolean calc_resid(TrkParams t)
+ {
+ double phi1 = phi1_at_radius(t,_hitrad);
+ double dphi=_hitphi-phi1;
+ if(dphi<-pi)dphi+=twopi;if(dphi>pi)dphi-=twopi;
+ _resid=dphi*_hitrad;
+ int prtflg=0;
+ if (prtflg!=0){
+ System.out.print("hittup ");System.out.print(_hitlay);
+ System.out.print(" ");System.out.print(_hitrad);
+ System.out.print(" ");System.out.print(_resid);
+ System.out.print(" ");System.out.print(t.getParam(0));
+ System.out.print(" ");System.out.print(t.getParam(1));
+ System.out.print(" ");System.out.println(t.getParam(2));
+ }//end printing
+ return true;
+ }
+ public double phi1_at_radius(TrkParams t, double rl)
+ {
+ double d0=t.getParam(0);
+ double phi0=t.getParam(1);
+ double omega=t.getParam(2);
+ double orcsq=(1.0+d0*omega)*(1.0+d0*omega);
+ double orlsq=(rl*omega)*(rl*omega);
+ double cphil=(1.0+orcsq-orlsq)/(2.0*(1.0+d0*omega));
+ double phil1=Math.acos(cphil);
+ double l1=phil1/omega; if (l1<0.0)l1=-1.0*l1;
+ double sp0=Math.sin(phi0); double cp0=Math.cos(phi0);
+ double xc=-sp0*(d0+1.0/omega);
+ double yc=cp0*(d0+1.0/omega);
+ double phit=phi0+phil1; if (omega<0.0)phit=phi0-phil1;
+ double xh=xc+Math.sin(phit)/omega;
+ double yh=yc-Math.cos(phit)/omega;
+ double rh=Math.sqrt(xh*xh+yh*yh);
+ double phih=Math.atan2(yh,xh);
+ return phih;
+ }
+ public double e(double pt){
+ double res = -0.00040625*pt+0.0208125;
+ if ((pt<2.0)&&(_hitlay==1))res=-0.01*pt+0.04;
+ if ((pt<2.0)&&(_hitlay==5))res=-0.18*pt+0.38;
+ if ((res<0.0130)&&(_hitlay==1))res=0.0130;
+ if ((res<0.0141)&&(_hitlay==5))res=0.0141;
+ return res;
+ }
+ public double e(){
+//fundamental res double res = 0.0007;
+// double res = 0.0140;
+//used for slac workshop double res = 0.0060;
+//
+//c for now change resolution by hand in SODHit e()
+//c to compensate for MS
+//
+// double res = 0.0020; //new012904 for 50 GeV 10 hits
+// double res = 0.0100; //new020204 for 1 GeV 10 hits
+// double res = 0.0040; //new081905 for 2 GeV 10 hits
+// double res = 0.0024; //new081905 for 5 GeV 10 hits
+// double res = 0.0021; //new081905 for 10 GeV 10 hits
+// double res = 0.0010; //new011006 for 20, 50, 100 GeV 10 hits
+ double res = sodres; if (_hitrad<10.0){res = vtxres;}
+ return res;
+ }
+ public boolean IsForegroundHit(){return _is_foreground_hit;}
+ public MCParticle getMCParticle(){return _mcParticle;}
+ public SimTrackerHit getHit(){return _trackerHit;}
+ private int _hitlay;
+ private double _hitx;
+ private double _hity;
+ private double _hitrad;
+ private double _hitphi;
+ private double _resid;
+ private boolean _is_foreground_hit;
+ public int _used_on_hel;
+ private static double pi=3.14159265359;
+ private static double twopi=6.283185307;
+ private static double mmtocm=0.1;
+ private static double vtxres=0.0005;
+ private static double sodres=0.0007;
+ private SimTrackerHit _trackerHit;
+ private MCParticle _mcParticle;
+}
lcsim/src/org/lcsim/contrib/SODTracker
diff -N SODHitAdder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SODHitAdder.java 10 Mar 2006 21:21:34 -0000 1.1
@@ -0,0 +1,104 @@
+package org.lcsim.contrib.SODTracker;
+
+import java.util.*;
+//import org.lcsim.util.aida.AIDA;
+
+public class SODHitAdder
+{
+// private AIDA aida = AIDA.defaultInstance();
+ public SODHitAdder()
+ {
+ out_fit = new SODFittedCir();
+ }
+ public SODHitAdder(SODFittedCir in_fit,
+ java.util.LinkedList l1shlist,
+ java.util.LinkedList l2shlist,
+ java.util.LinkedList l3shlist,
+ java.util.LinkedList l4shlist,
+ java.util.LinkedList l5shlist) {
+ SODFittedCir new_fit = new SODFittedCir();
+ SODL1Path firstpa = new SODL1Path(1,in_fit,l1shlist,l2shlist,l3shlist,l4shlist,l5shlist);
+ int best_nhits=firstpa.Nhits();
+ double best_rcs=firstpa.Rcs();
+ int best_iflag=firstpa.Complete_trk();
+ double best_pt=firstpa.Pt();
+ double previous_rcs=1000000.0;
+ int previous_iflag=-1;
+ new_fit = firstpa.out_fit();
+ SODL1Path nextpa = new SODL1Path(2,in_fit,l1shlist,l2shlist,l3shlist,l4shlist,l5shlist);
+ double next_rcs=nextpa.Rcs();
+ if ((next_rcs<best_rcs)&&(nextpa.Nhits()>=best_nhits)){
+ previous_rcs=best_rcs; previous_iflag=best_iflag;
+ best_pt=nextpa.Pt();
+ best_nhits=nextpa.Nhits();
+ best_rcs=next_rcs;
+ best_iflag=nextpa.Complete_trk();
+ new_fit = nextpa.out_fit();
+ }else{
+ previous_rcs=next_rcs;
+ previous_iflag=nextpa.Complete_trk();
+ }
+ SODL1Path thirdpa = new SODL1Path(3,in_fit,l1shlist,l2shlist,l3shlist,l4shlist,l5shlist);
+ double third_rcs=thirdpa.Rcs();
+ if ((third_rcs<best_rcs)&&(thirdpa.Nhits()>=best_nhits)){
+ previous_rcs=best_rcs; previous_iflag=best_iflag;
+ best_pt=thirdpa.Pt();
+ best_nhits=thirdpa.Nhits();
+ best_rcs=third_rcs;
+ best_iflag=thirdpa.Complete_trk();
+ new_fit = thirdpa.out_fit();
+ }else{
+ if (third_rcs<previous_rcs){
+ previous_rcs=third_rcs;
+ previous_iflag=thirdpa.Complete_trk();
+ }
+ }
+// SODL1Path fourthpa = new SODL1Path(4,in_fit,l1shlist,l2shlist,l3shlist,l4shlist,l5shlist);
+// double fourthpa_rcs=fourthpa.Rcs();
+// if ((fourthpa_rcs<best_rcs)&&(fourthpa.Nhits()>=best_nhits)){
+// previous_rcs=best_rcs; previous_iflag=best_iflag;
+// best_pt=fourthpa.Pt();
+// best_nhits=fourthpa.Nhits();
+// best_rcs=fourthpa_rcs;
+// best_iflag=fourthpa.Complete_trk();
+// new_fit = fourthpa.out_fit();
+// }else{
+// if (fourthpa_rcs<previous_rcs){
+// previous_rcs=fourthpa_rcs;
+// previous_iflag=fourthpa.Complete_trk();
+// }
+// }
+// SODL1Path fifthpa = new SODL1Path(3,in_fit,l1shlist,l2shlist,l3shlist,l4shlist,l5shlist);
+// double fifthpa_rcs=fifthpa.Rcs();
+// if ((fifthpa_rcs<best_rcs)&&(fifthpa.Nhits()>=best_nhits)){
+// previous_rcs=best_rcs; previous_iflag=best_iflag;
+// best_pt=fifthpa.Pt();
+// best_nhits=fifthpa.Nhits();
+// best_rcs=fifthpa_rcs;
+// best_iflag=fifthpa.Complete_trk();
+// new_fit = fifthpa.out_fit();
+// }else{
+// if (fifthpa_rcs<previous_rcs){
+// previous_rcs=fifthpa_rcs;
+// previous_iflag=fifthpa.Complete_trk();
+// }
+// }
+ int prtflg=0;
+ if (prtflg!=0){
+ System.out.print("SODHitAdder");
+ System.out.print(" ");System.out.print(best_rcs);
+ System.out.print(" ");System.out.print(best_nhits);
+ System.out.print(" ");System.out.print(best_iflag);
+ System.out.print(" ");System.out.print(best_pt);
+ System.out.print(" ");System.out.print(previous_rcs);
+ System.out.print(" ");System.out.print(previous_iflag);
+ System.out.println(" ");
+ }
+// aida.cloud1D("Mom").fill(best_pt);
+ out_fit=new_fit;
+ }//end ctor from 5 hit lists
+ public SODFittedCir get_new_fit(){return out_fit;}
+ private SODFittedCir out_fit;
+ private static double pi=3.14159265359;
+ private static double twopi=6.283185307;
+}
lcsim/src/org/lcsim/contrib/SODTracker
diff -N SODL1Path.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SODL1Path.java 10 Mar 2006 21:21:34 -0000 1.1
@@ -0,0 +1,125 @@
+package org.lcsim.contrib.SODTracker;
+
+import java.util.*;
+
+public class SODL1Path
+{
+ public SODL1Path(int npath, SODFittedCir in_fit,
+ java.util.LinkedList l1shlist,
+ java.util.LinkedList l2shlist,
+ java.util.LinkedList l3shlist,
+ java.util.LinkedList l4shlist,
+ java.util.LinkedList l5shlist)
+ {
+ test_fit = new SODFittedCir();
+ complete_trk=0;
+ int prtflg=0;
+ double rmin=10.0;
+ test_fit.clone_fitted_cir(in_fit);
+// see if hits in layers 1 5 are with tolerance
+ double l1_min_resid=rmin,l2_min_resid=rmin;
+ double l3_min_resid=rmin,l4_min_resid=rmin;
+ double l5_min_resid=rmin;
+ SODHit sh1_min=(SODHit)l1shlist.getFirst();
+ SODHit sh2_min=(SODHit)l2shlist.getFirst();
+ SODHit sh3_min=(SODHit)l3shlist.getFirst();
+ SODHit sh4_min=(SODHit)l4shlist.getFirst();
+ SODHit sh5_min=(SODHit)l5shlist.getFirst();
+//layer 1
+ int l1_correct_rank=0;
+ for (ListIterator ii = l1shlist.listIterator();ii.hasNext();){
+ SODHit sh1=(SODHit)ii.next();
+ if ((sh1.GetUsedOnHel()==0)&&(sh1.calc_resid(test_fit))){
+ double resid= sh1.get_resid() ;
+ if ( Math.abs(resid)<Math.abs(l1_min_resid) ){
+ l1_min_resid=resid;
+ sh1_min=sh1;
+ l1_correct_rank++;
+ if (sh1.IsForegroundHit())l1_correct_rank=0;
+ }
+ }
+ }
+ if ( Math.abs(l1_min_resid)<rmin ){
+ test_fit.add_hit_to_fit(sh1_min);
+ sh1_min.SetUsedOnHel(npath); //try this
+ if(sh1_min.IsForegroundHit())complete_trk+=1;
+ }
+//layer 2
+ for (ListIterator ii = l2shlist.listIterator();ii.hasNext();){
+ SODHit sh2=(SODHit)ii.next();
+ if ((sh2.GetUsedOnHel()==0)&&(sh2.calc_resid(test_fit))){
+ double resid= sh2.get_resid() ;
+ if ( Math.abs(resid)<Math.abs(l2_min_resid) ){l2_min_resid=resid; sh2_min=sh2;}
+ }
+ }
+ if ( Math.abs(l2_min_resid)<rmin ){
+ test_fit.add_hit_to_fit(sh2_min);
+ if(sh2_min.IsForegroundHit())complete_trk+=2;
+ }
+//layer 3
+ for (ListIterator ii = l3shlist.listIterator();ii.hasNext();){
+ SODHit sh3=(SODHit)ii.next();
+ if ((sh3.GetUsedOnHel()==0)&&(sh3.calc_resid(test_fit))){
+ double resid= sh3.get_resid() ;
+ if ( Math.abs(resid)<Math.abs(l3_min_resid) ){l3_min_resid=resid; sh3_min=sh3;}
+ }
+ }
+ if ( Math.abs(l3_min_resid)<rmin ){
+ test_fit.add_hit_to_fit(sh3_min);
+ if(sh3_min.IsForegroundHit())complete_trk+=4;
+ }
+//layer 4
+ for (ListIterator ii = l4shlist.listIterator();ii.hasNext();){
+ SODHit sh4=(SODHit)ii.next();
+ if ((sh4.GetUsedOnHel()==0)&&(sh4.calc_resid(test_fit))){
+ double resid= sh4.get_resid() ;
+ if ( Math.abs(resid)<Math.abs(l4_min_resid) ){l4_min_resid=resid; sh4_min=sh4;}
+ }
+ }
+ if ( Math.abs(l4_min_resid)<rmin ){
+ test_fit.add_hit_to_fit(sh4_min);
+ if(sh4_min.IsForegroundHit())complete_trk+=8;
+ }
+//layer 5
+ for (ListIterator ii = l5shlist.listIterator();ii.hasNext();){
+ SODHit sh5=(SODHit)ii.next();
+ if ((sh5.GetUsedOnHel()==0)&&(sh5.calc_resid(test_fit))){
+ double resid= sh5.get_resid() ;
+ if ( Math.abs(resid)<Math.abs(l5_min_resid) ){l5_min_resid=resid; sh5_min=sh5;}
+ }
+ }
+ if ( Math.abs(l5_min_resid)<rmin ){
+ test_fit.add_hit_to_fit(sh5_min);
+//was a bug if(sh3_min.IsForegroundHit())complete_trk+=16;
+ if(sh5_min.IsForegroundHit())complete_trk+=16;
+ }
+// if trk OK, add it to trk list
+// if (1==npath)System.out.print("SODL1Path11111");
+// if (2==npath)System.out.print("SODL1Path21111");
+// if (3==npath)System.out.print("SODL1Path31111");
+// System.out.print(" ");System.out.print(l1_min_resid);//usually
+// if (1==npath){System.out.print(" ");System.out.print(l1_correct_rank);}
+// System.out.print(" ");System.out.print(l2_min_resid);//usually
+// System.out.print(" ");System.out.print(l3_min_resid);//usually
+// System.out.print(" ");System.out.print(l4_min_resid);//usually
+// System.out.print(" ");System.out.print(l5_min_resid);//usually
+// System.out.print(" ");System.out.print(test_fit.Chisq());
+// System.out.print(" ");System.out.print(test_fit.Nhits());
+// System.out.print(" ");System.out.print(complete_trk);
+// System.out.println(" ");
+ }//end ctor from 5 hit lists
+ public double Rcs(){
+ double rcs=1000000.0; double xhits=(double)test_fit.Nhits();
+ if (xhits>3.0){rcs=test_fit.Chisq()/(xhits-3.0);}
+ return rcs;
+ }
+ public double Chisq(){return test_fit.Chisq();}
+ public int Nhits(){return test_fit.Nhits();}
+ public double Pt(){return test_fit.Pt();}
+ public int Complete_trk(){return complete_trk;}
+ public SODFittedCir out_fit(){return test_fit;}
+ private SODFittedCir test_fit;
+ private int complete_trk;
+ private static double pi=3.14159265359;
+ private static double twopi=6.283185307;
+}
lcsim/src/org/lcsim/contrib/SODTracker
diff -N SODNoiseMaker.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SODNoiseMaker.java 10 Mar 2006 21:21:34 -0000 1.1
@@ -0,0 +1,34 @@
+package org.lcsim.contrib.SODTracker;
+
+import java.util.*;
+import java.io.*;
+public class SODNoiseMaker{
+ public SODNoiseMaker(java.util.LinkedList l1list, java.util.LinkedList l2list,
+ java.util.LinkedList l3list, java.util.LinkedList l4list, java.util.LinkedList l5list) {
+ double rl1=20.05,rl2=46.30,rl3=72.55,rl4=98.80,rl5=125.05;
+ //split sod
+ int nl1=106,nl2=91,nl3=92,nl4=98,nl5=109;
+ //tiled sod int nl1=36,nl2=5,nl3=1,nl4=0,nl5=1;
+ boolean success;
+ for(int i=0; i<nl1; i++) {
+ double myphi = Math.random()*Math.PI*2.0;
+ SODHit sh1=new SODHit(); sh1.make_hit(rl1,myphi); success=l1list.add(sh1);
+ }
+ for(int i=0; i<nl2; i++) {
+ double myphi = Math.random()*Math.PI*2.0;
+ SODHit sh2=new SODHit(); sh2.make_hit(rl2,myphi); success=l2list.add(sh2);
+ }
+ for(int i=0; i<nl3; i++) {
+ double myphi = Math.random()*Math.PI*2.0;
+ SODHit sh3=new SODHit(); sh3.make_hit(rl3,myphi); success=l3list.add(sh3);
+ }
+ for(int i=0; i<nl4; i++) {
+ double myphi = Math.random()*Math.PI*2.0;
+ SODHit sh4=new SODHit(); sh4.make_hit(rl4,myphi); success=l4list.add(sh4);
+ }
+ for(int i=0; i<nl5; i++) {
+ double myphi = Math.random()*Math.PI*2.0;
+ SODHit sh5=new SODHit(); sh5.make_hit(rl5,myphi); success=l5list.add(sh5);
+ }
+ }
+}
lcsim/src/org/lcsim/contrib/SODTracker
diff -N SODReadGeom.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SODReadGeom.java 10 Mar 2006 21:21:35 -0000 1.1
@@ -0,0 +1,106 @@
+package org.lcsim.contrib.SODTracker;
+
+import java.util.*;
+import java.io.*;
+public class SODReadGeom
+{
+ public SODReadGeom(BufferedReader bread)
+ {
+ System.out.println("read in geometry with SODReadGeom ");
+ boolean success;
+ String input;
+ _radvtxlay = new double[10];
+ _radsodlay = new double[10];
+ _radsidlay = new double[20];
+ _resvtxlay = new double[10];
+ _ressodlay = new double[10];
+ _ressidlay = new double[20];
+ try{input = bread.readLine();
+ int[] values = new int[13];
+ int j=0;
+ StringTokenizer st = new StringTokenizer(input);
+ while(st.hasMoreElements()){
+ if (j<13){values[j] = Float.valueOf(st.nextToken()).intValue();}
+ j++;
+ }
+ if (1==j){
+ _nvtxlay=values[0];
+ System.out.print(values[0]);
+ System.out.println(" vtx layers in SODReadGeom ");
+ }else{System.out.println("BIG FUCKUP IN SODReadGeom");}
+
+ int nl1=values[0];
+ if (nl1>0){
+ input = bread.readLine(); st = new StringTokenizer(input); j=0;
+ while(st.hasMoreElements()){
+ double myphi=Double.valueOf(st.nextToken()).doubleValue();
+ _radvtxlay[j]=myphi;
+ j++;
+ }
+ if (nl1==j){// System.out.println(" ");
+ }else{System.out.println("BIG FUCKUP IN SODReadGeom vtx rad");}
+ }
+
+ int nl2=values[0]; if (nl2>0){
+ input = bread.readLine(); st = new StringTokenizer(input); j=0;
+ while(st.hasMoreElements()){
+ double myphi=Double.valueOf(st.nextToken()).doubleValue();
+ _resvtxlay[j]=myphi;
+ j++;
+ }
+ if (nl2==j){// System.out.println(" ");
+ }else{System.out.println("BIG FUCKUP IN SODReadGeom vtx res");}
+ }
+
+ input = bread.readLine(); st = new StringTokenizer(input); j=0;
+ while(st.hasMoreElements()){
+ if (j<13){values[j] = Float.valueOf(st.nextToken()).intValue();}
+ j++;
+ }
+ if (1==j){
+ _nsodlay=values[0];
+ System.out.print(values[0]);
+ System.out.println(" sod layers in SODReadGeom ");
+ }else{System.out.println("BIG FUCKUP IN SODReadGeom");}
+
+ nl1=values[0]; if (nl1>0){
+ input = bread.readLine(); st = new StringTokenizer(input); j=0;
+ while(st.hasMoreElements()){
+ double myphi=Double.valueOf(st.nextToken()).doubleValue();
+ _radsodlay[j]=myphi;
+ int k=_nvtxlay+j;
+ _radsidlay[k]=myphi;
+ j++;
+ }
+ if (nl1==j){// System.out.println(" ");
+ }else{System.out.println("BIG FUCKUP IN SODReadGeom vtx rad");}
+ }
+
+ nl2=values[0]; if (nl2>0){
+ input = bread.readLine(); st = new StringTokenizer(input); j=0;
+ while(st.hasMoreElements()){
+ double myphi=Double.valueOf(st.nextToken()).doubleValue();
+ _ressodlay[j]=myphi;
+ int k=_nvtxlay+j;
+ _ressidlay[k]=myphi;
+ j++;
+ }
+ if (nl2==j){// System.out.println(" ");
+ }else{System.out.println("BIG FUCKUP IN SODReadGeom sod res");}
+ }
+
+ }catch(Exception e){System.out.println("blow me!! readLine failed");e.printStackTrace();}
+ }
+ public double rad_at_vtx_layer(int i){return _radvtxlay[i-1];}
+ public double rad_at_sod_layer(int i){return _radsodlay[i-1];}
+ public double res_at_vtx_layer(int i){return _resvtxlay[i-1];}
+ public double res_at_sod_layer(int i){return _ressodlay[i-1];}
+ private int _nvtxlay;
+ private int _nsodlay;
+ private double[] _radvtxlay;
+ private double[] _radsodlay;
+ private double[] _radsidlay;
+ private double[] _resvtxlay;
+ private double[] _ressodlay;
+ private double[] _ressidlay;
+}
lcsim/src/org/lcsim/contrib/SODTracker
diff -N SODTrack.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SODTrack.java 10 Mar 2006 21:21:35 -0000 1.1
@@ -0,0 +1,83 @@
+package org.lcsim.contrib.SODTracker;
+
+/* SODTrack.java
+
+*/
+
+import org.lcsim.event.Track;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.event.MCParticle;
+
+import java.util.ArrayList;
+import java.util.List;
+
+
+public class SODTrack implements Track
+{
+
+ List<SimTrackerHit> list = new ArrayList();
+ double p = 0., px = 0., py = 0., pz = 0.;
+ double[] momentum;
+ int charge = 0;
+
+ double d0 = 0., phi0 = 0., omega = 0., z0 = 0., s = 0.;
+ double[] parameters = new double[] {0.,0.,0.,0.,0.};
+
+ SODTrack(int ch, double pX, double pY, double pZ, double dzero, double phizero, double omeg, double zzero, double spar)
+ {
+ charge = ch;
+ px = pX;
+ py = pY;
+ pz = pZ;
+ p = Math.sqrt(px*px + py*py + pz*pz);
+ momentum = new double[] {px, py, pz};
+ d0 = dzero;
+ phi0 = phizero;
+ // CAUTION: different definition of omega in SOD and lcsim!
+ omega = -0.1*omeg;
+ z0 = zzero;
+ s = spar;
+ parameters = new double[] {d0, phi0, omega,z0,s};
+ }
+
+ private MCParticle mcParticle;
+ public void setMCParticle(MCParticle particle) { this.mcParticle = particle; };
+ public MCParticle getMCParticle() { return mcParticle; }
+
+ public void addHit(SimTrackerHit hit) { list.add(hit); }
+ public int getNPoints() { return list.size(); }
+ public List<SimTrackerHit> getHits() { return list; }
+
+ public int getCharge() { return charge; }
+ public double getPX() { return px; }
+ public double getPY() { return py; }
+ public double getPT() { return Math.sqrt(px*px+py*py); }
+ public double getPtot() { return p; }
+ public double getPZ() { return pz; }
+ public double[] getMomentum() { return momentum; }
+
+ public double[] getReferencePoint() { double[] point = new double[] {0.,0.,0.}; return point; }
+ public double getReferencePointX() { return 0.; }
+ public double getReferencePointY() { return 0.; }
+ public double getReferencePointZ() { return 0.; }
+ public boolean isReferencePointPCA() { return false; }
+
+ public boolean fitSuccess() { return true; }
+ public double getTrackParameter(int i) { return parameters[i]; }
+ public double[] getTrackParameters() { return parameters; }
+ public double getErrorMatrixElement(int i, int j) { return 0.; }
+ public double[][] getErrorMatrix() { double[][] matrix = new double[][] {{0.},{0.}}; return matrix; }
+ public double getChi2() { return 1.; }
+ public int getNDF() { return 1; }
+
+ public double getdEdx() { return 0.; }
+ public double getdEdxError() { return 0.; }
+ public double getRadiusOfInnermostHit() { return 0.; }
+ public int[] getSubdetectorHitNumbers() { int[] hitNumbers = new int[]{0}; return hitNumbers; }
+ public List<Track> getTracks() { List<Track> list = new ArrayList<Track>(); return list; }
+ public List<TrackerHit> getTrackerHits() { List<TrackerHit> list = new ArrayList<TrackerHit>(); return list; }
+ public int getType() { return 0; }
+
+}
+
lcsim/src/org/lcsim/contrib/SODTracker
diff -N SODTrackFinder.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SODTrackFinder.java 10 Mar 2006 21:21:35 -0000 1.1
@@ -0,0 +1,266 @@
+package org.lcsim.contrib.SODTracker;
+
+import java.util.*;
+import org.lcsim.event.MCParticle;
+public class SODTrackFinder
+{
+ public SODTrackFinder()
+ {
+ _trklist = new java.util.LinkedList();
+ }
+
+ public SODTrackFinder(java.util.LinkedList l1shlist,java.util.LinkedList l2shlist,
+ java.util.LinkedList l3shlist,java.util.LinkedList l4shlist,
+ java.util.LinkedList l5shlist,boolean useMCTruth) {
+ if(useMCTruth) {
+ int prtflg=0;
+ _trklist = new java.util.LinkedList(); int ncombo=0;
+ if ((l1shlist.size()>0)&&(l2shlist.size()>0)&&(l3shlist.size()>0)&&
+ (l4shlist.size()>0)&&(l5shlist.size()>0)){
+ for (ListIterator j1 = l1shlist.listIterator();j1.hasNext();){
+ SODHit sh1=(SODHit)j1.next();
+ MCParticle mcp1 = sh1.getMCParticle();
+ for (ListIterator j2 = l2shlist.listIterator();j2.hasNext();){
+ SODHit sh2=(SODHit)j2.next();
+ MCParticle mcp2 = sh2.getMCParticle();
+ if(mcp1==mcp2) {
+ for (ListIterator j3 = l3shlist.listIterator();j3.hasNext();){
+ SODHit sh3=(SODHit)j3.next();
+ MCParticle mcp3 = sh3.getMCParticle();
+ if(mcp2==mcp3) {
+ for (ListIterator j4 = l4shlist.listIterator();j4.hasNext();){
+ SODHit sh4=(SODHit)j4.next();
+ MCParticle mcp4 = sh4.getMCParticle();
+ if(mcp3==mcp4) {
+ for (ListIterator j5 = l5shlist.listIterator();j5.hasNext();){
+ SODHit sh5=(SODHit)j5.next();
+ MCParticle mcp5 = sh5.getMCParticle();
+ if(mcp4==mcp5) {
+ ncombo++;
+ // TrkParams test_helix = new TrkParams();
+ SODHel test_helix = new SODHel();
+ test_helix = make_test_helix(sh2,sh3,sh4);
+
+ java.util.LinkedList shitlist = new java.util.LinkedList();
+ shitlist.add(sh1);shitlist.add(sh2);shitlist.add(sh3);
+ shitlist.add(sh4);shitlist.add(sh5);
+ SODFittedCir fit_helix = new SODFittedCir();
+ fit_helix.make_fitted_cir(shitlist,test_helix);
+ boolean sflag=false;
+ if (0==fit_helix.Fail())sflag=_trklist.add(fit_helix);
+ //
+ if (0==fit_helix.Fail()){
+ sh1.SetUsedOnHel(1);sh2.SetUsedOnHel(1);sh3.SetUsedOnHel(1);
+ sh4.SetUsedOnHel(1);sh5.SetUsedOnHel(1);
+ }
+ } // if(mcp4==mcp5) {
+ } // for l5
+ } // if(mcp3==mcp4) {
+ } // for l4
+ } // if(mcp2==mcp3) {
+ } // for l3
+ }// if(mcp1==mcp2) {
+ } // for l2
+ } // for l1
+ } // if lists not empty
+
+ } else { // if(useMCTruth) {
+
+ int prtflg=0;
+ _trklist = new java.util.LinkedList(); int ncombo=0;
+ if ((l1shlist.size()>0)&&(l2shlist.size()>0)&&(l3shlist.size()>0)&&
+ (l4shlist.size()>0)&&(l5shlist.size()>0)){
+ for (ListIterator j1 = l2shlist.listIterator();j1.hasNext();){
+ SODHit sh2=(SODHit)j1.next();
+ for (ListIterator j2 = l3shlist.listIterator();j2.hasNext();){
+ SODHit sh3=(SODHit)j2.next();
+ for (ListIterator j3 = l4shlist.listIterator();j3.hasNext();){
+ SODHit sh4=(SODHit)j3.next();
+ ncombo++;
+ // make test helix from hits in layers 2 3 4, if not already used
+ if ((sh2.GetUsedOnHel()==0)&&(sh3.GetUsedOnHel()==0)&&
+ (sh4.GetUsedOnHel()==0)){
+ SODHel test_helix = new SODHel();
+ test_helix = make_test_helix(sh2,sh3,sh4);
+ // see if hits in layers 1 5 are with tolerance
+ double l1_min_resid=1000.0,l5_min_resid=1000.0;
+ SODHit sh1_min=(SODHit)l1shlist.getFirst();
+ SODHit sh5_min=(SODHit)l5shlist.getFirst();
+ for (ListIterator ii = l1shlist.listIterator();ii.hasNext();){
+ SODHit sh1=(SODHit)ii.next();
+ if ((sh1.GetUsedOnHel()==0)&&(sh1.calc_resid(test_helix))){
+ double resid= sh1.get_resid() ;
+ if ( Math.abs(resid)<Math.abs(l1_min_resid) ){l1_min_resid=resid; sh1_min=sh1;}
+ }
+ }
+ for (ListIterator ii = l5shlist.listIterator();ii.hasNext();){
+ SODHit sh5=(SODHit)ii.next();
+ if ((sh5.GetUsedOnHel()==0)&&(sh5.calc_resid(test_helix))){
+ double resid= sh5.get_resid() ;
+ if ( Math.abs(resid)<Math.abs(l5_min_resid) ){l5_min_resid=resid; sh5_min=sh5;}
+ }
+ }
+ // if trk OK, add it to trk list
+ // double pt=0.015/Math.abs( test_helix.getParam(2) );
+ double pt=test_helix.Pt();
+ int ndof=2;
+ double search1=search_fcn(pt,1); double search5=search_fcn(pt,5);
+ if ((Math.abs(l1_min_resid)<5.0*search1)&&(Math.abs(l5_min_resid)<5.0*search5)){
+ double chisq=(l1_min_resid/search1)*(l1_min_resid/search1);
+ chisq+=(l5_min_resid/search5)*(l5_min_resid/search5);
+ // boolean sflag=_trklist.add(test_helix);
+ // no, do a fit now
+ java.util.LinkedList shitlist = new java.util.LinkedList();
+ shitlist.add(sh1_min);shitlist.add(sh2);shitlist.add(sh3);shitlist.add(sh4);shitlist.add(sh5_min);
+ // test_helix.printHelix();
+ double smear_d0=0.0;
+ if (Math.abs(l1_min_resid)<Math.abs(l5_min_resid)){
+ smear_d0=l5_min_resid/2.0;
+ }else{
+ smear_d0=l1_min_resid/2.0;
+ }
+ SODFittedCir fit_helix = new SODFittedCir();
+ fit_helix.make_fitted_cir(shitlist,test_helix);
+ //smear d0 to see if that help inverting
+ // SODHel smear_hel = new SODHel();
+ // smear_hel.make_hel(test_helix.D0()+smear_d0,test_helix.Phi0(),
+ // test_helix.Omega(),0.0,0.0);
+ // fit_helix.make_fitted_cir(shitlist,smear_hel);
+ //end smearing d0
+ // int sflg=fit_helix.DoFit();
+ // int sflg=fit_helix.IterateFit();
+ // move IterateFit into make_fitted_cir
+ boolean sflag=false;
+ if (0==fit_helix.Fail())sflag=_trklist.add(fit_helix);
+ //see if move ref pt helps
+ // test_helix.change_ref(sh1_min.x(),sh1_min.y());
+ // SODFittedCir move_helix = new SODFittedCir(); move_helix.make_fitted_cir(shitlist,test_helix);
+ // sflg=move_helix.DoFit();
+ prtflg=0;
+ if (prtflg!=0){
+ System.out.print("chisq Chisq ");System.out.print(chisq);
+ System.out.print(" ");System.out.println( fit_helix.Chisq() );
+ }//end prtflg
+
+ if (0==fit_helix.Fail()){
+ sh1_min.SetUsedOnHel(1);sh2.SetUsedOnHel(1);sh3.SetUsedOnHel(1);
+ sh4.SetUsedOnHel(1);sh5_min.SetUsedOnHel(1);
+ }
+ prtflg=0;
+ if (prtflg!=0){
+ System.out.print("sodtftup "); //System.out.print(ncombo);
+ System.out.print(" ");System.out.print(test_helix.D0());
+ System.out.print(" ");System.out.print(test_helix.Phi0());
+ System.out.print(" ");System.out.print(test_helix.Omega());
+ System.out.print(" ");System.out.print(test_helix.Z0());
+ System.out.print(" ");System.out.print(test_helix.Tanl());
+ System.out.print(" ");System.out.print(pt);
+ System.out.print(" ");System.out.print(ndof);
+ System.out.print(" ");System.out.println(chisq);
+ }//done printing
+ prtflg=0;
+ if (prtflg!=0){
+ System.out.print("restup ");System.out.print(l1_min_resid);
+ sh2.calc_resid(fit_helix);
+ System.out.print(" ");System.out.print(sh2.get_resid());
+ sh3.calc_resid(fit_helix);
+ System.out.print(" ");System.out.print(sh3.get_resid());
+ sh4.calc_resid(fit_helix);
+ System.out.print(" ");System.out.print(sh4.get_resid());
+ System.out.print(" ");System.out.println(l5_min_resid);
+ }//done printing
+ }//end it's a track
+ }//seed hits not already used
+ }//loop over l4 hits
+ }//loop over l3 hits
+ }//loop over l2 hits
+ }//(unindented) all list have hits
+ } // if(useMCTruth) { } // contructor
+
+
+// public TrkParams make_test_helix(SODHit sh1, SODHit sh2, SODHit sh3)
+ public SODHel make_test_helix(SODHit sh1, SODHit sh2, SODHit sh3)
+ {
+ double x1=sh1.x(),x2=sh2.x(),x3=sh3.x(),y1=sh1.y(),y2=sh2.y(),y3=sh3.y();
+ double d12x=x1-x2,d12y=y1-y2,d13x=x1-x3,d13y=y1-y3,d23x=x2-x3,d23y=y2-y3;
+ double d12=Math.sqrt(d12x*d12x+d12y*d12y);
+ double d13=Math.sqrt(d13x*d13x+d13y*d13y);
+ double d23=Math.sqrt(d23x*d23x+d23y*d23y);
+ double s=0.5*(d12+d13+d23);
+ double r_denom=4.0*Math.sqrt(s*(s-d12)*(s-d13)*(s-d23));
+ double aom=r_denom/(d12*d13*d23); double r=1/aom;
+ double hphi12=Math.asin(d12/(2.0*r));
+ double phi12=2.0*hphi12,cp12=Math.cos(phi12),sp12=Math.sin(phi12);
+ double xm=0.5*(x1+x2),ym=0.5*(y1+y2),dm=r*Math.cos(hphi12);
+ double cal=d12y/d12,sal=d12x/d12;
+ double rm=Math.sqrt(xm*xm+ym*ym);
+ double xc1=xm+dm*cal;
+ double yc1=ym-dm*sal;
+ double xc2=xm-dm*cal;
+ double yc2=ym+dm*sal;
+ double rc1=Math.sqrt(xc1*xc1+yc1*yc1),d01=rc1-r;
+ double rc2=Math.sqrt(xc2*xc2+yc2*yc2),d02=rc2-r;
+
+// double e11=Math.sqrt((x1-xc1)*(x1-xc1)+(y1-yc1)*(y1-yc1))-r;
+// double e21=Math.sqrt((x2-xc1)*(x2-xc1)+(y2-yc1)*(y2-yc1))-r;
+ double e31=Math.sqrt((x3-xc1)*(x3-xc1)+(y3-yc1)*(y3-yc1))-r;
+// double e12=Math.sqrt((x1-xc2)*(x1-xc2)+(y1-yc2)*(y1-yc2))-r;
+// double e22=Math.sqrt((x2-xc2)*(x2-xc2)+(y2-yc2)*(y2-yc2))-r;
+ double e32=Math.sqrt((x3-xc2)*(x3-xc2)+(y3-yc2)*(y3-yc2))-r;
+ if (e31<0.0)e31=-1.0*e31;if (e32<0.0)e32=-1.0*e32; if (d01<0.0)d01=-1.0*d01;if (d02<0.0)d02=-1.0*d02;
+ double xc=0.0,yc=0.0,rc=0.0,d0=0.0,phi0=0.0,omega=0.0;
+ if (d01<d02){
+ if (e31<e32){//System.out.println("d01 e31"); //right
+ xc=xc1;yc=yc1;rc=rc1;
+// cosmetic d0=d01; phi0=-1.0*Math.atan2(xc,yc); omega=1.0/r;if (r>rc)d0=-1.0*d0;
+ d0=d01; phi0=-Math.atan2(xc,yc); omega=1.0/r; if (r>rc)d0=-d0;
+ }else{//System.out.println("d01 e32"); //fixed;was bad wrong
+ xc=xc2;yc=yc2;rc=rc2;
+// d0=d01; phi0=-1.0*Math.atan2(xc,yc); omega=1.0/r;if (r>rc)d0=-1.0*d0;
+ d0=-d02; phi0=pi-Math.atan2(xc,yc); omega=-1.0/r;if (r>rc)d0=-1.0*d0;
+ }
+ }else{
+ if (e31<e32){//System.out.println("d02 e31"); //fixed;was kinda wrong
+// xc=xc2;yc=yc2;rc=rc2;
+ xc=xc1;yc=yc1;rc=rc1;
+// d0=d02; phi0=pi-Math.atan2(xc,yc); omega=-1.0/r;if (r<rc)d0=-1.0*d0;
+ d0=-d01; phi0=-Math.atan2(xc,yc); omega=1.0/r;if (r<rc)d0=-1.0*d0;
+ }else{//System.out.println("d02 e32"); //fixed;was bad wrong
+// xc=xc1;yc=yc1;rc=rc1; xc=xc2;yc=yc2;rc=rc2;
+// d0=d02; phi0=pi-Math.atan2(xc,yc); omega=-1.0/r;if (r<rc)d0=-1.0*d0;
+ d0=d02; phi0=pi-Math.atan2(xc,yc); omega=-1.0/r;if (r<rc)d0=-1.0*d0;
+ }
+ }
+ int prtflg=0;
+ if (prtflg!=0){
+ System.out.print("test_helix ");System.out.print(d0);
+ System.out.print(" ");System.out.print(phi0);
+ System.out.print(" ");System.out.print(omega);
+ System.out.print(" ");System.out.print(d01);
+ System.out.print(" ");System.out.print(d02);
+ System.out.print(" ");System.out.print(r);
+ System.out.print(" ");System.out.println(rc);
+// System.out.print("test_helix ");System.out.print(e11);
+// System.out.print(" ");System.out.print(e21);
+ System.out.print(" ");System.out.print(e31);
+// System.out.print(" ");System.out.print(e12);
+// System.out.print(" ");System.out.print(e22);
+ System.out.print(" ");System.out.println(e32);
+ }
+// TrkParams test_helix = new TrkParams(d0,phi0,omega,0.0,0.0);
+ SODHel test_helix=new SODHel();test_helix.make_hel(d0,phi0,omega,0.0,0.0);
+ return test_helix;
+ }
+ public double search_fcn(double pt, int lay){
+ double search = -0.00040625*pt+0.0208125;
+ if ((pt<2.0)&&(lay==1))search=-0.01*pt+0.04;
+ if ((pt<2.0)&&(lay==5))search=-0.18*pt+0.38;
+ if ((search<0.0130)&&(lay==1))search=0.0130;
+ if ((search<0.0141)&&(lay==5))search=0.0141;
+ return search;
+ }
+ public java.util.LinkedList trklist(){return _trklist;}
+ private java.util.LinkedList _trklist;
+ private static double pi=3.14159265359;
+ private static double twopi=6.283185307;
+}
lcsim/src/org/lcsim/contrib/SODTracker
diff -N SODTrackFinderDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SODTrackFinderDriver.java 10 Mar 2006 21:21:35 -0000 1.1
@@ -0,0 +1,203 @@
+package org.lcsim.contrib.SODTracker;
+
+import java.util.*;
+import java.io.*;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.event.EventHeader.LCMetaData;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.MCParticle;
+import org.lcsim.geometry.TrackerIDDecoder;
+
+public class SODTrackFinderDriver extends Driver
+{
+ public SODTrackFinderDriver()
+ {
+ boolean dooutfile=false; if (dooutfile){
+ try{
+ FileOutputStream fos = new FileOutputStream("output.dat");
+ PrintStream ps = new PrintStream(fos);
+ System.setOut(ps);
+ }
+ catch(FileNotFoundException x){
+ System.out.println("blow me!! output.dat doesn't exist");
+ x.printStackTrace();
+ }
+ }//end dooutfile
+ generator = new java.util.Random(10101);
+ }
+
+ protected void process(EventHeader event)
+ {
+ int ntrk=0,prtflg=0,clnflg=0;
+ java.util.LinkedList l1shlist = new java.util.LinkedList();
+ java.util.LinkedList l2shlist = new java.util.LinkedList();
+ java.util.LinkedList l3shlist = new java.util.LinkedList();
+ java.util.LinkedList l4shlist = new java.util.LinkedList();
+ java.util.LinkedList l5shlist = new java.util.LinkedList();
+ java.util.LinkedList l1vhlist = new java.util.LinkedList();
+ java.util.LinkedList l2vhlist = new java.util.LinkedList();
+ java.util.LinkedList l3vhlist = new java.util.LinkedList();
+ java.util.LinkedList l4vhlist = new java.util.LinkedList();
+ java.util.LinkedList l5vhlist = new java.util.LinkedList();
+ int collection=0;
+ List<List<SimTrackerHit>> simTrackerHitCollections = event.get(SimTrackerHit.class);
+ for ( List<SimTrackerHit> simTrackerHits : simTrackerHitCollections ) {
+ collection++;
+ LCMetaData meta = event.getMetaData(simTrackerHits);
+ TrackerIDDecoder decoder = (TrackerIDDecoder) meta.getIDDecoder();
+ for (SimTrackerHit trackerHit : simTrackerHits) {
+ decoder.setID(trackerHit.getCellID() );
+ int layer = decoder.getLayer();
+ double[] pos = trackerHit.getPoint();
+ // get the MC particle corresponding to this hit
+ MCParticle mcp = trackerHit.getMCParticle();
+ double rgauss = generator.nextGaussian();
+ SODHit sh=new SODHit();
+ if (collection==2){
+ // the second collection corresponds to the barrel vertex tracker
+ sh.make_hit(layer,pos[0],pos[1],rgauss,trackerHit,mcp);
+ } else {
+ sh.make_hit(layer,pos[0],pos[1],rgauss,trackerHit);
+ }
+ double res=-(0.1*pos[0]-sh.x())*Math.sin(sh.phi())
+ +(0.1*pos[1]-sh.y())*Math.cos(sh.phi());
+ if (collection==2){
+ if (layer==0){boolean success=l1vhlist.add(sh);}
+ if (layer==1){boolean success=l2vhlist.add(sh);}
+ if (layer==2){boolean success=l3vhlist.add(sh);}
+ if (layer==3){boolean success=l4vhlist.add(sh);}
+ if (layer==4){boolean success=l5vhlist.add(sh);}
+ }
+ if (collection==3){
+ if (layer==0){boolean success=l1shlist.add(sh);}
+ if (layer==1){boolean success=l2shlist.add(sh);}
+ if (layer==2){boolean success=l3shlist.add(sh);}
+ if (layer==3){boolean success=l4shlist.add(sh);}
+ if (layer==4){boolean success=l5shlist.add(sh);}
+ }
+ clnflg=1;
+ // if ((l1vhlist.size()!=1)||(l2vhlist.size()!=1)||(l3vhlist.size()!=1)||(l4vhlist.size()!=1)||(l5vhlist.size()!=1))clnflg=0;
+ }
+ }
+ // SODNoiseMaker nrdr = new SODNoiseMaker(l1shlist,l2shlist,l3shlist,l4shlist,l5shlist);
+ //
+ //c no SOD tracking this version; find probe trk (and maybe bgtrk) in vtx and SODHitAdder
+ //
+ ntrk=0;
+ java.util.LinkedList mytrkl = new java.util.LinkedList();
+ java.util.LinkedList bgtrkl = new java.util.LinkedList();
+ java.util.LinkedList shatrkl = new java.util.LinkedList();
+
+ // if (clnflg==1){
+ boolean useMCTruth = true;
+ SODTrackFinder sodtf = new SODTrackFinder(l1vhlist,l2vhlist,l3vhlist,
+ l4vhlist,l5vhlist,useMCTruth);
+
+ // get list of tracks found in vertexer by SODTrackFinder
+ mytrkl = sodtf.trklist();
+
+// System.out.println(event.getEventNumber());
+
+ //
+ //c now run SODHitAdder, which really is the outer trker trkfinder
+ //
+ ntrk=0;
+ List<SODTrack> listTrk = new ArrayList<SODTrack>();
+ //c for now only work if hits are present in all 5 layers
+ if ((l1shlist.size()!=0)&&(l2shlist.size()!=0)&&(l3shlist.size()!=0)&&
+ (l4shlist.size()!=0)&&(l5shlist.size()!=0)) {
+ for (ListIterator jj = mytrkl.listIterator();jj.hasNext();){
+ SODFittedCir toast_helix=(SODFittedCir)jj.next();
+ SODHitAdder sodha = new SODHitAdder(toast_helix,l1shlist,l2shlist,l3shlist,l4shlist,l5shlist);
+ free_hit_lists(l1shlist,l2shlist,l3shlist,l4shlist,l5shlist);
+ SODFittedCir out_helix= new SODFittedCir();
+ out_helix=sodha.get_new_fit();
+ boolean success=shatrkl.add(out_helix);
+ SODHel test_helix=(SODHel)out_helix;
+ ntrk++;
+ // create new Tracks
+ double l = 0;
+ SODTrack trk = new SODTrack(test_helix.Q(), test_helix.Px(l),
+ test_helix.Py(l), test_helix.Pz(),
+ test_helix.D0(), test_helix.Phi0(),
+ test_helix.Omega(), test_helix.Z0(),
+ test_helix.Tanl());
+ java.util.LinkedList hitlist = out_helix.hitlist();
+ for ( ListIterator hitItr = hitlist.listIterator() ; hitItr.hasNext();){
+ SODHit hit = (SODHit)hitItr.next();
+ SimTrackerHit h = hit.getHit();
+ trk.addHit(h);
+ }
+ listTrk.add(trk);
+ }//end loop over found trks
+
+ // write list of found SODTracks in the event
+// if(ntrk>0) { event.put("SODTracks", listTrk, SODTrack.class, 0); }
+
+ }//end lists not empty
+ event.put("SODTracks", listTrk, SODTrack.class, 0);
+
+ //
+ //print out evt info
+ prtflg=0;
+ if (prtflg!=0){
+ System.out.print("evttup");
+ System.out.print(" ");System.out.print(event.getRunNumber());
+ System.out.print(" ");System.out.print(event.getEventNumber());
+ System.out.print(" ");System.out.print(ntrk);
+ System.out.print(" ");System.out.print(l1vhlist.size());
+ System.out.print(" ");System.out.print(l2vhlist.size());
+ System.out.print(" ");System.out.print(l3vhlist.size());
+ System.out.print(" ");System.out.print(l4vhlist.size());
+ System.out.print(" ");System.out.print(l5vhlist.size());
+ System.out.print(" ");System.out.print(l1shlist.size());
+ System.out.print(" ");System.out.print(l2shlist.size());
+ System.out.print(" ");System.out.print(l3shlist.size());
+ System.out.print(" ");System.out.print(l4shlist.size());
+ System.out.print(" ");System.out.print(l5shlist.size());
+ System.out.println(" ");
+ }
+ }//end process
+//
+//
+//
+ public void free_hit_lists(java.util.LinkedList l1list, java.util.LinkedList l2list,
+ java.util.LinkedList l3list, java.util.LinkedList l4list,
+ java.util.LinkedList l5list)
+ {
+ if ((l1list.size()>0)){
+ for (ListIterator ii = l1list.listIterator();ii.hasNext();){
+ SODHit sh=(SODHit)ii.next(); sh.SetUsedOnHel(0);
+ }
+ }
+ if ((l2list.size()>0)){
+ for (ListIterator ii = l2list.listIterator();ii.hasNext();){
+ SODHit sh=(SODHit)ii.next(); sh.SetUsedOnHel(0);
+ }
+ }
+ if ((l3list.size()>0)){
+ for (ListIterator ii = l3list.listIterator();ii.hasNext();){
+ SODHit sh=(SODHit)ii.next(); sh.SetUsedOnHel(0);
+ }
+ }
+ if ((l4list.size()>0)){
+ for (ListIterator ii = l4list.listIterator();ii.hasNext();){
+ SODHit sh=(SODHit)ii.next(); sh.SetUsedOnHel(0);
+ }
+ }
+ if ((l5list.size()>0)){
+ for (ListIterator ii = l5list.listIterator();ii.hasNext();){
+ SODHit sh=(SODHit)ii.next(); sh.SetUsedOnHel(0);
+ }
+ }
+ }
+
+ FileInputStream fis6;
+ BufferedReader gread;
+ SODReadGeom geom;
+ private static double pi=3.14159265359;
+ private static double twopi=6.283185307;
+ private int cntbg=0;
+ java.util.Random generator;
+}
lcsim/src/org/lcsim/contrib/SODTracker
diff -N SODmatinv.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SODmatinv.java 10 Mar 2006 21:21:35 -0000 1.1
@@ -0,0 +1,197 @@
+package org.lcsim.contrib.SODTracker;
+
+import java.util.*;
+
+public class SODmatinv
+ {
+ public SODmatinv()
+ {
+ ik=new int[10];
+ jk=new int[10];
+ array=new double[121];
+ for (int isrw=0; isrw<121; ++isrw){array[isrw]=0.0;}
+ }
+// public int matinv(double array, int norder, double det){
+ public int matinv(double[][] matrix, int snorder, double sdet){
+
+ norder=snorder; det=sdet;
+ /* Parameter adjustments */
+ for (i = 0; i < nmax; ++i){
+ for (j = 0; j < nmax; ++j){
+ array[(i+1)+(j+1)*nmax]=matrix[i][j];
+ }//end j
+ }//end i
+
+ /* Function Body */
+ det = (double)1.;
+ for (k = 1; k <= norder; ++k) {
+ amax = (double)0.;
+ L21();
+ if (amax == 0.) {sdet = (double)0.; return 1001;}
+ }
+ L100();
+ sdet=det;
+// for (int isrw=0; isrw<121; ++isrw){
+// System.out.print("array ");System.out.print(isrw);
+// System.out.print(" ");System.out.println(array[isrw]);
+// }
+//impedence match
+//weird,weird,weird
+ matrix[0][0]=-array[13];matrix[0][1]= array[12];matrix[1][1]= array[22];
+ matrix[0][2]=-array[33];matrix[1][2]= array[21];matrix[2][2]= array[31];
+ matrix[1][0]= array[12];matrix[2][0]=-array[33];matrix[2][1]= array[21];
+// for (i = 0; i < nmax; ++i){
+// for (j = 0; j < nmax; ++j){
+// matrix[i][j]=array[(i+1)+(j+1)*nmax];
+// }//end j
+// }//end i
+ return 0;
+} /* Dcxmatinv */
+//thanks java; no goto's _and_ no f2c
+ private void L21(){
+/* FIND LARGEST ELEMENT ARRAY(I, J) IN REST OF MATRIX */
+ for (i = k; i <= norder; ++i) {
+ for (j = k; j <= norder; ++j) {
+ d__1 = array[i + j * nmax];
+ if ((Math.abs(amax)-Math.abs(d__1)) <= 0.) {
+ amax = array[i + j * nmax];
+ ik[k - 1] = i;
+ jk[k - 1] = j;
+ }
+ }
+ }
+ if (amax == 0.) {return;}
+ L31();
+ }
+ private void L31(){
+/* INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K, K) */
+
+
+ i = ik[k - 1];
+ if ((i__3 = i - k) < 0) {
+ L21();
+ } else if (i__3 == 0) {
+ L51();
+ } else {
+ L43();
+ L51();//new, would have dropped thru to L51
+ }
+ }
+ private void L43(){
+ for (j = 1; j <= norder; ++j) {
+ save = array[k + j * nmax];
+ array[k + j * nmax] = array[i + j * nmax];
+ array[i + j * nmax] = -save;
+ }
+ }
+ private void L51(){
+ j = jk[k - 1];
+ if ((i__3 = j - k) < 0) {
+ L21();
+ } else if (i__3 == 0) {
+ L61();
+ } else {
+ L53();
+ L61();//new, would have dropped thru to L61
+ }
+ }
+ private void L53(){
+ for (i = 1; i <= norder; ++i) {
+ save = array[i + k * nmax];
+ array[i + k * nmax] = array[i + j * nmax];
+ array[i + j * nmax] = -save;
+ }
+ }
+ private void L61(){
+/* ACCUMULATE ELEMENTS OF INVERSE MATRIX */
+
+ for (i = 1; i <= norder; ++i) {
+ if (i - k != 0) {
+ array[i + k * nmax] = -array[i + k * nmax] / amax;
+ }
+ }
+ for (i = 1; i <= norder; ++i) {
+ for (j = 1; j <= norder; ++j) {
+ if (i - k != 0) {
+ L74();
+ } else {
+ L80();
+ }
+ }
+ }
+ for (j = 1; j <= norder; ++j) {
+ if (j - k != 0) {
+ array[k + j * nmax] /= amax;
+ }
+ }
+ array[k + k * nmax] = (double)1. / amax;
+ det *= amax;
+ }
+ private void L74(){
+ if (j - k != 0) {
+ L75();
+ } else {
+ L80();
+ }
+ }
+ private void L75(){
+ array[i+j*nmax] += array[i+k*nmax] * array[k+j*nmax];
+ }
+ private void L80(){
+ ;
+ }
+ private void L100(){
+/* RESTORE ORDERING OF MATRIX */
+
+ for (l = 1; l <= norder; ++l) {
+ k = norder - l + 1;
+ j = ik[k - 1];
+ if (j - k <= 0) {
+ L111();
+ } else {
+ L105();
+ }
+ }
+ }
+ private void L105(){
+ for (i = 1; i <= norder; ++i) {
+ save = array[i + k * nmax];
+ array[i + k * nmax] = -array[i + j * nmax];
+ array[i + j * nmax] = save;
+ }
+ }
+ private void L111(){
+ i = jk[k - 1];
+ if (i - k <= 0) {
+ L130();
+ } else {
+ L113();
+ }
+ }
+ private void L113(){
+ for (j = 1; j <= norder; ++j) {
+ save = array[k + j * nmax];
+ array[k + j * nmax] = -array[i + j * nmax];
+ array[i + j * nmax] = save;
+ }
+ }
+ private void L130(){
+ ;
+ }
+ /* System generated locals */
+ private int i__3;
+ private double d__1;
+ private double amax;
+ private double save;
+ private int i;
+ private int j;
+ private int k;
+ private int l;
+ private int[] ik;
+ private int[] jk;
+ private double[] array;
+ private int norder;
+ private double det;
+ /* Local variables */
+ private static int nmax=10;
+}
lcsim/src/org/lcsim/contrib/SODTracker/test
diff -N TestSOD.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TestSOD.java 10 Mar 2006 21:21:36 -0000 1.1
@@ -0,0 +1,26 @@
+package org.lcsim.contrib.SODTracker;
+import java.util.List;
+
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+public class TestSOD extends Driver {
+ private AIDA aida = AIDA.defaultInstance();
+ public TestSOD(){
+ add(new SODTrackFinderDriver());
+ };
+ public void process(EventHeader event) {
+ super.process(event); // this takes care that the child Drivers are loaded and processed.
+ List<SODTrack> SODTrackList = (List<SODTrack>) event.get("SODTracks");
+ aida.cloud1D("nTracks").fill(SODTrackList.size());
+ for ( SODTrack SODtrk : SODTrackList ) {
+ aida.cloud1D("Pt").fill(SODtrk.getPT());
+ aida.cloud1D("ptot").fill(SODtrk.getPtot());
+// List<SimTrackerHit> hitList = (List<SimTrackerHit>) SODtrk.getHits();
+ aida.cloud1D("nHits").fill(SODtrk.getNPoints());
+ }
+ }
+}
CVSspam 0.2.8