Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AbsKalmanFitter.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file AbsKalmanFitter.cc
1 /* Copyright 2013, Ludwig-Maximilians Universität München,
2  Authors: Tobias Schlüter & Johannes Rauch
3 
4  This file is part of GENFIT.
5 
6  GENFIT is free software: you can redistribute it and/or modify
7  it under the terms of the GNU Lesser General Public License as published
8  by the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  GENFIT is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU Lesser General Public License for more details.
15 
16  You should have received a copy of the GNU Lesser General Public License
17  along with GENFIT. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 /* This implements the simple Kalman fitter with no reference track
21  that uses the stateSeed only until it forgets about it after the
22  first few hits. */
23 
24 #include "Track.h"
25 #include "TrackPoint.h"
26 #include "Exception.h"
27 #include "KalmanFitterInfo.h"
28 #include "AbsMeasurement.h"
29 
30 #include "AbsKalmanFitter.h"
31 #include <Math/ProbFunc.h>
32 
33 
34 //#define DEBUG
35 
36 using namespace genfit;
37 
38 
39 
40 void AbsKalmanFitter::getChiSquNdf(const Track* tr, const AbsTrackRep* rep,
41  double& bChi2, double& fChi2,
42  double& bNdf, double& fNdf) const {
43 
44  bChi2 = 0;
45  fChi2 = 0;
46  bNdf = -1. * rep->getDim();
47  fNdf = -1. * rep->getDim();
48 
49  const std::vector<TrackPoint*>& pointsWM = tr->getPointsWithMeasurement();
50  for (std::vector<TrackPoint*>::const_iterator tpIter = pointsWM.begin(), endIter = pointsWM.end(); tpIter != endIter; ++tpIter) {
51  if (! (*tpIter)->hasFitterInfo(rep))
52  continue;
53 
54  AbsFitterInfo* afi = (*tpIter)->getFitterInfo(rep);
55  KalmanFitterInfo* fi = dynamic_cast<KalmanFitterInfo*>(afi);
56  if (!fi) {
57  Exception exc("AbsKalmanFitter::getChiSqu(): fitterInfo is not a KalmanFitterInfo", __LINE__,__FILE__);
58  throw exc;
59  }
60 
63 
64  if (fup == nullptr || bup == nullptr) {
65  Exception exc("AbsKalmanFitter::getChiSqu(): fup == nullptr || bup == nullptr", __LINE__,__FILE__);
66  throw exc;
67  }
68 
69  bChi2 += bup->getChiSquareIncrement();
70  fChi2 += fup->getChiSquareIncrement();
71 
72  bNdf += bup->getNdf();
73  fNdf += fup->getNdf();
74  }
75 
76  if (bNdf < 0)
77  bNdf = 0;
78 
79  if (fNdf < 0)
80  fNdf = 0;
81 }
82 
83 
84 double AbsKalmanFitter::getChiSqu(const Track* tr, const AbsTrackRep* rep, int direction) const {
85  double bChi2(0), fChi2(0), bNdf(0), fNdf(0);
86 
87  getChiSquNdf(tr, rep, bChi2, fChi2, bNdf, fNdf);
88 
89  if (direction < 0)
90  return bChi2;
91  return fChi2;
92 }
93 
94 double AbsKalmanFitter::getNdf(const Track* tr, const AbsTrackRep* rep, int direction) const {
95  double bChi2(0), fChi2(0), bNdf(0), fNdf(0);
96 
97  getChiSquNdf(tr, rep, bChi2, fChi2, bNdf, fNdf);
98 
99  if (direction < 0)
100  return bNdf;
101  return fNdf;
102 }
103 
104 double AbsKalmanFitter::getRedChiSqu(const Track* tr, const AbsTrackRep* rep, int direction) const {
105  double bChi2(0), fChi2(0), bNdf(0), fNdf(0);
106 
107  getChiSquNdf(tr, rep, bChi2, fChi2, bNdf, fNdf);
108 
109  if (direction < 0)
110  return bChi2/bNdf;
111  return fChi2/fNdf;
112 }
113 
114 double AbsKalmanFitter::getPVal(const Track* tr, const AbsTrackRep* rep, int direction) const {
115  double bChi2(0), fChi2(0), bNdf(0), fNdf(0);
116 
117  getChiSquNdf(tr, rep, bChi2, fChi2, bNdf, fNdf);
118 
119  if (direction < 0)
120  return std::max(0.,ROOT::Math::chisquared_cdf_c(bChi2, bNdf));
121  return std::max(0.,ROOT::Math::chisquared_cdf_c(fChi2, fNdf));
122 }
123 
124 
125 bool AbsKalmanFitter::isTrackPrepared(const Track* tr, const AbsTrackRep* rep) const {
126  const std::vector<TrackPoint*>& points = tr->getPointsWithMeasurement();
127 
128  if (points.size() == 0)
129  return true;
130 
131  for (std::vector<TrackPoint*>::const_iterator pIt = points.begin(), pEnd = points.end(); pIt != pEnd; ++pIt) {
132  KalmanFitterInfo* fi = dynamic_cast<KalmanFitterInfo*>((*pIt)->getFitterInfo(rep));
133 
134  if (!fi)
135  continue;
136 
137  if (!(fi->checkConsistency()))
138  return false;
139 
140  if (!(fi->hasReferenceState()))
141  return false;
142  }
143 
144  return true;
145 }
146 
147 bool AbsKalmanFitter::isTrackFitted(const Track* tr, const AbsTrackRep* rep) const {
148  if (! tr->getFitStatus(rep)->isFitted())
149  return false;
150 
151  // in depth check
152 
153  const std::vector< TrackPoint* >& points = tr->getPointsWithMeasurement();
154 
155  if (points.size() == 0)
156  return true;
157 
158  for (std::vector<TrackPoint*>::const_iterator pIt = points.begin(), pEnd = points.end(); pIt != pEnd; ++pIt) {
159  KalmanFitterInfo* fi = dynamic_cast<KalmanFitterInfo*>((*pIt)->getFitterInfo(rep));
160  if (!fi)
161  return false;
162 
163  if (!(fi->checkConsistency()))
164  return false;
165 
166  if (!(fi->hasForwardUpdate()))
167  return false;
168 
169  if (!(fi->hasBackwardUpdate()))
170  return false;
171  }
172 
173  return true;
174 }
175 
176 
177 const std::vector<MeasurementOnPlane *> AbsKalmanFitter::getMeasurements(const KalmanFitterInfo* fi, const TrackPoint* tp, int direction) const {
178 
180  case weightedAverage :
181  case unweightedAverage :
182  return fi->getMeasurementsOnPlane();
183 
186  {
187  if (!fi->hasReferenceState()) {
188  Exception e("AbsKalmanFitter::getMeasurement: no ReferenceState.", __LINE__,__FILE__);
189  e.setFatal();
190  throw e;
191  }
192  std::vector<MeasurementOnPlane *> retVal;
193  retVal.push_back(fi->getClosestMeasurementOnPlane(fi->getReferenceState()));
194  return retVal;
195  }
196 
199  {
200  if (!fi->hasPrediction(direction)) {
201  Exception e("AbsKalmanFitter::getMeasurement: no prediction.", __LINE__,__FILE__);
202  e.setFatal();
203  throw e;
204  }
205  std::vector<MeasurementOnPlane *> retVal;
206  retVal.push_back(fi->getClosestMeasurementOnPlane(fi->getPrediction(direction)));
207  return retVal;
208  }
209 
210 
213  {
215  if (!fi->hasReferenceState()) {
216  Exception e("AbsKalmanFitter::getMeasurement: no ReferenceState.", __LINE__,__FILE__);
217  e.setFatal();
218  throw e;
219  }
220  std::vector<MeasurementOnPlane *> retVal;
221  retVal.push_back(fi->getClosestMeasurementOnPlane(fi->getReferenceState()));
222  return retVal;
223  }
224  else
225  return fi->getMeasurementsOnPlane();
226  }
227 
230  {
232  if (!fi->hasPrediction(direction)) {
233  Exception e("AbsKalmanFitter::getMeasurement: no prediction.", __LINE__,__FILE__);
234  e.setFatal();
235  throw e;
236  }
237  std::vector<MeasurementOnPlane *> retVal;
238  retVal.push_back(fi->getClosestMeasurementOnPlane(fi->getPrediction(direction)));
239  return retVal;
240  }
241  else
242  return fi->getMeasurementsOnPlane();
243  }
244 
245 
246  default:
247  {
248  Exception e("AbsKalmanFitter::getMeasurement: choice not valid.", __LINE__,__FILE__);
249  e.setFatal();
250  throw e;
251  }
252  }
253 }
254 
255 
258  case unweightedAverage :
263  return true;
264 
265  case weightedAverage :
270  return false;
271 
272  default:
273  {
274  Exception e("AbsKalmanFitter::canIgnoreWeights: choice not valid.", __LINE__,__FILE__);
275  e.setFatal();
276  throw e;
277  }
278  }
279 }