Belle II Software prerelease-10-00-00a
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
73
74
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
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
183
184
185bool
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
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
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: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
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.
enum KFitError::ECode setVertex(const HepPoint3D &v)
Set an initial vertex position for the mass-constraint fit.
bool m_FlagTrackVertexError
Flag to indicate if the vertex error matrix of the child particle is preset.
bool m_FlagAtDecayPoint
Flag controlled by setFlagAtDecayPoint().
bool getFlagAtDecayPoint(void) const
Get a flag if to constraint at the decay point in the mass-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 mass-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.
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 mass-constraint fit.
double m_InvariantMass
Invariant mass.
enum KFitError::ECode setInvariantMass(const double m)
Set an invariant mass for the mass-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.
~MassFitKFit(void)
Destruct the object.
enum KFitError::ECode setVertexError(const CLHEP::HepSymMatrix &e)
Set an initial vertex error matrix for the mass-constraint fit.
enum KFitError::ECode setTrackVertexError(const CLHEP::HepMatrix &e)
Set a vertex error matrix of the child particle in the addTrack'ed order.
enum KFitError::ECode prepareCorrelation(void) override
Build a grand correlation matrix from input-track properties.
HepPoint3D m_AfterVertex
Vertex position after the fit.
MassFitKFit(void)
Construct an object with no argument.
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.
double getInvariantMass(void) const
Get an invariant mass.
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: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