00001 #ifndef _ZAV_H 00002 #define _ZAV_H 00003 //#include "mscl.h" 00004 #include <iostream> 00005 00006 #include "CLHEP/Matrix/Matrix.h" 00007 #include "CLHEP/Matrix/SymMatrix.h" 00008 //#include "zhit.h" 00009 00010 using namespace CLHEP; 00011 00012 class zav { 00013 public: 00014 zav(); 00015 zav(int) {} // dummy for one of the chain constructor 00016 void set(const zav*); 00017 double chisq() const { return _chisq; } 00018 void add(double , double, double); 00019 //void add(const zhit&); 00020 double calculate(); 00021 double a() const { return _a; } 00022 double b() const { return _b; } 00023 double z(double s) const { return _a * s + _b; } 00024 double d(double s, double z) const { return z - _a * s - _b; } 00025 int nc() const { return _nc; } 00026 void clear(void); 00027 inline friend std::ostream&operator<<(std::ostream&, const zav&); 00028 //HepSymMatrix cov() const; 00029 //friend class chain; 00030 private: 00031 double _a; 00032 double _b; 00033 double _w; 00034 double _sav; 00035 double _ssav; 00036 double _zav; 00037 double _szav; 00038 double _zzav; 00039 double _chisq; 00040 double _sig_inv; 00041 double _c11; 00042 double _c21; 00043 double _c22; 00044 int _nc; 00045 }; 00046 00047 inline zav::zav() { 00048 _a = _b = _w = _sav = _ssav = _zav = _szav = _zzav = 0; 00049 _chisq = -1; 00050 _c22 = _c21 = _c11 = _sig_inv = 0; 00051 _nc = 0; 00052 } 00053 00054 inline void zav::clear(void){ 00055 _w = _sav = _ssav = _zav = _szav = _zzav = 0; 00056 _chisq = -1; 00057 _c22 = _c21 = _c11 = _sig_inv = 0; 00058 _nc = 0; 00059 } 00060 00061 inline void zav::set(const zav *c) { 00062 if(c) { 00063 _w = c->_w; 00064 _sav = c->_sav; 00065 _ssav = c->_ssav; 00066 _zav = c->_zav; 00067 _szav = c->_szav; 00068 _zzav = c->_zzav; 00069 _sig_inv = c->_sig_inv; 00070 _c11 = c->_c11; 00071 _c21 = c->_c21; 00072 _c22 = c->_c22; 00073 _nc = c->_nc; 00074 } else { 00075 _w = _sav = _ssav = _zav = _szav = _zzav = 00076 _sig_inv = _c11 = _c21 = _c22 = 0; 00077 _nc = 0; 00078 } 00079 _a = _b = 0; 00080 _chisq = -1; 00081 } 00082 00083 inline void zav::add( double s, double z, double w) { 00084 _w += w; 00085 double sw = s * w; 00086 _sav += sw; 00087 _ssav += sw * s; 00088 double zw = z * w; 00089 _zav += zw; 00090 _szav += zw * s; 00091 _zzav += zw * z; 00092 _chisq = -1; 00093 _nc++; 00094 } 00095 /* 00096 inline void zav::add(const zhit & h) { 00097 double errsq = h._z.error(); 00098 if (errsq==0) return; 00099 errsq *= errsq; 00100 double w = 1/errsq; 00101 _w += w; 00102 double s = (double)h._s; 00103 double sw = s * w; 00104 _sav += sw; 00105 _ssav += sw * s; 00106 double z = (double)h._z; 00107 double zw = z * w; 00108 _zav += zw; 00109 _szav += zw * s; 00110 _zzav += zw * z; 00111 _chisq = -1; 00112 _nc++; 00113 } 00114 */ 00115 inline double zav::calculate() { 00116 double sig = _ssav * _w - _sav * _sav; 00117 if (sig!=0) { 00118 _sig_inv = 1/sig; 00119 _a = ( _szav * _w - _sav * _zav ) * _sig_inv; 00120 _b = ( _ssav * _zav - _sav * _szav ) * _sig_inv; 00121 _chisq = _zzav - 2 * _a * _szav - 2 * _b * _zav + _a * _a * _ssav 00122 + _b * _b * _w + 2 * _a * _b * _sav; 00123 _c11 = _w * _sig_inv; 00124 _c21 = - _sav * _sig_inv; 00125 _c22 = _ssav * _sig_inv; 00126 } else { 00127 _sig_inv = 0; 00128 _c11 = _c21 = _c22 = 0; 00129 _chisq = -1; 00130 } 00131 if(_nc==2) { 00132 _chisq = 0; 00133 } 00134 return _chisq; 00135 } 00136 00137 inline std::ostream&operator<<(std::ostream&o,const zav&z) { 00138 o << " zav::w=" << z._w << " sav=" << z._sav << " zav=" << z._zav 00139 << " nc=" << z._nc << " chisq=" << z._chisq << " a=" << z._a 00140 << " b=" << z._b << " c11=" << z._c11 << " c21=" << z._c21 00141 << " c22=" << z._c22 << " sig_inv=" << z._sig_inv << std::endl; 00142 return o; 00143 } 00144 00145 /*inline HepSymMatrix zav::cov() const 00146 { 00147 #ifndef __GNUG__ 00148 HepSymMatrix vret(2); 00149 #endif 00150 // vret(1,1) = _ssav; 00151 // vret(2,1) = _sav; 00152 // vret(2,2) = _w; 00153 // int iret = invert(vret); 00154 // if (iret) cerr << "zav::cov(): cannot invert the matrix" << endl; 00155 vret(1,1) = _ssav; 00156 vret(2,1) = _sav; 00157 vret(2,2) = _w; 00158 return vret; 00159 }*/ 00160 #endif /* ZAV */ 00161