GaussianSumVertexFit/src/MultiRefittedTS.cc

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    * Transformation into a FreeTrajectoryState
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    * Vector containing the refitted track parameters. <br>
0052    * These are (signed transverse curvature, theta, phi,
0053    *  (signed) transverse , longitudinal impact parameter)
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    * The covariance matrix
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    * Position at which the momentum is defined.
0070    */
0071
0072 GlobalPoint MultiRefittedTS::position() const {
0073   throw VertexException("MultiRefittedTS::freeTrajectoryState(): Don't know how to do that yet...");
0074 }
0075
0076 /**
0077    * Vector containing the parameters describing the momentum as the vertex.
0078    * These are (signed transverse curvature, theta, phi)
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    * Transformation into a TSOS at a given surface
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 {  //fixme... is the propagation done correctly? Is there a gsf propagator?
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 }