File indexing completed on 2024-04-06 12:31:30
0001 #include "TrackingTools/GsfTools/interface/BasicMultiTrajectoryState.h"
0002
0003 #include "DataFormats/GeometrySurface/interface/BoundPlane.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 using namespace SurfaceSideDefinition;
0007
0008 BasicMultiTrajectoryState::BasicMultiTrajectoryState(const std::vector<TSOS>& tsvec)
0009 : BasicTrajectoryState(tsvec.front().surface()), theStates(tsvec) {
0010
0011 const BoundPlane* bp = dynamic_cast<const BoundPlane*>(&tsvec.begin()->surface());
0012 if UNLIKELY (bp == nullptr)
0013 throw cms::Exception("LogicError") << "MultiTrajectoryState constructed on cylinder";
0014
0015 for (auto i = tsvec.begin(); i != tsvec.end(); i++) {
0016 if UNLIKELY (!i->isValid()) {
0017 throw cms::Exception("LogicError") << "MultiTrajectoryState constructed with invalid state";
0018 }
0019 if UNLIKELY (i->hasError() != tsvec.front().hasError()) {
0020 throw cms::Exception("LogicError") << "MultiTrajectoryState mixes states with and without errors";
0021 }
0022 if UNLIKELY (&i->surface() != &tsvec.front().surface()) {
0023 throw cms::Exception("LogicError") << "MultiTrajectoryState mixes states with different surfaces";
0024 }
0025 if UNLIKELY (i->surfaceSide() != tsvec.front().surfaceSide()) {
0026 throw cms::Exception("LogicError") << "MultiTrajectoryState mixes states defined before and after material";
0027 }
0028 if UNLIKELY (i->localParameters().pzSign() * tsvec.front().localParameters().pzSign() < 0.) {
0029 throw cms::Exception("LogicError") << "MultiTrajectoryState mixes states with different signs of local p_z";
0030 }
0031 }
0032
0033 combine();
0034 }
0035
0036 void BasicMultiTrajectoryState::rescaleError(double factor) {
0037 if UNLIKELY (theStates.empty()) {
0038 edm::LogError("BasicMultiTrajectoryState") << "Trying to rescale errors of empty MultiTrajectoryState!";
0039 return;
0040 }
0041
0042 for (auto& is : theStates)
0043 is.rescaleError(factor);
0044 combine();
0045 }
0046
0047 void BasicMultiTrajectoryState::combine() {
0048 const std::vector<TrajectoryStateOnSurface>& tsos = theStates;
0049
0050 if UNLIKELY (tsos.empty()) {
0051 edm::LogError("MultiTrajectoryStateCombiner") << "Trying to collapse empty set of trajectory states!";
0052 return;
0053 }
0054
0055 double pzSign = tsos.front().localParameters().pzSign();
0056 for (std::vector<TrajectoryStateOnSurface>::const_iterator it = tsos.begin(); it != tsos.end(); it++) {
0057 if UNLIKELY (it->localParameters().pzSign() != pzSign) {
0058 edm::LogError("MultiTrajectoryStateCombiner")
0059 << "Trying to collapse trajectory states with different signs on p_z!";
0060 return;
0061 }
0062 }
0063
0064 if UNLIKELY (tsos.size() == 1) {
0065 BasicTrajectoryState::update(tsos.front().weight(),
0066 tsos.front().localParameters(),
0067 tsos.front().localError(),
0068 tsos.front().surface(),
0069 tsos.front().magneticField(),
0070 tsos.front().surfaceSide());
0071 return;
0072 }
0073
0074 double sumw = 0.;
0075
0076 AlgebraicVector5 mean;
0077 AlgebraicSymMatrix55 covarPart1, covarPart2, covtmp;
0078 for (auto it1 = tsos.begin(); it1 != tsos.end(); it1++) {
0079 auto weight = it1->weight();
0080 auto const& param = it1->localParameters().vector();
0081 sumw += weight;
0082 mean += weight * param;
0083 covarPart1 += weight * it1->localError().matrix();
0084 for (auto it2 = it1 + 1; it2 != tsos.end(); it2++) {
0085 AlgebraicVector5 diff = param - it2->localParameters().vector();
0086 ROOT::Math::AssignSym::Evaluate(covtmp, ROOT::Math::TensorProd(diff, diff));
0087 covarPart2 += (weight * it2->weight()) * covtmp;
0088 }
0089 }
0090 double sumwI = 1.0 / sumw;
0091 mean *= sumwI;
0092 covarPart1 *= sumwI;
0093 covarPart2 *= (sumwI * sumwI);
0094 AlgebraicSymMatrix55 covar = covarPart1 + covarPart2;
0095
0096 BasicTrajectoryState::update(sumw,
0097 LocalTrajectoryParameters(mean, pzSign),
0098 LocalTrajectoryError(covar),
0099 tsos.front().surface(),
0100 tsos.front().magneticField(),
0101 tsos.front().surfaceSide());
0102 }
0103
0104 void BasicMultiTrajectoryState::update(const LocalTrajectoryParameters& p,
0105 const Surface& aSurface,
0106 const MagneticField* field,
0107 const SurfaceSide side) {
0108 throw cms::Exception("LogicError",
0109 "BasicMultiTrajectoryState::update(LocalTrajectoryParameters, Surface, ...) called even if "
0110 "canUpdateLocalParameters() is false");
0111 }
0112
0113 void BasicMultiTrajectoryState::update(double weight,
0114 const LocalTrajectoryParameters& p,
0115 const LocalTrajectoryError& err,
0116 const Surface& aSurface,
0117 const MagneticField* field,
0118 const SurfaceSide side) {
0119 throw cms::Exception("LogicError",
0120 "BasicMultiTrajectoryState::update(LocalTrajectoryParameters, LocalTrajectoryError, ...) called "
0121 "even if canUpdateLocalParameters() is false");
0122 }