KinematicFit/src/FinalTreeBuilder.cc

0001 #include "RecoVertex/KinematicFit/interface/FinalTreeBuilder.h"
0002 #include "RecoVertex/KinematicFitPrimitives/interface/KinematicRefittedTrackState.h"
0003 #include "RecoVertex/KinematicFitPrimitives/interface/KinematicStatePropagator.h"
0004 #include "RecoVertex/KinematicFitPrimitives/interface/KinematicState.h"
0005 //#include "Vertex/KinematicFitPrimitives/interface/KinematicLinearizedTrackState.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 FinalTreeBuilder::FinalTreeBuilder() {
0009   kvFactory = new KinematicVertexFactory();
0010   KinematicStatePropagator* ksp = nullptr;
0011   pFactory = new VirtualKinematicParticleFactory(ksp);
0012 }
0013
0014 FinalTreeBuilder::~FinalTreeBuilder() {
0015   delete kvFactory;
0016   delete pFactory;
0017 }
0018
0019 RefCountedKinematicTree FinalTreeBuilder::buildTree(const CachingVertex<6>& vtx,
0020                                                     const std::vector<RefCountedKinematicParticle>& input) const {
0021   //creating resulting kinematic particle
0022   AlgebraicVector7 par;
0023   AlgebraicMatrix cov(7, 7, 0);
0024   par(0) = vtx.position().x();
0025   par(1) = vtx.position().y();
0026   par(2) = vtx.position().z();
0027   double en = 0.;
0028   int ch = 0;
0029
0030   //new particle momentum calculation and refitted kinematic states
0031   std::vector<KinematicRefittedTrackState*> rStates;
0032   std::vector<RefCountedVertexTrack> refTracks = vtx.tracks();
0033   for (std::vector<RefCountedVertexTrack>::const_iterator i = refTracks.begin(); i != refTracks.end(); ++i) {
0034     KinematicRefittedTrackState* rs = dynamic_cast<KinematicRefittedTrackState*>(&(*((*i)->refittedState())));
0035     AlgebraicVector4 f_mom = rs->kinematicMomentumVector();
0036     par(3) += f_mom(0);
0037     par(4) += f_mom(1);
0038     par(5) += f_mom(2);
0039     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));
0040     ch += (*i)->linearizedTrack()->charge();
0041     rStates.push_back(rs);
0042   }
0043
0044   //math precision check (numerical stability)
0045   double differ = en * en - (par(3) * par(3) + par(4) * par(4) + par(5) * par(5));
0046   if (differ > 0.) {
0047     par(6) = sqrt(differ);
0048   } else {
0049     LogDebug("FinalTreeBuilder") << "Fit failed: Current precision does not allow to calculate the mass\n";
0050     return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
0051   }
0052
0053   // covariance matrix calculation: momentum-momentum components part (4x4)
0054   // and position-momentum components part:
0055   AlgebraicMatrix m_all = momentumPart(vtx, par);
0056
0057   //position-position components part (3x3)
0058   // AlgebraicMatrix x_X = vtx.error().matrix();
0059
0060   //making new matrix itself
0061   cov.sub(1, 1, m_all);
0062
0063   //covariance sym matrix
0064   AlgebraicSymMatrix77 sCov;
0065   for (int i = 1; i < 8; i++) {
0066     for (int j = 1; j < 8; j++) {
0067       if (i <= j)
0068         sCov(i - 1, j - 1) = cov(i, j);
0069     }
0070   }
0071
0072   //valid decay vertex for our resulting particle
0073   RefCountedKinematicVertex dVrt = kvFactory->vertex(vtx);
0074
0075   //invalid production vertex for our resulting particle
0076   RefCountedKinematicVertex pVrt = kvFactory->vertex();
0077
0078   //new born kinematic particle
0079   KinematicParameters kPar(par);
0080   KinematicParametersError kEr(sCov);
0081   const MagneticField* field = input.front()->magneticField();
0082   KinematicState nState(kPar, kEr, ch, field);
0083
0084   //invalid previous particle and empty constraint:
0085   KinematicParticle* zp = nullptr;
0086   RefCountedKinematicParticle pPart = ReferenceCountingPointer<KinematicParticle>(zp);
0087
0088   float vChi = vtx.totalChiSquared();
0089   float vNdf = vtx.degreesOfFreedom();
0090   RefCountedKinematicParticle nPart = pFactory->particle(nState, vChi, vNdf, pPart);
0091
0092   //adding top particle to the tree
0093   RefCountedKinematicTree resTree = ReferenceCountingPointer<KinematicTree>(new KinematicTree());
0094   resTree->addParticle(pVrt, dVrt, nPart);
0095
0096   //making refitted kinematic particles and putting them to the tree
0097   std::vector<RefCountedKinematicParticle> rrP;
0098
0099   std::vector<RefCountedKinematicParticle>::const_iterator j;
0100   std::vector<RefCountedVertexTrack>::const_iterator i;
0101   for (j = input.begin(), i = refTracks.begin(); j != input.end() && i != refTracks.end(); ++j, ++i) {
0102     RefCountedLinearizedTrackState lT = (*i)->linearizedTrack();
0103     KinematicRefittedTrackState* rS = dynamic_cast<KinematicRefittedTrackState*>(&(*((*i)->refittedState())));
0104
0105     //   RefCountedRefittedTrackState rS = (*i)->refittedState();
0106     AlgebraicVector7 lPar = rS->kinematicParameters();
0107     KinematicParameters lkPar(lPar);
0108     AlgebraicSymMatrix77 lCov = rS->kinematicParametersCovariance();
0109     KinematicParametersError lkCov(lCov);
0110     TrackCharge lch = lT->charge();
0111     KinematicState nState(lkPar, lkCov, lch, field);
0112     RefCountedKinematicParticle nPart = (*j)->refittedParticle(nState, vChi, vNdf);
0113     rrP.push_back(nPart);
0114     if ((*j)->correspondingTree() != nullptr) {
0115       //here are the particles having trees after them
0116       KinematicTree* tree = (*j)->correspondingTree();
0117       tree->movePointerToTheTop();
0118       tree->replaceCurrentParticle(nPart);
0119       RefCountedKinematicVertex cdVertex = resTree->currentDecayVertex();
0120       resTree->addTree(cdVertex, tree);
0121     } else {
0122       //here are just particles fitted to this tree
0123       RefCountedKinematicVertex nV = kvFactory->vertex();
0124       resTree->addParticle(dVrt, nV, nPart);
0125     }
0126   }
0127   return resTree;
0128 }
0129
0130 //method returning the full covariance matrix
0131 //of new born kinematic particle
0132 AlgebraicMatrix FinalTreeBuilder::momentumPart(const CachingVertex<6>& vtx, const AlgebraicVector7& par) const {
0133   std::vector<RefCountedVertexTrack> refTracks = vtx.tracks();
0134   int size = refTracks.size();
0135   AlgebraicMatrix cov(7 + 4 * (size - 1), 7 + 4 * (size - 1));
0136   AlgebraicMatrix jac(7, 7 + 4 * (size - 1));
0137   jac(1, 1) = 1.;
0138   jac(2, 2) = 1.;
0139   jac(3, 3) = 1.;
0140
0141   double energy_total = sqrt(par(3) * par(3) + par(6) * par(6) + par(5) * par(5) + par(4) * par(4));
0142
0143   std::vector<RefCountedVertexTrack>::const_iterator rt_i;
0144   int i_int = 0;
0145   for (rt_i = refTracks.begin(); rt_i != refTracks.end(); rt_i++) {
0146     double a;
0147     AlgebraicVector6 param = (**rt_i).refittedState()->parameters();  // rho, theta, phi,tr_im, z_im, mass
0148     double rho = param[0];
0149     double theta = param[1];
0150     double phi = param[2];
0151     double mass = param[5];
0152
0153     if ((**rt_i).linearizedTrack()->charge() != 0) {
0154       a = -(**rt_i).refittedState()->freeTrajectoryState().parameters().magneticFieldInInverseGeV(vtx.position()).z() *
0155           (**rt_i).refittedState()->freeTrajectoryState().parameters().charge();
0156       if (a == 0.)
0157         throw cms::Exception("FinalTreeBuilder", "Field is 0");
0158     } else {
0159       a = 1;
0160     }
0161
0162     AlgebraicMatrix jc_el(4, 4, 0);
0163     jc_el(1, 1) = -a * cos(phi) / (rho * rho);    //dpx/d rho
0164     jc_el(2, 1) = -a * sin(phi) / (rho * rho);    //dpy/d rho
0165     jc_el(3, 1) = -a / (rho * rho * tan(theta));  //dpz/d rho
0166
0167     jc_el(3, 2) = -a / (rho * sin(theta) * sin(theta));  //dpz/d theta
0168
0169     jc_el(1, 3) = -a * sin(phi) / rho;  //dpx/d phi
0170     jc_el(2, 3) = a * cos(phi) / rho;   //dpy/d
0171
0172     //non-trival elements: mass correlations:
0173     double energy_local = sqrt(a * a / (rho * rho) * (1 + 1 / (tan(theta) * tan(theta))) + mass * mass);
0174
0175     jc_el(4, 4) = energy_total * mass / (par(6) * energy_local);  // dm/dm
0176
0177     jc_el(4, 1) = (-(energy_total / energy_local * a * a / (rho * rho * rho * sin(theta) * sin(theta))) +
0178                    par(3) * a / (rho * rho) * cos(phi) + par(4) * a / (rho * rho) * sin(phi) +
0179                    par(5) * a / (rho * rho * tan(theta))) /
0180                   par(6);  //dm / drho
0181
0182     jc_el(4, 2) = (-(energy_total / energy_local * a * a / (rho * rho * sin(theta) * sin(theta) * tan(theta))) +
0183                    par(5) * a / (rho * sin(theta) * sin(theta))) /
0184                   par(6);  //dm d theta
0185
0186     jc_el(4, 3) = (par(3) * sin(phi) - par(4) * cos(phi)) * a / (rho * par(6));  //dm/dphi
0187
0188     jac.sub(4, i_int * 4 + 4, jc_el);
0189
0190     //top left corner elements
0191     if (i_int == 0) {
0192       cov.sub(1, 1, asHepMatrix<7>((**rt_i).fullCovariance()));
0193     } else {
0194       //4-momentum corellatons: diagonal elements of the matrix
0195       AlgebraicMatrix fullCovMatrix(asHepMatrix<7>((**rt_i).fullCovariance()));
0196       AlgebraicMatrix m_m_cov = fullCovMatrix.sub(4, 7, 4, 7);
0197       AlgebraicMatrix x_p_cov = fullCovMatrix.sub(1, 3, 4, 7);
0198       AlgebraicMatrix p_x_cov = fullCovMatrix.sub(4, 7, 1, 3);
0199
0200       // cout << "Full covariance: \n"<< (**rt_i).fullCovariance()<<endl;
0201       // cout << "Full m_m_cov: "<< m_m_cov<<endl;
0202       //  cout<<"p_x_cov"<< p_x_cov<<endl;
0203       //  cout<<"x_p_cov"<< x_p_cov<<endl;
0204
0205       //putting everything to the joint covariance matrix:
0206       //diagonal momentum-momentum elements:
0207       cov.sub(i_int * 4 + 4, i_int * 4 + 4, m_m_cov);
0208
0209       //position momentum elements:
0210       cov.sub(1, i_int * 4 + 4, x_p_cov);
0211       cov.sub(i_int * 4 + 4, 1, p_x_cov);
0212
0213       //off diagonal elements: corellations
0214       // track momentum - track momentum
0215     }
0216     int j_int = 0;
0217     for (std::vector<RefCountedVertexTrack>::const_iterator rt_j = refTracks.begin(); rt_j != refTracks.end(); rt_j++) {
0218       if (i_int < j_int) {
0219         AlgebraicMatrix i_k_cov_m = asHepMatrix<4, 4>(vtx.tkToTkCovariance((*rt_i), (*rt_j)));
0220         //     cout<<"i_k_cov_m"<<i_k_cov_m <<endl;
0221         cov.sub(i_int * 4 + 4, j_int * 4 + 4, i_k_cov_m);
0222         cov.sub(j_int * 4 + 4, i_int * 4 + 4, i_k_cov_m.T());
0223       }
0224       j_int++;
0225     }
0226     i_int++;
0227   }
0228   // cout<<"jac"<<jac<<endl;
0229   // cout<<"cov"<<cov<<endl;
0230   //  cout << "final result new"<<jac*cov*jac.T()<<endl;
0231
0232   return jac * cov * jac.T();
0233 }