Belle II Software development
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 <analysis/VertexFitting/KFit/MassFourCFitKFit.h>
10
11#include <analysis/dataobjects/Particle.h>
12#include <analysis/VertexFitting/KFit/MakeMotherKFit.h>
13#include <analysis/utility/CLHEPToROOT.h>
14#include <framework/gearbox/Const.h>
15
16#include <cstdio>
17
18#include <TMath.h>
19
20using namespace std;
21using namespace Belle2;
22using namespace Belle2::analysis;
23using namespace CLHEP;
24using namespace ROOT::Math;
25
39
40
42
44MassFourCFitKFit::addMassConstraint(const double m, std::vector<unsigned>& childTrackId) {
45 if ((childTrackId.back() - childTrackId.front()) != childTrackId.size() - 1)
46 {
48 }
49 m_ConstraintMassList.push_back(m);
50 m_ConstraintMassChildLists.push_back(std::make_pair(childTrackId.front(), childTrackId.back()));
52}
53
55MassFourCFitKFit::setVertex(const HepPoint3D& v) {
57
59}
60
61
63MassFourCFitKFit::setVertexError(const HepSymMatrix& e) {
64 if (e.num_row() != 3)
65 {
67 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
68 return m_ErrorCode;
69 }
70
73
75}
76
77
84
85
87MassFourCFitKFit::setFourMomentum(const PxPyPzEVector& m) {
89
91}
92
93
100
101
104 m_IsFixMass.push_back(true);
105
107}
108
109
112 m_IsFixMass.push_back(false);
113
115}
116
117
120 if (e.num_row() != 3 || e.num_col() != KFitConst::kNumber7)
121 {
123 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
124 return m_ErrorCode;
125 }
126
127 m_BeforeTrackVertexError.push_back(e);
130
132}
133
134
137 HepMatrix zero(3, KFitConst::kNumber7, 0);
138
139 return this->setTrackVertexError(zero);
140}
141
142
145 return KFitBase::setCorrelation(m);
146}
147
148
153
154
155const HepPoint3D
156MassFourCFitKFit::getVertex(const int flag) const
157{
158 if (flag == KFitConst::kAfterFit && !isFitted()) return HepPoint3D();
159
160 switch (flag) {
162 return m_BeforeVertex;
163
165 return m_AfterVertex;
166
167 default:
168 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
169 return HepPoint3D();
170 }
171}
172
173
174const HepSymMatrix
176{
177 if (flag == KFitConst::kAfterFit && !isFitted()) return HepSymMatrix(3, 0);
178
179 if (flag == KFitConst::kBeforeFit)
180 return m_BeforeVertexError;
182 return m_AfterVertexError;
183 else {
184 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
185 return HepSymMatrix(3, 0);
186 }
187}
188
189
190double
195
196
197bool
202
203
204bool
209
210
211double
213{
214 return m_CHIsq;
215}
216
217
218const HepMatrix
219MassFourCFitKFit::getTrackVertexError(const int id, const int flag) const
220{
221 if (flag == KFitConst::kAfterFit && !isFitted()) return HepMatrix(3, KFitConst::kNumber7, 0);
222 if (!isTrackIDInRange(id)) return HepMatrix(3, KFitConst::kNumber7, 0);
223
224 if (flag == KFitConst::kBeforeFit)
225 return m_BeforeTrackVertexError[id];
227 return m_AfterTrackVertexError[id];
228 else {
229 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
230 return HepMatrix(3, KFitConst::kNumber7, 0);
231 }
232}
233
234
235double
237{
238 if (!isFitted()) return -1;
239 if (!isTrackIDInRange(id)) return -1;
240
241 if (m_IsFixMass[id]) {
242
243 HepMatrix da(m_Tracks[id].getFitParameter(KFitConst::kBeforeFit) - m_Tracks[id].getFitParameter(KFitConst::kAfterFit));
244 int err_inverse = 0;
245 const double chisq = (da.T() * (m_Tracks[id].getFitError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
246
247 if (err_inverse) {
248 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
249 return -1;
250 }
251
252 return chisq;
253
254 } else {
255
256 HepMatrix da(m_Tracks[id].getMomPos(KFitConst::kBeforeFit) - m_Tracks[id].getMomPos(KFitConst::kAfterFit));
257 int err_inverse = 0;
258 const double chisq = (da.T() * (m_Tracks[id].getError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
259
260 if (err_inverse) {
261 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
262 return -1;
263 }
264
265 return chisq;
266
267 }
268}
269
270
271const HepMatrix
272MassFourCFitKFit::getCorrelation(const int id1, const int id2, const int flag) const
273{
274 if (flag == KFitConst::kAfterFit && !isFitted()) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
275 if (!isTrackIDInRange(id1)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
276 if (!isTrackIDInRange(id2)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
277
278 switch (flag) {
280 return KFitBase::getCorrelation(id1, id2, flag);
281
283 return makeError3(
284 this->getTrackMomentum(id1),
285 this->getTrackMomentum(id2),
286 m_V_al_1.sub(KFitConst::kNumber7 * id1 + 1, KFitConst::kNumber7 * (id1 + 1), KFitConst::kNumber7 * id2 + 1,
287 KFitConst::kNumber7 * (id2 + 1)),
288 m_IsFixMass[id1],
289 m_IsFixMass[id2]);
290
291 default:
292 KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
293 return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
294 }
295}
296
297
302
303
307 m_d = HepMatrix(4 + m_ConstraintMassCount, 1, 0);
308 m_V_D = HepMatrix(4 + m_ConstraintMassCount, 4 + m_ConstraintMassCount, 0);
309 m_lam = HepMatrix(4 + m_ConstraintMassCount, 1, 0);
310
312 {
314 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
315 return m_ErrorCode;
316 }
317
318 if (m_IsFixMass.size() == 0)
319 {
320 // If no fix_mass flag at all,
321 // all tracks are considered to be fixed at mass.
322 for (int i = 0; i < m_TrackCount; i++) this->fixMass();
323 } else if (m_IsFixMass.size() != (unsigned int)m_TrackCount)
324 {
326 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
327 return m_ErrorCode;
328 }
329
331 {
332 int index = 0;
333 m_al_0 = HepMatrix(KFitConst::kNumber7 * m_TrackCount, 1, 0);
334 m_property = HepMatrix(m_TrackCount, 3, 0);
335 m_V_al_0 = HepSymMatrix(KFitConst::kNumber7 * m_TrackCount, 0);
336
337 for (auto& track : m_Tracks) {
338 // momentum x,y,z and position x,y,z
339 m_al_0[index * KFitConst::kNumber7 + 0][0] = track.getMomentum(KFitConst::kBeforeFit).x();
340 m_al_0[index * KFitConst::kNumber7 + 1][0] = track.getMomentum(KFitConst::kBeforeFit).y();
341 m_al_0[index * KFitConst::kNumber7 + 2][0] = track.getMomentum(KFitConst::kBeforeFit).z();
342 m_al_0[index * KFitConst::kNumber7 + 3][0] = track.getMomentum(KFitConst::kBeforeFit).t();
343 m_al_0[index * KFitConst::kNumber7 + 4][0] = track.getPosition(KFitConst::kBeforeFit).x();
344 m_al_0[index * KFitConst::kNumber7 + 5][0] = track.getPosition(KFitConst::kBeforeFit).y();
345 m_al_0[index * KFitConst::kNumber7 + 6][0] = track.getPosition(KFitConst::kBeforeFit).z();
346 // these error
347 m_V_al_0.sub(index * KFitConst::kNumber7 + 1, track.getError(KFitConst::kBeforeFit));
348 // charge, mass, a
349 m_property[index][0] = track.getCharge();
350 m_property[index][1] = track.getMass();
351 const double c = Const::speedOfLight * 1e-4;
352 m_property[index][2] = -c * m_MagneticField * track.getCharge();
353 index++;
354 }
355
356 // error between track and track
357 if (m_FlagCorrelation) {
358 this->prepareCorrelation();
360 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
361 return m_ErrorCode;
362 }
363 }
364
365 // set member matrix
366 m_al_1 = m_al_0;
367
368 // define size of matrix
371
372 } else
373 {
374 //TODO: Not Implemented
376 // m_FlagFitIncludingVertex == true
377 int index = 0;
378 m_al_0 = HepMatrix(KFitConst::kNumber7 * m_TrackCount + 3, 1, 0);
379 m_property = HepMatrix(m_TrackCount, 3, 0);
380 m_V_al_0 = HepSymMatrix(KFitConst::kNumber7 * m_TrackCount + 3, 0);
381
382 for (auto& track : m_Tracks) {
383 // momentum x,y,z and position x,y,z
384 m_al_0[index * KFitConst::kNumber7 + 0][0] = track.getMomentum(KFitConst::kBeforeFit).x();
385 m_al_0[index * KFitConst::kNumber7 + 1][0] = track.getMomentum(KFitConst::kBeforeFit).y();
386 m_al_0[index * KFitConst::kNumber7 + 2][0] = track.getMomentum(KFitConst::kBeforeFit).z();
387 m_al_0[index * KFitConst::kNumber7 + 3][0] = track.getMomentum(KFitConst::kBeforeFit).t();
388 m_al_0[index * KFitConst::kNumber7 + 4][0] = track.getPosition(KFitConst::kBeforeFit).x();
389 m_al_0[index * KFitConst::kNumber7 + 5][0] = track.getPosition(KFitConst::kBeforeFit).y();
390 m_al_0[index * KFitConst::kNumber7 + 6][0] = track.getPosition(KFitConst::kBeforeFit).z();
391 // these error
392 m_V_al_0.sub(index * KFitConst::kNumber7 + 1, track.getError(KFitConst::kBeforeFit));
393 // charge, mass, a
394 m_property[index][0] = track.getCharge();
395 m_property[index][1] = track.getMass();
396 const double c = Const::speedOfLight * 1e-4;
397 m_property[index][2] = -c * m_MagneticField * track.getCharge();
398 index++;
399 }
400
401 // vertex
406
407 // error between track and track
408 if (m_FlagCorrelation) {
409 this->prepareCorrelation();
411 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
412 return m_ErrorCode;
413 }
414 }
415
416 // set member matrix
417 m_al_1 = m_al_0;
418
419 // define size of matrix
421 m_D = m_V_al_1.sub(1, 4, 1, KFitConst::kNumber7 * m_TrackCount + 3);
422 }
423
425}
426
427
430 char buf[1024];
431 sprintf(buf, "%s:%s(): internal error; this function should never be called", __FILE__, __func__);
432 B2FATAL(buf);
433
434 /* NEVER REACHEd HERE */
436}
437
438
441 if (m_BeforeCorrelation.size() != static_cast<unsigned int>(m_TrackCount * (m_TrackCount - 1) / 2))
442 {
444 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
445 return m_ErrorCode;
446 }
447
448 int row = 0, col = 0;
449
450 for (auto& hm : m_BeforeCorrelation)
451 {
452 // counter
453 row++;
454 if (row == m_TrackCount) {
455 col++;
456 row = col + 1;
457 }
458
459 int ii = 0, jj = 0;
460 for (int i = KFitConst::kNumber7 * row; i < KFitConst::kNumber7 * (row + 1); i++) {
461 for (int j = KFitConst::kNumber7 * col; j < KFitConst::kNumber7 * (col + 1); j++) {
462 m_V_al_0[i][j] = hm[ii][jj];
463 jj++;
464 }
465 jj = 0;
466 ii++;
467 }
468 }
469
471 {
472 //TODO: Not Implemented
474
475 // ...error of vertex
477
478 // ...error matrix between vertex and tracks
480 if (m_BeforeTrackVertexError.size() != (unsigned int)m_TrackCount) {
482 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
483 return m_ErrorCode;
484 }
485
486 int i = 0;
487 for (auto& hm : m_BeforeTrackVertexError) {
488 for (int j = 0; j < 3; j++) for (int k = 0; k < KFitConst::kNumber7; k++) {
490 }
491 i++;
492 }
493 }
494 }
495
497}
498
499
502 Hep3Vector h3v;
503 int index = 0;
504 for (auto& pdata : m_Tracks)
505 {
506 // tracks
507 // momentum
508 h3v.setX(m_al_1[index * KFitConst::kNumber7 + 0][0]);
509 h3v.setY(m_al_1[index * KFitConst::kNumber7 + 1][0]);
510 h3v.setZ(m_al_1[index * KFitConst::kNumber7 + 2][0]);
511 if (m_IsFixMass[index])
512 pdata.setMomentum(HepLorentzVector(h3v, sqrt(h3v.mag2() + pdata.getMass()*pdata.getMass())), KFitConst::kAfterFit);
513 else
514 pdata.setMomentum(HepLorentzVector(h3v, m_al_1[index * KFitConst::kNumber7 + 3][0]), KFitConst::kAfterFit);
515 // position
516 pdata.setPosition(HepPoint3D(
517 m_al_1[index * KFitConst::kNumber7 + 4][0],
518 m_al_1[index * KFitConst::kNumber7 + 5][0],
520 // error of the tracks
521 pdata.setError(this->makeError3(pdata.getMomentum(),
522 m_V_al_1.sub(
523 index * KFitConst::kNumber7 + 1,
524 (index + 1)*KFitConst::kNumber7,
525 index * KFitConst::kNumber7 + 1,
526 (index + 1)*KFitConst::kNumber7), m_IsFixMass[index]),
528 if (m_ErrorCode != KFitError::kNoError) break;
529 index++;
530 }
531
533 {
534 //TODO: Not Implemented
536 // vertex
540 // error of the vertex
541 for (int i = 0; i < 3; i++) for (int j = i; j < 3; j++) {
543 }
544 // error between vertex and tracks
545 for (int i = 0; i < m_TrackCount; i++) {
546 HepMatrix hm(3, KFitConst::kNumber7, 0);
547 for (int j = 0; j < 3; j++) for (int k = 0; k < KFitConst::kNumber7; k++) {
549 }
550 if (m_IsFixMass[i])
551 m_AfterTrackVertexError.push_back(this->makeError4(m_Tracks[i].getMomentum(), hm));
552 else
553 m_AfterTrackVertexError.push_back(hm);
554 }
555 } else
556 {
557 // not fit
559 }
560
562}
563
564
568 {
569
570 HepMatrix al_1_prime(m_al_1);
571 HepMatrix Sum_al_1(4, 1, 0);
572 HepMatrix Sum_child_al_1(4 * m_ConstraintMassCount, 1, 0);
573 std::vector<double> energy(m_TrackCount);
574 double a;
575
576 for (int i = 0; i < m_TrackCount; i++) {
577 a = m_property[i][2];
578 if (!m_FlagAtDecayPoint) a = 0.;
579 al_1_prime[i * KFitConst::kNumber7 + 0][0] -= a * (m_BeforeVertex.y() - al_1_prime[i * KFitConst::kNumber7 + 5][0]);
580 al_1_prime[i * KFitConst::kNumber7 + 1][0] += a * (m_BeforeVertex.x() - al_1_prime[i * KFitConst::kNumber7 + 4][0]);
581 energy[i] = sqrt(al_1_prime[i * KFitConst::kNumber7 + 0][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] +
582 al_1_prime[i * KFitConst::kNumber7 + 1][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] +
583 al_1_prime[i * KFitConst::kNumber7 + 2][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] +
584 m_property[i][1] * m_property[i][1]);
585 if (m_IsFixMass[i])
586 Sum_al_1[3][0] += energy[i];
587 else
588 Sum_al_1[3][0] += al_1_prime[i * KFitConst::kNumber7 + 3][0];
589 }
590
591 for (int i = 0; i < m_TrackCount; i++) {
592 for (int j = 0; j < 3; j++) Sum_al_1[j][0] += al_1_prime[i * KFitConst::kNumber7 + j][0];
593 }
594
595 for (int i = 0; i < m_ConstraintMassCount; i++) {
596 for (int k = m_ConstraintMassChildLists[i].first; k <= m_ConstraintMassChildLists[i].second; k++) {
597 if (m_IsFixMass[k])
598 Sum_child_al_1[i * 4 + 3][0] += energy[k];
599 else
600 Sum_child_al_1[i * 4 + 3][0] += al_1_prime[k * KFitConst::kNumber7 + 3][0];
601 for (int j = 0; j < 3; j++) Sum_child_al_1[i * 4 + j][0] += al_1_prime[k * KFitConst::kNumber7 + j][0];
602 }
603 }
604
605 m_d[0][0] = Sum_al_1[0][0] - m_FourMomentum.Px();
606 m_d[1][0] = Sum_al_1[1][0] - m_FourMomentum.Py();
607 m_d[2][0] = Sum_al_1[2][0] - m_FourMomentum.Pz();
608 m_d[3][0] = Sum_al_1[3][0] - m_FourMomentum.E();
609
610 for (int i = 0; i < m_ConstraintMassCount; i++) {
611 m_d[4 + i][0] =
612 + 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]
613 - 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]
615 }
616
617 for (int i = 0; i < m_TrackCount; i++) {
618 if (energy[i] == 0) {
620 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
621 break;
622 }
623
624 a = m_property[i][2];
625 if (!m_FlagAtDecayPoint) a = 0.;
626
627 // four-momentum conservation constraint
628 for (int l = 0; l < 4; l++) {
629 for (int n = 0; n < 6; n++) {
630 m_D[l][i * KFitConst::kNumber7 + n] = 0;
631 }
632 }
633 if (m_IsFixMass[i]) {
634 double invE = 1. / energy[i];
635 m_D[0][i * KFitConst::kNumber7 + 0] = 1;
636 m_D[0][i * KFitConst::kNumber7 + 5] = -a;
637 m_D[1][i * KFitConst::kNumber7 + 1] = 1;
638 m_D[1][i * KFitConst::kNumber7 + 4] = a;
639 m_D[2][i * KFitConst::kNumber7 + 2] = 1;
640 m_D[3][i * KFitConst::kNumber7 + 0] = al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE;
641 m_D[3][i * KFitConst::kNumber7 + 1] = al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE;
642 m_D[3][i * KFitConst::kNumber7 + 2] = al_1_prime[i * KFitConst::kNumber7 + 2][0] * invE;
643 m_D[3][i * KFitConst::kNumber7 + 4] = -al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE * a;
644 m_D[3][i * KFitConst::kNumber7 + 5] = al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE * a;
645 } else {
646 m_D[0][i * KFitConst::kNumber7 + 0] = 1;
647 m_D[1][i * KFitConst::kNumber7 + 1] = 1;
648 m_D[2][i * KFitConst::kNumber7 + 2] = 1;
649 m_D[3][i * KFitConst::kNumber7 + 3] = 1;
650 }
651
652 // invariant mass constraint
653 for (int l = 0; l < m_ConstraintMassCount; l++) {
654 if (i >= m_ConstraintMassChildLists[l].first && i <= m_ConstraintMassChildLists[l].second) {
655 double invE = 1. / energy[i];
656 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 -
657 Sum_child_al_1[l * 4 + 0][0]);
658 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 -
659 Sum_child_al_1[l * 4 + 1][0]);
660 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 -
661 Sum_child_al_1[l * 4 + 2][0]);
662 m_D[4 + l][i * KFitConst::kNumber7 + 3] = 0.;
663 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 -
664 Sum_child_al_1[l * 4 + 1][0]) * a;
665 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 -
666 Sum_child_al_1[l * 4 + 0][0]) * a;
667 m_D[4 + l][i * KFitConst::kNumber7 + 6] = 0.;
668 } else {
669 for (int n = 0; n < 6; n++) {
670 m_D[4 + l][i * KFitConst::kNumber7 + n] = 0;
671 }
672 }
673 }
674 }
675
676 } else
677 {
678 //TODO: Not Implemented
680
681 // m_FlagFitIncludingVertex == true
682 HepMatrix al_1_prime(m_al_1);
683 HepMatrix Sum_al_1(7, 1, 0);
684 double energy[KFitConst::kMaxTrackCount2];
685
686 for (int i = 0; i < m_TrackCount; i++) {
687 const double a = m_property[i][2];
688 al_1_prime[i * KFitConst::kNumber7 + 0][0] -= a * (al_1_prime[KFitConst::kNumber7 * m_TrackCount + 1][0] - al_1_prime[i *
689 KFitConst::kNumber7 + 5][0]);
690 al_1_prime[i * KFitConst::kNumber7 + 1][0] += a * (al_1_prime[KFitConst::kNumber7 * m_TrackCount + 0][0] - al_1_prime[i *
691 KFitConst::kNumber7 + 4][0]);
692 energy[i] = sqrt(al_1_prime[i * KFitConst::kNumber7 + 0][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] +
693 al_1_prime[i * KFitConst::kNumber7 + 1][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] +
694 al_1_prime[i * KFitConst::kNumber7 + 2][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] +
695 m_property[i][1] * m_property[i][1]);
696 Sum_al_1[6][0] = + a;
697 }
698
699 for (int i = 0; i < m_TrackCount; i++) {
700 if (energy[i] == 0) {
702 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
703 break;
704 }
705
706 if (m_IsFixMass[i]) {
707 double invE = 1. / energy[i];
708 Sum_al_1[3][0] += energy[i];
709 Sum_al_1[4][0] += al_1_prime[i * KFitConst::kNumber7 + 1][0] * m_property[i][2] * invE;
710 Sum_al_1[5][0] += al_1_prime[i * KFitConst::kNumber7 + 0][0] * m_property[i][2] * invE;
711 } else {
712 Sum_al_1[3][0] += al_1_prime[i * KFitConst::kNumber7 + 3][0];
713 }
714
715 for (int j = 0; j < 3; j++) Sum_al_1[j][0] += al_1_prime[i * KFitConst::kNumber7 + j][0];
716 }
717
718 m_d[0][0] = Sum_al_1[0][0] - m_FourMomentum.Px();
719 m_d[1][0] = Sum_al_1[1][0] - m_FourMomentum.Py();
720 m_d[2][0] = Sum_al_1[2][0] - m_FourMomentum.Pz();
721 m_d[3][0] = Sum_al_1[3][0] - m_FourMomentum.E();
722
723 for (int i = 0; i < m_TrackCount; i++) {
724 if (energy[i] == 0) {
726 KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
727 break;
728 }
729
730 for (int l = 0; l < 4; l++) {
731 for (int n = 0; n < 6; n++) {
732 if (l == n) m_D[l][i * KFitConst::kNumber7 + n] = 1;
733 else m_D[l][i * KFitConst::kNumber7 + n] = 0;
734 }
735 }
736 }
737
738 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]);
739 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]);
740 m_D[0][KFitConst::kNumber7 * m_TrackCount + 2] = 0.;
741 }
742
744}
745
746
753
755{
756 MakeMotherKFit kmm;
758 unsigned n = getTrackCount();
759 for (unsigned i = 0; i < n; ++i) {
761 getTrack(i).getCharge());
764 for (unsigned j = i + 1; j < n; ++j) {
766 }
767 }
768 kmm.setVertex(getVertex());
771 m_ErrorCode = kmm.doMake();
773 return m_ErrorCode;
774 double chi2 = getCHIsq();
775 int ndf = getNDF();
776 double prob = TMath::Prob(chi2, ndf);
777 //
778 bool haschi2 = mother->hasExtraInfo("chiSquared");
779 if (haschi2) {
780 mother->setExtraInfo("chiSquared", chi2);
781 mother->setExtraInfo("ndf", ndf);
782 } else {
783 mother->addExtraInfo("chiSquared", chi2);
784 mother->addExtraInfo("ndf", ndf);
785 }
786
787 mother->updateMomentum(
788 CLHEPToROOT::getLorentzVector(kmm.getMotherMomentum()),
789 CLHEPToROOT::getXYZVector(kmm.getMotherPosition()),
790 CLHEPToROOT::getTMatrixFSym(kmm.getMotherError()),
791 prob);
793 return m_ErrorCode;
794}
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
@ kUnimplemented
Unprepared.
Definition KFitError.h:43
@ 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-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.
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.
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 kMaxTrackCount2
Maximum track size (internal use)
Definition KFitConst.h:40
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