Belle II Software  light-2212-foldex
ParticleVertexFitterModule.cc
1 /**************************************************************************
2  * basf2 (Belle II Analysis Software Framework) *
3  * Author: The Belle II Collaboration *
4  * *
5  * See git log for contributors and copyright holders. *
6  * This file is licensed under LGPL-3.0, see LICENSE.md. *
7  **************************************************************************/
8 
9 // Own include
10 #include <analysis/modules/ParticleVertexFitter/ParticleVertexFitterModule.h>
11 
12 // framework aux
13 #include <framework/gearbox/Unit.h>
14 #include <framework/gearbox/Const.h>
15 #include <framework/geometry/B2Vector3.h>
16 #include <framework/logging/Logger.h>
17 #include <framework/particledb/EvtGenDatabasePDG.h>
18 
19 // dataobjects
20 #include <analysis/dataobjects/Particle.h>
21 #include <analysis/dataobjects/Btube.h>
22 
23 // utilities
24 #include <analysis/utility/CLHEPToROOT.h>
25 #include <analysis/utility/PCmsLabTransform.h>
26 #include <analysis/utility/ParticleCopy.h>
27 #include <analysis/utility/ROOTToCLHEP.h>
28 
29 // Magnetic field
30 #include <framework/geometry/BFieldManager.h>
31 
32 // KFit
33 #include <analysis/VertexFitting/KFit/KFitConst.h>
34 
35 #include <TVector.h>
36 #include <TRotation.h>
37 #include <TMath.h>
38 
39 using namespace std;
40 using namespace Belle2;
41 
42 
43 //-----------------------------------------------------------------
44 // Register module
45 //-----------------------------------------------------------------
46 
47 REG_MODULE(ParticleVertexFitter);
48 
49 //-----------------------------------------------------------------
50 // Implementation
51 //-----------------------------------------------------------------
52 
53 ParticleVertexFitterModule::ParticleVertexFitterModule() : Module(),
54  m_Bfield(0)
55 {
56  // set module description (e.g. insert text)
57  setDescription("Vertex fitter for modular analysis");
59 
60  // Add parameters
61  addParam("listName", m_listName, "name of particle list", string(""));
62  addParam("confidenceLevel", m_confidenceLevel,
63  "Confidence level to accept the fit. Particle candidates with "
64  "p-value less than confidenceLevel are removed from the particle "
65  "list. If set to -1, all candidates are kept; if set to 0, "
66  "the candidates failing the fit are removed.",
67  0.001);
68  addParam("vertexFitter", m_vertexFitter, "KFit or Rave", string("KFit"));
69  addParam("fitType", m_fitType, "type of the kinematic fit (vertex, massvertex, mass)", string("vertex"));
70  addParam("withConstraint", m_withConstraint,
71  "additional constraint on vertex: ipprofile, iptube, mother, iptubecut, pointing, btube",
72  string(""));
73  addParam("decayString", m_decayString, "specifies which daughter particles are included in the kinematic fit", string(""));
74  addParam("updateDaughters", m_updateDaughters, "true: update the daughters after the vertex fit", false);
75  addParam("smearing", m_smearing, "smear IP tube width by given length", 0.002);
76  addParam("recoilMass", m_recoilMass, "recoil invariant mass (GeV)", 0.);
77  addParam("massConstraintList", m_massConstraintList,
78  "Type::[int]. List of daughter particles to mass constrain with int = pdg code. (only for MassFourCKFit)", {});
79  addParam("massConstraintListParticlename", m_massConstraintListParticlename,
80  "Type::[string]. List of daughter particles to mass constrain with string = particle name. (only for MassFourCKFit)", {});
81 }
82 
84 {
85  // Particle list with name m_listName has to exist
86  m_plist.isRequired(m_listName);
87 
88  // magnetic field
89  m_Bfield = BFieldManager::getFieldInTesla(ROOT::Math::XYZVector(0, 0, 0)).Z();
90 
91  // RAVE setup
92  if (m_vertexFitter == "Rave")
94 
95  B2DEBUG(1, "ParticleVertexFitterModule : magnetic field = " << m_Bfield);
96 
97 
98  if (m_decayString != "")
100 
101  if ((m_massConstraintList.size()) == 0 && (m_massConstraintListParticlename.size()) > 0) {
102  for (auto& containedParticle : m_massConstraintListParticlename) {
103  TParticlePDG* particletemp = TDatabasePDG::Instance()->GetParticle((containedParticle).c_str());
104  m_massConstraintList.push_back(particletemp->PdgCode());
105  }
106  }
107 
108  B2INFO("ParticleVertexFitter: Performing " << m_fitType << " fit on " << m_listName << " using " << m_vertexFitter);
109  if (m_decayString != "")
110  B2INFO("ParticleVertexFitter: Using specified decay string: " << m_decayString);
111  if (m_withConstraint != "")
112  B2INFO("ParticleVertexFitter: Additional " << m_withConstraint << " will be applied");
113 
114 }
115 
117 {
118  //TODO: set magnetic field for each run
119  //m_Bfield = BFieldMap::Instance().getBField(B2Vector3D(0,0,0)).Z();
120  //TODO: set IP spot size for each run
121 }
122 
124 {
125  if (m_vertexFitter == "Rave")
127 
128  m_BeamSpotCenter = m_beamSpotDB->getIPPosition();
129  m_beamSpotCov.ResizeTo(3, 3);
130  TMatrixDSym beamSpotCov(3);
131  if (m_withConstraint == "ipprofile") m_beamSpotCov = m_beamSpotDB->getCovVertex();
132  if (m_withConstraint == "iptube") {
133  if (m_smearing > 0 && m_vertexFitter == "KFit") {
135  } else {
137  }
138  }
139  if (m_withConstraint == "iptubecut") { // for development purpose only
140  m_BeamSpotCenter = B2Vector3D(0.001, 0., .013);
141  findConstraintBoost(0.03);
142  }
143  if ((m_vertexFitter == "Rave") && (m_withConstraint == "ipprofile" || m_withConstraint == "iptube"
144  || m_withConstraint == "mother" || m_withConstraint == "iptubecut" || m_withConstraint == "btube"))
146 
147  std::vector<unsigned int> toRemove;
148  unsigned int nParticles = m_plist->getListSize();
149 
150  for (unsigned iPart = 0; iPart < nParticles; iPart++) {
151  Particle* particle = m_plist->getParticle(iPart);
152  m_hasCovMatrix = false;
153  if (m_updateDaughters == true) {
154  if (m_decayString.empty() || m_vertexFitter == "KFit")
155  ParticleCopy::copyDaughters(particle);
156  else B2ERROR("Daughters update works only when all daughters are selected. Daughters will not be updated");
157  }
158 
159  if (m_withConstraint == "mother") {
160  m_BeamSpotCenter = B2Vector3D(particle->getVertex().x(), particle->getVertex().y(), particle->getVertex().z());
161  m_beamSpotCov = particle->getVertexErrorMatrix();
162  }
163 
164  TMatrixFSym mother_errMatrix(7);
165  mother_errMatrix = particle->getMomentumVertexErrorMatrix();
166  for (int k = 0; k < 7; k++) {
167  for (int j = 0; j < 7; j++) {
168  if (mother_errMatrix[k][j] > 0) {
169  m_hasCovMatrix = true;
170  }
171  }
172  }
173 
174  bool hasTube = true;
175  if (m_withConstraint == "btube") {
176  Btube* Ver = particle->getRelatedTo<Btube>();
177  if (!Ver) {
178  hasTube = false;
179  toRemove.push_back(particle->getArrayIndex());
180  } else {
181  m_BeamSpotCenter.SetXYZ(Ver->getTubeCenter()(0, 0), Ver->getTubeCenter()(1, 0), Ver->getTubeCenter()(2, 0));
182  m_beamSpotCov = Ver->getTubeMatrix();
183  }
184  }
185  bool ok = false;
186  if (hasTube) {
187  ok = doVertexFit(particle);
188  }
189  if (!ok)
190  particle->setPValue(-1);
191  if (particle->getPValue() < m_confidenceLevel)
192  toRemove.push_back(particle->getArrayIndex());
193 
194  }
195  m_plist->removeParticles(toRemove);
196 
197  //free memory allocated by rave. initialize() would be enough, except that we must clean things up before program end...
198  if (m_vertexFitter == "Rave")
200 }
201 
203 {
204  // steering starts here
205 
206  if (m_Bfield == 0) {
207  B2FATAL("ParticleVertexFitter: No magnetic field");
208  }
209 
210  if (m_withConstraint != "ipprofile" &&
211  m_withConstraint != "iptube" &&
212  m_withConstraint != "mother" &&
213  m_withConstraint != "iptubecut" &&
214  m_withConstraint != "pointing" &&
215  m_withConstraint != "btube" &&
216  m_withConstraint != "")
217  B2FATAL("ParticleVertexFitter: " << m_withConstraint << " ***invalid Constraint ");
218 
219  bool ok = false;
220  // fits with KFit
221  if (m_vertexFitter == "KFit") {
222 
223  if (m_decayString != "" and m_fitType != "vertex")
224  B2FATAL("ParticleVertexFitter: KFit does not support yet selection of daughters via decay string except for vertex fit!");
225 
226  // vertex fit
227  if (m_fitType == "vertex") {
228  if (m_withConstraint == "ipprofile") {
229  ok = doKVertexFit(mother, true, false);
230  } else if (m_withConstraint == "iptube") {
231  ok = doKVertexFit(mother, false, true);
232  } else {
233  ok = doKVertexFit(mother, false, false);
234  }
235  }
236 
237  // mass-constrained vertex fit
238  if (m_fitType == "massvertex") {
239  if (m_withConstraint == "ipprofile" || m_withConstraint == "iptube" || m_withConstraint == "iptubecut") {
240  B2FATAL("ParticleVertexFitter: Invalid options - mass-constrained fit using KFit does not work with iptube or ipprofile constraint.");
241  } else if (m_withConstraint == "pointing") {
242  ok = doKMassPointingVertexFit(mother);
243  } else {
244  ok = doKMassVertexFit(mother);
245  }
246  }
247 
248  // mass fit
249  if (m_fitType == "mass") {
250  if (m_withConstraint == "ipprofile" || m_withConstraint == "iptube" || m_withConstraint == "iptubecut") {
251  B2FATAL("ParticleVertexFitter: Invalid options - mass fit using KFit does not work with iptube or ipprofile constraint.");
252  } else {
253  ok = doKMassFit(mother);
254  }
255  }
256 
257  // four C fit
258  if (m_fitType == "fourC") {
259  if (m_withConstraint == "ipprofile" || m_withConstraint == "iptube" || m_withConstraint == "iptubecut") {
260  B2FATAL("ParticleVertexFitter: Invalid options - four C fit using KFit does not work with iptube or ipprofile constraint.");
261  } else {
262  ok = doKFourCFit(mother);
263  }
264  }
265 
266  // four mass C fit
267  if (m_fitType == "massfourC") {
268  if (m_withConstraint == "ipprofile" || m_withConstraint == "iptube" || m_withConstraint == "iptubecut") {
269  B2FATAL("ParticleVertexFitter: Invalid options - four C fit using KFit does not work with iptube or ipprofile constraint.");
270  } else {
271  ok = doKMassFourCFit(mother);
272  }
273  }
274 
275  // recoil mass C fit
276  if (m_fitType == "recoilmass") {
277  if (m_withConstraint == "ipprofile" || m_withConstraint == "iptube" || m_withConstraint == "iptubecut") {
278  B2FATAL("ParticleVertexFitter: Invalid options - recoil mass fit using KFit does not work with iptube or ipprofile constraint.");
279  } else {
280  ok = doKRecoilMassFit(mother);
281  }
282  }
283 
284  // invalid KFit fit type
285  if (m_fitType != "vertex"
286  && m_fitType != "massvertex"
287  && m_fitType != "mass"
288  && m_fitType != "fourC"
289  && m_fitType != "massfourC"
290  && m_fitType != "recoilmass")
291  B2FATAL("ParticleVertexFitter: " << m_fitType << " ***invalid fit type for the vertex fitter ");
292  }
293 
294  // fits using Rave
295  if (m_vertexFitter == "Rave") {
296  try {
297  ok = doRaveFit(mother);
298  } catch (const rave::CheckedFloatException&) {
299  B2ERROR("Invalid inputs (nan/inf)?");
300  ok = false;
301  }
302  }
303 
304  // invalid fitter
305  if (m_vertexFitter != "KFit" && m_vertexFitter != "Rave")
306  B2FATAL("ParticleVertexFitter: " << m_vertexFitter << " ***invalid vertex fitter ");
307 
308  if (!ok) return false;
309 
310  // steering ends here
311 
312  //if (mother->getPValue() < m_confidenceLevel) return false;
313  return true;
314 
315 }
316 
317 bool ParticleVertexFitterModule::fillFitParticles(const Particle* mother, std::vector<const Particle*>& fitChildren,
318  std::vector<const Particle*>& twoPhotonChildren)
319 {
320  if (m_decayString.empty()) {
321  // if decayString is empty, just use all primary daughters
322  for (unsigned ichild = 0; ichild < mother->getNDaughters(); ichild++) {
323  const Particle* child = mother->getDaughter(ichild);
324  // This if allows to skip the daughters, which cannot be used in the fits, particularly K_L0 from KLM.
325  // Useful for fully-inclusive particles.
326  if (mother->getProperty() == Particle::PropertyFlags::c_IsUnspecified and child->getPValue() < 0) {
327  continue;
328  }
329  fitChildren.push_back(child);
330  }
331  } else {
332  fitChildren = m_decaydescriptor.getSelectionParticles(mother);
333  }
334 
335  auto itr = fitChildren.begin();
336  while (itr != fitChildren.end()) {
337  const Particle* child = *itr;
338 
339  if (child->getPValue() < 0) {
340  B2WARNING("Daughter with PDG code " << child->getPDGCode() << " does not have a valid error matrix.");
341  return false; // error matrix not valid
342  }
343  bool isTwoPhotonParticle = false;
344  if (m_hasCovMatrix == false) {
345  if (child->getPDGCode() == Const::pi0.getPDGCode() or child->getPDGCode() == 221) { // pi0 or eta
346  if (child->getNDaughters() == 2) {
347  if (child->getDaughter(0)->getPDGCode() == Const::photon.getPDGCode()
348  && child->getDaughter(1)->getPDGCode() == Const::photon.getPDGCode()) {
349  isTwoPhotonParticle = true;
350  }
351  }
352  }
353  }
354  if (isTwoPhotonParticle) {
355  // move children from fitChildren to twoPhotonChildren
356  twoPhotonChildren.push_back(child);
357  itr = fitChildren.erase(itr);
358  } else {
359  itr++;
360  }
361  }
362 
363  return true;
364 }
365 
366 bool ParticleVertexFitterModule::fillNotFitParticles(const Particle* mother, std::vector<const Particle*>& notFitChildren,
367  const std::vector<const Particle*>& fitChildren)
368 {
369  if (fitChildren.empty())
370  B2WARNING("[ParticleVertexFitterModule::fillNotFitParticles] fitChildren is empty! Please call fillFitParticles firstly");
371  if (!notFitChildren.empty())
372  B2WARNING("[ParticleVertexFitterModule::fillNotFitParticles] notFitChildren is NOT empty!"
373  << " The function should be called only once");
374 
375  if (m_decayString.empty())
376  // if decayString is empty, just use all primary daughters
377  return true;
378 
379  std::function<bool(const Particle*)> funcCheckInFit =
380  [&funcCheckInFit, &notFitChildren, fitChildren](const Particle * part) {
381 
382  // check if the given particle in fitChildren
383  // if it is included, return true
384  if (std::find(fitChildren.begin(), fitChildren.end(), part) != fitChildren.end())
385  return true;
386 
387  // if not, firstly check if particle has children
388  if (part->getNDaughters() == 0)
389  // if it has no children (=final-state-particle), return false
390  return false;
391 
392  // here, the given particle is not in fitChildren and has children
393  bool isAnyChildrenInFit = false;
394  vector<const Particle*> notFitChildren_tmp;
395  for (unsigned ichild = 0; ichild < part->getNDaughters(); ichild++) {
396  // call funcCheckInFit recursively for all children
397  const Particle* child = part->getDaughter(ichild);
398  bool isChildrenInFit = funcCheckInFit(child);
399  isAnyChildrenInFit = isChildrenInFit or isAnyChildrenInFit;
400 
401  // if the child is not in fitChildren, fill the child in a temporary vector
402  if (!isChildrenInFit)
403  notFitChildren_tmp.push_back(child);
404  }
405 
406  // if there are a sister in fitChildren, the children in the temporary vector will be filled in notFitChildren
407  if (isAnyChildrenInFit)
408  notFitChildren.insert(notFitChildren.end(), notFitChildren_tmp.begin(), notFitChildren_tmp.end());
409 
410  // if no children in fitChildren, the given particle should be filled instead of all children.
411 
412  return isAnyChildrenInFit;
413  };
414 
415 
416  // call funcCheckInFit for all primary children
417  for (unsigned ichild = 0; ichild < mother->getNDaughters(); ichild++) {
418  const Particle* child = mother->getDaughter(ichild);
419  bool isGivenParticleOrAnyChildrenInFit = funcCheckInFit(child);
420  if (!isGivenParticleOrAnyChildrenInFit)
421  notFitChildren.push_back(child);
422  }
423 
424  return true;
425 }
426 
428  const analysis::VertexFitKFit& kv)
429 {
430  // TODO: something like setGammaError is necessary
431  // this is just workaround for the moment
432 
433  const Particle* g1Orig = preFit->getDaughter(0);
434  const Particle* g2Orig = preFit->getDaughter(1);
435  Particle g1Temp(g1Orig->get4Vector(), 22);
436  Particle g2Temp(g2Orig->get4Vector(), 22);
437 
438  TMatrixFSym g1ErrMatrix = g1Orig->getMomentumVertexErrorMatrix();
439  TMatrixFSym g2ErrMatrix = g2Orig->getMomentumVertexErrorMatrix();
440 
441  ROOT::Math::XYZVector pos(kv.getVertex().x(), kv.getVertex().y(), kv.getVertex().z());
442  CLHEP::HepSymMatrix posErrorMatrix = kv.getVertexError();
443 
444  TMatrixFSym errMatrix(3);
445  for (int i = 0; i < 3; i++)
446  for (int j = 0; j < 3; j++)
447  errMatrix(i, j) = posErrorMatrix[i][j];
448 
449  g1ErrMatrix.SetSub(4, errMatrix);
450  g2ErrMatrix.SetSub(4, errMatrix);
451 
452  g1Temp.updateMomentum(g1Orig->get4Vector(), pos, g1ErrMatrix, 1.0);
453  g2Temp.updateMomentum(g2Orig->get4Vector(), pos, g2ErrMatrix, 1.0);
454 
455  // perform the mass fit for the two-photon particle
458 
459  km.addParticle(&g1Temp);
460  km.addParticle(&g2Temp);
461 
462  km.setVertex(kv.getVertex());
464  km.setInvariantMass(preFit->getPDGMass());
465 
466  int err = km.doFit();
467  if (err != 0) {
468  return false;
469  }
470 
471  // The update of the daughters is disabled for this mass fit.
472  bool updateDaughters = m_updateDaughters;
473  m_updateDaughters = false;
474  bool ok = makeKMassMother(km, postFit);
475  m_updateDaughters = updateDaughters;
476 
477  return ok;
478 }
479 
480 bool ParticleVertexFitterModule::doKVertexFit(Particle* mother, bool ipProfileConstraint, bool ipTubeConstraint)
481 {
482  if ((mother->getNDaughters() < 2 && !ipTubeConstraint) || mother->getNDaughters() < 1) return false;
483 
484  std::vector<const Particle*> fitChildren;
485  std::vector<const Particle*> twoPhotonChildren;
486  bool validChildren = fillFitParticles(mother, fitChildren, twoPhotonChildren);
487 
488  if (!validChildren)
489  return false;
490 
491  std::vector<const Particle*> notFitChildren;
492  fillNotFitParticles(mother, notFitChildren, fitChildren);
493 
494 
495  if (twoPhotonChildren.size() > 1) {
496  B2FATAL("[ParticleVertexFitterModule::doKVertexFit] Vertex fit using KFit does not support fit with multiple particles decaying to two photons like pi0 (yet).");
497  }
498 
499  if ((fitChildren.size() < 2 && !ipTubeConstraint) || fitChildren.size() < 1) {
500  B2WARNING("[ParticleVertexFitterModule::doKVertexFit] Number of particles with valid error matrix entering the vertex fit using KFit is too low.");
501  return false;
502  }
503 
504  // Initialise the Fitter
507 
508  for (auto& child : fitChildren)
509  kv.addParticle(child);
510 
511  if (ipProfileConstraint)
512  addIPProfileToKFit(kv);
513 
514  if (ipTubeConstraint)
515  addIPTubeToKFit(kv);
516 
517  // Perform vertex fit using only the particles with valid error matrices
518  int err = kv.doFit();
519  if (err != 0)
520  return false;
521 
522  double chi2_track = getChi2TracksLBoost(kv);
523  unsigned track_count = kv.getTrackCount();
524  mother->writeExtraInfo("chiSquared_trackL", chi2_track);
525  mother->writeExtraInfo("kFit_nTracks", track_count);
526 
527  bool ok = false;
528  if (twoPhotonChildren.size() == 0)
529  // in the case daughters do not include pi0 - this is it (fit done)
530  ok = makeKVertexMother(kv, mother);
531  else if (twoPhotonChildren.size() == 1) {
532  // there is a daughter reconstructed from two photons so without position information
533  // 1. determine vertex based on all other valid daughters
534  // 2. set position and error matrix of two-photon daughter to previously determined vertex
535  // 3. redo the fit using all particles (including two-photon particle this time)
536 
537  const Particle* twoPhotonDaughter = twoPhotonChildren[0];
538  Particle fixedTwoPhotonDaughter(twoPhotonDaughter->get4Vector(), twoPhotonDaughter->getPDGCode());
539  ok = redoTwoPhotonDaughterMassFit(&fixedTwoPhotonDaughter, twoPhotonDaughter, kv);
540  if (!ok)
541  return false;
542 
543  // finally perform the fit using all daughter particles
546 
547  for (auto& child : fitChildren)
548  kv2.addParticle(child);
549 
550  kv2.addParticle(&fixedTwoPhotonDaughter);
551 
552  if (ipProfileConstraint)
553  addIPProfileToKFit(kv2);
554 
555  err = kv2.doFit();
556 
557  if (err != 0)
558  return false;
559 
560  ok = makeKVertexMother(kv2, mother);
561  }
562 
563  // update 4-vector using not-fit-particles
564  ROOT::Math::PxPyPzEVector total4Vector(mother->get4Vector());
565  for (auto& child : notFitChildren)
566  total4Vector += child->get4Vector();
567  mother->set4Vector(total4Vector);
568 
569  return ok;
570 }
571 
573 {
574  if (mother->getNDaughters() < 2) return false;
575 
576  std::vector<const Particle*> fitChildren;
577  std::vector<const Particle*> twoPhotonChildren;
578  bool validChildren = fillFitParticles(mother, fitChildren, twoPhotonChildren);
579 
580  if (!validChildren)
581  return false;
582 
583  if (twoPhotonChildren.size() > 1) {
584  B2FATAL("[ParticleVertexFitterModule::doKVertexFit] MassVertex fit using KFit does not support fit with multiple particles decaying to two photons like pi0 (yet).");
585  }
586 
587  if (fitChildren.size() < 2) {
588  B2WARNING("[ParticleVertexFitterModule::doKVertexFit] Number of particles with valid error matrix entering the vertex fit using KFit is less than 2.");
589  return false;
590  }
591 
592  bool ok = false;
593  if (twoPhotonChildren.size() == 0) {
594  // Initialise the Fitter
597 
598  for (auto child : fitChildren)
599  kmv.addParticle(child);
600 
601  kmv.setInvariantMass(mother->getPDGMass());
602  int err = kmv.doFit();
603  if (err != 0)
604  return false;
605 
606  // in the case daughters do not include particles with two photon daughters like pi0 - this is it (fit done)
607  ok = makeKMassVertexMother(kmv, mother);
608  } else if (twoPhotonChildren.size() == 1) {
609  // there is a daughter reconstructed from two photons so without position information
610  // 1. determine vertex based on all other valid daughters
611  // 2. set position and error matrix of two-photon daughter to previously determined vertex
612  // 3. redo the fit using all particles (including two-photon particle this time)
613 
616 
617  for (auto child : fitChildren)
618  kv.addParticle(child);
619 
620  // Perform vertex fit using only the particles with valid error matrices
621  int err = kv.doFit();
622  if (err != 0)
623  return false;
624 
625  const Particle* twoPhotonDaughter = twoPhotonChildren[0];
626  Particle fixedTwoPhotonDaughter(twoPhotonDaughter->get4Vector(), twoPhotonDaughter->getPDGCode());
627  ok = redoTwoPhotonDaughterMassFit(&fixedTwoPhotonDaughter, twoPhotonDaughter, kv);
628  if (!ok)
629  return false;
630 
631  // finally perform the fit using all daughter particles
634 
635  for (auto child : fitChildren)
636  kmv2.addParticle(child);
637  kmv2.addParticle(&fixedTwoPhotonDaughter);
638 
639  kmv2.setInvariantMass(mother->getPDGMass());
640  err = kmv2.doFit();
641 
642  if (err != 0)
643  return false;
644 
645  ok = makeKMassVertexMother(kmv2, mother);
646  }
647 
648  return ok;
649 
650 }
651 
653 {
654  if (!(mother->hasExtraInfo("prodVertX") && mother->hasExtraInfo("prodVertY") && mother->hasExtraInfo("prodVertZ"))) {
655  return false;
656  }
657 
658  if (mother->getNDaughters() < 2) return false;
659 
660  std::vector<const Particle*> fitChildren;
661  std::vector<const Particle*> twoPhotonChildren;
662  bool validChildren = fillFitParticles(mother, fitChildren, twoPhotonChildren);
663 
664  if (!validChildren)
665  return false;
666 
667  if (twoPhotonChildren.size() > 0) {
668  B2FATAL("[ParticleVertexFitterModule::doKMassPointingVertexFit] MassPointingVertex fit using KFit does not support fit with two-photon daughters (yet).");
669  }
670 
671  if (fitChildren.size() < 2) {
672  B2WARNING("[ParticleVertexFitterModule::doKMassPointingVertexFit] Number of particles with valid error matrix entering the vertex fit using KFit is less than 2.");
673  return false;
674  }
675 
676  bool ok = false;
677  // Initialise the Fitter
680 
681  for (auto child : fitChildren)
682  kmpv.addParticle(child);
683 
684  kmpv.setInvariantMass(mother->getPDGMass());
685  HepPoint3D productionVertex(mother->getExtraInfo("prodVertX"),
686  mother->getExtraInfo("prodVertY"),
687  mother->getExtraInfo("prodVertZ"));
688  kmpv.setProductionVertex(productionVertex);
689  int err = kmpv.doFit();
690  if (err != 0) return false;
691 
692  ok = makeKMassPointingVertexMother(kmpv, mother);
693 
694  return ok;
695 }
696 
698 {
699  if (mother->getNDaughters() < 2) return false;
700 
703 
704  for (unsigned ichild = 0; ichild < mother->getNDaughters(); ichild++) {
705  const Particle* child = mother->getDaughter(ichild);
706 
707  if (child->getPValue() < 0) return false; // error matrix not valid
708 
709  km.addParticle(child);
710  }
711 
712  // apply mass constraint
713  km.setInvariantMass(mother->getPDGMass());
714 
715  int err = km.doFit();
716 
717  if (err != 0) return false;
718 
719  bool ok = makeKMassMother(km, mother);
720 
721  return ok;
722 }
723 
725 {
726  if (mother->getNDaughters() < 2) return false;
727 
730 
731  for (unsigned ichild = 0; ichild < mother->getNDaughters(); ichild++) {
732  const Particle* child = mother->getDaughter(ichild);
733 
734  if (child->getNDaughters() > 0) {
735  bool err = addChildofParticletoKFit(kf, child);
736  if (!err) return false;
737  } else {
738  if (child->getPValue() < 0) return false; // error matrix not valid
739 
740  kf.addParticle(child);
741  }
742  }
743 
744  // apply four momentum constraint
747 
748  int err = kf.doFit();
749 
750  if (err != 0) return false;
751 
752  bool ok = makeKFourCMother(kf, mother);
753 
754  return ok;
755 }
756 
758 {
759  if (mother->getNDaughters() < 2) return false;
760 
763 
764  for (unsigned ichild = 0; ichild < mother->getNDaughters(); ichild++) {
765  const Particle* child = mother->getDaughter(ichild);
766 
767  if (child->getNDaughters() > 0) {
768  bool massconstraint = std::find(m_massConstraintList.begin(), m_massConstraintList.end(),
769  std::abs(child->getPDGCode())) != m_massConstraintList.end();
770  std::vector<unsigned> childId;
771  bool err = addChildofParticletoMassKFit(kf, child, childId);
772  if (massconstraint) kf.addMassConstraint(child->getPDGMass(), childId);
773  if (!err) return false;
774  } else {
775  if (child->getPValue() < 0) return false; // error matrix not valid
776  kf.addParticle(child);
777  }
778  }
779 
780  // apply four momentum constraint
783 
784  int err = kf.doFit();
785 
786  if (err != 0) return false;
787 
788  bool ok = makeMassKFourCMother(kf, mother);
789 
790  return ok;
791 }
792 
794 {
797 
798  for (unsigned ichild = 0; ichild < mother->getNDaughters(); ichild++) {
799  const Particle* child = mother->getDaughter(ichild);
800 
801  if (child->getPValue() < 0) return false; // error matrix not valid
802 
803  kf.addParticle(child);
804  }
805 
806  // apply four momentum constraint
809 
810  // apply recoil mass constraint
812 
813  int err = kf.doFit();
814 
815  if (err != 0) return false;
816 
817  bool ok = makeKRecoilMassMother(kf, mother);
818 
819  return ok;
820 }
821 
823  Particle* mother)
824 {
825  enum analysis::KFitError::ECode fitError;
826  fitError = kv.updateMother(mother);
827  if (fitError != analysis::KFitError::kNoError)
828  return false;
829  if (m_decayString.empty() && m_updateDaughters == true) {
830  // update daughter momenta as well
831  // the order of daughters in the *fitter is the same as in the mother Particle
832 
833  std::vector<Particle*> daughters = mother->getDaughters();
834 
835  unsigned track_count = kv.getTrackCount();
836  if (daughters.size() != track_count)
837  return false;
838 
839  for (unsigned iChild = 0; iChild < track_count; iChild++) {
840  double a = -1 * Belle2::Const::speedOfLight * 1e-4 * m_Bfield * daughters[iChild]->getCharge();
841  double dx = kv.getVertex().x() - kv.getTrackPosition(iChild).x();
842  double dy = kv.getVertex().y() - kv.getTrackPosition(iChild).y();
843 
844  ROOT::Math::PxPyPzEVector i4Vector(kv.getTrackMomentum(iChild).x() - a * dy,
845  kv.getTrackMomentum(iChild).y() + a * dx,
846  kv.getTrackMomentum(iChild).z(),
847  kv.getTrackMomentum(iChild).t());
848  daughters[iChild]->set4VectorDividingByMomentumScaling(i4Vector);
849 
850  daughters[iChild]->setVertex(
851  CLHEPToROOT::getXYZVector(kv.getTrackPosition(iChild)));
852  daughters[iChild]->setMomentumVertexErrorMatrix(
853  CLHEPToROOT::getTMatrixFSym(kv.getTrackError(iChild)));
854  }
855 
856  } else if (m_updateDaughters == true) { // if decayString is not empty
857  // first, update only the fit children
858  std::vector<const Particle*> fitChildren = m_decaydescriptor.getSelectionParticles(mother);
859 
860  unsigned track_count = kv.getTrackCount();
861  if (fitChildren.size() != track_count)
862  return false;
863 
864  for (unsigned iChild = 0; iChild < track_count; iChild++) {
865  auto daughter = const_cast<Particle*>(fitChildren[iChild]);
866 
867  double a = -1 * Belle2::Const::speedOfLight * 1e-4 * m_Bfield * daughter->getCharge();
868  double dx = kv.getVertex().x() - kv.getTrackPosition(iChild).x();
869  double dy = kv.getVertex().y() - kv.getTrackPosition(iChild).y();
870 
871  ROOT::Math::PxPyPzEVector i4Vector(kv.getTrackMomentum(iChild).x() - a * dy,
872  kv.getTrackMomentum(iChild).y() + a * dx,
873  kv.getTrackMomentum(iChild).z(),
874  kv.getTrackMomentum(iChild).t());
875  daughter->set4VectorDividingByMomentumScaling(i4Vector);
876 
877  daughter->setVertex(CLHEPToROOT::getXYZVector(kv.getTrackPosition(iChild)));
878  daughter->setMomentumVertexErrorMatrix(CLHEPToROOT::getTMatrixFSym(kv.getTrackError(iChild)));
879  }
880 
881  // then, update other particles that have a fit-child in decay
882  std::function<bool(Particle*)> funcUpdateMomentum =
883  [&funcUpdateMomentum, fitChildren](Particle * part) {
884 
885  if (part->getNDaughters() == 0) {
886  // check if part is included in fitChildren
887  if (std::find(fitChildren.begin(), fitChildren.end(), part) != fitChildren.end())
888  return true;
889  else
890  return false;
891  }
892 
893  bool includeFitChildren = false;
894 
895  // Update daughters' momentum
896  for (auto daughter : part->getDaughters())
897  includeFitChildren = funcUpdateMomentum(daughter) || includeFitChildren;
898 
899  if (includeFitChildren) {
900  // Using updated daughters, update part's momentum
901  ROOT::Math::PxPyPzEVector sum4Vector;
902  for (auto daughter : part->getDaughters())
903  sum4Vector += daughter->get4Vector();
904 
905  part->set4VectorDividingByMomentumScaling(sum4Vector);
906  }
907 
908  return includeFitChildren;
909  };
910 
911  // Update all daughters
912  for (auto daughter : mother->getDaughters())
913  funcUpdateMomentum(daughter);
914 
915  }
916 
917 
918  return true;
919 }
920 
922  Particle* mother)
923 {
924  enum analysis::KFitError::ECode fitError;
925  fitError = kmv.updateMother(mother);
926  if (fitError != analysis::KFitError::kNoError)
927  return false;
928  if (m_decayString.empty() && m_updateDaughters == true) {
929  // update daughter momenta as well
930  // the order of daughters in the *fitter is the same as in the mother Particle
931 
932  std::vector<Particle*> daughters = mother->getDaughters();
933 
934  unsigned track_count = kmv.getTrackCount();
935  if (daughters.size() != track_count)
936  return false;
937 
938  for (unsigned iChild = 0; iChild < track_count; iChild++) {
939  double a = -1 * Belle2::Const::speedOfLight * 1e-4 * m_Bfield * daughters[iChild]->getCharge();
940  double dx = kmv.getVertex().x() - kmv.getTrackPosition(iChild).x();
941  double dy = kmv.getVertex().y() - kmv.getTrackPosition(iChild).y();
942 
943  ROOT::Math::PxPyPzEVector i4Vector(kmv.getTrackMomentum(iChild).x() - a * dy,
944  kmv.getTrackMomentum(iChild).y() + a * dx,
945  kmv.getTrackMomentum(iChild).z(),
946  kmv.getTrackMomentum(iChild).t());
947  daughters[iChild]->set4VectorDividingByMomentumScaling(i4Vector);
948 
949  daughters[iChild]->setVertex(
950  CLHEPToROOT::getXYZVector(kmv.getTrackPosition(iChild)));
951  daughters[iChild]->setMomentumVertexErrorMatrix(
952  CLHEPToROOT::getTMatrixFSym(kmv.getTrackError(iChild)));
953  }
954  }
955 
956  return true;
957 }
958 
960  Particle* mother)
961 {
962  enum analysis::KFitError::ECode fitError;
963  fitError = kmpv.updateMother(mother);
964  if (fitError != analysis::KFitError::kNoError) {
965  return false;
966  }
967 
968  if (m_decayString.empty() && m_updateDaughters == true) {
969  // update daughter momenta as well
970  // the order of daughters in the *fitter is the same as in the mother Particle
971 
972  std::vector<Particle*> daughters = mother->getDaughters();
973 
974  unsigned track_count = kmpv.getTrackCount();
975  if (daughters.size() != track_count)
976  return false;
977 
978  for (unsigned iChild = 0; iChild < track_count; iChild++) {
979  double a = -1 * Belle2::Const::speedOfLight * 1e-4 * m_Bfield * daughters[iChild]->getCharge();
980  double dx = kmpv.getVertex().x() - kmpv.getTrackPosition(iChild).x();
981  double dy = kmpv.getVertex().y() - kmpv.getTrackPosition(iChild).y();
982 
983  ROOT::Math::PxPyPzEVector i4Vector(kmpv.getTrackMomentum(iChild).x() - a * dy,
984  kmpv.getTrackMomentum(iChild).y() + a * dx,
985  kmpv.getTrackMomentum(iChild).z(),
986  kmpv.getTrackMomentum(iChild).t());
987  daughters[iChild]->set4VectorDividingByMomentumScaling(i4Vector);
988 
989  daughters[iChild]->setVertex(
990  CLHEPToROOT::getXYZVector(kmpv.getTrackPosition(iChild)));
991  daughters[iChild]->setMomentumVertexErrorMatrix(
992  CLHEPToROOT::getTMatrixFSym(kmpv.getTrackError(iChild)));
993  }
994  }
995 
996  return true;
997 }
998 
999 
1001  Particle* mother)
1002 {
1003  enum analysis::KFitError::ECode fitError;
1004  fitError = km.updateMother(mother);
1005  if (fitError != analysis::KFitError::kNoError)
1006  return false;
1007  if (m_decayString.empty() && m_updateDaughters == true) {
1008  // update daughter momenta as well
1009  // the order of daughters in the *fitter is the same as in the mother Particle
1010 
1011  std::vector<Particle*> daughters = mother->getDaughters();
1012 
1013  unsigned track_count = km.getTrackCount();
1014  if (daughters.size() != track_count)
1015  return false;
1016 
1017  for (unsigned iChild = 0; iChild < track_count; iChild++) {
1018  double a = -1 * Belle2::Const::speedOfLight * 1e-4 * m_Bfield * daughters[iChild]->getCharge();
1019  double dx = km.getVertex().x() - km.getTrackPosition(iChild).x();
1020  double dy = km.getVertex().y() - km.getTrackPosition(iChild).y();
1021 
1022  ROOT::Math::PxPyPzEVector i4Vector(km.getTrackMomentum(iChild).x() - a * dy,
1023  km.getTrackMomentum(iChild).y() + a * dx,
1024  km.getTrackMomentum(iChild).z(),
1025  km.getTrackMomentum(iChild).t());
1026  daughters[iChild]->set4VectorDividingByMomentumScaling(i4Vector);
1027 
1028  daughters[iChild]->setVertex(
1029  CLHEPToROOT::getXYZVector(km.getTrackPosition(iChild)));
1030  daughters[iChild]->setMomentumVertexErrorMatrix(
1031  CLHEPToROOT::getTMatrixFSym(km.getTrackError(iChild)));
1032  }
1033  }
1034 
1035  return true;
1036 }
1037 
1038 
1039 
1041 {
1042  enum analysis::KFitError::ECode fitError;
1043  fitError = kf.updateMother(mother);
1044  if (fitError != analysis::KFitError::kNoError)
1045  return false;
1046  mother->addExtraInfo("FourCFitProb", kf.getCHIsq());
1047  mother->addExtraInfo("FourCFitChi2", kf.getNDF());
1048  if (m_decayString.empty() && m_updateDaughters == true) {
1049  // update daughter momenta as well
1050  // the order of daughters in the *fitter is the same as in the mother Particle
1051 
1052  std::vector<Particle*> daughters = mother->getDaughters();
1053 
1054  const unsigned nd = daughters.size();
1055  unsigned l = 0;
1056  std::vector<std::vector<unsigned>> pars;
1057  std::vector<Particle*> allparticles;
1058  for (unsigned ichild = 0; ichild < nd; ichild++) {
1059  const Particle* daughter = mother->getDaughter(ichild);
1060  std::vector<unsigned> pard;
1061  if (daughter->getNDaughters() > 0) {
1062  updateMapOfTrackAndDaughter(l, pars, pard, allparticles, daughter);
1063  pars.push_back(pard);
1064  allparticles.push_back(daughters[ichild]);
1065  } else {
1066  pard.push_back(l);
1067  pars.push_back(pard);
1068  allparticles.push_back(daughters[ichild]);
1069  l++;
1070  }
1071  }
1072 
1073  unsigned track_count = kf.getTrackCount();
1074  if (l != track_count)
1075  return false;
1076 
1077  for (unsigned iDaug = 0; iDaug < allparticles.size(); iDaug++) {
1078  ROOT::Math::PxPyPzEVector childMoms;
1079  ROOT::Math::XYZVector childPoss;
1080  TMatrixFSym childErrMatrixs(7);
1081  for (unsigned int iChild : pars[iDaug]) {
1082  childMoms = childMoms +
1083  CLHEPToROOT::getLorentzVector(
1084  kf.getTrackMomentum(iChild));
1085  childPoss = childPoss +
1086  CLHEPToROOT::getXYZVector(
1087  kf.getTrackPosition(iChild));
1088  TMatrixFSym childErrMatrix =
1089  CLHEPToROOT::getTMatrixFSym(kf.getTrackError(iChild));
1090  childErrMatrixs = childErrMatrixs + childErrMatrix;
1091  }
1092  allparticles[iDaug]->set4Vector(childMoms);
1093  allparticles[iDaug]->setVertex(childPoss);
1094  allparticles[iDaug]->setMomentumVertexErrorMatrix(childErrMatrixs);
1095  }
1096  }
1097 
1098  return true;
1099 }
1100 
1102 {
1103  enum analysis::KFitError::ECode fitError;
1104  fitError = kf.updateMother(mother);
1105  if (fitError != analysis::KFitError::kNoError)
1106  return false;
1107  mother->addExtraInfo("MassFourCFitProb", TMath::Prob(kf.getCHIsq(), kf.getNDF()));
1108  mother->addExtraInfo("MassFourCFitChi2", kf.getCHIsq());
1109  mother->addExtraInfo("MassFourCFitNDF", kf.getNDF());
1110  if (m_decayString.empty() && m_updateDaughters == true) {
1111  // update daughter momenta as well
1112  // the order of daughters in the *fitter is the same as in the mother Particle
1113 
1114  std::vector<Particle*> daughters = mother->getDaughters();
1115 
1116  const unsigned nd = daughters.size();
1117  unsigned l = 0;
1118  std::vector<std::vector<unsigned>> pars;
1119  std::vector<Particle*> allparticles;
1120  for (unsigned ichild = 0; ichild < nd; ichild++) {
1121  const Particle* daughter = mother->getDaughter(ichild);
1122  std::vector<unsigned> pard;
1123  if (daughter->getNDaughters() > 0) {
1124  updateMapOfTrackAndDaughter(l, pars, pard, allparticles, daughter);
1125  pars.push_back(pard);
1126  allparticles.push_back(daughters[ichild]);
1127  } else {
1128  pard.push_back(l);
1129  pars.push_back(pard);
1130  allparticles.push_back(daughters[ichild]);
1131  l++;
1132  }
1133  }
1134 
1135  unsigned track_count = kf.getTrackCount();
1136  if (l != track_count)
1137  return false;
1138 
1139  for (unsigned iDaug = 0; iDaug < allparticles.size(); iDaug++) {
1140  ROOT::Math::PxPyPzEVector childMoms;
1141  ROOT::Math::XYZVector childPoss;
1142  TMatrixFSym childErrMatrixs(7);
1143  for (unsigned int iChild : pars[iDaug]) {
1144  childMoms = childMoms +
1145  CLHEPToROOT::getLorentzVector(
1146  kf.getTrackMomentum(iChild));
1147  childPoss = childPoss +
1148  CLHEPToROOT::getXYZVector(
1149  kf.getTrackPosition(iChild));
1150  TMatrixFSym childErrMatrix =
1151  CLHEPToROOT::getTMatrixFSym(kf.getTrackError(iChild));
1152  childErrMatrixs = childErrMatrixs + childErrMatrix;
1153  }
1154  allparticles[iDaug]->set4Vector(childMoms);
1155  allparticles[iDaug]->setVertex(childPoss);
1156  allparticles[iDaug]->setMomentumVertexErrorMatrix(childErrMatrixs);
1157  }
1158  }
1159 
1160  return true;
1161 }
1162 
1164 {
1165  enum analysis::KFitError::ECode fitError;
1166  fitError = kf.updateMother(mother);
1167  if (fitError != analysis::KFitError::kNoError)
1168  return false;
1169  mother->addExtraInfo("RecoilMassFitProb", TMath::Prob(kf.getCHIsq(), kf.getNDF()));
1170  mother->addExtraInfo("RecoilMassFitChi2", kf.getCHIsq());
1171  mother->addExtraInfo("RecoilMassFitNDF", kf.getNDF());
1172  if (m_decayString.empty() && m_updateDaughters == true) {
1173  // update daughter momenta as well
1174  // the order of daughters in the *fitter is the same as in the mother Particle
1175 
1176  std::vector<Particle*> daughters = mother->getDaughters();
1177 
1178  unsigned track_count = kf.getTrackCount();
1179  if (daughters.size() != track_count)
1180  return false;
1181 
1182  for (unsigned iChild = 0; iChild < track_count; iChild++) {
1183  double a = -1 * Belle2::Const::speedOfLight * 1e-4 * m_Bfield * daughters[iChild]->getCharge();
1184  double dx = kf.getVertex().x() - kf.getTrackPosition(iChild).x();
1185  double dy = kf.getVertex().y() - kf.getTrackPosition(iChild).y();
1186 
1187  ROOT::Math::PxPyPzEVector i4Vector(kf.getTrackMomentum(iChild).x() - a * dy,
1188  kf.getTrackMomentum(iChild).y() + a * dx,
1189  kf.getTrackMomentum(iChild).z(),
1190  kf.getTrackMomentum(iChild).t());
1191  daughters[iChild]->set4VectorDividingByMomentumScaling(i4Vector);
1192 
1193  daughters[iChild]->setVertex(
1194  CLHEPToROOT::getXYZVector(kf.getTrackPosition(iChild)));
1195  daughters[iChild]->setMomentumVertexErrorMatrix(
1196  CLHEPToROOT::getTMatrixFSym(kf.getTrackError(iChild)));
1197  }
1198  }
1199 
1200  return true;
1201 }
1202 
1203 
1204 void ParticleVertexFitterModule::updateMapOfTrackAndDaughter(unsigned& l, std::vector<std::vector<unsigned>>& pars,
1205  std::vector<unsigned>& parm, std::vector<Particle*>& allparticles, const Particle* daughter)
1206 {
1207  std::vector <Belle2::Particle*> childs = daughter->getDaughters();
1208  for (unsigned ichild = 0; ichild < daughter->getNDaughters(); ichild++) {
1209  const Particle* child = daughter->getDaughter(ichild);
1210  std::vector<unsigned> pard;
1211  if (child->getNDaughters() > 0) {
1212  updateMapOfTrackAndDaughter(l, pars, pard, allparticles, child);
1213  parm.insert(parm.end(), pard.begin(), pard.end());
1214  pars.push_back(pard);
1215  allparticles.push_back(childs[ichild]);
1216  } else {
1217  pard.push_back(l);
1218  parm.push_back(l);
1219  pars.push_back(pard);
1220  allparticles.push_back(childs[ichild]);
1221  l++;
1222  }
1223  }
1224 }
1225 
1226 
1228 {
1229  if ((m_decayString.empty() ||
1230  (m_withConstraint == "" && m_fitType != "mass")) && mother->getNDaughters() < 2) return false;
1232  if (m_withConstraint == "ipprofile" || m_withConstraint == "iptube" || m_withConstraint == "mother"
1233  || m_withConstraint == "iptubecut" || m_withConstraint == "btube")
1235 
1237  if (m_fitType == "mass") rf.setVertFit(false);
1238 
1239  if (m_decayString.empty()) {
1240  rf.addMother(mother);
1241  } else {
1242  std::vector<const Particle*> tracksVertex = m_decaydescriptor.getSelectionParticles(mother);
1243  std::vector<std::string> tracksName = m_decaydescriptor.getSelectionNames();
1244 
1245  if (allSelectedDaughters(mother, tracksVertex)) {
1246  for (auto& itrack : tracksVertex) {
1247  if (itrack != mother) rf.addTrack(itrack);
1248  }
1249  rf.setMother(mother);
1250  } else {
1251 
1253  bool mothSel = false;
1254  int nTrk = 0;
1255  for (unsigned itrack = 0; itrack < tracksVertex.size(); itrack++) {
1256  if (tracksVertex[itrack] != mother) {
1257  rsf.addTrack(tracksVertex[itrack]);
1258  B2DEBUG(1, "ParticleVertexFitterModule: Adding particle " << tracksName[itrack] << " to vertex fit ");
1259  nTrk++;
1260  }
1261  if (tracksVertex[itrack] == mother) mothSel = true;
1262  }
1263 
1264 
1265  // Fit one particle constrained to originate from the beam spot
1266  bool mothIPfit = false;
1267  if (tracksVertex.size() == 1 && mothSel == true && m_withConstraint != "" && nTrk == 0) {
1268  rsf.addTrack(tracksVertex[0]);
1269  if (tracksVertex[0] != mother)
1270  B2FATAL("ParticleVertexFitterModule: FATAL Error in IP constrained mother fit");
1271  nTrk++;
1272  mothIPfit = true;
1273  }
1274 
1275 
1276  ROOT::Math::XYZVector pos;
1277  TMatrixDSym RerrMatrix(3);
1278  int nvert = 0;
1279 
1280  // one track fit is not kinematic
1281  if (nTrk == 1) {
1283  for (auto& itrack : tracksVertex) {
1284  rsg.addTrack(itrack);
1285  nvert = rsg.fit("kalman");
1286  if (nvert > 0) {
1287  pos = rsg.getPos(0);
1288  RerrMatrix = rsg.getCov(0);
1289  double prob = rsg.getPValue(0);
1290  ROOT::Math::PxPyPzEVector mom(mother->get4Vector());
1291  TMatrixDSym errMatrix(7);
1292  for (int i = 0; i < 7; i++) {
1293  for (int j = 0; j < 7; j++) {
1294  if (i > 3 && j > 3) {errMatrix[i][j] = RerrMatrix[i - 4][j - 4];}
1295  else {errMatrix[i][j] = 0;}
1296  }
1297  }
1298  if (mothIPfit) {
1299  mother->writeExtraInfo("prodVertX", pos.X());
1300  mother->writeExtraInfo("prodVertY", pos.Y());
1301  mother->writeExtraInfo("prodVertZ", pos.Z());
1302  mother->writeExtraInfo("prodVertSxx", RerrMatrix[0][0]);
1303  mother->writeExtraInfo("prodVertSxy", RerrMatrix[0][1]);
1304  mother->writeExtraInfo("prodVertSxz", RerrMatrix[0][2]);
1305  mother->writeExtraInfo("prodVertSyx", RerrMatrix[1][0]);
1306  mother->writeExtraInfo("prodVertSyy", RerrMatrix[1][1]);
1307  mother->writeExtraInfo("prodVertSyz", RerrMatrix[1][2]);
1308  mother->writeExtraInfo("prodVertSzx", RerrMatrix[2][0]);
1309  mother->writeExtraInfo("prodVertSzy", RerrMatrix[2][1]);
1310  mother->writeExtraInfo("prodVertSzz", RerrMatrix[2][2]);
1311  } else {
1312  mother->updateMomentum(mom, pos, errMatrix, prob);
1313  }
1314  return true;
1315  } else {return false;}
1316  }
1317  } else {
1318  nvert = rsf.fit();
1319  }
1320 
1321  if (nvert > 0) {
1322  pos = rsf.getPos();
1323  RerrMatrix = rsf.getVertexErrorMatrix();
1324  double prob = rsf.getPValue();
1325  ROOT::Math::PxPyPzEVector mom(mother->get4Vector());
1326  TMatrixDSym errMatrix(7);
1327  for (int i = 0; i < 7; i++) {
1328  for (int j = 0; j < 7; j++) {
1329  if (i > 3 && j > 3) {errMatrix[i][j] = RerrMatrix[i - 4][j - 4];}
1330  else {errMatrix[i][j] = 0;}
1331  }
1332  }
1333  mother->updateMomentum(mom, pos, errMatrix, prob);
1334  } else {return false;}
1335 
1336 
1337  if (mothSel && nTrk > 1) {
1338  analysis::RaveSetup::getInstance()->setBeamSpot(B2Vector3D(pos.x(), pos.y(), pos.z()), RerrMatrix);
1339  rf.addMother(mother);
1340  int nKfit = rf.fit();
1341  rf.updateMother();
1343 
1344  if (nKfit > 0) {return true;}
1345  else return false;
1346  } else return true;
1347  }
1348  }
1349 
1350  bool okFT = false;
1351  if (m_fitType == "vertex") {
1352  okFT = true;
1353  int nVert = rf.fit();
1354  rf.updateMother();
1355  if (m_decayString.empty() && m_updateDaughters == true) rf.updateDaughters();
1356  if (nVert != 1) return false;
1357  }
1358  if (m_fitType == "mass") {
1359  // add protection
1360  okFT = true;
1361  rf.setMassConstFit(true);
1362  rf.setVertFit(false);
1363  int nVert = rf.fit();
1364  rf.updateMother();
1365  if (nVert != 1) return false;
1366  };
1367  if (m_fitType == "massvertex") {
1368  okFT = true;
1369  rf.setMassConstFit(true);
1370  int nVert = rf.fit();
1371  rf.updateMother();
1372  if (m_decayString.empty() && m_updateDaughters == true) rf.updateDaughters();
1373  if (nVert != 1) return false;
1374  };
1375  if (!okFT) {
1376  B2FATAL("fitType : " << m_fitType << " ***invalid fit type ");
1377  }
1378 
1379  return true;
1380 }
1381 
1383  const std::vector<const Particle*>& tracksVertex)
1384 {
1385  bool isAll = false;
1386  if (mother->getNDaughters() == 0) return false;
1387 
1388  int nNotIncluded = mother->getNDaughters();
1389 
1390  for (unsigned i = 0; i < mother->getNDaughters(); i++) {
1391  bool dauOk = false;
1392  for (auto& vi : tracksVertex) {
1393  if (vi == mother->getDaughter(i)) {
1394  nNotIncluded = nNotIncluded - 1;
1395  dauOk = true;
1396  }
1397  }
1398  if (!dauOk) {
1399  if (allSelectedDaughters(mother->getDaughter(i), tracksVertex)) nNotIncluded--;
1400  }
1401  }
1402  if (nNotIncluded == 0) isAll = true;
1403  return isAll;
1404 }
1405 
1407 {
1408  for (unsigned ichild = 0; ichild < particle->getNDaughters(); ichild++) {
1409  const Particle* child = particle->getDaughter(ichild);
1410  if (child->getNDaughters() > 0) addChildofParticletoKFit(kf, child);
1411  else {
1412  if (child->getPValue() < 0) return false; // error matrix not valid
1413 
1414  kf.addParticle(child);
1415  }
1416  }
1417  return true;
1418 }
1419 
1421  std::vector<unsigned>& particleId)
1422 {
1423  for (unsigned ichild = 0; ichild < particle->getNDaughters(); ichild++) {
1424  const Particle* child = particle->getDaughter(ichild);
1425  if (child->getNDaughters() > 0) {
1426  bool massconstraint = std::find(m_massConstraintList.begin(), m_massConstraintList.end(),
1427  std::abs(child->getPDGCode())) != m_massConstraintList.end();
1428  std::vector<unsigned> childId;
1429  addChildofParticletoMassKFit(kf, child, childId);
1430  if (massconstraint) kf.addMassConstraint(child->getPDGMass(), childId);
1431  particleId.insert(particleId.end(), childId.begin(), childId.end());
1432  } else {
1433  if (child->getPValue() < 0) return false; // error matrix not valid
1434  kf.addParticle(child);
1435  particleId.push_back(kf.getTrackCount() - 1);
1436  }
1437  }
1438  return true;
1439 }
1440 
1442 {
1443  HepPoint3D pos(0.0, 0.0, 0.0);
1444  CLHEP::HepSymMatrix covMatrix(3, 0);
1445 
1446  for (int i = 0; i < 3; i++) {
1447  pos[i] = m_BeamSpotCenter(i);
1448  for (int j = 0; j < 3; j++) {
1449  covMatrix[i][j] = m_beamSpotCov(i, j);
1450  }
1451  }
1452 
1453  kv.setIpProfile(pos, covMatrix);
1454 }
1455 
1457 {
1458  CLHEP::HepSymMatrix err(7, 0);
1459 
1460  for (int i = 0; i < 3; i++) {
1461  for (int j = 0; j < 3; j++) {
1462  err[i + 4][j + 4] = m_beamSpotCov(i, j);
1463  }
1464  }
1465 
1466  PCmsLabTransform T;
1467  ROOT::Math::PxPyPzEVector iptube_mom = T.getBeamFourMomentum();
1468 
1469  kv.setIpTubeProfile(
1470  ROOTToCLHEP::getHepLorentzVector(iptube_mom),
1471  ROOTToCLHEP::getPoint3DFromB2Vector(m_BeamSpotCenter),
1472  err,
1473  0.);
1474 }
1475 
1477 {
1478  PCmsLabTransform T;
1479 
1480  B2Vector3D boost = T.getBoostVector();
1481  B2Vector3D boostDir = boost.Unit();
1482 
1483  TMatrixDSym beamSpotCov = m_beamSpotDB->getCovVertex();
1484  beamSpotCov(2, 2) = cut * cut;
1485  double thetab = boostDir.Theta();
1486  double phib = boostDir.Phi();
1487 
1488  double stb = TMath::Sin(thetab);
1489  double ctb = TMath::Cos(thetab);
1490  double spb = TMath::Sin(phib);
1491  double cpb = TMath::Cos(phib);
1492 
1493 
1494  TMatrix rz(3, 3); rz(2, 2) = 1;
1495  rz(0, 0) = cpb; rz(0, 1) = spb;
1496  rz(1, 0) = -1 * spb; rz(1, 1) = cpb;
1497 
1498  TMatrix ry(3, 3); ry(1, 1) = 1;
1499  ry(0, 0) = ctb; ry(0, 2) = -1 * stb;
1500  ry(2, 0) = stb; ry(2, 2) = ctb;
1501 
1502  TMatrix r(3, 3); r.Mult(rz, ry);
1503  TMatrix rt(3, 3); rt.Transpose(r);
1504 
1505  TMatrix TubePart(3, 3); TubePart.Mult(rt, beamSpotCov);
1506  TMatrix Tube(3, 3); Tube.Mult(TubePart, r);
1507 
1508  m_beamSpotCov(0, 0) = Tube(0, 0); m_beamSpotCov(0, 1) = Tube(0, 1); m_beamSpotCov(0, 2) = Tube(0, 2);
1509  m_beamSpotCov(1, 0) = Tube(1, 0); m_beamSpotCov(1, 1) = Tube(1, 1); m_beamSpotCov(1, 2) = Tube(1, 2);
1510  m_beamSpotCov(2, 0) = Tube(2, 0); m_beamSpotCov(2, 1) = Tube(2, 1); m_beamSpotCov(2, 2) = Tube(2, 2);
1511 }
1512 
1514 {
1515  TMatrixDSym beamSpotCov = m_beamSpotDB->getCovVertex();
1516  for (int i = 0; i < 3; i++)
1517  beamSpotCov(i, i) += width * width;
1518 
1519  m_beamSpotCov = beamSpotCov;
1520 }
1521 
1523 {
1524  double chi2TrackL = 0;
1525 
1526  for (int iTrack = 0; iTrack < kv.getTrackCount(); iTrack++) {
1527 
1528  analysis::KFitTrack trk_i = kv.getTrack(iTrack); // KFitTrack contains parameters before/after fit.
1529 
1530  TMatrixFSym err = CLHEPToROOT::getTMatrixFSym(trk_i.getError(analysis::KFitConst::kBeforeFit)); // px, py, pz, E, x, y, z
1531 
1532  B2Vector3D x_before = CLHEPToROOT::getXYZVector(trk_i.getPosition(analysis::KFitConst::kBeforeFit));
1533  B2Vector3D x_after = CLHEPToROOT::getXYZVector(trk_i.getPosition());
1534  B2Vector3D dPos = x_after - x_before;
1535 
1536  PCmsLabTransform T;
1537  B2Vector3D boost3 = T.getBoostVector().Unit();
1538  TVectorD boostD(0, 6, 0., 0., 0., 0., boost3.X(), boost3.Y(), boost3.Z(), "END");
1539 
1540  double dLBoost = dPos.Dot(boost3);
1541 
1542  chi2TrackL += TMath::Power(dLBoost, 2) / err.Similarity(boostD);
1543  }
1544  return chi2TrackL;
1545 }
DataType Phi() const
The azimuth angle.
Definition: B2Vector3.h:151
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:435
DataType Theta() const
The polar angle.
Definition: B2Vector3.h:153
DataType X() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:431
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:433
B2Vector3< DataType > Unit() const
Unit vector parallel to this.
Definition: B2Vector3.h:269
DataType Dot(const B2Vector3< DataType > &p) const
Scalar product.
Definition: B2Vector3.h:290
void SetXYZ(DataType x, DataType y, DataType z)
set all coordinates using data type
Definition: B2Vector3.h:464
static ROOT::Math::XYZVector getFieldInTesla(const ROOT::Math::XYZVector &pos)
return the magnetic field at a given position in Tesla.
Definition: BFieldManager.h:61
For each MCParticle with hits in the CDC, this class stores some summarising information on those hit...
Definition: Btube.h:27
TMatrixFSym getTubeMatrix() const
Returns Btube matrix.
Definition: Btube.h:78
Eigen::Matrix< double, 3, 1 > getTubeCenter() const
Returns Btube center.
Definition: Btube.h:64
int getPDGCode() const
PDG code.
Definition: Const.h:464
static const ParticleType pi0
neutral pion particle
Definition: Const.h:665
static const double speedOfLight
[cm/ns]
Definition: Const.h:685
static const ParticleType photon
photon particle
Definition: Const.h:664
bool init(const std::string &str)
Initialise the DecayDescriptor from given string.
std::vector< std::string > getSelectionNames()
Return list of human readable names of selected particles.
std::vector< const Particle * > getSelectionParticles(const Particle *particle)
Get a vector of pointers with selected daughters in the decay tree.
Base class for Modules.
Definition: Module.h:72
void setDescription(const std::string &description)
Sets the description of the module.
Definition: Module.cc:214
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition: Module.cc:208
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition: Module.h:80
Class to hold Lorentz transformations from/to CMS and boost vector.
ROOT::Math::PxPyPzEVector getBeamFourMomentum() const
Returns LAB four-momentum of e+e-.
B2Vector3D getBoostVector() const
Returns boost vector (beta=p/E)
TMatrixDSym m_beamSpotCov
Beam spot covariance matrix.
bool makeKMassMother(analysis::MassFitKFit &kv, Particle *p)
Update mother particle after mass fit using KFit.
std::string m_withConstraint
additional constraint on vertex
bool doKMassPointingVertexFit(Particle *p)
Mass-constrained vertex fit with additional pointing constraint using KFit.
void smearBeamSpot(double width)
smear beam spot covariance
bool makeMassKFourCMother(analysis::MassFourCFitKFit &kv, Particle *p)
Update mother particle after MassFourC fit using KFit.
bool doKVertexFit(Particle *p, bool ipProfileConstraint, bool ipTubeConstraint)
Unconstrained vertex fit using KFit.
bool doKMassVertexFit(Particle *p)
Mass-constrained vertex fit using KFit.
virtual void initialize() override
Initialize the Module.
void addIPTubeToKFit(analysis::VertexFitKFit &kv)
Adds IPTube constraint to the vertex fit using KFit.
bool addChildofParticletoMassKFit(analysis::MassFourCFitKFit &kf, const Particle *particle, std::vector< unsigned > &particleId)
Adds given particle's child to the MassFourCFitKFit.
virtual void event() override
Event processor.
bool makeKMassPointingVertexMother(analysis::MassPointingVertexFitKFit &kv, Particle *p)
Update mother particle after mass-constrained vertex fit with additional pointing constraint using KF...
bool m_updateDaughters
flag for daughters update
std::string m_decayString
daughter particles selection
std::vector< std::string > m_massConstraintListParticlename
Name of the particles to be mass constraint (massfourC)
std::string m_listName
particle list name
std::vector< int > m_massConstraintList
PDG codes of the particles to be mass constraint (massfourC)
bool m_hasCovMatrix
flag for mother covariance matrix (PseudoFitter)
bool makeKVertexMother(analysis::VertexFitKFit &kv, Particle *p)
Update mother particle after unconstrained vertex fit using KFit.
bool redoTwoPhotonDaughterMassFit(Particle *postFit, const Particle *preFit, const analysis::VertexFitKFit &kv)
Combines preFit particle and vertex information from vertex fit kv to create new postFit particle.
bool makeKMassVertexMother(analysis::MassVertexFitKFit &kv, Particle *p)
Update mother particle after mass-constrained vertex fit using KFit.
bool makeKFourCMother(analysis::FourCFitKFit &kv, Particle *p)
Update mother particle after FourC fit using KFit.
bool fillFitParticles(const Particle *mother, std::vector< const Particle * > &fitChildren, std::vector< const Particle * > &twoPhotonChildren)
Fills valid particle's children (with valid error matrix) in the vector of Particles that will enter ...
bool doKFourCFit(Particle *p)
FourC fit using KFit.
B2Vector3D m_BeamSpotCenter
Beam spot position.
virtual void beginRun() override
Called when entering a new run.
bool doKRecoilMassFit(Particle *p)
RecoilMass fit using KFit.
void updateMapOfTrackAndDaughter(unsigned &l, std::vector< std::vector< unsigned >> &pars, std::vector< unsigned > &pard, std::vector< Particle * > &allparticles, const Particle *daughter)
update the map of daughter and tracks, find out which tracks belong to each daughter.
DBObjPtr< BeamSpot > m_beamSpotDB
Beam spot database object.
bool doVertexFit(Particle *p)
Main steering routine.
std::string m_vertexFitter
Vertex Fitter name.
bool doRaveFit(Particle *mother)
Fit using Rave.
bool makeKRecoilMassMother(analysis::RecoilMassKFit &kf, Particle *p)
Update mother particle after RecoilMass fit using KFit.
bool doKMassFit(Particle *p)
Mass fit using KFit.
double m_recoilMass
recoil mass for constraint
double m_confidenceLevel
required fit confidence level
bool doKMassFourCFit(Particle *p)
MassFourC fit using KFit.
DecayDescriptor m_decaydescriptor
Decay descriptor of decays to look for.
bool fillNotFitParticles(const Particle *mother, std::vector< const Particle * > &notFitChildren, const std::vector< const Particle * > &fitChildren)
Fills valid particle's children (with valid error matrix) in the vector of Particles that will not en...
void findConstraintBoost(double cut)
calculate iptube constraint (quasi cylinder along boost direction) for RAVE fit
double m_Bfield
magnetic field from data base
bool addChildofParticletoKFit(analysis::FourCFitKFit &kv, const Particle *particle)
Adds given particle's child to the FourCFitKFit.
double getChi2TracksLBoost(const analysis::VertexFitKFit &kv)
calculate the chi2 using only lboost information of tracks
double m_smearing
smearing width applied to IP tube
std::string m_fitType
type of the kinematic fit
void addIPProfileToKFit(analysis::VertexFitKFit &kv)
Adds IPProfile constraint to the vertex fit using KFit.
StoreObjPtr< ParticleList > m_plist
particle list
bool allSelectedDaughters(const Particle *mother, const std::vector< const Particle * > &tracksVertex)
check if all the Daughters (o grand-daughters) are selected for the vertex fit
Class to store reconstructed particles.
Definition: Particle.h:74
TMatrixFSym getVertexErrorMatrix() const
Returns the 3x3 position error sub-matrix.
Definition: Particle.cc:479
void writeExtraInfo(const std::string &name, const double value)
Sets the user defined extraInfo.
Definition: Particle.cc:1335
double getPValue() const
Returns chi^2 probability of fit if done or -1.
Definition: Particle.h:640
ROOT::Math::XYZVector getVertex() const
Returns vertex position (POCA for charged, IP for neutral FS particles)
Definition: Particle.h:604
bool hasExtraInfo(const std::string &name) const
Return whether the extra info with the given name is set.
Definition: Particle.cc:1293
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:441
int getProperty() const
Returns particle property as a bit pattern The values are defined in the PropertyFlags enum and descr...
Definition: Particle.h:485
unsigned getNDaughters(void) const
Returns number of daughter particles.
Definition: Particle.h:700
std::vector< Belle2::Particle * > getDaughters() const
Returns a vector of pointers to daughter particles.
Definition: Particle.cc:669
double getPDGMass(void) const
Returns uncertainty on the invariant mass (requires valid momentum error matrix)
Definition: Particle.cc:636
ROOT::Math::PxPyPzEVector get4Vector() const
Returns Lorentz vector.
Definition: Particle.h:532
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 set4Vector(const ROOT::Math::PxPyPzEVector &p4)
Sets Lorentz vector.
Definition: Particle.h:286
void setPValue(double pValue)
Sets chi^2 probability of fit.
Definition: Particle.h:353
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
TMatrixFSym getMomentumVertexErrorMatrix() const
Returns 7x7 error matrix.
Definition: Particle.cc:452
const Particle * getDaughter(unsigned i) const
Returns a pointer to the i-th daughter particle.
Definition: Particle.cc:663
double getExtraInfo(const std::string &name) const
Return given value if set.
Definition: Particle.cc:1316
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
TO * getRelatedTo(const std::string &name="", const std::string &namedRelation="") const
Get the object to which this object has a relation.
FourCFitKFit is a derived class from KFitBase to perform 4 momentum-constraint kinematical fit.
Definition: FourCFitKFit.h:30
enum KFitError::ECode setFourMomentum(const ROOT::Math::PxPyPzEVector &m)
Set an 4 Momentum for the FourC-constraint fit.
Definition: FourCFitKFit.cc:76
double getCHIsq(void) const override
Get a chi-square of the fit.
enum KFitError::ECode updateMother(Particle *mother)
Update mother particle.
enum KFitError::ECode doFit(void)
Perform a four momentum-constraint fit.
const CLHEP::HepSymMatrix getTrackError(const int id) const
Get an error matrix of the track.
Definition: KFitBase.cc:168
const CLHEP::HepLorentzVector getTrackMomentum(const int id) const
Get a Lorentz vector of the track.
Definition: KFitBase.cc:154
const HepPoint3D getTrackPosition(const int id) const
Get a position of the track.
Definition: KFitBase.cc:161
virtual int getNDF(void) const
Get an NDF of the fit.
Definition: KFitBase.cc:114
enum KFitError::ECode setMagneticField(const double mf)
Change a magnetic field from the default value KFitConst::kDefaultMagneticField.
Definition: KFitBase.cc:93
const KFitTrack getTrack(const int id) const
Get a specified track object.
Definition: KFitBase.cc:175
enum KFitError::ECode addParticle(const Particle *particle)
Add a particle to the fitter.
Definition: KFitBase.cc:59
int getTrackCount(void) const
Get the number of added tracks.
Definition: KFitBase.cc:107
ECode
ECode is a error code enumerate.
Definition: KFitError.h:34
KFitTrack is a container of the track information (Lorentz vector, position, and error matrix),...
Definition: KFitTrack.h:38
const CLHEP::HepSymMatrix getError(const int flag=KFitConst::kAfterFit) const
Get an error matrix of the track.
Definition: KFitTrack.cc:172
const HepPoint3D getPosition(const int flag=KFitConst::kAfterFit) const
Get a position of the track.
Definition: KFitTrack.cc:164
MassFitKFit is a derived class from KFitBase to perform mass-constraint kinematical fit.
Definition: MassFitKFit.h:33
enum KFitError::ECode setVertex(const HepPoint3D &v)
Set an initial vertex position for the mass-constraint fit.
Definition: MassFitKFit.cc:43
enum KFitError::ECode updateMother(Particle *mother)
Update mother particle.
Definition: MassFitKFit.cc:691
enum KFitError::ECode doFit(void)
Perform a mass-constraint fit.
Definition: MassFitKFit.cc:279
enum KFitError::ECode setInvariantMass(const double m)
Set an invariant mass for the mass-constraint fit.
Definition: MassFitKFit.cc:67
enum KFitError::ECode setVertexError(const CLHEP::HepSymMatrix &e)
Set an initial vertex error matrix for the mass-constraint fit.
Definition: MassFitKFit.cc:51
const HepPoint3D getVertex(const int flag=KFitConst::kAfterFit) const
Get a vertex position.
Definition: MassFitKFit.cc:136
MassFourCFitKFit is a derived class from KFitBase to perform mass and 4 momentum-constraint kinematic...
enum KFitError::ECode setFourMomentum(const ROOT::Math::PxPyPzEVector &m)
Set an 4 Momentum for the mass-four-constraint fit.
enum KFitError::ECode addMassConstraint(const double m, std::vector< unsigned > &childTrackId)
Set an invariant mass of daughter particle for the mass-four-momentum-constraint fit.
double getCHIsq(void) const override
Get a chi-square of the fit.
enum KFitError::ECode updateMother(Particle *mother)
Update mother particle.
enum KFitError::ECode doFit(void)
Perform a mass-four-momentum-constraint fit.
MassPointingVertexFitKFit is a derived class from KFitBase It performs a kinematical fit with three c...
enum KFitError::ECode updateMother(Particle *mother)
Update mother particle.
enum KFitError::ECode doFit(void)
Perform a mass-vertex-pointing constraint fit.
enum KFitError::ECode setInvariantMass(const double m)
Set an invariant mass for the mass-vertex-pointing constraint fit.
enum KFitError::ECode setProductionVertex(const HepPoint3D &v)
Set the production vertex of the particle.
const HepPoint3D getVertex(const int flag=KFitConst::kAfterFit) const
Get a vertex position.
MassVertexFitKFit is a derived class from KFitBase to perform mass-vertex-constraint kinematical fit.
enum KFitError::ECode updateMother(Particle *mother)
Update mother particle.
enum KFitError::ECode doFit(void)
Perform a mass-vertex-constraint fit.
enum KFitError::ECode setInvariantMass(const double m)
Set an invariant mass for the mass-vertex constraint fit.
const HepPoint3D getVertex(const int flag=KFitConst::kAfterFit) const
Get a vertex position.
The RaveKinematicVertexFitter class is part of the RaveInterface together with RaveSetup.
void addMother(const Particle *aMotherParticlePtr)
All daughters of the argument of this function will be used as input for the vertex fit.
ROOT::Math::XYZVector getPos()
get the position of the fitted vertex.
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
void setMother(const Particle *aMotherParticlePtr)
Set Mother particle for Vertex/momentum update.
int fit()
do the kinematic vertex fit with all tracks previously added with the addTrack or addMother function.
void setVertFit(bool isVertFit=true)
Set vertex fit: set false in case of mass fit only.
void setMassConstFit(bool isConstFit=true)
Set mass constrained fit
double getPValue()
get the p value of the fitted vertex.
TMatrixDSym getVertexErrorMatrix()
get the covariance matrix (3x3) of the of the fitted vertex position.
void updateDaughters()
update the Daughters particles
void unsetBeamSpot()
unset beam spot constraint
Definition: RaveSetup.cc:80
static void initialize(int verbosity=1, double MagneticField=1.5)
Set everything up so everything needed for vertex fitting is there.
Definition: RaveSetup.cc:33
static RaveSetup * getInstance()
get the pointer to the instance to get/set any of options stored in RaveSetup
Definition: RaveSetup.h:43
void setBeamSpot(const B2Vector3D &beamSpot, const TMatrixDSym &beamSpotCov)
The beam spot position and covariance is known you can set it here so that and a vertex in the beam s...
Definition: RaveSetup.cc:72
void reset()
frees memory allocated by initialize().
Definition: RaveSetup.cc:58
The RaveVertexFitter class is part of the RaveInterface together with RaveSetup.
TMatrixDSym getCov(VecSize vertexId=0) const
get the covariance matrix (3x3) of the of the fitted vertex position.
int fit(std::string options="default")
do the vertex fit with all tracks previously added with the addTrack or addMother function.
void addTrack(const Particle *const aParticlePtr)
add a track (in the format of a Belle2::Particle) to set of tracks that should be fitted to a vertex
B2Vector3D getPos(VecSize vertexId=0) const
get the position of the fitted vertex.
double getPValue(VecSize vertexId=0) const
get the p value of the fitted vertex.
RecoilMassKFit is a derived class from KFitBase to perform a kinematical fit with a recoil mass const...
enum KFitError::ECode setFourMomentum(const ROOT::Math::PxPyPzEVector &m)
Set a recoil mass .
double getCHIsq(void) const override
Get a chi-square of the fit.
enum KFitError::ECode updateMother(Particle *mother)
Update mother particle.
enum KFitError::ECode setRecoilMass(const double m)
Set an invariant mass for the four momentum-constraint fit.
enum KFitError::ECode doFit(void)
Perform a recoil-mass constraint fit.
const HepPoint3D getVertex(const int flag=KFitConst::kAfterFit) const
Get a vertex position.
VertexFitKFit is a derived class from KFitBase to perform vertex-constraint kinematical fit.
Definition: VertexFitKFit.h:34
enum KFitError::ECode setIpTubeProfile(const CLHEP::HepLorentzVector &p, const HepPoint3D &x, const CLHEP::HepSymMatrix &e, const double q)
Set a virtual IP-tube track for the vertex constraint fit.
enum KFitError::ECode updateMother(Particle *mother)
Update mother particle.
const CLHEP::HepSymMatrix getVertexError(void) const
Get a fitted vertex error matrix.
enum KFitError::ECode doFit(void)
Perform a vertex-constraint fit.
enum KFitError::ECode setIpProfile(const HepPoint3D &ip, const CLHEP::HepSymMatrix &ipe)
Set an IP-ellipsoid shape for the vertex constraint fit.
const HepPoint3D getVertex(const int flag=KFitConst::kAfterFit) const
Get a vertex position.
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition: Module.h:560
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition: Module.h:650
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
void copyDaughters(Particle *mother)
Function copies all (grand-)^n-daughter particles of the argument mother Particle.
Definition: ParticleCopy.cc:55
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:23
static const int kBeforeFit
Input parameter to specify before-fit when setting/getting a track attribute.
Definition: KFitConst.h:35