Belle II Software development
FourCFitKFit.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9#include <analysis/VertexFitting/KFit/FourCFitKFit.h>
10
11#include <analysis/dataobjects/Particle.h>
12#include <analysis/VertexFitting/KFit/MakeMotherKFit.h>
13#include <analysis/utility/CLHEPToROOT.h>
14#include <framework/gearbox/Const.h>
15
16#include <cstdio>
17
18#include <TMath.h>
19
20using namespace std;
21using namespace Belle2;
22using namespace Belle2::analysis;
23using namespace CLHEP;
24using namespace ROOT::Math;
25
27{
28 m_FlagFitted = false;
31 m_FlagAtDecayPoint = true;
33 m_d = HepMatrix(4, 1, 0);
34 m_V_D = HepMatrix(4, 4, 0);
35 m_lam = HepMatrix(4, 1, 0);
36 m_AfterVertexError = HepSymMatrix(3, 0);
37 m_InvariantMass = -1.0;
38 m_FourMomentum = PxPyPzEVector();
39}
40
41
43
44
46FourCFitKFit::setVertex(const HepPoint3D& v) {
48
50}
51
52
54FourCFitKFit::setVertexError(const HepSymMatrix& e) {
55 if (e.num_row() != 3)
56 {
58 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
59 return m_ErrorCode;
60 }
61
64
66}
67
68
75
76
78FourCFitKFit::setFourMomentum(const PxPyPzEVector& m) {
80
82}
83
84
91
92
95 m_IsFixMass.push_back(true);
96
98}
99
100
103 m_IsFixMass.push_back(false);
104
106}
107
108
111 if (e.num_row() != 3 || e.num_col() != KFitConst::kNumber7)
112 {
114 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
115 return m_ErrorCode;
116 }
117
118 m_BeforeTrackVertexError.push_back(e);
121
123}
124
125
128 HepMatrix zero(3, KFitConst::kNumber7, 0);
129
130 return this->setTrackVertexError(zero);
131}
132
133
135FourCFitKFit::setCorrelation(const HepMatrix& m) {
136 return KFitBase::setCorrelation(m);
137}
138
139
144
145
146const HepPoint3D
147FourCFitKFit::getVertex(const int flag) const
148{
149 if (flag == KFitConst::kAfterFit && !isFitted()) return HepPoint3D();
150
151 switch (flag) {
153 return m_BeforeVertex;
154
156 return m_AfterVertex;
157
158 default:
159 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
160 return HepPoint3D();
161 }
162}
163
164
165const HepSymMatrix
166FourCFitKFit::getVertexError(const int flag) const
167{
168 if (flag == KFitConst::kAfterFit && !isFitted()) return HepSymMatrix(3, 0);
169
170 if (flag == KFitConst::kBeforeFit)
171 return m_BeforeVertexError;
173 return m_AfterVertexError;
174 else {
175 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
176 return HepSymMatrix(3, 0);
177 }
178}
179
180
181double
186
187
188bool
193
194
195bool
200
201
202double
204{
205 return m_CHIsq;
206}
207
208
209const HepMatrix
210FourCFitKFit::getTrackVertexError(const int id, const int flag) const
211{
212 if (flag == KFitConst::kAfterFit && !isFitted()) return HepMatrix(3, KFitConst::kNumber7, 0);
213 if (!isTrackIDInRange(id)) return HepMatrix(3, KFitConst::kNumber7, 0);
214
215 if (flag == KFitConst::kBeforeFit)
216 return m_BeforeTrackVertexError[id];
218 return m_AfterTrackVertexError[id];
219 else {
220 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
221 return HepMatrix(3, KFitConst::kNumber7, 0);
222 }
223}
224
225
226double
228{
229 if (!isFitted()) return -1;
230 if (!isTrackIDInRange(id)) return -1;
231
232 if (m_IsFixMass[id]) {
233
234 HepMatrix da(m_Tracks[id].getFitParameter(KFitConst::kBeforeFit) - m_Tracks[id].getFitParameter(KFitConst::kAfterFit));
235 int err_inverse = 0;
236 const double chisq = (da.T() * (m_Tracks[id].getFitError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
237
238 if (err_inverse) {
239 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
240 return -1;
241 }
242
243 return chisq;
244
245 } else {
246
247 HepMatrix da(m_Tracks[id].getMomPos(KFitConst::kBeforeFit) - m_Tracks[id].getMomPos(KFitConst::kAfterFit));
248 int err_inverse = 0;
249 const double chisq = (da.T() * (m_Tracks[id].getError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
250
251 if (err_inverse) {
252 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
253 return -1;
254 }
255
256 return chisq;
257
258 }
259}
260
261
262const HepMatrix
263FourCFitKFit::getCorrelation(const int id1, const int id2, const int flag) const
264{
265 if (flag == KFitConst::kAfterFit && !isFitted()) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
266 if (!isTrackIDInRange(id1)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
267 if (!isTrackIDInRange(id2)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
268
269 switch (flag) {
271 return KFitBase::getCorrelation(id1, id2, flag);
272
274 return makeError3(
275 this->getTrackMomentum(id1),
276 this->getTrackMomentum(id2),
277 m_V_al_1.sub(KFitConst::kNumber7 * id1 + 1, KFitConst::kNumber7 * (id1 + 1), KFitConst::kNumber7 * id2 + 1,
278 KFitConst::kNumber7 * (id2 + 1)),
279 m_IsFixMass[id1],
280 m_IsFixMass[id2]);
281
282 default:
283 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
284 return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
285 }
286}
287
288
293
294
298 {
300 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
301 return m_ErrorCode;
302 }
303
304
305 if (m_IsFixMass.size() == 0)
306 {
307 // If no fix_mass flag at all,
308 // all tracks are considered to be fixed at mass.
309 for (int i = 0; i < m_TrackCount; i++) this->fixMass();
310 } else if (m_IsFixMass.size() != (unsigned int)m_TrackCount)
311 {
313 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
314 return m_ErrorCode;
315 }
316
317
319 {
320 int index = 0;
321 m_al_0 = HepMatrix(KFitConst::kNumber7 * m_TrackCount, 1, 0);
322 m_property = HepMatrix(m_TrackCount, 3, 0);
323 m_V_al_0 = HepSymMatrix(KFitConst::kNumber7 * m_TrackCount, 0);
324
325 for (auto& track : m_Tracks) {
326 // momentum x,y,z and position x,y,z
327 m_al_0[index * KFitConst::kNumber7 + 0][0] = track.getMomentum(KFitConst::kBeforeFit).x();
328 m_al_0[index * KFitConst::kNumber7 + 1][0] = track.getMomentum(KFitConst::kBeforeFit).y();
329 m_al_0[index * KFitConst::kNumber7 + 2][0] = track.getMomentum(KFitConst::kBeforeFit).z();
330 m_al_0[index * KFitConst::kNumber7 + 3][0] = track.getMomentum(KFitConst::kBeforeFit).t();
331 m_al_0[index * KFitConst::kNumber7 + 4][0] = track.getPosition(KFitConst::kBeforeFit).x();
332 m_al_0[index * KFitConst::kNumber7 + 5][0] = track.getPosition(KFitConst::kBeforeFit).y();
333 m_al_0[index * KFitConst::kNumber7 + 6][0] = track.getPosition(KFitConst::kBeforeFit).z();
334 // these error
335 m_V_al_0.sub(index * KFitConst::kNumber7 + 1, track.getError(KFitConst::kBeforeFit));
336 // charge, mass, a
337 m_property[index][0] = track.getCharge();
338 m_property[index][1] = track.getMass();
339 const double c = Const::speedOfLight * 1e-4;
340 m_property[index][2] = -c * m_MagneticField * track.getCharge();
341 index++;
342 }
343
344 // error between track and track
345 if (m_FlagCorrelation) {
346 this->prepareCorrelation();
348 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
349 return m_ErrorCode;
350 }
351 }
352
353 // set member matrix
354 m_al_1 = m_al_0;
355
356 // define size of matrix
359
360 } else
361 {
362 // m_FlagFitIncludingVertex == true
363 int index = 0;
364 m_al_0 = HepMatrix(KFitConst::kNumber7 * m_TrackCount + 3, 1, 0);
365 m_property = HepMatrix(m_TrackCount, 3, 0);
366 m_V_al_0 = HepSymMatrix(KFitConst::kNumber7 * m_TrackCount + 3, 0);
367
368 for (auto& track : m_Tracks) {
369 // momentum x,y,z and position x,y,z
370 m_al_0[index * KFitConst::kNumber7 + 0][0] = track.getMomentum(KFitConst::kBeforeFit).x();
371 m_al_0[index * KFitConst::kNumber7 + 1][0] = track.getMomentum(KFitConst::kBeforeFit).y();
372 m_al_0[index * KFitConst::kNumber7 + 2][0] = track.getMomentum(KFitConst::kBeforeFit).z();
373 m_al_0[index * KFitConst::kNumber7 + 3][0] = track.getMomentum(KFitConst::kBeforeFit).t();
374 m_al_0[index * KFitConst::kNumber7 + 4][0] = track.getPosition(KFitConst::kBeforeFit).x();
375 m_al_0[index * KFitConst::kNumber7 + 5][0] = track.getPosition(KFitConst::kBeforeFit).y();
376 m_al_0[index * KFitConst::kNumber7 + 6][0] = track.getPosition(KFitConst::kBeforeFit).z();
377 // these error
378 m_V_al_0.sub(index * KFitConst::kNumber7 + 1, track.getError(KFitConst::kBeforeFit));
379 // charge, mass, a
380 m_property[index][0] = track.getCharge();
381 m_property[index][1] = track.getMass();
382 const double c = Const::speedOfLight * 1e-4;
383 m_property[index][2] = -c * m_MagneticField * track.getCharge();
384 index++;
385 }
386
387 // vertex
392
393 // error between track and track
394 if (m_FlagCorrelation) {
395 this->prepareCorrelation();
397 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
398 return m_ErrorCode;
399 }
400 }
401
402 // set member matrix
403 m_al_1 = m_al_0;
404
405 // define size of matrix
407 m_D = m_V_al_1.sub(1, 4, 1, KFitConst::kNumber7 * m_TrackCount + 3);
408 }
409
411}
412
413
416 char buf[1024];
417 sprintf(buf, "%s:%s(): internal error; this function should never be called", __FILE__, __func__);
418 B2FATAL(buf);
419
420 /* NEVER REACHEd HERE */
422}
423
424
427 if (m_BeforeCorrelation.size() != static_cast<unsigned int>(m_TrackCount * (m_TrackCount - 1) / 2))
428 {
430 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
431 return m_ErrorCode;
432 }
433
434 int row = 0, col = 0;
435
436 for (auto& hm : m_BeforeCorrelation)
437 {
438 // counter
439 row++;
440 if (row == m_TrackCount) {
441 col++;
442 row = col + 1;
443 }
444
445 int ii = 0, jj = 0;
446 for (int i = KFitConst::kNumber7 * row; i < KFitConst::kNumber7 * (row + 1); i++) {
447 for (int j = KFitConst::kNumber7 * col; j < KFitConst::kNumber7 * (col + 1); j++) {
448 m_V_al_0[i][j] = hm[ii][jj];
449 jj++;
450 }
451 jj = 0;
452 ii++;
453 }
454 }
455
457 {
458 // ...error of vertex
460
461 // ...error matrix between vertex and tracks
463 if (m_BeforeTrackVertexError.size() != (unsigned int)m_TrackCount) {
465 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
466 return m_ErrorCode;
467 }
468
469 int i = 0;
470 for (auto& hm : m_BeforeTrackVertexError) {
471 for (int j = 0; j < 3; j++) for (int k = 0; k < KFitConst::kNumber7; k++) {
473 }
474 i++;
475 }
476 }
477 }
478
480}
481
482
485 Hep3Vector h3v;
486 int index = 0;
487 for (auto& pdata : m_Tracks)
488 {
489 // tracks
490 // momentum
491 h3v.setX(m_al_1[index * KFitConst::kNumber7 + 0][0]);
492 h3v.setY(m_al_1[index * KFitConst::kNumber7 + 1][0]);
493 h3v.setZ(m_al_1[index * KFitConst::kNumber7 + 2][0]);
494 if (m_IsFixMass[index])
495 pdata.setMomentum(HepLorentzVector(h3v, sqrt(h3v.mag2() + pdata.getMass()*pdata.getMass())), KFitConst::kAfterFit);
496 else
497 pdata.setMomentum(HepLorentzVector(h3v, m_al_1[index * KFitConst::kNumber7 + 3][0]), KFitConst::kAfterFit);
498 // position
499 pdata.setPosition(HepPoint3D(
500 m_al_1[index * KFitConst::kNumber7 + 4][0],
501 m_al_1[index * KFitConst::kNumber7 + 5][0],
503 // error of the tracks
504 pdata.setError(this->makeError3(pdata.getMomentum(),
505 m_V_al_1.sub(
506 index * KFitConst::kNumber7 + 1,
507 (index + 1)*KFitConst::kNumber7,
508 index * KFitConst::kNumber7 + 1,
509 (index + 1)*KFitConst::kNumber7), m_IsFixMass[index]),
511 if (m_ErrorCode != KFitError::kNoError) break;
512 index++;
513 }
514
516 {
517 // vertex
521 // error of the vertex
522 for (int i = 0; i < 3; i++) for (int j = i; j < 3; j++) {
524 }
525 // error between vertex and tracks
526 for (int i = 0; i < m_TrackCount; i++) {
527 HepMatrix hm(3, KFitConst::kNumber7, 0);
528 for (int j = 0; j < 3; j++) for (int k = 0; k < KFitConst::kNumber7; k++) {
530 }
531 if (m_IsFixMass[i])
532 m_AfterTrackVertexError.push_back(this->makeError4(m_Tracks[i].getMomentum(), hm));
533 else
534 m_AfterTrackVertexError.push_back(hm);
535 }
536 } else
537 {
538 // not fit
540 }
541
543}
544
545
549 {
550
551 HepMatrix al_1_prime(m_al_1);
552 HepMatrix Sum_al_1(4, 1, 0);
553 std::vector<double> energy(m_TrackCount);
554 double a;
555
556 for (int i = 0; i < m_TrackCount; i++) {
557 a = m_property[i][2];
558 if (!m_FlagAtDecayPoint) a = 0.;
559 al_1_prime[i * KFitConst::kNumber7 + 0][0] -= a * (m_BeforeVertex.y() - al_1_prime[i * KFitConst::kNumber7 + 5][0]);
560 al_1_prime[i * KFitConst::kNumber7 + 1][0] += a * (m_BeforeVertex.x() - al_1_prime[i * KFitConst::kNumber7 + 4][0]);
561 energy[i] = sqrt(al_1_prime[i * KFitConst::kNumber7 + 0][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] +
562 al_1_prime[i * KFitConst::kNumber7 + 1][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] +
563 al_1_prime[i * KFitConst::kNumber7 + 2][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] +
564 m_property[i][1] * m_property[i][1]);
565 if (m_IsFixMass[i])
566 Sum_al_1[3][0] += energy[i];
567 else
568 Sum_al_1[3][0] += al_1_prime[i * KFitConst::kNumber7 + 3][0];
569 }
570
571 for (int i = 0; i < m_TrackCount; i++) {
572 for (int j = 0; j < 3; j++) Sum_al_1[j][0] += al_1_prime[i * KFitConst::kNumber7 + j][0];
573 }
574
575 m_d[0][0] = Sum_al_1[0][0] - m_FourMomentum.Px();
576 m_d[1][0] = Sum_al_1[1][0] - m_FourMomentum.Py();
577 m_d[2][0] = Sum_al_1[2][0] - m_FourMomentum.Pz();
578 m_d[3][0] = Sum_al_1[3][0] - m_FourMomentum.E();
579
580 for (int i = 0; i < m_TrackCount; i++) {
581 if (energy[i] == 0) {
583 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
584 break;
585 }
586
587 a = m_property[i][2];
588 if (!m_FlagAtDecayPoint) a = 0.;
589
590 if (m_IsFixMass[i]) {
591 double invE = 1. / energy[i];
592 for (int l = 0; l < 4; l++) {
593 for (int n = 0; n < 6; n++) {
594 m_D[l][i * KFitConst::kNumber7 + n] = 0;
595 }
596 }
597 m_D[0][i * KFitConst::kNumber7 + 0] = 1;
598 m_D[0][i * KFitConst::kNumber7 + 5] = -a;
599 m_D[1][i * KFitConst::kNumber7 + 1] = 1;
600 m_D[1][i * KFitConst::kNumber7 + 4] = a;
601 m_D[2][i * KFitConst::kNumber7 + 2] = 1;
602 m_D[3][i * KFitConst::kNumber7 + 0] = al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE;
603 m_D[3][i * KFitConst::kNumber7 + 1] = al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE;
604 m_D[3][i * KFitConst::kNumber7 + 2] = al_1_prime[i * KFitConst::kNumber7 + 2][0] * invE;
605 m_D[3][i * KFitConst::kNumber7 + 4] = -al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE * a;
606 m_D[3][i * KFitConst::kNumber7 + 5] = al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE * a;
607 } else {
608 m_D[0][i * KFitConst::kNumber7 + 0] = 1;
609 m_D[1][i * KFitConst::kNumber7 + 1] = 1;
610 m_D[2][i * KFitConst::kNumber7 + 2] = 1;
611 m_D[3][i * KFitConst::kNumber7 + 3] = 1;
612 }
613 }
614
615 } else
616 {
617
618 // m_FlagFitIncludingVertex == true
619 HepMatrix al_1_prime(m_al_1);
620 HepMatrix Sum_al_1(7, 1, 0);
621 std::vector<double> energy(m_TrackCount);
622 double a;
623
624 for (int i = 0; i < m_TrackCount; i++) {
625 a = m_property[i][2];
626 al_1_prime[i * KFitConst::kNumber7 + 0][0] -= a * (al_1_prime[KFitConst::kNumber7 * m_TrackCount + 1][0] - al_1_prime[i *
627 KFitConst::kNumber7 + 5][0]);
628 al_1_prime[i * KFitConst::kNumber7 + 1][0] += a * (al_1_prime[KFitConst::kNumber7 * m_TrackCount + 0][0] - al_1_prime[i *
629 KFitConst::kNumber7 + 4][0]);
630 energy[i] = sqrt(al_1_prime[i * KFitConst::kNumber7 + 0][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] +
631 al_1_prime[i * KFitConst::kNumber7 + 1][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] +
632 al_1_prime[i * KFitConst::kNumber7 + 2][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] +
633 m_property[i][1] * m_property[i][1]);
634 Sum_al_1[6][0] = + a;
635 }
636
637 for (int i = 0; i < m_TrackCount; i++) {
638 if (energy[i] == 0) {
640 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
641 break;
642 }
643
644 if (m_IsFixMass[i]) {
645 double invE = 1. / energy[i];
646 Sum_al_1[3][0] += energy[i];
647 Sum_al_1[4][0] += al_1_prime[i * KFitConst::kNumber7 + 1][0] * m_property[i][2] * invE;
648 Sum_al_1[5][0] += al_1_prime[i * KFitConst::kNumber7 + 0][0] * m_property[i][2] * invE;
649 } else {
650 Sum_al_1[3][0] += al_1_prime[i * KFitConst::kNumber7 + 3][0];
651 }
652
653 for (int j = 0; j < 3; j++) Sum_al_1[j][0] += al_1_prime[i * KFitConst::kNumber7 + j][0];
654 }
655
656 m_d[0][0] = Sum_al_1[0][0] - m_FourMomentum.Px();
657 m_d[1][0] = Sum_al_1[1][0] - m_FourMomentum.Py();
658 m_d[2][0] = Sum_al_1[2][0] - m_FourMomentum.Pz();
659 m_d[3][0] = Sum_al_1[3][0] - m_FourMomentum.E();
660
661 for (int i = 0; i < m_TrackCount; i++) {
662 if (energy[i] == 0) {
664 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
665 break;
666 }
667
668 a = m_property[i][2];
669 if (!m_FlagAtDecayPoint) a = 0.;
670
671 if (m_IsFixMass[i]) {
672 double invE = 1. / energy[i];
673 for (int l = 0; l < 4; l++) {
674 for (int n = 0; n < 6; n++) {
675 m_D[l][i * KFitConst::kNumber7 + n] = 0;
676 }
677 }
678 m_D[0][i * KFitConst::kNumber7 + 0] = 1;
679 m_D[0][i * KFitConst::kNumber7 + 5] = -a;
680 m_D[1][i * KFitConst::kNumber7 + 1] = 1;
681 m_D[1][i * KFitConst::kNumber7 + 4] = a;
682 m_D[2][i * KFitConst::kNumber7 + 2] = 1;
683 m_D[3][i * KFitConst::kNumber7 + 0] = al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE;
684 m_D[3][i * KFitConst::kNumber7 + 1] = al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE;
685 m_D[3][i * KFitConst::kNumber7 + 2] = al_1_prime[i * KFitConst::kNumber7 + 2][0] * invE;
686 m_D[3][i * KFitConst::kNumber7 + 4] = -al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE * a;
687 m_D[3][i * KFitConst::kNumber7 + 5] = al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE * a;
688 } else {
689 m_D[0][i * KFitConst::kNumber7 + 0] = 1;
690 m_D[1][i * KFitConst::kNumber7 + 1] = 1;
691 m_D[2][i * KFitConst::kNumber7 + 2] = 1;
692 m_D[3][i * KFitConst::kNumber7 + 3] = 1;
693 }
694 }
695
696 m_D[0][KFitConst::kNumber7 * m_TrackCount + 0] = 2.*(Sum_al_1[3][0] * Sum_al_1[4][0] - Sum_al_1[1][0] * Sum_al_1[6][0]);
697 m_D[0][KFitConst::kNumber7 * m_TrackCount + 1] = -2.*(Sum_al_1[3][0] * Sum_al_1[5][0] - Sum_al_1[0][0] * Sum_al_1[6][0]);
698 m_D[0][KFitConst::kNumber7 * m_TrackCount + 2] = 0.;
699 }
700
702}
703
704
711
713{
714 MakeMotherKFit kmm;
716 unsigned n = getTrackCount();
717 for (unsigned i = 0; i < n; ++i) {
719 getTrack(i).getCharge());
722 for (unsigned j = i + 1; j < n; ++j) {
724 }
725 }
726 kmm.setVertex(getVertex());
729 m_ErrorCode = kmm.doMake();
731 return m_ErrorCode;
732 double chi2 = getCHIsq();
733 int ndf = getNDF();
734 double prob = TMath::Prob(chi2, ndf);
735 //
736 bool haschi2 = mother->hasExtraInfo("chiSquared");
737 if (haschi2) {
738 mother->setExtraInfo("chiSquared", chi2);
739 mother->setExtraInfo("ndf", ndf);
740 } else {
741 mother->addExtraInfo("chiSquared", chi2);
742 mother->addExtraInfo("ndf", ndf);
743 }
744
745 mother->updateMomentum(
746 CLHEPToROOT::getLorentzVector(kmm.getMotherMomentum()),
747 CLHEPToROOT::getXYZVector(kmm.getMotherPosition()),
748 CLHEPToROOT::getTMatrixFSym(kmm.getMotherError()),
749 prob);
751 return m_ErrorCode;
752}
static const double speedOfLight
[cm/ns]
Definition Const.h:695
Class to store reconstructed particles.
Definition Particle.h:76
void setExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
Definition Particle.cc:1402
bool hasExtraInfo(const std::string &name) const
Return whether the extra info with the given name is set.
Definition Particle.cc:1351
void addExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
Definition Particle.cc:1421
void updateMomentum(const ROOT::Math::PxPyPzEVector &p4, const ROOT::Math::XYZVector &vertex, const TMatrixFSym &errMatrix, double pValue)
Sets Lorentz vector, position, 7x7 error matrix and p-value.
Definition Particle.h:397
enum KFitError::ECode setZeroCorrelation(void) override
Indicate no correlation between tracks.
enum KFitError::ECode setVertex(const HepPoint3D &v)
Set an initial vertex position for the four momentum-constraint fit.
bool m_FlagTrackVertexError
Flag to indicate if the vertex error matrix of the child particle is preset.
enum KFitError::ECode setFourMomentum(const ROOT::Math::PxPyPzEVector &m)
Set an 4 Momentum for the FourC-constraint fit.
bool m_FlagAtDecayPoint
Flag controlled by setFlagAtDecayPoint().
bool getFlagAtDecayPoint(void) const
Get a flag if to constraint at the decay point in the four momentum-constraint fit.
enum KFitError::ECode prepareInputMatrix(void) override
Build grand matrices for minimum search from input-track properties.
bool getFlagFitWithVertex(void) const
Get a flag if the fit is allowed with moving the vertex position.
enum KFitError::ECode calculateNDF(void) override
Calculate an NDF of the fit.
double getCHIsq(void) const override
Get a chi-square of the fit.
enum KFitError::ECode setFlagAtDecayPoint(const bool flag)
Set a flag if to constraint at the decay point in the four momentum-constraint fit.
enum KFitError::ECode setTrackZeroVertexError(void)
Indicate no vertex uncertainty in the child particle in the addTrack'ed order.
std::vector< int > m_IsFixMass
Array of flags whether the track property is fixed at the mass.
enum KFitError::ECode updateMother(Particle *mother)
Update mother particle.
enum KFitError::ECode unfixMass(void)
Tell the object to unfix the last added track property at the invariant mass.
FourCFitKFit(void)
Construct an object with no argument.
enum KFitError::ECode prepareInputSubMatrix(void) override
Build sub-matrices for minimum search from input-track properties.
bool m_FlagFitIncludingVertex
Flag to indicate if the fit is allowed with moving the vertex position.
enum KFitError::ECode doFit(void)
Perform a four momentum-constraint fit.
double m_InvariantMass
Invariant mass.
enum KFitError::ECode setInvariantMass(const double m)
Set an invariant mass for the four momentum-constraint fit.
std::vector< CLHEP::HepMatrix > m_BeforeTrackVertexError
array of vertex error matrices before the fit.
enum KFitError::ECode makeCoreMatrix(void) override
Build matrices using the kinematical constraint.
enum KFitError::ECode setVertexError(const CLHEP::HepSymMatrix &e)
Set an initial vertex error matrix for the four momentum-constraint fit.
enum KFitError::ECode setTrackVertexError(const CLHEP::HepMatrix &e)
Set a vertex error matrix of the child particle in the addTrack'ed order.
~FourCFitKFit(void)
Destruct the object.
enum KFitError::ECode prepareCorrelation(void) override
Build a grand correlation matrix from input-track properties.
HepPoint3D m_AfterVertex
Vertex position after the fit.
enum KFitError::ECode setCorrelation(const CLHEP::HepMatrix &m) override
Set a correlation matrix.
std::vector< CLHEP::HepMatrix > m_AfterTrackVertexError
array of vertex error matrices after the fit.
const CLHEP::HepMatrix getCorrelation(const int id1, const int id2, const int flag=KFitConst::kAfterFit) const override
Get a correlation matrix between two tracks.
double getTrackCHIsq(const int id) const override
Get a chi-square of the track.
const CLHEP::HepSymMatrix getVertexError(const int flag=KFitConst::kAfterFit) const
Get a vertex error matrix.
const HepPoint3D getVertex(const int flag=KFitConst::kAfterFit) const
Get a vertex position.
const CLHEP::HepMatrix getTrackVertexError(const int id, const int flag=KFitConst::kAfterFit) const
Get a vertex error matrix of the track.
CLHEP::HepSymMatrix m_AfterVertexError
Vertex error matrix after the fit.
enum KFitError::ECode prepareOutputMatrix(void) override
Build an output error matrix.
CLHEP::HepSymMatrix m_BeforeVertexError
Vertex error matrix before the fit.
HepPoint3D m_BeforeVertex
Vertex position before the fit.
enum KFitError::ECode fixMass(void)
Tell the object to fix the last added track property at the invariant mass.
ROOT::Math::PxPyPzEVector m_FourMomentum
Four Momentum.
double getInvariantMass(void) const
Get an invariant mass.
int m_NecessaryTrackCount
Number needed tracks to perform fit.
Definition KFitBase.h:303
double m_MagneticField
Magnetic field.
Definition KFitBase.h:311
CLHEP::HepMatrix m_al_1
See J.Tanaka Ph.D (2001) p136 for definition.
Definition KFitBase.h:259
virtual enum KFitError::ECode setCorrelation(const CLHEP::HepMatrix &c)
Set a correlation matrix.
Definition KFitBase.cc:70
const CLHEP::HepSymMatrix makeError3(const CLHEP::HepLorentzVector &p, const CLHEP::HepMatrix &e, const bool is_fix_mass) const
Rebuild an error matrix from a Lorentz vector and an error matrix.
Definition KFitBase.cc:320
const CLHEP::HepSymMatrix getTrackError(const int id) const
Get an error matrix of the track.
Definition KFitBase.cc:168
const CLHEP::HepLorentzVector getTrackMomentum(const int id) const
Get a Lorentz vector of the track.
Definition KFitBase.cc:154
CLHEP::HepMatrix m_lam
See J.Tanaka Ph.D (2001) p137 for definition.
Definition KFitBase.h:276
const HepPoint3D getTrackPosition(const int id) const
Get a position of the track.
Definition KFitBase.cc:161
CLHEP::HepMatrix m_property
Container of charges and masses.
Definition KFitBase.h:263
enum KFitError::ECode m_ErrorCode
Error code.
Definition KFitBase.h:243
virtual enum KFitError::ECode setZeroCorrelation(void)
Indicate no correlation between tracks.
Definition KFitBase.cc:85
CLHEP::HepMatrix m_V_al_1
See J.Tanaka Ph.D (2001) p138 for definition.
Definition KFitBase.h:274
virtual int getNDF(void) const
Get an NDF of the fit.
Definition KFitBase.cc:114
CLHEP::HepMatrix m_d
See J.Tanaka Ph.D (2001) p137 for definition.
Definition KFitBase.h:268
bool isFitted(void) const
Return false if fit is not performed yet or performed fit is failed; otherwise true.
Definition KFitBase.cc:728
CLHEP::HepMatrix m_D
See J.Tanaka Ph.D (2001) p137 for definition.
Definition KFitBase.h:266
CLHEP::HepMatrix m_V_D
See J.Tanaka Ph.D (2001) p138 for definition.
Definition KFitBase.h:271
bool isTrackIDInRange(const int id) const
Check if the id is in the range.
Definition KFitBase.cc:739
virtual const CLHEP::HepMatrix getCorrelation(const int id1, const int id2, const int flag=KFitConst::kAfterFit) const
Get a correlation matrix between two tracks.
Definition KFitBase.cc:183
bool m_FlagCorrelation
Flag whether a correlation among tracks exists.
Definition KFitBase.h:306
CLHEP::HepSymMatrix m_V_al_0
See J.Tanaka Ph.D (2001) p137 for definition.
Definition KFitBase.h:255
const KFitTrack getTrack(const int id) const
Get a specified track object.
Definition KFitBase.cc:175
std::vector< CLHEP::HepMatrix > m_BeforeCorrelation
Container of input correlation matrices.
Definition KFitBase.h:251
bool m_FlagFitted
Flag to indicate if the fit is performed and succeeded.
Definition KFitBase.h:245
double m_CHIsq
chi-square of the fit.
Definition KFitBase.h:297
int getTrackCount(void) const
Get the number of added tracks.
Definition KFitBase.cc:107
int m_NDF
NDF of the fit.
Definition KFitBase.h:295
std::vector< KFitTrack > m_Tracks
Container of input tracks.
Definition KFitBase.h:249
const CLHEP::HepMatrix makeError4(const CLHEP::HepLorentzVector &p, const CLHEP::HepMatrix &e) const
Rebuild an error matrix from a Lorentz vector and an error matrix.
Definition KFitBase.cc:439
int m_TrackCount
Number of tracks.
Definition KFitBase.h:301
CLHEP::HepMatrix m_al_0
See J.Tanaka Ph.D (2001) p136 for definition.
Definition KFitBase.h:257
enum KFitError::ECode doFit1(void)
Perform a fit (used in MassFitKFit::doFit()).
Definition KFitBase.cc:502
static void displayError(const char *file, const int line, const char *func, const enum ECode code)
Display a description of error and its location.
Definition KFitError.h:71
ECode
ECode is a error code enumerate.
Definition KFitError.h:33
@ kCannotGetMatrixInverse
Cannot calculate matrix inverse (bad track property or internal error)
Definition KFitError.h:57
@ kOutOfRange
Specified track-id out of range.
Definition KFitError.h:41
@ kDivisionByZero
Division by zero (bad track property or internal error)
Definition KFitError.h:55
@ kBadTrackSize
Track count too small to perform fit.
Definition KFitError.h:46
@ kBadMatrixSize
Wrong correlation matrix size.
Definition KFitError.h:48
@ kBadCorrelationSize
Wrong correlation matrix size (internal error)
Definition KFitError.h:50
MakeMotherKFit is a class to build mother particle from kinematically fitted daughters.
enum KFitError::ECode setVertex(const HepPoint3D &v)
Set a vertex position of the mother particle.
enum KFitError::ECode addTrack(const KFitTrack &kp)
Add a track to the make-mother object.
enum KFitError::ECode doMake(void)
Perform a reconstruction of mother particle.
const CLHEP::HepSymMatrix getMotherError(void) const
Get an error matrix of the mother particle.
enum KFitError::ECode setCorrelation(const CLHEP::HepMatrix &e)
Set a correlation matrix.
const HepPoint3D getMotherPosition(void) const
Get a position of the mother particle.
enum KFitError::ECode setVertexError(const CLHEP::HepSymMatrix &e)
Set a vertex error matrix of the mother particle.
enum KFitError::ECode setTrackVertexError(const CLHEP::HepMatrix &e)
Set a vertex error matrix of the child particle in the addTrack'ed order.
const CLHEP::HepLorentzVector getMotherMomentum(void) const
Get a Lorentz vector of the mother particle.
enum KFitError::ECode setMagneticField(const double mf)
Change a magnetic field from the default value KFitConst::kDefaultMagneticField.
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
Abstract base class for different kinds of events.
STL namespace.
static const int kMaxTrackCount
Maximum track size.
Definition KFitConst.h:38
static const int kAfterFit
Input parameter to specify after-fit when setting/getting a track attribute.
Definition KFitConst.h:35
static const int kBeforeFit
Input parameter to specify before-fit when setting/getting a track attribute.
Definition KFitConst.h:33
static const int kNumber7
Constant 7 to check matrix size (internal use)
Definition KFitConst.h:30