GsfTools/src/BasicMultiTrajectoryState.cc

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   // only accept planes!!
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   //int dim = tsos.front().localParameters().vector().num_row();
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 }