Belle II Software  release-05-01-25
FourCFitKFit.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2017 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributor: Yu Hu *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 
12 #include <cstdio>
13 
14 #include <TMatrixFSym.h>
15 
16 #include <analysis/VertexFitting/KFit/FourCFitKFit.h>
17 #include <analysis/VertexFitting/KFit/MakeMotherKFit.h>
18 #include <analysis/utility/CLHEPToROOT.h>
19 #include <TLorentzVector.h>
20 
21 
22 using namespace std;
23 using namespace Belle2;
24 using namespace Belle2::analysis;
25 using namespace CLHEP;
26 
27 FourCFitKFit::FourCFitKFit()
28 {
29  m_FlagFitted = false;
30  m_FlagTrackVertexError = false;
31  m_FlagFitIncludingVertex = false;
32  m_FlagAtDecayPoint = true;
33  m_NecessaryTrackCount = 2;
34  m_d = HepMatrix(4, 1, 0);
35  m_V_D = HepMatrix(4, 4, 0);
36  m_lam = HepMatrix(4, 1, 0);
37  m_AfterVertexError = HepSymMatrix(3, 0);
38  m_InvariantMass = -1.0;
39  m_FourMomentum = TLorentzVector();
40 }
41 
42 
43 FourCFitKFit::~FourCFitKFit() = default;
44 
45 
46 enum KFitError::ECode
47 FourCFitKFit::setVertex(const HepPoint3D& v) {
48  m_BeforeVertex = v;
49 
50  return m_ErrorCode = KFitError::kNoError;
51 }
52 
53 
54 enum KFitError::ECode
55 FourCFitKFit::setVertexError(const HepSymMatrix& e) {
56  if (e.num_row() != 3)
57  {
58  m_ErrorCode = KFitError::kBadMatrixSize;
59  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
60  return m_ErrorCode;
61  }
62 
63  m_BeforeVertexError = e;
64  m_FlagFitIncludingVertex = true;
65 
66  return m_ErrorCode = KFitError::kNoError;
67 }
68 
69 
70 enum KFitError::ECode
71 FourCFitKFit::setInvariantMass(const double m) {
72  m_InvariantMass = m;
73 
74  return m_ErrorCode = KFitError::kNoError;
75 }
76 
77 
78 enum KFitError::ECode
79 FourCFitKFit::setFourMomentum(const TLorentzVector& m) {
80  m_FourMomentum = m;
81 
82  return m_ErrorCode = KFitError::kNoError;
83 }
84 
85 
86 enum KFitError::ECode
87 FourCFitKFit::setFlagAtDecayPoint(const bool flag) {
88  m_FlagAtDecayPoint = flag;
89 
90  return m_ErrorCode = KFitError::kNoError;
91 }
92 
93 
94 enum KFitError::ECode
95 FourCFitKFit::fixMass() {
96  m_IsFixMass.push_back(true);
97 
98  return m_ErrorCode = KFitError::kNoError;
99 }
100 
101 
102 enum KFitError::ECode
103 FourCFitKFit::unfixMass() {
104  m_IsFixMass.push_back(false);
105 
106  return m_ErrorCode = KFitError::kNoError;
107 }
108 
109 
110 enum KFitError::ECode
111 FourCFitKFit::setTrackVertexError(const HepMatrix& e) {
112  if (e.num_row() != 3 || e.num_col() != KFitConst::kNumber7)
113  {
114  m_ErrorCode = KFitError::kBadMatrixSize;
115  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
116  return m_ErrorCode;
117  }
118 
119  m_BeforeTrackVertexError.push_back(e);
120  m_FlagTrackVertexError = true;
121  m_FlagFitIncludingVertex = true;
122 
123  return m_ErrorCode = KFitError::kNoError;
124 }
125 
126 
127 enum KFitError::ECode
128 FourCFitKFit::setTrackZeroVertexError() {
129  HepMatrix zero(3, KFitConst::kNumber7, 0);
130 
131  return this->setTrackVertexError(zero);
132 }
133 
134 
135 enum KFitError::ECode
136 FourCFitKFit::setCorrelation(const HepMatrix& m) {
137  return KFitBase::setCorrelation(m);
138 }
139 
140 
141 enum KFitError::ECode
142 FourCFitKFit::setZeroCorrelation() {
143  return KFitBase::setZeroCorrelation();
144 }
145 
146 
147 const HepPoint3D
148 FourCFitKFit::getVertex(const int flag) const
149 {
150  if (flag == KFitConst::kAfterFit && !isFitted()) return HepPoint3D();
151 
152  switch (flag) {
153  case KFitConst::kBeforeFit:
154  return m_BeforeVertex;
155 
156  case KFitConst::kAfterFit:
157  return m_AfterVertex;
158 
159  default:
160  KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
161  return HepPoint3D();
162  }
163 }
164 
165 
166 const HepSymMatrix
167 FourCFitKFit::getVertexError(const int flag) const
168 {
169  if (flag == KFitConst::kAfterFit && !isFitted()) return HepSymMatrix(3, 0);
170 
171  if (flag == KFitConst::kBeforeFit)
172  return m_BeforeVertexError;
173  else if (flag == KFitConst::kAfterFit && m_FlagFitIncludingVertex)
174  return m_AfterVertexError;
175  else {
176  KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
177  return HepSymMatrix(3, 0);
178  }
179 }
180 
181 
182 double
183 FourCFitKFit::getInvariantMass() const
184 {
185  return m_InvariantMass;
186 }
187 
188 
189 bool
190 FourCFitKFit::getFlagAtDecayPoint() const
191 {
192  return m_FlagAtDecayPoint;
193 }
194 
195 
196 bool
197 FourCFitKFit::getFlagFitWithVertex() const
198 {
199  return m_FlagFitIncludingVertex;
200 }
201 
202 
203 double
204 FourCFitKFit::getCHIsq() const
205 {
206  return m_CHIsq;
207 }
208 
209 
210 const HepMatrix
211 FourCFitKFit::getTrackVertexError(const int id, const int flag) const
212 {
213  if (flag == KFitConst::kAfterFit && !isFitted()) return HepMatrix(3, KFitConst::kNumber7, 0);
214  if (!isTrackIDInRange(id)) return HepMatrix(3, KFitConst::kNumber7, 0);
215 
216  if (flag == KFitConst::kBeforeFit)
217  return m_BeforeTrackVertexError[id];
218  else if (flag == KFitConst::kAfterFit && m_FlagFitIncludingVertex)
219  return m_AfterTrackVertexError[id];
220  else {
221  KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
222  return HepMatrix(3, KFitConst::kNumber7, 0);
223  }
224 }
225 
226 
227 double
228 FourCFitKFit::getTrackCHIsq(const int id) const
229 {
230  if (!isFitted()) return -1;
231  if (!isTrackIDInRange(id)) return -1;
232 
233  if (m_IsFixMass[id]) {
234 
235  HepMatrix da(m_Tracks[id].getFitParameter(KFitConst::kBeforeFit) - m_Tracks[id].getFitParameter(KFitConst::kAfterFit));
236  int err_inverse = 0;
237  const double chisq = (da.T() * (m_Tracks[id].getFitError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
238 
239  if (err_inverse) {
240  KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
241  return -1;
242  }
243 
244  return chisq;
245 
246  } else {
247 
248  HepMatrix da(m_Tracks[id].getMomPos(KFitConst::kBeforeFit) - m_Tracks[id].getMomPos(KFitConst::kAfterFit));
249  int err_inverse = 0;
250  const double chisq = (da.T() * (m_Tracks[id].getError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
251 
252  if (err_inverse) {
253  KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
254  return -1;
255  }
256 
257  return chisq;
258 
259  }
260 }
261 
262 
263 const HepMatrix
264 FourCFitKFit::getCorrelation(const int id1, const int id2, const int flag) const
265 {
266  if (flag == KFitConst::kAfterFit && !isFitted()) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
267  if (!isTrackIDInRange(id1)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
268  if (!isTrackIDInRange(id2)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
269 
270  switch (flag) {
271  case KFitConst::kBeforeFit:
272  return KFitBase::getCorrelation(id1, id2, flag);
273 
274  case KFitConst::kAfterFit:
275  return makeError3(
276  this->getTrackMomentum(id1),
277  this->getTrackMomentum(id2),
278  m_V_al_1.sub(KFitConst::kNumber7 * id1 + 1, KFitConst::kNumber7 * (id1 + 1), KFitConst::kNumber7 * id2 + 1,
279  KFitConst::kNumber7 * (id2 + 1)),
280  m_IsFixMass[id1],
281  m_IsFixMass[id2]);
282 
283  default:
284  KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
285  return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
286  }
287 }
288 
289 
290 enum KFitError::ECode
291 FourCFitKFit::doFit() {
292  return KFitBase::doFit1();
293 }
294 
295 
296 enum KFitError::ECode
297 FourCFitKFit::prepareInputMatrix() {
298  if (m_TrackCount > KFitConst::kMaxTrackCount)
299  {
300  m_ErrorCode = KFitError::kBadTrackSize;
301  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
302  return m_ErrorCode;
303  }
304 
305 
306  if (m_IsFixMass.size() == 0)
307  {
308  // If no fix_mass flag at all,
309  // all tracks are considered to be fixed at mass.
310  for (int i = 0; i < m_TrackCount; i++) this->fixMass();
311  } else if (m_IsFixMass.size() != (unsigned int)m_TrackCount)
312  {
313  m_ErrorCode = KFitError::kBadTrackSize;
314  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
315  return m_ErrorCode;
316  }
317 
318 
319  if (!m_FlagFitIncludingVertex)
320  {
321  int index = 0;
322  m_al_0 = HepMatrix(KFitConst::kNumber7 * m_TrackCount, 1, 0);
323  m_property = HepMatrix(m_TrackCount, 3, 0);
324  m_V_al_0 = HepSymMatrix(KFitConst::kNumber7 * m_TrackCount, 0);
325 
326  for (auto& track : m_Tracks) {
327  // momentum x,y,z and position x,y,z
328  m_al_0[index * KFitConst::kNumber7 + 0][0] = track.getMomentum(KFitConst::kBeforeFit).x();
329  m_al_0[index * KFitConst::kNumber7 + 1][0] = track.getMomentum(KFitConst::kBeforeFit).y();
330  m_al_0[index * KFitConst::kNumber7 + 2][0] = track.getMomentum(KFitConst::kBeforeFit).z();
331  m_al_0[index * KFitConst::kNumber7 + 3][0] = track.getMomentum(KFitConst::kBeforeFit).t();
332  m_al_0[index * KFitConst::kNumber7 + 4][0] = track.getPosition(KFitConst::kBeforeFit).x();
333  m_al_0[index * KFitConst::kNumber7 + 5][0] = track.getPosition(KFitConst::kBeforeFit).y();
334  m_al_0[index * KFitConst::kNumber7 + 6][0] = track.getPosition(KFitConst::kBeforeFit).z();
335  // these error
336  m_V_al_0.sub(index * KFitConst::kNumber7 + 1, track.getError(KFitConst::kBeforeFit));
337  // charge, mass, a
338  m_property[index][0] = track.getCharge();
339  m_property[index][1] = track.getMass();
340  const double c = KFitConst::kLightSpeed; // C++ bug?
341  // m_property[index][2] = -KFitConst::kLightSpeed * m_MagneticField * it->getCharge();
342  m_property[index][2] = -c * m_MagneticField * track.getCharge();
343  index++;
344  }
345 
346  // error between track and track
347  if (m_FlagCorrelation) {
348  this->prepareCorrelation();
349  if (m_ErrorCode != KFitError::kNoError) {
350  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
351  return m_ErrorCode;
352  }
353  }
354 
355  // set member matrix
356  m_al_1 = m_al_0;
357 
358  // define size of matrix
359  m_V_al_1 = HepMatrix(KFitConst::kNumber7 * m_TrackCount, KFitConst::kNumber7 * m_TrackCount, 0);
360  m_D = m_V_al_1.sub(1, 4, 1, KFitConst::kNumber7 * m_TrackCount);
361 
362  } else {
363  // m_FlagFitIncludingVertex == true
364  int index = 0;
365  m_al_0 = HepMatrix(KFitConst::kNumber7 * m_TrackCount + 3, 1, 0);
366  m_property = HepMatrix(m_TrackCount, 3, 0);
367  m_V_al_0 = HepSymMatrix(KFitConst::kNumber7 * m_TrackCount + 3, 0);
368 
369  for (auto& track : m_Tracks)
370  {
371  // momentum x,y,z and position x,y,z
372  m_al_0[index * KFitConst::kNumber7 + 0][0] = track.getMomentum(KFitConst::kBeforeFit).x();
373  m_al_0[index * KFitConst::kNumber7 + 1][0] = track.getMomentum(KFitConst::kBeforeFit).y();
374  m_al_0[index * KFitConst::kNumber7 + 2][0] = track.getMomentum(KFitConst::kBeforeFit).z();
375  m_al_0[index * KFitConst::kNumber7 + 3][0] = track.getMomentum(KFitConst::kBeforeFit).t();
376  m_al_0[index * KFitConst::kNumber7 + 4][0] = track.getPosition(KFitConst::kBeforeFit).x();
377  m_al_0[index * KFitConst::kNumber7 + 5][0] = track.getPosition(KFitConst::kBeforeFit).y();
378  m_al_0[index * KFitConst::kNumber7 + 6][0] = track.getPosition(KFitConst::kBeforeFit).z();
379  // these error
380  m_V_al_0.sub(index * KFitConst::kNumber7 + 1, track.getError(KFitConst::kBeforeFit));
381  // charge, mass, a
382  m_property[index][0] = track.getCharge();
383  m_property[index][1] = track.getMass();
384  const double c = KFitConst::kLightSpeed; // C++ bug?
385  // m_property[index][2] = -KFitConst::kLightSpeed * m_MagneticField * it->getCharge();
386  m_property[index][2] = -c * m_MagneticField * track.getCharge();
387  index++;
388  }
389 
390  // vertex
391  m_al_0[KFitConst::kNumber7 * m_TrackCount + 0][0] = m_BeforeVertex.x();
392  m_al_0[KFitConst::kNumber7 * m_TrackCount + 1][0] = m_BeforeVertex.y();
393  m_al_0[KFitConst::kNumber7 * m_TrackCount + 2][0] = m_BeforeVertex.z();
394  m_V_al_0.sub(KFitConst::kNumber7 * m_TrackCount + 1, m_BeforeVertexError);
395 
396  // error between track and track
397  if (m_FlagCorrelation)
398  {
399  this->prepareCorrelation();
400  if (m_ErrorCode != KFitError::kNoError) {
401  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
402  return m_ErrorCode;
403  }
404  }
405 
406  // set member matrix
407  m_al_1 = m_al_0;
408 
409  // define size of matrix
410  m_V_al_1 = HepMatrix(KFitConst::kNumber7 * m_TrackCount + 3, KFitConst::kNumber7 * m_TrackCount + 3, 0);
411  m_D = m_V_al_1.sub(1, 4, 1, KFitConst::kNumber7 * m_TrackCount + 3);
412  }
413 
414  return m_ErrorCode = KFitError::kNoError;
415 }
416 
417 
418 enum KFitError::ECode
419 FourCFitKFit::prepareInputSubMatrix() { // unused
420  char buf[1024];
421  sprintf(buf, "%s:%s(): internal error; this function should never be called", __FILE__, __func__);
422  B2FATAL(buf);
423 
424  /* NEVER REACHEd HERE */
425  return KFitError::kOutOfRange;
426 }
427 
428 
429 enum KFitError::ECode
430 FourCFitKFit::prepareCorrelation() {
431  if (m_BeforeCorrelation.size() != static_cast<unsigned int>(m_TrackCount * (m_TrackCount - 1) / 2))
432  {
433  m_ErrorCode = KFitError::kBadCorrelationSize;
434  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
435  return m_ErrorCode;
436  }
437 
438  int row = 0, col = 0;
439 
440  for (auto& hm : m_BeforeCorrelation)
441  {
442  // counter
443  row++;
444  if (row == m_TrackCount) {
445  col++;
446  row = col + 1;
447  }
448 
449  int ii = 0, jj = 0;
450  for (int i = KFitConst::kNumber7 * row; i < KFitConst::kNumber7 * (row + 1); i++) {
451  for (int j = KFitConst::kNumber7 * col; j < KFitConst::kNumber7 * (col + 1); j++) {
452  m_V_al_0[i][j] = hm[ii][jj];
453  jj++;
454  }
455  jj = 0;
456  ii++;
457  }
458  }
459 
460  if (m_FlagFitIncludingVertex)
461  {
462  // ...error of vertex
463  m_V_al_0.sub(KFitConst::kNumber7 * m_TrackCount + 1, m_BeforeVertexError);
464 
465  // ...error matrix between vertex and tracks
466  if (m_FlagTrackVertexError) {
467  if (m_BeforeTrackVertexError.size() != (unsigned int)m_TrackCount) {
468  m_ErrorCode = KFitError::kBadCorrelationSize;
469  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
470  return m_ErrorCode;
471  }
472 
473  int i = 0;
474  for (auto& hm : m_BeforeTrackVertexError) {
475  for (int j = 0; j < 3; j++) for (int k = 0; k < KFitConst::kNumber7; k++) {
476  m_V_al_0[j + KFitConst::kNumber7 * m_TrackCount][k + i * KFitConst::kNumber7] = hm[j][k];
477  }
478  i++;
479  }
480  }
481  }
482 
483  return m_ErrorCode = KFitError::kNoError;
484 }
485 
486 
487 enum KFitError::ECode
488 FourCFitKFit::prepareOutputMatrix() {
489  Hep3Vector h3v;
490  int index = 0;
491  for (auto& pdata : m_Tracks)
492  {
493  // tracks
494  // momentum
495  h3v.setX(m_al_1[index * KFitConst::kNumber7 + 0][0]);
496  h3v.setY(m_al_1[index * KFitConst::kNumber7 + 1][0]);
497  h3v.setZ(m_al_1[index * KFitConst::kNumber7 + 2][0]);
498  pdata.setMomentum(HepLorentzVector(h3v, m_al_1[index * KFitConst::kNumber7 + 3][0]), KFitConst::kAfterFit);
499  // position
500  pdata.setPosition(HepPoint3D(
501  m_al_1[index * KFitConst::kNumber7 + 4][0],
502  m_al_1[index * KFitConst::kNumber7 + 5][0],
503  m_al_1[index * KFitConst::kNumber7 + 6][0]), KFitConst::kAfterFit);
504  // error of the tracks
505  pdata.setError(this->makeError3(pdata.getMomentum(),
506  m_V_al_1.sub(
507  index * KFitConst::kNumber7 + 1,
508  (index + 1)*KFitConst::kNumber7,
509  index * KFitConst::kNumber7 + 1,
510  (index + 1)*KFitConst::kNumber7), m_IsFixMass[index]),
511  KFitConst::kAfterFit);
512  if (m_ErrorCode != KFitError::kNoError) break;
513  index++;
514  }
515 
516  if (m_FlagFitIncludingVertex)
517  {
518  // vertex
519  m_AfterVertex.setX(m_al_1[KFitConst::kNumber7 * m_TrackCount + 0][0]);
520  m_AfterVertex.setY(m_al_1[KFitConst::kNumber7 * m_TrackCount + 1][0]);
521  m_AfterVertex.setZ(m_al_1[KFitConst::kNumber7 * m_TrackCount + 2][0]);
522  // error of the vertex
523  for (int i = 0; i < 3; i++) for (int j = i; j < 3; j++) {
524  m_AfterVertexError[i][j] = m_V_al_1[KFitConst::kNumber7 * m_TrackCount + i][KFitConst::kNumber7 * m_TrackCount + j];
525  }
526  // error between vertex and tracks
527  for (int i = 0; i < m_TrackCount; i++) {
528  HepMatrix hm(3, KFitConst::kNumber7, 0);
529  for (int j = 0; j < 3; j++) for (int k = 0; k < KFitConst::kNumber7; k++) {
530  hm[j][k] = m_V_al_1[KFitConst::kNumber7 * m_TrackCount + j][KFitConst::kNumber7 * i + k];
531  }
532  if (m_IsFixMass[i])
533  m_AfterTrackVertexError.push_back(this->makeError4(m_Tracks[i].getMomentum(), hm));
534  else
535  m_AfterTrackVertexError.push_back(hm);
536  }
537  } else {
538  // not fit
539  m_AfterVertex = m_BeforeVertex;
540  }
541 
542  return m_ErrorCode = KFitError::kNoError;
543 }
544 
545 
546 enum KFitError::ECode
547 FourCFitKFit::makeCoreMatrix() {
548  if (!m_FlagFitIncludingVertex)
549  {
550 
551  HepMatrix al_1_prime(m_al_1);
552  HepMatrix Sum_al_1(4, 1, 0);
553  double energy[KFitConst::kMaxTrackCount2];
554 
555  for (int i = 0; i < m_TrackCount; i++) {
556  energy[i] = sqrt(al_1_prime[i * KFitConst::kNumber7 + 0][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] +
557  al_1_prime[i * KFitConst::kNumber7 + 1][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] +
558  al_1_prime[i * KFitConst::kNumber7 + 2][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] +
559  m_property[i][1] * m_property[i][1]);
560  }
561 
562  for (int i = 0; i < m_TrackCount; i++) {
563  // 3->4
564  for (int j = 0; j < 4; j++) Sum_al_1[j][0] += al_1_prime[i * KFitConst::kNumber7 + j][0];
565  }
566 
567  m_d[0][0] = Sum_al_1[0][0] - m_FourMomentum.Px();
568  m_d[1][0] = Sum_al_1[1][0] - m_FourMomentum.Py();
569  m_d[2][0] = Sum_al_1[2][0] - m_FourMomentum.Pz();
570  m_d[3][0] = Sum_al_1[3][0] - m_FourMomentum.E();
571 
572  for (int i = 0; i < m_TrackCount; i++) {
573  if (energy[i] == 0) {
574  m_ErrorCode = KFitError::kDivisionByZero;
575  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
576  break;
577  }
578 
579  for (int l = 0; l < 4; l++) {
580  for (int n = 0; n < 6; n++) {
581  if (l == n) m_D[l][i * KFitConst::kNumber7 + n] = 1;
582  else m_D[l][i * KFitConst::kNumber7 + n] = 0;
583  }
584  }
585  }
586 
587  } else {
588 
589  // m_FlagFitIncludingVertex == true
590  HepMatrix al_1_prime(m_al_1);
591  HepMatrix Sum_al_1(7, 1, 0);
592  double energy[KFitConst::kMaxTrackCount2];
593 
594  for (int i = 0; i < m_TrackCount; i++)
595  {
596  const double a = m_property[i][2];
597  al_1_prime[i * KFitConst::kNumber7 + 0][0] -= a * (al_1_prime[KFitConst::kNumber7 * m_TrackCount + 1][0] - al_1_prime[i *
598  KFitConst::kNumber7 + 5][0]);
599  al_1_prime[i * KFitConst::kNumber7 + 1][0] += a * (al_1_prime[KFitConst::kNumber7 * m_TrackCount + 0][0] - al_1_prime[i *
600  KFitConst::kNumber7 + 4][0]);
601  energy[i] = sqrt(al_1_prime[i * KFitConst::kNumber7 + 0][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] +
602  al_1_prime[i * KFitConst::kNumber7 + 1][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] +
603  al_1_prime[i * KFitConst::kNumber7 + 2][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] +
604  m_property[i][1] * m_property[i][1]);
605  Sum_al_1[6][0] = + a;
606  }
607 
608  for (int i = 0; i < m_TrackCount; i++)
609  {
610  if (energy[i] == 0) {
611  m_ErrorCode = KFitError::kDivisionByZero;
612  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
613  break;
614  }
615 
616  if (m_IsFixMass[i]) {
617  double invE = 1. / energy[i];
618  Sum_al_1[3][0] += energy[i];
619  Sum_al_1[4][0] += al_1_prime[i * KFitConst::kNumber7 + 1][0] * m_property[i][2] * invE;
620  Sum_al_1[5][0] += al_1_prime[i * KFitConst::kNumber7 + 0][0] * m_property[i][2] * invE;
621  } else {
622  Sum_al_1[3][0] += al_1_prime[i * KFitConst::kNumber7 + 3][0];
623  }
624 
625  for (int j = 0; j < 3; j++) Sum_al_1[j][0] += al_1_prime[i * KFitConst::kNumber7 + j][0];
626  }
627 
628  m_d[0][0] = Sum_al_1[0][0] - m_FourMomentum.Px();
629  m_d[1][0] = Sum_al_1[1][0] - m_FourMomentum.Py();
630  m_d[2][0] = Sum_al_1[2][0] - m_FourMomentum.Pz();
631  m_d[3][0] = Sum_al_1[3][0] - m_FourMomentum.E();
632 
633  for (int i = 0; i < m_TrackCount; i++)
634  {
635  if (energy[i] == 0) {
636  m_ErrorCode = KFitError::kDivisionByZero;
637  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
638  break;
639  }
640 
641  for (int l = 0; l < 4; l++) {
642  for (int n = 0; n < 6; n++) {
643  if (l == n) m_D[l][i * KFitConst::kNumber7 + n] = 1;
644  else m_D[l][i * KFitConst::kNumber7 + n] = 0;
645  }
646  }
647  }
648 
649  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]);
650  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]);
651  m_D[0][KFitConst::kNumber7 * m_TrackCount + 2] = 0.;
652  }
653 
654  return m_ErrorCode = KFitError::kNoError;
655 }
656 
657 
658 enum KFitError::ECode
659 FourCFitKFit::calculateNDF() {
660  m_NDF = 4;
661 
662  return m_ErrorCode = KFitError::kNoError;
663 }
664 
665 enum KFitError::ECode FourCFitKFit::updateMother(Particle* mother)
666 {
667  MakeMotherKFit kmm;
668  kmm.setMagneticField(m_MagneticField);
669  unsigned n = getTrackCount();
670  for (unsigned i = 0; i < n; ++i) {
671  kmm.addTrack(getTrackMomentum(i), getTrackPosition(i), getTrackError(i),
672  getTrack(i).getCharge());
673  if (getFlagFitWithVertex())
674  kmm.setTrackVertexError(getTrackVertexError(i));
675  for (unsigned j = i + 1; j < n; ++j) {
676  kmm.setCorrelation(getCorrelation(i, j));
677  }
678  }
679  kmm.setVertex(getVertex());
680  if (getFlagFitWithVertex())
681  kmm.setVertexError(getVertexError());
682  m_ErrorCode = kmm.doMake();
683  if (m_ErrorCode != KFitError::kNoError)
684  return m_ErrorCode;
685  double chi2 = getCHIsq();
686  int ndf = getNDF();
687  double prob = TMath::Prob(chi2, ndf);
688  //
689  bool haschi2 = mother->hasExtraInfo("chiSquared");
690  if (haschi2) {
691  mother->setExtraInfo("chiSquared", chi2);
692  mother->setExtraInfo("ndf", ndf);
693  } else {
694  mother->addExtraInfo("chiSquared", chi2);
695  mother->addExtraInfo("ndf", ndf);
696  }
697 
698  mother->updateMomentum(
699  CLHEPToROOT::getTLorentzVector(kmm.getMotherMomentum()),
700  CLHEPToROOT::getTVector3(kmm.getMotherPosition()),
701  CLHEPToROOT::getTMatrixFSym(kmm.getMotherError()),
702  prob);
703  m_ErrorCode = KFitError::kNoError;
704  return m_ErrorCode;
705 }
Belle2::analysis::MakeMotherKFit::setTrackVertexError
enum KFitError::ECode setTrackVertexError(const CLHEP::HepMatrix &e)
Set a vertex error matrix of the child particle in the addTrack'ed order.
Definition: MakeMotherKFit.cc:95
Belle2::analysis::MakeMotherKFit::setVertex
enum KFitError::ECode setVertex(const HepPoint3D &v)
Set a vertex position of the mother particle.
Definition: MakeMotherKFit.cc:71
Belle2::analysis::MakeMotherKFit::getMotherError
const CLHEP::HepSymMatrix getMotherError(void) const
Get an error matrix of the mother particle.
Definition: MakeMotherKFit.cc:185
Belle2::analysis::MakeMotherKFit
MakeMotherKFit is a class to build mother particle from kinematically fitted daughters.
Definition: MakeMotherKFit.h:41
Belle2::analysis::MakeMotherKFit::getMotherMomentum
const CLHEP::HepLorentzVector getMotherMomentum(void) const
Get a Lorentz vector of the mother particle.
Definition: MakeMotherKFit.cc:171
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::analysis::MakeMotherKFit::setMagneticField
enum KFitError::ECode setMagneticField(const double mf)
Change a magnetic field from the default value KFitConst::kDefaultMagneticField.
Definition: MakeMotherKFit.cc:63
Belle2::analysis::MakeMotherKFit::setCorrelation
enum KFitError::ECode setCorrelation(const CLHEP::HepMatrix &e)
Set a correlation matrix.
Definition: MakeMotherKFit.cc:119
Belle2::analysis::MakeMotherKFit::setVertexError
enum KFitError::ECode setVertexError(const CLHEP::HepSymMatrix &e)
Set a vertex error matrix of the mother particle.
Definition: MakeMotherKFit.cc:79
Belle2::analysis::MakeMotherKFit::getMotherPosition
const HepPoint3D getMotherPosition(void) const
Get a position of the mother particle.
Definition: MakeMotherKFit.cc:178
Belle2::Particle
Class to store reconstructed particles.
Definition: Particle.h:77
Belle2::analysis::MakeMotherKFit::addTrack
enum KFitError::ECode addTrack(const KFitTrack &kp)
Add a track to the make-mother object.
Definition: MakeMotherKFit.cc:40
Belle2::analysis::KFitError::ECode
ECode
ECode is a error code enumerate.
Definition: KFitError.h:43
Belle2::analysis::MakeMotherKFit::doMake
enum KFitError::ECode doMake(void)
Perform a reconstruction of mother particle.
Definition: MakeMotherKFit.cc:192
HepGeom::Point3D< double >