Belle II Software development
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 <analysis/VertexFitting/KFit/MassFitKFit.h>
11
12#include <analysis/dataobjects/Particle.h>
13#include <analysis/VertexFitting/KFit/MakeMotherKFit.h>
14#include <analysis/utility/CLHEPToROOT.h>
15#include <framework/gearbox/Const.h>
16
17#include <cstdio>
18
19#include <TMath.h>
20
21using namespace std;
22using namespace Belle2;
23using namespace Belle2::analysis;
24using namespace CLHEP;
25
27{
28 m_FlagFitted = false;
31 m_FlagAtDecayPoint = true;
33 m_d = HepMatrix(1, 1, 0);
34 m_V_D = HepMatrix(1, 1, 0);
35 m_lam = HepMatrix(1, 1, 0);
36 m_AfterVertexError = HepSymMatrix(3, 0);
37 m_InvariantMass = -1.0;
38}
39
40
42
43
45MassFitKFit::setVertex(const HepPoint3D& v) {
47
49}
50
51
53MassFitKFit::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
74
75
82
83
86 m_IsFixMass.push_back(true);
87
89}
90
91
94 m_IsFixMass.push_back(false);
95
97}
98
99
102 if (e.num_row() != 3 || e.num_col() != KFitConst::kNumber7)
103 {
105 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
106 return m_ErrorCode;
107 }
108
109 m_BeforeTrackVertexError.push_back(e);
112
114}
115
116
119 HepMatrix zero(3, KFitConst::kNumber7, 0);
120
121 return this->setTrackVertexError(zero);
122}
123
124
126MassFitKFit::setCorrelation(const HepMatrix& m) {
127 return KFitBase::setCorrelation(m);
128}
129
130
135
136
137const HepPoint3D
138MassFitKFit::getVertex(const int flag) const
139{
140 if (flag == KFitConst::kAfterFit && !isFitted()) return HepPoint3D();
141
142 switch (flag) {
144 return m_BeforeVertex;
145
147 return m_AfterVertex;
148
149 default:
150 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
151 return HepPoint3D();
152 }
153}
154
155
156const HepSymMatrix
157MassFitKFit::getVertexError(const int flag) const
158{
159 if (flag == KFitConst::kAfterFit && !isFitted()) return HepSymMatrix(3, 0);
160
161 if (flag == KFitConst::kBeforeFit)
162 return m_BeforeVertexError;
164 return m_AfterVertexError;
165 else {
166 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
167 return HepSymMatrix(3, 0);
168 }
169}
170
171
172double
174{
175 return m_InvariantMass;
176}
177
178
179bool
184
185
186bool
191
192
193double
195{
196 return m_CHIsq;
197}
198
199
200const HepMatrix
201MassFitKFit::getTrackVertexError(const int id, const int flag) const
202{
203 if (flag == KFitConst::kAfterFit && !isFitted()) return HepMatrix(3, KFitConst::kNumber7, 0);
204 if (!isTrackIDInRange(id)) return HepMatrix(3, KFitConst::kNumber7, 0);
205
206 if (flag == KFitConst::kBeforeFit)
207 return m_BeforeTrackVertexError[id];
209 return m_AfterTrackVertexError[id];
210 else {
211 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
212 return HepMatrix(3, KFitConst::kNumber7, 0);
213 }
214}
215
216
217double
218MassFitKFit::getTrackCHIsq(const int id) const
219{
220 if (!isFitted()) return -1;
221 if (!isTrackIDInRange(id)) return -1;
222
223 if (m_IsFixMass[id]) {
224
225 HepMatrix da(m_Tracks[id].getFitParameter(KFitConst::kBeforeFit) - m_Tracks[id].getFitParameter(KFitConst::kAfterFit));
226 int err_inverse = 0;
227 const double chisq = (da.T() * (m_Tracks[id].getFitError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
228
229 if (err_inverse) {
230 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
231 return -1;
232 }
233
234 return chisq;
235
236 } else {
237
238 HepMatrix da(m_Tracks[id].getMomPos(KFitConst::kBeforeFit) - m_Tracks[id].getMomPos(KFitConst::kAfterFit));
239 int err_inverse = 0;
240 const double chisq = (da.T() * (m_Tracks[id].getError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
241
242 if (err_inverse) {
243 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
244 return -1;
245 }
246
247 return chisq;
248
249 }
250}
251
252
253const HepMatrix
254MassFitKFit::getCorrelation(const int id1, const int id2, const int flag) const
255{
256 if (flag == KFitConst::kAfterFit && !isFitted()) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
257 if (!isTrackIDInRange(id1)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
258 if (!isTrackIDInRange(id2)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
259
260 switch (flag) {
262 return KFitBase::getCorrelation(id1, id2, flag);
263
265 return makeError3(
266 this->getTrackMomentum(id1),
267 this->getTrackMomentum(id2),
268 m_V_al_1.sub(KFitConst::kNumber7 * id1 + 1, KFitConst::kNumber7 * (id1 + 1), KFitConst::kNumber7 * id2 + 1,
269 KFitConst::kNumber7 * (id2 + 1)),
270 m_IsFixMass[id1],
271 m_IsFixMass[id2]);
272
273 default:
274 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
275 return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
276 }
277}
278
279
284
285
289 {
291 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
292 return m_ErrorCode;
293 }
294
295
296 if (m_IsFixMass.size() == 0)
297 {
298 // If no fix_mass flag at all,
299 // all tracks are considered to be fixed at mass.
300 for (int i = 0; i < m_TrackCount; i++) this->fixMass();
301 } else if (m_IsFixMass.size() != (unsigned int)m_TrackCount)
302 {
304 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
305 return m_ErrorCode;
306 }
307
308
310 {
311 int index = 0;
312 m_al_0 = HepMatrix(KFitConst::kNumber7 * m_TrackCount, 1, 0);
313 m_property = HepMatrix(m_TrackCount, 3, 0);
314 m_V_al_0 = HepSymMatrix(KFitConst::kNumber7 * m_TrackCount, 0);
315
316 for (auto& track : m_Tracks) {
317 // momentum x,y,z and position x,y,z
318 m_al_0[index * KFitConst::kNumber7 + 0][0] = track.getMomentum(KFitConst::kBeforeFit).x();
319 m_al_0[index * KFitConst::kNumber7 + 1][0] = track.getMomentum(KFitConst::kBeforeFit).y();
320 m_al_0[index * KFitConst::kNumber7 + 2][0] = track.getMomentum(KFitConst::kBeforeFit).z();
321 m_al_0[index * KFitConst::kNumber7 + 3][0] = track.getMomentum(KFitConst::kBeforeFit).t();
322 m_al_0[index * KFitConst::kNumber7 + 4][0] = track.getPosition(KFitConst::kBeforeFit).x();
323 m_al_0[index * KFitConst::kNumber7 + 5][0] = track.getPosition(KFitConst::kBeforeFit).y();
324 m_al_0[index * KFitConst::kNumber7 + 6][0] = track.getPosition(KFitConst::kBeforeFit).z();
325 // these error
326 m_V_al_0.sub(index * KFitConst::kNumber7 + 1, track.getError(KFitConst::kBeforeFit));
327 // charge, mass, a
328 m_property[index][0] = track.getCharge();
329 m_property[index][1] = track.getMass();
330 const double c = Const::speedOfLight * 1e-4;
331 m_property[index][2] = -c * m_MagneticField * track.getCharge();
332 index++;
333 }
334
335 // error between track and track
336 if (m_FlagCorrelation) {
337 this->prepareCorrelation();
339 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
340 return m_ErrorCode;
341 }
342 }
343
344 // set member matrix
345 m_al_1 = m_al_0;
346
347 // define size of matrix
350
351 } else
352 {
353 // m_FlagFitIncludingVertex == true
354 int index = 0;
355 m_al_0 = HepMatrix(KFitConst::kNumber7 * m_TrackCount + 3, 1, 0);
356 m_property = HepMatrix(m_TrackCount, 3, 0);
357 m_V_al_0 = HepSymMatrix(KFitConst::kNumber7 * m_TrackCount + 3, 0);
358
359 for (auto& track : m_Tracks) {
360 // momentum x,y,z and position x,y,z
361 m_al_0[index * KFitConst::kNumber7 + 0][0] = track.getMomentum(KFitConst::kBeforeFit).x();
362 m_al_0[index * KFitConst::kNumber7 + 1][0] = track.getMomentum(KFitConst::kBeforeFit).y();
363 m_al_0[index * KFitConst::kNumber7 + 2][0] = track.getMomentum(KFitConst::kBeforeFit).z();
364 m_al_0[index * KFitConst::kNumber7 + 3][0] = track.getMomentum(KFitConst::kBeforeFit).t();
365 m_al_0[index * KFitConst::kNumber7 + 4][0] = track.getPosition(KFitConst::kBeforeFit).x();
366 m_al_0[index * KFitConst::kNumber7 + 5][0] = track.getPosition(KFitConst::kBeforeFit).y();
367 m_al_0[index * KFitConst::kNumber7 + 6][0] = track.getPosition(KFitConst::kBeforeFit).z();
368 // these error
369 m_V_al_0.sub(index * KFitConst::kNumber7 + 1, track.getError(KFitConst::kBeforeFit));
370 // charge, mass, a
371 m_property[index][0] = track.getCharge();
372 m_property[index][1] = track.getMass();
373 const double c = Const::speedOfLight * 1e-4;
374 m_property[index][2] = -c * m_MagneticField * track.getCharge();
375 index++;
376 }
377
378 // vertex
383
384 // error between track and track
385 if (m_FlagCorrelation) {
386 this->prepareCorrelation();
388 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
389 return m_ErrorCode;
390 }
391 }
392
393 // set member matrix
394 m_al_1 = m_al_0;
395
396 // define size of matrix
398 m_D = m_V_al_1.sub(1, 1, 1, KFitConst::kNumber7 * m_TrackCount + 3);
399 }
400
402}
403
404
407 char buf[1024];
408 sprintf(buf, "%s:%s(): internal error; this function should never be called", __FILE__, __func__);
409 B2FATAL(buf);
410
411 /* NEVER REACHEd HERE */
413}
414
415
418 if (m_BeforeCorrelation.size() != static_cast<unsigned int>(m_TrackCount * (m_TrackCount - 1) / 2))
419 {
421 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
422 return m_ErrorCode;
423 }
424
425 int row = 0, col = 0;
426
427 for (auto& hm : m_BeforeCorrelation)
428 {
429 // counter
430 row++;
431 if (row == m_TrackCount) {
432 col++;
433 row = col + 1;
434 }
435
436 int ii = 0, jj = 0;
437 for (int i = KFitConst::kNumber7 * row; i < KFitConst::kNumber7 * (row + 1); i++) {
438 for (int j = KFitConst::kNumber7 * col; j < KFitConst::kNumber7 * (col + 1); j++) {
439 m_V_al_0[i][j] = hm[ii][jj];
440 jj++;
441 }
442 jj = 0;
443 ii++;
444 }
445 }
446
448 {
449 // ...error of vertex
451
452 // ...error matrix between vertex and tracks
454 if (m_BeforeTrackVertexError.size() != (unsigned int)m_TrackCount) {
456 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
457 return m_ErrorCode;
458 }
459
460 int i = 0;
461 for (auto& hm : m_BeforeTrackVertexError) {
462 for (int j = 0; j < 3; j++) for (int k = 0; k < KFitConst::kNumber7; k++) {
464 }
465 i++;
466 }
467 }
468 }
469
471}
472
473
476 Hep3Vector h3v;
477 int index = 0;
478 for (auto& pdata : m_Tracks)
479 {
480 // tracks
481 // momentum
482 h3v.setX(m_al_1[index * KFitConst::kNumber7 + 0][0]);
483 h3v.setY(m_al_1[index * KFitConst::kNumber7 + 1][0]);
484 h3v.setZ(m_al_1[index * KFitConst::kNumber7 + 2][0]);
485 if (m_IsFixMass[index])
486 pdata.setMomentum(HepLorentzVector(h3v, sqrt(h3v.mag2() + pdata.getMass()*pdata.getMass())), KFitConst::kAfterFit);
487 else
488 pdata.setMomentum(HepLorentzVector(h3v, m_al_1[index * KFitConst::kNumber7 + 3][0]), KFitConst::kAfterFit);
489 // position
490 pdata.setPosition(HepPoint3D(
491 m_al_1[index * KFitConst::kNumber7 + 4][0],
492 m_al_1[index * KFitConst::kNumber7 + 5][0],
494 // error of the tracks
495 pdata.setError(this->makeError3(pdata.getMomentum(),
496 m_V_al_1.sub(
497 index * KFitConst::kNumber7 + 1,
498 (index + 1)*KFitConst::kNumber7,
499 index * KFitConst::kNumber7 + 1,
500 (index + 1)*KFitConst::kNumber7), m_IsFixMass[index]),
502 if (m_ErrorCode != KFitError::kNoError) break;
503 index++;
504 }
505
507 {
508 // vertex
512 // error of the vertex
513 for (int i = 0; i < 3; i++) for (int j = i; j < 3; j++) {
515 }
516 // error between vertex and tracks
517 for (int i = 0; i < m_TrackCount; i++) {
518 HepMatrix hm(3, KFitConst::kNumber7, 0);
519 for (int j = 0; j < 3; j++) for (int k = 0; k < KFitConst::kNumber7; k++) {
521 }
522 if (m_IsFixMass[i])
523 m_AfterTrackVertexError.push_back(this->makeError4(m_Tracks[i].getMomentum(), hm));
524 else
525 m_AfterTrackVertexError.push_back(hm);
526 }
527 } else
528 {
529 // not fit
531 }
532
534}
535
536
540 {
541
542 HepMatrix al_1_prime(m_al_1);
543 HepMatrix Sum_al_1(4, 1, 0);
544 std::vector<double> energy(m_TrackCount);
545 double a;
546
547 for (int i = 0; i < m_TrackCount; i++) {
548 a = m_property[i][2];
549 if (!m_FlagAtDecayPoint) a = 0.;
550 al_1_prime[i * KFitConst::kNumber7 + 0][0] -= a * (m_BeforeVertex.y() - al_1_prime[i * KFitConst::kNumber7 + 5][0]);
551 al_1_prime[i * KFitConst::kNumber7 + 1][0] += a * (m_BeforeVertex.x() - al_1_prime[i * KFitConst::kNumber7 + 4][0]);
552 energy[i] = sqrt(al_1_prime[i * KFitConst::kNumber7 + 0][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] +
553 al_1_prime[i * KFitConst::kNumber7 + 1][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] +
554 al_1_prime[i * KFitConst::kNumber7 + 2][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] +
555 m_property[i][1] * m_property[i][1]);
556 if (m_IsFixMass[i])
557 Sum_al_1[3][0] += energy[i];
558 else
559 Sum_al_1[3][0] += al_1_prime[i * KFitConst::kNumber7 + 3][0];
560 }
561
562 for (int i = 0; i < m_TrackCount; i++) {
563 for (int j = 0; j < 3; j++) Sum_al_1[j][0] += al_1_prime[i * KFitConst::kNumber7 + j][0];
564 }
565
566 m_d[0][0] =
567 + Sum_al_1[3][0] * Sum_al_1[3][0] - Sum_al_1[0][0] * Sum_al_1[0][0]
568 - Sum_al_1[1][0] * Sum_al_1[1][0] - Sum_al_1[2][0] * Sum_al_1[2][0]
570
571 for (int i = 0; i < m_TrackCount; i++) {
572 if (energy[i] == 0) {
574 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
575 break;
576 }
577
578 a = m_property[i][2];
579 if (!m_FlagAtDecayPoint) a = 0.;
580
581 if (m_IsFixMass[i]) {
582 double invE = 1. / energy[i];
583 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]);
584 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]);
585 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]);
586 m_D[0][i * KFitConst::kNumber7 + 3] = 0.;
587 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;
588 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;
589 m_D[0][i * KFitConst::kNumber7 + 6] = 0.;
590 } else {
591 m_D[0][i * KFitConst::kNumber7 + 0] = -2.*Sum_al_1[0][0];
592 m_D[0][i * KFitConst::kNumber7 + 1] = -2.*Sum_al_1[1][0];
593 m_D[0][i * KFitConst::kNumber7 + 2] = -2.*Sum_al_1[2][0];
594 m_D[0][i * KFitConst::kNumber7 + 3] = 2.*Sum_al_1[3][0];
595 m_D[0][i * KFitConst::kNumber7 + 4] = 2.*Sum_al_1[1][0] * a;
596 m_D[0][i * KFitConst::kNumber7 + 5] = -2.*Sum_al_1[0][0] * a;
597 m_D[0][i * KFitConst::kNumber7 + 6] = 0.;
598 }
599 }
600
601 } else
602 {
603
604 // m_FlagFitIncludingVertex == true
605 HepMatrix al_1_prime(m_al_1);
606 HepMatrix Sum_al_1(7, 1, 0);
607 std::vector<double> energy(m_TrackCount);
608 double a;
609
610 for (int i = 0; i < m_TrackCount; i++) {
611 a = m_property[i][2];
612 al_1_prime[i * KFitConst::kNumber7 + 0][0] -= a * (al_1_prime[KFitConst::kNumber7 * m_TrackCount + 1][0] - al_1_prime[i *
613 KFitConst::kNumber7 + 5][0]);
614 al_1_prime[i * KFitConst::kNumber7 + 1][0] += a * (al_1_prime[KFitConst::kNumber7 * m_TrackCount + 0][0] - al_1_prime[i *
615 KFitConst::kNumber7 + 4][0]);
616 energy[i] = sqrt(al_1_prime[i * KFitConst::kNumber7 + 0][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] +
617 al_1_prime[i * KFitConst::kNumber7 + 1][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] +
618 al_1_prime[i * KFitConst::kNumber7 + 2][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] +
619 m_property[i][1] * m_property[i][1]);
620 Sum_al_1[6][0] = + a;
621 }
622
623 for (int i = 0; i < m_TrackCount; i++) {
624 if (energy[i] == 0) {
626 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
627 break;
628 }
629
630 if (m_IsFixMass[i]) {
631 double invE = 1. / energy[i];
632 Sum_al_1[3][0] += energy[i];
633 Sum_al_1[4][0] += al_1_prime[i * KFitConst::kNumber7 + 1][0] * m_property[i][2] * invE;
634 Sum_al_1[5][0] += al_1_prime[i * KFitConst::kNumber7 + 0][0] * m_property[i][2] * invE;
635 } else {
636 Sum_al_1[3][0] += al_1_prime[i * KFitConst::kNumber7 + 3][0];
637 }
638
639 for (int j = 0; j < 3; j++) Sum_al_1[j][0] += al_1_prime[i * KFitConst::kNumber7 + j][0];
640 }
641
642 m_d[0][0] =
643 + Sum_al_1[3][0] * Sum_al_1[3][0] - Sum_al_1[0][0] * Sum_al_1[0][0]
644 - Sum_al_1[1][0] * Sum_al_1[1][0] - Sum_al_1[2][0] * Sum_al_1[2][0]
646
647 for (int i = 0; i < m_TrackCount; i++) {
648 if (energy[i] == 0) {
650 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
651 break;
652 }
653
654 a = m_property[i][2];
655
656 if (m_IsFixMass[i]) {
657 double invE = 1. / energy[i];
658 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]);
659 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]);
660 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]);
661 m_D[0][i * KFitConst::kNumber7 + 3] = 0.;
662 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;
663 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;
664 m_D[0][i * KFitConst::kNumber7 + 6] = 0.;
665 } else {
666 m_D[0][i * KFitConst::kNumber7 + 0] = -2.*Sum_al_1[0][0];
667 m_D[0][i * KFitConst::kNumber7 + 1] = -2.*Sum_al_1[1][0];
668 m_D[0][i * KFitConst::kNumber7 + 2] = -2.*Sum_al_1[2][0];
669 m_D[0][i * KFitConst::kNumber7 + 3] = 2.*Sum_al_1[3][0];
670 m_D[0][i * KFitConst::kNumber7 + 4] = 2.*Sum_al_1[1][0] * a;
671 m_D[0][i * KFitConst::kNumber7 + 5] = -2.*Sum_al_1[0][0] * a;
672 m_D[0][i * KFitConst::kNumber7 + 6] = 0.;
673 }
674 }
675
676 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]);
677 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]);
678 m_D[0][KFitConst::kNumber7 * m_TrackCount + 2] = 0.;
679 }
680
682}
683
684
691
693{
694 MakeMotherKFit kmm;
696 unsigned n = getTrackCount();
697 for (unsigned i = 0; i < n; ++i) {
699 getTrack(i).getCharge());
702 for (unsigned j = i + 1; j < n; ++j) {
704 }
705 }
706 kmm.setVertex(getVertex());
709 m_ErrorCode = kmm.doMake();
711 return m_ErrorCode;
712 double chi2 = getCHIsq();
713 int ndf = getNDF();
714 double prob = TMath::Prob(chi2, ndf);
715 //
716 bool haschi2 = mother->hasExtraInfo("chiSquared");
717 if (haschi2) {
718 mother->setExtraInfo("chiSquared", chi2);
719 mother->setExtraInfo("ndf", ndf);
720 } else {
721 mother->addExtraInfo("chiSquared", chi2);
722 mother->addExtraInfo("ndf", ndf);
723 }
724
725 mother->updateMomentum(
726 CLHEPToROOT::getLorentzVector(kmm.getMotherMomentum()),
727 CLHEPToROOT::getXYZVector(kmm.getMotherPosition()),
728 CLHEPToROOT::getTMatrixFSym(kmm.getMotherError()),
729 prob);
731 return m_ErrorCode;
732}
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:71
ECode
ECode is a error code enumerate.
Definition KFitError.h:33
@ kCannotGetMatrixInverse
Cannot calculate matrix inverse (bad track property or internal error)
Definition KFitError.h:57
@ kOutOfRange
Specified track-id out of range.
Definition KFitError.h:41
@ kDivisionByZero
Division by zero (bad track property or internal error)
Definition KFitError.h:55
@ kBadTrackSize
Track count too small to perform fit.
Definition KFitError.h:46
@ kBadMatrixSize
Wrong correlation matrix size.
Definition KFitError.h:48
@ kBadCorrelationSize
Wrong correlation matrix size (internal error)
Definition KFitError.h:50
MakeMotherKFit is a class to build mother particle from kinematically fitted daughters.
enum KFitError::ECode setVertex(const HepPoint3D &v)
Set a vertex position of the mother particle.
enum KFitError::ECode addTrack(const KFitTrack &kp)
Add a track to the make-mother object.
enum KFitError::ECode doMake(void)
Perform a reconstruction of mother particle.
const CLHEP::HepSymMatrix getMotherError(void) const
Get an error matrix of the mother particle.
enum KFitError::ECode setCorrelation(const CLHEP::HepMatrix &e)
Set a correlation matrix.
const HepPoint3D getMotherPosition(void) const
Get a position of the mother particle.
enum KFitError::ECode setVertexError(const CLHEP::HepSymMatrix &e)
Set a vertex error matrix of the mother particle.
enum KFitError::ECode setTrackVertexError(const CLHEP::HepMatrix &e)
Set a vertex error matrix of the child particle in the addTrack'ed order.
const CLHEP::HepLorentzVector getMotherMomentum(void) const
Get a Lorentz vector of the mother particle.
enum KFitError::ECode setMagneticField(const double mf)
Change a magnetic field from the default value KFitConst::kDefaultMagneticField.
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:38
static const int kAfterFit
Input parameter to specify after-fit when setting/getting a track attribute.
Definition KFitConst.h:35
static const int kBeforeFit
Input parameter to specify before-fit when setting/getting a track attribute.
Definition KFitConst.h:33
static const int kNumber7
Constant 7 to check matrix size (internal use)
Definition KFitConst.h:30