13 #include <TDatabasePDG.h>
15 #include <analysis/VertexFitting/TreeFitter/ParticleBase.h>
16 #include <analysis/VertexFitting/TreeFitter/InternalParticle.h>
17 #include <analysis/VertexFitting/TreeFitter/RecoComposite.h>
18 #include <analysis/VertexFitting/TreeFitter/RecoResonance.h>
19 #include <analysis/VertexFitting/TreeFitter/RecoTrack.h>
20 #include <analysis/VertexFitting/TreeFitter/FeedthroughParticle.h>
21 #include <analysis/VertexFitting/TreeFitter/InternalTrack.h>
23 #include <analysis/VertexFitting/TreeFitter/RecoPhoton.h>
24 #include <analysis/VertexFitting/TreeFitter/RecoKlong.h>
25 #include <analysis/VertexFitting/TreeFitter/Resonance.h>
26 #include <analysis/VertexFitting/TreeFitter/Origin.h>
27 #include <analysis/VertexFitting/TreeFitter/FitParams.h>
29 #include <framework/geometry/BFieldManager.h>
31 namespace TreeFitter {
36 m_isStronglyDecayingResonance(false),
39 m_pdgMass(particle->getPDGMass()),
41 m_pdgLifeTime(TDatabasePDG::Instance()->GetParticle(particle->getPDGCode())->Lifetime() * 1e9),
46 m_isStronglyDecayingResonance = isAResonance(particle);
47 const int pdgcode = particle->getPDGCode();
50 double fltcharge = particle->getCharge();
53 m_charge = fltcharge < 0 ? int(fltcharge - 0.5) : int(fltcharge + 0.5);
54 m_name = particle->getName();
56 m_charge = particle->getCharge() > 0 ? 1 : (particle->getCharge() < 0 ? -1 : 0);
64 m_isStronglyDecayingResonance(false),
67 m_pdgMass(particle->getPDGMass()),
69 m_pdgLifeTime(TDatabasePDG::Instance()->GetParticle(particle->getPDGCode())->Lifetime() * 1e9),
74 m_isStronglyDecayingResonance = isAResonance(particle);
75 const int pdgcode = particle->getPDGCode();
78 double fltcharge = particle->getCharge();
81 m_charge = fltcharge < 0 ? int(fltcharge - 0.5) : int(fltcharge + 0.5);
82 m_name = particle->getName();
84 m_charge = particle->getCharge() > 0 ? 1 : (particle->getCharge() < 0 ? -1 : 0);
92 m_isStronglyDecayingResonance(false),
126 B2ERROR(
"Cannot remove particle, because not found ...");
133 daughter->updateIndex(offset);
145 return new Origin(daughter, config, forceFitAll);
165 rc =
new FeedthroughParticle(
particle,
mother, config, forceFitAll);
167 rc =
new InternalTrack(
particle,
mother, config, forceFitAll,
true,
true);
220 lifetime = TDatabasePDG::Instance()->GetParticle(pdgcode)->Lifetime() * 1e9;
252 particles.push_back(
this);
256 daughter->collectVertexDaughters(particles, posindex);
267 for (
int i = 0; i < 3; ++i) {
268 fitparams.getCovariance()(posindex + i, posindex + i) = 1;
275 for (
int i = 0; i < maxrow; ++i) {
276 fitparams.getCovariance()(momindex + i, momindex + i) = 10.;
282 fitparams.getCovariance()(tauindex, tauindex) = 1.;
290 while (rc && rc->type() == kFeedthroughParticle) {
298 std::string rc =
name();
300 case 0: rc +=
"_x ";
break;
301 case 1: rc +=
"_y ";
break;
302 case 2: rc +=
"_z ";
break;
303 case 3: rc +=
"_tau";
break;
304 case 4: rc +=
"_px ";
break;
305 case 5: rc +=
"_py ";
break;
306 case 6: rc +=
"_pz ";
break;
307 case 7: rc +=
"_E ";
break;
327 map.push_back(std::pair<const ParticleBase*, int>(
this,
index()));
329 daughter->retrieveIndexMap(map);
337 rc += daughter->chiSquare(fitparams);
346 rc += daughter->nFinalChargedCandidates();
362 Eigen::Matrix < double, 1, -1, 1, 1, 3 > x_vec = fitparams.
getStateVector().segment(posindex,
dim);
363 Eigen::Matrix < double, 1, -1, 1, 1, 3 > x_m = fitparams.
getStateVector().segment(posindexmother,
dim);
364 Eigen::Matrix < double, 1, -1, 1, 1, 3 > p_vec = fitparams.
getStateVector().segment(momindex,
dim);
365 const double mom = p_vec.norm();
366 const double mom3 = mom * mom * mom;
371 p.getH()(0, momindex) = tau * (p_vec(1) * p_vec(1) + p_vec(2) * p_vec(2)) / mom3 ;
372 p.getH()(1, momindex + 1) = tau * (p_vec(0) * p_vec(0) + p_vec(2) * p_vec(2)) / mom3 ;
373 p.getH()(2, momindex + 2) = tau * (p_vec(0) * p_vec(0) + p_vec(1) * p_vec(1)) / mom3 ;
376 p.getH()(0, momindex + 1) = - tau * p_vec(0) * p_vec(1) / mom3 ;
377 p.getH()(0, momindex + 2) = - tau * p_vec(0) * p_vec(2) / mom3 ;
379 p.getH()(1, momindex + 0) = - tau * p_vec(1) * p_vec(0) / mom3 ;
380 p.getH()(1, momindex + 2) = - tau * p_vec(1) * p_vec(2) / mom3 ;
382 p.getH()(2, momindex + 0) = - tau * p_vec(2) * p_vec(0) / mom3 ;
383 p.getH()(2, momindex + 1) = - tau * p_vec(2) * p_vec(1) / mom3 ;
385 }
else if (2 ==
dim) {
388 p.getH()(0, momindex) = tau * (p_vec(1) * p_vec(1)) / mom3 ;
389 p.getH()(1, momindex + 1) = tau * (p_vec(0) * p_vec(0)) / mom3 ;
392 p.getH()(0, momindex + 1) = - tau * p_vec(0) * p_vec(1) / mom3 ;
393 p.getH()(1, momindex + 0) = - tau * p_vec(1) * p_vec(0) / mom3 ;
395 B2FATAL(
"Dimension of Geometric constraint is not 2 or 3. This will crash many things. You should feel bad.");
398 for (
int row = 0; row <
dim; ++row) {
400 double posxmother = x_m(row);
401 double posx = x_vec(row);
402 double momx = p_vec(row);
407 p.getResiduals()(row) = posxmother + tau * momx / mom - posx ;
408 p.getH()(row, posindexmother + row) = 1;
409 p.getH()(row, posindex + row) = -1;
410 p.getH()(row, tauindex) = momx / mom;
413 return ErrCode(ErrCode::Status::success);
416 void inline setExtraInfo(
Belle2::Particle* part,
const std::string& name,
const double value)
419 if (part->hasExtraInfo(name)) {
420 part->setExtraInfo(name, value);
422 part->addExtraInfo(name, value);
431 const double mass2 = mass * mass;
439 const int momindex = daughter->momIndex();
441 const double px_daughter = fitparams.getStateVector()(momindex);
442 const double py_daughter = fitparams.getStateVector()(momindex + 1);
443 const double pz_daughter = fitparams.getStateVector()(momindex + 2);
448 if (daughter->hasEnergy()) {
449 E += fitparams.getStateVector()(momindex + 3);
452 const double m = daughter->pdgMass();
453 E += std::sqrt(m * m + px_daughter * px_daughter + py_daughter * py_daughter + pz_daughter * pz_daughter);
460 p.getResiduals()(0) = mass2 - E * E + px * px + py * py + pz * pz;
464 const int momindex = daughter->momIndex();
465 p.getH()(0, momindex) = 2.0 * px;
466 p.getH()(0, momindex + 1) = 2.0 * py;
467 p.getH()(0, momindex + 2) = 2.0 * pz;
469 if (daughter->hasEnergy()) {
470 p.getH()(0, momindex + 3) = -2.0 * E;
472 const double px_daughter = fitparams.getStateVector()(momindex);
473 const double py_daughter = fitparams.getStateVector()(momindex + 1);
474 const double pz_daughter = fitparams.getStateVector()(momindex + 2);
475 const double m = daughter->pdgMass();
477 const double E_daughter = std::sqrt(m * m + px_daughter * px_daughter + py_daughter * py_daughter + pz_daughter * pz_daughter);
478 const double E_by_E_daughter = E / E_daughter;
479 p.getH()(0, momindex) -= 2.0 * E_by_E_daughter * px_daughter;
480 p.getH()(0, momindex + 1) -= 2.0 * E_by_E_daughter * py_daughter;
481 p.getH()(0, momindex + 2) -= 2.0 * E_by_E_daughter * pz_daughter;
485 return ErrCode(ErrCode::Status::success);
492 const double mass2 = mass * mass;
494 const double px = fitparams.getStateVector()(momindex);
495 const double py = fitparams.getStateVector()(momindex + 1);
496 const double pz = fitparams.getStateVector()(momindex + 2);
497 const double E = fitparams.getStateVector()(momindex + 3);
502 p.getResiduals()(0) = mass2 - E * E + px * px + py * py + pz * pz;
504 p.getH()(0, momindex) = 2.0 * px;
505 p.getH()(0, momindex + 1) = 2.0 * py;
506 p.getH()(0, momindex + 2) = 2.0 * pz;
507 p.getH()(0, momindex + 3) = -2.0 * E;
515 return ErrCode(ErrCode::Status::success);
532 if (
type == Constraint::mass) {
535 B2FATAL(
"Trying to project constraint of ParticleBase type. This is undefined.");
556 const Eigen::Matrix < double, 1, -1, 1, 1, 3 > vertex_dist =
557 fitparams.getStateVector().segment(posindex,
dim) - fitparams.getStateVector().segment(mother_ps_index,
dim);
558 const Eigen::Matrix < double, 1, -1, 1, 1, 3 >
559 mom = fitparams.getStateVector().segment(posindex,
dim) - fitparams.getStateVector().segment(mother_ps_index,
dim);
564 const double mom_norm = mom.norm();
565 const double dot = std::abs(vertex_dist.dot(mom));
566 const double tau = dot / mom_norm;
567 if (0 == mom_norm || 0 == dot) {
570 fitparams.getStateVector()(tauindex) = tau;
574 return ErrCode(ErrCode::Status::success);