VertexTools/src/InvariantMassFromVertex.cc

0001 #include "RecoVertex/VertexTools/interface/InvariantMassFromVertex.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 InvariantMassFromVertex::LorentzVector InvariantMassFromVertex::p4(const CachingVertex<5>& vertex,
0005                                                                    const double mass) const {
0006   return p4(vertex, std::vector<double>(vertex.tracks().size(), mass));
0007 }
0008
0009 InvariantMassFromVertex::LorentzVector InvariantMassFromVertex::p4(const CachingVertex<5>& vertex,
0010                                                                    const std::vector<double>& masses) const {
0011   LorentzVector totalP4;
0012
0013   // Check that tkToTkCovarianceIsAvailable
0014   if (!vertex.tkToTkCovarianceIsAvailable()) {
0015     LogDebug("InvariantMassFromVertex") << "Fit failed: vertex has not been smoothed\n";
0016     return totalP4;
0017   }
0018
0019   if (vertex.tracks().size() != masses.size()) {
0020     LogDebug("InvariantMassFromVertex") << "Vector of masses does not have the same size as tracks in vertex\n";
0021     return totalP4;
0022   }
0023
0024   std::vector<RefCountedVertexTrack> refTracks = vertex.tracks();
0025   std::vector<RefCountedVertexTrack>::const_iterator i_s = refTracks.begin();
0026   std::vector<double>::const_iterator i_m = masses.begin();
0027
0028   for (; i_s != refTracks.end() && i_m != masses.end(); ++i_s, ++i_m) {
0029     GlobalVector momentum = (**i_s).refittedState()->freeTrajectoryState().momentum();
0030     totalP4 += LorentzVector(momentum.x(), momentum.y(), momentum.z(), *i_m);
0031   }
0032   return totalP4;
0033 }
0034
0035 GlobalVector InvariantMassFromVertex::momentum(const CachingVertex<5>& vertex) const {
0036   GlobalVector momentum_;
0037
0038   // Check that tkToTkCovarianceIsAvailable
0039   if (!vertex.tkToTkCovarianceIsAvailable()) {
0040     LogDebug("InvariantMassFromVertex") << "Fit failed: vertex has not been smoothed\n";
0041     return momentum_;
0042   }
0043
0044   std::vector<RefCountedVertexTrack> refTracks = vertex.tracks();
0045   std::vector<RefCountedVertexTrack>::const_iterator i_s = refTracks.begin();
0046
0047   for (; i_s != refTracks.end(); ++i_s) {
0048     momentum_ += (**i_s).refittedState()->freeTrajectoryState().momentum();
0049   }
0050   return momentum_;
0051 }
0052
0053 Measurement1D InvariantMassFromVertex::invariantMass(const CachingVertex<5>& vertex, const double mass) const {
0054   return invariantMass(vertex, std::vector<double>(vertex.tracks().size(), mass));
0055 }
0056
0057 Measurement1D InvariantMassFromVertex::invariantMass(const CachingVertex<5>& vertex,
0058                                                      const std::vector<double>& masses) const {
0059   // Check that tkToTkCovarianceIsAvailable
0060   if (!vertex.tkToTkCovarianceIsAvailable()) {
0061     LogDebug("InvariantMassFromVertex") << "Fit failed: vertex has not been smoothed\n";
0062     return Measurement1D(0., 0.);
0063   }
0064   if (vertex.tracks().size() != masses.size()) {
0065     LogDebug("InvariantMassFromVertex") << "Vector of masses does not have the same size as tracks in vertex\n";
0066     return Measurement1D(0., 0.);
0067   }
0068
0069   LorentzVector totalP4 = p4(vertex, masses);
0070   double u = uncertainty(totalP4, vertex, masses);
0071   // std::cout << u<<std::endl;
0072   return Measurement1D(totalP4.M(), u);
0073 }
0074
0075 //method returning the full covariance matrix
0076 //of new born kinematic particle
0077 double InvariantMassFromVertex::uncertainty(const LorentzVector& totalP4,
0078                                             const CachingVertex<5>& vertex,
0079                                             const std::vector<double>& masses) const {
0080   std::vector<RefCountedVertexTrack> refTracks = vertex.tracks();
0081   int size = refTracks.size();
0082   AlgebraicMatrix cov(3 * size, 3 * size);
0083   AlgebraicMatrix jac(1, 3 * size);
0084
0085   double energy_total = totalP4.E();
0086
0087   std::vector<RefCountedVertexTrack>::const_iterator rt_i = refTracks.begin();
0088   std::vector<double>::const_iterator i_m = masses.begin();
0089
0090   int i_int = 0;
0091   for (; rt_i != refTracks.end() && i_m != masses.end(); ++rt_i, ++i_m) {
0092     double a;
0093     AlgebraicVector5 param = (**rt_i).refittedState()->parameters();  // rho, theta, phi,tr_im, z_im
0094     double rho = param[0];
0095     double theta = param[1];
0096     double phi = param[2];
0097     double mass = *i_m;
0098
0099     if ((**rt_i).linearizedTrack()->charge() != 0) {
0100       a = -(**rt_i).refittedState()->freeTrajectoryState().parameters().magneticFieldInInverseGeV(vertex.position()).z() *
0101           (**rt_i).refittedState()->freeTrajectoryState().parameters().charge();
0102       if (a == 0.)
0103         throw cms::Exception("InvariantMassFromVertex", "Field is 0");
0104     } else {
0105       a = 1;
0106     }
0107
0108     double energy_local = sqrt(a * a / (rho * rho) * (1 + 1 / (tan(theta) * tan(theta))) + mass * mass);
0109
0110     jac(1, i_int * 3 + 1) = (-(energy_total / energy_local * a * a / (rho * rho * rho * sin(theta) * sin(theta))) +
0111                              totalP4.X() * a / (rho * rho) * cos(phi) + totalP4.Y() * a / (rho * rho) * sin(phi) +
0112                              totalP4.Z() * a / (rho * rho * tan(theta))) /
0113                             totalP4.M();  //dm / drho
0114
0115     jac(1, i_int * 3 + 2) =
0116         (-(energy_total / energy_local * a * a / (rho * rho * sin(theta) * sin(theta) * tan(theta))) +
0117          totalP4.Z() * a / (rho * sin(theta) * sin(theta))) /
0118         totalP4.M();  //dm d theta
0119
0120     jac(1, i_int * 3 + 3) = (totalP4.X() * sin(phi) - totalP4.Y() * cos(phi)) * a / (rho * totalP4.M());  //dm/dphi
0121
0122     // momentum corellatons: diagonal elements of the matrix
0123     cov.sub(i_int * 3 + 1, i_int * 3 + 1, asHepMatrix<6>((**rt_i).fullCovariance()).sub(4, 6));
0124
0125     //off diagonal elements: track momentum - track momentum corellations
0126
0127     int j_int = 0;
0128     for (std::vector<RefCountedVertexTrack>::const_iterator rt_j = refTracks.begin(); rt_j != refTracks.end(); rt_j++) {
0129       if (i_int < j_int) {
0130         AlgebraicMatrix i_k_cov_m = asHepMatrix<3, 3>(vertex.tkToTkCovariance((*rt_i), (*rt_j)));
0131         cov.sub(i_int * 3 + 1, j_int * 3 + 1, i_k_cov_m);
0132         cov.sub(j_int * 3 + 1, i_int * 3 + 1, i_k_cov_m.T());
0133       }
0134       j_int++;
0135     }
0136     i_int++;
0137   }
0138   //   std::cout<<"jac"<<jac<<std::endl;
0139   //   std::cout<<"cov"<<cov<<std::endl;
0140   //   std::cout << "final result"<<(jac*cov*jac.T())<<std::endl;
0141
0142   return sqrt((jac * cov * jac.T())(1, 1));
0143 }