Belle II Software light-2406-ragdoll
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 <cstdio>
10
11#include <analysis/VertexFitting/KFit/FourCFitKFit.h>
12#include <analysis/VertexFitting/KFit/MakeMotherKFit.h>
13#include <analysis/utility/CLHEPToROOT.h>
14#include <framework/gearbox/Const.h>
15
16#include <TMath.h>
17#include <TMatrixFSym.h>
18
19using namespace std;
20using namespace Belle2;
21using namespace Belle2::analysis;
22using namespace CLHEP;
23using namespace ROOT::Math;
24
26{
27 m_FlagFitted = false;
30 m_FlagAtDecayPoint = true;
32 m_d = HepMatrix(4, 1, 0);
33 m_V_D = HepMatrix(4, 4, 0);
34 m_lam = HepMatrix(4, 1, 0);
35 m_AfterVertexError = HepSymMatrix(3, 0);
36 m_InvariantMass = -1.0;
37 m_FourMomentum = PxPyPzEVector();
38}
39
40
42
43
45FourCFitKFit::setVertex(const HepPoint3D& v) {
47
49}
50
51
53FourCFitKFit::setVertexError(const HepSymMatrix& e) {
54 if (e.num_row() != 3)
55 {
57 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
58 return m_ErrorCode;
59 }
60
63
65}
66
67
71
73}
74
75
77FourCFitKFit::setFourMomentum(const PxPyPzEVector& m) {
79
81}
82
83
86 m_FlagAtDecayPoint = flag;
87
89}
90
91
94 m_IsFixMass.push_back(true);
95
97}
98
99
102 m_IsFixMass.push_back(false);
103
105}
106
107
110 if (e.num_row() != 3 || e.num_col() != KFitConst::kNumber7)
111 {
113 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
114 return m_ErrorCode;
115 }
116
117 m_BeforeTrackVertexError.push_back(e);
120
122}
123
124
127 HepMatrix zero(3, KFitConst::kNumber7, 0);
128
129 return this->setTrackVertexError(zero);
130}
131
132
134FourCFitKFit::setCorrelation(const HepMatrix& m) {
135 return KFitBase::setCorrelation(m);
136}
137
138
142}
143
144
145const HepPoint3D
146FourCFitKFit::getVertex(const int flag) const
147{
148 if (flag == KFitConst::kAfterFit && !isFitted()) return HepPoint3D();
149
150 switch (flag) {
152 return m_BeforeVertex;
153
155 return m_AfterVertex;
156
157 default:
158 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
159 return HepPoint3D();
160 }
161}
162
163
164const HepSymMatrix
165FourCFitKFit::getVertexError(const int flag) const
166{
167 if (flag == KFitConst::kAfterFit && !isFitted()) return HepSymMatrix(3, 0);
168
169 if (flag == KFitConst::kBeforeFit)
170 return m_BeforeVertexError;
172 return m_AfterVertexError;
173 else {
174 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
175 return HepSymMatrix(3, 0);
176 }
177}
178
179
180double
182{
183 return m_InvariantMass;
184}
185
186
187bool
189{
190 return m_FlagAtDecayPoint;
191}
192
193
194bool
196{
198}
199
200
201double
203{
204 return m_CHIsq;
205}
206
207
208const HepMatrix
209FourCFitKFit::getTrackVertexError(const int id, const int flag) const
210{
211 if (flag == KFitConst::kAfterFit && !isFitted()) return HepMatrix(3, KFitConst::kNumber7, 0);
212 if (!isTrackIDInRange(id)) return HepMatrix(3, KFitConst::kNumber7, 0);
213
214 if (flag == KFitConst::kBeforeFit)
215 return m_BeforeTrackVertexError[id];
217 return m_AfterTrackVertexError[id];
218 else {
219 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
220 return HepMatrix(3, KFitConst::kNumber7, 0);
221 }
222}
223
224
225double
227{
228 if (!isFitted()) return -1;
229 if (!isTrackIDInRange(id)) return -1;
230
231 if (m_IsFixMass[id]) {
232
233 HepMatrix da(m_Tracks[id].getFitParameter(KFitConst::kBeforeFit) - m_Tracks[id].getFitParameter(KFitConst::kAfterFit));
234 int err_inverse = 0;
235 const double chisq = (da.T() * (m_Tracks[id].getFitError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
236
237 if (err_inverse) {
238 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
239 return -1;
240 }
241
242 return chisq;
243
244 } else {
245
246 HepMatrix da(m_Tracks[id].getMomPos(KFitConst::kBeforeFit) - m_Tracks[id].getMomPos(KFitConst::kAfterFit));
247 int err_inverse = 0;
248 const double chisq = (da.T() * (m_Tracks[id].getError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
249
250 if (err_inverse) {
251 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
252 return -1;
253 }
254
255 return chisq;
256
257 }
258}
259
260
261const HepMatrix
262FourCFitKFit::getCorrelation(const int id1, const int id2, const int flag) const
263{
264 if (flag == KFitConst::kAfterFit && !isFitted()) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
265 if (!isTrackIDInRange(id1)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
266 if (!isTrackIDInRange(id2)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
267
268 switch (flag) {
270 return KFitBase::getCorrelation(id1, id2, flag);
271
273 return makeError3(
274 this->getTrackMomentum(id1),
275 this->getTrackMomentum(id2),
276 m_V_al_1.sub(KFitConst::kNumber7 * id1 + 1, KFitConst::kNumber7 * (id1 + 1), KFitConst::kNumber7 * id2 + 1,
277 KFitConst::kNumber7 * (id2 + 1)),
278 m_IsFixMass[id1],
279 m_IsFixMass[id2]);
280
281 default:
282 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
283 return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
284 }
285}
286
287
290 return KFitBase::doFit1();
291}
292
293
297 {
299 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
300 return m_ErrorCode;
301 }
302
303
304 if (m_IsFixMass.size() == 0)
305 {
306 // If no fix_mass flag at all,
307 // all tracks are considered to be fixed at mass.
308 for (int i = 0; i < m_TrackCount; i++) this->fixMass();
309 } else if (m_IsFixMass.size() != (unsigned int)m_TrackCount)
310 {
312 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
313 return m_ErrorCode;
314 }
315
316
318 {
319 int index = 0;
320 m_al_0 = HepMatrix(KFitConst::kNumber7 * m_TrackCount, 1, 0);
321 m_property = HepMatrix(m_TrackCount, 3, 0);
322 m_V_al_0 = HepSymMatrix(KFitConst::kNumber7 * m_TrackCount, 0);
323
324 for (auto& track : m_Tracks) {
325 // momentum x,y,z and position x,y,z
326 m_al_0[index * KFitConst::kNumber7 + 0][0] = track.getMomentum(KFitConst::kBeforeFit).x();
327 m_al_0[index * KFitConst::kNumber7 + 1][0] = track.getMomentum(KFitConst::kBeforeFit).y();
328 m_al_0[index * KFitConst::kNumber7 + 2][0] = track.getMomentum(KFitConst::kBeforeFit).z();
329 m_al_0[index * KFitConst::kNumber7 + 3][0] = track.getMomentum(KFitConst::kBeforeFit).t();
330 m_al_0[index * KFitConst::kNumber7 + 4][0] = track.getPosition(KFitConst::kBeforeFit).x();
331 m_al_0[index * KFitConst::kNumber7 + 5][0] = track.getPosition(KFitConst::kBeforeFit).y();
332 m_al_0[index * KFitConst::kNumber7 + 6][0] = track.getPosition(KFitConst::kBeforeFit).z();
333 // these error
334 m_V_al_0.sub(index * KFitConst::kNumber7 + 1, track.getError(KFitConst::kBeforeFit));
335 // charge, mass, a
336 m_property[index][0] = track.getCharge();
337 m_property[index][1] = track.getMass();
338 const double c = Belle2::Const::speedOfLight * 1e-4;
339 m_property[index][2] = -c * m_MagneticField * track.getCharge();
340 index++;
341 }
342
343 // error between track and track
344 if (m_FlagCorrelation) {
345 this->prepareCorrelation();
347 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
348 return m_ErrorCode;
349 }
350 }
351
352 // set member matrix
353 m_al_1 = m_al_0;
354
355 // define size of matrix
358
359 } else
360 {
361 // m_FlagFitIncludingVertex == true
362 int index = 0;
363 m_al_0 = HepMatrix(KFitConst::kNumber7 * m_TrackCount + 3, 1, 0);
364 m_property = HepMatrix(m_TrackCount, 3, 0);
365 m_V_al_0 = HepSymMatrix(KFitConst::kNumber7 * m_TrackCount + 3, 0);
366
367 for (auto& track : m_Tracks) {
368 // momentum x,y,z and position x,y,z
369 m_al_0[index * KFitConst::kNumber7 + 0][0] = track.getMomentum(KFitConst::kBeforeFit).x();
370 m_al_0[index * KFitConst::kNumber7 + 1][0] = track.getMomentum(KFitConst::kBeforeFit).y();
371 m_al_0[index * KFitConst::kNumber7 + 2][0] = track.getMomentum(KFitConst::kBeforeFit).z();
372 m_al_0[index * KFitConst::kNumber7 + 3][0] = track.getMomentum(KFitConst::kBeforeFit).t();
373 m_al_0[index * KFitConst::kNumber7 + 4][0] = track.getPosition(KFitConst::kBeforeFit).x();
374 m_al_0[index * KFitConst::kNumber7 + 5][0] = track.getPosition(KFitConst::kBeforeFit).y();
375 m_al_0[index * KFitConst::kNumber7 + 6][0] = track.getPosition(KFitConst::kBeforeFit).z();
376 // these error
377 m_V_al_0.sub(index * KFitConst::kNumber7 + 1, track.getError(KFitConst::kBeforeFit));
378 // charge, mass, a
379 m_property[index][0] = track.getCharge();
380 m_property[index][1] = track.getMass();
381 const double c = Belle2::Const::speedOfLight * 1e-4;
382 m_property[index][2] = -c * m_MagneticField * track.getCharge();
383 index++;
384 }
385
386 // vertex
391
392 // error between track and track
393 if (m_FlagCorrelation) {
394 this->prepareCorrelation();
396 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
397 return m_ErrorCode;
398 }
399 }
400
401 // set member matrix
402 m_al_1 = m_al_0;
403
404 // define size of matrix
406 m_D = m_V_al_1.sub(1, 4, 1, KFitConst::kNumber7 * m_TrackCount + 3);
407 }
408
410}
411
412
415 char buf[1024];
416 sprintf(buf, "%s:%s(): internal error; this function should never be called", __FILE__, __func__);
417 B2FATAL(buf);
418
419 /* NEVER REACHEd HERE */
421}
422
423
426 if (m_BeforeCorrelation.size() != static_cast<unsigned int>(m_TrackCount * (m_TrackCount - 1) / 2))
427 {
429 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
430 return m_ErrorCode;
431 }
432
433 int row = 0, col = 0;
434
435 for (auto& hm : m_BeforeCorrelation)
436 {
437 // counter
438 row++;
439 if (row == m_TrackCount) {
440 col++;
441 row = col + 1;
442 }
443
444 int ii = 0, jj = 0;
445 for (int i = KFitConst::kNumber7 * row; i < KFitConst::kNumber7 * (row + 1); i++) {
446 for (int j = KFitConst::kNumber7 * col; j < KFitConst::kNumber7 * (col + 1); j++) {
447 m_V_al_0[i][j] = hm[ii][jj];
448 jj++;
449 }
450 jj = 0;
451 ii++;
452 }
453 }
454
456 {
457 // ...error of vertex
459
460 // ...error matrix between vertex and tracks
462 if (m_BeforeTrackVertexError.size() != (unsigned int)m_TrackCount) {
464 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
465 return m_ErrorCode;
466 }
467
468 int i = 0;
469 for (auto& hm : m_BeforeTrackVertexError) {
470 for (int j = 0; j < 3; j++) for (int k = 0; k < KFitConst::kNumber7; k++) {
472 }
473 i++;
474 }
475 }
476 }
477
479}
480
481
484 Hep3Vector h3v;
485 int index = 0;
486 for (auto& pdata : m_Tracks)
487 {
488 // tracks
489 // momentum
490 h3v.setX(m_al_1[index * KFitConst::kNumber7 + 0][0]);
491 h3v.setY(m_al_1[index * KFitConst::kNumber7 + 1][0]);
492 h3v.setZ(m_al_1[index * KFitConst::kNumber7 + 2][0]);
493 if (m_IsFixMass[index])
494 pdata.setMomentum(HepLorentzVector(h3v, sqrt(h3v.mag2() + pdata.getMass()*pdata.getMass())), KFitConst::kAfterFit);
495 else
496 pdata.setMomentum(HepLorentzVector(h3v, m_al_1[index * KFitConst::kNumber7 + 3][0]), KFitConst::kAfterFit);
497 // position
498 pdata.setPosition(HepPoint3D(
499 m_al_1[index * KFitConst::kNumber7 + 4][0],
500 m_al_1[index * KFitConst::kNumber7 + 5][0],
502 // error of the tracks
503 pdata.setError(this->makeError3(pdata.getMomentum(),
504 m_V_al_1.sub(
505 index * KFitConst::kNumber7 + 1,
506 (index + 1)*KFitConst::kNumber7,
507 index * KFitConst::kNumber7 + 1,
508 (index + 1)*KFitConst::kNumber7), m_IsFixMass[index]),
510 if (m_ErrorCode != KFitError::kNoError) break;
511 index++;
512 }
513
515 {
516 // vertex
520 // error of the vertex
521 for (int i = 0; i < 3; i++) for (int j = i; j < 3; j++) {
523 }
524 // error between vertex and tracks
525 for (int i = 0; i < m_TrackCount; i++) {
526 HepMatrix hm(3, KFitConst::kNumber7, 0);
527 for (int j = 0; j < 3; j++) for (int k = 0; k < KFitConst::kNumber7; k++) {
529 }
530 if (m_IsFixMass[i])
531 m_AfterTrackVertexError.push_back(this->makeError4(m_Tracks[i].getMomentum(), hm));
532 else
533 m_AfterTrackVertexError.push_back(hm);
534 }
535 } else
536 {
537 // not fit
539 }
540
542}
543
544
548 {
549
550 HepMatrix al_1_prime(m_al_1);
551 HepMatrix Sum_al_1(4, 1, 0);
552 std::vector<double> energy(m_TrackCount);
553 double a;
554
555 for (int i = 0; i < m_TrackCount; i++) {
556 a = m_property[i][2];
557 if (!m_FlagAtDecayPoint) a = 0.;
558 al_1_prime[i * KFitConst::kNumber7 + 0][0] -= a * (m_BeforeVertex.y() - al_1_prime[i * KFitConst::kNumber7 + 5][0]);
559 al_1_prime[i * KFitConst::kNumber7 + 1][0] += a * (m_BeforeVertex.x() - al_1_prime[i * KFitConst::kNumber7 + 4][0]);
560 energy[i] = sqrt(al_1_prime[i * KFitConst::kNumber7 + 0][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] +
561 al_1_prime[i * KFitConst::kNumber7 + 1][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] +
562 al_1_prime[i * KFitConst::kNumber7 + 2][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] +
563 m_property[i][1] * m_property[i][1]);
564 if (m_IsFixMass[i])
565 Sum_al_1[3][0] += energy[i];
566 else
567 Sum_al_1[3][0] += al_1_prime[i * KFitConst::kNumber7 + 3][0];
568 }
569
570 for (int i = 0; i < m_TrackCount; i++) {
571 for (int j = 0; j < 3; j++) Sum_al_1[j][0] += al_1_prime[i * KFitConst::kNumber7 + j][0];
572 }
573
574 m_d[0][0] = Sum_al_1[0][0] - m_FourMomentum.Px();
575 m_d[1][0] = Sum_al_1[1][0] - m_FourMomentum.Py();
576 m_d[2][0] = Sum_al_1[2][0] - m_FourMomentum.Pz();
577 m_d[3][0] = Sum_al_1[3][0] - m_FourMomentum.E();
578
579 for (int i = 0; i < m_TrackCount; i++) {
580 if (energy[i] == 0) {
582 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
583 break;
584 }
585
586 a = m_property[i][2];
587 if (!m_FlagAtDecayPoint) a = 0.;
588
589 if (m_IsFixMass[i]) {
590 double invE = 1. / energy[i];
591 for (int l = 0; l < 4; l++) {
592 for (int n = 0; n < 6; n++) {
593 m_D[l][i * KFitConst::kNumber7 + n] = 0;
594 }
595 }
596 m_D[0][i * KFitConst::kNumber7 + 0] = 1;
597 m_D[0][i * KFitConst::kNumber7 + 5] = -a;
598 m_D[1][i * KFitConst::kNumber7 + 1] = 1;
599 m_D[1][i * KFitConst::kNumber7 + 4] = a;
600 m_D[2][i * KFitConst::kNumber7 + 2] = 1;
601 m_D[3][i * KFitConst::kNumber7 + 0] = al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE;
602 m_D[3][i * KFitConst::kNumber7 + 1] = al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE;
603 m_D[3][i * KFitConst::kNumber7 + 2] = al_1_prime[i * KFitConst::kNumber7 + 2][0] * invE;
604 m_D[3][i * KFitConst::kNumber7 + 4] = -al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE * a;
605 m_D[3][i * KFitConst::kNumber7 + 5] = al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE * a;
606 } else {
607 m_D[0][i * KFitConst::kNumber7 + 0] = 1;
608 m_D[1][i * KFitConst::kNumber7 + 1] = 1;
609 m_D[2][i * KFitConst::kNumber7 + 2] = 1;
610 m_D[3][i * KFitConst::kNumber7 + 3] = 1;
611 }
612 }
613
614 } else
615 {
616
617 // m_FlagFitIncludingVertex == true
618 HepMatrix al_1_prime(m_al_1);
619 HepMatrix Sum_al_1(7, 1, 0);
620 std::vector<double> energy(m_TrackCount);
621 double a;
622
623 for (int i = 0; i < m_TrackCount; i++) {
624 a = m_property[i][2];
625 al_1_prime[i * KFitConst::kNumber7 + 0][0] -= a * (al_1_prime[KFitConst::kNumber7 * m_TrackCount + 1][0] - al_1_prime[i *
626 KFitConst::kNumber7 + 5][0]);
627 al_1_prime[i * KFitConst::kNumber7 + 1][0] += a * (al_1_prime[KFitConst::kNumber7 * m_TrackCount + 0][0] - al_1_prime[i *
628 KFitConst::kNumber7 + 4][0]);
629 energy[i] = sqrt(al_1_prime[i * KFitConst::kNumber7 + 0][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] +
630 al_1_prime[i * KFitConst::kNumber7 + 1][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] +
631 al_1_prime[i * KFitConst::kNumber7 + 2][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] +
632 m_property[i][1] * m_property[i][1]);
633 Sum_al_1[6][0] = + a;
634 }
635
636 for (int i = 0; i < m_TrackCount; i++) {
637 if (energy[i] == 0) {
639 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
640 break;
641 }
642
643 if (m_IsFixMass[i]) {
644 double invE = 1. / energy[i];
645 Sum_al_1[3][0] += energy[i];
646 Sum_al_1[4][0] += al_1_prime[i * KFitConst::kNumber7 + 1][0] * m_property[i][2] * invE;
647 Sum_al_1[5][0] += al_1_prime[i * KFitConst::kNumber7 + 0][0] * m_property[i][2] * invE;
648 } else {
649 Sum_al_1[3][0] += al_1_prime[i * KFitConst::kNumber7 + 3][0];
650 }
651
652 for (int j = 0; j < 3; j++) Sum_al_1[j][0] += al_1_prime[i * KFitConst::kNumber7 + j][0];
653 }
654
655 m_d[0][0] = Sum_al_1[0][0] - m_FourMomentum.Px();
656 m_d[1][0] = Sum_al_1[1][0] - m_FourMomentum.Py();
657 m_d[2][0] = Sum_al_1[2][0] - m_FourMomentum.Pz();
658 m_d[3][0] = Sum_al_1[3][0] - m_FourMomentum.E();
659
660 for (int i = 0; i < m_TrackCount; i++) {
661 if (energy[i] == 0) {
663 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
664 break;
665 }
666
667 a = m_property[i][2];
668 if (!m_FlagAtDecayPoint) a = 0.;
669
670 if (m_IsFixMass[i]) {
671 double invE = 1. / energy[i];
672 for (int l = 0; l < 4; l++) {
673 for (int n = 0; n < 6; n++) {
674 m_D[l][i * KFitConst::kNumber7 + n] = 0;
675 }
676 }
677 m_D[0][i * KFitConst::kNumber7 + 0] = 1;
678 m_D[0][i * KFitConst::kNumber7 + 5] = -a;
679 m_D[1][i * KFitConst::kNumber7 + 1] = 1;
680 m_D[1][i * KFitConst::kNumber7 + 4] = a;
681 m_D[2][i * KFitConst::kNumber7 + 2] = 1;
682 m_D[3][i * KFitConst::kNumber7 + 0] = al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE;
683 m_D[3][i * KFitConst::kNumber7 + 1] = al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE;
684 m_D[3][i * KFitConst::kNumber7 + 2] = al_1_prime[i * KFitConst::kNumber7 + 2][0] * invE;
685 m_D[3][i * KFitConst::kNumber7 + 4] = -al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE * a;
686 m_D[3][i * KFitConst::kNumber7 + 5] = al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE * a;
687 } else {
688 m_D[0][i * KFitConst::kNumber7 + 0] = 1;
689 m_D[1][i * KFitConst::kNumber7 + 1] = 1;
690 m_D[2][i * KFitConst::kNumber7 + 2] = 1;
691 m_D[3][i * KFitConst::kNumber7 + 3] = 1;
692 }
693 }
694
695 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]);
696 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]);
697 m_D[0][KFitConst::kNumber7 * m_TrackCount + 2] = 0.;
698 }
699
701}
702
703
706 m_NDF = 4;
707
709}
710
712{
713 MakeMotherKFit kmm;
715 unsigned n = getTrackCount();
716 for (unsigned i = 0; i < n; ++i) {
718 getTrack(i).getCharge());
721 for (unsigned j = i + 1; j < n; ++j) {
723 }
724 }
725 kmm.setVertex(getVertex());
728 m_ErrorCode = kmm.doMake();
730 return m_ErrorCode;
731 double chi2 = getCHIsq();
732 int ndf = getNDF();
733 double prob = TMath::Prob(chi2, ndf);
734 //
735 bool haschi2 = mother->hasExtraInfo("chiSquared");
736 if (haschi2) {
737 mother->setExtraInfo("chiSquared", chi2);
738 mother->setExtraInfo("ndf", ndf);
739 } else {
740 mother->addExtraInfo("chiSquared", chi2);
741 mother->addExtraInfo("ndf", ndf);
742 }
743
744 mother->updateMomentum(
745 CLHEPToROOT::getLorentzVector(kmm.getMotherMomentum()),
746 CLHEPToROOT::getXYZVector(kmm.getMotherPosition()),
747 CLHEPToROOT::getTMatrixFSym(kmm.getMotherError()),
748 prob);
750 return m_ErrorCode;
751}
static const double speedOfLight
[cm/ns]
Definition: Const.h:695
Class to store reconstructed particles.
Definition: Particle.h:75
void setExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1317
bool hasExtraInfo(const std::string &name) const
Return whether the extra info with the given name is set.
Definition: Particle.cc:1266
void addExtraInfo(const std::string &name, double value)
Sets the user-defined data of given name to the given value.
Definition: Particle.cc:1336
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:386
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.
Definition: FourCFitKFit.cc:45
bool m_FlagTrackVertexError
Flag to indicate if the vertex error matrix of the child particle is preset.
Definition: FourCFitKFit.h:158
enum KFitError::ECode setFourMomentum(const ROOT::Math::PxPyPzEVector &m)
Set an 4 Momentum for the FourC-constraint fit.
Definition: FourCFitKFit.cc:77
bool m_FlagAtDecayPoint
Flag controlled by setFlagAtDecayPoint().
Definition: FourCFitKFit.h:162
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.
Definition: FourCFitKFit.cc:85
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.
Definition: FourCFitKFit.h:171
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.
Definition: FourCFitKFit.cc:25
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.
Definition: FourCFitKFit.h:160
enum KFitError::ECode doFit(void)
Perform a four momentum-constraint fit.
double m_InvariantMass
Invariant mass.
Definition: FourCFitKFit.h:165
enum KFitError::ECode setInvariantMass(const double m)
Set an invariant mass for the four momentum-constraint fit.
Definition: FourCFitKFit.cc:69
std::vector< CLHEP::HepMatrix > m_BeforeTrackVertexError
array of vertex error matrices before the fit.
Definition: FourCFitKFit.h:148
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.
Definition: FourCFitKFit.cc:53
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.
Definition: FourCFitKFit.h:151
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.
Definition: FourCFitKFit.h:155
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.
Definition: FourCFitKFit.h:153
enum KFitError::ECode prepareOutputMatrix(void) override
Build an output error matrix.
CLHEP::HepSymMatrix m_BeforeVertexError
Vertex error matrix before the fit.
Definition: FourCFitKFit.h:146
HepPoint3D m_BeforeVertex
Vertex position before the fit.
Definition: FourCFitKFit.h:144
enum KFitError::ECode fixMass(void)
Tell the object to fix the last added track property at the invariant mass.
Definition: FourCFitKFit.cc:93
ROOT::Math::PxPyPzEVector m_FourMomentum
Four Momentum.
Definition: FourCFitKFit.h:168
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:72
ECode
ECode is a error code enumerate.
Definition: KFitError.h:34
@ kCannotGetMatrixInverse
Cannot calculate matrix inverse (bad track property or internal error)
Definition: KFitError.h:58
@ kOutOfRange
Specified track-id out of range.
Definition: KFitError.h:42
@ kDivisionByZero
Division by zero (bad track property or internal error)
Definition: KFitError.h:56
@ kBadTrackSize
Track count too small to perform fit.
Definition: KFitError.h:47
@ kBadMatrixSize
Wrong correlation matrix size.
Definition: KFitError.h:49
@ kBadCorrelationSize
Wrong correlation matrix size (internal error)
Definition: KFitError.h:51
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.
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24
STL namespace.
static const int kMaxTrackCount
Maximum track size.
Definition: KFitConst.h:40
static const int kAfterFit
Input parameter to specify after-fit when setting/getting a track attribute.
Definition: KFitConst.h:37
static const int kBeforeFit
Input parameter to specify before-fit when setting/getting a track attribute.
Definition: KFitConst.h:35
static const int kNumber7
Constant 7 to check matrix size (internal use)
Definition: KFitConst.h:32