LArSoft  v09_90_00
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:449
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:417
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 354 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().

355 {
356  int nreps = trk->getNumReps();
357  for (int irep = 0; irep < nreps; ++irep) {
358  GFAbsTrackRep* arep = trk->getTrackRep(irep);
359  //dont do it for already compromsied reps, since they wont be fitted anyway
360  if (arep->getStatusFlag() == 0) {
361  TMatrixT<Double_t> cov = arep->getCov();
362  for (int i = 0; i < cov.GetNrows(); ++i) {
363  for (int j = 0; j < cov.GetNcols(); ++j) {
364  if (i != j) { //off diagonal
365  cov[i][j] = 0.;
366  }
367  else { //diagonal
368  cov[i][j] = cov[i][j] * fBlowUpFactor;
369  }
370  }
371  }
372  arep->setCov(cov);
373  }
374  }
375 }
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 395 of file GFDaf.cxx.

References invertMatrix().

Referenced by processTrack(), and setBlowUpFactor().

399 {
400 
401  //get C^-1
402  TMatrixT<Double_t> Cinv;
403  invertMatrix(C, Cinv);
404  TMatrixT<Double_t> Vinv;
405  invertMatrix(V, Vinv);
406  //get H^T
407  TMatrixT<Double_t> Htransp(H);
408  Htransp.T();
409 
410  TMatrixT<Double_t> covsum = Cinv + (p * Htransp * Vinv * H);
411  TMatrixT<Double_t> covsumInv;
412  invertMatrix(covsum, covsumInv);
413 
414  return (covsumInv * Htransp * Vinv);
415 }
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:377
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 377 of file GFDaf.cxx.

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

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

378 {
379  inv.ResizeTo(mat);
380  inv = (mat);
381  double det = 0;
382  inv.Invert(&det);
383  if (TMath::IsNaN(det)) {
384  GFException e("Daf Gain: det of matrix is nan", __LINE__, __FILE__);
385  e.setFatal();
386  throw e;
387  }
388  if (det == 0) {
389  GFException e("cannot invert matrix in Daf - det=0", __LINE__, __FILE__);
390  e.setFatal();
391  throw e;
392  }
393 }
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  hits.at(0);
112  pl = hits.at(0)->getDetPlane(rep);
113  //do the extrapolation
114  rep->extrapolate(pl, state, cov);
115 
116  if (cov[0][0] < 1.E-50) {
117  GFException exc(COVEXC, __LINE__, __FILE__);
118  throw exc;
119  }
120  }
121  catch (GFException& exc) {
122  std::cerr << exc.what();
123  exc.info();
124  for (unsigned int i = 0; i < planes.at(ipl)->size(); ++i)
125  trk->addFailedHit(irep, planes.at(ipl)->at(i));
126  if (exc.isFatal()) {
127  rep->setStatusFlag(1);
128  //abort();
129  }
130  continue; //go to next rep immediately
131  }
132  }
133  else {
134  pl = rep->getReferencePlane();
135  state = rep->getState();
136  cov = rep->getCov();
137  }
138  //the H matrix has to be identical for all hits in the plane
139  TMatrixT<Double_t> H = hits.at(0)->getHMatrix(rep);
140  bk->setMatrix("fPredSt", planes.at(ipl)->at(0), state);
141  bk->setMatrix("fPredCov", planes.at(ipl)->at(0), cov);
142 
143  TMatrixT<Double_t> covInv;
144  invertMatrix(cov, covInv);
145 
146  TMatrixT<Double_t> Htransp(H);
147  Htransp.T();
148 
149  bk->setMatrix("H", planes.at(ipl)->at(0), H);
150  bk->setDetPlane("pl", planes.at(ipl)->at(0), pl);
151 
152  TMatrixT<Double_t> stMod(state.GetNrows(), 1);
153 
154  double sumPk(0);
155  for (unsigned int ihit = 0; ihit < nhitsInPlane; ++ihit) {
156  double pki;
157  bk->getNumber("p", planes.at(ipl)->at(ihit), pki);
158  sumPk += pki;
159  }
160  TMatrixT<Double_t> V = hits.at(0)->getHitCov(pl);
161  TMatrixT<Double_t> Vinv;
162  invertMatrix(V, Vinv);
163  bk->setMatrix("V", planes.at(ipl)->at(0), V);
164  for (unsigned int ihit = 0; ihit < nhitsInPlane; ++ihit) {
165  //std::cout << "%%%% forward hit " << planes.at(ipl)->at(ihit) << std::endl;
166 
167  TMatrixT<Double_t> m = hits.at(ihit)->getHitCoord(pl);
168  bk->setMatrix("m", planes.at(ipl)->at(ihit), m);
169 
170  double pki;
171  bk->getNumber("p", planes.at(ipl)->at(ihit), pki);
172  TMatrixT<Double_t> Gain = calcGain(cov, V, H, sumPk);
173  //std::cout << "using weight " << pki << std::endl;
174  stMod += pki * Gain * (m - H * state);
175  }
176  state += stMod;
177  invertMatrix(covInv + sumPk * Htransp * Vinv * H, cov);
178  rep->setData(state, pl, &cov);
179 
180  } //loop over reps
181  } //loop over planes for forward fit
182  blowUpCovs(trk);
183  //now do the loop over the planes in backward direction
184  for (int ipl = nPlanes - 1; ipl >= 0; --ipl) {
185  //std::cout << "@bw@ plane " << ipl << std::endl;
186  std::vector<GFAbsRecoHit*> hits;
187  unsigned int nhitsInPlane = planes.at(ipl)->size();
188  for (unsigned int ihitInPl = 0; ihitInPl < nhitsInPlane; ++ihitInPl) {
189  hits.push_back(trk->getHit(planes.at(ipl)->at(ihitInPl)));
190  }
191  //now hits contains pointers to all hits in the next plane
192  //for non-planar detectors this will always be just one hit
193  for (unsigned int irep = 0; irep < nreps; ++irep) {
194  GFAbsTrackRep* rep = trk->getTrackRep(irep);
195  if (rep->getStatusFlag() != 0) continue;
196  //std::cout << "@bw@ rep " << irep << std::endl;
197  GFBookkeeping* bk = trk->getBK(irep);
198  if (bk->hitFailed(planes.at(ipl)->at(0)) > 0) continue;
199  // get prototypes for matrices
200  int repDim = rep->getDim();
201  TMatrixT<Double_t> state(repDim, 1);
202  TMatrixT<Double_t> cov(repDim, repDim);
203  ;
204 
205  TMatrixT<Double_t> H;
206  bk->getMatrix("H", planes.at(ipl)->at(0), H);
207  TMatrixT<Double_t> Htransp(H);
208  Htransp.T();
209  GFDetPlane pl;
210  bk->getDetPlane("pl", planes.at(ipl)->at(0), pl);
211  //std::cout << "backward ### iBeta, ipl, irep " << iBeta<< " " << ipl << " " << irep << std::endl;
212  if (ipl < (nPlanes - 1)) {
213  try {
214  //do the extrapolation
215  rep->extrapolate(pl, state, cov);
216  if (cov[0][0] < 1.E-50) {
217  GFException exc(COVEXC, __LINE__, __FILE__);
218  throw exc;
219  }
220  }
221  catch (GFException& exc) {
222  //std::cout << __FILE__ << " " << __LINE__ << std::endl;
223  std::cerr << exc.what();
224  exc.info();
225  //TAG
226  for (unsigned int i = 0; i < planes.at(ipl)->size(); ++i)
227  trk->addFailedHit(irep, planes.at(ipl)->at(i));
228  if (exc.isFatal()) { rep->setStatusFlag(1); }
229  continue; //go to next rep immediately
230  }
231  }
232  else {
233  state = rep->getState();
234  cov = rep->getCov();
235  }
236  TMatrixT<Double_t> covInv;
237  invertMatrix(cov, covInv);
238 
239  bk->setMatrix("bPredSt", planes.at(ipl)->at(0), state);
240  bk->setMatrix("bPredCov", planes.at(ipl)->at(0), cov);
241 
242  TMatrixT<Double_t> stMod(state.GetNrows(), 1);
243 
244  double sumPk(0);
245  for (unsigned int ihit = 0; ihit < nhitsInPlane; ++ihit) {
246  double pki;
247  bk->getNumber("p", planes.at(ipl)->at(ihit), pki);
248  sumPk += pki;
249  }
250 
251  TMatrixT<Double_t> V;
252  bk->getMatrix("V", planes.at(ipl)->at(0), V);
253  TMatrixT<Double_t> Vinv;
254  invertMatrix(V, Vinv);
255 
256  for (unsigned int ihit = 0; ihit < nhitsInPlane; ++ihit) {
257  //std::cout << "%%bw%% hit " << planes.at(ipl)->at(ihit) << std::endl;
258  TMatrixT<Double_t> m;
259  bk->getMatrix("m", planes.at(ipl)->at(ihit), m);
260 
261  double pki;
262  bk->getNumber("p", planes.at(ipl)->at(ihit), pki);
263  TMatrixT<Double_t> Gain = calcGain(cov, V, H, sumPk);
264  //std::cout << "using pki " << pki << std::endl;
265 
266  stMod += pki * Gain * (m - H * state);
267  }
268  state += stMod;
269  invertMatrix(covInv + sumPk * Htransp * Vinv * H, cov);
270 
271  rep->setData(state, pl, &cov);
272  } //done with loop over reps
273  } //loop over planes for backward
274  //calculate smoothed states and covs
275  TMatrixT<Double_t> fSt, fCov, bSt, bCov, smooSt, smooCov, fCovInv, bCovInv, m, H, smooResid;
276  for (int ipl = 0; ipl < nPlanes; ++ipl) {
277  for (unsigned int irep = 0; irep < nreps; ++irep) {
278  if (trk->getTrackRep(irep)->getStatusFlag() != 0) continue;
279  GFBookkeeping* bk = trk->getBK(irep);
280  if (bk->hitFailed(planes.at(ipl)->at(0)) > 0) continue;
281  bk->getMatrix("fPredSt", planes.at(ipl)->at(0), fSt);
282  bk->getMatrix("bPredSt", planes.at(ipl)->at(0), bSt);
283  bk->getMatrix("fPredCov", planes.at(ipl)->at(0), fCov);
284  bk->getMatrix("bPredCov", planes.at(ipl)->at(0), bCov);
285  fCovInv.ResizeTo(fCov);
286  bCovInv.ResizeTo(bCov);
287  invertMatrix(fCov, fCovInv);
288  invertMatrix(bCov, bCovInv);
289  smooCov.ResizeTo(fCov);
290  invertMatrix(fCovInv + bCovInv, smooCov);
291  smooSt.ResizeTo(fSt.GetNrows(), 1);
292  smooSt = smooCov * (fCovInv * fSt + bCovInv * bSt);
293  //std::cout << "#@#@#@#@#@# smooResid ipl, irep " << ipl << " " << irep << std::endl;
294  bk->getMatrix("H", planes.at(ipl)->at(0), H);
295  for (unsigned int ihit = 0; ihit < planes.at(ipl)->size(); ++ihit) {
296  //TAG
297  bk->getMatrix("m", planes.at(ipl)->at(ihit), m);
298  smooResid.ResizeTo(m.GetNrows(), 1);
299  smooResid = m - H * smooSt;
300  //std::cout << "%%%% hit " << planes.at(ipl)->at(ihit) << std::endl;
301  //smooResid.Print();
302  bk->setMatrix("smooResid", planes.at(ipl)->at(ihit), smooResid);
303  }
304  }
305  } //end loop over planes for smoothed state calculation
306 
307  //calculate the new weights
308  if (iBeta != fBeta.size() - 1) { //dont do it for the last beta, ause fitting will stop
309  for (unsigned int irep = 0; irep < nreps; ++irep) {
310  GFBookkeeping* bk = trk->getBK(irep);
311  for (int ipl = 0; ipl < nPlanes; ++ipl) {
312  if (trk->getTrackRep(irep)->getStatusFlag() != 0) continue;
313  if (bk->hitFailed(planes.at(ipl)->at(0)) > 0) continue;
314  double sumPhi = 0.;
315  double sumPhiCut = 0.;
316  std::vector<double> phi;
317  TMatrixT<Double_t> V;
318  bk->getMatrix("V", planes.at(ipl)->at(0), V);
319  V *= fBeta.at(iBeta);
320  TMatrixT<Double_t> Vinv;
321  invertMatrix(V, Vinv);
322  for (unsigned int ihit = 0; ihit < planes.at(ipl)->size(); ++ihit) {
323  //TMatrixT<Double_t> smooResid;
324  bk->getMatrix("smooResid", planes.at(ipl)->at(ihit), smooResid);
325  TMatrixT<Double_t> smooResidTransp(smooResid);
326  smooResidTransp.T();
327  int dimV = V.GetNrows();
328  double detV = V.Determinant();
329  TMatrixT<Double_t> expArg = smooResidTransp * Vinv * smooResid;
330  if (expArg.GetNrows() != 1)
331  throw std::runtime_error("genf::GFDaf::processTrack(): wrong matrix dimensions!");
332  double thisPhi =
333  1. / (pow(2. * TMath::Pi(), dimV / 2.) * sqrt(detV)) * exp(-0.5 * expArg[0][0]);
334  phi.push_back(thisPhi);
335  sumPhi += thisPhi;
336 
337  double cutVal = chi2Cuts[dimV];
338  if (cutVal <= 1.E-6)
339  throw std::runtime_error("genf::GFDaf::processTrack(): chi2 cut too low");
340  sumPhiCut += 1. / (2 * TMath::Pi() * sqrt(detV)) * exp(-0.5 * cutVal / fBeta.at(iBeta));
341  }
342  for (unsigned int ihit = 0; ihit < planes.at(ipl)->size(); ++ihit) {
343  bk->setNumber("p", planes.at(ipl)->at(ihit), phi.at(ihit) / (sumPhi + sumPhiCut));
344  }
345  }
346  }
347  }
348 
349  } //loop over betas
350 
351  return;
352 }
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:377
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:395
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:354
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 449 of file GFDaf.cxx.

References fBeta.

Referenced by GFDaf(), and setBlowUpFactor().

460 {
461 
462  /* // asserting version
463  assert(b1>0);fBeta.push_back(b1);
464  assert(b2>0 && b2<=b1);fBeta.push_back(b2);
465  if(b3>=0.) {
466  assert(b3<=b2);fBeta.push_back(b3);
467  if(b4>=0.) {
468  assert(b4<=b3);fBeta.push_back(b4);
469  if(b5>=0.) {
470  assert(b5<=b4);fBeta.push_back(b5);
471  if(b6>=0.) {
472  assert(b6<=b5);fBeta.push_back(b6);
473  if(b7>=0.) {
474  assert(b7<=b6);fBeta.push_back(b7);
475  if(b8>=0.) {
476  assert(b8<=b7);fBeta.push_back(b8);
477  if(b9>=0.) {
478  assert(b9<=b8);fBeta.push_back(b9);
479  if(b10>=0.) {
480  assert(b10<=b9);fBeta.push_back(b10);
481  }
482  }
483  }
484  }
485  }
486  }
487  }
488  }
489  */
490  if (b1 <= 0) throw std::runtime_error("genf::GFDaf::setBetas(): b1 <= 0");
491  fBeta.push_back(b1);
492 
493  if (b2 <= 0 || b2 > b1) throw std::runtime_error("genf::GFDaf::setBetas(): b2 in wrong range");
494  fBeta.push_back(b2);
495 
496  if (b3 < 0.) return;
497  if (b3 > b2) throw std::runtime_error("genf::GFDaf::setBetas(): b3 > b2");
498  fBeta.push_back(b3);
499 
500  if (b4 < 0.) return;
501  if (b4 > b3) throw std::runtime_error("genf::GFDaf::setBetas(): b4 > b3");
502  fBeta.push_back(b4);
503 
504  if (b5 < 0.) return;
505  if (b5 > b4) throw std::runtime_error("genf::GFDaf::setBetas(): b5 > b4");
506  fBeta.push_back(b5);
507 
508  if (b6 < 0.) return;
509  if (b6 > b5) throw std::runtime_error("genf::GFDaf::setBetas(): b6 > b5");
510  fBeta.push_back(b6);
511 
512  if (b7 < 0.) return;
513  if (b7 > b6) throw std::runtime_error("genf::GFDaf::setBetas(): b7 > b6");
514  fBeta.push_back(b7);
515 
516  if (b8 < 0.) return;
517  if (b8 > b7) throw std::runtime_error("genf::GFDaf::setBetas(): b8 > b7");
518  fBeta.push_back(b8);
519 
520  if (b9 < 0.) return;
521  if (b9 > b8) throw std::runtime_error("genf::GFDaf::setBetas(): b9 > b8");
522  fBeta.push_back(b9);
523 
524  if (b10 < 0.) return;
525  if (b10 > b9) throw std::runtime_error("genf::GFDaf::setBetas(): b10 > b9");
526  fBeta.push_back(b10);
527 }
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 417 of file GFDaf.cxx.

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

Referenced by GFDaf(), and setBlowUpFactor().

418 {
419 
420  if (fabs(val - 0.01) < 1.E-10) {
421  chi2Cuts[1] = 6.63;
422  chi2Cuts[2] = 9.21;
423  chi2Cuts[3] = 11.34;
424  chi2Cuts[4] = 13.23;
425  chi2Cuts[5] = 15.09;
426  }
427  else if (fabs(val - 0.005) < 1.E-10) {
428  chi2Cuts[1] = 7.88;
429  chi2Cuts[2] = 10.60;
430  chi2Cuts[3] = 12.84;
431  chi2Cuts[4] = 14.86;
432  chi2Cuts[5] = 16.75;
433  }
434  else if (fabs(val - 0.001) < 1.E-10) {
435  chi2Cuts[1] = 10.83;
436  chi2Cuts[2] = 13.82;
437  chi2Cuts[3] = 16.27;
438  chi2Cuts[4] = 18.47;
439  chi2Cuts[5] = 20.51;
440  }
441  else {
442  throw GFException(
443  "genf::GFTrackCand::GFDafsetProbCut(): value not supported", __LINE__, __FILE__)
444  .setFatal()
445  .setNumbers("value", {val});
446  }
447 }
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: