this class More...

#include <FitManager.h>

Collaboration diagram for FitManager:

[legend]

Public Types
enum	VertexStatus { Success = 0 , NonConverged , BadInput , Failed , UnFitted }
	status flag of the fit-itereation (the step in the newton method) More...

Public Member Functions
	FitManager ()
	constructor

	FitManager (Belle2::Particle *particle, const ConstraintConfiguration &config, double prec=0.01, bool updateDaughters=false)
	constructor

	FitManager (const FitManager &other)=delete
	use default copy constructor

FitManager &	operator= (const FitManager &other)=delete
	use default assignment op

	~FitManager ()
	destructor does stuff

bool	fit ()
	main fit function that uses the kalman filter

bool	updateCand (Belle2::Particle &particle, const bool isTreeHead) const
	update particles parameters with the fit results

void	updateCand (const ParticleBase &pb, Belle2::Particle &cand, const bool isTreeHead) const
	locate particle base for a belle2 particle and update the particle with the values from particle base

void	updateTree (Belle2::Particle &particle, const bool isTreeHead) const
	update the Belle2::Particles with the fit results

void	getCovFromPB (const ParticleBase *pb, TMatrixFSym &returncov) const
	extract cov from particle base

std::tuple< double, double >	getLifeTime (Belle2::Particle &cand) const
	get lifetime

std::tuple< double, double >	getDecayLength (const ParticleBase *pb) const
	get decay length

std::tuple< double, double >	getDecayLength (Belle2::Particle &cand) const
	get decay length

Belle2::Particle *	particle ()
	getter for the head of the tree

Static Public Member Functions
static std::tuple< double, double >	getDecayLength (const ParticleBase *pb, const FitParams &fitparams)
	get decay length

Private Attributes
Belle2::Particle *	m_particle
	head of the tree

DecayChain *	m_decaychain
	the decay tree

int	m_status
	status of the current iteration

double	m_chiSquare
	chi2 of the current iteration

double	m_prec
	precision that is needed for status:converged (delta chi2)

ErrCode	m_errCode
	errorcode

const bool	m_updateDaugthers
	if this is set all daughters will be updated otherwise only the head of the tree

int	m_ndf
	number of degrees of freedom for this topology

FitParams *	m_fitparams
	parameters to be fitted

const ConstraintConfiguration	m_config
	config container

Detailed Description

this class

Definition at line 29 of file FitManager.h.

Member Enumeration Documentation

◆ VertexStatus

enum VertexStatus

status flag of the fit-itereation (the step in the newton method)

Definition at line 33 of file FitManager.h.

33{ Success = 0, NonConverged, BadInput, Failed, UnFitted };

Constructor & Destructor Documentation

◆ FitManager() [1/2]

FitManager ( )

inline

constructor

Definition at line 36 of file FitManager.h.

                 : m_particle(0), m_decaychain(0), m_status(VertexStatus::UnFitted),
      m_chiSquare(-1), m_prec(0.01), m_updateDaugthers(false), m_ndf(0),
      m_fitparams(0)
    {}

◆ FitManager() [2/2]

FitManager	(	Belle2::Particle *	particle,
		const ConstraintConfiguration &	config,
		double	prec = 0.01,
		bool	updateDaughters = false )

constructor

Definition at line 24 of file FitManager.cc.

                          :
    m_particle(particle),
    m_decaychain(nullptr),
    m_status(VertexStatus::UnFitted),
    m_chiSquare(-1),
    m_prec(prec),
    m_updateDaugthers(updateDaughters),
    m_ndf(0),
    m_fitparams(nullptr),
    m_config(config)
  {
    m_decaychain = new DecayChain(particle, config, false);
    m_fitparams  = new FitParams(m_decaychain->dim());
  }

◆ ~FitManager()

~FitManager ( )

destructor does stuff

Definition at line 43 of file FitManager.cc.

  {
    delete m_decaychain;
    delete m_fitparams;
  }

Member Function Documentation

◆ fit()

bool fit ( )

main fit function that uses the kalman filter

Definition at line 49 of file FitManager.cc.

  {
    const int nitermax = 100;
    const int maxndiverging = 3;
    const double dChisqConv = m_prec;
    m_chiSquare = -1;
    m_errCode.reset();
 
    if (m_status == VertexStatus::UnFitted) {
      m_errCode = m_decaychain->initialize(*m_fitparams);
    }
 
    if (m_errCode.failure()) {
      m_status = VertexStatus::BadInput;
    } else {
      m_status = VertexStatus::UnFitted;
      int ndiverging = 0;
      bool finished = false;
      int niter = 0;
      for (niter = 0; niter < nitermax && !finished; ++niter) {
        if (niter == 0) {
          m_errCode = m_decaychain->filter(*m_fitparams);
        } else {
          auto* tempState = new FitParams(*m_fitparams);
          m_errCode = m_decaychain->filterWithReference(*m_fitparams, *tempState);
          delete tempState;
        }
        m_ndf = m_fitparams->nDof();
        double chisq = m_fitparams->chiSquare();
        double deltachisq = chisq - m_chiSquare;
        if (m_errCode.failure()) {
          finished = true ;
          m_status = VertexStatus::Failed;
          m_particle->writeExtraInfo("failed", 1);
        } else {
          if (niter > 0) {
            if ((std::abs(deltachisq) / m_chiSquare < dChisqConv)) {
              m_chiSquare = chisq;
              m_status = VertexStatus::Success;
              finished = true ;
              m_particle->writeExtraInfo("failed", 0);
            } else if (deltachisq > 0 && ++ndiverging >= maxndiverging) {
              m_particle->writeExtraInfo("failed", 2);
              m_status = VertexStatus::NonConverged;
              m_errCode = ErrCode(ErrCode::Status::slowdivergingfit);
              finished = true ;
            }
          }
          if (deltachisq < 0) {
            ndiverging = 0;
          }
          m_chiSquare = chisq;
        }
      }
      if (niter == nitermax && m_status != VertexStatus::Success) {
        m_particle->writeExtraInfo("failed", 3);
        m_status = VertexStatus::NonConverged;
      }
      if (!(m_fitparams->testCovariance())) {
        m_particle->writeExtraInfo("failed", 4);
        m_status = VertexStatus::Failed;
      }
    }
 
    if (m_status == VertexStatus::Success) {
      // mass constraints comes after kine so we have to
      // update the mothers with the values set by the mass constraint
      if (m_config.m_massConstraintListPDG.size() != 0) {
        m_decaychain->locate(m_particle)->forceP4Sum(*m_fitparams);
      }
      updateTree(*m_particle, true);
    }
 
    return (m_status == VertexStatus::Success);
  }

◆ getCovFromPB()

void getCovFromPB	(	const ParticleBase *	pb,
		TMatrixFSym &	returncov ) const

extract cov from particle base

Definition at line 128 of file FitManager.cc.

  {
 
    Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> cov = m_fitparams->getCovariance().selfadjointView<Eigen::Lower>();
    int posindex = pb->posIndex();
    // hack: for tracks and photons, use the production vertex
    if (posindex < 0 && pb->mother()) {
      posindex = pb->mother()->posIndex();
    }
    int momindex = pb->momIndex();
    if (pb->hasEnergy()) {
      // if particle has energy, get full p4 from fitparams and put them directly in the return type
      // very important! Belle2 uses p,E,x! Change order here!
      for (int row = 0; row < 4; ++row) {
        for (int col = 0; col < 4; ++col) {
          returncov(row, col) = cov(momindex + row, momindex + col);
        }
      }
 
      for (int row = 0; row < 3; ++row) {
        for (int col = 0; col < 3; ++col) {
          returncov(row + 4, col + 4) = cov(posindex + row, posindex + col);
        }
      }
 
    } else {
      Eigen::Matrix<double, 6, 6> cov6 =
        Eigen::Matrix<double, 6, 6>::Zero(6, 6);
 
      for (int row = 0; row < 3; ++row) {
        for (int col = 0; col < 3; ++col) {
          cov6(row, col) = cov(momindex + row, momindex + col);
          cov6(row + 3, col + 3) = cov(posindex + row, posindex + col);
        }
      }
 
      double mass = 0;
      if (pb->particle()->hasExtraInfo("treeFitterMassConstraintValue")) {
        mass = pb->particle()->getExtraInfo("treeFitterMassConstraintValue");
      } else mass = pb->particle()->getPDGMass();
      Eigen::Matrix<double, 3, 1> momVec =
        m_fitparams->getStateVector().segment(momindex, 3);
 
      double energy2 = momVec.transpose() * momVec;
      energy2 += mass * mass;
      double energy = sqrt(energy2);
 
      Eigen::Matrix<double, 7, 6> jacobian =
        Eigen::Matrix<double, 7, 6>::Zero(7, 6);
 
      for (int col = 0; col < 3; ++col) {
        jacobian(col, col) = 1;
        jacobian(3, col) = m_fitparams->getStateVector()(momindex + col) / energy;
        jacobian(col + 4, col + 3) = 1;
      }
 
      Eigen::Matrix<double, 7, 7> cov7
        = jacobian * cov6.selfadjointView<Eigen::Lower>() * jacobian.transpose();
 
      for (int row = 0; row < 7; ++row) {
        for (int col = 0; col < 7; ++col) {
          returncov(row, col) = cov7(row, col);
        }
      }
    } // else
  }

◆ getDecayLength() [1/3]

std::tuple< double, double > getDecayLength ( Belle2::Particle & cand ) const

get decay length

Definition at line 361 of file FitManager.cc.

  {
    std::tuple<double, double> rc = std::make_tuple(-999, -999);
    const ParticleBase* pb = m_decaychain->locate(&cand) ;
    if (pb && pb->tauIndex() >= 0 && pb->mother()) {
      rc = getDecayLength(pb);
    }
    return rc;
  }

◆ getDecayLength() [2/3]

std::tuple< double, double > getDecayLength ( const ParticleBase * pb ) const

get decay length

Definition at line 344 of file FitManager.cc.

  {
    // returns the decaylength in the lab frame
    return getDecayLength(pb, *m_fitparams);
  }

◆ getDecayLength() [3/3]

std::tuple< double, double > getDecayLength	(	const ParticleBase *	pb,
		const FitParams &	fitparams )

static

get decay length

Definition at line 350 of file FitManager.cc.

  {
    if (pb->tauIndex() >= 0 && pb->mother()) {
      const int tauindex = pb->tauIndex();
      const double len = fitparams.getStateVector()(tauindex);
      const double lenErr2 = fitparams.getCovariance()(tauindex, tauindex);
      return std::make_tuple(len, std::sqrt(lenErr2));
    }
    return std::make_tuple(-999, -999);
  }

◆ getLifeTime()

std::tuple< double, double > getLifeTime ( Belle2::Particle & cand ) const

get lifetime

Definition at line 298 of file FitManager.cc.

  {
    const ParticleBase* pb = m_decaychain->locate(&cand);
 
    if (pb && pb->tauIndex() >= 0 && pb->mother()) {
      const int momindex = pb->momIndex();
      const int tauIndex = pb->tauIndex();
      const Eigen::Matrix<double, 1, 3> mom_vec = m_fitparams->getStateVector().segment(momindex, 3);
 
      const Eigen::Matrix<double, 3, 3> mom_cov = m_fitparams->getCovariance().block<3, 3>(momindex, momindex);
      Eigen::Matrix<double, 4, 4> comb_cov =  Eigen::Matrix<double, 4, 4>::Zero(4, 4);
 
      const std::tuple<double, double> lenTuple  = getDecayLength(cand);
 
      const double lenErr = std::get<1>(lenTuple);
      comb_cov(0, 0) = lenErr * lenErr;
      comb_cov(1, 0) = m_fitparams->getCovariance()(momindex, tauIndex);
      comb_cov(2, 0) = m_fitparams->getCovariance()(momindex + 1, tauIndex);
      comb_cov(3, 0) = m_fitparams->getCovariance()(momindex + 2, tauIndex);
 
      comb_cov.block<3, 3>(1, 1) = mom_cov;
 
      double mass = 0;
      if (pb->particle()->hasExtraInfo("treeFitterMassConstraintValue")) {
        mass = pb->particle()->getExtraInfo("treeFitterMassConstraintValue");
      } else mass = pb->particle()->getPDGMass();
      const double mBYc = mass / Belle2::Const::speedOfLight;
      const double mom = mom_vec.norm();
      const double mom3 = mom * mom * mom;
 
      const double len = std::get<0>(lenTuple);
      const double t = len / mom * mBYc;
 
      Eigen::Matrix<double, 1, 4> jac = Eigen::Matrix<double, 1, 4>::Zero();
      jac(0) =  1. / mom * mBYc;
      jac(1) = -1. * len * mom_vec(0) / mom3 * mBYc;
      jac(2) = -1. * len * mom_vec(1) / mom3 * mBYc;
      jac(3) = -1. * len * mom_vec(2) / mom3 * mBYc;
 
      const double tErr2 = jac * comb_cov.selfadjointView<Eigen::Lower>() * jac.transpose();
      // time in nanosec
      return std::make_tuple(t, std::sqrt(tErr2));
    }
    return std::make_tuple(-999, -999);
  }

◆ particle()

Belle2::Particle * particle ( )

inline

getter for the head of the tree

Definition at line 85 of file FitManager.h.

85{ return m_particle; }

◆ updateCand() [1/2]

bool updateCand	(	Belle2::Particle &	particle,
		const bool	isTreeHead ) const

update particles parameters with the fit results

Definition at line 195 of file FitManager.cc.

  {
    const ParticleBase* pb = m_decaychain->locate(&cand);
    if (pb) {
      updateCand(*pb, cand, isTreeHead);
    } else {
      B2ERROR("Can't find candidate " << cand.getName() << " in tree " << m_particle->getName());
    }
    return pb != nullptr;
  }

◆ updateCand() [2/2]

void updateCand	(	const ParticleBase &	pb,
		Belle2::Particle &	cand,
		const bool	isTreeHead ) const

locate particle base for a belle2 particle and update the particle with the values from particle base

Definition at line 206 of file FitManager.cc.

  {
    int posindex = pb.posIndex();
    if (posindex < 0 && pb.mother()) {
      posindex = pb.mother()->posIndex();
    }
 
    if (m_updateDaugthers || isTreeHead) {
      TMatrixFSym cov7b2(7);
      if (posindex >= 0) {
        const ROOT::Math::XYZVector pos(m_fitparams->getStateVector()(posindex),
                                        m_fitparams->getStateVector()(posindex + 1),
                                        m_fitparams->getStateVector()(posindex + 2));
        cand.setVertex(pos);
        if (&pb == m_decaychain->cand()) { // if head
          const double fitparchi2 = m_fitparams->chiSquare();
          cand.setPValue(TMath::Prob(fitparchi2, m_ndf));//if m_ndf<1, this is 0.
          cand.writeExtraInfo("chiSquared", fitparchi2);
          cand.writeExtraInfo("modifiedPValue", TMath::Prob(fitparchi2, 3));
          cand.writeExtraInfo("ndf", m_ndf);
        }
        if (pb.mother()) {
          int motherPosIndex = pb.mother()->posIndex();
          if (motherPosIndex >= 0) {
            cand.writeExtraInfo("prodVertexX", m_fitparams->getStateVector()(motherPosIndex));
            cand.writeExtraInfo("prodVertexY", m_fitparams->getStateVector()(motherPosIndex + 1));
            if (pb.mother()->dim() > 2)
              cand.writeExtraInfo("prodVertexZ", m_fitparams->getStateVector()(motherPosIndex + 2));
            if (not isTreeHead) {
              getCovFromPB(pb.mother(), cov7b2);
              cand.writeExtraInfo("prodVertSxx", cov7b2[4][4]);
              cand.writeExtraInfo("prodVertSxy", cov7b2[4][5]);
              cand.writeExtraInfo("prodVertSyx", cov7b2[5][4]);
              cand.writeExtraInfo("prodVertSyy", cov7b2[5][5]);
              if (pb.mother()->dim() > 2) {
                cand.writeExtraInfo("prodVertexZ", m_fitparams->getStateVector()(motherPosIndex + 2));
                cand.writeExtraInfo("prodVertSxz", cov7b2[4][6]);
                cand.writeExtraInfo("prodVertSyz", cov7b2[5][6]);
                cand.writeExtraInfo("prodVertSzx", cov7b2[6][4]);
                cand.writeExtraInfo("prodVertSzy", cov7b2[6][5]);
                cand.writeExtraInfo("prodVertSzz", cov7b2[6][6]);
              }
            }
          }
        }
      }
 
      const int momindex = pb.momIndex();
      ROOT::Math::PxPyPzEVector p;
      p.SetPx(m_fitparams->getStateVector()(momindex));
      p.SetPy(m_fitparams->getStateVector()(momindex + 1));
      p.SetPz(m_fitparams->getStateVector()(momindex + 2));
      if (pb.hasEnergy()) {
        p.SetE(m_fitparams->getStateVector()(momindex + 3));
        cand.set4VectorDividingByMomentumScaling(p);
      } else {
        double mass = 0;
        if (cand.hasExtraInfo("treeFitterMassConstraintValue")) {
          mass = cand.getExtraInfo("treeFitterMassConstraintValue");
        } else mass = cand.getPDGMass();
        p.SetE(std::sqrt(p.P2() + mass * mass));
        cand.set4VectorDividingByMomentumScaling(p);
      }
      getCovFromPB(&pb, cov7b2);
      cand.setMomentumVertexErrorMatrix(cov7b2);
    }
 
    if (pb.tauIndex() > 0) {
      const std::tuple<double, double>tau  = getDecayLength(cand);
      const std::tuple<double, double>life = getLifeTime(cand);
      cand.writeExtraInfo("decayLength", std::get<0>(tau));
      cand.writeExtraInfo("decayLengthErr", std::get<1>(tau));
      cand.writeExtraInfo("lifeTime", std::get<0>(life));
      cand.writeExtraInfo("lifeTimeErr", std::get<1>(life));
    }
  }

◆ updateTree()

void updateTree	(	Belle2::Particle &	particle,
		const bool	isTreeHead ) const

update the Belle2::Particles with the fit results

Definition at line 284 of file FitManager.cc.

  {
    const bool updateableMother = updateCand(cand, isTreeHead);
 
    if (updateableMother and not cand.hasExtraInfo("bremsCorrected") and
        not(cand.hasExtraInfo("treeFitterTreatMeAsInvisible") and cand.getExtraInfo("treeFitterTreatMeAsInvisible") == 1)) {
      const int ndaughters = cand.getNDaughters();
      for (int i = 0; i < ndaughters; i++) {
        auto* daughter = const_cast<Belle2::Particle*>(cand.getDaughter(i));
        updateTree(*daughter, false);
      }
    }
  }

Member Data Documentation

◆ m_chiSquare

double m_chiSquare

private

chi2 of the current iteration

Definition at line 98 of file FitManager.h.

◆ m_config

const ConstraintConfiguration m_config

private

config container

Definition at line 116 of file FitManager.h.

◆ m_decaychain

DecayChain* m_decaychain

private

the decay tree

Definition at line 92 of file FitManager.h.

◆ m_errCode

ErrCode m_errCode

private

errorcode

Definition at line 104 of file FitManager.h.

◆ m_fitparams

FitParams* m_fitparams

private

parameters to be fitted

Definition at line 113 of file FitManager.h.

◆ m_ndf

int m_ndf

private

number of degrees of freedom for this topology

Definition at line 110 of file FitManager.h.

◆ m_particle

Belle2::Particle* m_particle

private

head of the tree

Definition at line 89 of file FitManager.h.

◆ m_prec

double m_prec

private

precision that is needed for status:converged (delta chi2)

Definition at line 101 of file FitManager.h.

◆ m_status

int m_status

private

status of the current iteration

Definition at line 95 of file FitManager.h.

◆ m_updateDaugthers

const bool m_updateDaugthers

private

if this is set all daughters will be updated otherwise only the head of the tree

Definition at line 107 of file FitManager.h.

The documentation for this class was generated from the following files:

analysis/VertexFitting/TreeFitter/include/FitManager.h
analysis/VertexFitting/TreeFitter/src/FitManager.cc

Public Types

Public Member Functions

Static Public Member Functions

Private Attributes

Detailed Description

Member Enumeration Documentation

◆ VertexStatus

Constructor & Destructor Documentation

◆ FitManager() [1/2]

◆ FitManager() [2/2]

◆ ~FitManager()

Member Function Documentation

◆ fit()

◆ getCovFromPB()

◆ getDecayLength() [1/3]

◆ getDecayLength() [2/3]

◆ getDecayLength() [3/3]

◆ getLifeTime()

◆ particle()

◆ updateCand() [1/2]

◆ updateCand() [2/2]

◆ updateTree()

Member Data Documentation

◆ m_chiSquare

◆ m_config

◆ m_decaychain

◆ m_errCode

◆ m_fitparams

◆ m_ndf

◆ m_particle

◆ m_prec

◆ m_status

◆ m_updateDaugthers