Belle II Software  release-05-01-25
RaveKinematicVertexFitter.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2013 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Moritz Nadler, Luigi Li Gioi *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 #include <analysis/VertexFitting/RaveInterface/RaveKinematicVertexFitter.h>
12 #include <analysis/VertexFitting/RaveInterface/RaveSetup.h>
13 
14 //root
15 #include <Math/ProbFunc.h>
16 
17 #include <rave/TransientTrackKinematicParticle.h>
18 #include <rave/impl/RaveBase/Converters/interface/RaveStreamers.h>
19 
20 #include <rave/KinematicTreeFactory.h>
21 #include <rave/KinematicConstraint.h>
22 #include <rave/KinematicConstraintBuilder.h>
23 #include <rave/VertexFactory.h>
24 
25 
26 //#include <analysis/dataobjects/Particle.h>
27 
28 //c++ std lib
29 using std::string;
30 #include <vector>
31 using std::vector;
32 #include <iostream>
33 
34 using namespace Belle2;
35 using namespace analysis;
36 
37 
38 
39 RaveKinematicVertexFitter::RaveKinematicVertexFitter(): m_useBeamSpot(false), m_motherParticlePtr(nullptr), m_raveAlgorithm(""),
40  m_massConstFit(false), m_vertFit(true)
41 {
42  if (RaveSetup::getRawInstance() == nullptr) {
43  B2FATAL("RaveSetup::initialize was not called. It has to be called before RaveSetup or RaveKinematicVertexFitter are used");
44  }
45  m_useBeamSpot = RaveSetup::getRawInstance()->m_useBeamSpot;
46 }
47 
48 IOIntercept::InterceptorScopeGuard<IOIntercept::OutputToLogMessages> RaveKinematicVertexFitter::captureOutput()
49 {
50  static IOIntercept::OutputToLogMessages s_captureOutput("Rave", LogConfig::c_Debug, LogConfig::c_Debug);
51  return IOIntercept::start_intercept(s_captureOutput);
52 }
53 
54 
55 RaveKinematicVertexFitter::~RaveKinematicVertexFitter() = default;
56 
57 
58 void RaveKinematicVertexFitter::setMassConstFit(bool isConstFit)
59 {
60  m_massConstFit = isConstFit;
61 }
62 
63 void RaveKinematicVertexFitter::setVertFit(bool isVertFit)
64 {
65  m_vertFit = isVertFit;
66 }
67 
68 
69 
70 void RaveKinematicVertexFitter::addTrack(const Particle* aParticlePtr)
71 {
72  rave::Vector7D raveState(aParticlePtr->getX(), aParticlePtr->getY(), aParticlePtr->getZ(), aParticlePtr->getPx(),
73  aParticlePtr->getPy(), aParticlePtr->getPz(), aParticlePtr->getMass());
74  TMatrixDSym covP = aParticlePtr->getMomentumVertexErrorMatrix();
75 
76  TMatrixDSym covE(7);
77  for (int i = 0; i < 7; i++) {
78  for (int j = 0; j < 7; j++) {
79  if (i < 3 && j < 3) covE(i, j) = covP(i + 4, j + 4);
80  if (i > 2 && j > 2) covE(i, j) = covP(i - 3, j - 3);
81  if (i < 3 && j > 2) covE(i, j) = covP(i + 4, j - 3);
82  if (i > 2 && j < 3) covE(i, j) = covP(i - 3, j + 4);
83  }
84  }
85 
86  TMatrixDSym cov = ErrorMatrixEnergyToMass(aParticlePtr->get4Vector(), covE);
87 
88  rave::Covariance7D raveCov(cov(0, 0), cov(0, 1), cov(0, 2), // x x , x y, x z
89  cov(1, 1), cov(1, 2), cov(2, 2), // y y , y z, z z
90  cov(0, 3), cov(0, 4), cov(0, 5), // x px , x py, x pz
91  cov(1, 3), cov(1, 4), cov(1, 5), // y px , y py, y pz
92  cov(2, 3), cov(2, 4), cov(2, 5), // z px , z py, z pz
93  cov(3, 3), cov(3, 4), cov(3, 5), // px px , px py, px pz
94  cov(4, 4), cov(4, 5), cov(5, 5), // py py , py pz, pz pz
95  cov(0, 6), cov(1, 6), cov(2, 6), // x m , y m, z m
96  cov(3, 6), cov(4, 6), cov(5, 6), // px m, py m, pz m
97  cov(6, 6)); // mm
98 
99  rave::TransientTrackKinematicParticle aRaveParticle(raveState, raveCov, rave::Charge(aParticlePtr->getCharge()), 1, 1);
100  // 1 and 1 are dummy values for chi2 and ndf. the are not used for the vertex fit
101 
102  m_inputParticles.push_back(aRaveParticle);
103  m_belleDaughters.push_back(const_cast<Particle*>(aParticlePtr));
104 
105 }
106 
107 void RaveKinematicVertexFitter::addMother(const Particle* aMotherParticlePtr)
108 {
109  vector<Particle*> daughters = aMotherParticlePtr->getDaughters();
110 
111  int nDaughters = daughters.size();
112  for (int i = 0; i not_eq nDaughters; ++i) {
113  addTrack(daughters[i]);
114  }
115 
116  //store input pointer so fit results can be written to the mother particle after the fit
117  auto* tmp = const_cast<Particle*>(aMotherParticlePtr);
118  m_motherParticlePtr = tmp;
119 
120 }
121 
122 
123 void RaveKinematicVertexFitter::setMother(const Particle* aMotherParticlePtr)
124 {
125  auto* tmp = const_cast<Particle*>(aMotherParticlePtr);
126  m_motherParticlePtr = tmp;
127 }
128 
129 
130 
131 int RaveKinematicVertexFitter::fit()
132 {
133  // make sure all output in this function is converted to log messages
134  auto output_capture = captureOutput();
135 
136  if (m_inputParticles.size() < 2 && m_vertFit) {
137  return -1;
138  }
139  int nOfVertices = -100;
140  if (m_useBeamSpot == true) {
141  const TVector3& bsPos = RaveSetup::getRawInstance()->m_beamSpot;
142  const TMatrixDSym& bsCov = RaveSetup::getRawInstance()->m_beamSpotCov;
143  const rave::Covariance3D bsCovRave(bsCov(0, 0), bsCov(0, 1), bsCov(0, 2), bsCov(1, 1), bsCov(1, 2), bsCov(2, 2));
144  RaveSetup::getRawInstance()->m_raveVertexFactory->setBeamSpot(rave::Ellipsoid3D(rave::Point3D(bsPos.X(), bsPos.Y(), bsPos.Z()),
145  bsCovRave));
146  }
147 
148  if (m_vertFit && m_massConstFit) {
149 
150  rave::KinematicConstraint cs2 = rave::KinematicConstraintBuilder().createMultiTrackMassKinematicConstraint(
151  m_motherParticlePtr->getPDGMass(), m_inputParticles.size());
152  try {
153  m_fittedResult = RaveSetup::getRawInstance()->m_raveKinematicTreeFactory->useVertexFitter(m_inputParticles, cs2, "", m_useBeamSpot);
154  m_fittedParticle = m_fittedResult.topParticle();
155  } catch (...) {
156  nOfVertices = 0;
157  }
158  } else {
159 
160  if (m_vertFit) {
161  if (!m_massConstFit) {
162  try {
163  m_fittedResult = RaveSetup::getRawInstance()->m_raveKinematicTreeFactory->useVertexFitter(m_inputParticles, "", m_useBeamSpot);
164  m_fittedParticle = m_fittedResult.topParticle();
165  } catch (...) {
166  nOfVertices = 0;
167  }
168  }
169 
170  }
171 
172  if (!m_vertFit && m_massConstFit) {
173 
174  try {
175 
176  Particle* aParticlePtr = m_motherParticlePtr;
177 
178  rave::Vector7D raveState(aParticlePtr->getX(), aParticlePtr->getY(), aParticlePtr->getZ(), aParticlePtr->getPx(),
179  aParticlePtr->getPy(), aParticlePtr->getPz(), aParticlePtr->getMass());
180  TMatrixFSym covP = aParticlePtr->getMomentumVertexErrorMatrix();
181 
182  float chi2 = 1;
183  float ndf = 1.;
184 
185  TMatrixDSym covE(7);
186  for (int i = 0; i < 7; i++) {
187  for (int j = 0; j < 7; j++) {
188  if (i < 3 && j < 3) covE(i, j) = covP(i + 4, j + 4);
189  if (i > 2 && j > 2) covE(i, j) = covP(i - 3, j - 3);
190  if (i < 3 && j > 2) covE(i, j) = covP(i + 4, j - 3);
191  if (i > 2 && j < 3) covE(i, j) = covP(i - 3, j + 4);
192  }
193  }
194 
195  TMatrixDSym cov = ErrorMatrixEnergyToMass(aParticlePtr->get4Vector(), covE);
196 
197  rave::Covariance7D raveCov(cov(0, 0), cov(0, 1), cov(0, 2), // x x , x y, x z
198  cov(1, 1), cov(1, 2), cov(2, 2), // y y , y z, z z
199  cov(0, 3), cov(0, 4), cov(0, 5), // x px , x py, x pz
200  cov(1, 3), cov(1, 4), cov(1, 5), // y px , y py, y pz
201  cov(2, 3), cov(2, 4), cov(2, 5), // z px , z py, z pz
202  cov(3, 3), cov(3, 4), cov(3, 5), // px px , px py, px pz
203  cov(4, 4), cov(4, 5), cov(5, 5), // py py , py pz, pz pz
204  cov(0, 6), cov(1, 6), cov(2, 6), // x m , y m, z m
205  cov(3, 6), cov(4, 6), cov(5, 6), // px m, py m, pz m
206  cov(6, 6)); // mm
207 
208  rave::TransientTrackKinematicParticle aRaveParticle(raveState, raveCov, rave::Charge(aParticlePtr->getCharge()), chi2, ndf);
209 
210  std::vector< rave::KinematicParticle > parts; parts.push_back(aRaveParticle);
211 
212  rave::KinematicConstraint constraint = rave::KinematicConstraintBuilder().createMassKinematicConstraint(
213  m_motherParticlePtr->getPDGMass(), 0.);
214 
215  if (m_motherParticlePtr->getMomentumVertexErrorMatrix().Determinant() != 0) {
216 
217  std::vector< rave::KinematicParticle > refitted = RaveSetup::getRawInstance()->m_raveKinematicTreeFactory->useParticleFitter(parts,
218  constraint, "ppf:lppf");
219 
220  m_fittedParticle = refitted[0];
221 
222  } else {
223  B2ERROR("[RaveKinematicVertexFitter]: VertexException saying ParentParticleFitter::error inverting covariance matrix occured");
224  nOfVertices = 0;
225  }
226 
227 
228  } catch (...) {
229  nOfVertices = 0;
230  }
231  }
232 
233  }
234 
235 
236  rave::Vector7D fittedState;
237  rave::Covariance7D fittedCov;
238 
239  try {
240  fittedState = m_fittedParticle.fullstate();
241  fittedCov = m_fittedParticle.fullerror();
242  } catch (...) {
243  nOfVertices = 0;
244  }
245 
246 
247  for (auto& i : m_inputParticles) i.unlink();
248  m_inputParticles.clear();
249  //
250 
251  if (nOfVertices == 0) { // vertex fit not successful
252  return 0;
253  }
254 
255  if (m_vertFit) {
256  rave::KinematicVertex fittedVertex = m_fittedResult.currentDecayVertex();
257 
258  m_fittedNdf = fittedVertex.ndf();
259  m_fittedChi2 = fittedVertex.chiSquared();
260  m_fittedPValue = ROOT::Math::chisquared_cdf_c(m_fittedChi2, m_fittedNdf);
261  m_fittedPos.SetXYZ(fittedVertex.position().x(), fittedVertex.position().y(), fittedVertex.position().z());
262 
263  } else {
264  m_fittedNdf = m_fittedParticle.ndof();
265  m_fittedChi2 = m_fittedParticle.chi2();
266  m_fittedPValue = ROOT::Math::chisquared_cdf_c(m_fittedChi2, m_fittedNdf);
267  m_fittedPos = m_motherParticlePtr->getVertex();
268  }
269 
270 
271  m_fitted4Vector.SetXYZT(fittedState.p4().p3().x(), fittedState.p4().p3().y(), fittedState.p4().p3().z(), fittedState.p4().energy());
272 
273  m_fitted7Cov.ResizeTo(7, 7);
274 
275  TMatrixDSym fitted7CovM(7);
276 
277  fitted7CovM(3, 6) = fittedCov.dpxm(); fitted7CovM(6, 3) = fitted7CovM(3, 6);
278  fitted7CovM(4, 6) = fittedCov.dpym(); fitted7CovM(6, 4) = fitted7CovM(4, 6);
279  fitted7CovM(5, 6) = fittedCov.dpzm(); fitted7CovM(6, 5) = fitted7CovM(5, 6);
280  fitted7CovM(6, 6) = fittedCov.dmm(); fitted7CovM(6, 6) = fitted7CovM(6, 6);
281  fitted7CovM(0, 6) = fittedCov.dxm(); fitted7CovM(6, 0) = fitted7CovM(0, 6);
282  fitted7CovM(1, 6) = fittedCov.dym(); fitted7CovM(6, 1) = fitted7CovM(1, 6);
283  fitted7CovM(2, 6) = fittedCov.dzm(); fitted7CovM(6, 2) = fitted7CovM(2, 6);
284 
285  fitted7CovM(3, 3) = fittedCov.dpxpx(); fitted7CovM(3, 3) = fitted7CovM(3, 3);
286  fitted7CovM(3, 4) = fittedCov.dpxpy(); fitted7CovM(4, 3) = fitted7CovM(3, 4);
287  fitted7CovM(3, 5) = fittedCov.dpxpz(); fitted7CovM(5, 3) = fitted7CovM(3, 5);
288  fitted7CovM(3, 0) = fittedCov.dxpx(); fitted7CovM(0, 3) = fitted7CovM(3, 0);
289  fitted7CovM(3, 1) = fittedCov.dypx(); fitted7CovM(1, 3) = fitted7CovM(3, 1);
290  fitted7CovM(3, 2) = fittedCov.dzpx(); fitted7CovM(2, 3) = fitted7CovM(3, 2);
291 
292  fitted7CovM(4, 4) = fittedCov.dpypy(); fitted7CovM(4, 4) = fitted7CovM(4, 4);
293  fitted7CovM(4, 5) = fittedCov.dpypz(); fitted7CovM(5, 4) = fitted7CovM(4, 5);
294  fitted7CovM(4, 0) = fittedCov.dxpy(); fitted7CovM(0, 4) = fitted7CovM(4, 0);
295  fitted7CovM(4, 1) = fittedCov.dypy(); fitted7CovM(1, 4) = fitted7CovM(4, 1);
296  fitted7CovM(4, 2) = fittedCov.dzpy(); fitted7CovM(2, 4) = fitted7CovM(4, 2);
297 
298  fitted7CovM(5, 5) = fittedCov.dpzpz(); fitted7CovM(5, 5) = fitted7CovM(5, 5);
299  fitted7CovM(5, 0) = fittedCov.dxpz(); fitted7CovM(0, 5) = fitted7CovM(5, 0);
300  fitted7CovM(5, 1) = fittedCov.dypz(); fitted7CovM(1, 5) = fitted7CovM(5, 1);
301  fitted7CovM(5, 2) = fittedCov.dzpz(); fitted7CovM(2, 5) = fitted7CovM(5, 2);
302 
303  fitted7CovM(0, 0) = fittedCov.dxx(); fitted7CovM(0, 0) = fitted7CovM(0, 0);
304  fitted7CovM(0, 1) = fittedCov.dxy(); fitted7CovM(1, 0) = fitted7CovM(0, 1);
305  fitted7CovM(0, 2) = fittedCov.dxz(); fitted7CovM(2, 0) = fitted7CovM(0, 2);
306 
307  fitted7CovM(1, 1) = fittedCov.dyy(); fitted7CovM(1, 1) = fitted7CovM(1, 1);
308  fitted7CovM(1, 2) = fittedCov.dyz(); fitted7CovM(2, 1) = fitted7CovM(1, 2);
309 
310  fitted7CovM(2, 2) = fittedCov.dzz(); fitted7CovM(2, 2) = fitted7CovM(2, 2);
311 
312  TMatrixDSym fitted7CovE = ErrorMatrixMassToEnergy(m_fitted4Vector, fitted7CovM);
313 
314  for (int i = 0; i < 7; i++) {
315  for (int j = 0; j < 7; j++) {
316  if (i < 4 && j < 4) m_fitted7Cov(i, j) = fitted7CovE(i + 3, j + 3);
317  if (i > 3 && j > 3) m_fitted7Cov(i, j) = fitted7CovE(i - 4, j - 4);
318  if (i < 4 && j > 3) m_fitted7Cov(i, j) = fitted7CovE(i + 3, j - 4);
319  if (i > 3 && j < 4) m_fitted7Cov(i, j) = fitted7CovE(i - 4, j + 3);
320  }
321  }
322 
323  return 1;
324 }
325 
326 void RaveKinematicVertexFitter::updateMother()
327 {
328  m_motherParticlePtr->updateMomentum(m_fitted4Vector, m_fittedPos, m_fitted7Cov, m_fittedPValue);
329 }
330 
331 void RaveKinematicVertexFitter::updateDaughters()
332 {
333 
334  m_fittedResult.topParticle();
335  std::vector< rave::KinematicParticle > rDau = m_fittedResult.daughterParticles();
336  std::vector<Belle2::Particle*> bDau = m_belleDaughters;
337  if (rDau.size() == bDau.size()) {
338  for (unsigned ii = 0; ii < bDau.size(); ii++) {
339  rave::Vector7D fittedState;
340  rave::Covariance7D fittedCov;
341  fittedState = rDau[ii].fullstate();
342  fittedCov = rDau[ii].fullerror();
343 
344  TLorentzVector p4;
345  p4.SetXYZT(fittedState.p4().p3().x(), fittedState.p4().p3().y(), fittedState.p4().p3().z(), fittedState.p4().energy());
346 
347  TVector3 x3(fittedState.x(), fittedState.y(), fittedState.z());
348 
349 
350  TMatrixDSym fitted7CovM(7);
351  fitted7CovM(3, 6) = fittedCov.dpxm(); fitted7CovM(6, 3) = fitted7CovM(3, 6);
352  fitted7CovM(4, 6) = fittedCov.dpym(); fitted7CovM(6, 4) = fitted7CovM(4, 6);
353  fitted7CovM(5, 6) = fittedCov.dpzm(); fitted7CovM(6, 5) = fitted7CovM(5, 6);
354  fitted7CovM(6, 6) = fittedCov.dmm(); fitted7CovM(6, 6) = fitted7CovM(6, 6);
355  fitted7CovM(0, 6) = fittedCov.dxm(); fitted7CovM(6, 0) = fitted7CovM(0, 6);
356  fitted7CovM(1, 6) = fittedCov.dym(); fitted7CovM(6, 1) = fitted7CovM(1, 6);
357  fitted7CovM(2, 6) = fittedCov.dzm(); fitted7CovM(6, 2) = fitted7CovM(2, 6);
358 
359  fitted7CovM(3, 3) = fittedCov.dpxpx(); fitted7CovM(3, 3) = fitted7CovM(3, 3);
360  fitted7CovM(3, 4) = fittedCov.dpxpy(); fitted7CovM(4, 3) = fitted7CovM(3, 4);
361  fitted7CovM(3, 5) = fittedCov.dpxpz(); fitted7CovM(5, 3) = fitted7CovM(3, 5);
362  fitted7CovM(3, 0) = fittedCov.dxpx(); fitted7CovM(0, 3) = fitted7CovM(3, 0);
363  fitted7CovM(3, 1) = fittedCov.dypx(); fitted7CovM(1, 3) = fitted7CovM(3, 1);
364  fitted7CovM(3, 2) = fittedCov.dzpx(); fitted7CovM(2, 3) = fitted7CovM(3, 2);
365 
366  fitted7CovM(4, 4) = fittedCov.dpypy(); fitted7CovM(4, 4) = fitted7CovM(4, 4);
367  fitted7CovM(4, 5) = fittedCov.dpypz(); fitted7CovM(5, 4) = fitted7CovM(4, 5);
368  fitted7CovM(4, 0) = fittedCov.dxpy(); fitted7CovM(0, 4) = fitted7CovM(4, 0);
369  fitted7CovM(4, 1) = fittedCov.dypy(); fitted7CovM(1, 4) = fitted7CovM(4, 1);
370  fitted7CovM(4, 2) = fittedCov.dzpy(); fitted7CovM(2, 4) = fitted7CovM(4, 2);
371 
372  fitted7CovM(5, 5) = fittedCov.dpzpz(); fitted7CovM(5, 5) = fitted7CovM(5, 5);
373  fitted7CovM(5, 0) = fittedCov.dxpz(); fitted7CovM(0, 5) = fitted7CovM(5, 0);
374  fitted7CovM(5, 1) = fittedCov.dypz(); fitted7CovM(1, 5) = fitted7CovM(5, 1);
375  fitted7CovM(5, 2) = fittedCov.dzpz(); fitted7CovM(2, 5) = fitted7CovM(5, 2);
376 
377  fitted7CovM(0, 0) = fittedCov.dxx(); fitted7CovM(0, 0) = fitted7CovM(0, 0);
378  fitted7CovM(0, 1) = fittedCov.dxy(); fitted7CovM(1, 0) = fitted7CovM(0, 1);
379  fitted7CovM(0, 2) = fittedCov.dxz(); fitted7CovM(2, 0) = fitted7CovM(0, 2);
380 
381  fitted7CovM(1, 1) = fittedCov.dyy(); fitted7CovM(1, 1) = fitted7CovM(1, 1);
382  fitted7CovM(1, 2) = fittedCov.dyz(); fitted7CovM(2, 1) = fitted7CovM(1, 2);
383 
384  fitted7CovM(2, 2) = fittedCov.dzz(); fitted7CovM(2, 2) = fitted7CovM(2, 2);
385 
386  TMatrixDSym fitted7CovE = ErrorMatrixMassToEnergy(m_fitted4Vector, fitted7CovM);
387 
388  TMatrixDSym fitted7CovDauM(7);
389  for (int i = 0; i < 7; i++) {
390  for (int j = 0; j < 7; j++) {
391  if (i < 4 && j < 4) fitted7CovDauM(i, j) = fitted7CovE(i + 3, j + 3);
392  if (i > 3 && j > 3) fitted7CovDauM(i, j) = fitted7CovE(i - 4, j - 4);
393  if (i < 4 && j > 3) fitted7CovDauM(i, j) = fitted7CovE(i + 3, j - 4);
394  if (i > 3 && j < 4) fitted7CovDauM(i, j) = fitted7CovE(i - 4, j + 3);
395  }
396  }
397 
398  float pValDau = rDau[ii].chi2();
399 
400  bDau[ii]->updateMomentum(p4, x3, fitted7CovDauM, pValDau);
401 
402  }
403 
404  } else B2ERROR("Error in Daughters update");
405 
406 }
407 
408 Particle* RaveKinematicVertexFitter::getMother()
409 {
410  return m_motherParticlePtr;
411 }
412 
413 TVector3 RaveKinematicVertexFitter::getPos()
414 {
415  return m_fittedPos;
416 }
417 
418 double RaveKinematicVertexFitter::getPValue()
419 {
420  return m_fittedPValue;
421 }
422 
423 double RaveKinematicVertexFitter::getNdf()
424 {
425  return m_fittedNdf;
426 }
427 
428 double RaveKinematicVertexFitter::getChi2()
429 {
430  return m_fittedChi2;
431 }
432 
433 TMatrixDSym RaveKinematicVertexFitter::getCov()
434 {
435  return m_fitted7Cov;
436 }
437 
438 TMatrixDSym RaveKinematicVertexFitter::getVertexErrorMatrix()
439 {
440  TMatrixDSym posErr(3);
441  posErr(0, 0) = m_fitted7Cov(4, 4);
442  posErr(0, 1) = m_fitted7Cov(4, 5);
443  posErr(0, 2) = m_fitted7Cov(4, 6);
444  posErr(1, 0) = m_fitted7Cov(5, 4);
445  posErr(1, 1) = m_fitted7Cov(5, 5);
446  posErr(1, 2) = m_fitted7Cov(5, 6);
447  posErr(2, 0) = m_fitted7Cov(6, 4);
448  posErr(2, 1) = m_fitted7Cov(6, 5);
449  posErr(2, 2) = m_fitted7Cov(6, 6);
450 
451  return posErr;
452 }
453 
454 
455 TMatrixDSym RaveKinematicVertexFitter::ErrorMatrixMassToEnergy(const TLorentzVector& p4, const TMatrixDSym& MassErr)
456 {
457 
458  TMatrix jac(7, 7);
459  for (int i = 0; i < 7; i++) {
460  for (int j = 0; j < 7; j++) {
461  if (i == j) jac(i, j) = 1;
462  else jac(i, j) = 0;
463  }
464  }
465  jac(6, 3) = p4.Px() / p4.E();
466  jac(6, 4) = p4.Py() / p4.E();
467  jac(6, 5) = p4.Pz() / p4.E();
468  jac(6, 6) = p4.M() / p4.E();
469 
470 
471  TMatrix jact(7, 7); jact.Transpose(jac);
472  TMatrix EnergyErrPart(7, 7); EnergyErrPart.Mult(jac, MassErr);
473  TMatrix EnergyErrTemp(7, 7); EnergyErrTemp.Mult(EnergyErrPart, jact);
474 
475  TMatrixDSym EnergyErr(7);
476  for (int i = 0; i < 7; i++) {
477  for (int j = 0; j < 7; j++) {
478  EnergyErr(i, j) = EnergyErrTemp(i, j);
479  }
480  }
481 
482  return EnergyErr;
483 }
484 
485 TMatrixDSym RaveKinematicVertexFitter::ErrorMatrixEnergyToMass(const TLorentzVector& p4, const TMatrixDSym& EnergyErr)
486 {
487 
488  TMatrix jac(7, 7);
489  for (int i = 0; i < 7; i++) {
490  for (int j = 0; j < 7; j++) {
491  if (i == j) jac(i, j) = 1;
492  else jac(i, j) = 0;
493  }
494  }
495  jac(6, 3) = -1 * p4.Px() / p4.M();
496  jac(6, 4) = -1 * p4.Py() / p4.M();
497  jac(6, 5) = -1 * p4.Pz() / p4.M();
498  jac(6, 6) = p4.E() / p4.M();
499 
500  TMatrix jact(7, 7); jact.Transpose(jac);
501  TMatrix MassErrPart(7, 7); MassErrPart.Mult(jac, EnergyErr);
502  TMatrix MassErrTemp(7, 7); MassErrTemp.Mult(MassErrPart, jact);
503 
504  TMatrixDSym MassErr(7);
505  for (int i = 0; i < 7; i++) {
506  for (int j = 0; j < 7; j++) {
507  MassErr(i, j) = MassErrTemp(i, j);
508  }
509  }
510 
511  return MassErr;
512 }
Belle2::Particle::getPx
float getPx() const
Returns x component of momentum.
Definition: Particle.h:502
Belle2::analysis::RaveKinematicVertexFitter::RaveKinematicVertexFitter
RaveKinematicVertexFitter()
The default constructor checks if RaveSetup was initialized and will set the attributes of RaveKinema...
Definition: RaveKinematicVertexFitter.cc:39
Belle2::analysis::RaveKinematicVertexFitter::m_massConstFit
bool m_massConstFit
flag determines if the mass fit is performed
Definition: RaveKinematicVertexFitter.h:167
Belle2::analysis::RaveKinematicVertexFitter::getChi2
double getChi2()
get the χ² of the fitted vertex.
Definition: RaveKinematicVertexFitter.cc:428
Belle2::analysis::RaveKinematicVertexFitter::m_fittedChi2
double m_fittedChi2
chi^2 of the vertex fit
Definition: RaveKinematicVertexFitter.h:175
Belle2::Particle::getZ
float getZ() const
Returns z component of vertex position.
Definition: Particle.h:556
Belle2::Particle::getPy
float getPy() const
Returns y component of momentum.
Definition: Particle.h:511
Belle2::Particle::getPz
float getPz() const
Returns z component of momentum.
Definition: Particle.h:520
Belle2::analysis::RaveKinematicVertexFitter::m_fittedParticle
rave::KinematicParticle m_fittedParticle
Particle fit output.
Definition: RaveKinematicVertexFitter.h:161
Belle2::analysis::RaveKinematicVertexFitter::getPValue
double getPValue()
get the p value of the fitted vertex.
Definition: RaveKinematicVertexFitter.cc:418
Belle2::Particle::getVertex
TVector3 getVertex() const
Returns vertex position (POCA for charged, IP for neutral FS particles)
Definition: Particle.h:529
Belle2::analysis::RaveKinematicVertexFitter::ErrorMatrixMassToEnergy
TMatrixDSym ErrorMatrixMassToEnergy(const TLorentzVector &p4, const TMatrixDSym &MassErr)
Convert the error matrix from P-M to P-E.
Definition: RaveKinematicVertexFitter.cc:455
Belle2::analysis::RaveKinematicVertexFitter::m_inputParticles
std::vector< rave::KinematicParticle > m_inputParticles
input particles for vertex fit in rave format
Definition: RaveKinematicVertexFitter.h:159
Belle2::analysis::RaveKinematicVertexFitter::m_fittedResult
rave::KinematicTree m_fittedResult
the output of the kinematic fit
Definition: RaveKinematicVertexFitter.h:150
Belle2::Particle::getY
float getY() const
Returns y component of vertex position.
Definition: Particle.h:547
Belle2::analysis::RaveKinematicVertexFitter::getNdf
double getNdf()
get the number of degrees of freedom (NDF) of the fitted vertex.
Definition: RaveKinematicVertexFitter.cc:423
Belle2::analysis::RaveKinematicVertexFitter::getCov
TMatrixDSym getCov()
get the covariance matrix (7x7).
Definition: RaveKinematicVertexFitter.cc:433
Belle2::analysis::RaveKinematicVertexFitter::getMother
Particle * getMother()
returns a pointer to the mother particle
Definition: RaveKinematicVertexFitter.cc:408
Belle2::analysis::RaveKinematicVertexFitter::getVertexErrorMatrix
TMatrixDSym getVertexErrorMatrix()
get the covariance matrix (3x3) of the of the fitted vertex position.
Definition: RaveKinematicVertexFitter.cc:438
Belle2::analysis::RaveKinematicVertexFitter::updateMother
void updateMother()
update the mother particle
Definition: RaveKinematicVertexFitter.cc:326
Belle2::analysis::RaveKinematicVertexFitter::m_fittedNdf
double m_fittedNdf
Ndf of the vertex fit.
Definition: RaveKinematicVertexFitter.h:171
Belle2::Particle::updateMomentum
void updateMomentum(const TLorentzVector &p4, const TVector3 &vertex, const TMatrixFSym &errMatrix, float pValue)
Sets Lorentz vector, position, 7x7 error matrix and p-value.
Definition: Particle.h:328
Belle2::analysis::RaveSetup::m_useBeamSpot
bool m_useBeamSpot
flag determines if beam spot information should be used for vertex fit.
Definition: RaveSetup.h:74
Belle2::analysis::RaveKinematicVertexFitter::captureOutput
IOIntercept::InterceptorScopeGuard< IOIntercept::OutputToLogMessages > captureOutput()
Start capturing the output of rave and divert it to log messages.
Definition: RaveKinematicVertexFitter.cc:48
Belle2::analysis::RaveKinematicVertexFitter::m_belleDaughters
std::vector< Particle * > m_belleDaughters
Belle Particle pointers input.
Definition: RaveKinematicVertexFitter.h:164
Belle2::analysis::RaveKinematicVertexFitter::addTrack
void addTrack(const Particle *aParticlePtr)
add a track (in the format of a Belle2::Particle) to set of tracks that should be fitted to a vertex
Definition: RaveKinematicVertexFitter.cc:70
Belle2::analysis::RaveKinematicVertexFitter::ErrorMatrixEnergyToMass
TMatrixDSym ErrorMatrixEnergyToMass(const TLorentzVector &p4, const TMatrixDSym &EnergyErr)
Convert the error matrix from P-E to P-M.
Definition: RaveKinematicVertexFitter.cc:485
Belle2::Particle::getCharge
float getCharge(void) const
Returns particle charge.
Definition: Particle.cc:586
Belle2::analysis::RaveSetup::m_raveVertexFactory
rave::VertexFactory * m_raveVertexFactory
The RAVE vertex factory is the principal interface offered by the RAVE vertex fitting library.
Definition: RaveSetup.h:79
Belle2::Particle::getDaughters
std::vector< Belle2::Particle * > getDaughters() const
Returns a vector of pointers to daughter particles.
Definition: Particle.cc:601
Belle2::Particle::getMass
float getMass() const
Returns invariant mass (= nominal for FS particles)
Definition: Particle.h:433
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::Particle::getMomentumVertexErrorMatrix
TMatrixFSym getMomentumVertexErrorMatrix() const
Returns 7x7 error matrix.
Definition: Particle.cc:386
Belle2::IOIntercept::start_intercept
InterceptorScopeGuard< T > start_intercept(T &interceptor)
Convenience wrapper to simplify use of InterceptorScopeGuard<T>.
Definition: IOIntercept.h:355
Belle2::analysis::RaveKinematicVertexFitter::m_vertFit
bool m_vertFit
flag determines if the vertex fit is performed
Definition: RaveKinematicVertexFitter.h:169
Belle2::analysis::RaveSetup::m_raveKinematicTreeFactory
rave::KinematicTreeFactory * m_raveKinematicTreeFactory
< The RAVE Kinematic Tree factory is the principal interface offered by the RAVE for kinematic vertex...
Definition: RaveSetup.h:83
Belle2::analysis::RaveKinematicVertexFitter::setVertFit
void setVertFit(bool isVertFit=true)
Set vertex fit: set false in case of mass fit only.
Definition: RaveKinematicVertexFitter.cc:63
Belle2::IOIntercept::InterceptorScopeGuard
Simple RAII guard for output interceptor.
Definition: IOIntercept.h:311
Belle2::analysis::RaveKinematicVertexFitter::fit
int fit()
do the kinematic vertex fit with all tracks previously added with the addTrack or addMother function.
Definition: RaveKinematicVertexFitter.cc:131
Belle2::analysis::RaveKinematicVertexFitter::m_motherParticlePtr
Particle * m_motherParticlePtr
pointer to the mother particle who's daughters will be used in the fit.
Definition: RaveKinematicVertexFitter.h:144
Belle2::analysis::RaveKinematicVertexFitter::addMother
void addMother(const Particle *aMotherParticlePtr)
All daughters of the argument of this function will be used as input for the vertex fit.
Definition: RaveKinematicVertexFitter.cc:107
Belle2::Particle::getX
float getX() const
Returns x component of vertex position.
Definition: Particle.h:538
Belle2::Particle
Class to store reconstructed particles.
Definition: Particle.h:77
Belle2::analysis::RaveKinematicVertexFitter::updateDaughters
void updateDaughters()
update the Daughters particles
Definition: RaveKinematicVertexFitter.cc:331
Belle2::analysis::RaveKinematicVertexFitter::getPos
TVector3 getPos()
get the position of the fitted vertex.
Definition: RaveKinematicVertexFitter.cc:413
Belle2::analysis::RaveKinematicVertexFitter::m_useBeamSpot
bool m_useBeamSpot
flag determines if the beam spot will be used or not.
Definition: RaveKinematicVertexFitter.h:141
Belle2::analysis::RaveKinematicVertexFitter::m_fitted4Vector
TLorentzVector m_fitted4Vector
4 momentum of the mother particle after the fit
Definition: RaveKinematicVertexFitter.h:179
Belle2::analysis::RaveSetup::getRawInstance
static RaveSetup * getRawInstance()
Same as getInstance(), but no check if the instance is initialised.
Definition: RaveSetup.cc:29
Belle2::IOIntercept::OutputToLogMessages
Capture stdout and stderr and convert into log messages.
Definition: IOIntercept.h:236
Belle2::Particle::getPDGMass
float getPDGMass(void) const
Returns uncertaint on the invariant mass (requiers valid momentum error matrix)
Definition: Particle.cc:577
Belle2::LogConfig::c_Debug
@ c_Debug
Debug: for code development.
Definition: LogConfig.h:36
Belle2::analysis::RaveKinematicVertexFitter::setMother
void setMother(const Particle *aMotherParticlePtr)
Set Mother particle for Vertex/momentum update.
Definition: RaveKinematicVertexFitter.cc:123
Belle2::analysis::RaveKinematicVertexFitter::setMassConstFit
void setMassConstFit(bool isConstFit=true)
Set mass constrained fit
Definition: RaveKinematicVertexFitter.cc:58
Belle2::analysis::RaveKinematicVertexFitter::~RaveKinematicVertexFitter
~RaveKinematicVertexFitter()
Destructor.
Belle2::Particle::get4Vector
TLorentzVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:464
Belle2::analysis::RaveSetup::m_beamSpot
TVector3 m_beamSpot
beam spot position.
Definition: RaveSetup.h:75
Belle2::Point3D
HepGeom::Point3D< double > Point3D
3D point
Definition: Cell.h:33
Belle2::analysis::RaveKinematicVertexFitter::m_fitted7Cov
TMatrixFSym m_fitted7Cov
7x7 errror matrix of the mother particle after the fit
Definition: RaveKinematicVertexFitter.h:181
Belle2::analysis::RaveKinematicVertexFitter::m_fittedPos
TVector3 m_fittedPos
Fitted vertex position.
Definition: RaveKinematicVertexFitter.h:177
Belle2::analysis::RaveSetup::m_beamSpotCov
TMatrixDSym m_beamSpotCov
beam spot position covariance matrix.
Definition: RaveSetup.h:76
Belle2::analysis::RaveKinematicVertexFitter::m_fittedPValue
double m_fittedPValue
Pvalue of the fit result.
Definition: RaveKinematicVertexFitter.h:173