View Javadoc

1   package org.lcsim.recon.cat;
2   
3   //import hep.physics.*;
4   import java.util.*;
5   
6   /**
7    * Simple circle from 2D point calculation for helix estimates.
8    *
9    * @author  E. von Toerne
10   * @version $Id: CircleFromPoints.java,v 1.1 2007/04/06 21:48:14 onoprien Exp $
11   */
12  final public class CircleFromPoints {
13    public double[] a,b,c;
14    private double[] ab,ac,cb,caperp;
15    private double sin_a, sin_ab;
16    public double[] center;
17    public double kappa;
18    public double[] dirAtC;
19    // constructors
20    
21    /**
22     * empty constructor
23     */
24    CircleFromPoints() {
25      a=new double[]{0.,0.};
26      b=new double[]{0.,0.};
27      c=new double[]{0.,0.};
28      ab=new double[]{0.,0.};
29      ac=new double[]{0.,0.};
30      cb=new double[]{0.,0.};
31      caperp=new double[]{0.,0.};
32      center=new double[]{0.,0.};
33      dirAtC=new double[]{0.,0.};
34    }
35    
36    public void calculate(double ax, double ay, double bx, double by, double cx, double cy) {
37      int i;
38      a[0]=ax;
39      a[1]=ay;
40      b[0]=bx;
41      b[1]=by;
42      c[0]=cx;
43      c[1]=cy;
44      for (i=0;i<2;i++){
45        ab[i]=b[i]-a[i];
46        ac[i]=c[i]-a[i];
47        cb[i]=c[i]-b[i];
48      }
49      
50      double lencb=leng(cb);
51      double lenac=leng(ac);
52      if (lenac<1.E-5){
53        kappa = 0.;
54        System.out.println("Problem in CircleFromPoints"+lencb);
55        return;
56      }
57      if (lencb<1.E-5){
58        kappa = 0.;
59        dirAtC[0]=(c[0]-a[0])/lenac;
60        dirAtC[1]=(c[1]-a[1])/lenac;
61        //System.out.println("Problem in CircleFromPoints"+lencb);
62        return;
63      }
64      
65      caperp[0]=ac[1]/lenac;
66      caperp[1]=-ac[0]/lenac;
67      
68      double sin_a=Math.abs(cdot(ab,cb)/leng(ab)/lencb);
69      double sin_b=Math.abs(cdot(ac,cb)/leng(ac)/lencb);
70      double cos_a=1.-sqr(sin_a);
71      cos_a = cos_a >0.? Math.sqrt(cos_a) : 0.;
72      cos_a = cos_a <1.? cos_a : 1.;
73      double cos_b=1.-sqr(sin_b);
74      cos_b = cos_b >0.? Math.sqrt(cos_b) : 0.;
75      cos_b = cos_b<1.? cos_b : 1.;
76      double sin_ab=sin_a * cos_b + sin_b * cos_a;
77      double theSign=1.;
78      if (cdot(caperp,cb)<0.) theSign = -1.;
79      
80      for (i=0;i<2;i++) center[i]=0.5*(a[i]+c[i])+theSign*0.5*sin_a/sin_ab*lencb*caperp[i];
81      kappa=1./Math.sqrt(sqr(center[0]-a[0])+sqr(center[1]-a[1]));
82      dirAtC[0]=(center[1]-c[1])*kappa;
83      dirAtC[1]=-(center[0]-c[0])*kappa;
84      if (dirAtC[0]*(c[0]-b[0])+dirAtC[1]*(c[1]-b[1])<0.){
85        dirAtC[0]=-dirAtC[0];
86        dirAtC[1]=-dirAtC[1];
87      } else { kappa=-kappa;}
88      
89    }
90    
91    private double sqr(double x){return x*x;}
92    private double cdot(double[] a1, double[] b1){return a1[0]*b1[0]+a1[1]*b1[1];}
93    private double leng(double[] a1){return Math.sqrt(a1[0]*a1[0]+a1[1]*a1[1]);}
94    public double getCenter(int i){return center[i];}
95    public double getKappa(){return kappa;}
96    public double getDirAtC(int i){return dirAtC[i];}
97    public void debug(){
98      double[] disa=new double[]{center[0]-a[0],center[1]-a[1]};
99      double[] disb=new double[]{center[0]-b[0],center[1]-b[1]};
100     double[] disc=new double[]{center[0]-c[0],center[1]-c[1]};
101     System.out.println("a="+a[0]+" "+a[1]);
102     System.out.println("b="+b[0]+" "+b[1]);
103     System.out.println("c="+c[0]+" "+c[1]);
104     System.out.println("center="+center[0]+" "+center[1]);
105     System.out.println("kappa="+kappa+" Radius="+1./kappa);
106     System.out.println("distance point A to center="+leng(disa));
107     System.out.println("distance point B to center="+leng(disb));
108     System.out.println("distance point C to center="+leng(disc));
109   }
110 }