12 #include <analysis/dataobjects/Particle.h>
14 #include <framework/logging/Logger.h>
15 #include <framework/gearbox/Const.h>
17 #include <analysis/VertexFitting/TreeFitter/FitManager.h>
18 #include <analysis/VertexFitting/TreeFitter/FitParams.h>
19 #include <analysis/VertexFitting/TreeFitter/DecayChain.h>
20 #include <analysis/VertexFitting/TreeFitter/ParticleBase.h>
22 namespace TreeFitter {
25 const ConstraintConfiguration& config,
28 const bool useReferencing
31 m_decaychain(nullptr),
36 m_updateDaugthers(updateDaughters),
39 m_useReferencing(useReferencing),
46 m_fitparams =
new FitParams(m_decaychain->dim());
58 if (part->hasExtraInfo(name)) {
59 part->setExtraInfo(name, value);
61 part->addExtraInfo(name, value);
68 const int nitermax = 100;
69 const int maxndiverging = 3;
70 double dChisqConv =
m_prec;
83 bool finished =
false;
101 if ((std::abs(deltachisq) /
m_chiSquare < dChisqConv)) {
106 }
else if (deltachisq > 0 && ++ndiverging >= maxndiverging) {
108 m_status = VertexStatus::NonConverged;
109 m_errCode = ErrCode(ErrCode::Status::slowdivergingfit);
111 }
else if (deltachisq > 0) {
114 if (deltachisq < 0) {
122 m_status = VertexStatus::NonConverged;
130 if (
m_status == VertexStatus::Success) {
139 return (
m_status == VertexStatus::Success);
166 Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> cov =
m_fitparams->
getCovariance().selfadjointView<Eigen::Lower>();
169 if (posindex < 0 && pb->mother()) {
176 for (
int row = 0; row < 4; ++row) {
177 for (
int col = 0; col < 4; ++col) {
178 returncov(row, col) = cov(momindex + row, momindex + col);
182 for (
int row = 0; row < 3; ++row) {
183 for (
int col = 0; col < 3; ++col) {
184 returncov(row + 4, col + 4) = cov(posindex + row, posindex + col);
189 Eigen::Matrix<double, 6, 6> cov6 =
190 Eigen::Matrix<double, 6, 6>::Zero(6, 6);
192 for (
int row = 0; row < 3; ++row) {
193 for (
int col = 0; col < 3; ++col) {
194 cov6(row, col) = cov(momindex + row, momindex + col);
195 cov6(row + 3, col + 3) = cov(posindex + row, posindex + col);
200 Eigen::Matrix<double, 3, 1> momVec =
203 double energy2 = momVec.transpose() * momVec;
204 energy2 += mass * mass;
205 double energy = sqrt(energy2);
207 Eigen::Matrix<double, 7, 6> jacobian =
208 Eigen::Matrix<double, 7, 6>::Zero(7, 6);
210 for (
int col = 0; col < 3; ++col) {
211 jacobian(col, col) = 1;
213 jacobian(col + 4, col + 3) = 1;
216 Eigen::Matrix<double, 7, 7> cov7
217 = jacobian * cov6.selfadjointView<Eigen::Lower>() * jacobian.transpose();
219 for (
int row = 0; row < 7; ++row) {
220 for (
int col = 0; col < 7; ++col) {
221 returncov(row, col) = cov7(row, col);
235 return pb !=
nullptr;
242 if (posindex < 0 && pb.
mother()) {
255 setExtraInfo(&cand,
"modifiedPValue", TMath::Prob(fitparchi2, 3));
270 p.SetE(std::sqrt(p.Vect()*p.Vect() + mass * mass));
273 TMatrixFSym cov7b2(7);
280 std::tuple<double, double>life =
getLifeTime(cand);
281 setExtraInfo(&cand, std::string(
"decayLength"), std::get<0>(tau));
282 setExtraInfo(&cand, std::string(
"decayLengthErr"), std::get<1>(tau));
283 setExtraInfo(&cand, std::string(
"lifeTime"), std::get<0>(life));
284 setExtraInfo(&cand, std::string(
"lifeTimeErr"), std::get<1>(life));
290 const bool updateableMother =
updateCand(cand, isTreeHead);
292 if (updateableMother) {
294 for (
int i = 0; i < ndaughters; i++) {
305 if (pb && pb->tauIndex() >= 0 && pb->mother()) {
306 const int momindex = pb->
momIndex();
307 const int tauIndex = pb->tauIndex();
311 Eigen::Matrix<double, 4, 4> comb_cov = Eigen::Matrix<double, 4, 4>::Zero(4, 4);
315 const double lenErr = std::get<1>(lenTuple);
316 comb_cov(0, 0) = lenErr * lenErr;
321 comb_cov.block<3, 3>(1, 1) = mom_cov;
323 const double mass = pb->pdgMass();
325 const double mom = mom_vec.norm();
326 const double mom3 = mom * mom * mom;
328 const double len = std::get<0>(lenTuple);
329 const double t = len / mom * mBYc;
331 Eigen::Matrix<double, 1, 4> jac = Eigen::Matrix<double, 1, 4>::Zero();
332 jac(0) = 1. / mom * mBYc;
333 jac(1) = -1. * len * mom_vec(0) / mom3 * mBYc;
334 jac(2) = -1. * len * mom_vec(1) / mom3 * mBYc;
335 jac(3) = -1. * len * mom_vec(2) / mom3 * mBYc;
337 const double tErr2 = jac * comb_cov.selfadjointView<Eigen::Lower>() * jac.transpose();
339 return std::make_tuple(t, std::sqrt(tErr2));
341 return std::make_tuple(-999, -999);
352 if (pb->tauIndex() >= 0 && pb->mother()) {
353 const int tauindex = pb->tauIndex();
354 const double len = fitparams.getStateVector()(tauindex);
355 const double lenErr2 = fitparams.getCovariance()(tauindex, tauindex);
356 return std::make_tuple(len, std::sqrt(lenErr2));
358 return std::make_tuple(-999, -999);
363 std::tuple<double, double> rc = std::make_tuple(-999, -999);