File indexing completed on 2024-04-06 12:31:30
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
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
0027
0028 GetComponents comps(tsos);
0029 auto const& components = comps();
0030 auto numb = components.size();
0031
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
0039 for (std::vector<TrajectoryStateOnSurface>::const_iterator ic = components.begin(); ic != components.end(); ++ic) {
0040
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
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
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
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
0074
0075
0076 GetComponents comps(tsos);
0077 auto const& components = comps();
0078 auto numb = components.size();
0079
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
0087 for (std::vector<TrajectoryStateOnSurface>::const_iterator ic = components.begin(); ic != components.end(); ++ic) {
0088
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
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
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
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
0122
0123 double qpMode(0);
0124 double dxdzMode(0);
0125 double dydzMode(0);
0126
0127
0128
0129 for (unsigned int iv = 0; iv < 3; ++iv) {
0130
0131 MultiGaussianState1D state1D = MultiGaussianStateTransform::multiState1D(tsos, iv);
0132 GaussianSumUtilities1D utils(state1D);
0133
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
0145 LocalVector localP(dxdzMode, dydzMode, 1.);
0146 localP *= tsos.localParameters().pzSign() / fabs(qpMode) / sqrt(dxdzMode * dxdzMode + dydzMode * dydzMode + 1.);
0147
0148 momentum = tsos.surface().toGlobal(localP);
0149 return true;
0150 }
0151
0152 bool momentumFromModeQP(TrajectoryStateOnSurface const& tsos, double& momentum) {
0153
0154
0155
0156 momentum = 0.;
0157 if (!tsos.isValid()) {
0158 edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
0159 return false;
0160 }
0161
0162
0163
0164 double qpMode(0);
0165
0166
0167
0168
0169 MultiGaussianState1D state1D = MultiGaussianStateTransform::multiState1D(tsos, 0);
0170 GaussianSumUtilities1D utils(state1D);
0171
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
0183
0184 momentum = 0.;
0185 if (!tsos.isValid()) {
0186 edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
0187 return false;
0188 }
0189
0190
0191
0192
0193 MultiGaussianState1D qpMultiState = MultiGaussianStateTransform::multiState1D(tsos, 0);
0194 std::vector<SingleGaussianState1D> states(qpMultiState.components());
0195
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
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
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
0226
0227 double xMode(0);
0228 double yMode(0);
0229
0230
0231
0232 for (unsigned int iv = 3; iv < 5; ++iv) {
0233
0234 MultiGaussianState1D state1D = MultiGaussianStateTransform::multiState1D(tsos, iv);
0235 GaussianSumUtilities1D utils(state1D);
0236
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
0246 LocalPoint localP(xMode, yMode, 0.);
0247
0248 position = tsos.surface().toGlobal(localP);
0249 return true;
0250 }
0251
0252 bool momentumFromModePPhiEta(TrajectoryStateOnSurface const& tsos, GlobalVector& momentum) {
0253
0254
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
0263
0264 GetComponents comps(tsos);
0265 auto const& components = comps();
0266 auto numb = components.size();
0267
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
0275 AlgebraicMatrix33 jacobian;
0276 AlgebraicSymMatrix33 covCart;
0277 AlgebraicSymMatrix33 covPPhiEta;
0278
0279 for (std::vector<TrajectoryStateOnSurface>::const_iterator ic = components.begin(); ic != components.end(); ++ic) {
0280
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
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
0300
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
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
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
0323 auto tanth2 = std::exp(-eta);
0324 auto pt = p * 2 * tanth2 / (1 + tanth2 * tanth2);
0325 auto pz = p * (1 - tanth2 * tanth2) / (1 + tanth2 * tanth2);
0326
0327 momentum = GlobalVector(pt * cos(phi), pt * sin(phi), pz);
0328 return true;
0329 }
0330
0331 int chargeFromMode(TrajectoryStateOnSurface const& tsos) {
0332
0333
0334
0335 if (!tsos.isValid()) {
0336 edm::LogInfo("multiTrajectoryStateMode") << "Cannot calculate mode from invalid TSOS";
0337 return 0;
0338 }
0339
0340
0341
0342 MultiGaussianState1D state1D = MultiGaussianStateTransform::multiState1D(tsos, 0);
0343 GaussianSumUtilities1D utils(state1D);
0344
0345 double result = utils.mode().mean();
0346 if (!utils.modeIsValid())
0347 result = utils.mean();
0348
0349 return result > 0. ? 1 : -1;
0350 }
0351
0352 }