Belle II Software light-2406-ragdoll
MassFitKFit.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * External Contributor: J. Tanaka *
5 * *
6 * See git log for contributors and copyright holders. *
7 * This file is licensed under LGPL-3.0, see LICENSE.md. *
8 **************************************************************************/
9
10#include <cstdio>
11
12#include <analysis/VertexFitting/KFit/MakeMotherKFit.h>
13#include <analysis/VertexFitting/KFit/MassFitKFit.h>
14#include <analysis/utility/CLHEPToROOT.h>
15#include <framework/gearbox/Const.h>
16
17#include <TMath.h>
18#include <TMatrixFSym.h>
19
20using namespace std;
21using namespace Belle2;
22using namespace Belle2::analysis;
23using namespace CLHEP;
24
26{
27 m_FlagFitted = false;
30 m_FlagAtDecayPoint = true;
32 m_d = HepMatrix(1, 1, 0);
33 m_V_D = HepMatrix(1, 1, 0);
34 m_lam = HepMatrix(1, 1, 0);
35 m_AfterVertexError = HepSymMatrix(3, 0);
36 m_InvariantMass = -1.0;
37}
38
39
41
42
44MassFitKFit::setVertex(const HepPoint3D& v) {
46
48}
49
50
52MassFitKFit::setVertexError(const HepSymMatrix& e) {
53 if (e.num_row() != 3)
54 {
56 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
57 return m_ErrorCode;
58 }
59
62
64}
65
66
70
72}
73
74
77 m_FlagAtDecayPoint = flag;
78
80}
81
82
85 m_IsFixMass.push_back(true);
86
88}
89
90
93 m_IsFixMass.push_back(false);
94
96}
97
98
101 if (e.num_row() != 3 || e.num_col() != KFitConst::kNumber7)
102 {
104 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
105 return m_ErrorCode;
106 }
107
108 m_BeforeTrackVertexError.push_back(e);
111
113}
114
115
118 HepMatrix zero(3, KFitConst::kNumber7, 0);
119
120 return this->setTrackVertexError(zero);
121}
122
123
125MassFitKFit::setCorrelation(const HepMatrix& m) {
126 return KFitBase::setCorrelation(m);
127}
128
129
133}
134
135
136const HepPoint3D
137MassFitKFit::getVertex(const int flag) const
138{
139 if (flag == KFitConst::kAfterFit && !isFitted()) return HepPoint3D();
140
141 switch (flag) {
143 return m_BeforeVertex;
144
146 return m_AfterVertex;
147
148 default:
149 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
150 return HepPoint3D();
151 }
152}
153
154
155const HepSymMatrix
156MassFitKFit::getVertexError(const int flag) const
157{
158 if (flag == KFitConst::kAfterFit && !isFitted()) return HepSymMatrix(3, 0);
159
160 if (flag == KFitConst::kBeforeFit)
161 return m_BeforeVertexError;
163 return m_AfterVertexError;
164 else {
165 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
166 return HepSymMatrix(3, 0);
167 }
168}
169
170
171double
173{
174 return m_InvariantMass;
175}
176
177
178bool
180{
181 return m_FlagAtDecayPoint;
182}
183
184
185bool
187{
189}
190
191
192double
194{
195 return m_CHIsq;
196}
197
198
199const HepMatrix
200MassFitKFit::getTrackVertexError(const int id, const int flag) const
201{
202 if (flag == KFitConst::kAfterFit && !isFitted()) return HepMatrix(3, KFitConst::kNumber7, 0);
203 if (!isTrackIDInRange(id)) return HepMatrix(3, KFitConst::kNumber7, 0);
204
205 if (flag == KFitConst::kBeforeFit)
206 return m_BeforeTrackVertexError[id];
208 return m_AfterTrackVertexError[id];
209 else {
210 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
211 return HepMatrix(3, KFitConst::kNumber7, 0);
212 }
213}
214
215
216double
217MassFitKFit::getTrackCHIsq(const int id) const
218{
219 if (!isFitted()) return -1;
220 if (!isTrackIDInRange(id)) return -1;
221
222 if (m_IsFixMass[id]) {
223
224 HepMatrix da(m_Tracks[id].getFitParameter(KFitConst::kBeforeFit) - m_Tracks[id].getFitParameter(KFitConst::kAfterFit));
225 int err_inverse = 0;
226 const double chisq = (da.T() * (m_Tracks[id].getFitError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
227
228 if (err_inverse) {
229 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
230 return -1;
231 }
232
233 return chisq;
234
235 } else {
236
237 HepMatrix da(m_Tracks[id].getMomPos(KFitConst::kBeforeFit) - m_Tracks[id].getMomPos(KFitConst::kAfterFit));
238 int err_inverse = 0;
239 const double chisq = (da.T() * (m_Tracks[id].getError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
240
241 if (err_inverse) {
242 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
243 return -1;
244 }
245
246 return chisq;
247
248 }
249}
250
251
252const HepMatrix
253MassFitKFit::getCorrelation(const int id1, const int id2, const int flag) const
254{
255 if (flag == KFitConst::kAfterFit && !isFitted()) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
256 if (!isTrackIDInRange(id1)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
257 if (!isTrackIDInRange(id2)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
258
259 switch (flag) {
261 return KFitBase::getCorrelation(id1, id2, flag);
262
264 return makeError3(
265 this->getTrackMomentum(id1),
266 this->getTrackMomentum(id2),
267 m_V_al_1.sub(KFitConst::kNumber7 * id1 + 1, KFitConst::kNumber7 * (id1 + 1), KFitConst::kNumber7 * id2 + 1,
268 KFitConst::kNumber7 * (id2 + 1)),
269 m_IsFixMass[id1],
270 m_IsFixMass[id2]);
271
272 default:
273 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
274 return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
275 }
276}
277
278
281 return KFitBase::doFit1();
282}
283
284
288 {
290 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
291 return m_ErrorCode;
292 }
293
294
295 if (m_IsFixMass.size() == 0)
296 {
297 // If no fix_mass flag at all,
298 // all tracks are considered to be fixed at mass.
299 for (int i = 0; i < m_TrackCount; i++) this->fixMass();
300 } else if (m_IsFixMass.size() != (unsigned int)m_TrackCount)
301 {
303 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
304 return m_ErrorCode;
305 }
306
307
309 {
310 int index = 0;
311 m_al_0 = HepMatrix(KFitConst::kNumber7 * m_TrackCount, 1, 0);
312 m_property = HepMatrix(m_TrackCount, 3, 0);
313 m_V_al_0 = HepSymMatrix(KFitConst::kNumber7 * m_TrackCount, 0);
314
315 for (auto& track : m_Tracks) {
316 // momentum x,y,z and position x,y,z
317 m_al_0[index * KFitConst::kNumber7 + 0][0] = track.getMomentum(KFitConst::kBeforeFit).x();
318 m_al_0[index * KFitConst::kNumber7 + 1][0] = track.getMomentum(KFitConst::kBeforeFit).y();
319 m_al_0[index * KFitConst::kNumber7 + 2][0] = track.getMomentum(KFitConst::kBeforeFit).z();
320 m_al_0[index * KFitConst::kNumber7 + 3][0] = track.getMomentum(KFitConst::kBeforeFit).t();
321 m_al_0[index * KFitConst::kNumber7 + 4][0] = track.getPosition(KFitConst::kBeforeFit).x();
322 m_al_0[index * KFitConst::kNumber7 + 5][0] = track.getPosition(KFitConst::kBeforeFit).y();
323 m_al_0[index * KFitConst::kNumber7 + 6][0] = track.getPosition(KFitConst::kBeforeFit).z();
324 // these error
325 m_V_al_0.sub(index * KFitConst::kNumber7 + 1, track.getError(KFitConst::kBeforeFit));
326 // charge, mass, a
327 m_property[index][0] = track.getCharge();
328 m_property[index][1] = track.getMass();
329 const double c = Belle2::Const::speedOfLight * 1e-4;
330 m_property[index][2] = -c * m_MagneticField * track.getCharge();
331 index++;
332 }
333
334 // error between track and track
335 if (m_FlagCorrelation) {
336 this->prepareCorrelation();
338 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
339 return m_ErrorCode;
340 }
341 }
342
343 // set member matrix
344 m_al_1 = m_al_0;
345
346 // define size of matrix
349
350 } else
351 {
352 // m_FlagFitIncludingVertex == true
353 int index = 0;
354 m_al_0 = HepMatrix(KFitConst::kNumber7 * m_TrackCount + 3, 1, 0);
355 m_property = HepMatrix(m_TrackCount, 3, 0);
356 m_V_al_0 = HepSymMatrix(KFitConst::kNumber7 * m_TrackCount + 3, 0);
357
358 for (auto& track : m_Tracks) {
359 // momentum x,y,z and position x,y,z
360 m_al_0[index * KFitConst::kNumber7 + 0][0] = track.getMomentum(KFitConst::kBeforeFit).x();
361 m_al_0[index * KFitConst::kNumber7 + 1][0] = track.getMomentum(KFitConst::kBeforeFit).y();
362 m_al_0[index * KFitConst::kNumber7 + 2][0] = track.getMomentum(KFitConst::kBeforeFit).z();
363 m_al_0[index * KFitConst::kNumber7 + 3][0] = track.getMomentum(KFitConst::kBeforeFit).t();
364 m_al_0[index * KFitConst::kNumber7 + 4][0] = track.getPosition(KFitConst::kBeforeFit).x();
365 m_al_0[index * KFitConst::kNumber7 + 5][0] = track.getPosition(KFitConst::kBeforeFit).y();
366 m_al_0[index * KFitConst::kNumber7 + 6][0] = track.getPosition(KFitConst::kBeforeFit).z();
367 // these error
368 m_V_al_0.sub(index * KFitConst::kNumber7 + 1, track.getError(KFitConst::kBeforeFit));
369 // charge, mass, a
370 m_property[index][0] = track.getCharge();
371 m_property[index][1] = track.getMass();
372 const double c = Belle2::Const::speedOfLight * 1e-4;
373 m_property[index][2] = -c * m_MagneticField * track.getCharge();
374 index++;
375 }
376
377 // vertex
382
383 // error between track and track
384 if (m_FlagCorrelation) {
385 this->prepareCorrelation();
387 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
388 return m_ErrorCode;
389 }
390 }
391
392 // set member matrix
393 m_al_1 = m_al_0;
394
395 // define size of matrix
397 m_D = m_V_al_1.sub(1, 1, 1, KFitConst::kNumber7 * m_TrackCount + 3);
398 }
399
401}
402
403
406 char buf[1024];
407 sprintf(buf, "%s:%s(): internal error; this function should never be called", __FILE__, __func__);
408 B2FATAL(buf);
409
410 /* NEVER REACHEd HERE */
412}
413
414
417 if (m_BeforeCorrelation.size() != static_cast<unsigned int>(m_TrackCount * (m_TrackCount - 1) / 2))
418 {
420 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
421 return m_ErrorCode;
422 }
423
424 int row = 0, col = 0;
425
426 for (auto& hm : m_BeforeCorrelation)
427 {
428 // counter
429 row++;
430 if (row == m_TrackCount) {
431 col++;
432 row = col + 1;
433 }
434
435 int ii = 0, jj = 0;
436 for (int i = KFitConst::kNumber7 * row; i < KFitConst::kNumber7 * (row + 1); i++) {
437 for (int j = KFitConst::kNumber7 * col; j < KFitConst::kNumber7 * (col + 1); j++) {
438 m_V_al_0[i][j] = hm[ii][jj];
439 jj++;
440 }
441 jj = 0;
442 ii++;
443 }
444 }
445
447 {
448 // ...error of vertex
450
451 // ...error matrix between vertex and tracks
453 if (m_BeforeTrackVertexError.size() != (unsigned int)m_TrackCount) {
455 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
456 return m_ErrorCode;
457 }
458
459 int i = 0;
460 for (auto& hm : m_BeforeTrackVertexError) {
461 for (int j = 0; j < 3; j++) for (int k = 0; k < KFitConst::kNumber7; k++) {
463 }
464 i++;
465 }
466 }
467 }
468
470}
471
472
475 Hep3Vector h3v;
476 int index = 0;
477 for (auto& pdata : m_Tracks)
478 {
479 // tracks
480 // momentum
481 h3v.setX(m_al_1[index * KFitConst::kNumber7 + 0][0]);
482 h3v.setY(m_al_1[index * KFitConst::kNumber7 + 1][0]);
483 h3v.setZ(m_al_1[index * KFitConst::kNumber7 + 2][0]);
484 if (m_IsFixMass[index])
485 pdata.setMomentum(HepLorentzVector(h3v, sqrt(h3v.mag2() + pdata.getMass()*pdata.getMass())), KFitConst::kAfterFit);
486 else
487 pdata.setMomentum(HepLorentzVector(h3v, m_al_1[index * KFitConst::kNumber7 + 3][0]), KFitConst::kAfterFit);
488 // position
489 pdata.setPosition(HepPoint3D(
490 m_al_1[index * KFitConst::kNumber7 + 4][0],
491 m_al_1[index * KFitConst::kNumber7 + 5][0],
493 // error of the tracks
494 pdata.setError(this->makeError3(pdata.getMomentum(),
495 m_V_al_1.sub(
496 index * KFitConst::kNumber7 + 1,
497 (index + 1)*KFitConst::kNumber7,
498 index * KFitConst::kNumber7 + 1,
499 (index + 1)*KFitConst::kNumber7), m_IsFixMass[index]),
501 if (m_ErrorCode != KFitError::kNoError) break;
502 index++;
503 }
504
506 {
507 // vertex
511 // error of the vertex
512 for (int i = 0; i < 3; i++) for (int j = i; j < 3; j++) {
514 }
515 // error between vertex and tracks
516 for (int i = 0; i < m_TrackCount; i++) {
517 HepMatrix hm(3, KFitConst::kNumber7, 0);
518 for (int j = 0; j < 3; j++) for (int k = 0; k < KFitConst::kNumber7; k++) {
520 }
521 if (m_IsFixMass[i])
522 m_AfterTrackVertexError.push_back(this->makeError4(m_Tracks[i].getMomentum(), hm));
523 else
524 m_AfterTrackVertexError.push_back(hm);
525 }
526 } else
527 {
528 // not fit
530 }
531
533}
534
535
539 {
540
541 HepMatrix al_1_prime(m_al_1);
542 HepMatrix Sum_al_1(4, 1, 0);
543 std::vector<double> energy(m_TrackCount);
544 double a;
545
546 for (int i = 0; i < m_TrackCount; i++) {
547 a = m_property[i][2];
548 if (!m_FlagAtDecayPoint) a = 0.;
549 al_1_prime[i * KFitConst::kNumber7 + 0][0] -= a * (m_BeforeVertex.y() - al_1_prime[i * KFitConst::kNumber7 + 5][0]);
550 al_1_prime[i * KFitConst::kNumber7 + 1][0] += a * (m_BeforeVertex.x() - al_1_prime[i * KFitConst::kNumber7 + 4][0]);
551 energy[i] = sqrt(al_1_prime[i * KFitConst::kNumber7 + 0][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] +
552 al_1_prime[i * KFitConst::kNumber7 + 1][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] +
553 al_1_prime[i * KFitConst::kNumber7 + 2][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] +
554 m_property[i][1] * m_property[i][1]);
555 if (m_IsFixMass[i])
556 Sum_al_1[3][0] += energy[i];
557 else
558 Sum_al_1[3][0] += al_1_prime[i * KFitConst::kNumber7 + 3][0];
559 }
560
561 for (int i = 0; i < m_TrackCount; i++) {
562 for (int j = 0; j < 3; j++) Sum_al_1[j][0] += al_1_prime[i * KFitConst::kNumber7 + j][0];
563 }
564
565 m_d[0][0] =
566 + Sum_al_1[3][0] * Sum_al_1[3][0] - Sum_al_1[0][0] * Sum_al_1[0][0]
567 - Sum_al_1[1][0] * Sum_al_1[1][0] - Sum_al_1[2][0] * Sum_al_1[2][0]
569
570 for (int i = 0; i < m_TrackCount; i++) {
571 if (energy[i] == 0) {
573 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
574 break;
575 }
576
577 a = m_property[i][2];
578 if (!m_FlagAtDecayPoint) a = 0.;
579
580 if (m_IsFixMass[i]) {
581 double invE = 1. / energy[i];
582 m_D[0][i * KFitConst::kNumber7 + 0] = 2.*(Sum_al_1[3][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE - Sum_al_1[0][0]);
583 m_D[0][i * KFitConst::kNumber7 + 1] = 2.*(Sum_al_1[3][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE - Sum_al_1[1][0]);
584 m_D[0][i * KFitConst::kNumber7 + 2] = 2.*(Sum_al_1[3][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] * invE - Sum_al_1[2][0]);
585 m_D[0][i * KFitConst::kNumber7 + 3] = 0.;
586 m_D[0][i * KFitConst::kNumber7 + 4] = -2.*(Sum_al_1[3][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE - Sum_al_1[1][0]) * a;
587 m_D[0][i * KFitConst::kNumber7 + 5] = 2.*(Sum_al_1[3][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE - Sum_al_1[0][0]) * a;
588 m_D[0][i * KFitConst::kNumber7 + 6] = 0.;
589 } else {
590 m_D[0][i * KFitConst::kNumber7 + 0] = -2.*Sum_al_1[0][0];
591 m_D[0][i * KFitConst::kNumber7 + 1] = -2.*Sum_al_1[1][0];
592 m_D[0][i * KFitConst::kNumber7 + 2] = -2.*Sum_al_1[2][0];
593 m_D[0][i * KFitConst::kNumber7 + 3] = 2.*Sum_al_1[3][0];
594 m_D[0][i * KFitConst::kNumber7 + 4] = 2.*Sum_al_1[1][0] * a;
595 m_D[0][i * KFitConst::kNumber7 + 5] = -2.*Sum_al_1[0][0] * a;
596 m_D[0][i * KFitConst::kNumber7 + 6] = 0.;
597 }
598 }
599
600 } else
601 {
602
603 // m_FlagFitIncludingVertex == true
604 HepMatrix al_1_prime(m_al_1);
605 HepMatrix Sum_al_1(7, 1, 0);
606 std::vector<double> energy(m_TrackCount);
607 double a;
608
609 for (int i = 0; i < m_TrackCount; i++) {
610 a = m_property[i][2];
611 al_1_prime[i * KFitConst::kNumber7 + 0][0] -= a * (al_1_prime[KFitConst::kNumber7 * m_TrackCount + 1][0] - al_1_prime[i *
612 KFitConst::kNumber7 + 5][0]);
613 al_1_prime[i * KFitConst::kNumber7 + 1][0] += a * (al_1_prime[KFitConst::kNumber7 * m_TrackCount + 0][0] - al_1_prime[i *
614 KFitConst::kNumber7 + 4][0]);
615 energy[i] = sqrt(al_1_prime[i * KFitConst::kNumber7 + 0][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] +
616 al_1_prime[i * KFitConst::kNumber7 + 1][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] +
617 al_1_prime[i * KFitConst::kNumber7 + 2][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] +
618 m_property[i][1] * m_property[i][1]);
619 Sum_al_1[6][0] = + a;
620 }
621
622 for (int i = 0; i < m_TrackCount; i++) {
623 if (energy[i] == 0) {
625 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
626 break;
627 }
628
629 if (m_IsFixMass[i]) {
630 double invE = 1. / energy[i];
631 Sum_al_1[3][0] += energy[i];
632 Sum_al_1[4][0] += al_1_prime[i * KFitConst::kNumber7 + 1][0] * m_property[i][2] * invE;
633 Sum_al_1[5][0] += al_1_prime[i * KFitConst::kNumber7 + 0][0] * m_property[i][2] * invE;
634 } else {
635 Sum_al_1[3][0] += al_1_prime[i * KFitConst::kNumber7 + 3][0];
636 }
637
638 for (int j = 0; j < 3; j++) Sum_al_1[j][0] += al_1_prime[i * KFitConst::kNumber7 + j][0];
639 }
640
641 m_d[0][0] =
642 + Sum_al_1[3][0] * Sum_al_1[3][0] - Sum_al_1[0][0] * Sum_al_1[0][0]
643 - Sum_al_1[1][0] * Sum_al_1[1][0] - Sum_al_1[2][0] * Sum_al_1[2][0]
645
646 for (int i = 0; i < m_TrackCount; i++) {
647 if (energy[i] == 0) {
649 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
650 break;
651 }
652
653 a = m_property[i][2];
654
655 if (m_IsFixMass[i]) {
656 double invE = 1. / energy[i];
657 m_D[0][i * KFitConst::kNumber7 + 0] = 2.*(Sum_al_1[3][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE - Sum_al_1[0][0]);
658 m_D[0][i * KFitConst::kNumber7 + 1] = 2.*(Sum_al_1[3][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE - Sum_al_1[1][0]);
659 m_D[0][i * KFitConst::kNumber7 + 2] = 2.*(Sum_al_1[3][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] * invE - Sum_al_1[2][0]);
660 m_D[0][i * KFitConst::kNumber7 + 3] = 0.;
661 m_D[0][i * KFitConst::kNumber7 + 4] = -2.*(Sum_al_1[3][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE - Sum_al_1[1][0]) * a;
662 m_D[0][i * KFitConst::kNumber7 + 5] = 2.*(Sum_al_1[3][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE - Sum_al_1[0][0]) * a;
663 m_D[0][i * KFitConst::kNumber7 + 6] = 0.;
664 } else {
665 m_D[0][i * KFitConst::kNumber7 + 0] = -2.*Sum_al_1[0][0];
666 m_D[0][i * KFitConst::kNumber7 + 1] = -2.*Sum_al_1[1][0];
667 m_D[0][i * KFitConst::kNumber7 + 2] = -2.*Sum_al_1[2][0];
668 m_D[0][i * KFitConst::kNumber7 + 3] = 2.*Sum_al_1[3][0];
669 m_D[0][i * KFitConst::kNumber7 + 4] = 2.*Sum_al_1[1][0] * a;
670 m_D[0][i * KFitConst::kNumber7 + 5] = -2.*Sum_al_1[0][0] * a;
671 m_D[0][i * KFitConst::kNumber7 + 6] = 0.;
672 }
673 }
674
675 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]);
676 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]);
677 m_D[0][KFitConst::kNumber7 * m_TrackCount + 2] = 0.;
678 }
679
681}
682
683
686 m_NDF = 1;
687
689}
690
692{
693 MakeMotherKFit kmm;
695 unsigned n = getTrackCount();
696 for (unsigned i = 0; i < n; ++i) {
698 getTrack(i).getCharge());
701 for (unsigned j = i + 1; j < n; ++j) {
703 }
704 }
705 kmm.setVertex(getVertex());
708 m_ErrorCode = kmm.doMake();
710 return m_ErrorCode;
711 double chi2 = getCHIsq();
712 int ndf = getNDF();
713 double prob = TMath::Prob(chi2, ndf);
714 //
715 bool haschi2 = mother->hasExtraInfo("chiSquared");
716 if (haschi2) {
717 mother->setExtraInfo("chiSquared", chi2);
718 mother->setExtraInfo("ndf", ndf);
719 } else {
720 mother->addExtraInfo("chiSquared", chi2);
721 mother->addExtraInfo("ndf", ndf);
722 }
723
724 mother->updateMomentum(
725 CLHEPToROOT::getLorentzVector(kmm.getMotherMomentum()),
726 CLHEPToROOT::getXYZVector(kmm.getMotherPosition()),
727 CLHEPToROOT::getTMatrixFSym(kmm.getMotherError()),
728 prob);
730 return m_ErrorCode;
731}
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
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.
enum KFitError::ECode setZeroCorrelation(void) override
Indicate no correlation between tracks.
Definition: MassFitKFit.cc:131
enum KFitError::ECode setVertex(const HepPoint3D &v)
Set an initial vertex position for the mass-constraint fit.
Definition: MassFitKFit.cc:44
bool m_FlagTrackVertexError
Flag to indicate if the vertex error matrix of the child particle is preset.
Definition: MassFitKFit.h:156
bool m_FlagAtDecayPoint
Flag controlled by setFlagAtDecayPoint().
Definition: MassFitKFit.h:160
bool getFlagAtDecayPoint(void) const
Get a flag if to constraint at the decay point in the mass-constraint fit.
Definition: MassFitKFit.cc:179
enum KFitError::ECode prepareInputMatrix(void) override
Build grand matrices for minimum search from input-track properties.
Definition: MassFitKFit.cc:286
bool getFlagFitWithVertex(void) const
Get a flag if the fit is allowed with moving the vertex position.
Definition: MassFitKFit.cc:186
enum KFitError::ECode calculateNDF(void) override
Calculate an NDF of the fit.
Definition: MassFitKFit.cc:685
double getCHIsq(void) const override
Get a chi-square of the fit.
Definition: MassFitKFit.cc:193
enum KFitError::ECode setFlagAtDecayPoint(const bool flag)
Set a flag if to constraint at the decay point in the mass-constraint fit.
Definition: MassFitKFit.cc:76
enum KFitError::ECode setTrackZeroVertexError(void)
Indicate no vertex uncertainty in the child particle in the addTrack'ed order.
Definition: MassFitKFit.cc:117
std::vector< int > m_IsFixMass
Array of flags whether the track property is fixed at the mass.
Definition: MassFitKFit.h:165
enum KFitError::ECode updateMother(Particle *mother)
Update mother particle.
Definition: MassFitKFit.cc:691
enum KFitError::ECode unfixMass(void)
Tell the object to unfix the last added track property at the invariant mass.
Definition: MassFitKFit.cc:92
enum KFitError::ECode prepareInputSubMatrix(void) override
Build sub-matrices for minimum search from input-track properties.
Definition: MassFitKFit.cc:405
bool m_FlagFitIncludingVertex
Flag to indicate if the fit is allowed with moving the vertex position.
Definition: MassFitKFit.h:158
enum KFitError::ECode doFit(void)
Perform a mass-constraint fit.
Definition: MassFitKFit.cc:280
double m_InvariantMass
Invariant mass.
Definition: MassFitKFit.h:163
enum KFitError::ECode setInvariantMass(const double m)
Set an invariant mass for the mass-constraint fit.
Definition: MassFitKFit.cc:68
std::vector< CLHEP::HepMatrix > m_BeforeTrackVertexError
array of vertex error matrices before the fit.
Definition: MassFitKFit.h:146
enum KFitError::ECode makeCoreMatrix(void) override
Build matrices using the kinematical constraint.
Definition: MassFitKFit.cc:537
~MassFitKFit(void)
Destruct the object.
enum KFitError::ECode setVertexError(const CLHEP::HepSymMatrix &e)
Set an initial vertex error matrix for the mass-constraint fit.
Definition: MassFitKFit.cc:52
enum KFitError::ECode setTrackVertexError(const CLHEP::HepMatrix &e)
Set a vertex error matrix of the child particle in the addTrack'ed order.
Definition: MassFitKFit.cc:100
enum KFitError::ECode prepareCorrelation(void) override
Build a grand correlation matrix from input-track properties.
Definition: MassFitKFit.cc:416
HepPoint3D m_AfterVertex
Vertex position after the fit.
Definition: MassFitKFit.h:149
MassFitKFit(void)
Construct an object with no argument.
Definition: MassFitKFit.cc:25
enum KFitError::ECode setCorrelation(const CLHEP::HepMatrix &m) override
Set a correlation matrix.
Definition: MassFitKFit.cc:125
std::vector< CLHEP::HepMatrix > m_AfterTrackVertexError
array of vertex error matrices after the fit.
Definition: MassFitKFit.h:153
const CLHEP::HepMatrix getCorrelation(const int id1, const int id2, const int flag=KFitConst::kAfterFit) const override
Get a correlation matrix between two tracks.
Definition: MassFitKFit.cc:253
double getTrackCHIsq(const int id) const override
Get a chi-square of the track.
Definition: MassFitKFit.cc:217
const CLHEP::HepSymMatrix getVertexError(const int flag=KFitConst::kAfterFit) const
Get a vertex error matrix.
Definition: MassFitKFit.cc:156
const HepPoint3D getVertex(const int flag=KFitConst::kAfterFit) const
Get a vertex position.
Definition: MassFitKFit.cc:137
const CLHEP::HepMatrix getTrackVertexError(const int id, const int flag=KFitConst::kAfterFit) const
Get a vertex error matrix of the track.
Definition: MassFitKFit.cc:200
CLHEP::HepSymMatrix m_AfterVertexError
Vertex error matrix after the fit.
Definition: MassFitKFit.h:151
enum KFitError::ECode prepareOutputMatrix(void) override
Build an output error matrix.
Definition: MassFitKFit.cc:474
CLHEP::HepSymMatrix m_BeforeVertexError
Vertex error matrix before the fit.
Definition: MassFitKFit.h:144
HepPoint3D m_BeforeVertex
Vertex position before the fit.
Definition: MassFitKFit.h:142
enum KFitError::ECode fixMass(void)
Tell the object to fix the last added track property at the invariant mass.
Definition: MassFitKFit.cc:84
double getInvariantMass(void) const
Get an invariant mass.
Definition: MassFitKFit.cc:172
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