00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "VertexFit/KalmanVertexFit.h"
00011 #include "VertexFit/HTrackParameter.h"
00012 #include "VertexFit/WTrackParameter.h"
00013 #include "math.h"
00014
00015 using namespace std;
00016
00017 KalmanVertexFit * KalmanVertexFit::m_pointer = 0;
00018
00019 KalmanVertexFit * KalmanVertexFit::instance() {
00020 if(m_pointer) return m_pointer;
00021 m_pointer = new KalmanVertexFit();
00022 return m_pointer;
00023 }
00024
00025 KalmanVertexFit::~KalmanVertexFit() {;}
00026
00027 KalmanVertexFit::KalmanVertexFit() {;}
00028
00029 void KalmanVertexFit::init() {
00030 m_x = HepVector(3, 0);
00031 m_C0 = HepSymMatrix(3, 0);
00032 m_C = HepSymMatrix(3, 0);
00033 m_chisq = 0;
00034 m_ndof = -3;
00035
00036 m_flag.clear();
00037 m_p.clear();
00038 m_hTrkOrigin.clear();
00039 m_hTrkInfit.clear();
00040
00041
00042 m_G.clear();
00043 m_A.clear();
00044 m_B.clear();
00045 m_c.clear();
00046 m_W.clear();
00047 m_GB.clear();
00048 m_chiF.clear();
00049
00050
00051
00052
00053 m_chisqCutforTrack = 50;
00054 m_maxVertexIteration = 3;
00055 m_maxTrackIteration = 5;
00056 m_chisqCutforTrackIteration = 0.1;
00057 }
00058
00059 void KalmanVertexFit::initVertex(const VertexParameter vtx) {
00060 int ifail;
00061 m_x = vtx.Vx();
00062 m_C0 = (vtx.Evx()).inverse(ifail);
00063 m_C = m_C0;
00064 }
00065
00066 VertexParameter KalmanVertexFit::vtx() const {
00067 int ifail;
00068 VertexParameter v;
00069 v.setVx(m_x);
00070
00071 std::vector<int> trackid;
00072 v.setEvx(m_C.inverse(ifail));
00073 return v;
00074 }
00075
00076 HepSymMatrix KalmanVertexFit::Ex() const {
00077 int ifail;
00078 return m_C.inverse(ifail);
00079 }
00080
00081 void KalmanVertexFit::addTrack(const HTrackParameter htrk) {
00082
00083 int ifail;
00084
00085 HepVector p(3, 0);
00086 p = htrk.p();
00087
00088 m_p.push_back(p);
00089 m_flag.push_back(0);
00090 m_chiF.push_back(999.);
00091
00092 HepMatrix A(5, 3, 0);
00093 HepMatrix B(5, 3, 0);
00094 HepSymMatrix G(5, 0);
00095
00096 G = (htrk.eHel()).inverse(ifail);
00097
00098 m_hTrkOrigin.push_back(htrk);
00099
00100
00101 m_hTrkInfit.push_back(htrk);
00102
00103
00104 m_G.push_back(G);
00105 m_A.push_back(A);
00106 m_B.push_back(B);
00107
00108 HepSymMatrix W(3, 0);
00109 HepSymMatrix GB(5, 0);
00110 HepVector c(5, 0);
00111
00112 m_W.push_back(W);
00113 m_GB.push_back(GB);
00114 m_c.push_back(c);
00115
00116 }
00117
00118 int KalmanVertexFit::numTrack() const {
00119 int num = 0;
00120 for(int k = 0; k < m_hTrkOrigin.size(); k++)
00121 if(m_flag[k] == 0) num = num +1;
00122
00123 return num;
00124 }
00125
00126 std::vector<int> KalmanVertexFit::trackIDList() const {
00127 std::vector<int> trackid;
00128 for(int k = 0; k < m_hTrkOrigin.size(); k++)
00129 if(m_flag[k] == 0) trackid.push_back(trackID(k));
00130
00131 return trackid;
00132 }
00133
00134
00135 int KalmanVertexFit::filter(const int k) {
00136 int ifail;
00137 double chisq[10];
00139
00141 HepVector pp(3, 0);
00142 HepVector xp(3, 0);
00143 HepSymMatrix CP(3, 0);
00144 xp = m_x;
00145 CP = m_C;
00146 pp = m_p[k];
00147
00148 for(int iter = 0; iter < m_maxTrackIteration; iter++) {
00149
00150
00151
00152 updateMatrices(k, pp, xp);
00153
00154
00155
00156 HepSymMatrix CK(3, 0);
00157 CK = CP + (m_GB[k].similarity(m_A[k].T()));
00158
00159
00160
00161 HepVector x(3, 0);
00162 HepVector p(3, 0);
00163 x = CK.inverse(ifail) * (CP * xp + m_A[k].T() * m_GB[k] *
00164 (m_hTrkOrigin[k].hel() - m_c[k]));
00165 p = m_W[k] * m_B[k].T() * m_G[k] *(m_hTrkOrigin[k].hel()- m_c[k]-
00166 m_A[k]*x);
00167
00168 HepSymMatrix chif(1, 0);
00169 chif = CP.similarity((x-xp).T()) +
00170 m_G[k].similarity((m_hTrkOrigin[k].hel()-
00171 m_c[k]- m_A[k]*x - m_B[k]*p).T());
00172
00173
00174
00175
00176
00177 pp = p;
00178
00179 xp = x;
00180 CP = CK;
00181 chisq[iter] = chif[0][0];
00182 m_chiF[k] = chif[0][0];
00183
00184 if(iter > 0) {
00185 double delchi = (chisq[iter]-chisq[iter-1]);
00186 if(fabs(delchi) < m_chisqCutforTrackIteration) break;
00187 }
00188 }
00189
00190 if(m_chiF[k] > m_chisqCutforTrack) {
00191
00192 m_flag[k]=1;
00193 return 1;
00194 } else {
00195
00196 m_x = xp;
00197 m_C = CP;
00198 m_p[k] = pp;
00199 return 0;
00200 }
00201 }
00202
00203 void KalmanVertexFit::inverse(const int k) {
00204
00205
00206
00207
00208 int ifail;
00209 HepVector x(3, 0);
00210 HepVector p(3, 0);
00211 HepSymMatrix C(3, 0);
00212 C = m_C - m_GB[k].similarity(m_A[k].T());
00213 x = C.inverse(ifail)*(m_C*m_x-m_A[k].T() * m_GB[k] *
00214 (m_hTrkOrigin[k].hel()-m_c[k]));
00215
00216
00217
00218 m_x = x;
00219 m_C = C;
00220 }
00221 void KalmanVertexFit::filter() {
00222
00223 for(int iter = 0; iter < m_maxVertexIteration; iter++) {
00224
00225
00226
00227
00228
00229 m_C = m_C0;
00230 m_chisq = 0;
00231 m_ndof = -3;
00232
00233 for(int k = 0; k < m_hTrkOrigin.size(); k++) {
00234 if(m_flag[k]==1) continue;
00235 if(filter(k)==0)
00236 m_ndof += 2;
00237
00238 }
00239
00240
00241
00242
00243
00244
00245 for(int k = 0; k < m_hTrkOrigin.size(); k++) {
00246 if(m_flag[k]==1) continue;
00247 double chi2 = chiS(k);
00248 if(chi2 < 0) {
00249 remove(k);
00250 continue;
00251 }
00252
00253 if(chi2 > m_chisqCutforTrack) {
00254 remove(k);
00255 continue;
00256 }
00257 m_chisq += chi2;
00258
00259 }
00260 }
00261 }
00262
00263 double KalmanVertexFit::chiS(const int k) const {
00264 if(m_flag[k]==1) return 999;
00265 int ifail;
00266 HepVector x(3, 0);
00267 HepVector p(3, 0);
00268 HepSymMatrix C(3, 0);
00269 C = m_C - m_GB[k].similarity(m_A[k].T());
00270 x = C.inverse(ifail)*(m_C*m_x-m_A[k].T() * m_GB[k] *
00271 (m_hTrkOrigin[k].hel()-m_c[k]));
00272 p = m_W[k]* m_B[k].T()*m_G[k]*
00273 (m_hTrkOrigin[k].hel()-m_c[k]-m_A[k]*m_x);
00274 HepSymMatrix chis(1, 0);
00275 chis = C.similarity((x-m_x).T())+
00276 m_G[k].similarity((m_hTrkOrigin[k].hel()-m_c[k]-
00277 m_A[k]*m_x-m_B[k]*p).T());
00278 return chis[0][0];
00279 }
00280
00281 void KalmanVertexFit::smooth(const int k) {
00282
00283
00284
00285
00286 if(m_flag[k] == 1) return;
00287 int ifail;
00288
00289 HepVector pp(3, 0);
00290 HepVector xp(3, 0);
00291 HepSymMatrix CP(3, 0);
00292 xp = m_x;
00293 CP = m_C;
00294 pp = m_p[k];
00295
00296 updateMatrices(k, pp, xp);
00297 pp = m_W[k] * m_B[k].T() * m_G[k] * (m_hTrkOrigin[k].hel() - m_c[k] - m_A[k] * xp);
00298
00299 updateMatrices(k, pp, xp);
00300
00301 HepVector helix(5, 0);
00302 helix = m_c[k] + m_A[k] * xp + m_B[k] * pp;
00303 HepMatrix E(3, 3, 0);
00304 E = -CP.inverse(ifail) * m_A[k].T() * m_G[k] * m_B[k] * m_W[k];
00305 HepSymMatrix D(3, 0);
00306 D = m_W[k] + CP.similarity(E.T());
00307
00308 HepMatrix middle(5, 5, 0);
00309 HepSymMatrix error(5, 0);
00310 middle = (CP.inverse(ifail)).similarity(m_A[k]) + m_A[k] * E * m_B[k].T() +
00311 (m_A[k] * E * m_B[k].T()).T() + D.similarity(m_B[k]);
00312 error.assign(middle);
00313
00314 m_hTrkInfit[k].setHel(helix);
00315 m_hTrkInfit[k].setEHel(error);
00316
00317 m_p[k] = pp;
00318
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329
00330
00331
00332
00333
00334
00335
00336
00337
00338
00339
00340
00341
00342
00343
00344
00345
00346
00347
00348
00349
00350
00351
00352
00353 }
00354
00355 void KalmanVertexFit::remove(const int k) {
00356
00357
00358
00359
00360
00361
00362
00363
00364
00365
00366
00367
00368
00369
00370
00371
00372 if(m_flag[k]==1) return;
00373 inverse(k);
00374 m_ndof -= 2;
00375 m_flag[k]=1;
00376
00377
00378 }
00379
00380
00381 HTrackParameter KalmanVertexFit::hTrack(const int k) const{
00382
00383
00384
00385 return m_hTrkInfit[k];
00386 }
00387
00388 WTrackParameter KalmanVertexFit::wTrack(const int k, const double mass) const {
00389 if(m_flag[k] == 1) return m_hTrkOrigin[k].wTrack();
00390 int ifail;
00391
00392 HepVector pp(3, 0);
00393 HepVector xp(3, 0);
00394 HepSymMatrix CP(3, 0);
00395 xp = m_x;
00396 CP = m_C;
00397 pp = m_p[k];
00398
00399 WTrackParameter wtrk;
00400 double Energy = sqrt(pp[0]*pp[0] + pp[1]*pp[1] + pp[2]*pp[2] + mass * mass);
00401
00402 HepMatrix Awtrk(7, 3, 0), Bwtrk(7, 3, 0);
00403 Awtrk[4][0] = Awtrk[5][1] = Awtrk[6][2] = 1;
00404 Bwtrk[0][0] = Bwtrk[1][1] = Bwtrk[2][2] = 1;
00405 Bwtrk[3][0] = pp[0] / Energy;
00406 Bwtrk[3][1] = pp[1] / Energy;
00407 Bwtrk[3][2] = pp[2] / Energy;
00408
00409 HepVector w(7, 0);
00410 HepSymMatrix Ew(7, 0);
00411 w[0] = pp[0]; w[1] = pp[1]; w[2] = pp[2]; w[3] = Energy;
00412 w[4] = xp[0]; w[5] = xp[1]; w[6] = xp[2];
00413
00414 HepSymMatrix Gwtrk(7, 0);
00415 Gwtrk = m_hTrkOrigin[k].wTrack().Ew().inverse(ifail);
00416 HepSymMatrix Wwtrk(3, 0);
00417 Wwtrk = Gwtrk.similarity(Bwtrk.T()).inverse(ifail);
00418
00419 HepMatrix Ewtrk(3, 3, 0);
00420 Ewtrk = -CP.inverse(ifail) * Awtrk.T() * Gwtrk * Bwtrk * Wwtrk;
00421 HepSymMatrix Dwtrk(3, 0);
00422 Dwtrk = Wwtrk + CP.similarity(Ewtrk.T());
00423
00424 HepMatrix Ewmiddle(7, 7, 0);
00425 Ewmiddle = (CP.inverse(ifail)).similarity(Awtrk) + Awtrk * Ewtrk * Bwtrk.T() +
00426 (Awtrk * Ewtrk * Bwtrk.T()).T() + Dwtrk.similarity(Bwtrk);
00427 Ew.assign(Ewmiddle);
00428
00429 wtrk.setCharge(m_hTrkOrigin[k].charge());
00430 wtrk.setW(w);
00431 wtrk.setEw(Ew);
00432
00433 return wtrk;
00434 }
00435
00436
00437
00438
00439 void KalmanVertexFit::updateMatrices(const int k, const HepVector p, const HepVector x) {
00440
00441 int ifail;
00442
00443
00444
00445 HTrackParameter he(m_hTrkOrigin[k].charge(), p, x);
00446
00447
00448
00449 m_A[k] = he.dHdx(p, x);
00450 m_B[k] = he.dHdp(p, x);
00451
00452
00453 m_W[k] = (m_G[k].similarity(m_B[k].T())).inverse(ifail);
00454 m_GB[k] = m_G[k] - (m_W[k].similarity(m_B[k])).similarity(m_G[k]);
00455 m_c[k] = he.hel() - m_A[k] * x - m_B[k] * p;
00456 }
00457
00458 void KalmanVertexFit::updateMatrices(const int k) {
00459
00460 int ifail;
00461 HepVector p(3, 0);
00462 HepVector x(3, 0);
00463 p = m_p[k];
00464 x = m_x;
00465 updateMatrices(k, p, x);
00466 }
00467
00468 HepVector KalmanVertexFit::pull(const int k) {
00469 HepVector n_pull(5, 0);
00470 n_pull[0] = n_pull[1] = n_pull[2] = n_pull[3] = n_pull[4] = 999.;
00471 if(m_flag[k] == 1) return n_pull;
00472 for (int i = 0 ; i < 5; i++) {
00473 double cov = m_hTrkOrigin[k].eHel()[i][i] - m_hTrkInfit[k].eHel()[i][i];
00474 if (cov == 0.) continue;
00475 n_pull[i] = (m_hTrkOrigin[k].hel()[i] - m_hTrkInfit[k].hel()[i]) / sqrt(cov);
00476 }
00477 return n_pull;
00478 }
00479
00480 double KalmanVertexFit::pullmomentum(const int k) {
00481 double pull = 999.;
00482 if(m_flag[k] == 1) return pull;
00483
00484 double kappa_origin = m_hTrkOrigin[k].kappa();
00485 double kappa_infit = m_hTrkInfit[k].kappa();
00486 double lamb_origin = m_hTrkOrigin[k].lambda();
00487 double lamb_infit = m_hTrkInfit[k].lambda();
00488
00489 double Vkappa_origin = m_hTrkOrigin[k].eHel()[2][2];
00490 double Vkappa_infit = m_hTrkInfit[k].eHel()[2][2];
00491 double Vlamb_origin = m_hTrkOrigin[k].eHel()[4][4];
00492 double Vlamb_infit = m_hTrkInfit[k].eHel()[4][4];
00493 double V_kappa_lamb_origin = m_hTrkOrigin[k].eHel()[4][2];
00494 double V_kappa_lamb_infit = m_hTrkInfit[k].eHel()[4][2];
00495
00496 double P_origin = calculationP(kappa_origin, lamb_origin);
00497
00498
00499
00500
00501
00502 double P_infit = calculationP(kappa_infit, lamb_infit);
00503
00504 double SigmaP_origin = calculationSigmaP(kappa_origin, lamb_origin, Vkappa_origin,
00505 Vlamb_origin, V_kappa_lamb_origin);
00506 double SigmaP_infit = calculationSigmaP(kappa_infit, lamb_infit, Vkappa_infit,
00507 Vlamb_infit, V_kappa_lamb_infit);
00508 if ((SigmaP_origin - SigmaP_infit) <= 0.) return pull;
00509 pull = (P_origin - P_infit) / sqrt (SigmaP_origin - SigmaP_infit);
00510
00511 return pull;
00512 }
00513
00514
00515 double KalmanVertexFit::calculationP(const double kappa, const double lamb) {
00516 double P = 0.0;
00517 P = sqrt(1 + lamb * lamb) / kappa;
00518
00519 return P;
00520 }
00521
00522 double KalmanVertexFit::calculationSigmaP(const double kappa, const double lamb,
00523 const double Vkappa, const double Vlamb,
00524 const double Vkappa_lamb) {
00525 double SigmaP = 0.0;
00526 double dp_dkappa = -sqrt(1 + lamb*lamb) /kappa/kappa;
00527 double dp_dlamb = lamb / kappa / sqrt(1 + lamb*lamb);
00528 SigmaP = dp_dkappa * dp_dkappa * Vkappa + dp_dlamb * dp_dlamb * Vlamb
00529 + 2 * dp_dkappa * dp_dlamb * Vkappa_lamb;
00530 return SigmaP;
00531 }