1 added + 8 modified, total 9 files
lcsim/src/org/lcsim/contrib/SODTracker
diff -N SODFittedHel.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SODFittedHel.java 28 Jul 2006 20:21:41 -0000 1.1
@@ -0,0 +1,269 @@
+package org.lcsim.contrib.SODTracker;
+
+import java.util.*;
+
+public class SODFittedHel extends SODHel
+ {
+ public SODFittedHel()
+ {
+// new
+ current_hel=new SODHel();
+ next_hel=new SODHel();
+ _hitlist=new java.util.LinkedList();
+ nhits=0;
+ ematrix = new double[5][5];
+ }
+ public void make_fitted_hel(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;
+ 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_hel(SODFittedHel 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 >= 2*nhits) {return ftemp;}
+//more involved now int ndof=nhits-nfree;
+ ndof=0; nphits=0; nzhits=0;
+ 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){
+ ndof=0; nphits=0; nzhits=0;
+ ftemp=DoFit(); chisq_mem[iloop]=chisq; ndof=nphits+nzhits-nfree;
+ 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.print(" ");System.out.print(ssqd0);
+ System.out.print(" ");System.out.print(ssqz0);
+ 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=5;
+ 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[5][5];
+ double[] B = new double[5];
+ double[] D = new double[5];
+ 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;}
+ }
+//first loop over hits for phi derivatives
+ for (ListIterator ii = _hitlist.listIterator(); ii.hasNext();){
+ SODHit shit=(SODHit)ii.next();
+ double [] numb_derivs=new double[6];
+ boolean whaa=shit.calc_p_derivatives(current_hel,numb_derivs);
+ prtflg=0;
+ if (prtflg!=0){
+ System.out.print("pderiv ");
+ for (int jj=0; jj<6; ++jj){System.out.print(" ");System.out.print(numb_derivs[jj]);}
+ System.out.println(" ");
+ }
+ chisq+=numb_derivs[0]*numb_derivs[0];
+ for(int ipar=0; ipar<norder; ipar++){
+ D[ipar]+=numb_derivs[0]*numb_derivs[ipar+1];
+ //inner parameter loop
+ for(int jpar=0; jpar<norder; jpar++){
+ A[ipar][jpar]+=numb_derivs[ipar+1]*numb_derivs[jpar+1];
+ }//endof inner parameter loop
+ }//endof outer parameter loop
+ nphits++; //now _count_ phi and z hits separately
+ }//endof pointloop
+//then loop over hits for z derivatives
+ for (ListIterator ii = _hitlist.listIterator(); ii.hasNext();){
+ SODHit shit=(SODHit)ii.next();
+ if (shit.has_z_info()){
+ double [] numb_derivs=new double[6];
+ boolean whaa=shit.calc_z_derivatives(current_hel,numb_derivs);
+ prtflg=0;
+ if (prtflg!=0){
+ System.out.print("zderiv ");
+ for (int jj=0; jj<6; ++jj){System.out.print(" ");System.out.print(numb_derivs[jj]);}
+ System.out.println(" ");
+ }
+ chisq+=numb_derivs[0]*numb_derivs[0];
+ for(int ipar=0; ipar<norder; ipar++){
+ D[ipar]+=numb_derivs[0]*numb_derivs[ipar+1];
+ //inner parameter loop
+ for(int jpar=0; jpar<norder; jpar++){
+ A[ipar][jpar]+=numb_derivs[ipar+1]*numb_derivs[jpar+1];
+ }//endof inner parameter loop
+ }//endof outer parameter loop
+ nzhits++; //now _count_ phi and z hits separately
+ }//end shit has_z_info
+ }//endof pointloop
+ if ((nphits<4)||(nzhits<3)){
+ System.out.print("DANGER - HELIX FIT UNSTABLE ");
+ System.out.print(nphits);System.out.print(" ");System.out.println(nzhits);
+ }
+ prtflg=0;
+ if (prtflg!=0){
+ System.out.print("chisq before fit ");System.out.print(chisq);
+// System.out.println(" ");
+ System.out.print(" "); current_hel.printHelix();
+ }
+ 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);
+//if matrix inverted, update parameters
+ 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()-0.5*B[3],current_hel.Tanl()-0.5*B[4]);
+ ssqd0=A[0][0]; ssqz0=A[3][3];
+ for(int ii=0;ii<norder;ii++){for(int jj=0;jj<norder;jj++){ematrix[ii][jj]=A[ii][jj];}}
+// new_hel.printHelix();
+ next_hel=new_hel;
+ }//end matrix inverted
+ return ftemp;
+}
+
+ boolean check_inversion( double[][] A, double[][] Z) {
+ int prtflg=0;
+ boolean rtflg=false; int norder=5;
+ double[][] I = new double[5][5];
+ 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]
+ +A[ipar][3]*Z[3][jpar]+A[ipar][4]*Z[4][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[3][3]-1.0)<epsilon)&&(Math.abs(I[4][4]-1.0)<epsilon)&&
+ (Math.abs(I[0][1] )<epsilon)&&(Math.abs(I[0][2] )<epsilon)&&(Math.abs(I[0][3] )<epsilon)&&
+ (Math.abs(I[0][4] )<epsilon)&&(Math.abs(I[1][2] )<epsilon)&&(Math.abs(I[1][3] )<epsilon)&&
+ (Math.abs(I[1][4] )<epsilon)&&(Math.abs(I[2][3] )<epsilon)&&(Math.abs(I[2][4] )<epsilon)&&
+ (Math.abs(I[3][4] )<epsilon))rtflg=true;
+ if (!rtflg)System.out.println("MATRIX INVERSION FAILED!!!!!");
+ 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.print(A[ii][2]);
+ System.out.print(" ");System.out.print(A[ii][3]);
+ System.out.print(" ");System.out.println(A[ii][4]);
+ }
+ 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.print(Z[ii][2]);
+ System.out.print(" ");System.out.print(Z[ii][3]);
+ System.out.print(" ");System.out.println(Z[ii][4]);
+ }
+ 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.print(I[ii][2]);
+ System.out.print(" ");System.out.print(I[ii][3]);
+ System.out.print(" ");System.out.println(I[ii][4]);
+ }
+ }//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;}
+ public double[][] get_ematrix(){return ematrix;}
+ private double chisq;
+ private double ssqd0;
+ private double ssqz0;
+ private double[][] ematrix;
+ private int norder;
+ private int fail;
+ private int nhits;
+ private int ndof;
+ private int nphits;
+ private int nzhits;
+ private static int nfree=5;
+ 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 -u -r1.1 -r1.2
--- History 10 Mar 2006 21:21:34 -0000 1.1
+++ History 28 Jul 2006 20:21:41 -0000 1.2
@@ -3,6 +3,36 @@
## Please summarize changes to the SODTracker code here.
## Most recent first please.
+28 July 2006 Fred Blanc V01-01-00
+Modifications from Steve Wagner to include helix fit.
+SODTrack now contains the fitted helix parameters and the
+error matrix. Further details of the modifications are
+given below.
+
+SODHel.java - 2 new fcns (l1_at_radius and z1_at_radius)
+
+SODHit.java - 6 new fcns (make_hit, z, has_z_info, calc_p_derivatives,
+ calc_z_derivatives, calc_z_resid)
+
+SODTrackFinder.java - 2 bug fixes, 1 new fcn (l1_at_radius), and
+ major changes to make_test_helix and MC version
+ of MC truth version of SODTrackFinder to get it
+ to use z info and SODFittedHel
+
+SODTrackFinderDriver.java - make SODHits with z info for vtx det hits
+
+SODmatinv.java - correctly handle 5x5 inversion. Really needs to be
+ fixed for the general case
+
+SODFittedHel.java - all new
+
+SODTrack.java - added capability to fill in and to extract the
+ error matrix.
+
+test/TestSOD.java - add a few hists
+
+
+
10 March 2006 Fred Blanc V01-00-00
- Create SODTracker package. Included files:
SODFittedCir.java
lcsim/src/org/lcsim/contrib/SODTracker
diff -u -r1.1 -r1.2
--- SODHel.java 10 Mar 2006 21:21:34 -0000 1.1
+++ SODHel.java 28 Jul 2006 20:21:41 -0000 1.2
@@ -178,6 +178,20 @@
// derivs[1]=derivs[1]/100.;derivs[2]=derivs[2]/10000.;derivs[3]=derivs[3]/10000.;
return true;
}
+ public double l1_at_radius(double rl)
+ {
+ 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;
+ return l1;
+ }
+ public double z1_at_radius(double rl)
+ {
+ double z1 = z0 + tanl*l1_at_radius(rl);
+ return z1;
+ }
private double d0;
private double phi0;
private double omega;
lcsim/src/org/lcsim/contrib/SODTracker
diff -u -r1.2 -r1.3
--- SODHit.java 12 Jul 2006 18:15:55 -0000 1.2
+++ SODHit.java 28 Jul 2006 20:21:41 -0000 1.3
@@ -23,6 +23,7 @@
double[] hitpos = h.getPosition();
_hitx = hitpos[0];
_hity = hitpos[1];
+ _has_z_info = false;
_hitrad = Math.sqrt(_hitx*_hitx+_hity*_hity);
_hitlay=0;
// if _hitlay=0 it's not a SODHIT and other data may be invalid,
@@ -50,6 +51,7 @@
_is_foreground_hit=true;
_hitx = xin*mmtocm;
_hity = yin*mmtocm;
+ _has_z_info = false;
_hitrad = Math.sqrt(_hitx*_hitx+_hity*_hity);
_hitlay=lay+1;
_hitphi = Math.atan2(_hity,_hitx);
@@ -60,12 +62,32 @@
_trackerHit = trackerHit;
_mcParticle = mcp;
}
+ public void make_hit(int lay, double xin, double yin, double zin, double rgauss, double rg2, SimTrackerHit trackerHit, MCParticle mcp)
+ {
+ _used_on_hel=0;
+ _is_foreground_hit=true;
+ _hitx = xin*mmtocm;
+ _hity = yin*mmtocm;
+ _hitz = zin*mmtocm;
+ _has_z_info = true;
+ _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);
+ _hitz+=rg2*myres;
+ _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;
+ _has_z_info = false;
_hitrad = Math.sqrt(_hitx*_hitx+_hity*_hity);
_hitlay=lay+1;
_hitphi = Math.atan2(_hity,_hitx);
@@ -81,6 +103,7 @@
_is_foreground_hit=false;
_hitx = xrad*Math.cos(xphi);
_hity = xrad*Math.sin(xphi);
+ _has_z_info = false;
_hitrad = xrad;
_hitlay=0;
// if _hitlay=0 it's not a SODHIT and other data may be invalid,
@@ -100,8 +123,10 @@
public int layer(){return _hitlay;}
public double x(){return _hitx;}
public double y(){return _hity;}
+ public double z(){return _hitz;}
public double r(){return _hitrad;}
public double phi(){return _hitphi;}
+ public boolean has_z_info(){return _has_z_info;}
// 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)
@@ -128,6 +153,39 @@
calc_resid(new_hel);derivs[3]=(get_resid()-resid_base)/(delta*error);
return true;
}
+ public boolean calc_p_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);
+ derivs[4]=0.0; derivs[5]=0.0;
+ return true;
+ }
+ public boolean calc_z_derivatives(SODHel t, double [] derivs)
+ {
+ double resid_base=calc_z_resid(t); 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());
+ derivs[1]=(calc_z_resid(new_hel)-resid_base)/(delta*error);
+ new_hel.make_hel(t.D0(),t.Phi0()+delta,t.Omega(),t.Z0(),t.Tanl());
+ derivs[2]=(calc_z_resid(new_hel)-resid_base)/(delta*error);
+ new_hel.make_hel(t.D0(),t.Phi0(),t.Omega(),t.Z0()+delta,t.Tanl());
+ derivs[4]=(calc_z_resid(new_hel)-resid_base)/(delta*error);
+ new_hel.make_hel(t.D0(),t.Phi0(),t.Omega(),t.Z0(),t.Tanl()+delta);
+ derivs[5]=(calc_z_resid(new_hel)-resid_base)/(delta*error);
+ delta=0.01*t.Omega();
+ new_hel.make_hel(t.D0(),t.Phi0(),t.Omega()+delta,t.Z0(),t.Tanl());
+ derivs[3]=(calc_z_resid(new_hel)-resid_base)/(delta*error);
+ return true;
+ }
public boolean calc_resid(SODHel t)
{
double phi1 = t.phi1_at_radius(_hitrad);
@@ -145,6 +203,11 @@
}//end printing
return true;
}
+ public double calc_z_resid(SODHel t)
+ {
+ double _z_resid = _hitz - t.z1_at_radius(_hitrad);
+ return _z_resid;
+ }
public boolean calc_resid(TrkParams t)
{
double phi1 = phi1_at_radius(t,_hitrad);
@@ -213,10 +276,12 @@
private int _hitlay;
private double _hitx;
private double _hity;
+ private double _hitz;
private double _hitrad;
private double _hitphi;
private double _resid;
private boolean _is_foreground_hit;
+ private boolean _has_z_info;
public int _used_on_hel;
private static double pi=3.14159265359;
private static double twopi=6.283185307;
lcsim/src/org/lcsim/contrib/SODTracker
diff -u -r1.1 -r1.2
--- SODTrack.java 10 Mar 2006 21:21:35 -0000 1.1
+++ SODTrack.java 28 Jul 2006 20:21:41 -0000 1.2
@@ -32,11 +32,11 @@
pz = pZ;
p = Math.sqrt(px*px + py*py + pz*pz);
momentum = new double[] {px, py, pz};
- d0 = dzero;
+ d0 = 10.*dzero;
phi0 = phizero;
// CAUTION: different definition of omega in SOD and lcsim!
omega = -0.1*omeg;
- z0 = zzero;
+ z0 = 10.*zzero;
s = spar;
parameters = new double[] {d0, phi0, omega,z0,s};
}
@@ -66,8 +66,20 @@
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; }
+
+ private double[][] ematrix;
+ public void add_ematrix(double[][] emat)
+ {
+ ematrix = new double[5][5];
+ double[] conv = new double[] {10.,1.,-0.1,10.,1.};
+ for(int i=0; i<5; i++) {
+ for(int j=0; j<5; j++) {
+ ematrix[i][j] = conv[i]*conv[j]*emat[i][j];
+ }
+ }
+ }
+ public double getErrorMatrixElement(int i, int j) { return ematrix[i][j]; }
+ public double[][] getErrorMatrix() { return ematrix; }
public double getChi2() { return 1.; }
public int getNDF() { return 1; }
lcsim/src/org/lcsim/contrib/SODTracker
diff -u -r1.3 -r1.4
--- SODTrackFinder.java 15 Jul 2006 23:04:15 -0000 1.3
+++ SODTrackFinder.java 28 Jul 2006 20:21:41 -0000 1.4
@@ -2,8 +2,11 @@
import java.util.*;
import org.lcsim.event.MCParticle;
+import org.lcsim.util.aida.AIDA;
+
public class SODTrackFinder
{
+ private AIDA aida = AIDA.defaultInstance();
public SODTrackFinder()
{
_trklist = new java.util.LinkedList();
@@ -44,8 +47,15 @@
java.util.LinkedList shitlist = new java.util.LinkedList();
shitlist.add(sh1);shitlist.add(sh2);shitlist.add(sh3);
shitlist.add(sh4);shitlist.add(sh5);
+//move to helix fit SODFittedCir fit_helix = new SODFittedCir();
+//move to helix fit fit_helix.make_fitted_cir(shitlist,test_helix);
+ SODFittedHel fit_real_helix = new SODFittedHel();
+ fit_real_helix.make_fitted_hel(shitlist,test_helix);
+ SODHel real_helix = (SODHel)fit_real_helix;
SODFittedCir fit_helix = new SODFittedCir();
- fit_helix.make_fitted_cir(shitlist,test_helix);
+ aida.cloud1D("stubpt").fill(real_helix.Pt());
+ fit_helix.make_fitted_cir(fit_real_helix.hitlist(),real_helix);
+//move to helix fit
boolean sflag=false;
if (0==fit_helix.Fail())sflag=_trklist.add(fit_helix);
//
@@ -65,7 +75,7 @@
} // if lists not empty
} else { // if(useMCTruth) {
-
+
// System.out.println("find seed track from combinatorics")
int prtflg=0;
_trklist = new java.util.LinkedList(); int ncombo=0;
@@ -177,12 +187,12 @@
}//loop over l2 hits
}//(unindented) all list have hits
} // if(useMCTruth) { } // contructor
- }
+ }//close bracket deleated by mistake?
-// 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 z0=0.0, tanl=0.0;
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);
@@ -201,35 +211,27 @@
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;
+ 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;
}
}
@@ -242,15 +244,20 @@
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);
+ if ((sh1.has_z_info())&&(sh2.has_z_info())){
+ double r1=Math.sqrt(x1*x1+y1*y1); double l1=l1_at_radius(d0,phi0,omega,r1);
+ double r2=Math.sqrt(x2*x2+y2*y2); double l2=l1_at_radius(d0,phi0,omega,r2);
+ tanl=(sh1.z()-sh2.z())/(l1-l2); z0=sh1.z()-tanl*l1;
+// System.out.print("test_helix ");System.out.print(l1);
+// System.out.print(" ");System.out.print(l2);
+// System.out.print(" ");System.out.print(z0);
+// System.out.print(" ");System.out.println(tanl);
+ }
+ SODHel test_helix=new SODHel();
+ test_helix.make_hel(d0,phi0,omega,z0,tanl);
return test_helix;
}
public double search_fcn(double pt, int lay){
@@ -261,6 +268,15 @@
if ((search<0.0141)&&(lay==5))search=0.0141;
return search;
}
+ public double l1_at_radius(double d0, double phi0, double omega, double rl)
+ {
+ 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;
+ return l1;
+ }
public java.util.LinkedList trklist(){return _trklist;}
private java.util.LinkedList _trklist;
private static double pi=3.14159265359;
lcsim/src/org/lcsim/contrib/SODTracker
diff -u -r1.1 -r1.2
--- SODTrackFinderDriver.java 10 Mar 2006 21:21:35 -0000 1.1
+++ SODTrackFinderDriver.java 28 Jul 2006 20:21:42 -0000 1.2
@@ -56,7 +56,9 @@
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);
+ // sh.make_hit(layer,pos[0],pos[1],rgauss,trackerHit,mcp);
+ double rg2=generator.nextGaussian();
+ sh.make_hit(layer,pos[0],pos[1],pos[2],rgauss,rg2,trackerHit,mcp);
} else {
sh.make_hit(layer,pos[0],pos[1],rgauss,trackerHit);
}
@@ -102,7 +104,7 @@
//
//c now run SODHitAdder, which really is the outer trker trkfinder
//
- ntrk=0;
+ ntrk=0; int nhots=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)&&
@@ -114,7 +116,25 @@
SODFittedCir out_helix= new SODFittedCir();
out_helix=sodha.get_new_fit();
boolean success=shatrkl.add(out_helix);
- SODHel test_helix=(SODHel)out_helix;
+// SODHel test_helix=(SODHel)out_helix;
+//do full helix fit now
+//resid SODHit reshit=new SODHit();
+ SODHel test_cir=(SODHel)out_helix; java.util.LinkedList shitlist = new java.util.LinkedList();
+ for ( ListIterator hitItr = out_helix.hitlist().listIterator() ; hitItr.hasNext();){
+ SODHit hit = (SODHit)hitItr.next();
+//resid if ((hit.layer()==1)&&(hit.r()>15.0)){
+//resid reshit=hit; shitlist.add(hit);
+//resid }else{
+ shitlist.add(hit);
+//resid }
+// System.out.print(hit.layer());System.out.print(" ");System.out.println(hit.r());
+ }
+ SODFittedHel test_helix= new SODFittedHel(); test_helix.make_fitted_hel(shitlist,test_cir);
+//resid SODHel temp=(SODHel)test_helix; reshit.calc_resid(temp);
+//resid System.out.print("hit residual ");System.out.println(reshit.get_resid());
+ double[][] ematrix=new double[5][5]; ematrix=test_helix.get_ematrix();
+ // System.out.print("ssqd0 ");System.out.print(ematrix[0][0]);
+ // System.out.print(" ssqz0 ");System.out.println(ematrix[3][3]);
ntrk++;
// create new Tracks
double l = 0;
@@ -123,11 +143,12 @@
test_helix.D0(), test_helix.Phi0(),
test_helix.Omega(), test_helix.Z0(),
test_helix.Tanl());
+ trk.add_ematrix(test_helix.get_ematrix());
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);
+ trk.addHit(h); nhots++;
}
listTrk.add(trk);
}//end loop over found trks
@@ -146,6 +167,7 @@
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(nhots);
System.out.print(" ");System.out.print(l1vhlist.size());
System.out.print(" ");System.out.print(l2vhlist.size());
System.out.print(" ");System.out.print(l3vhlist.size());
lcsim/src/org/lcsim/contrib/SODTracker
diff -u -r1.1 -r1.2
--- SODmatinv.java 10 Mar 2006 21:21:35 -0000 1.1
+++ SODmatinv.java 28 Jul 2006 20:21:42 -0000 1.2
@@ -25,6 +25,7 @@
/* Function Body */
det = (double)1.;
for (k = 1; k <= norder; ++k) {
+/* FIND LARGEST ELEMENT ARRAY(I, J) IN REST OF MATRIX */
amax = (double)0.;
L21();
if (amax == 0.) {sdet = (double)0.; return 1001;}
@@ -37,19 +38,32 @@
// }
//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
+//obviously some small but not insignificant part of matinv didn't get translated properly
+//from C++ to java. Will have to sort it out someday but for now just use the patches derived
+//from stand-alone C++ version. At least starting to see a pattern now.
+// if (3==norder){
+// matrix[0][0]= array[11];matrix[0][1]= array[12];matrix[0][2]= array[13];
+// matrix[1][0]= array[21];matrix[1][1]= array[22];matrix[1][2]= array[23];
+// matrix[2][0]= array[31];matrix[2][1]= array[32];matrix[2][2]= array[33];
+// }
+// if (5==norder){
+// matrix[0][0]= array[11];matrix[0][1]= array[12];matrix[0][2]= array[13];matrix[0][3]= array[14];matrix[0][4]= array[15];
+// matrix[1][0]= array[21];matrix[1][1]= array[22];matrix[1][2]= array[23];matrix[1][3]= array[24];matrix[1][4]= array[25];
+// matrix[2][0]= array[31];matrix[2][1]= array[32];matrix[2][2]= array[33];matrix[2][3]= array[34];matrix[2][4]= array[35];
+// matrix[3][0]= array[41];matrix[3][1]= array[42];matrix[3][2]= array[43];matrix[3][3]= array[44];matrix[3][4]= array[45];
+// matrix[4][0]= array[51];matrix[4][1]= array[52];matrix[4][2]= array[53];matrix[4][3]= array[54];matrix[4][4]= array[55];
+// }
+//should be all fixed now
+ 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 */
+//System.out.println("L21");
for (i = k; i <= norder; ++i) {
for (j = k; j <= norder; ++j) {
d__1 = array[i + j * nmax];
@@ -60,12 +74,9 @@
}
}
}
- if (amax == 0.) {return;}
- L31();
- }
- private void L31(){
/* INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K, K) */
-
+//this return just exits L21 - next line in Function Body is real return
+ if (amax == 0.) {return;}
i = ik[k - 1];
if ((i__3 = i - k) < 0) {
@@ -74,17 +85,19 @@
L51();
} else {
L43();
- L51();//new, would have dropped thru to L51
}
}
private void L43(){
+//System.out.println("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;
}
+ L51();//new, need to drop thru to L51
}
private void L51(){
+//System.out.println("L51");
j = jk[k - 1];
if ((i__3 = j - k) < 0) {
L21();
@@ -92,17 +105,19 @@
L61();
} else {
L53();
- L61();//new, would have dropped thru to L61
}
}
private void L53(){
+//System.out.println("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;
}
+ L61();//new, need to drop thru to L61
}
private void L61(){
+//System.out.println("L61");
/* ACCUMULATE ELEMENTS OF INVERSE MATRIX */
for (i = 1; i <= norder; ++i) {
@@ -128,6 +143,7 @@
det *= amax;
}
private void L74(){
+//System.out.println("L74");
if (j - k != 0) {
L75();
} else {
@@ -135,12 +151,17 @@
}
}
private void L75(){
+//System.out.println("L75");
array[i+j*nmax] += array[i+k*nmax] * array[k+j*nmax];
+ L80();
}
private void L80(){
+//System.out.println("L80");
;
}
+
private void L100(){
+//System.out.println("L100");
/* RESTORE ORDERING OF MATRIX */
for (l = 1; l <= norder; ++l) {
@@ -153,14 +174,19 @@
}
}
}
+
private void L105(){
+//System.out.println("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;
}
+ L111();//new, need to drop thru to L111
}
+
private void L111(){
+//System.out.println("L111");
i = jk[k - 1];
if (i - k <= 0) {
L130();
@@ -168,14 +194,19 @@
L113();
}
}
+
private void L113(){
+//System.out.println("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;
}
+ L130();
}
+
private void L130(){
+//System.out.println("L130");
;
}
/* System generated locals */
lcsim/src/org/lcsim/contrib/SODTracker/test
diff -u -r1.2 -r1.3
--- TestSOD.java 10 Mar 2006 21:40:09 -0000 1.2
+++ TestSOD.java 28 Jul 2006 20:21:43 -0000 1.3
@@ -1,5 +1,7 @@
import java.util.List;
+import java.util.*;
+import java.io.*;
import org.lcsim.event.EventHeader;
import org.lcsim.event.MCParticle;
import org.lcsim.event.TrackerHit;
@@ -10,16 +12,54 @@
private AIDA aida = AIDA.defaultInstance();
public TestSOD(){
add(new SODTrackFinderDriver());
+ // 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();
+ // }
};
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 ) {
+ if (SODTrackList.size()<100){
+ aida.cloud1D("nTracks").fill(SODTrackList.size());
+ for ( SODTrack SODtrk : SODTrackList ) {
aida.cloud1D("Pt").fill(SODtrk.getPT());
aida.cloud1D("ptot").fill(SODtrk.getPtot());
+ aida.cloud1D("d0 ").fill(SODtrk.getTrackParameter(0));
+ aida.cloud1D("phi0").fill(SODtrk.getTrackParameter(1));
+ aida.cloud1D("z0 ").fill(SODtrk.getTrackParameter(3));
+ aida.cloud1D("tanl").fill(SODtrk.getTrackParameter(4));
// List<SimTrackerHit> hitList = (List<SimTrackerHit>) SODtrk.getHits();
aida.cloud1D("nHits").fill(SODtrk.getNPoints());
- }
+ // System.out.print("pttup ");System.out.print(SODtrk.getNPoints());
+ // System.out.print(" ");System.out.println(SODtrk.getPT());
+
+// double[][] emat = SODtrk.getErrorMatrix();
+// aida.cloud1D("d0Pull").fill(SODtrk.getTrackParameter(0)/Math.sqrt(emat[0][0]));
+// aida.cloud1D("phiErr").fill(Math.sqrt(emat[1][1]));
+// aida.cloud1D("omega").fill(SODtrk.getTrackParameter(2));
+// aida.cloud1D("omegaPull").fill((Math.abs(SODtrk.getTrackParameter(2))-0.00015)/Math.sqrt(emat[2][2]));
+// aida.cloud1D("omegaErr").fill(Math.sqrt(emat[2][2]));
+// aida.cloud1D("z0Pull").fill(SODtrk.getTrackParameter(3)/Math.sqrt(emat[3][3]));
+// aida.cloud1D("tanlPull").fill(SODtrk.getTrackParameter(4)/Math.sqrt(emat[4][4]));
+
+// for(int i=0; i<5; i++) {
+// System.out.print(" ");
+// System.out.print(SODtrk.getTrackParameter(i));
+// System.out.print(" | ");
+// for(int j=0; j<5; j++) {
+// System.out.print(" ");
+// System.out.print(emat[i][j]);
+// }
+// System.out.println("");
+// }
+
+ }
+ }
}
}
CVSspam 0.2.8