LArSoft  v06_85_00
Liquid Argon Software toolkit - http://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 52 of file GFDaf.h.

Constructor & Destructor Documentation

genf::GFDaf::GFDaf ( )

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

Definition at line 45 of file GFDaf.cxx.

References setBetas(), and setProbCut().

45  :fBlowUpFactor(500.){
46  setProbCut(0.01);
47  setBetas(81,8,4,1,1,1);
48 }
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:452
genf::GFDaf::fBlowUpFactor
double fBlowUpFactor
Definition: GFDaf.h:108
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:422
genf::GFDaf::~GFDaf ( )

Definition at line 51 of file GFDaf.cxx.

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 356 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().

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

References invertMatrix().

Referenced by processTrack(), and setBlowUpFactor().

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

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

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

381  {
382  inv.ResizeTo(mat);
383  inv = (mat);
384  double det=0;
385  inv.Invert(&det);
386  if(TMath::IsNaN(det)) {
387  GFException e("Daf Gain: det of matrix is nan",__LINE__,__FILE__);
388  e.setFatal();
389  throw e;
390  }
391  if(det==0){
392  GFException e("cannot invert matrix in Daf - det=0",
393  __LINE__,__FILE__);
394  e.setFatal();
395  throw e;
396  }
397 
398 }
mat
Float_t mat
Definition: plot.C:40
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:50
e
Float_t e
Definition: plot.C:34
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().

53  {
54  trk->clearBookkeeping();
55 
56 
57 
58  unsigned int nreps=trk->getNumReps();
59  unsigned int nhits=trk->getNumHits();
60 
61  for(unsigned int i=0; i<nreps; ++i){
62  GFBookkeeping *bk = trk->getBK(i);
63  bk->setNhits(nhits);
64  bk->bookMatrices("fPredSt");
65  bk->bookMatrices("fPredCov");
66  bk->bookMatrices("bPredSt");
67  bk->bookMatrices("bPredCov");
68  bk->bookMatrices("H");
69  bk->bookMatrices("m");
70  bk->bookMatrices("V");
71  bk->bookGFDetPlanes("pl");
72  bk->bookMatrices("smooResid");
73  bk->bookNumbers("p",1.);
74  }
75 
76  //plane grouping
77  std::vector< std::vector<int>* > planes;
78  trk->getHitsByPlane(planes);
79  int nPlanes = planes.size();
80 
81  trk->clearFailedHits();
82  for(unsigned int iBeta=0; iBeta<fBeta.size(); iBeta++){
83  //std::cout << "@@ beta " << fBeta.at(iBeta) << std::endl;
84  if(iBeta>0) blowUpCovs(trk);
85  for(int ipl=0; ipl<nPlanes; ++ipl){
86  //std::cout << "@@ plane " << ipl << std::endl;
87  std::vector< GFAbsRecoHit* > hits;
88  unsigned int nhitsInPlane = planes.at(ipl)->size();
89  for(unsigned int ihitInPl=0;ihitInPl<nhitsInPlane;++ihitInPl){
90  hits.push_back( trk->getHit(planes.at(ipl)->at(ihitInPl)) );
91  }
92  //now hits contains pointers to all hits in the next plane
93  //for non-planar detectors this will always be just one hit
94  for(unsigned int irep=0; irep<nreps; ++irep){
95  GFAbsTrackRep* rep=trk->getTrackRep(irep);
96  if(rep->getStatusFlag()!=0) continue;
97  //std::cout << "@@ rep " << irep << std::endl;
98  GFBookkeeping *bk = trk->getBK(irep);
99  if(bk->hitFailed(planes.at(ipl)->at(0))>0) continue;
100  GFDetPlane pl;
101  // get prototypes for matrices
102  int repDim=rep->getDim();
103  TMatrixT<Double_t> state(repDim,1);
104  TMatrixT<Double_t> cov(repDim,repDim);;
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) 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  TMatrixT<Double_t> H;
204  bk->getMatrix("H",planes.at(ipl)->at(0),H);
205  TMatrixT<Double_t> Htransp(H);
206  Htransp.T();
207  GFDetPlane pl;
208  bk->getDetPlane("pl",planes.at(ipl)->at(0),pl);
209  //std::cout << "backward ### iBeta, ipl, irep " << iBeta<< " " << ipl << " " << irep << std::endl;
210  if(ipl<(nPlanes-1)){
211  try{
212  //do the extrapolation
213  rep->extrapolate(pl,state,cov);
214  if(cov[0][0]<1.E-50){
215  GFException exc(COVEXC,__LINE__,__FILE__);
216  throw exc;
217  }
218  }
219  catch(GFException& exc){
220  //std::cout << __FILE__ << " " << __LINE__ << std::endl;
221  std::cerr << exc.what();
222  exc.info();
223  //TAG
224  for(unsigned int i=0;i<planes.at(ipl)->size();++i) trk->addFailedHit(irep,planes.at(ipl)->at(i));
225  if(exc.isFatal()){
226  rep->setStatusFlag(1);
227  }
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 
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 
257  for(unsigned int ihit=0;ihit<nhitsInPlane;++ihit){
258  //std::cout << "%%bw%% hit " << planes.at(ipl)->at(ihit) << std::endl;
259  TMatrixT<Double_t> m;
260  bk->getMatrix("m",planes.at(ipl)->at(ihit),m);
261 
262  double pki;
263  bk->getNumber("p",planes.at(ipl)->at(ihit),pki);
264  TMatrixT<Double_t> Gain = calcGain(cov,V,H,sumPk);
265  //std::cout << "using pki " << pki << std::endl;
266 
267  stMod += pki*Gain*(m-H*state);
268  }
269  state += stMod;
270  invertMatrix( covInv + sumPk*Htransp*Vinv*H , cov);
271 
272  rep->setData(state,pl,&cov);
273  }//done with loop over reps
274  }//loop over planes for backward
275  //calculate smoothed states and covs
276  TMatrixT<Double_t> fSt,fCov,bSt,bCov,smooSt,smooCov,fCovInv,bCovInv,m,H,smooResid;
277  for(int ipl=0; ipl<nPlanes; ++ipl){
278  for(unsigned int irep=0; irep<nreps; ++irep){
279  if(trk->getTrackRep(irep)->getStatusFlag()!=0) continue;
280  GFBookkeeping *bk = trk->getBK(irep);
281  if(bk->hitFailed(planes.at(ipl)->at(0))>0) continue;
282  bk->getMatrix("fPredSt",planes.at(ipl)->at(0),fSt);
283  bk->getMatrix("bPredSt",planes.at(ipl)->at(0),bSt);
284  bk->getMatrix("fPredCov",planes.at(ipl)->at(0),fCov);
285  bk->getMatrix("bPredCov",planes.at(ipl)->at(0),bCov);
286  fCovInv.ResizeTo(fCov);
287  bCovInv.ResizeTo(bCov);
288  invertMatrix(fCov,fCovInv);
289  invertMatrix(bCov,bCovInv);
290  smooCov.ResizeTo(fCov);
291  invertMatrix( fCovInv + bCovInv , smooCov);
292  smooSt.ResizeTo(fSt.GetNrows(),1);
293  smooSt = smooCov * (fCovInv*fSt + bCovInv*bSt);
294  //std::cout << "#@#@#@#@#@# smooResid ipl, irep " << ipl << " " << irep << std::endl;
295  bk->getMatrix("H",planes.at(ipl)->at(0),H);
296  for(unsigned int ihit=0;ihit<planes.at(ipl)->size();++ihit){
297  //TAG
298  bk->getMatrix("m",planes.at(ipl)->at(ihit),m);
299  smooResid.ResizeTo(m.GetNrows(),1);
300  smooResid = m - H*smooSt;
301  //std::cout << "%%%% hit " << planes.at(ipl)->at(ihit) << std::endl;
302  //smooResid.Print();
303  bk->setMatrix("smooResid",planes.at(ipl)->at(ihit),smooResid);
304  }
305  }
306  }//end loop over planes for smoothed state calculation
307 
308  //calculate the new weights
309  if(iBeta!=fBeta.size()-1){//dont do it for the last beta, ause fitting will stop
310  for(unsigned int irep=0; irep<nreps; ++irep){
311  GFBookkeeping *bk = trk->getBK(irep);
312  for(int ipl=0; ipl<nPlanes; ++ipl){
313  if(trk->getTrackRep(irep)->getStatusFlag()!=0) continue;
314  if(bk->hitFailed(planes.at(ipl)->at(0))>0) continue;
315  double sumPhi=0.;
316  double sumPhiCut=0.;
317  std::vector<double> phi;
318  TMatrixT<Double_t> V;
319  bk->getMatrix("V",planes.at(ipl)->at(0),V);
320  V *= fBeta.at(iBeta);
321  TMatrixT<Double_t> Vinv;
322  invertMatrix(V,Vinv);
323  for(unsigned int ihit=0;ihit<planes.at(ipl)->size();++ihit){
324  //TMatrixT<Double_t> smooResid;
325  bk->getMatrix("smooResid",planes.at(ipl)->at(ihit),smooResid);
326  TMatrixT<Double_t> smooResidTransp(smooResid);
327  smooResidTransp.T();
328  int dimV = V.GetNrows();
329  double detV = V.Determinant();
330  TMatrixT<Double_t> expArg = smooResidTransp*Vinv*smooResid;
331  if (expArg.GetNrows() != 1)
332  throw std::runtime_error("genf::GFDaf::processTrack(): wrong matrix dimensions!");
333  double thisPhi = 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 
352  return;
353 
354  }
E
Float_t E
Definition: plot.C:23
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:381
GFException::what
virtual const char * what() const
standard error message handling for exceptions. use like "std::cerr << e.what();" ...
Definition: GFException.cxx:46
genf::GFDaf::chi2Cuts
std::map< int, double > chi2Cuts
Definition: GFDaf.h:110
hits
void hits()
Definition: readHits.C:15
genf::GFDaf::fBeta
std::vector< double > fBeta
Definition: GFDaf.h:109
GFException::isFatal
bool isFatal()
get fatal flag.
Definition: GFException.h:82
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:400
GFException::info
void info()
print information in the exception object
Definition: GFException.cxx:56
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:50
COVEXC
#define COVEXC
Definition: GFDaf.cxx:42
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:356
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 452 of file GFDaf.cxx.

References fBeta.

Referenced by GFDaf(), and setBlowUpFactor().

462  {
463 
464  /* // asserting version
465  assert(b1>0);fBeta.push_back(b1);
466  assert(b2>0 && b2<=b1);fBeta.push_back(b2);
467  if(b3>=0.) {
468  assert(b3<=b2);fBeta.push_back(b3);
469  if(b4>=0.) {
470  assert(b4<=b3);fBeta.push_back(b4);
471  if(b5>=0.) {
472  assert(b5<=b4);fBeta.push_back(b5);
473  if(b6>=0.) {
474  assert(b6<=b5);fBeta.push_back(b6);
475  if(b7>=0.) {
476  assert(b7<=b6);fBeta.push_back(b7);
477  if(b8>=0.) {
478  assert(b8<=b7);fBeta.push_back(b8);
479  if(b9>=0.) {
480  assert(b9<=b8);fBeta.push_back(b9);
481  if(b10>=0.) {
482  assert(b10<=b9);fBeta.push_back(b10);
483  }
484  }
485  }
486  }
487  }
488  }
489  }
490  }
491  */
492  if (b1 <= 0)
493  throw std::runtime_error("genf::GFDaf::setBetas(): b1 <= 0");
494  fBeta.push_back(b1);
495 
496  if (b2 <= 0 || b2 > b1)
497  throw std::runtime_error("genf::GFDaf::setBetas(): b2 in wrong range");
498  fBeta.push_back(b2);
499 
500  if(b3 < 0.) return;
501  if (b3 > b2)
502  throw std::runtime_error("genf::GFDaf::setBetas(): b3 > b2");
503  fBeta.push_back(b3);
504 
505  if(b4 < 0.) return;
506  if (b4 > b3)
507  throw std::runtime_error("genf::GFDaf::setBetas(): b4 > b3");
508  fBeta.push_back(b4);
509 
510  if(b5 < 0.) return;
511  if (b5 > b4)
512  throw std::runtime_error("genf::GFDaf::setBetas(): b5 > b4");
513  fBeta.push_back(b5);
514 
515  if(b6 < 0.) return;
516  if (b6 > b5)
517  throw std::runtime_error("genf::GFDaf::setBetas(): b6 > b5");
518  fBeta.push_back(b6);
519 
520  if(b7 < 0.) return;
521  if (b7 > b6)
522  throw std::runtime_error("genf::GFDaf::setBetas(): b7 > b6");
523  fBeta.push_back(b7);
524 
525  if(b8 < 0.) return;
526  if (b8 > b7)
527  throw std::runtime_error("genf::GFDaf::setBetas(): b8 > b7");
528  fBeta.push_back(b8);
529 
530  if(b9 < 0.) return;
531  if (b9 > b8)
532  throw std::runtime_error("genf::GFDaf::setBetas(): b9 > b8");
533  fBeta.push_back(b9);
534 
535  if(b10 < 0.) return;
536  if (b10 > b9)
537  throw std::runtime_error("genf::GFDaf::setBetas(): b10 > b9");
538  fBeta.push_back(b10);
539 
540 }
genf::GFDaf::fBeta
std::vector< double > fBeta
Definition: GFDaf.h:109
void genf::GFDaf::setBlowUpFactor ( double  f)
inline

Set the blowup factor (see blowUpCovs() )

Definition at line 73 of file GFDaf.h.

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

73 {fBlowUpFactor=f;}
genf::GFDaf::fBlowUpFactor
double fBlowUpFactor
Definition: GFDaf.h:108
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 422 of file GFDaf.cxx.

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

Referenced by GFDaf(), and setBlowUpFactor().

422  {
423 
424  if(fabs(val-0.01)<1.E-10){
425  chi2Cuts[1] = 6.63;
426  chi2Cuts[2] = 9.21;
427  chi2Cuts[3] = 11.34;
428  chi2Cuts[4] = 13.23;
429  chi2Cuts[5] = 15.09;
430  }
431  else if(fabs(val-0.005)<1.E-10){
432  chi2Cuts[1] = 7.88;
433  chi2Cuts[2] = 10.60;
434  chi2Cuts[3] = 12.84;
435  chi2Cuts[4] = 14.86;
436  chi2Cuts[5] = 16.75;
437  }
438  else if(fabs(val-0.001)<1.E-10){
439  chi2Cuts[1] = 10.83;
440  chi2Cuts[2] = 13.82;
441  chi2Cuts[3] = 16.27;
442  chi2Cuts[4] = 18.47;
443  chi2Cuts[5] = 20.51;
444  }
445  else{
446  throw GFException("genf::GFTrackCand::GFDafsetProbCut(): value not supported", __LINE__, __FILE__)
447  .setFatal().setNumbers("value", { val });
448  }
449 
450 }
GFException::setNumbers
GFException & setNumbers(std::string, const std::vector< double > &)
set list of numbers with description
Definition: GFException.cxx:32
E
Float_t E
Definition: plot.C:23
genf::GFDaf::chi2Cuts
std::map< int, double > chi2Cuts
Definition: GFDaf.h:110
GFException
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: GFException.h:50
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:80

Member Data Documentation

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

Definition at line 110 of file GFDaf.h.

Referenced by processTrack(), and setProbCut().

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

Definition at line 109 of file GFDaf.h.

Referenced by processTrack(), and setBetas().

double genf::GFDaf::fBlowUpFactor
private

Definition at line 108 of file GFDaf.h.

Referenced by blowUpCovs(), and setBlowUpFactor().

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