FinalTreeBuilder Class Reference

#include <FinalTreeBuilder.h>

Public Member Functions
RefCountedKinematicTree	buildTree (const CachingVertex< 6 > &vtx, const std::vector< RefCountedKinematicParticle > &input) const
	FinalTreeBuilder ()
	~FinalTreeBuilder ()

Private Types
typedef ReferenceCountingPointer< LinearizedTrackState< 6 > >	RefCountedLinearizedTrackState
typedef ReferenceCountingPointer< RefittedTrackState< 6 > >	RefCountedRefittedTrackState
typedef ReferenceCountingPointer< VertexTrack< 6 > >	RefCountedVertexTrack

Private Member Functions
AlgebraicMatrix	momentumPart (const CachingVertex< 6 > &vtx, const AlgebraicVector7 &par) const

Private Attributes
KinematicVertexFactory *	kvFactory
VirtualKinematicParticleFactory *	pFactory

Detailed Description

Class building a resulting output tree out of the information provided by KinematicParticleVertexFitter.

Definition at line 16 of file FinalTreeBuilder.h.

Member Typedef Documentation

RefCountedLinearizedTrackState

typedef ReferenceCountingPointer<LinearizedTrackState<6> > FinalTreeBuilder::RefCountedLinearizedTrackState

private

Definition at line 27 of file FinalTreeBuilder.h.

RefCountedRefittedTrackState

typedef ReferenceCountingPointer<RefittedTrackState<6> > FinalTreeBuilder::RefCountedRefittedTrackState

private

Definition at line 28 of file FinalTreeBuilder.h.

RefCountedVertexTrack

typedef ReferenceCountingPointer<VertexTrack<6> > FinalTreeBuilder::RefCountedVertexTrack

private

Definition at line 26 of file FinalTreeBuilder.h.

Constructor & Destructor Documentation

FinalTreeBuilder()

FinalTreeBuilder::FinalTreeBuilder

(

)

Definition at line 8 of file FinalTreeBuilder.cc.

References kvFactory, and pFactory.

                                   {
  kvFactory = new KinematicVertexFactory();
  KinematicStatePropagator* ksp = nullptr;
  pFactory = new VirtualKinematicParticleFactory(ksp);
}

~FinalTreeBuilder()

FinalTreeBuilder::~FinalTreeBuilder

(

)

Definition at line 14 of file FinalTreeBuilder.cc.

References kvFactory, and pFactory.

                                    {
  delete kvFactory;
  delete pFactory;
}

Member Function Documentation

buildTree()

RefCountedKinematicTree FinalTreeBuilder::buildTree	(	const CachingVertex< 6 > &	vtx,
		const std::vector< RefCountedKinematicParticle > &	input
	)		const

Definition at line 19 of file FinalTreeBuilder.cc.

References mps_fire::i, input, dqmiolumiharvest::j, KinematicRefittedTrackState::kinematicMomentumVector(), KinematicRefittedTrackState::kinematicParameters(), KinematicRefittedTrackState::kinematicParametersCovariance(), kvFactory, LogDebug, momentumPart(), nPart(), hcaldqm::flag::nState, VirtualKinematicParticleFactory::particle(), pFactory, mathSSE::sqrt(), KinematicVertexFactory::vertex(), and extraflags_cff::vtx.

Referenced by KinematicParticleVertexFitter::fit().

                                                                                                               {
  //creating resulting kinematic particle
  AlgebraicVector7 par;
  AlgebraicMatrix cov(7, 7, 0);
  par(0) = vtx.position().x();
  par(1) = vtx.position().y();
  par(2) = vtx.position().z();
  double en = 0.;
  int ch = 0;

  //new particle momentum calculation and refitted kinematic states
  std::vector<KinematicRefittedTrackState*> rStates;
  std::vector<RefCountedVertexTrack> refTracks = vtx.tracks();
  for (std::vector<RefCountedVertexTrack>::const_iterator i = refTracks.begin(); i != refTracks.end(); ++i) {
    KinematicRefittedTrackState* rs = dynamic_cast<KinematicRefittedTrackState*>(&(*((*i)->refittedState())));
    AlgebraicVector4 f_mom = rs->kinematicMomentumVector();
    par(3) += f_mom(0);
    par(4) += f_mom(1);
    par(5) += f_mom(2);
    en += sqrt(f_mom(1) * f_mom(1) + f_mom(2) * f_mom(2) + f_mom(3) * f_mom(3) + f_mom(0) * f_mom(0));
    ch += (*i)->linearizedTrack()->charge();
    rStates.push_back(rs);
  }

  //math precision check (numerical stability)
  double differ = en * en - (par(3) * par(3) + par(4) * par(4) + par(5) * par(5));
  if (differ > 0.) {
    par(6) = sqrt(differ);
  } else {
    LogDebug("FinalTreeBuilder") << "Fit failed: Current precision does not allow to calculate the mass\n";
    return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
  }

  // covariance matrix calculation: momentum-momentum components part (4x4)
  // and position-momentum components part:
  AlgebraicMatrix m_all = momentumPart(vtx, par);

  //position-position components part (3x3)
  // AlgebraicMatrix x_X = vtx.error().matrix();

  //making new matrix itself
  cov.sub(1, 1, m_all);

  //covariance sym matrix
  AlgebraicSymMatrix77 sCov;
  for (int i = 1; i < 8; i++) {
    for (int j = 1; j < 8; j++) {
      if (i <= j)
        sCov(i - 1, j - 1) = cov(i, j);
    }
  }

  //valid decay vertex for our resulting particle
  RefCountedKinematicVertex dVrt = kvFactory->vertex(vtx);

  //invalid production vertex for our resulting particle
  RefCountedKinematicVertex pVrt = kvFactory->vertex();

  //new born kinematic particle
  KinematicParameters kPar(par);
  KinematicParametersError kEr(sCov);
  const MagneticField* field = input.front()->magneticField();
  KinematicState nState(kPar, kEr, ch, field);

  //invalid previous particle and empty constraint:
  KinematicParticle* zp = nullptr;
  RefCountedKinematicParticle pPart = ReferenceCountingPointer<KinematicParticle>(zp);

  float vChi = vtx.totalChiSquared();
  float vNdf = vtx.degreesOfFreedom();
  RefCountedKinematicParticle nPart = pFactory->particle(nState, vChi, vNdf, pPart);

  //adding top particle to the tree
  RefCountedKinematicTree resTree = ReferenceCountingPointer<KinematicTree>(new KinematicTree());
  resTree->addParticle(pVrt, dVrt, nPart);

  //making refitted kinematic particles and putting them to the tree
  std::vector<RefCountedKinematicParticle> rrP;

  std::vector<RefCountedKinematicParticle>::const_iterator j;
  std::vector<RefCountedVertexTrack>::const_iterator i;
  for (j = input.begin(), i = refTracks.begin(); j != input.end() && i != refTracks.end(); ++j, ++i) {
    RefCountedLinearizedTrackState lT = (*i)->linearizedTrack();
    KinematicRefittedTrackState* rS = dynamic_cast<KinematicRefittedTrackState*>(&(*((*i)->refittedState())));

    //   RefCountedRefittedTrackState rS = (*i)->refittedState();
    AlgebraicVector7 lPar = rS->kinematicParameters();
    KinematicParameters lkPar(lPar);
    AlgebraicSymMatrix77 lCov = rS->kinematicParametersCovariance();
    KinematicParametersError lkCov(lCov);
    TrackCharge lch = lT->charge();
    KinematicState nState(lkPar, lkCov, lch, field);
    RefCountedKinematicParticle nPart = (*j)->refittedParticle(nState, vChi, vNdf);
    rrP.push_back(nPart);
    if ((*j)->correspondingTree() != nullptr) {
      //here are the particles having trees after them
      KinematicTree* tree = (*j)->correspondingTree();
      tree->movePointerToTheTop();
      tree->replaceCurrentParticle(nPart);
      RefCountedKinematicVertex cdVertex = resTree->currentDecayVertex();
      resTree->addTree(cdVertex, tree);
    } else {
      //here are just particles fitted to this tree
      RefCountedKinematicVertex nV = kvFactory->vertex();
      resTree->addParticle(dVrt, nV, nPart);
    }
  }
  return resTree;
}

momentumPart()

AlgebraicMatrix FinalTreeBuilder::momentumPart ( const CachingVertex< 6 > &  vtx, const AlgebraicVector7 &  par  ) const

private

Definition at line 132 of file FinalTreeBuilder.cc.

References a, funct::cos(), Exception, EgHLTOffHistBins_cfi::mass, phi, rho, funct::sin(), findQualityFiles::size, mathSSE::sqrt(), funct::tan(), theta(), extraflags_cff::vtx, and z.

Referenced by buildTree().

                                                                                                             {
  std::vector<RefCountedVertexTrack> refTracks = vtx.tracks();
  int size = refTracks.size();
  AlgebraicMatrix cov(7 + 4 * (size - 1), 7 + 4 * (size - 1));
  AlgebraicMatrix jac(7, 7 + 4 * (size - 1));
  jac(1, 1) = 1.;
  jac(2, 2) = 1.;
  jac(3, 3) = 1.;

  double energy_total = sqrt(par(3) * par(3) + par(6) * par(6) + par(5) * par(5) + par(4) * par(4));

  std::vector<RefCountedVertexTrack>::const_iterator rt_i;
  int i_int = 0;
  for (rt_i = refTracks.begin(); rt_i != refTracks.end(); rt_i++) {
    double a;
    AlgebraicVector6 param = (**rt_i).refittedState()->parameters();  // rho, theta, phi,tr_im, z_im, mass
    double rho = param[0];
    double theta = param[1];
    double phi = param[2];
    double mass = param[5];

    if ((**rt_i).linearizedTrack()->charge() != 0) {
      a = -(**rt_i).refittedState()->freeTrajectoryState().parameters().magneticFieldInInverseGeV(vtx.position()).z() *
          (**rt_i).refittedState()->freeTrajectoryState().parameters().charge();
      if (a == 0.)
        throw cms::Exception("FinalTreeBuilder", "Field is 0");
    } else {
      a = 1;
    }

    AlgebraicMatrix jc_el(4, 4, 0);
    jc_el(1, 1) = -a * cos(phi) / (rho * rho);    //dpx/d rho
    jc_el(2, 1) = -a * sin(phi) / (rho * rho);    //dpy/d rho
    jc_el(3, 1) = -a / (rho * rho * tan(theta));  //dpz/d rho

    jc_el(3, 2) = -a / (rho * sin(theta) * sin(theta));  //dpz/d theta

    jc_el(1, 3) = -a * sin(phi) / rho;  //dpx/d phi
    jc_el(2, 3) = a * cos(phi) / rho;   //dpy/d

    //non-trival elements: mass correlations:
    double energy_local = sqrt(a * a / (rho * rho) * (1 + 1 / (tan(theta) * tan(theta))) + mass * mass);

    jc_el(4, 4) = energy_total * mass / (par(6) * energy_local);  // dm/dm

    jc_el(4, 1) = (-(energy_total / energy_local * a * a / (rho * rho * rho * sin(theta) * sin(theta))) +
                   par(3) * a / (rho * rho) * cos(phi) + par(4) * a / (rho * rho) * sin(phi) +
                   par(5) * a / (rho * rho * tan(theta))) /
                  par(6);  //dm / drho

    jc_el(4, 2) = (-(energy_total / energy_local * a * a / (rho * rho * sin(theta) * sin(theta) * tan(theta))) +
                   par(5) * a / (rho * sin(theta) * sin(theta))) /
                  par(6);  //dm d theta

    jc_el(4, 3) = (par(3) * sin(phi) - par(4) * cos(phi)) * a / (rho * par(6));  //dm/dphi

    jac.sub(4, i_int * 4 + 4, jc_el);

    //top left corner elements
    if (i_int == 0) {
      cov.sub(1, 1, asHepMatrix<7>((**rt_i).fullCovariance()));
    } else {
      //4-momentum corellatons: diagonal elements of the matrix
      AlgebraicMatrix fullCovMatrix(asHepMatrix<7>((**rt_i).fullCovariance()));
      AlgebraicMatrix m_m_cov = fullCovMatrix.sub(4, 7, 4, 7);
      AlgebraicMatrix x_p_cov = fullCovMatrix.sub(1, 3, 4, 7);
      AlgebraicMatrix p_x_cov = fullCovMatrix.sub(4, 7, 1, 3);

      // cout << "Full covariance: \n"<< (**rt_i).fullCovariance()<<endl;
      // cout << "Full m_m_cov: "<< m_m_cov<<endl;
      //  cout<<"p_x_cov"<< p_x_cov<<endl;
      //  cout<<"x_p_cov"<< x_p_cov<<endl;

      //putting everything to the joint covariance matrix:
      //diagonal momentum-momentum elements:
      cov.sub(i_int * 4 + 4, i_int * 4 + 4, m_m_cov);

      //position momentum elements:
      cov.sub(1, i_int * 4 + 4, x_p_cov);
      cov.sub(i_int * 4 + 4, 1, p_x_cov);

      //off diagonal elements: corellations
      // track momentum - track momentum
    }
    int j_int = 0;
    for (std::vector<RefCountedVertexTrack>::const_iterator rt_j = refTracks.begin(); rt_j != refTracks.end(); rt_j++) {
      if (i_int < j_int) {
        AlgebraicMatrix i_k_cov_m = asHepMatrix<4, 4>(vtx.tkToTkCovariance((*rt_i), (*rt_j)));
        //     cout<<"i_k_cov_m"<<i_k_cov_m <<endl;
        cov.sub(i_int * 4 + 4, j_int * 4 + 4, i_k_cov_m);
        cov.sub(j_int * 4 + 4, i_int * 4 + 4, i_k_cov_m.T());
      }
      j_int++;
    }
    i_int++;
  }
  // cout<<"jac"<<jac<<endl;
  // cout<<"cov"<<cov<<endl;
  //  cout << "final result new"<<jac*cov*jac.T()<<endl;

  return jac * cov * jac.T();
}

Member Data Documentation

kvFactory

KinematicVertexFactory* FinalTreeBuilder::kvFactory

private

Definition at line 33 of file FinalTreeBuilder.h.

Referenced by buildTree(), FinalTreeBuilder(), and ~FinalTreeBuilder().

pFactory

VirtualKinematicParticleFactory* FinalTreeBuilder::pFactory

private

Definition at line 34 of file FinalTreeBuilder.h.

Referenced by buildTree(), FinalTreeBuilder(), and ~FinalTreeBuilder().