KinematicFit/src/ConstrainedTreeBuilder.cc

0001 #include "RecoVertex/KinematicFit/interface/ConstrainedTreeBuilder.h"
0002 #include "DataFormats/CLHEP/interface/Migration.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 ConstrainedTreeBuilder::ConstrainedTreeBuilder() {
0006   pFactory = new VirtualKinematicParticleFactory();
0007   vFactory = new KinematicVertexFactory();
0008 }
0009
0010 ConstrainedTreeBuilder::~ConstrainedTreeBuilder() {
0011   delete pFactory;
0012   delete vFactory;
0013 }
0014
0015 RefCountedKinematicTree ConstrainedTreeBuilder::buildTree(
0016     const std::vector<RefCountedKinematicParticle>& initialParticles,
0017     const std::vector<KinematicState>& finalStates,
0018     const RefCountedKinematicVertex vertex,
0019     const AlgebraicMatrix& fullCov) const {
0020   if (!vertex->vertexIsValid()) {
0021     LogDebug("ConstrainedTreeBuilder") << "Vertex is invalid\n";
0022     return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
0023   }
0024   AlgebraicVector3 vtx;
0025   vtx(0) = vertex->position().x();
0026   vtx(1) = vertex->position().y();
0027   vtx(2) = vertex->position().z();
0028   AlgebraicMatrix33 vertexCov = asSMatrix<3, 3>(fullCov.sub(1, 3, 1, 3));
0029
0030   // cout << fullCov<<endl;
0031   //  cout << "RecoVertex/ConstrainedTreeBuilder"<<vtx<<endl;
0032
0033   double ent = 0.;
0034   int charge = 0;
0035   AlgebraicVector7 virtualPartPar;
0036   virtualPartPar(0) = vertex->position().x();
0037   virtualPartPar(1) = vertex->position().y();
0038   virtualPartPar(2) = vertex->position().z();
0039
0040   //making refitted particles out of refitted states.
0041   //none of the operations above violates the order of particles
0042
0043   ROOT::Math::SMatrix<double, 7, 7, ROOT::Math::MatRepStd<double, 7, 7> > aMatrix;
0044   aMatrix(3, 3) = 1;
0045   aMatrix(4, 4) = 1;
0046   aMatrix(5, 5) = 1;
0047   aMatrix(6, 6) = 1;
0048   ROOT::Math::SMatrix<double, 7, 3, ROOT::Math::MatRepStd<double, 7, 3> > bMatrix;
0049   bMatrix(0, 0) = 1;
0050   bMatrix(1, 1) = 1;
0051   bMatrix(2, 2) = 1;
0052   AlgebraicMatrix77 trackParCov;
0053   ROOT::Math::SMatrix<double, 3, 7, ROOT::Math::MatRepStd<double, 3, 7> > vtxTrackCov;
0054   AlgebraicMatrix77 nCovariance;
0055
0056   std::vector<RefCountedKinematicParticle>::const_iterator i = initialParticles.begin();
0057   std::vector<KinematicState>::const_iterator iStates = finalStates.begin();
0058   std::vector<RefCountedKinematicParticle> rParticles;
0059   int n = 0;
0060   for (; i != initialParticles.end() && iStates != finalStates.end(); ++i, ++iStates) {
0061     AlgebraicVector7 p = iStates->kinematicParameters().vector();
0062     double a = -iStates->particleCharge() * iStates->magneticField()->inInverseGeV(iStates->globalPosition()).z();
0063
0064     aMatrix(4, 0) = -a;
0065     aMatrix(3, 1) = a;
0066     bMatrix(4, 0) = a;
0067     bMatrix(3, 1) = -a;
0068
0069     AlgebraicVector7 par = aMatrix * p + bMatrix * vtx;
0070
0071     trackParCov = asSMatrix<7, 7>(fullCov.sub(4 + n * 7, 10 + n * 7, 4 + n * 7, 10 + n * 7));
0072     vtxTrackCov = asSMatrix<3, 7>(fullCov.sub(1, 3, 4 + n * 7, 10 + n * 7));
0073     ;
0074     nCovariance = aMatrix * trackParCov * ROOT::Math::Transpose(aMatrix) +
0075                   aMatrix * ROOT::Math::Transpose(vtxTrackCov) * ROOT::Math::Transpose(bMatrix) +
0076                   bMatrix * vtxTrackCov * ROOT::Math::Transpose(aMatrix) +
0077                   bMatrix * vertexCov * ROOT::Math::Transpose(bMatrix);
0078
0079     KinematicState stateAtVertex(KinematicParameters(par),
0080                                  KinematicParametersError(AlgebraicSymMatrix77(nCovariance.LowerBlock())),
0081                                  iStates->particleCharge(),
0082                                  iStates->magneticField());
0083     rParticles.push_back((*i)->refittedParticle(stateAtVertex, vertex->chiSquared(), vertex->degreesOfFreedom()));
0084
0085     virtualPartPar(3) += par(3);
0086     virtualPartPar(4) += par(4);
0087     virtualPartPar(5) += par(5);
0088     ent += sqrt(par(6) * par(6) + par(3) * par(3) + par(4) * par(4) + par(5) * par(5));
0089     charge += iStates->particleCharge();
0090
0091     ++n;
0092   }
0093
0094   //total reconstructed mass
0095   double differ = ent * ent - (virtualPartPar(3) * virtualPartPar(3) + virtualPartPar(5) * virtualPartPar(5) +
0096                                virtualPartPar(4) * virtualPartPar(4));
0097   if (differ > 0.) {
0098     virtualPartPar(6) = sqrt(differ);
0099   } else {
0100     LogDebug("ConstrainedTreeBuilder") << "Fit failed: Current precision does not allow to calculate the mass\n";
0101     return ReferenceCountingPointer<KinematicTree>(new KinematicTree());
0102   }
0103
0104   // covariance matrix:
0105
0106   AlgebraicMatrix77 cov = asSMatrix<7, 7>(covarianceMatrix(rParticles, virtualPartPar, fullCov));
0107
0108   KinematicState nState(KinematicParameters(virtualPartPar),
0109                         KinematicParametersError(AlgebraicSymMatrix77(cov.LowerBlock())),
0110                         charge,
0111                         initialParticles[0]->magneticField());
0112
0113   //newborn kinematic particle
0114   float chi2 = vertex->chiSquared();
0115   float ndf = vertex->degreesOfFreedom();
0116   KinematicParticle* zp = nullptr;
0117   RefCountedKinematicParticle virtualParticle = pFactory->particle(nState, chi2, ndf, zp);
0118
0119   return buildTree(virtualParticle, vertex, rParticles);
0120 }
0121
0122 RefCountedKinematicTree ConstrainedTreeBuilder::buildTree(
0123     const RefCountedKinematicParticle virtualParticle,
0124     const RefCountedKinematicVertex vtx,
0125     const std::vector<RefCountedKinematicParticle>& particles) const {
0126   //making a resulting tree:
0127   RefCountedKinematicTree resTree = ReferenceCountingPointer<KinematicTree>(new KinematicTree());
0128
0129   //fake production vertex:
0130   RefCountedKinematicVertex fVertex = vFactory->vertex();
0131   resTree->addParticle(fVertex, vtx, virtualParticle);
0132
0133   //adding final state
0134   for (std::vector<RefCountedKinematicParticle>::const_iterator il = particles.begin(); il != particles.end(); il++) {
0135     if ((*il)->previousParticle()->correspondingTree() != nullptr) {
0136       KinematicTree* tree = (*il)->previousParticle()->correspondingTree();
0137       tree->movePointerToTheTop();
0138       tree->replaceCurrentParticle(*il);
0139       RefCountedKinematicVertex cdVertex = resTree->currentDecayVertex();
0140       resTree->addTree(cdVertex, tree);
0141     } else {
0142       RefCountedKinematicVertex ffVertex = vFactory->vertex();
0143       resTree->addParticle(vtx, ffVertex, *il);
0144     }
0145   }
0146   return resTree;
0147 }
0148
0149 AlgebraicMatrix ConstrainedTreeBuilder::covarianceMatrix(const std::vector<RefCountedKinematicParticle>& rPart,
0150                                                          const AlgebraicVector7& newPar,
0151                                                          const AlgebraicMatrix& fitCov) const {
0152   //constructing the full matrix using the simple fact
0153   //that we have never broken the order of tracks
0154   // during our fit.
0155
0156   int size = rPart.size();
0157   //global propagation to the vertex position
0158   //Jacobian is done for all the parameters together
0159   AlgebraicMatrix jac(3 + 7 * size, 3 + 7 * size, 0);
0160   jac(1, 1) = 1;
0161   jac(2, 2) = 1;
0162   jac(3, 3) = 1;
0163   int i_int = 0;
0164   for (std::vector<RefCountedKinematicParticle>::const_iterator i = rPart.begin(); i != rPart.end(); i++) {
0165     //vertex position related components of the matrix
0166     double a_i = -(*i)->currentState().particleCharge() *
0167                  (*i)->magneticField()->inInverseGeV((*i)->currentState().globalPosition()).z();
0168
0169     AlgebraicMatrix upper(3, 7, 0);
0170     AlgebraicMatrix diagonal(7, 7, 0);
0171     upper(1, 1) = 1;
0172     upper(2, 2) = 1;
0173     upper(3, 3) = 1;
0174     upper(2, 4) = -a_i;
0175     upper(1, 5) = a_i;
0176     jac.sub(1, 4 + i_int * 7, upper);
0177
0178     diagonal(4, 4) = 1;
0179     diagonal(5, 5) = 1;
0180     diagonal(6, 6) = 1;
0181     diagonal(7, 7) = 1;
0182     diagonal(2, 4) = a_i;
0183     diagonal(1, 5) = -a_i;
0184     jac.sub(4 + i_int * 7, 4 + i_int * 7, diagonal);
0185     i_int++;
0186   }
0187
0188   // jacobian is constructed in such a way, that
0189   // right operation for transformation will be
0190   // fitCov.similarityT(jac)
0191   // WARNING: normal similarity operation is
0192   // not valid in this case
0193   //now making reduced matrix:
0194   int vSize = rPart.size();
0195   AlgebraicSymMatrix fit_cov_sym(7 * vSize + 3, 0);
0196   for (int i = 1; i < 7 * vSize + 4; ++i) {
0197     for (int j = 1; j < 7 * vSize + 4; ++j) {
0198       if (i <= j)
0199         fit_cov_sym(i, j) = fitCov(i, j);
0200     }
0201   }
0202
0203   AlgebraicMatrix reduced(3 + 4 * size, 3 + 4 * size, 0);
0204   AlgebraicMatrix transform = fit_cov_sym.similarityT(jac);
0205
0206   //jacobian to add matrix components
0207   AlgebraicMatrix jac_t(7, 7 + 4 * (size - 1));
0208   jac_t(1, 1) = 1.;
0209   jac_t(2, 2) = 1.;
0210   jac_t(3, 3) = 1.;
0211
0212   //debug code:
0213   //CLHEP bugs: avoiding the
0214   // HepMatrix::sub method use
0215   AlgebraicMatrix reduced_tm(7, 7, 0);
0216   for (int i = 1; i < 8; ++i) {
0217     for (int j = 1; j < 8; ++j) {
0218       reduced_tm(i, j) = transform(i + 3, j + 3);
0219     }
0220   }
0221
0222   //left top corner
0223   // reduced.sub(1,1,transform.sub(4,10,4,10));
0224
0225   //debug code:
0226   //CLHEP bugs: avoiding the
0227   // HepMatrix::sub method use
0228   reduced.sub(1, 1, reduced_tm);
0229
0230   // cout<<"reduced matrix"<<reduced<<endl;
0231   int il_int = 0;
0232   double energy_global =
0233       sqrt(newPar(3) * newPar(3) + newPar(4) * newPar(4) + newPar(5) * newPar(5) + newPar(6) * newPar(6));
0234   for (std::vector<RefCountedKinematicParticle>::const_iterator rs = rPart.begin(); rs != rPart.end(); rs++) {
0235     //jacobian components:
0236     AlgebraicMatrix jc_el(4, 4, 0);
0237     jc_el(1, 1) = 1.;
0238     jc_el(2, 2) = 1.;
0239     jc_el(3, 3) = 1.;
0240
0241     //non-trival elements: mass correlations:
0242     AlgebraicVector7 l_Par = (*rs)->currentState().kinematicParameters().vector();
0243     double energy_local = sqrt(l_Par(6) * l_Par(6) + l_Par(3) * l_Par(3) + l_Par(4) * l_Par(4) + l_Par(5) * l_Par(5));
0244     jc_el(4, 4) = energy_global * l_Par(6) / (newPar(6) * energy_local);
0245     jc_el(4, 1) = ((energy_global * l_Par(3) / energy_local) - newPar(3)) / newPar(6);
0246     jc_el(4, 2) = ((energy_global * l_Par(4) / energy_local) - newPar(4)) / newPar(6);
0247     jc_el(4, 3) = ((energy_global * l_Par(5) / energy_local) - newPar(5)) / newPar(6);
0248     jac_t.sub(4, il_int * 4 + 4, jc_el);
0249     il_int++;
0250   }
0251   // cout<<"jac_t"<<jac_t<<endl;
0252   //debug code
0253   //CLHEP bugs workaround
0254   // cout<<"Transform matrix"<< transform<<endl;
0255
0256   for (int i = 1; i < size; i++) {
0257     //non-trival elements: mass correlations:
0258     //upper row and right column
0259     AlgebraicMatrix transform_sub1(3, 4, 0);
0260     AlgebraicMatrix transform_sub2(4, 3, 0);
0261     for (int l1 = 1; l1 < 4; ++l1) {
0262       for (int l2 = 1; l2 < 5; ++l2) {
0263         transform_sub1(l1, l2) = transform(3 + l1, 6 + 7 * i + l2);
0264       }
0265     }
0266
0267     for (int l1 = 1; l1 < 5; ++l1) {
0268       for (int l2 = 1; l2 < 4; ++l2) {
0269         transform_sub2(l1, l2) = transform(6 + 7 * i + l1, 3 + l2);
0270       }
0271     }
0272
0273     AlgebraicMatrix transform_sub3(4, 4, 0);
0274     for (int l1 = 1; l1 < 5; ++l1) {
0275       for (int l2 = 1; l2 < 5; ++l2) {
0276         transform_sub3(l1, l2) = transform(6 + 7 * i + l1, 6 + 7 * i + l2);
0277       }
0278     }
0279
0280     reduced.sub(1, 4 + 4 * i, transform_sub1);
0281     reduced.sub(4 + 4 * i, 1, transform_sub2);
0282
0283     //diagonal elements
0284     reduced.sub(4 + 4 * i, 4 + 4 * i, transform_sub3);
0285
0286     //off diagonal elements
0287     for (int j = 1; j < size; j++) {
0288       AlgebraicMatrix transform_sub4(4, 4, 0);
0289       AlgebraicMatrix transform_sub5(4, 4, 0);
0290       for (int l1 = 1; l1 < 5; ++l1) {
0291         for (int l2 = 1; l2 < 5; ++l2) {
0292           transform_sub4(l1, l2) = transform(6 + 7 * (i - 1) + l1, 6 + 7 * j + l2);
0293           transform_sub5(l1, l2) = transform(6 + 7 * j + l1, 6 + 7 * (i - 1) + l2);
0294         }
0295       }
0296       reduced.sub(4 + 4 * (i - 1), 4 + 4 * j, transform_sub4);
0297       reduced.sub(4 + 4 * j, 4 + 4 * (i - 1), transform_sub5);
0298     }
0299   }
0300
0301   return jac_t * reduced * jac_t.T();
0302 }