File indexing completed on 2024-04-06 12:29:07
0001 #include "RecoVertex/GaussianSumVertexFit/interface/MultiRefittedTS.h"
0002 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0003 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0004 #include "TrackingTools/GsfTools/interface/BasicMultiTrajectoryState.h"
0005 #include "RecoVertex/VertexPrimitives/interface/VertexException.h"
0006 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0007 #include "TrackingTools/TransientTrack/interface/TransientTrackFromFTSFactory.h"
0008 #include <cfloat>
0009 using namespace std;
0010
0011 MultiRefittedTS::MultiRefittedTS(const std::vector<RefCountedRefittedTrackState>& prtsComp,
0012 const Surface& referenceSurface)
0013 : theComponents(prtsComp), ftsAvailable(false), refSurface(&referenceSurface), surf(true) {}
0014
0015 MultiRefittedTS::MultiRefittedTS(const std::vector<RefCountedRefittedTrackState>& prtsComp,
0016 const GlobalPoint& referencePosition)
0017 : theComponents(prtsComp), ftsAvailable(false), refPosition(referencePosition), surf(false) {}
0018
0019
0020
0021
0022
0023 FreeTrajectoryState MultiRefittedTS::freeTrajectoryState() const {
0024 if (!ftsAvailable)
0025 computeFreeTrajectoryState();
0026 return fts;
0027 }
0028
0029 void MultiRefittedTS::computeFreeTrajectoryState() const {
0030 if (surf) {
0031 fts = *(trajectoryStateOnSurface(*refSurface).freeTrajectoryState());
0032 } else {
0033 double maxWeight = -1.;
0034 RTSvector::const_iterator maxIt;
0035 for (RTSvector::const_iterator it = theComponents.begin(); it != theComponents.end(); it++) {
0036 if ((**it).weight() > maxWeight) {
0037 maxWeight = (**it).weight();
0038 maxIt = it;
0039 }
0040 }
0041
0042 TransverseImpactPointExtrapolator tipe(&((**maxIt).freeTrajectoryState().parameters().magneticField()));
0043 TrajectoryStateOnSurface initialTSOS = tipe.extrapolate((**maxIt).freeTrajectoryState(), refPosition);
0044
0045 fts = *(trajectoryStateOnSurface(initialTSOS.surface()).freeTrajectoryState());
0046 }
0047 ftsAvailable = true;
0048 }
0049
0050
0051
0052
0053
0054
0055
0056 MultiRefittedTS::AlgebraicVectorN MultiRefittedTS::parameters() const {
0057 throw VertexException("MultiRefittedTS::freeTrajectoryState(): Don't know how to do that yet...");
0058 }
0059
0060
0061
0062
0063
0064 MultiRefittedTS::AlgebraicSymMatrixNN MultiRefittedTS::covariance() const {
0065 throw VertexException("MultiRefittedTS::freeTrajectoryState(): Don't know how to do that yet...");
0066 }
0067
0068
0069
0070
0071
0072 GlobalPoint MultiRefittedTS::position() const {
0073 throw VertexException("MultiRefittedTS::freeTrajectoryState(): Don't know how to do that yet...");
0074 }
0075
0076
0077
0078
0079
0080
0081 MultiRefittedTS::AlgebraicVectorM MultiRefittedTS::momentumVector() const {
0082 throw VertexException("MultiRefittedTS::freeTrajectoryState(): Don't know how to do that yet...");
0083 }
0084
0085 double MultiRefittedTS::weight() const {
0086 if (!totalWeightAvailable) {
0087 totalWeight = 0.;
0088 if (theComponents.empty()) {
0089 cout << "Asking for weight of empty MultiRefittedTS, returning zero!" << endl;
0090 }
0091 for (RTSvector::const_iterator it = theComponents.begin(); it != theComponents.end(); it++) {
0092 totalWeight += (**it).weight();
0093 }
0094 }
0095 return totalWeight;
0096 }
0097
0098 ReferenceCountingPointer<RefittedTrackState<5> > MultiRefittedTS::stateWithNewWeight(const double newWeight) const {
0099 if (weight() < DBL_MIN) {
0100 throw VertexException(
0101 "MultiRefittedTS::stateWithNewWeight(): Can not reweight multi-state with total weight < DBL_MIN");
0102 }
0103 double factor = newWeight / weight();
0104
0105 RTSvector reWeightedRTSC;
0106 reWeightedRTSC.reserve(theComponents.size());
0107
0108 for (RTSvector::const_iterator it = theComponents.begin(); it != theComponents.end(); it++) {
0109 reWeightedRTSC.push_back((**it).stateWithNewWeight((**it).weight() * factor));
0110 }
0111 if (surf) {
0112 return RefCountedRefittedTrackState(new MultiRefittedTS(reWeightedRTSC, *refSurface));
0113 } else {
0114 return RefCountedRefittedTrackState(new MultiRefittedTS(reWeightedRTSC, refPosition));
0115 }
0116 }
0117
0118
0119
0120
0121 TrajectoryStateOnSurface MultiRefittedTS::trajectoryStateOnSurface(const Surface& surface) const {
0122 vector<TrajectoryStateOnSurface> tsosComponents;
0123 tsosComponents.reserve(theComponents.size());
0124 for (RTSvector::const_iterator it = theComponents.begin(); it != theComponents.end(); it++) {
0125 tsosComponents.push_back((**it).trajectoryStateOnSurface(surface));
0126 }
0127 return TrajectoryStateOnSurface((BasicTrajectoryState*)new BasicMultiTrajectoryState(tsosComponents));
0128 }
0129
0130 TrajectoryStateOnSurface MultiRefittedTS::trajectoryStateOnSurface(const Surface& surface, const Propagator& propagator)
0131 const {
0132 vector<TrajectoryStateOnSurface> tsosComponents;
0133 tsosComponents.reserve(theComponents.size());
0134 for (RTSvector::const_iterator it = theComponents.begin(); it != theComponents.end(); it++) {
0135 tsosComponents.push_back((**it).trajectoryStateOnSurface(surface, propagator));
0136 }
0137 return TrajectoryStateOnSurface((BasicTrajectoryState*)new BasicMultiTrajectoryState(tsosComponents));
0138 }
0139
0140 reco::TransientTrack MultiRefittedTS::transientTrack() const {
0141 TransientTrackFromFTSFactory factory;
0142 return factory.build(freeTrajectoryState());
0143 }