LArSoft  v10_04_05
Liquid Argon Software toolkit - https://larsoft.org/
Public Member Functions | Private Member Functions | Private Attributes | List of all members
genf::GFDaf Class Reference

#include "GFDaf.h"

Public Member Functions

 GFDaf ()
 Standard CTOR. Sets default values for annealing scheme and probablity cut. More...
 
 ~GFDaf ()
 
void processTrack (GFTrack *)
 Performs DAF fit on all track representations in a GFTrack. More...
 
void setBlowUpFactor (double f)
 Set the blowup factor (see blowUpCovs() ) More...
 
void setProbCut (double val)
 Set the probabilty cut for the weight calculation for the hits. Currently supported are the values 0.01 0.005, and 0.001. The corresponding chi2 cuts for different hits dimensionalities are hardcoded in the implementation because I did not yet figure out how to calculate them. Please feel very welcome to change the implementtion if you know how to do it. More...
 
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 temperatures. The maximum would ten tempertatures. More...
 

Private Member Functions

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. More...
 
void blowUpCovs (GFTrack *trk)
 This is needed to blow up the covariance matrix before a fitting pass. The method drops off-diagonal elements and blows up diagonal by blowUpFactor. More...
 
void invertMatrix (const TMatrixT< Double_t > &, TMatrixT< Double_t > &)
 invert a matrix. First argument is matrix to be inverted, second is return by ref. More...
 

Private Attributes

double fBlowUpFactor
 
std::vector< double > fBeta
 
std::map< int, double > chi2Cuts
 

Detailed Description

Definition at line 46 of file GFDaf.h.

Constructor & Destructor Documentation

genf::GFDaf::GFDaf ( )

Standard CTOR. Sets default values for annealing scheme and probablity cut.

Definition at line 42 of file GFDaf.cxx.

References setBetas(), and setProbCut().

42  : fBlowUpFactor(500.)
43 {
44  setProbCut(0.01);
45  setBetas(81, 8, 4, 1, 1, 1);
46 }
genf::GFDaf::setBetas
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
genf::GFDaf::fBlowUpFactor
double fBlowUpFactor
Definition: GFDaf.h:103
genf::GFDaf::setProbCut
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
genf::GFDaf::~GFDaf ( )

Definition at line 48 of file GFDaf.cxx.

49 {
50  ;
51 }

Member Function Documentation

void genf::GFDaf::blowUpCovs ( GFTrack *  trk)
private

This is needed to blow up the covariance matrix before a fitting pass. The method drops off-diagonal elements and blows up diagonal by blowUpFactor.

Definition at line 353 of file GFDaf.cxx.

References fBlowUpFactor, genf::GFAbsTrackRep::getCov(), genf::GFTrack::getNumReps(), genf::GFAbsTrackRep::getStatusFlag(), genf::GFTrack::getTrackRep(), and genf::GFAbsTrackRep::setCov().

Referenced by processTrack(), and setBlowUpFactor().

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 }
genf::GFDaf::fBlowUpFactor
double fBlowUpFactor
Definition: GFDaf.h:103
TMatrixT< Double_t > genf::GFDaf::calcGain ( const TMatrixT< Double_t > &  cov,
const TMatrixT< Double_t > &  HitCov,
const TMatrixT< Double_t > &  H,
const double &  p 
)
private

Calculate Kalman Gain.

Definition at line 394 of file GFDaf.cxx.

References invertMatrix().

Referenced by processTrack(), and setBlowUpFactor().

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 }
genf::GFDaf::invertMatrix
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 genf::GFDaf::invertMatrix ( const TMatrixT< Double_t > &  mat,
TMatrixT< Double_t > &  inv 
)
private

invert a matrix. First argument is matrix to be inverted, second is return by ref.

Definition at line 376 of file GFDaf.cxx.

References e, mat, and GFException::setFatal().

Referenced by calcGain(), processTrack(), and setBlowUpFactor().

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 }
mat
Float_t mat
Definition: plot.C:38
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:47
e
Float_t e
Definition: plot.C:35
void genf::GFDaf::processTrack ( GFTrack *  trk)

Performs DAF fit on all track representations in a GFTrack.

Definition at line 53 of file GFDaf.cxx.

References genf::GFTrack::addFailedHit(), blowUpCovs(), genf::GFBookkeeping::bookGFDetPlanes(), genf::GFBookkeeping::bookMatrices(), genf::GFBookkeeping::bookNumbers(), calcGain(), chi2Cuts, genf::GFTrack::clearBookkeeping(), genf::GFTrack::clearFailedHits(), COVEXC, E, genf::GFAbsTrackRep::extrapolate(), fBeta, genf::GFTrack::getBK(), genf::GFAbsTrackRep::getCov(), genf::GFBookkeeping::getDetPlane(), genf::GFAbsTrackRep::getDim(), genf::GFTrack::getHit(), genf::GFTrack::getHitsByPlane(), genf::GFBookkeeping::getMatrix(), genf::GFBookkeeping::getNumber(), genf::GFTrack::getNumHits(), genf::GFTrack::getNumReps(), genf::GFAbsTrackRep::getReferencePlane(), genf::GFAbsTrackRep::getState(), genf::GFAbsTrackRep::getStatusFlag(), genf::GFTrack::getTrackRep(), genf::GFBookkeeping::hitFailed(), hits(), GFException::info(), invertMatrix(), GFException::isFatal(), genf::GFAbsTrackRep::setData(), genf::GFBookkeeping::setDetPlane(), genf::GFBookkeeping::setMatrix(), genf::GFBookkeeping::setNhits(), genf::GFBookkeeping::setNumber(), genf::GFAbsTrackRep::setStatusFlag(), and GFException::what().

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 }
genf::GFDaf::invertMatrix
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
hits
void hits()
Definition: readHits.C:15
E
Float_t E
Definition: plot.C:20
genf::GFDaf::fBeta
std::vector< double > fBeta
Definition: GFDaf.h:104
GFException::isFatal
bool isFatal()
get fatal flag.
Definition: GFException.h:81
genf::GFDaf::calcGain
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
genf::GFDaf::chi2Cuts
std::map< int, double > chi2Cuts
Definition: GFDaf.h:105
GFException::info
void info()
print information in the exception object
Definition: GFException.cxx:61
GFException::what
const char * what() const noexcept override
standard error message handling for exceptions. use like "std::cerr << e.what();" ...
Definition: GFException.cxx:55
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:47
COVEXC
#define COVEXC
Definition: GFDaf.cxx:40
genf::GFDaf::blowUpCovs
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
void genf::GFDaf::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 temperatures. The maximum would ten tempertatures.

Definition at line 448 of file GFDaf.cxx.

References fBeta.

Referenced by GFDaf(), and setBlowUpFactor().

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 }
genf::GFDaf::fBeta
std::vector< double > fBeta
Definition: GFDaf.h:104
void genf::GFDaf::setBlowUpFactor ( double  f)
inline

Set the blowup factor (see blowUpCovs() )

Definition at line 61 of file GFDaf.h.

References blowUpCovs(), calcGain(), f, fBlowUpFactor, invertMatrix(), setBetas(), and setProbCut().

61 { fBlowUpFactor = f; }
genf::GFDaf::fBlowUpFactor
double fBlowUpFactor
Definition: GFDaf.h:103
f
TFile f
Definition: plotHisto.C:6
void genf::GFDaf::setProbCut ( double  val)

Set the probabilty cut for the weight calculation for the hits. Currently supported are the values 0.01 0.005, and 0.001. The corresponding chi2 cuts for different hits dimensionalities are hardcoded in the implementation because I did not yet figure out how to calculate them. Please feel very welcome to change the implementtion if you know how to do it.

Definition at line 416 of file GFDaf.cxx.

References chi2Cuts, E, GFException::setFatal(), and GFException::setNumbers().

Referenced by GFDaf(), and setBlowUpFactor().

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 }
GFException::setNumbers
GFException & setNumbers(std::string, const std::vector< double > &)
set list of numbers with description
Definition: GFException.cxx:40
E
Float_t E
Definition: plot.C:20
genf::GFDaf::chi2Cuts
std::map< int, double > chi2Cuts
Definition: GFDaf.h:105
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:47
GFException::setFatal
GFException & setFatal(bool b=true)
set fatal flag. if this is true, the fit stops for this current track repr.
Definition: GFException.h:75

Member Data Documentation

std::map<int, double> genf::GFDaf::chi2Cuts
private

Definition at line 105 of file GFDaf.h.

Referenced by processTrack(), and setProbCut().

std::vector<double> genf::GFDaf::fBeta
private

Definition at line 104 of file GFDaf.h.

Referenced by processTrack(), and setBetas().

double genf::GFDaf::fBlowUpFactor
private

Definition at line 103 of file GFDaf.h.

Referenced by blowUpCovs(), and setBlowUpFactor().

The documentation for this class was generated from the following files: