GsfTools/src/MultiTrajectoryStateMode.cc

0001 #include "TrackingTools/GsfTools/interface/MultiTrajectoryStateMode.h"
0002 #include "TrackingTools/GsfTools/interface/GetComponents.h"
0003
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0007 #include "DataFormats/GeometrySurface/interface/Surface.h"
0008 #include "TrackingTools/GsfTools/interface/MultiGaussianStateTransform.h"
0009 #include "TrackingTools/GsfTools/interface/MultiGaussianState1D.h"
0010 #include "TrackingTools/GsfTools/interface/GaussianSumUtilities1D.h"
0011
0012 #include <iostream>
0013
0014 namespace multiTrajectoryStateMode {
0015
0016   bool momentumFromModeCartesian(TrajectoryStateOnSurface const& tsos, GlobalVector& momentum) {
0017     //
0018     // clear result vector and check validity of the TSOS
0019     //
0020     momentum = GlobalVector(0., 0., 0.);
0021     if (!tsos.isValid()) {
0022       edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
0023       return false;
0024     }
0025     //
0026     // 1D mode computation for px, py and pz
0027     //
0028     GetComponents comps(tsos);
0029     auto const& components = comps();
0030     auto numb = components.size();
0031     // vectors of components in x, y and z
0032     std::vector<SingleGaussianState1D> pxStates;
0033     pxStates.reserve(numb);
0034     std::vector<SingleGaussianState1D> pyStates;
0035     pyStates.reserve(numb);
0036     std::vector<SingleGaussianState1D> pzStates;
0037     pzStates.reserve(numb);
0038     // iteration over components
0039     for (std::vector<TrajectoryStateOnSurface>::const_iterator ic = components.begin(); ic != components.end(); ++ic) {
0040       // extraction of parameters and variances
0041       GlobalVector mom(ic->globalMomentum());
0042       AlgebraicSymMatrix66 cov(ic->cartesianError().matrix());
0043       pxStates.push_back(SingleGaussianState1D(mom.x(), cov(3, 3), ic->weight()));
0044       pyStates.push_back(SingleGaussianState1D(mom.y(), cov(4, 4), ic->weight()));
0045       pzStates.push_back(SingleGaussianState1D(mom.z(), cov(5, 5), ic->weight()));
0046     }
0047     //
0048     // transformation in 1D multi-states and creation of utility classes
0049     //
0050     MultiGaussianState1D pxState(pxStates);
0051     MultiGaussianState1D pyState(pyStates);
0052     MultiGaussianState1D pzState(pzStates);
0053     GaussianSumUtilities1D pxUtils(pxState);
0054     GaussianSumUtilities1D pyUtils(pyState);
0055     GaussianSumUtilities1D pzUtils(pzState);
0056     //
0057     // cartesian momentum vector from modes
0058     //
0059     momentum = GlobalVector(pxUtils.mode().mean(), pyUtils.mode().mean(), pzUtils.mode().mean());
0060     return true;
0061   }
0062
0063   bool positionFromModeCartesian(TrajectoryStateOnSurface const& tsos, GlobalPoint& position) {
0064     //
0065     // clear result vector and check validity of the TSOS
0066     //
0067     position = GlobalPoint(0., 0., 0.);
0068     if (!tsos.isValid()) {
0069       edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
0070       return false;
0071     }
0072     //
0073     // 1D mode computation for x, y and z
0074     //
0075
0076     GetComponents comps(tsos);
0077     auto const& components = comps();
0078     auto numb = components.size();
0079     // vectors of components in x, y and z
0080     std::vector<SingleGaussianState1D> xStates;
0081     xStates.reserve(numb);
0082     std::vector<SingleGaussianState1D> yStates;
0083     yStates.reserve(numb);
0084     std::vector<SingleGaussianState1D> zStates;
0085     zStates.reserve(numb);
0086     // iteration over components
0087     for (std::vector<TrajectoryStateOnSurface>::const_iterator ic = components.begin(); ic != components.end(); ++ic) {
0088       // extraction of parameters and variances
0089       GlobalPoint pos(ic->globalPosition());
0090       AlgebraicSymMatrix66 cov(ic->cartesianError().matrix());
0091       xStates.push_back(SingleGaussianState1D(pos.x(), cov(0, 0), ic->weight()));
0092       yStates.push_back(SingleGaussianState1D(pos.y(), cov(1, 1), ic->weight()));
0093       zStates.push_back(SingleGaussianState1D(pos.z(), cov(2, 2), ic->weight()));
0094     }
0095     //
0096     // transformation in 1D multi-states and creation of utility classes
0097     //
0098     MultiGaussianState1D xState(xStates);
0099     MultiGaussianState1D yState(yStates);
0100     MultiGaussianState1D zState(zStates);
0101     GaussianSumUtilities1D xUtils(xState);
0102     GaussianSumUtilities1D yUtils(yState);
0103     GaussianSumUtilities1D zUtils(zState);
0104     //
0105     // cartesian position vector from modes
0106     //
0107     position = GlobalPoint(xUtils.mode().mean(), yUtils.mode().mean(), zUtils.mode().mean());
0108     return true;
0109   }
0110
0111   bool momentumFromModeLocal(TrajectoryStateOnSurface const& tsos, GlobalVector& momentum) {
0112     //
0113     // clear result vector and check validity of the TSOS
0114     //
0115     momentum = GlobalVector(0., 0., 0.);
0116     if (!tsos.isValid()) {
0117       edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
0118       return false;
0119     }
0120     //
0121     // mode computation for local co-ordinates q/p, dx/dz, dy/dz
0122     //
0123     double qpMode(0);
0124     double dxdzMode(0);
0125     double dydzMode(0);
0126     //
0127     // first 3 elements of local parameters = q/p, dx/dz, dy/dz
0128     //
0129     for (unsigned int iv = 0; iv < 3; ++iv) {
0130       // extraction of multi-state using helper class
0131       MultiGaussianState1D state1D = MultiGaussianStateTransform::multiState1D(tsos, iv);
0132       GaussianSumUtilities1D utils(state1D);
0133       // mode (in case of failure: mean)
0134       double result = utils.mode().mean();
0135       if (!utils.modeIsValid())
0136         result = utils.mean();
0137       if (iv == 0)
0138         qpMode = result;
0139       else if (iv == 1)
0140         dxdzMode = result;
0141       else
0142         dydzMode = result;
0143     }
0144     // local momentum vector from dx/dz, dy/dz and q/p + sign of local pz
0145     LocalVector localP(dxdzMode, dydzMode, 1.);
0146     localP *= tsos.localParameters().pzSign() / fabs(qpMode) / sqrt(dxdzMode * dxdzMode + dydzMode * dydzMode + 1.);
0147     // conversion to global coordinates
0148     momentum = tsos.surface().toGlobal(localP);
0149     return true;
0150   }
0151
0152   bool momentumFromModeQP(TrajectoryStateOnSurface const& tsos, double& momentum) {
0153     //
0154     // clear result vector and check validity of the TSOS
0155     //
0156     momentum = 0.;
0157     if (!tsos.isValid()) {
0158       edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
0159       return false;
0160     }
0161     //
0162     // mode computation for local co-ordinates q/p, dx/dz, dy/dz
0163     //
0164     double qpMode(0);
0165     //
0166     // first element of local parameters = q/p
0167     //
0168     // extraction of multi-state using helper class
0169     MultiGaussianState1D state1D = MultiGaussianStateTransform::multiState1D(tsos, 0);
0170     GaussianSumUtilities1D utils(state1D);
0171     // mode (in case of failure: mean)
0172     qpMode = utils.mode().mean();
0173     if (!utils.modeIsValid())
0174       qpMode = utils.mean();
0175
0176     momentum = 1. / fabs(qpMode);
0177     return true;
0178   }
0179
0180   bool momentumFromModeP(TrajectoryStateOnSurface const& tsos, double& momentum) {
0181     //
0182     // clear result vector and check validity of the TSOS
0183     //
0184     momentum = 0.;
0185     if (!tsos.isValid()) {
0186       edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
0187       return false;
0188     }
0189     //
0190     // first element of local parameters = q/p
0191     //
0192     // extraction of multi-state using helper class
0193     MultiGaussianState1D qpMultiState = MultiGaussianStateTransform::multiState1D(tsos, 0);
0194     std::vector<SingleGaussianState1D> states(qpMultiState.components());
0195     // transform from q/p to p
0196     for (unsigned int i = 0; i < states.size(); ++i) {
0197       SingleGaussianState1D& qpState = states[i];
0198       double wgt = qpState.weight();
0199       double qp = qpState.mean();
0200       double varQp = qpState.variance();
0201       double p = 1. / fabs(qp);
0202       double varP = p * p * p * p * varQp;
0203       states[i] = SingleGaussianState1D(p, varP, wgt);
0204     }
0205     MultiGaussianState1D pMultiState(states);
0206     GaussianSumUtilities1D utils(pMultiState);
0207     // mode (in case of failure: mean)
0208     momentum = utils.mode().mean();
0209     if (!utils.modeIsValid())
0210       momentum = utils.mean();
0211
0212     return true;
0213   }
0214
0215   bool positionFromModeLocal(TrajectoryStateOnSurface const& tsos, GlobalPoint& position) {
0216     //
0217     // clear result vector and check validity of the TSOS
0218     //
0219     position = GlobalPoint(0., 0., 0.);
0220     if (!tsos.isValid()) {
0221       edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
0222       return false;
0223     }
0224     //
0225     // mode computation for local co-ordinates x, y
0226     //
0227     double xMode(0);
0228     double yMode(0);
0229     //
0230     // last 2 elements of local parameters = x, y
0231     //
0232     for (unsigned int iv = 3; iv < 5; ++iv) {
0233       // extraction of multi-state using helper class
0234       MultiGaussianState1D state1D = MultiGaussianStateTransform::multiState1D(tsos, iv);
0235       GaussianSumUtilities1D utils(state1D);
0236       // mode (in case of failure: mean)
0237       double result = utils.mode().mean();
0238       if (!utils.modeIsValid())
0239         result = utils.mean();
0240       if (iv == 3)
0241         xMode = result;
0242       else
0243         yMode = result;
0244     }
0245     // local position vector from x, y
0246     LocalPoint localP(xMode, yMode, 0.);
0247     // conversion to global coordinates
0248     position = tsos.surface().toGlobal(localP);
0249     return true;
0250   }
0251
0252   bool momentumFromModePPhiEta(TrajectoryStateOnSurface const& tsos, GlobalVector& momentum) {
0253     //
0254     // clear result vector and check validity of the TSOS
0255     //
0256     momentum = GlobalVector(0., 0., 0.);
0257     if (!tsos.isValid()) {
0258       edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
0259       return false;
0260     }
0261     //
0262     // 1D mode computation for p, phi, eta
0263     //
0264     GetComponents comps(tsos);
0265     auto const& components = comps();
0266     auto numb = components.size();
0267     // vectors of components in p, phi and eta
0268     std::vector<SingleGaussianState1D> pStates;
0269     pStates.reserve(numb);
0270     std::vector<SingleGaussianState1D> phiStates;
0271     phiStates.reserve(numb);
0272     std::vector<SingleGaussianState1D> etaStates;
0273     etaStates.reserve(numb);
0274     // covariances in cartesian and p-phi-eta and jacobian
0275     AlgebraicMatrix33 jacobian;
0276     AlgebraicSymMatrix33 covCart;
0277     AlgebraicSymMatrix33 covPPhiEta;
0278     // iteration over components
0279     for (std::vector<TrajectoryStateOnSurface>::const_iterator ic = components.begin(); ic != components.end(); ++ic) {
0280       // parameters
0281       GlobalVector mom(ic->globalMomentum());
0282       auto px = mom.x();
0283       auto py = mom.y();
0284       auto pz = mom.z();
0285       auto op = 1. / mom.mag();
0286       auto opt2 = 1. / mom.perp2();
0287       auto phi = mom.phi();
0288       auto eta = mom.eta();
0289       // jacobian
0290       jacobian(0, 0) = px * op;
0291       jacobian(0, 1) = py * op;
0292       jacobian(0, 2) = pz * op;
0293       jacobian(1, 0) = py * opt2;
0294       jacobian(1, 1) = -px * opt2;
0295       jacobian(1, 2) = 0;
0296       jacobian(2, 0) = px * pz * opt2 * op;
0297       jacobian(2, 1) = py * pz * opt2 * op;
0298       jacobian(2, 2) = -op;
0299       // extraction of the momentum part from the 6x6 cartesian error matrix
0300       // and conversion to p-phi-eta
0301       covCart = ic->cartesianError().matrix().Sub<AlgebraicSymMatrix33>(3, 3);
0302       covPPhiEta = ROOT::Math::Similarity(jacobian, covCart);
0303       pStates.push_back(SingleGaussianState1D(1 / op, covPPhiEta(0, 0), ic->weight()));
0304       phiStates.push_back(SingleGaussianState1D(phi, covPPhiEta(1, 1), ic->weight()));
0305       etaStates.push_back(SingleGaussianState1D(eta, covPPhiEta(2, 2), ic->weight()));
0306     }
0307     //
0308     // transformation in 1D multi-states and creation of utility classes
0309     //
0310     MultiGaussianState1D pState(pStates);
0311     MultiGaussianState1D phiState(phiStates);
0312     MultiGaussianState1D etaState(etaStates);
0313     GaussianSumUtilities1D pUtils(pState);
0314     GaussianSumUtilities1D phiUtils(phiState);
0315     GaussianSumUtilities1D etaUtils(etaState);
0316     //
0317     // parameters from mode (in case of failure: mean)
0318     //
0319     auto p = pUtils.modeIsValid() ? pUtils.mode().mean() : pUtils.mean();
0320     auto phi = phiUtils.modeIsValid() ? phiUtils.mode().mean() : phiUtils.mean();
0321     auto eta = etaUtils.modeIsValid() ? etaUtils.mode().mean() : etaUtils.mean();
0322     //   double theta = 2*atan(exp(-eta));
0323     auto tanth2 = std::exp(-eta);
0324     auto pt = p * 2 * tanth2 / (1 + tanth2 * tanth2);             // p*sin(theta)
0325     auto pz = p * (1 - tanth2 * tanth2) / (1 + tanth2 * tanth2);  // p*cos(theta)
0326     // conversion to a cartesian momentum vector
0327     momentum = GlobalVector(pt * cos(phi), pt * sin(phi), pz);
0328     return true;
0329   }
0330
0331   int chargeFromMode(TrajectoryStateOnSurface const& tsos) {
0332     //
0333     // clear result vector and check validity of the TSOS
0334     //
0335     if (!tsos.isValid()) {
0336       edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
0337       return 0;
0338     }
0339     //
0340     // mode computation for local co-ordinates q/p
0341     // extraction of multi-state using helper class
0342     MultiGaussianState1D state1D = MultiGaussianStateTransform::multiState1D(tsos, 0);
0343     GaussianSumUtilities1D utils(state1D);
0344     // mode (in case of failure: mean)
0345     double result = utils.mode().mean();
0346     if (!utils.modeIsValid())
0347       result = utils.mean();
0348
0349     return result > 0. ? 1 : -1;
0350   }
0351
0352 }  // namespace multiTrajectoryStateMode