LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
GFDaf.cxx
Go to the documentation of this file.
1 /* Copyright 2008-2010, Technische Universitaet Muenchen,
2  Authors: Christian Hoeppner & Sebastian Neubert
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 /*
21  The DAF algorithm stems from the following sources:
22  R. Fruehwirth & A. Strandlie,
23  Computer Physics Communications 120 (199) 197-214
24  &
25  CERN thesis: Dissertation by Matthias Winkler
26  */
27 
28 #include "larreco/Genfit/GFDaf.h"
32 
33 #include "TMath.h"
34 #include <math.h>
35 
38 #include "larreco/Genfit/GFTrack.h"
39 
40 #define COVEXC "cov_is_zero"
41 
42 genf::GFDaf::GFDaf() : fBlowUpFactor(500.)
43 {
44  setProbCut(0.01);
45  setBetas(81, 8, 4, 1, 1, 1);
46 }
47 
49 {
50  ;
51 }
52 
54 {
55  trk->clearBookkeeping();
56 
57  unsigned int nreps = trk->getNumReps();
58  unsigned int nhits = trk->getNumHits();
59 
60  for (unsigned int i = 0; i < nreps; ++i) {
61  GFBookkeeping* bk = trk->getBK(i);
62  bk->setNhits(nhits);
63  bk->bookMatrices("fPredSt");
64  bk->bookMatrices("fPredCov");
65  bk->bookMatrices("bPredSt");
66  bk->bookMatrices("bPredCov");
67  bk->bookMatrices("H");
68  bk->bookMatrices("m");
69  bk->bookMatrices("V");
70  bk->bookGFDetPlanes("pl");
71  bk->bookMatrices("smooResid");
72  bk->bookNumbers("p", 1.);
73  }
74 
75  //plane grouping
76  std::vector<std::vector<int>*> planes;
77  trk->getHitsByPlane(planes);
78  int nPlanes = planes.size();
79 
80  trk->clearFailedHits();
81  for (unsigned int iBeta = 0; iBeta < fBeta.size(); iBeta++) {
82  //std::cout << "@@ beta " << fBeta.at(iBeta) << std::endl;
83  if (iBeta > 0) blowUpCovs(trk);
84  for (int ipl = 0; ipl < nPlanes; ++ipl) {
85  //std::cout << "@@ plane " << ipl << std::endl;
86  std::vector<GFAbsRecoHit*> hits;
87  unsigned int nhitsInPlane = planes.at(ipl)->size();
88  for (unsigned int ihitInPl = 0; ihitInPl < nhitsInPlane; ++ihitInPl) {
89  hits.push_back(trk->getHit(planes.at(ipl)->at(ihitInPl)));
90  }
91  //now hits contains pointers to all hits in the next plane
92  //for non-planar detectors this will always be just one hit
93  for (unsigned int irep = 0; irep < nreps; ++irep) {
94  GFAbsTrackRep* rep = trk->getTrackRep(irep);
95  if (rep->getStatusFlag() != 0) continue;
96  //std::cout << "@@ rep " << irep << std::endl;
97  GFBookkeeping* bk = trk->getBK(irep);
98  if (bk->hitFailed(planes.at(ipl)->at(0)) > 0) continue;
99  GFDetPlane pl;
100  // get prototypes for matrices
101  int repDim = rep->getDim();
102  TMatrixT<Double_t> state(repDim, 1);
103  TMatrixT<Double_t> cov(repDim, repDim);
104  ;
105 
106  if (ipl > 0 || (ipl == 0 && iBeta == 0)) {
107  // get the (virtual) detector plane
108  //std::cout << "forward ### iBeta, ipl, irep " << iBeta<< " " << ipl << " " << irep << std::endl;
109  // rep->getState().Print();
110  try {
111  pl = hits.at(0)->getDetPlane(rep);
112  //do the extrapolation
113  rep->extrapolate(pl, state, cov);
114 
115  if (cov[0][0] < 1.E-50) {
116  GFException exc(COVEXC, __LINE__, __FILE__);
117  throw exc;
118  }
119  }
120  catch (GFException& exc) {
121  std::cerr << exc.what();
122  exc.info();
123  for (unsigned int i = 0; i < planes.at(ipl)->size(); ++i)
124  trk->addFailedHit(irep, planes.at(ipl)->at(i));
125  if (exc.isFatal()) {
126  rep->setStatusFlag(1);
127  //abort();
128  }
129  continue; //go to next rep immediately
130  }
131  }
132  else {
133  pl = rep->getReferencePlane();
134  state = rep->getState();
135  cov = rep->getCov();
136  }
137  //the H matrix has to be identical for all hits in the plane
138  TMatrixT<Double_t> H = hits.at(0)->getHMatrix(rep);
139  bk->setMatrix("fPredSt", planes.at(ipl)->at(0), state);
140  bk->setMatrix("fPredCov", planes.at(ipl)->at(0), cov);
141 
142  TMatrixT<Double_t> covInv;
143  invertMatrix(cov, covInv);
144 
145  TMatrixT<Double_t> Htransp(H);
146  Htransp.T();
147 
148  bk->setMatrix("H", planes.at(ipl)->at(0), H);
149  bk->setDetPlane("pl", planes.at(ipl)->at(0), pl);
150 
151  TMatrixT<Double_t> stMod(state.GetNrows(), 1);
152 
153  double sumPk(0);
154  for (unsigned int ihit = 0; ihit < nhitsInPlane; ++ihit) {
155  double pki;
156  bk->getNumber("p", planes.at(ipl)->at(ihit), pki);
157  sumPk += pki;
158  }
159  TMatrixT<Double_t> V = hits.at(0)->getHitCov(pl);
160  TMatrixT<Double_t> Vinv;
161  invertMatrix(V, Vinv);
162  bk->setMatrix("V", planes.at(ipl)->at(0), V);
163  for (unsigned int ihit = 0; ihit < nhitsInPlane; ++ihit) {
164  //std::cout << "%%%% forward hit " << planes.at(ipl)->at(ihit) << std::endl;
165 
166  TMatrixT<Double_t> m = hits.at(ihit)->getHitCoord(pl);
167  bk->setMatrix("m", planes.at(ipl)->at(ihit), m);
168 
169  double pki;
170  bk->getNumber("p", planes.at(ipl)->at(ihit), pki);
171  TMatrixT<Double_t> Gain = calcGain(cov, V, H, sumPk);
172  //std::cout << "using weight " << pki << std::endl;
173  stMod += pki * Gain * (m - H * state);
174  }
175  state += stMod;
176  invertMatrix(covInv + sumPk * Htransp * Vinv * H, cov);
177  rep->setData(state, pl, &cov);
178 
179  } //loop over reps
180  } //loop over planes for forward fit
181  blowUpCovs(trk);
182  //now do the loop over the planes in backward direction
183  for (int ipl = nPlanes - 1; ipl >= 0; --ipl) {
184  //std::cout << "@bw@ plane " << ipl << std::endl;
185  std::vector<GFAbsRecoHit*> hits;
186  unsigned int nhitsInPlane = planes.at(ipl)->size();
187  for (unsigned int ihitInPl = 0; ihitInPl < nhitsInPlane; ++ihitInPl) {
188  hits.push_back(trk->getHit(planes.at(ipl)->at(ihitInPl)));
189  }
190  //now hits contains pointers to all hits in the next plane
191  //for non-planar detectors this will always be just one hit
192  for (unsigned int irep = 0; irep < nreps; ++irep) {
193  GFAbsTrackRep* rep = trk->getTrackRep(irep);
194  if (rep->getStatusFlag() != 0) continue;
195  //std::cout << "@bw@ rep " << irep << std::endl;
196  GFBookkeeping* bk = trk->getBK(irep);
197  if (bk->hitFailed(planes.at(ipl)->at(0)) > 0) continue;
198  // get prototypes for matrices
199  int repDim = rep->getDim();
200  TMatrixT<Double_t> state(repDim, 1);
201  TMatrixT<Double_t> cov(repDim, repDim);
202  ;
203 
204  TMatrixT<Double_t> H;
205  bk->getMatrix("H", planes.at(ipl)->at(0), H);
206  TMatrixT<Double_t> Htransp(H);
207  Htransp.T();
208  GFDetPlane pl;
209  bk->getDetPlane("pl", planes.at(ipl)->at(0), pl);
210  //std::cout << "backward ### iBeta, ipl, irep " << iBeta<< " " << ipl << " " << irep << std::endl;
211  if (ipl < (nPlanes - 1)) {
212  try {
213  //do the extrapolation
214  rep->extrapolate(pl, state, cov);
215  if (cov[0][0] < 1.E-50) {
216  GFException exc(COVEXC, __LINE__, __FILE__);
217  throw exc;
218  }
219  }
220  catch (GFException& exc) {
221  //std::cout << __FILE__ << " " << __LINE__ << std::endl;
222  std::cerr << exc.what();
223  exc.info();
224  //TAG
225  for (unsigned int i = 0; i < planes.at(ipl)->size(); ++i)
226  trk->addFailedHit(irep, planes.at(ipl)->at(i));
227  if (exc.isFatal()) { rep->setStatusFlag(1); }
228  continue; //go to next rep immediately
229  }
230  }
231  else {
232  state = rep->getState();
233  cov = rep->getCov();
234  }
235  TMatrixT<Double_t> covInv;
236  invertMatrix(cov, covInv);
237 
238  bk->setMatrix("bPredSt", planes.at(ipl)->at(0), state);
239  bk->setMatrix("bPredCov", planes.at(ipl)->at(0), cov);
240 
241  TMatrixT<Double_t> stMod(state.GetNrows(), 1);
242 
243  double sumPk(0);
244  for (unsigned int ihit = 0; ihit < nhitsInPlane; ++ihit) {
245  double pki;
246  bk->getNumber("p", planes.at(ipl)->at(ihit), pki);
247  sumPk += pki;
248  }
249 
250  TMatrixT<Double_t> V;
251  bk->getMatrix("V", planes.at(ipl)->at(0), V);
252  TMatrixT<Double_t> Vinv;
253  invertMatrix(V, Vinv);
254 
255  for (unsigned int ihit = 0; ihit < nhitsInPlane; ++ihit) {
256  //std::cout << "%%bw%% hit " << planes.at(ipl)->at(ihit) << std::endl;
257  TMatrixT<Double_t> m;
258  bk->getMatrix("m", planes.at(ipl)->at(ihit), m);
259 
260  double pki;
261  bk->getNumber("p", planes.at(ipl)->at(ihit), pki);
262  TMatrixT<Double_t> Gain = calcGain(cov, V, H, sumPk);
263  //std::cout << "using pki " << pki << std::endl;
264 
265  stMod += pki * Gain * (m - H * state);
266  }
267  state += stMod;
268  invertMatrix(covInv + sumPk * Htransp * Vinv * H, cov);
269 
270  rep->setData(state, pl, &cov);
271  } //done with loop over reps
272  } //loop over planes for backward
273  //calculate smoothed states and covs
274  TMatrixT<Double_t> fSt, fCov, bSt, bCov, smooSt, smooCov, fCovInv, bCovInv, m, H, smooResid;
275  for (int ipl = 0; ipl < nPlanes; ++ipl) {
276  for (unsigned int irep = 0; irep < nreps; ++irep) {
277  if (trk->getTrackRep(irep)->getStatusFlag() != 0) continue;
278  GFBookkeeping* bk = trk->getBK(irep);
279  if (bk->hitFailed(planes.at(ipl)->at(0)) > 0) continue;
280  bk->getMatrix("fPredSt", planes.at(ipl)->at(0), fSt);
281  bk->getMatrix("bPredSt", planes.at(ipl)->at(0), bSt);
282  bk->getMatrix("fPredCov", planes.at(ipl)->at(0), fCov);
283  bk->getMatrix("bPredCov", planes.at(ipl)->at(0), bCov);
284  fCovInv.ResizeTo(fCov);
285  bCovInv.ResizeTo(bCov);
286  invertMatrix(fCov, fCovInv);
287  invertMatrix(bCov, bCovInv);
288  smooCov.ResizeTo(fCov);
289  invertMatrix(fCovInv + bCovInv, smooCov);
290  smooSt.ResizeTo(fSt.GetNrows(), 1);
291  smooSt = smooCov * (fCovInv * fSt + bCovInv * bSt);
292  //std::cout << "#@#@#@#@#@# smooResid ipl, irep " << ipl << " " << irep << std::endl;
293  bk->getMatrix("H", planes.at(ipl)->at(0), H);
294  for (unsigned int ihit = 0; ihit < planes.at(ipl)->size(); ++ihit) {
295  //TAG
296  bk->getMatrix("m", planes.at(ipl)->at(ihit), m);
297  smooResid.ResizeTo(m.GetNrows(), 1);
298  smooResid = m - H * smooSt;
299  //std::cout << "%%%% hit " << planes.at(ipl)->at(ihit) << std::endl;
300  //smooResid.Print();
301  bk->setMatrix("smooResid", planes.at(ipl)->at(ihit), smooResid);
302  }
303  }
304  } //end loop over planes for smoothed state calculation
305 
306  //calculate the new weights
307  if (iBeta != fBeta.size() - 1) { //dont do it for the last beta, ause fitting will stop
308  for (unsigned int irep = 0; irep < nreps; ++irep) {
309  GFBookkeeping* bk = trk->getBK(irep);
310  for (int ipl = 0; ipl < nPlanes; ++ipl) {
311  if (trk->getTrackRep(irep)->getStatusFlag() != 0) continue;
312  if (bk->hitFailed(planes.at(ipl)->at(0)) > 0) continue;
313  double sumPhi = 0.;
314  double sumPhiCut = 0.;
315  std::vector<double> phi;
316  TMatrixT<Double_t> V;
317  bk->getMatrix("V", planes.at(ipl)->at(0), V);
318  V *= fBeta.at(iBeta);
319  TMatrixT<Double_t> Vinv;
320  invertMatrix(V, Vinv);
321  for (unsigned int ihit = 0; ihit < planes.at(ipl)->size(); ++ihit) {
322  //TMatrixT<Double_t> smooResid;
323  bk->getMatrix("smooResid", planes.at(ipl)->at(ihit), smooResid);
324  TMatrixT<Double_t> smooResidTransp(smooResid);
325  smooResidTransp.T();
326  int dimV = V.GetNrows();
327  double detV = V.Determinant();
328  TMatrixT<Double_t> expArg = smooResidTransp * Vinv * smooResid;
329  if (expArg.GetNrows() != 1)
330  throw std::runtime_error("genf::GFDaf::processTrack(): wrong matrix dimensions!");
331  double thisPhi =
332  1. / (pow(2. * TMath::Pi(), dimV / 2.) * sqrt(detV)) * exp(-0.5 * expArg[0][0]);
333  phi.push_back(thisPhi);
334  sumPhi += thisPhi;
335 
336  double cutVal = chi2Cuts[dimV];
337  if (cutVal <= 1.E-6)
338  throw std::runtime_error("genf::GFDaf::processTrack(): chi2 cut too low");
339  sumPhiCut += 1. / (2 * TMath::Pi() * sqrt(detV)) * exp(-0.5 * cutVal / fBeta.at(iBeta));
340  }
341  for (unsigned int ihit = 0; ihit < planes.at(ipl)->size(); ++ihit) {
342  bk->setNumber("p", planes.at(ipl)->at(ihit), phi.at(ihit) / (sumPhi + sumPhiCut));
343  }
344  }
345  }
346  }
347 
348  } //loop over betas
349 
350  return;
351 }
352 
354 {
355  int nreps = trk->getNumReps();
356  for (int irep = 0; irep < nreps; ++irep) {
357  GFAbsTrackRep* arep = trk->getTrackRep(irep);
358  //dont do it for already compromsied reps, since they wont be fitted anyway
359  if (arep->getStatusFlag() == 0) {
360  TMatrixT<Double_t> cov = arep->getCov();
361  for (int i = 0; i < cov.GetNrows(); ++i) {
362  for (int j = 0; j < cov.GetNcols(); ++j) {
363  if (i != j) { //off diagonal
364  cov[i][j] = 0.;
365  }
366  else { //diagonal
367  cov[i][j] = cov[i][j] * fBlowUpFactor;
368  }
369  }
370  }
371  arep->setCov(cov);
372  }
373  }
374 }
375 
376 void genf::GFDaf::invertMatrix(const TMatrixT<Double_t>& mat, TMatrixT<Double_t>& inv)
377 {
378  inv.ResizeTo(mat);
379  inv = (mat);
380  double det = 0;
381  inv.Invert(&det);
382  if (TMath::IsNaN(det)) {
383  GFException e("Daf Gain: det of matrix is nan", __LINE__, __FILE__);
384  e.setFatal();
385  throw e;
386  }
387  if (det == 0) {
388  GFException e("cannot invert matrix in Daf - det=0", __LINE__, __FILE__);
389  e.setFatal();
390  throw e;
391  }
392 }
393 
394 TMatrixT<Double_t> genf::GFDaf::calcGain(const TMatrixT<Double_t>& C,
395  const TMatrixT<Double_t>& V,
396  const TMatrixT<Double_t>& H,
397  const double& p)
398 {
399 
400  //get C^-1
401  TMatrixT<Double_t> Cinv;
402  invertMatrix(C, Cinv);
403  TMatrixT<Double_t> Vinv;
404  invertMatrix(V, Vinv);
405  //get H^T
406  TMatrixT<Double_t> Htransp(H);
407  Htransp.T();
408 
409  TMatrixT<Double_t> covsum = Cinv + (p * Htransp * Vinv * H);
410  TMatrixT<Double_t> covsumInv;
411  invertMatrix(covsum, covsumInv);
412 
413  return (covsumInv * Htransp * Vinv);
414 }
415 
416 void genf::GFDaf::setProbCut(double val)
417 {
418 
419  if (fabs(val - 0.01) < 1.E-10) {
420  chi2Cuts[1] = 6.63;
421  chi2Cuts[2] = 9.21;
422  chi2Cuts[3] = 11.34;
423  chi2Cuts[4] = 13.23;
424  chi2Cuts[5] = 15.09;
425  }
426  else if (fabs(val - 0.005) < 1.E-10) {
427  chi2Cuts[1] = 7.88;
428  chi2Cuts[2] = 10.60;
429  chi2Cuts[3] = 12.84;
430  chi2Cuts[4] = 14.86;
431  chi2Cuts[5] = 16.75;
432  }
433  else if (fabs(val - 0.001) < 1.E-10) {
434  chi2Cuts[1] = 10.83;
435  chi2Cuts[2] = 13.82;
436  chi2Cuts[3] = 16.27;
437  chi2Cuts[4] = 18.47;
438  chi2Cuts[5] = 20.51;
439  }
440  else {
441  throw GFException(
442  "genf::GFTrackCand::GFDafsetProbCut(): value not supported", __LINE__, __FILE__)
443  .setFatal()
444  .setNumbers("value", {val});
445  }
446 }
447 
448 void genf::GFDaf::setBetas(double b1,
449  double b2,
450  double b3 /* = -1. */,
451  double b4 /* = -1. */,
452  double b5 /* = -1. */,
453  double b6 /* = -1. */,
454  double b7 /* = -1. */,
455  double b8 /* = -1. */,
456  double b9 /* = -1. */,
457  double b10 /* = -1. */
458 )
459 {
460 
461  /* // asserting version
462  assert(b1>0);fBeta.push_back(b1);
463  assert(b2>0 && b2<=b1);fBeta.push_back(b2);
464  if(b3>=0.) {
465  assert(b3<=b2);fBeta.push_back(b3);
466  if(b4>=0.) {
467  assert(b4<=b3);fBeta.push_back(b4);
468  if(b5>=0.) {
469  assert(b5<=b4);fBeta.push_back(b5);
470  if(b6>=0.) {
471  assert(b6<=b5);fBeta.push_back(b6);
472  if(b7>=0.) {
473  assert(b7<=b6);fBeta.push_back(b7);
474  if(b8>=0.) {
475  assert(b8<=b7);fBeta.push_back(b8);
476  if(b9>=0.) {
477  assert(b9<=b8);fBeta.push_back(b9);
478  if(b10>=0.) {
479  assert(b10<=b9);fBeta.push_back(b10);
480  }
481  }
482  }
483  }
484  }
485  }
486  }
487  }
488  */
489  if (b1 <= 0) throw std::runtime_error("genf::GFDaf::setBetas(): b1 <= 0");
490  fBeta.push_back(b1);
491 
492  if (b2 <= 0 || b2 > b1) throw std::runtime_error("genf::GFDaf::setBetas(): b2 in wrong range");
493  fBeta.push_back(b2);
494 
495  if (b3 < 0.) return;
496  if (b3 > b2) throw std::runtime_error("genf::GFDaf::setBetas(): b3 > b2");
497  fBeta.push_back(b3);
498 
499  if (b4 < 0.) return;
500  if (b4 > b3) throw std::runtime_error("genf::GFDaf::setBetas(): b4 > b3");
501  fBeta.push_back(b4);
502 
503  if (b5 < 0.) return;
504  if (b5 > b4) throw std::runtime_error("genf::GFDaf::setBetas(): b5 > b4");
505  fBeta.push_back(b5);
506 
507  if (b6 < 0.) return;
508  if (b6 > b5) throw std::runtime_error("genf::GFDaf::setBetas(): b6 > b5");
509  fBeta.push_back(b6);
510 
511  if (b7 < 0.) return;
512  if (b7 > b6) throw std::runtime_error("genf::GFDaf::setBetas(): b7 > b6");
513  fBeta.push_back(b7);
514 
515  if (b8 < 0.) return;
516  if (b8 > b7) throw std::runtime_error("genf::GFDaf::setBetas(): b8 > b7");
517  fBeta.push_back(b8);
518 
519  if (b9 < 0.) return;
520  if (b9 > b8) throw std::runtime_error("genf::GFDaf::setBetas(): b9 > b8");
521  fBeta.push_back(b9);
522 
523  if (b10 < 0.) return;
524  if (b10 > b9) throw std::runtime_error("genf::GFDaf::setBetas(): b10 > b9");
525  fBeta.push_back(b10);
526 }
void setDetPlane(std::string key, unsigned int index, const GFDetPlane &pl)
GFException & setNumbers(std::string, const std::vector< double > &)
set list of numbers with description
Definition: GFException.cxx:40
GFDaf()
Standard CTOR. Sets default values for annealing scheme and probablity cut.
Definition: GFDaf.cxx:42
const GFDetPlane & getReferencePlane() const
void invertMatrix(const TMatrixT< Double_t > &, TMatrixT< Double_t > &)
invert a matrix. First argument is matrix to be inverted, second is return by ref.
Definition: GFDaf.cxx:376
void setBetas(double b1, double b2, double b3=-1., double b4=-1., double b5=-1., double b6=-1., double b7=-1., double b8=-1., double b9=-1., double b10=-1.)
Configure the annealing scheme. In the current implementation you need to provide at least two temper...
Definition: GFDaf.cxx:448
GFAbsTrackRep * getTrackRep(int id) const
Accessor for track representations.
Definition: GFTrack.h:193
void setNumber(std::string key, unsigned int index, const double &num)
const TMatrixT< Double_t > & getState() const
void setNhits(int n)
Definition: GFBookkeeping.h:47
void processTrack(GFTrack *)
Performs DAF fit on all track representations in a GFTrack.
Definition: GFDaf.cxx:53
double fBlowUpFactor
Definition: GFDaf.h:103
Base Class for genfit track representations. Defines interface for track parameterizations.
Definition: GFAbsTrackRep.h:81
void setCov(const TMatrixT< Double_t > &aCov)
void setMatrix(std::string key, unsigned int index, const TMatrixT< Double_t > &mat)
void hits()
Definition: readHits.C:15
void bookGFDetPlanes(std::string key)
Float_t E
Definition: plot.C:20
unsigned int getNumReps() const
Get number of track represenatations.
Definition: GFTrack.h:200
std::vector< double > fBeta
Definition: GFDaf.h:104
virtual void setData(const TMatrixT< Double_t > &st, const GFDetPlane &pl, const TMatrixT< Double_t > *cov=NULL)
void addFailedHit(unsigned int irep, unsigned int id)
Definition: GFTrack.h:294
bool isFatal()
get fatal flag.
Definition: GFException.h:81
GFAbsRecoHit * getHit(int id) const
Definition: GFTrack.h:161
TMatrixT< Double_t > calcGain(const TMatrixT< Double_t > &cov, const TMatrixT< Double_t > &HitCov, const TMatrixT< Double_t > &H, const double &p)
Calculate Kalman Gain.
Definition: GFDaf.cxx:394
virtual double extrapolate(const GFDetPlane &plane, TMatrixT< Double_t > &statePred)
returns the tracklength spanned in this extrapolation
Float_t mat
Definition: plot.C:38
unsigned int hitFailed(unsigned int)
const TMatrixT< Double_t > & getCov() const
GFBookkeeping * getBK(int index=-1)
get GFBookKeeping object for particular track rep (default is cardinal rep)
Definition: GFTrack.h:334
std::map< int, double > chi2Cuts
Definition: GFDaf.h:105
void setProbCut(double val)
Set the probabilty cut for the weight calculation for the hits. Currently supported are the values 0...
Definition: GFDaf.cxx:416
void info()
print information in the exception object
Definition: GFException.cxx:61
const char * what() const noexcept override
standard error message handling for exceptions. use like "std::cerr << e.what();" ...
Definition: GFException.cxx:55
bool getNumber(std::string key, unsigned int index, double &num) const
unsigned int getNumHits() const
Definition: GFTrack.h:163
void clearFailedHits()
Definition: GFTrack.h:432
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:47
void setStatusFlag(int _val)
void bookMatrices(std::string key)
bool getMatrix(std::string key, unsigned int index, TMatrixT< Double_t > &mat) const
void bookNumbers(std::string key, double val=0.)
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:75
void getHitsByPlane(std::vector< std::vector< int > * > &retVal)
Definition: GFTrack.cxx:215
#define COVEXC
Definition: GFDaf.cxx:40
unsigned int getDim() const
returns dimension of state vector
Float_t e
Definition: plot.C:35
void blowUpCovs(GFTrack *trk)
This is needed to blow up the covariance matrix before a fitting pass. The method drops off-diagonal ...
Definition: GFDaf.cxx:353
bool getDetPlane(std::string key, unsigned int index, GFDetPlane &pl) const
void clearBookkeeping()
Definition: GFTrack.h:425