Belle II Software  release-05-01-25
MassPointingVertexFitKFit.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2019 - Belle II Collaboration *
4  * Contributors: Frank Meier *
5  * *
6  * This software is provided "as is" without any warranty. *
7  **************************************************************************/
8 
9 #include <analysis/VertexFitting/KFit/MakeMotherKFit.h>
10 #include <analysis/VertexFitting/KFit/MassPointingVertexFitKFit.h>
11 #include <analysis/utility/CLHEPToROOT.h>
12 
13 using namespace std;
14 using namespace Belle2;
15 using namespace Belle2::analysis;
16 using namespace CLHEP;
17 
18 MassPointingVertexFitKFit::MassPointingVertexFitKFit():
19  m_BeforeVertex(HepPoint3D(0., 0., 0.)),
20  m_AfterVertexError(HepSymMatrix(3, 0)),
21  m_ProductionVertex(HepPoint3D(0., 0., 0.))
22 {
23  m_FlagFitted = false;
25  m_V_E = HepMatrix(3, 3, 0);
26  m_v = HepMatrix(3, 1, 0);
27  m_v_a = HepMatrix(3, 1, 0);
28  m_InvariantMass = -1.0;
29 }
30 
31 
33 
34 
37  m_BeforeVertex = v;
38 
40 }
41 
42 
45  m_InvariantMass = m;
46 
48 }
49 
50 
54 
56 }
57 
58 
61  m_IsFixMass.push_back(true);
62 
64 }
65 
66 
69  m_IsFixMass.push_back(false);
70 
72 }
73 
74 
77  return KFitBase::setCorrelation(m);
78 }
79 
80 
84 }
85 
86 
87 const HepPoint3D
89 {
90  if (flag == KFitConst::kAfterFit && !isFitted()) return HepPoint3D();
91 
92  switch (flag) {
94  return m_BeforeVertex;
95 
97  return m_AfterVertex;
98 
99  default:
100  KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
101  return HepPoint3D();
102  }
103 }
104 
105 
106 const HepSymMatrix
108 {
109  return m_AfterVertexError;
110 }
111 
112 
113 double
115 {
116  return m_InvariantMass;
117 }
118 
119 
120 double
122 {
123  return m_CHIsq;
124 }
125 
126 
127 const HepMatrix
129 {
130  if (!isTrackIDInRange(id)) return HepMatrix(3, KFitConst::kNumber7, 0);
131 
132  return m_AfterTrackVertexError[id];
133 }
134 
135 
136 double
138 {
139  if (!isFitted()) return -1;
140  if (!isTrackIDInRange(id)) return -1;
141 
142  if (m_IsFixMass[id]) {
143 
144  HepMatrix da(m_Tracks[id].getFitParameter(KFitConst::kBeforeFit) - m_Tracks[id].getFitParameter(KFitConst::kAfterFit));
145  int err_inverse = 0;
146  const double chisq = (da.T() * (m_Tracks[id].getFitError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
147 
148  if (err_inverse) {
149  KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
150  return -1;
151  }
152 
153  return chisq;
154 
155  } else {
156 
157  HepMatrix da(m_Tracks[id].getMomPos(KFitConst::kBeforeFit) - m_Tracks[id].getMomPos(KFitConst::kAfterFit));
158  int err_inverse = 0;
159  const double chisq = (da.T() * (m_Tracks[id].getError(KFitConst::kBeforeFit).inverse(err_inverse)) * da)[0][0];
160 
161  if (err_inverse) {
162  KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kCannotGetMatrixInverse);
163  return -1;
164  }
165 
166  return chisq;
167  }
168 }
169 
170 
171 const HepMatrix
172 MassPointingVertexFitKFit::getCorrelation(const int id1, const int id2, const int flag) const
173 {
174  if (flag == KFitConst::kAfterFit && !isFitted()) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
175  if (!isTrackIDInRange(id1)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
176  if (!isTrackIDInRange(id2)) return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
177 
178  switch (flag) {
180  return KFitBase::getCorrelation(id1, id2, flag);
181 
183  return makeError3(
184  this->getTrackMomentum(id1),
185  this->getTrackMomentum(id2),
186  m_V_al_1.sub(KFitConst::kNumber7 * id1 + 1, KFitConst::kNumber7 * (id1 + 1), KFitConst::kNumber7 * id2 + 1,
187  KFitConst::kNumber7 * (id2 + 1)),
188  m_IsFixMass[id1],
189  m_IsFixMass[id2]);
190 
191  default:
192  KFitError::displayError(__FILE__, __LINE__, __func__, KFitError::kOutOfRange);
193  return HepMatrix(KFitConst::kNumber7, KFitConst::kNumber7, 0);
194  }
195 }
196 
197 
198 enum KFitError::ECode
200  return KFitBase::doFit2();
201 }
202 
203 
204 enum KFitError::ECode
207  {
209  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
210  return m_ErrorCode;
211  }
212 
213 
214  if (m_IsFixMass.size() == 0)
215  {
216  // If no fix_mass flag at all,
217  // all tracks are considered to be fixed at mass.
218  for (int i = 0; i < m_TrackCount; i++) this->fixMass();
219  } else if (m_IsFixMass.size() != (unsigned int)m_TrackCount)
220  {
222  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
223  return m_ErrorCode;
224  }
225 
226 
227  int index = 0;
228  m_al_0 = HepMatrix(KFitConst::kNumber7 * m_TrackCount, 1, 0);
229  m_property = HepMatrix(m_TrackCount, 3, 0);
230  m_V_al_0 = HepSymMatrix(KFitConst::kNumber7 * m_TrackCount, 0);
231 
232  for (auto& track : m_Tracks)
233  {
234  // momentum x,y,z,E and position x,y,z
235  m_al_0[index * KFitConst::kNumber7 + 0][0] = track.getMomentum(KFitConst::kBeforeFit).x();
236  m_al_0[index * KFitConst::kNumber7 + 1][0] = track.getMomentum(KFitConst::kBeforeFit).y();
237  m_al_0[index * KFitConst::kNumber7 + 2][0] = track.getMomentum(KFitConst::kBeforeFit).z();
238  m_al_0[index * KFitConst::kNumber7 + 3][0] = track.getMomentum(KFitConst::kBeforeFit).t();
239  m_al_0[index * KFitConst::kNumber7 + 4][0] = track.getPosition(KFitConst::kBeforeFit).x();
240  m_al_0[index * KFitConst::kNumber7 + 5][0] = track.getPosition(KFitConst::kBeforeFit).y();
241  m_al_0[index * KFitConst::kNumber7 + 6][0] = track.getPosition(KFitConst::kBeforeFit).z();
242  // these error
243  m_V_al_0.sub(index * KFitConst::kNumber7 + 1, track.getError(KFitConst::kBeforeFit));
244  // charge, mass, a
245  m_property[index][0] = track.getCharge();
246  m_property[index][1] = track.getMass();
247  const double c = KFitConst::kLightSpeed; // C++ bug?
248  // m_property[index][2] = - it->getCharge() * KFitConst::kLightSpeed * m_MagneticField;
249  m_property[index][2] = -c * m_MagneticField * track.getCharge();
250  index++;
251  }
252 
253  // error between track and track
254  if (m_FlagCorrelation)
255  {
256  this->prepareCorrelation();
258  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
259  return m_ErrorCode;
260  }
261  }
262 
263  // vertex
264  m_v_a[0][0] = m_BeforeVertex.x();
265  m_v_a[1][0] = m_BeforeVertex.y();
266  m_v_a[2][0] = m_BeforeVertex.z();
267 
268  // set member matrix
269  m_al_1 = m_al_0;
270 
272  m_D = m_V_al_1.sub(1, m_TrackCount * 2 + 3, 1, KFitConst::kNumber7 * m_TrackCount);
273  m_E = m_V_al_1.sub(1, m_TrackCount * 2 + 3, 1, 3);
274  m_d = m_V_al_1.sub(1, m_TrackCount * 2 + 3, 1, 1);
275  m_V_D = m_V_al_1.sub(1, m_TrackCount * 2 + 3, 1, m_TrackCount * 2 + 3);
276  m_lam = m_V_al_1.sub(1, m_TrackCount * 2 + 3, 1, 1);
277  m_lam0 = m_V_al_1.sub(1, m_TrackCount * 2 + 3, 1, 1);
278  m_V_Dt = m_V_al_1.sub(1, m_TrackCount * 2 + 3, 1, m_TrackCount * 2 + 3);
280 
282 }
283 
284 
285 enum KFitError::ECode
287  // vertex
288  for (int i = 0; i < 3; i++) m_v[i][0] = m_v_a[i][0];
289 
291 }
292 
293 
294 enum KFitError::ECode
296  if (m_BeforeCorrelation.size() != static_cast<unsigned int>(m_TrackCount * (m_TrackCount - 1) / 2))
297  {
299  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
300  return m_ErrorCode;
301  }
302 
303  int row = 0, col = 0;
304 
305  for (auto& hm : m_BeforeCorrelation)
306  {
307  // counter
308  row++;
309  if (row == m_TrackCount) {
310  col++;
311  row = col + 1;
312  }
313 
314  int ii = 0, jj = 0;
315  for (int i = KFitConst::kNumber7 * row; i < KFitConst::kNumber7 * (row + 1); i++) {
316  for (int j = KFitConst::kNumber7 * col; j < KFitConst::kNumber7 * (col + 1); j++) {
317  m_V_al_0[i][j] = hm[ii][jj];
318  jj++;
319  }
320  jj = 0;
321  ii++;
322  }
323  }
324 
326 }
327 
328 
329 enum KFitError::ECode
331  Hep3Vector h3v;
332  int index = 0;
333  for (auto& pdata : m_Tracks)
334  {
335  // tracks
336  // momentum
337  h3v.setX(m_al_1[index * KFitConst::kNumber7 + 0][0]);
338  h3v.setY(m_al_1[index * KFitConst::kNumber7 + 1][0]);
339  h3v.setZ(m_al_1[index * KFitConst::kNumber7 + 2][0]);
340  if (m_IsFixMass[index])
341  pdata.setMomentum(HepLorentzVector(h3v, sqrt(h3v.mag2() + pdata.getMass()*pdata.getMass())), KFitConst::kAfterFit);
342  else
343  pdata.setMomentum(HepLorentzVector(h3v, m_al_1[index * KFitConst::kNumber7 + 3][0]), KFitConst::kAfterFit);
344  // position
345  pdata.setPosition(HepPoint3D(
346  m_al_1[index * KFitConst::kNumber7 + 4][0],
347  m_al_1[index * KFitConst::kNumber7 + 5][0],
349  // error of the tracks
350  pdata.setError(this->makeError3(pdata.getMomentum(),
351  m_V_al_1.sub(
352  index * KFitConst::kNumber7 + 1,
353  (index + 1)*KFitConst::kNumber7,
354  index * KFitConst::kNumber7 + 1,
355  (index + 1)*KFitConst::kNumber7), m_IsFixMass[index]),
357  if (m_ErrorCode != KFitError::kNoError) break;
358  index++;
359  }
360 
361  // vertex
362  m_AfterVertex.setX(m_v_a[0][0]);
363  m_AfterVertex.setY(m_v_a[1][0]);
364  m_AfterVertex.setZ(m_v_a[2][0]);
365  // error of the vertex
366  for (int i = 0; i < 3; i++) for (int j = i; j < 3; j++)
367  {
368  m_AfterVertexError[i][j] = m_V_E[i][j];
369  }
370  // error between vertex and tracks
371  for (int i = 0; i < m_TrackCount; i++)
372  {
373  HepMatrix hm(3, KFitConst::kNumber7, 0);
374  for (int j = 0; j < 3; j++) for (int k = 0; k < KFitConst::kNumber7; k++) {
375  hm[j][k] = m_Cov_v_al_1[j][KFitConst::kNumber7 * i + k];
376  }
377  if (m_IsFixMass[i])
378  m_AfterTrackVertexError.push_back(this->makeError4(m_Tracks[i].getMomentum(), hm));
379  else
380  m_AfterTrackVertexError.push_back(hm);
381  }
382 
384 }
385 
386 
387 enum KFitError::ECode
389  // Mass Constraint
390  HepMatrix al_1_prime(m_al_1);
391  HepMatrix Sum_al_1(4, 1, 0);
392  double energy[KFitConst::kMaxTrackCount2];
393  double a;
394 
395  for (int i = 0; i < m_TrackCount; i++)
396  {
397  a = m_property[i][2];
398  al_1_prime[i * KFitConst::kNumber7 + 0][0] -= a * (m_v_a[1][0] - al_1_prime[i * KFitConst::kNumber7 + 5][0]);
399  al_1_prime[i * KFitConst::kNumber7 + 1][0] += a * (m_v_a[0][0] - al_1_prime[i * KFitConst::kNumber7 + 4][0]);
400  energy[i] = sqrt(al_1_prime[i * KFitConst::kNumber7 + 0][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] +
401  al_1_prime[i * KFitConst::kNumber7 + 1][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] +
402  al_1_prime[i * KFitConst::kNumber7 + 2][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] +
403  m_property[i][1] * m_property[i][1]);
404  }
405 
406  for (int i = 0; i < m_TrackCount; i++)
407  {
408  if (m_IsFixMass[i])
409  Sum_al_1[3][0] += energy[i];
410  else
411  Sum_al_1[3][0] += al_1_prime[i * KFitConst::kNumber7 + 3][0];
412 
413  for (int j = 0; j < 3; j++) Sum_al_1[j][0] += al_1_prime[i * KFitConst::kNumber7 + j][0];
414  }
415 
416  double vtx = sqrt((m_v_a[0][0] - m_ProductionVertex.x()) * (m_v_a[0][0] - m_ProductionVertex.x())
417  + (m_v_a[1][0] - m_ProductionVertex.y()) * (m_v_a[1][0] - m_ProductionVertex.y())
418  + (m_v_a[2][0] - m_ProductionVertex.z()) * (m_v_a[2][0] - m_ProductionVertex.z()));
419  double vtxpt2 = (m_ProductionVertex.x() - m_v_a[0][0]) * (m_ProductionVertex.x() - m_v_a[0][0]) + (m_ProductionVertex.y() - m_v_a[1][0]) * (m_ProductionVertex.y() - m_v_a[1][0]);
420  double mom = sqrt(Sum_al_1[0][0] * Sum_al_1[0][0] + Sum_al_1[1][0] * Sum_al_1[1][0] + Sum_al_1[2][0] * Sum_al_1[2][0]);
421  double Pt2 = Sum_al_1[0][0] * Sum_al_1[0][0] + Sum_al_1[1][0] * Sum_al_1[1][0];
422 
423  m_d[2 * m_TrackCount][0] =
424  + Sum_al_1[3][0] * Sum_al_1[3][0] - Sum_al_1[0][0] * Sum_al_1[0][0]
425  - Sum_al_1[1][0] * Sum_al_1[1][0] - Sum_al_1[2][0] * Sum_al_1[2][0]
427 
428  double phiPointingConstraint = atan2(m_v_a[1][0] - m_ProductionVertex.y(), m_v_a[0][0] - m_ProductionVertex.x()) - atan2(Sum_al_1[1][0], Sum_al_1[0][0]);
429  phiPointingConstraint = std::fmod(phiPointingConstraint + TMath::Pi(), TMath::TwoPi());
430  if (phiPointingConstraint < 0) phiPointingConstraint += TMath::TwoPi();
431  phiPointingConstraint -= TMath::Pi();
432 
433  m_d[2 * m_TrackCount + 1][0] = phiPointingConstraint;
434  m_d[2 * m_TrackCount + 2][0] = acos((m_v_a[2][0] - m_ProductionVertex.z()) / vtx) - acos(Sum_al_1[2][0] / mom);
435 
436  double Sum_a = 0., Sum_tmpx = 0., Sum_tmpy = 0.;
437  for (int i = 0; i < m_TrackCount; i++)
438  {
439  if (energy[i] == 0) {
441  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
442  return m_ErrorCode;
443  }
444 
445  a = m_property[i][2];
446 
447  if (m_IsFixMass[i]) {
448  double invE = 1. / energy[i];
449  m_D[2 * m_TrackCount][i * KFitConst::kNumber7 + 0] = 2.*(Sum_al_1[3][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE -
450  Sum_al_1[0][0]);
451  m_D[2 * m_TrackCount][i * KFitConst::kNumber7 + 1] = 2.*(Sum_al_1[3][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE -
452  Sum_al_1[1][0]);
453  m_D[2 * m_TrackCount][i * KFitConst::kNumber7 + 2] = 2.*(Sum_al_1[3][0] * al_1_prime[i * KFitConst::kNumber7 + 2][0] * invE -
454  Sum_al_1[2][0]);
455  m_D[2 * m_TrackCount][i * KFitConst::kNumber7 + 3] = 0.;
456  m_D[2 * m_TrackCount][i * KFitConst::kNumber7 + 4] = -2.*(Sum_al_1[3][0] * al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE -
457  Sum_al_1[1][0]) * a;
458  m_D[2 * m_TrackCount][i * KFitConst::kNumber7 + 5] = 2.*(Sum_al_1[3][0] * al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE -
459  Sum_al_1[0][0]) * a;
460  m_D[2 * m_TrackCount][i * KFitConst::kNumber7 + 6] = 0.;
461  Sum_tmpx += al_1_prime[i * KFitConst::kNumber7 + 0][0] * invE * a;
462  Sum_tmpy += al_1_prime[i * KFitConst::kNumber7 + 1][0] * invE * a;
463  } else {
464  m_D[2 * m_TrackCount][i * KFitConst::kNumber7 + 0] = -2.*Sum_al_1[0][0];
465  m_D[2 * m_TrackCount][i * KFitConst::kNumber7 + 1] = -2.*Sum_al_1[1][0];
466  m_D[2 * m_TrackCount][i * KFitConst::kNumber7 + 2] = -2.*Sum_al_1[2][0];
467  m_D[2 * m_TrackCount][i * KFitConst::kNumber7 + 3] = 2.*Sum_al_1[3][0];
468  m_D[2 * m_TrackCount][i * KFitConst::kNumber7 + 4] = 2.*Sum_al_1[1][0] * a;
469  m_D[2 * m_TrackCount][i * KFitConst::kNumber7 + 5] = -2.*Sum_al_1[0][0] * a;
470  m_D[2 * m_TrackCount][i * KFitConst::kNumber7 + 6] = 0.;
471  }
472  m_D[2 * m_TrackCount + 1][i * KFitConst::kNumber7 + 0] = Sum_al_1[1][0] / Pt2;
473  m_D[2 * m_TrackCount + 1][i * KFitConst::kNumber7 + 1] = - Sum_al_1[0][0] / Pt2;
474  m_D[2 * m_TrackCount + 1][i * KFitConst::kNumber7 + 2] = 0.;
475  m_D[2 * m_TrackCount + 1][i * KFitConst::kNumber7 + 3] = 0.;
476  m_D[2 * m_TrackCount + 1][i * KFitConst::kNumber7 + 4] = a * Sum_al_1[0][0] / Pt2;
477  m_D[2 * m_TrackCount + 1][i * KFitConst::kNumber7 + 5] = a * Sum_al_1[1][0] / Pt2;
478  m_D[2 * m_TrackCount + 1][i * KFitConst::kNumber7 + 6] = 0.;
479  m_D[2 * m_TrackCount + 2][i * KFitConst::kNumber7 + 0] = - Sum_al_1[0][0] * Sum_al_1[2][0] / (sqrt(Pt2) * mom * mom);
480  m_D[2 * m_TrackCount + 2][i * KFitConst::kNumber7 + 1] = - Sum_al_1[1][0] * Sum_al_1[2][0] / (sqrt(Pt2) * mom * mom);
481  m_D[2 * m_TrackCount + 2][i * KFitConst::kNumber7 + 2] = sqrt(Pt2) / (mom * mom);
482  m_D[2 * m_TrackCount + 2][i * KFitConst::kNumber7 + 3] = 0.;
483  m_D[2 * m_TrackCount + 2][i * KFitConst::kNumber7 + 4] = a * Sum_al_1[1][0] * Sum_al_1[2][0] / (sqrt(Pt2) * mom * mom);
484  m_D[2 * m_TrackCount + 2][i * KFitConst::kNumber7 + 5] = - a * Sum_al_1[0][0] * Sum_al_1[2][0] / (sqrt(Pt2) * mom * mom);
485  m_D[2 * m_TrackCount + 2][i * KFitConst::kNumber7 + 6] = 0.;
486  Sum_a += a;
487  }
488 
489  // m_E
490  m_E[2 * m_TrackCount][0] = -2.*Sum_al_1[1][0] * Sum_a + 2.*Sum_al_1[3][0] * Sum_tmpy;
491  m_E[2 * m_TrackCount][1] = 2.*Sum_al_1[0][0] * Sum_a - 2.*Sum_al_1[3][0] * Sum_tmpx;
492  m_E[2 * m_TrackCount][2] = 0.;
493  m_E[2 * m_TrackCount + 1][0] = (m_ProductionVertex.y() - m_v_a[1][0]) / vtxpt2 - Sum_a * Sum_al_1[0][0] / Pt2;
494  m_E[2 * m_TrackCount + 1][1] = (m_v_a[0][0] - m_ProductionVertex.x()) / vtxpt2 - Sum_a * Sum_al_1[1][0] / Pt2;
495  m_E[2 * m_TrackCount + 1][2] = 0.;
496  m_E[2 * m_TrackCount + 2][0] = (m_v_a[0][0] - m_ProductionVertex.x()) * (m_v_a[2][0] - m_ProductionVertex.z()) / (vtx* vtx * sqrt(vtxpt2)) - Sum_a * Sum_al_1[1][0] * Sum_al_1[2][0] / (sqrt(Pt2) * mom * mom);
497  m_E[2 * m_TrackCount + 2][1] = (m_v_a[1][0] - m_ProductionVertex.y()) * (m_v_a[2][0] - m_ProductionVertex.z()) / (vtx* vtx * sqrt(vtxpt2)) + Sum_a * Sum_al_1[0][0] * Sum_al_1[2][0] / (sqrt(Pt2) * mom * mom);
498  m_E[2 * m_TrackCount + 2][2] = - sqrt(vtxpt2) / (vtx * vtx);
499 
500  for (int i = 0; i < m_TrackCount; i++)
501  {
502  double B, S, U;
503  double sininv, sqrtag;
504 
505  double px = m_al_1[i * KFitConst::kNumber7 + 0][0];
506  double py = m_al_1[i * KFitConst::kNumber7 + 1][0];
507  double pz = m_al_1[i * KFitConst::kNumber7 + 2][0];
508  double x = m_al_1[i * KFitConst::kNumber7 + 4][0];
509  double y = m_al_1[i * KFitConst::kNumber7 + 5][0];
510  double z = m_al_1[i * KFitConst::kNumber7 + 6][0];
511  a = m_property[i][2];
512 
513  double pt = sqrt(px * px + py * py);
514 
515  if (pt == 0) {
517  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
518  return m_ErrorCode;
519  }
520 
521  double invPt = 1. / pt;
522  double invPt2 = invPt * invPt;
523  double dlx = m_v_a[0][0] - x;
524  double dly = m_v_a[1][0] - y;
525  double dlz = m_v_a[2][0] - z;
526  double a1 = -dlx * py + dly * px;
527  double a2 = dlx * px + dly * py;
528  double r2d2 = dlx * dlx + dly * dly;
529  double Rx = dlx - 2.*px * a2 * invPt2;
530  double Ry = dly - 2.*py * a2 * invPt2;
531 
532  if (a != 0.) { // charged
533 
534  B = a * a2 * invPt2;
535  if (fabs(B) > 1.) {
537  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
538  return m_ErrorCode;
539  }
540  // sin^(-1)(B)
541  sininv = asin(B);
542  double tmp0 = 1.0 - B * B;
543  if (tmp0 == 0) {
545  KFitError::displayError(__FILE__, __LINE__, __func__, m_ErrorCode);
546  return m_ErrorCode;
547  }
548  // 1/sqrt(1-B^2)
549  sqrtag = 1.0 / sqrt(tmp0);
550  S = sqrtag * invPt2;
551  U = dlz - pz * sininv / a;
552 
553  } else { // neutral
554 
555  B = 0.0;
556  sininv = 0.0;
557  sqrtag = 1.0;
558  S = invPt2;
559  U = dlz - pz * a2 * invPt2;
560 
561  }
562 
563  // d
564  m_d[i * 2 + 0][0] = a1 - 0.5 * a * r2d2;
565  m_d[i * 2 + 1][0] = U * pt;
566 
567  // D
568  m_D[i * 2 + 0][i * KFitConst::kNumber7 + 0] = dly;
569  m_D[i * 2 + 0][i * KFitConst::kNumber7 + 1] = -dlx;
570  m_D[i * 2 + 0][i * KFitConst::kNumber7 + 2] = 0.0;
571  m_D[i * 2 + 0][i * KFitConst::kNumber7 + 4] = py + a * dlx;
572  m_D[i * 2 + 0][i * KFitConst::kNumber7 + 5] = -px + a * dly;
573  m_D[i * 2 + 0][i * KFitConst::kNumber7 + 6] = 0.0;
574  m_D[i * 2 + 1][i * KFitConst::kNumber7 + 0] = -pz * pt * S * Rx + U * px * invPt;
575  m_D[i * 2 + 1][i * KFitConst::kNumber7 + 1] = -pz * pt * S * Ry + U * py * invPt;
576  if (a != 0.)
577  m_D[i * 2 + 1][i * KFitConst::kNumber7 + 2] = -sininv * pt / a;
578  else
579  m_D[i * 2 + 1][i * KFitConst::kNumber7 + 2] = -a2 * invPt;
580  m_D[i * 2 + 1][i * KFitConst::kNumber7 + 4] = px * pz * pt * S;
581  m_D[i * 2 + 1][i * KFitConst::kNumber7 + 5] = py * pz * pt * S;
582  m_D[i * 2 + 1][i * KFitConst::kNumber7 + 6] = -pt;
583 
584  // E
585  m_E[i * 2 + 0][0] = -py - a * dlx;
586  m_E[i * 2 + 0][1] = px - a * dly;
587  m_E[i * 2 + 0][2] = 0.0;
588  m_E[i * 2 + 1][0] = -px * pz * pt * S;
589  m_E[i * 2 + 1][1] = -py * pz * pt * S;
590  m_E[i * 2 + 1][2] = pt;
591  }
592 
594 }
595 
596 
597 enum KFitError::ECode
599  m_NDF = 2 * m_TrackCount - 3 + 3;
600 
602 }
603 
605 {
606  MakeMotherKFit kmm;
608  unsigned n = getTrackCount();
609  for (unsigned i = 0; i < n; ++i) {
611  getTrack(i).getCharge());
613  for (unsigned j = i + 1; j < n; ++j) {
614  kmm.setCorrelation(getCorrelation(i, j));
615  }
616  }
617  kmm.setVertex(getVertex());
619  m_ErrorCode = kmm.doMake();
621  B2DEBUG(2, "Error in doMake");
622  return m_ErrorCode;
623  }
624  double chi2 = getCHIsq();
625  int ndf = getNDF();
626  double prob = TMath::Prob(chi2, ndf);
627  mother->addExtraInfo("ndf", ndf);
628  mother->addExtraInfo("chiSquared", chi2);
629  mother->updateMomentum(
630  CLHEPToROOT::getTLorentzVector(kmm.getMotherMomentum()),
631  CLHEPToROOT::getTVector3(kmm.getMotherPosition()),
632  CLHEPToROOT::getTMatrixFSym(kmm.getMotherError()),
633  prob);
635  return m_ErrorCode;
636 }
Belle2::analysis::KFitBase::m_MagneticField
double m_MagneticField
Magnetic field.
Definition: KFitBase.h:311
Belle2::analysis::MassPointingVertexFitKFit::m_BeforeVertex
HepPoint3D m_BeforeVertex
Vertex position before the fit.
Definition: MassPointingVertexFitKFit.h:121
Belle2::analysis::KFitBase::m_v
CLHEP::HepMatrix m_v
See J.Tanaka Ph.D (2001) p137 for definition.
Definition: KFitBase.h:285
Belle2::analysis::MassPointingVertexFitKFit::setProductionVertex
enum KFitError::ECode setProductionVertex(const HepPoint3D &v)
Set the production vertex of the particle.
Definition: MassPointingVertexFitKFit.cc:52
Belle2::analysis::KFitBase::makeError3
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.
Belle2::analysis::KFitBase::m_CHIsq
double m_CHIsq
chi-square of the fit.
Definition: KFitBase.h:297
Belle2::analysis::KFitBase::m_NDF
int m_NDF
NDF of the fit.
Definition: KFitBase.h:295
Belle2::analysis::MassPointingVertexFitKFit::m_AfterTrackVertexError
std::vector< CLHEP::HepMatrix > m_AfterTrackVertexError
array of vertex error matrices after the fit.
Definition: MassPointingVertexFitKFit.h:127
Belle2::analysis::KFitBase::m_ErrorCode
enum KFitError::ECode m_ErrorCode
Error code.
Definition: KFitBase.h:243
Belle2::analysis::MassPointingVertexFitKFit::~MassPointingVertexFitKFit
~MassPointingVertexFitKFit(void)
Destruct the object.
Belle2::analysis::KFitBase::m_FlagCorrelation
bool m_FlagCorrelation
Flag whether a correlation among tracks exists.
Definition: KFitBase.h:306
Belle2::analysis::KFitBase::m_D
CLHEP::HepMatrix m_D
See J.Tanaka Ph.D (2001) p137 for definition.
Definition: KFitBase.h:266
Belle2::analysis::MassPointingVertexFitKFit::getCorrelation
const CLHEP::HepMatrix getCorrelation(const int id1, const int id2, const int flag=KFitConst::kAfterFit) const override
Get a correlation matrix between two tracks.
Definition: MassPointingVertexFitKFit.cc:172
Belle2::analysis::KFitConst::kBeforeFit
static const int kBeforeFit
Input parameter to specify before-fit when setting/getting a track attribute.
Definition: KFitConst.h:43
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::KFitBase::getCorrelation
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:184
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::KFitBase::isFitted
bool isFitted(void) const
Return false if fit is not performed yet or performed fit is failed; otherwise true.
Definition: KFitBase.cc:728
Belle2::analysis::KFitBase::doFit2
enum KFitError::ECode doFit2(void)
Perform a fit (used in VertexFitKFit::doFit() and MassVertexFitKFit::doFit()).
Definition: KFitBase.cc:578
Belle2::analysis::KFitBase::m_NecessaryTrackCount
int m_NecessaryTrackCount
Number needed tracks to perform fit.
Definition: KFitBase.h:303
Belle2::analysis::MassPointingVertexFitKFit::prepareCorrelation
enum KFitError::ECode prepareCorrelation(void) override
Build a grand correlation matrix from input-track properties.
Definition: MassPointingVertexFitKFit.cc:295
Belle2::analysis::KFitBase::getTrackError
const CLHEP::HepSymMatrix getTrackError(const int id) const
Get an error matrix of the track.
Definition: KFitBase.cc:169
Belle2::analysis::KFitBase::m_lam0
CLHEP::HepMatrix m_lam0
See J.Tanaka Ph.D (2001) p138 for definition.
Definition: KFitBase.h:283
Belle2::analysis::KFitBase::m_Cov_v_al_1
CLHEP::HepMatrix m_Cov_v_al_1
See J.Tanaka Ph.D (2001) p137 for definition.
Definition: KFitBase.h:291
Belle2::analysis::KFitBase::m_Tracks
std::vector< KFitTrack > m_Tracks
Container of input tracks.
Definition: KFitBase.h:249
Belle2::analysis::KFitBase::m_property
CLHEP::HepMatrix m_property
Container of charges and masses.
Definition: KFitBase.h:263
Belle2::analysis::MassPointingVertexFitKFit::setInvariantMass
enum KFitError::ECode setInvariantMass(const double m)
Set an invariant mass for the mass-vertex-pointing constraint fit.
Definition: MassPointingVertexFitKFit.cc:44
Belle2::analysis::KFitBase::m_lam
CLHEP::HepMatrix m_lam
See J.Tanaka Ph.D (2001) p137 for definition.
Definition: KFitBase.h:276
Belle2::analysis::KFitBase::m_V_D
CLHEP::HepMatrix m_V_D
See J.Tanaka Ph.D (2001) p138 for definition.
Definition: KFitBase.h:271
Belle2::analysis::KFitBase::m_al_1
CLHEP::HepMatrix m_al_1
See J.Tanaka Ph.D (2001) p136 for definition.
Definition: KFitBase.h:259
Belle2::analysis::MassPointingVertexFitKFit::prepareInputSubMatrix
enum KFitError::ECode prepareInputSubMatrix(void) override
Build sub-matrices for minimum search from input-track properties.
Definition: MassPointingVertexFitKFit.cc:286
Belle2::analysis::KFitError::displayError
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:79
Belle2::analysis::KFitError::kBadCorrelationSize
@ kBadCorrelationSize
Wrong correlation matrix size (internal error)
Definition: KFitError.h:58
Belle2::analysis::KFitBase::m_E
CLHEP::HepMatrix m_E
See J.Tanaka Ph.D (2001) p137 for definition.
Definition: KFitBase.h:279
Belle2::analysis::MakeMotherKFit::getMotherError
const CLHEP::HepSymMatrix getMotherError(void) const
Get an error matrix of the mother particle.
Definition: MakeMotherKFit.cc:185
Belle2::analysis::KFitBase::setZeroCorrelation
virtual enum KFitError::ECode setZeroCorrelation(void)
Indicate no correlation between tracks.
Definition: KFitBase.cc:86
Belle2::analysis::MassPointingVertexFitKFit::calculateNDF
enum KFitError::ECode calculateNDF(void) override
Calculate an NDF of the fit.
Definition: MassPointingVertexFitKFit.cc:598
Belle2::analysis::KFitBase::m_FlagFitted
bool m_FlagFitted
Flag to indicate if the fit is performed and succeeded.
Definition: KFitBase.h:245
Belle2::analysis::MassPointingVertexFitKFit::updateMother
enum KFitError::ECode updateMother(Particle *mother)
Update mother particle.
Definition: MassPointingVertexFitKFit.cc:604
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::analysis::KFitError::kNoError
@ kNoError
No error.
Definition: KFitError.h:45
Belle2::analysis::MassPointingVertexFitKFit::m_AfterVertex
HepPoint3D m_AfterVertex
Vertex position after the fit.
Definition: MassPointingVertexFitKFit.h:123
Belle2::analysis::KFitBase::getTrackCount
int getTrackCount(void) const
Get the number of added tracks.
Definition: KFitBase.cc:108
Belle2::analysis::KFitConst::kAfterFit
static const int kAfterFit
Input parameter to specify after-fit when setting/getting a track attribute.
Definition: KFitConst.h:45
Belle2::analysis::KFitBase::m_v_a
CLHEP::HepMatrix m_v_a
See J.Tanaka Ph.D (2001) p137 for definition.
Definition: KFitBase.h:287
Belle2::analysis::KFitConst::kMaxTrackCount
static const int kMaxTrackCount
Maximum track size.
Definition: KFitConst.h:48
Belle2::analysis::KFitConst::kMaxTrackCount2
static const int kMaxTrackCount2
Maximum track size (internal use)
Definition: KFitConst.h:50
Belle2::analysis::MassPointingVertexFitKFit::setZeroCorrelation
enum KFitError::ECode setZeroCorrelation(void) override
Indicate no correlation between tracks.
Definition: MassPointingVertexFitKFit.cc:82
Belle2::analysis::KFitBase::m_V_al_1
CLHEP::HepMatrix m_V_al_1
See J.Tanaka Ph.D (2001) p138 for definition.
Definition: KFitBase.h:274
Belle2::analysis::MassPointingVertexFitKFit::getTrackVertexError
const CLHEP::HepMatrix getTrackVertexError(const int id) const
Get a vertex error matrix of the track.
Definition: MassPointingVertexFitKFit.cc:128
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::analysis::MassPointingVertexFitKFit::setCorrelation
enum KFitError::ECode setCorrelation(const CLHEP::HepMatrix &m) override
Set a correlation matrix.
Definition: MassPointingVertexFitKFit.cc:76
Belle2::analysis::KFitBase::getNDF
virtual int getNDF(void) const
Get an NDF of the fit.
Definition: KFitBase.cc:115
Belle2::analysis::KFitError::kDivisionByZero
@ kDivisionByZero
Division by zero (bad track property or internal error)
Definition: KFitError.h:63
Belle2::analysis::KFitBase::setCorrelation
virtual enum KFitError::ECode setCorrelation(const CLHEP::HepMatrix &c)
Set a correlation matrix.
Definition: KFitBase.cc:71
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::KFitBase::m_TrackCount
int m_TrackCount
Number of tracks.
Definition: KFitBase.h:301
Belle2::analysis::MassPointingVertexFitKFit::prepareOutputMatrix
enum KFitError::ECode prepareOutputMatrix(void) override
Build an output error matrix.
Definition: MassPointingVertexFitKFit.cc:330
Belle2::analysis::KFitError::kCannotGetMatrixInverse
@ kCannotGetMatrixInverse
Cannot calculate matrix inverse (bad track property or internal error)
Definition: KFitError.h:65
Belle2::analysis::MassPointingVertexFitKFit::m_ProductionVertex
HepPoint3D m_ProductionVertex
Production vertex position.
Definition: MassPointingVertexFitKFit.h:129
Belle2::analysis::KFitBase::m_al_0
CLHEP::HepMatrix m_al_0
See J.Tanaka Ph.D (2001) p136 for definition.
Definition: KFitBase.h:257
Belle2::analysis::KFitConst::kLightSpeed
static constexpr double kLightSpeed
Speed of light.
Definition: KFitConst.h:62
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::MassPointingVertexFitKFit::doFit
enum KFitError::ECode doFit(void)
Perform a mass-vertex-pointing constraint fit.
Definition: MassPointingVertexFitKFit.cc:199
Belle2::analysis::KFitBase::getTrack
const KFitTrack getTrack(const int id) const
Get a specified track object.
Definition: KFitBase.cc:176
Belle2::analysis::KFitBase::m_d
CLHEP::HepMatrix m_d
See J.Tanaka Ph.D (2001) p137 for definition.
Definition: KFitBase.h:268
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::MassPointingVertexFitKFit::m_IsFixMass
std::vector< int > m_IsFixMass
Array of flags whether the track property is fixed at the mass.
Definition: MassPointingVertexFitKFit.h:133
Belle2::analysis::KFitBase::m_V_Dt
CLHEP::HepMatrix m_V_Dt
See J.Tanaka Ph.D (2001) p138 for definition.
Definition: KFitBase.h:289
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::KFitBase::isTrackIDInRange
bool isTrackIDInRange(const int id) const
Check if the id is in the range.
Definition: KFitBase.cc:739
Belle2::analysis::MassPointingVertexFitKFit::unfixMass
enum KFitError::ECode unfixMass(void)
Tell the object to unfix the last added track property at the invariant mass.
Definition: MassPointingVertexFitKFit.cc:68
Belle2::analysis::KFitBase::getTrackPosition
const HepPoint3D getTrackPosition(const int id) const
Get a position of the track.
Definition: KFitBase.cc:162
Belle2::analysis::KFitError::ECode
ECode
ECode is a error code enumerate.
Definition: KFitError.h:43
Belle2::analysis::MassPointingVertexFitKFit::m_AfterVertexError
CLHEP::HepSymMatrix m_AfterVertexError
Vertex error matrix after the fit.
Definition: MassPointingVertexFitKFit.h:125
Belle2::analysis::MakeMotherKFit::doMake
enum KFitError::ECode doMake(void)
Perform a reconstruction of mother particle.
Definition: MakeMotherKFit.cc:192
Belle2::analysis::KFitError::kOutOfRange
@ kOutOfRange
Specified track-id out of range.
Definition: KFitError.h:51
Belle2::analysis::MassPointingVertexFitKFit::getInvariantMass
double getInvariantMass(void) const
Get an invariant mass.
Definition: MassPointingVertexFitKFit.cc:114
HepGeom::Point3D< double >
Belle2::analysis::KFitBase::makeError4
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
Belle2::analysis::MassPointingVertexFitKFit::getVertexError
const CLHEP::HepSymMatrix getVertexError(void) const
Get a fitted vertex error matrix.
Definition: MassPointingVertexFitKFit.cc:107
Belle2::analysis::MassPointingVertexFitKFit::prepareInputMatrix
enum KFitError::ECode prepareInputMatrix(void) override
Build grand matrices for minimum search from input-track properties.
Definition: MassPointingVertexFitKFit.cc:205
Belle2::analysis::MassPointingVertexFitKFit::fixMass
enum KFitError::ECode fixMass(void)
Tell the object to fix the last added track property at the invariant mass.
Definition: MassPointingVertexFitKFit.cc:60
Belle2::analysis::KFitConst::kNumber7
static const int kNumber7
Constant 7 to check matrix size (internal use)
Definition: KFitConst.h:40
Belle2::analysis::MassPointingVertexFitKFit::getTrackCHIsq
double getTrackCHIsq(const int id) const override
Get a chi-square of the track.
Definition: MassPointingVertexFitKFit.cc:137
Belle2::analysis::KFitBase::m_V_E
CLHEP::HepMatrix m_V_E
See J.Tanaka Ph.D (2001) p138 for definition.
Definition: KFitBase.h:281
Belle2::analysis::KFitError::kBadTrackSize
@ kBadTrackSize
Track count too small to perform fit.
Definition: KFitError.h:54
Belle2::analysis::KFitBase::m_V_al_0
CLHEP::HepSymMatrix m_V_al_0
See J.Tanaka Ph.D (2001) p137 for definition.
Definition: KFitBase.h:255
Belle2::analysis::MassPointingVertexFitKFit::getVertex
const HepPoint3D getVertex(const int flag=KFitConst::kAfterFit) const
Get a vertex position.
Definition: MassPointingVertexFitKFit.cc:88
Belle2::analysis::MassPointingVertexFitKFit::getCHIsq
double getCHIsq(void) const override
Get a chi-square of the fit.
Definition: MassPointingVertexFitKFit.cc:121
Belle2::analysis::KFitBase::m_BeforeCorrelation
std::vector< CLHEP::HepMatrix > m_BeforeCorrelation
Container of input correlation matrices.
Definition: KFitBase.h:251
Belle2::analysis::MassPointingVertexFitKFit::m_InvariantMass
double m_InvariantMass
Invariant mass.
Definition: MassPointingVertexFitKFit.h:131
Belle2::analysis::MassPointingVertexFitKFit::setInitialVertex
enum KFitError::ECode setInitialVertex(const HepPoint3D &v)
Set an initial vertex point for the mass-vertex-pointing constraint fit.
Definition: MassPointingVertexFitKFit.cc:36
Belle2::analysis::KFitBase::getTrackMomentum
const CLHEP::HepLorentzVector getTrackMomentum(const int id) const
Get a Lorentz vector of the track.
Definition: KFitBase.cc:155
Belle2::analysis::MassPointingVertexFitKFit::makeCoreMatrix
enum KFitError::ECode makeCoreMatrix(void) override
Build matrices using the kinematical constraint.
Definition: MassPointingVertexFitKFit.cc:388