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