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