fitters/src/AbsKalmanFitter.cc

0001 /* Copyright 2013, Ludwig-Maximilians Universität München,
0002    Authors: Tobias Schlüter & Johannes Rauch
0003
0004    This file is part of GENFIT.
0005
0006    GENFIT is free software: you can redistribute it and/or modify
0007    it under the terms of the GNU Lesser General Public License as published
0008    by the Free Software Foundation, either version 3 of the License, or
0009    (at your option) any later version.
0010
0011    GENFIT is distributed in the hope that it will be useful,
0012    but WITHOUT ANY WARRANTY; without even the implied warranty of
0013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0014    GNU Lesser General Public License for more details.
0015
0016    You should have received a copy of the GNU Lesser General Public License
0017    along with GENFIT.  If not, see <http://www.gnu.org/licenses/>.
0018 */
0019
0020 /* This implements the simple Kalman fitter with no reference track
0021    that uses the stateSeed only until it forgets about it after the
0022    first few hits.  */
0023
0024 #include "Track.h"
0025 #include "TrackPoint.h"
0026 #include "Exception.h"
0027 #include "KalmanFitterInfo.h"
0028 #include "AbsMeasurement.h"
0029
0030 #include "AbsKalmanFitter.h"
0031 #include <Math/ProbFunc.h>
0032
0033
0034 //#define DEBUG
0035
0036 using namespace genfit;
0037
0038
0039
0040 void AbsKalmanFitter::getChiSquNdf(const Track* tr, const AbsTrackRep* rep,
0041     double& bChi2, double& fChi2,
0042     double& bNdf,  double& fNdf) const {
0043
0044   bChi2 = 0;
0045   fChi2 = 0;
0046   bNdf = -1. * rep->getDim();
0047   fNdf = -1. * rep->getDim();
0048
0049   const std::vector<TrackPoint*>& pointsWM = tr->getPointsWithMeasurement();
0050   for (std::vector<TrackPoint*>::const_iterator tpIter = pointsWM.begin(), endIter = pointsWM.end(); tpIter != endIter; ++tpIter) {
0051     if (! (*tpIter)->hasFitterInfo(rep))
0052       continue;
0053
0054     AbsFitterInfo* afi = (*tpIter)->getFitterInfo(rep);
0055     KalmanFitterInfo* fi = dynamic_cast<KalmanFitterInfo*>(afi);
0056     if (!fi) {
0057       Exception exc("AbsKalmanFitter::getChiSqu(): fitterInfo is not a KalmanFitterInfo", __LINE__,__FILE__);
0058       throw exc;
0059     }
0060
0061     KalmanFittedStateOnPlane* fup = fi->getForwardUpdate();
0062     KalmanFittedStateOnPlane* bup = fi->getBackwardUpdate();
0063
0064     if (fup == nullptr || bup == nullptr) {
0065       Exception exc("AbsKalmanFitter::getChiSqu(): fup == nullptr || bup == nullptr", __LINE__,__FILE__);
0066       throw exc;
0067     }
0068
0069     bChi2 += bup->getChiSquareIncrement();
0070     fChi2 += fup->getChiSquareIncrement();
0071
0072     bNdf += bup->getNdf();
0073     fNdf += fup->getNdf();
0074   }
0075
0076   if (bNdf < 0)
0077     bNdf = 0;
0078
0079   if (fNdf < 0)
0080     fNdf = 0;
0081 }
0082
0083
0084 double AbsKalmanFitter::getChiSqu(const Track* tr, const AbsTrackRep* rep, int direction) const {
0085   double bChi2(0), fChi2(0), bNdf(0), fNdf(0);
0086
0087   getChiSquNdf(tr, rep, bChi2, fChi2, bNdf, fNdf);
0088
0089   if (direction < 0)
0090     return bChi2;
0091   return fChi2;
0092 }
0093
0094 double AbsKalmanFitter::getNdf(const Track* tr, const AbsTrackRep* rep, int direction) const {
0095   double bChi2(0), fChi2(0), bNdf(0), fNdf(0);
0096
0097   getChiSquNdf(tr, rep, bChi2, fChi2, bNdf, fNdf);
0098
0099   if (direction < 0)
0100     return bNdf;
0101   return fNdf;
0102 }
0103
0104 double AbsKalmanFitter::getRedChiSqu(const Track* tr, const AbsTrackRep* rep, int direction) const {
0105   double bChi2(0), fChi2(0), bNdf(0), fNdf(0);
0106
0107   getChiSquNdf(tr, rep, bChi2, fChi2, bNdf, fNdf);
0108
0109   if (direction < 0)
0110     return bChi2/bNdf;
0111   return fChi2/fNdf;
0112 }
0113
0114 double AbsKalmanFitter::getPVal(const Track* tr, const AbsTrackRep* rep, int direction) const {
0115   double bChi2(0), fChi2(0), bNdf(0), fNdf(0);
0116
0117   getChiSquNdf(tr, rep, bChi2, fChi2, bNdf, fNdf);
0118
0119   if (direction < 0)
0120     return std::max(0.,ROOT::Math::chisquared_cdf_c(bChi2, bNdf));
0121   return std::max(0.,ROOT::Math::chisquared_cdf_c(fChi2, fNdf));
0122 }
0123
0124
0125 bool AbsKalmanFitter::isTrackPrepared(const Track* tr, const AbsTrackRep* rep) const {
0126   const std::vector<TrackPoint*>& points = tr->getPointsWithMeasurement();
0127
0128   if (points.size() == 0)
0129     return true;
0130
0131   for (std::vector<TrackPoint*>::const_iterator pIt = points.begin(), pEnd = points.end(); pIt != pEnd; ++pIt) {
0132     KalmanFitterInfo* fi = dynamic_cast<KalmanFitterInfo*>((*pIt)->getFitterInfo(rep));
0133
0134     if (!fi)
0135       continue;
0136
0137     if (!(fi->checkConsistency()))
0138       return false;
0139
0140     if (!(fi->hasReferenceState()))
0141       return false;
0142   }
0143
0144   return true;
0145 }
0146
0147 bool AbsKalmanFitter::isTrackFitted(const Track* tr, const AbsTrackRep* rep) const {
0148   if (! tr->getFitStatus(rep)->isFitted())
0149     return false;
0150
0151   // in depth check
0152
0153   const std::vector< TrackPoint* >& points = tr->getPointsWithMeasurement();
0154
0155   if (points.size() == 0)
0156     return true;
0157
0158   for (std::vector<TrackPoint*>::const_iterator pIt = points.begin(), pEnd = points.end(); pIt != pEnd; ++pIt) {
0159     KalmanFitterInfo* fi = dynamic_cast<KalmanFitterInfo*>((*pIt)->getFitterInfo(rep));
0160     if (!fi)
0161       return false;
0162
0163     if (!(fi->checkConsistency()))
0164       return false;
0165
0166     if (!(fi->hasForwardUpdate()))
0167       return false;
0168
0169     if (!(fi->hasBackwardUpdate()))
0170       return false;
0171   }
0172
0173   return true;
0174 }
0175
0176
0177 const std::vector<MeasurementOnPlane *> AbsKalmanFitter::getMeasurements(const KalmanFitterInfo* fi, const TrackPoint* tp, int direction) const {
0178
0179   switch (multipleMeasurementHandling_) {
0180     case weightedAverage :
0181     case unweightedAverage :
0182       return fi->getMeasurementsOnPlane();
0183
0184     case weightedClosestToReference :
0185     case unweightedClosestToReference :
0186     {
0187       if (!fi->hasReferenceState()) {
0188         Exception e("AbsKalmanFitter::getMeasurement: no ReferenceState.", __LINE__,__FILE__);
0189         e.setFatal();
0190         throw e;
0191       }
0192       std::vector<MeasurementOnPlane *> retVal;
0193       retVal.push_back(fi->getClosestMeasurementOnPlane(fi->getReferenceState()));
0194       return retVal;
0195     }
0196
0197     case weightedClosestToPrediction :
0198     case unweightedClosestToPrediction :
0199     {
0200       if (!fi->hasPrediction(direction)) {
0201         Exception e("AbsKalmanFitter::getMeasurement: no prediction.", __LINE__,__FILE__);
0202         e.setFatal();
0203         throw e;
0204       }
0205       std::vector<MeasurementOnPlane *> retVal;
0206       retVal.push_back(fi->getClosestMeasurementOnPlane(fi->getPrediction(direction)));
0207       return retVal;
0208     }
0209
0210
0211     case weightedClosestToReferenceWire :
0212     case unweightedClosestToReferenceWire :
0213     {
0214       if (tp->getNumRawMeasurements() == 1 && tp->getRawMeasurement()->isLeftRightMeasurement()) {
0215         if (!fi->hasReferenceState()) {
0216           Exception e("AbsKalmanFitter::getMeasurement: no ReferenceState.", __LINE__,__FILE__);
0217           e.setFatal();
0218           throw e;
0219         }
0220         std::vector<MeasurementOnPlane *> retVal;
0221         retVal.push_back(fi->getClosestMeasurementOnPlane(fi->getReferenceState()));
0222         return retVal;
0223       }
0224       else
0225         return fi->getMeasurementsOnPlane();
0226     }
0227
0228     case weightedClosestToPredictionWire :
0229     case unweightedClosestToPredictionWire :
0230     {
0231       if (tp->getNumRawMeasurements() == 1 && tp->getRawMeasurement()->isLeftRightMeasurement()) {
0232         if (!fi->hasPrediction(direction)) {
0233           Exception e("AbsKalmanFitter::getMeasurement: no prediction.", __LINE__,__FILE__);
0234           e.setFatal();
0235           throw e;
0236         }
0237         std::vector<MeasurementOnPlane *> retVal;
0238         retVal.push_back(fi->getClosestMeasurementOnPlane(fi->getPrediction(direction)));
0239         return retVal;
0240       }
0241       else
0242         return fi->getMeasurementsOnPlane();
0243     }
0244
0245
0246     default:
0247     {
0248       Exception e("AbsKalmanFitter::getMeasurement: choice not valid.", __LINE__,__FILE__);
0249       e.setFatal();
0250       throw e;
0251     }
0252   }
0253 }
0254
0255
0256 bool AbsKalmanFitter::canIgnoreWeights() const {
0257   switch (multipleMeasurementHandling_) {
0258     case unweightedAverage :
0259     case unweightedClosestToReference :
0260     case unweightedClosestToPrediction :
0261     case unweightedClosestToReferenceWire :
0262     case unweightedClosestToPredictionWire :
0263       return true;
0264
0265     case weightedAverage :
0266     case weightedClosestToReference :
0267     case weightedClosestToPrediction :
0268     case weightedClosestToReferenceWire :
0269     case weightedClosestToPredictionWire :
0270       return false;
0271
0272     default:
0273     {
0274       Exception e("AbsKalmanFitter::canIgnoreWeights: choice not valid.", __LINE__,__FILE__);
0275       e.setFatal();
0276       throw e;
0277     }
0278   }
0279 }