KinematicFit/interface/ConstrainedTreeBuilderT.h

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