00001 #ifndef VertexFit_KalmanKinematicFit_H
00002 #define VertexFit_KalmanKinematicFit_H
00003
00004 #include <vector>
00005 #include "VertexFit/WTrackParameter.h"
00006 #include "VertexFit/KinematicConstraints.h"
00007 #include "VertexFit/TrackPool.h"
00008 #include "VertexFit/GammaShape.h"
00009 #include "TGraph2D.h"
00010
00011 class KalmanKinematicFit : public TrackPool{
00012
00013 public:
00014
00015
00016 static KalmanKinematicFit * instance();
00017 ~KalmanKinematicFit();
00018
00019
00020
00021 void AddResonance(int number, double mres, std::vector<int> tlis);
00022 void AddResonance(int number, double mres, int n1);
00023 void AddResonance(int number, double mres, int n1, int n2);
00024 void AddResonance(int number, double mres, int n1, int n2, int n3);
00025 void AddResonance(int number, double mres, int n1, int n2, int n3, int n4);
00026 void AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
00027 int n5);
00028 void AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
00029 int n5, int n6);
00030 void AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
00031 int n5, int n6, int n7);
00032 void AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
00033 int n5, int n6, int n7, int n8);
00034 void AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
00035 int n5, int n6, int n7, int n8, int n9);
00036 void AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
00037 int n5, int n6, int n7, int n8, int n9,
00038 int n10);
00039 void AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
00040 int n5, int n6, int n7, int n8, int n9,
00041 int n10, int n11);
00042 void AddResonance(int number, double mres, int n1, int n2, int n3, int n4,
00043 int n5, int n6, int n7, int n8, int n9,
00044 int n10, int n11, int n12);
00045
00046
00047
00048 void AddTotalEnergy(int number, double etot, std::vector<int> lis);
00049 void AddTotalEnergy(int number, double etot, int n1);
00050 void AddTotalEnergy(int number, double etot, int n1, int n2);
00051 void AddTotalEnergy(int number, double etot, int n1, int n2, int n3);
00052 void AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4);
00053 void AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
00054 int n5);
00055 void AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
00056 int n5, int n6);
00057 void AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
00058 int n5, int n6, int n7);
00059 void AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
00060 int n5, int n6, int n7, int n8);
00061 void AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
00062 int n5, int n6, int n7, int n8, int n9);
00063 void AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
00064 int n5, int n6, int n7, int n8, int n9, int n10);
00065 void AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
00066 int n5, int n6, int n7, int n8, int n9, int n10, int n11);
00067 void AddTotalEnergy(int number, double etot, int n1, int n2, int n3, int n4,
00068 int n5, int n6, int n7, int n8, int n9,
00069 int n10, int n11, int n12);
00070
00071
00072
00073
00074 void AddTotalMomentum(int number, double ptot, std::vector<int> lis);
00075 void AddTotalMomentum(int number, double ptot, int n1);
00076 void AddTotalMomentum(int number, double ptot, int n1, int n2);
00077 void AddTotalMomentum(int number, double ptot, int n1, int n2, int n3);
00078 void AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4);
00079 void AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
00080 int n5);
00081 void AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
00082 int n5, int n6);
00083 void AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
00084 int n5, int n6, int n7);
00085 void AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
00086 int n5, int n6, int n7, int n8);
00087 void AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
00088 int n5, int n6, int n7, int n8, int n9);
00089 void AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
00090 int n5, int n6, int n7, int n8, int n9,
00091 int n10);
00092 void AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
00093 int n5, int n6, int n7, int n8, int n9,
00094 int n10, int n11);
00095 void AddTotalMomentum(int number, double ptot, int n1, int n2, int n3, int n4,
00096 int n5, int n6, int n7, int n8, int n9,
00097 int n10, int n11, int n12);
00098
00099
00100
00101 void AddThreeMomentum(int number, Hep3Vector p3);
00102
00103
00104
00105 void AddFourMomentum(int number, HepLorentzVector p4);
00106 void AddFourMomentum(int number, double etot);
00107
00108
00109
00110
00111
00112 void AddEqualMass(int number, std::vector<int> tlis1, std::vector<int> tlis2);
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122 void BuildVirtualParticle(int number);
00123
00124
00125
00126 void init();
00127
00128
00129
00130 void setFlag(const bool flag = 1) {m_flag = flag;}
00131 void setIterNumber(const int niter = 5) {m_niter = niter;}
00132 void setChisqCut(const double chicut = 200, const double chiter=0.05) {m_chicut = chicut;m_chiter=chiter;}
00133
00134
00135
00136 void setEspread (const double espread = 0.0009) {m_espread = espread;}
00137 void setCollideangle (const double collideangle = 11e-3) {m_collideangle = collideangle;}
00138 void setDynamicerror (const bool dynamicerror = 1) {m_dynamicerror = dynamicerror;}
00139 void setTgraph ( TGraph2D* graph2d) {m_graph2d = graph2d;}
00140
00141
00142
00143 bool Fit();
00144 bool Fit(int n);
00145
00146
00147
00148
00149
00150 double chisq() const {return m_chi;}
00151 double chisq(int n) const {return m_chisq[n];}
00152
00153
00154 HepLorentzVector pfit(int n) const {return p4Infit(n);}
00155
00156
00157 HepLorentzVector pfit1(int n) {return p4Origin(n);}
00158 HepVector xfit() {return m_q.sub(1,3);}
00159
00160 WTrackParameter origin(int n) const {return wTrackOrigin(n);}
00161 WTrackParameter infit(int n) const {return wTrackInfit(n);}
00162
00163 HepVector pull(int n) ;
00164
00165
00166 double espread() const {return m_espread;}
00167 double collideangle() const {return m_collideangle;}
00168 bool dynamicerror() const {return m_dynamicerror;}
00169
00170 HepVector cpu() const {return m_cpu;}
00171 HepSymMatrix getCOrigin(int i) const;
00172 HepSymMatrix getCInfit(int i) const;
00173
00174 WTrackParameter wVirtualTrack(int n) const {return m_virtual_wtrk[n];}
00175
00176 private:
00177
00178 std::vector<WTrackParameter> m_virtual_wtrk;
00179
00180 private:
00181 void updateConstraints(KinematicConstraints kc);
00182 void fit();
00183
00184 void upCovmtx();
00185 void upTrkpar();
00186 void clearDDT();
00187 void fit(int n);
00188 void covMatrix(int n);
00189 void gda();
00190 private:
00191 std::vector<KinematicConstraints> m_kc;
00192 std::vector<double> m_chisq;
00193 double m_chi;
00194 HepSymMatrix m_Vm;
00195 HepMatrix m_A;
00196 void setA(int ic, int itk, const HepMatrix &p) {m_A.sub(ic+1, itk, p);}
00197 HepMatrix m_AT;
00198 void setAT(int itk, int ic, const HepMatrix &p) { m_AT.sub(itk, ic+1, p);}
00199 HepVector m_G;
00200 HepSymMatrix m_W;
00201 HepMatrix m_KP;
00202
00203
00204 HepMatrix m_B;
00205 void setB(int ic, int itk, const HepMatrix &p) {m_B.sub(ic+1, itk, p);}
00206 HepMatrix m_BT;
00207 void setBT(int itk, int ic, const HepMatrix &p) { m_BT.sub(itk, ic+1, p);}
00208
00209
00210 HepMatrix m_KQ;
00211 int m_nc;
00212 int m_nktrk;
00213
00214
00215 HepVector m_p0;
00216 HepVector m_p;
00217 HepSymMatrix m_C0;
00218 HepSymMatrix m_C;
00219 HepVector pOrigin(int i) const ;
00220 HepLorentzVector p4Origin(int i) const { HepVector p(4, 0); p = pOrigin(i); return HepLorentzVector(p[0], p[1], p[2], p[3]);}
00221 HepVector pInfit(int i) const ;
00222 HepLorentzVector p4Infit(int i) const { HepVector p(4, 0); p = pInfit(i); return HepLorentzVector(p[0], p[1], p[2], p[3]); }
00223
00224
00225 void setPOrigin(int i, const HepVector &p) { m_p0.sub(i, p);}
00226 void setPInfit(int i, const HepVector &p) {m_p.sub(i, p);}
00227 void setCOrigin(int i, const HepSymMatrix &D) {m_C0.sub(i, D);}
00228 void setCInfit(int i, const HepSymMatrix &D) {m_C.sub(i,D);}
00229
00230 HepVector m_q0;
00231 HepVector m_q;
00232 HepSymMatrix m_D0;
00233 HepSymMatrix m_D;
00234 HepSymMatrix m_D0inv;
00235 HepSymMatrix m_Dinv;
00236
00237 void setQOrigin(int i, const HepVector &q) { m_q0.sub(i, q);}
00238 void setQInfit(int i, const HepVector &q) {m_q.sub(i, q);}
00239 void setDOrigin(int i, const HepSymMatrix &D) {m_D0.sub(i, D);}
00240 void setDInfit(int i, const HepSymMatrix &D) {m_D.sub(i,D);}
00241 void setDOriginInv(int i, const HepSymMatrix &Dinv) {m_D0inv.sub(i,Dinv);}
00242
00243
00244
00245 private:
00246 KalmanKinematicFit();
00247 static KalmanKinematicFit * m_pointer;
00248 private:
00249 int m_niter;
00250 bool m_flag;
00251 double m_chicut;
00252 double m_chiter;
00253 private:
00254 double m_espread;
00255 double m_collideangle;
00256 private:
00257 HepVector m_cpu;
00258 private:
00259 bool m_dynamicerror;
00260 private:
00261 static const int NTRKPAR;
00262 static const int NKFPAR;
00263 static const int Resonance;
00264 static const int TotalEnergy;
00265 static const int TotalMomentum;
00266 static const int ThreeMomentum;
00267 static const int FourMomentum;
00268 static const int EqualMass;
00269 static const int Position;
00270 private:
00271 TGraph2D* m_graph2d;
00272
00273
00274 };
00275 #endif
00276