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