Belle II Software light-2406-ragdoll
MassFourCFitKFit.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/MassFourCFitKFit.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
25MassFourCFitKFit::MassFourCFitKFit() : m_AfterVertexError(HepSymMatrix(3, 0)),
26 m_FourMomentum(PxPyPzEVector())
27{
28 m_FlagFitted = false;
31 m_FlagAtDecayPoint = true;
33 m_InvariantMass = -1.0;
37}
38
39
41
43MassFourCFitKFit::addMassConstraint(const double m, std::vector<unsigned>& childTrackId) {
44 if ((childTrackId.back() - childTrackId.front()) != childTrackId.size() - 1)
45 {
47 }
48 m_ConstraintMassList.push_back(m);
49 m_ConstraintMassChildLists.push_back(std::make_pair(childTrackId.front(), childTrackId.back()));
51}
52
54MassFourCFitKFit::setVertex(const HepPoint3D& v) {
56
58}
59
60
62MassFourCFitKFit::setVertexError(const HepSymMatrix& e) {
63 if (e.num_row() != 3)
64 {
66 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
67 return m_ErrorCode;
68 }
69
72
74}
75
76
80
82}
83
84
86MassFourCFitKFit::setFourMomentum(const PxPyPzEVector& m) {
88
90}
91
92
95 m_FlagAtDecayPoint = flag;
96
98}
99
100
103 m_IsFixMass.push_back(true);
104
106}
107
108
111 m_IsFixMass.push_back(false);
112
114}
115
116
119 if (e.num_row() != 3 || e.num_col() != KFitConst::kNumber7)
120 {
122 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
123 return m_ErrorCode;
124 }
125
126 m_BeforeTrackVertexError.push_back(e);
129
131}
132
133
136 HepMatrix zero(3, KFitConst::kNumber7, 0);
137
138 return this->setTrackVertexError(zero);
139}
140
141
144 return KFitBase::setCorrelation(m);
145}
146
147
151}
152
153
154const HepPoint3D
155MassFourCFitKFit::getVertex(const int flag) const
156{
157 if (flag == KFitConst::kAfterFit && !isFitted()) return HepPoint3D();
158
159 switch (flag) {
161 return m_BeforeVertex;
162
164 return m_AfterVertex;
165
166 default:
167 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
168 return HepPoint3D();
169 }
170}
171
172
173const HepSymMatrix
175{
176 if (flag == KFitConst::kAfterFit && !isFitted()) return HepSymMatrix(3, 0);
177
178 if (flag == KFitConst::kBeforeFit)
179 return m_BeforeVertexError;
181 return m_AfterVertexError;
182 else {
183 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
184 return HepSymMatrix(3, 0);
185 }
186}
187
188
189double
191{
192 return m_InvariantMass;
193}
194
195
196bool
198{
199 return m_FlagAtDecayPoint;
200}
201
202
203bool
205{
207}
208
209
210double
212{
213 return m_CHIsq;
214}
215
216
217const HepMatrix
218MassFourCFitKFit::getTrackVertexError(const int id, const int flag) const
219{
220 if (flag == KFitConst::kAfterFit && !isFitted()) return HepMatrix(3, KFitConst::kNumber7, 0);
221 if (!isTrackIDInRange(id)) return HepMatrix(3, KFitConst::kNumber7, 0);
222
223 if (flag == KFitConst::kBeforeFit)
224 return m_BeforeTrackVertexError[id];
226 return m_AfterTrackVertexError[id];
227 else {
228 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
229 return HepMatrix(3, KFitConst::kNumber7, 0);
230 }
231}
232
233
234double
236{
237 if (!isFitted()) return -1;
238 if (!isTrackIDInRange(id)) return -1;
239
240 if (m_IsFixMass[id]) {
241
242 HepMatrix da(m_Tracks[id].getFitParameter(KFitConst::kBeforeFit) - m_Tracks[id].getFitParameter(KFitConst::kAfterFit));
243 int err_inverse = 0;
244 const double chisq = (da.T() * (m_Tracks[id].getFitError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
245
246 if (err_inverse) {
247 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
248 return -1;
249 }
250
251 return chisq;
252
253 } else {
254
255 HepMatrix da(m_Tracks[id].getMomPos(KFitConst::kBeforeFit) - m_Tracks[id].getMomPos(KFitConst::kAfterFit));
256 int err_inverse = 0;
257 const double chisq = (da.T() * (m_Tracks[id].getError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
258
259 if (err_inverse) {
260 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
261 return -1;
262 }
263
264 return chisq;
265
266 }
267}
268
269
270const HepMatrix
271MassFourCFitKFit::getCorrelation(const int id1, const int id2, const int flag) const
272{
273 if (flag == KFitConst::kAfterFit && !isFitted()) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
274 if (!isTrackIDInRange(id1)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
275 if (!isTrackIDInRange(id2)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
276
277 switch (flag) {
279 return KFitBase::getCorrelation(id1, id2, flag);
280
282 return makeError3(
283 this->getTrackMomentum(id1),
284 this->getTrackMomentum(id2),
285 m_V_al_1.sub(KFitConst::kNumber7 * id1 + 1, KFitConst::kNumber7 * (id1 + 1), KFitConst::kNumber7 * id2 + 1,
286 KFitConst::kNumber7 * (id2 + 1)),
287 m_IsFixMass[id1],
288 m_IsFixMass[id2]);
289
290 default:
291 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
292 return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
293 }
294}
295
296
299 return KFitBase::doFit1();
300}
301
302
306 m_d = HepMatrix(4 + m_ConstraintMassCount, 1, 0);
307 m_V_D = HepMatrix(4 + m_ConstraintMassCount, 4 + m_ConstraintMassCount, 0);
308 m_lam = HepMatrix(4 + m_ConstraintMassCount, 1, 0);
309
311 {
313 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
314 return m_ErrorCode;
315 }
316
317 if (m_IsFixMass.size() == 0)
318 {
319 // If no fix_mass flag at all,
320 // all tracks are considered to be fixed at mass.
321 for (int i = 0; i < m_TrackCount; i++) this->fixMass();
322 } else if (m_IsFixMass.size() != (unsigned int)m_TrackCount)
323 {
325 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
326 return m_ErrorCode;
327 }
328
330 {
331 int index = 0;
332 m_al_0 = HepMatrix(KFitConst::kNumber7 * m_TrackCount, 1, 0);
333 m_property = HepMatrix(m_TrackCount, 3, 0);
334 m_V_al_0 = HepSymMatrix(KFitConst::kNumber7 * m_TrackCount, 0);
335
336 for (auto& track : m_Tracks) {
337 // momentum x,y,z and position x,y,z
338 m_al_0[index * KFitConst::kNumber7 + 0][0] = track.getMomentum(KFitConst::kBeforeFit).x();
339 m_al_0[index * KFitConst::kNumber7 + 1][0] = track.getMomentum(KFitConst::kBeforeFit).y();
340 m_al_0[index * KFitConst::kNumber7 + 2][0] = track.getMomentum(KFitConst::kBeforeFit).z();
341 m_al_0[index * KFitConst::kNumber7 + 3][0] = track.getMomentum(KFitConst::kBeforeFit).t();
342 m_al_0[index * KFitConst::kNumber7 + 4][0] = track.getPosition(KFitConst::kBeforeFit).x();
343 m_al_0[index * KFitConst::kNumber7 + 5][0] = track.getPosition(KFitConst::kBeforeFit).y();
344 m_al_0[index * KFitConst::kNumber7 + 6][0] = track.getPosition(KFitConst::kBeforeFit).z();
345 // these error
346 m_V_al_0.sub(index * KFitConst::kNumber7 + 1, track.getError(KFitConst::kBeforeFit));
347 // charge, mass, a
348 m_property[index][0] = track.getCharge();
349 m_property[index][1] = track.getMass();
350 const double c = Belle2::Const::speedOfLight * 1e-4;
351 m_property[index][2] = -c * m_MagneticField * track.getCharge();
352 index++;
353 }
354
355 // error between track and track
356 if (m_FlagCorrelation) {
357 this->prepareCorrelation();
359 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
360 return m_ErrorCode;
361 }
362 }
363
364 // set member matrix
365 m_al_1 = m_al_0;
366
367 // define size of matrix
370
371 } else
372 {
373 //TODO: Not Implemented
375 // m_FlagFitIncludingVertex == true
376 int index = 0;
377 m_al_0 = HepMatrix(KFitConst::kNumber7 * m_TrackCount + 3, 1, 0);
378 m_property = HepMatrix(m_TrackCount, 3, 0);
379 m_V_al_0 = HepSymMatrix(KFitConst::kNumber7 * m_TrackCount + 3, 0);
380
381 for (auto& track : m_Tracks) {
382 // momentum x,y,z and position x,y,z
383 m_al_0[index * KFitConst::kNumber7 + 0][0] = track.getMomentum(KFitConst::kBeforeFit).x();
384 m_al_0[index * KFitConst::kNumber7 + 1][0] = track.getMomentum(KFitConst::kBeforeFit).y();
385 m_al_0[index * KFitConst::kNumber7 + 2][0] = track.getMomentum(KFitConst::kBeforeFit).z();
386 m_al_0[index * KFitConst::kNumber7 + 3][0] = track.getMomentum(KFitConst::kBeforeFit).t();
387 m_al_0[index * KFitConst::kNumber7 + 4][0] = track.getPosition(KFitConst::kBeforeFit).x();
388 m_al_0[index * KFitConst::kNumber7 + 5][0] = track.getPosition(KFitConst::kBeforeFit).y();
389 m_al_0[index * KFitConst::kNumber7 + 6][0] = track.getPosition(KFitConst::kBeforeFit).z();
390 // these error
391 m_V_al_0.sub(index * KFitConst::kNumber7 + 1, track.getError(KFitConst::kBeforeFit));
392 // charge, mass, a
393 m_property[index][0] = track.getCharge();
394 m_property[index][1] = track.getMass();
395 const double c = Belle2::Const::speedOfLight * 1e-4;
396 m_property[index][2] = -c * m_MagneticField * track.getCharge();
397 index++;
398 }
399
400 // vertex
405
406 // error between track and track
407 if (m_FlagCorrelation) {
408 this->prepareCorrelation();
410 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
411 return m_ErrorCode;
412 }
413 }
414
415 // set member matrix
416 m_al_1 = m_al_0;
417
418 // define size of matrix
420 m_D = m_V_al_1.sub(1, 4, 1, KFitConst::kNumber7 * m_TrackCount + 3);
421 }
422
424}
425
426
429 char buf[1024];
430 sprintf(buf, "%s:%s(): internal error; this function should never be called", __FILE__, __func__);
431 B2FATAL(buf);
432
433 /* NEVER REACHEd HERE */
435}
436
437
440 if (m_BeforeCorrelation.size() != static_cast<unsigned int>(m_TrackCount * (m_TrackCount - 1) / 2))
441 {
443 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
444 return m_ErrorCode;
445 }
446
447 int row = 0, col = 0;
448
449 for (auto& hm : m_BeforeCorrelation)
450 {
451 // counter
452 row++;
453 if (row == m_TrackCount) {
454 col++;
455 row = col + 1;
456 }
457
458 int ii = 0, jj = 0;
459 for (int i = KFitConst::kNumber7 * row; i < KFitConst::kNumber7 * (row + 1); i++) {
460 for (int j = KFitConst::kNumber7 * col; j < KFitConst::kNumber7 * (col + 1); j++) {
461 m_V_al_0[i][j] = hm[ii][jj];
462 jj++;
463 }
464 jj = 0;
465 ii++;
466 }
467 }
468
470 {
471 //TODO: Not Implemented
473
474 // ...error of vertex
476
477 // ...error matrix between vertex and tracks
479 if (m_BeforeTrackVertexError.size() != (unsigned int)m_TrackCount) {
481 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
482 return m_ErrorCode;
483 }
484
485 int i = 0;
486 for (auto& hm : m_BeforeTrackVertexError) {
487 for (int j = 0; j < 3; j++) for (int k = 0; k < KFitConst::kNumber7; k++) {
489 }
490 i++;
491 }
492 }
493 }
494
496}
497
498
501 Hep3Vector h3v;
502 int index = 0;
503 for (auto& pdata : m_Tracks)
504 {
505 // tracks
506 // momentum
507 h3v.setX(m_al_1[index * KFitConst::kNumber7 + 0][0]);
508 h3v.setY(m_al_1[index * KFitConst::kNumber7 + 1][0]);
509 h3v.setZ(m_al_1[index * KFitConst::kNumber7 + 2][0]);
510 if (m_IsFixMass[index])
511 pdata.setMomentum(HepLorentzVector(h3v, sqrt(h3v.mag2() + pdata.getMass()*pdata.getMass())), KFitConst::kAfterFit);
512 else
513 pdata.setMomentum(HepLorentzVector(h3v, m_al_1[index * KFitConst::kNumber7 + 3][0]), KFitConst::kAfterFit);
514 // position
515 pdata.setPosition(HepPoint3D(
516 m_al_1[index * KFitConst::kNumber7 + 4][0],
517 m_al_1[index * KFitConst::kNumber7 + 5][0],
519 // error of the tracks
520 pdata.setError(this->makeError3(pdata.getMomentum(),
521 m_V_al_1.sub(
522 index * KFitConst::kNumber7 + 1,
523 (index + 1)*KFitConst::kNumber7,
524 index * KFitConst::kNumber7 + 1,
525 (index + 1)*KFitConst::kNumber7), m_IsFixMass[index]),
527 if (m_ErrorCode != KFitError::kNoError) break;
528 index++;
529 }
530
532 {
533 //TODO: Not Implemented
535 // vertex
539 // error of the vertex
540 for (int i = 0; i < 3; i++) for (int j = i; j < 3; j++) {
542 }
543 // error between vertex and tracks
544 for (int i = 0; i < m_TrackCount; i++) {
545 HepMatrix hm(3, KFitConst::kNumber7, 0);
546 for (int j = 0; j < 3; j++) for (int k = 0; k < KFitConst::kNumber7; k++) {
548 }
549 if (m_IsFixMass[i])
550 m_AfterTrackVertexError.push_back(this->makeError4(m_Tracks[i].getMomentum(), hm));
551 else
552 m_AfterTrackVertexError.push_back(hm);
553 }
554 } else
555 {
556 // not fit
558 }
559
561}
562
563
567 {
568
569 HepMatrix al_1_prime(m_al_1);
570 HepMatrix Sum_al_1(4, 1, 0);
571 HepMatrix Sum_child_al_1(4 * m_ConstraintMassCount, 1, 0);
572 std::vector<double> energy(m_TrackCount);
573 double a;
574
575 for (int i = 0; i < m_TrackCount; i++) {
576 a = m_property[i][2];
577 if (!m_FlagAtDecayPoint) a = 0.;
578 al_1_prime[i * KFitConst::kNumber7 + 0][0] -= a * (m_BeforeVertex.y() - al_1_prime[i * KFitConst::kNumber7 + 5][0]);
579 al_1_prime[i * KFitConst::kNumber7 + 1][0] += a * (m_BeforeVertex.x() - al_1_prime[i * KFitConst::kNumber7 + 4][0]);
580 energy[i] = sqrt(al_1_prime[i * KFitConst::kNumber7 + 0][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] +
581 al_1_prime[i * KFitConst::kNumber7 + 1][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] +
582 al_1_prime[i * KFitConst::kNumber7 + 2][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] +
583 m_property[i][1] * m_property[i][1]);
584 if (m_IsFixMass[i])
585 Sum_al_1[3][0] += energy[i];
586 else
587 Sum_al_1[3][0] += al_1_prime[i * KFitConst::kNumber7 + 3][0];
588 }
589
590 for (int i = 0; i < m_TrackCount; i++) {
591 for (int j = 0; j < 3; j++) Sum_al_1[j][0] += al_1_prime[i * KFitConst::kNumber7 + j][0];
592 }
593
594 for (int i = 0; i < m_ConstraintMassCount; i++) {
595 for (int k = m_ConstraintMassChildLists[i].first; k <= m_ConstraintMassChildLists[i].second; k++) {
596 if (m_IsFixMass[k])
597 Sum_child_al_1[i * 4 + 3][0] += energy[k];
598 else
599 Sum_child_al_1[i * 4 + 3][0] += al_1_prime[k * KFitConst::kNumber7 + 3][0];
600 for (int j = 0; j < 3; j++) Sum_child_al_1[i * 4 + j][0] += al_1_prime[k * KFitConst::kNumber7 + j][0];
601 }
602 }
603
604 m_d[0][0] = Sum_al_1[0][0] - m_FourMomentum.Px();
605 m_d[1][0] = Sum_al_1[1][0] - m_FourMomentum.Py();
606 m_d[2][0] = Sum_al_1[2][0] - m_FourMomentum.Pz();
607 m_d[3][0] = Sum_al_1[3][0] - m_FourMomentum.E();
608
609 for (int i = 0; i < m_ConstraintMassCount; i++) {
610 m_d[4 + i][0] =
611 + Sum_child_al_1[i * 4 + 3][0] * Sum_child_al_1[i * 4 + 3][0] - Sum_child_al_1[i * 4 + 0][0] * Sum_child_al_1[i * 4 + 0][0]
612 - Sum_child_al_1[i * 4 + 1][0] * Sum_child_al_1[i * 4 + 1][0] - Sum_child_al_1[i * 4 + 2][0] * Sum_child_al_1[i * 4 + 2][0]
614 }
615
616 for (int i = 0; i < m_TrackCount; i++) {
617 if (energy[i] == 0) {
619 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
620 break;
621 }
622
623 a = m_property[i][2];
624 if (!m_FlagAtDecayPoint) a = 0.;
625
626 // four-momentum conservation constraint
627 for (int l = 0; l < 4; l++) {
628 for (int n = 0; n < 6; n++) {
629 m_D[l][i * KFitConst::kNumber7 + n] = 0;
630 }
631 }
632 if (m_IsFixMass[i]) {
633 double invE = 1. / energy[i];
634 m_D[0][i * KFitConst::kNumber7 + 0] = 1;
635 m_D[0][i * KFitConst::kNumber7 + 5] = -a;
636 m_D[1][i * KFitConst::kNumber7 + 1] = 1;
637 m_D[1][i * KFitConst::kNumber7 + 4] = a;
638 m_D[2][i * KFitConst::kNumber7 + 2] = 1;
639 m_D[3][i * KFitConst::kNumber7 + 0] = al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE;
640 m_D[3][i * KFitConst::kNumber7 + 1] = al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE;
641 m_D[3][i * KFitConst::kNumber7 + 2] = al_1_prime[i * KFitConst::kNumber7 + 2][0] * invE;
642 m_D[3][i * KFitConst::kNumber7 + 4] = -al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE * a;
643 m_D[3][i * KFitConst::kNumber7 + 5] = al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE * a;
644 } else {
645 m_D[0][i * KFitConst::kNumber7 + 0] = 1;
646 m_D[1][i * KFitConst::kNumber7 + 1] = 1;
647 m_D[2][i * KFitConst::kNumber7 + 2] = 1;
648 m_D[3][i * KFitConst::kNumber7 + 3] = 1;
649 }
650
651 // invariant mass constraint
652 for (int l = 0; l < m_ConstraintMassCount; l++) {
653 if (i >= m_ConstraintMassChildLists[l].first && i <= m_ConstraintMassChildLists[l].second) {
654 double invE = 1. / energy[i];
655 m_D[4 + l][i * KFitConst::kNumber7 + 0] = 2.*(Sum_child_al_1[l * 4 + 3][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE -
656 Sum_child_al_1[l * 4 + 0][0]);
657 m_D[4 + l][i * KFitConst::kNumber7 + 1] = 2.*(Sum_child_al_1[l * 4 + 3][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE -
658 Sum_child_al_1[l * 4 + 1][0]);
659 m_D[4 + l][i * KFitConst::kNumber7 + 2] = 2.*(Sum_child_al_1[l * 4 + 3][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] * invE -
660 Sum_child_al_1[l * 4 + 2][0]);
661 m_D[4 + l][i * KFitConst::kNumber7 + 3] = 0.;
662 m_D[4 + l][i * KFitConst::kNumber7 + 4] = -2.*(Sum_child_al_1[l * 4 + 3][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE -
663 Sum_child_al_1[l * 4 + 1][0]) * a;
664 m_D[4 + l][i * KFitConst::kNumber7 + 5] = 2.*(Sum_child_al_1[l * 4 + 3][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE -
665 Sum_child_al_1[l * 4 + 0][0]) * a;
666 m_D[4 + l][i * KFitConst::kNumber7 + 6] = 0.;
667 } else {
668 for (int n = 0; n < 6; n++) {
669 m_D[4 + l][i * KFitConst::kNumber7 + n] = 0;
670 }
671 }
672 }
673 }
674
675 } else
676 {
677 //TODO: Not Implemented
679
680 // m_FlagFitIncludingVertex == true
681 HepMatrix al_1_prime(m_al_1);
682 HepMatrix Sum_al_1(7, 1, 0);
683 double energy[KFitConst::kMaxTrackCount2];
684
685 for (int i = 0; i < m_TrackCount; i++) {
686 const double a = m_property[i][2];
687 al_1_prime[i * KFitConst::kNumber7 + 0][0] -= a * (al_1_prime[KFitConst::kNumber7 * m_TrackCount + 1][0] - al_1_prime[i *
688 KFitConst::kNumber7 + 5][0]);
689 al_1_prime[i * KFitConst::kNumber7 + 1][0] += a * (al_1_prime[KFitConst::kNumber7 * m_TrackCount + 0][0] - al_1_prime[i *
690 KFitConst::kNumber7 + 4][0]);
691 energy[i] = sqrt(al_1_prime[i * KFitConst::kNumber7 + 0][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] +
692 al_1_prime[i * KFitConst::kNumber7 + 1][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] +
693 al_1_prime[i * KFitConst::kNumber7 + 2][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] +
694 m_property[i][1] * m_property[i][1]);
695 Sum_al_1[6][0] = + a;
696 }
697
698 for (int i = 0; i < m_TrackCount; i++) {
699 if (energy[i] == 0) {
701 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
702 break;
703 }
704
705 if (m_IsFixMass[i]) {
706 double invE = 1. / energy[i];
707 Sum_al_1[3][0] += energy[i];
708 Sum_al_1[4][0] += al_1_prime[i * KFitConst::kNumber7 + 1][0] * m_property[i][2] * invE;
709 Sum_al_1[5][0] += al_1_prime[i * KFitConst::kNumber7 + 0][0] * m_property[i][2] * invE;
710 } else {
711 Sum_al_1[3][0] += al_1_prime[i * KFitConst::kNumber7 + 3][0];
712 }
713
714 for (int j = 0; j < 3; j++) Sum_al_1[j][0] += al_1_prime[i * KFitConst::kNumber7 + j][0];
715 }
716
717 m_d[0][0] = Sum_al_1[0][0] - m_FourMomentum.Px();
718 m_d[1][0] = Sum_al_1[1][0] - m_FourMomentum.Py();
719 m_d[2][0] = Sum_al_1[2][0] - m_FourMomentum.Pz();
720 m_d[3][0] = Sum_al_1[3][0] - m_FourMomentum.E();
721
722 for (int i = 0; i < m_TrackCount; i++) {
723 if (energy[i] == 0) {
725 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
726 break;
727 }
728
729 for (int l = 0; l < 4; l++) {
730 for (int n = 0; n < 6; n++) {
731 if (l == n) m_D[l][i * KFitConst::kNumber7 + n] = 1;
732 else m_D[l][i * KFitConst::kNumber7 + n] = 0;
733 }
734 }
735 }
736
737 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]);
738 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]);
739 m_D[0][KFitConst::kNumber7 * m_TrackCount + 2] = 0.;
740 }
741
743}
744
745
749
751}
752
754{
755 MakeMotherKFit kmm;
757 unsigned n = getTrackCount();
758 for (unsigned i = 0; i < n; ++i) {
760 getTrack(i).getCharge());
763 for (unsigned j = i + 1; j < n; ++j) {
765 }
766 }
767 kmm.setVertex(getVertex());
770 m_ErrorCode = kmm.doMake();
772 return m_ErrorCode;
773 double chi2 = getCHIsq();
774 int ndf = getNDF();
775 double prob = TMath::Prob(chi2, ndf);
776 //
777 bool haschi2 = mother->hasExtraInfo("chiSquared");
778 if (haschi2) {
779 mother->setExtraInfo("chiSquared", chi2);
780 mother->setExtraInfo("ndf", ndf);
781 } else {
782 mother->addExtraInfo("chiSquared", chi2);
783 mother->addExtraInfo("ndf", ndf);
784 }
785
786 mother->updateMomentum(
787 CLHEPToROOT::getLorentzVector(kmm.getMotherMomentum()),
788 CLHEPToROOT::getXYZVector(kmm.getMotherPosition()),
789 CLHEPToROOT::getTMatrixFSym(kmm.getMotherError()),
790 prob);
792 return m_ErrorCode;
793}
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
@ kUnimplemented
Unprepared.
Definition: KFitError.h:44
@ 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-four-momentum-constraint fit.
bool m_FlagTrackVertexError
Flag to indicate if the vertex error matrix of the child particle is preset.
enum KFitError::ECode setFourMomentum(const ROOT::Math::PxPyPzEVector &m)
Set an 4 Momentum for the mass-four-constraint fit.
bool m_FlagAtDecayPoint
Flag controlled by setFlagAtDecayPoint().
bool getFlagAtDecayPoint(void) const
Get a flag if to constraint at the decay point in the mass-four-momentum-constraint fit.
enum KFitError::ECode addMassConstraint(const double m, std::vector< unsigned > &childTrackId)
Set an invariant mass of daughter particle for the mass-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 mass-four-momentum-constraint fit.
enum KFitError::ECode setTrackZeroVertexError(void)
Indicate no vertex uncertainty in the child particle in the addTrack'ed order.
std::vector< int > m_IsFixMass
Array of flags whether the track property is fixed at the mass.
enum KFitError::ECode updateMother(Particle *mother)
Update mother particle.
enum KFitError::ECode unfixMass(void)
Tell the object to unfix the last added track property at the invariant mass.
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.
~MassFourCFitKFit(void)
Destruct the object.
enum KFitError::ECode doFit(void)
Perform a mass-four-momentum-constraint fit.
double m_InvariantMass
Invariant mass.
enum KFitError::ECode setInvariantMass(const double m)
Set an invariant mass for the mass-four-momentum-constraint fit.
std::vector< CLHEP::HepMatrix > m_BeforeTrackVertexError
array of vertex error matrices before the fit.
enum KFitError::ECode makeCoreMatrix(void) override
Build matrices using the kinematical constraint.
enum KFitError::ECode setVertexError(const CLHEP::HepSymMatrix &e)
Set an initial vertex error matrix for the mass-four-momentum-constraint fit.
enum KFitError::ECode setTrackVertexError(const CLHEP::HepMatrix &e)
Set a vertex error matrix of the child particle in the addTrack'ed order.
enum KFitError::ECode prepareCorrelation(void) override
Build a grand correlation matrix from input-track properties.
HepPoint3D m_AfterVertex
Vertex position after the fit.
enum KFitError::ECode setCorrelation(const CLHEP::HepMatrix &m) override
Set a correlation matrix.
std::vector< std::pair< int, int > > m_ConstraintMassChildLists
Daughter track id of constrained particle.
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.
MassFourCFitKFit()
Construct an object with no argument.
const HepPoint3D getVertex(const int flag=KFitConst::kAfterFit) const
Get a vertex position.
std::vector< double > m_ConstraintMassList
constrained mass
const CLHEP::HepMatrix getTrackVertexError(const int id, const int flag=KFitConst::kAfterFit) const
Get a vertex error matrix of the track.
CLHEP::HepSymMatrix m_AfterVertexError
Vertex error matrix after the fit.
enum KFitError::ECode prepareOutputMatrix(void) override
Build an output error matrix.
CLHEP::HepSymMatrix m_BeforeVertexError
Vertex error matrix before the fit.
HepPoint3D m_BeforeVertex
Vertex position before the fit.
enum KFitError::ECode fixMass(void)
Tell the object to fix the last added track property at the invariant mass.
ROOT::Math::PxPyPzEVector m_FourMomentum
Four Momentum.
double getInvariantMass(void) const
Get an invariant mass.
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 kMaxTrackCount2
Maximum track size (internal use)
Definition: KFitConst.h:42
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