00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #if !defined(PACKAGE_LPAR_H_INCLUDED)
00016 #define PACKAGE_LPAR_H_INCLUDED
00017
00018
00019 #include <iosfwd>
00020 #include <cmath>
00021 #include <iostream>
00022
00023
00024 #include "CLHEP/Matrix/Vector.h"
00025 #include "CLHEP/Matrix/Matrix.h"
00026
00027 #include "CLHEP/Geometry/Point3D.h"
00028 #ifndef ENABLE_BACKWARDS_COMPATIBILITY
00029 typedef HepGeom::Point3D<double> HepPoint3D;
00030 #endif
00031
00032
00033 using CLHEP::HepVector;
00034 using CLHEP::Hep3Vector;
00035 using CLHEP::HepMatrix;
00036 using CLHEP::HepSymMatrix;
00037
00038 using namespace CLHEP;
00039
00040
00041
00042 class Lpar
00043 {
00044
00045
00046 public:
00047
00048
00049
00050 Lpar();
00051
00052 virtual ~Lpar();
00053
00054
00055 inline const Lpar& operator=( const Lpar& );
00056
00057
00058 inline void neg();
00059 void circle(double x1, double y1, double x2, double y2,
00060 double x3, double y3);
00061
00062
00063 double kappa() const { return m_kappa; }
00064 double radius() const { return 0.5/std::fabs(m_kappa);}
00065 HepVector center() const;
00066 double s(double x, double y) const;
00067 inline double d(double x, double y) const;
00068 inline double dr(double x, double y) const;
00069 double s(double r, int dir=0) const;
00070 double phi(double r, int dir=0) const;
00071 inline int sd(double r, double x, double y,
00072 double limit, double & s, double & d) const;
00073 inline HepVector Hpar(const HepPoint3D & pivot) const;
00074
00075
00076
00077
00078 friend class Lpav;
00079 friend std::ostream & operator<<(std::ostream &o, Lpar &);
00080 friend int intersect(const Lpar&, const Lpar&, HepVector&, HepVector&);
00081
00082 protected:
00083
00084
00085
00086
00087 private:
00088
00089 class Cpar {
00090 public:
00091 Cpar(const Lpar&);
00092 double xi() const{ return 1 + 2 * m_cu * m_da; }
00093 double sfi() const { return m_sfi; }
00094 double cfi() const { return m_cfi; }
00095 double da() const { return m_da; }
00096 double cu() const { return m_cu; }
00097 double fi() const { return m_fi; }
00098 private:
00099 double m_cu;
00100 double m_fi;
00101 double m_da;
00102 double m_sfi;
00103 double m_cfi;
00104 };
00105 friend class Lpar::Cpar;
00106
00107 inline Lpar( const Lpar& );
00108
00109
00110 bool operator==( const Lpar& ) const;
00111 bool operator!=( const Lpar& ) const;
00112
00113
00114 void scale(double s) { m_kappa /= s; m_gamma *= s; }
00115 inline void rotate(double c, double s);
00116 inline void move (double x, double y);
00117
00118
00119 double alpha() const { return m_alpha; }
00120 double beta() const { return m_beta; }
00121 double gamma() const { return m_gamma; }
00122 inline double check() const;
00123 HepMatrix dldc() const;
00124 inline double d0(double x, double y) const;
00125 double kr2g(double r) const { return m_kappa * r * r + m_gamma; }
00126 double x(double r) const;
00127 double y(double r) const;
00128 void xhyh(double x, double y, double&xh, double&yh) const;
00129 double xi2() const { return 1 + 4 * m_kappa * m_gamma; }
00130 bool xy(double, double&, double&, int dir=0) const;
00131 inline double r_max() const;
00132 double xc() const { return - m_alpha/2/m_kappa; }
00133 double yc() const { return - m_beta/2/m_kappa; }
00134 double da() const { return 2 * gamma() / (std::sqrt(xi2()) + 1); }
00135 inline double arcfun(double xh, double yh) const;
00136
00137
00138 double m_alpha;
00139 double m_beta;
00140 double m_gamma;
00141 double m_kappa;
00142
00143
00144 static const double BELLE_ALPHA;
00145
00146
00147
00148 };
00149
00150
00151
00152
00153
00154
00155
00156 inline Lpar::Lpar() {
00157 m_alpha = 0;
00158 m_beta = 1;
00159 m_gamma = 0;
00160 m_kappa = 0;
00161 }
00162
00163 inline Lpar::Lpar(const Lpar&l) {
00164 m_alpha = l.m_alpha;
00165 m_beta = l.m_beta;
00166 m_gamma = l.m_gamma;
00167 m_kappa = l.m_kappa;
00168 }
00169
00170 inline const Lpar&Lpar::operator=(const Lpar&l) {
00171 if (this != &l) {
00172 m_alpha = l.m_alpha;
00173 m_beta = l.m_beta;
00174 m_gamma = l.m_gamma;
00175 m_kappa = l.m_kappa;
00176 }
00177 return *this;
00178 }
00179
00180 inline void Lpar::rotate(double c, double s) {
00181 double aLpar = c * m_alpha + s * m_beta;
00182 double betar = -s * m_alpha + c * m_beta;
00183 m_alpha = aLpar;
00184 m_beta = betar;
00185 }
00186
00187 inline void Lpar::move (double x, double y) {
00188 m_gamma += m_kappa * ( x * x + y * y ) + m_alpha * x + m_beta * y;
00189 m_alpha += 2 * m_kappa * x;
00190 m_beta += 2 * m_kappa * y;
00191 }
00192
00193 inline double Lpar::check() const {
00194 return m_alpha * m_alpha + m_beta * m_beta - 4 * m_kappa * m_gamma - 1;
00195 }
00196
00197 inline void Lpar::neg() {
00198 m_alpha = -m_alpha;
00199 m_beta = -m_beta;
00200 m_gamma = -m_gamma;
00201 m_kappa = -m_kappa;
00202 }
00203
00204 inline double Lpar::d0(double x, double y) const {
00205 return m_alpha * x + m_beta * y + m_gamma + m_kappa * ( x * x + y * y);
00206 }
00207
00208 inline double Lpar::d(double x, double y) const {
00209 double dd = d0(x,y);
00210 const double approx_limit = 0.2;
00211 if(std::fabs(m_kappa*dd)>approx_limit) return -1;
00212 return dd * ( 1 - m_kappa * dd );
00213 }
00214
00215 inline double Lpar::dr(double x, double y) const {
00216 double dx = xc() - x;
00217 double dy = yc() - y;
00218 double r = 0.5/std::fabs(m_kappa);
00219 return std::fabs(std::sqrt(dx * dx + dy * dy) - r);
00220 }
00221
00222 inline double Lpar::r_max() const {
00223 if (m_kappa==0) return 100000000.0;
00224 if (m_gamma==0) return 1/std::fabs(m_kappa);
00225 return std::fabs(2*m_gamma/(std::sqrt(1+4*m_gamma*m_kappa)-1));
00226 }
00227
00228 inline double Lpar::arcfun(double xh, double yh) const {
00229
00230
00231
00232 double r2kap = 2.0 * m_kappa;
00233 double xi = std::sqrt(xi2());
00234 double xinv = 1.0 / xi;
00235 double ar2kap = std::fabs(r2kap);
00236 double cross = m_alpha * yh - m_beta * xh;
00237 double a1 = ar2kap * cross * xinv;
00238 double a2 = r2kap * (m_alpha * xh + m_beta * yh) * xinv + xi;
00239 if (a1>=0 && a2>0 && a1<0.3) {
00240 double arg2 = a1*a1;
00241 return cross * ( 1.0 + arg2 * (1./6. + arg2 * (3./40.))) * xinv;
00242 } else {
00243 double at2 = std::atan2(a1,a2);
00244 if (at2<0) at2 += (2*M_PI);
00245 return at2/ar2kap;
00246 }
00247 }
00248
00249 inline int Lpar::sd(double r, double x, double y,
00250 double limit, double & s, double & d) const{
00251 if ((x*yc()-y*xc())*m_kappa<0) return 0;
00252 double dd = d0(x,y);
00253 d = dd * ( 1 - m_kappa * dd );
00254 double d_cross_limit = d*limit;
00255 if (d_cross_limit < 0 || d_cross_limit > limit*limit) return 0;
00256
00257 double rc = std::sqrt(m_alpha*m_alpha+m_beta*m_beta)/(2*m_kappa);
00258 double rho = 1./(-2*m_kappa);
00259 double cosPhi = (rc*rc + rho*rho - r*r)/(-2*rc*rho);
00260 double phi = std::acos(cosPhi);
00261 s = std::fabs(rho)*phi;
00262 d *= r/(std::fabs(rc)*std::sin(phi));
00263 if (abs(d) > abs(limit)) return 0;
00264 d_cross_limit = d*limit;
00265 if (d_cross_limit > limit*limit) return 0;
00266 return 1;
00267 }
00268
00269 inline HepVector Lpar::Hpar(const HepPoint3D & pivot) const{
00270 HepVector a(5);
00271 double dd = d0(pivot.x(),pivot.y());
00272 a(1) = dd * ( m_kappa * dd - 1 );
00273 a(2) = (m_kappa>0) ? std::atan2(yc() - pivot.y(), xc() - pivot.x()) + M_PI
00274 : std::atan2(pivot.y() - yc(), pivot.x() - xc()) - M_PI;
00275 a(3) = -2.0*BELLE_ALPHA*m_kappa;
00276 a(4) = 0;
00277 a(5) = 0;
00278 return a;
00279 }
00280
00281 #endif
00282