GENFIT  Rev:NoNumberAvailable
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules
DAF.cc
Go to the documentation of this file.
1 /* Copyright 2008-2013, Technische Universitaet Muenchen, Ludwig-Maximilians-Universität München
2  Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch & Tobias Schlüter
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 #include "DAF.h"
21 #include "Exception.h"
22 #include "IO.h"
23 #include "KalmanFitterInfo.h"
24 #include "KalmanFitter.h"
25 #include "KalmanFitterRefTrack.h"
26 #include "KalmanFitStatus.h"
27 #include "Tools.h"
28 #include "Track.h"
29 #include "TrackPoint.h"
30 
31 #include <assert.h>
32 #include <cmath>
33 
34 //root stuff
35 #include <TBuffer.h>
36 #include <TMath.h>
37 #include <Math/QuantFuncMathCore.h>
38 
39 
40 
41 namespace genfit {
42 
43 DAF::DAF(bool useRefKalman, double deltaPval, double deltaWeight)
44  : AbsKalmanFitter(10, deltaPval), deltaWeight_(deltaWeight)
45 {
46  if (useRefKalman) {
47  kalman_.reset(new KalmanFitterRefTrack());
48  static_cast<KalmanFitterRefTrack*>(kalman_.get())->setRefitAll();
49  }
50  else
51  kalman_.reset(new KalmanFitter());
52 
53  kalman_->setMultipleMeasurementHandling(weightedAverage);
54  kalman_->setMaxIterations(1);
55 
56  setAnnealingScheme(100, 0.1, 5); // also sets maxIterations_
57  setProbCut(0.001);
58 }
59 
60 DAF::DAF(AbsKalmanFitter* kalman, double deltaPval, double deltaWeight)
61  : AbsKalmanFitter(10, deltaPval), deltaWeight_(deltaWeight)
62 {
63  kalman_.reset(kalman);
64  kalman_->setMultipleMeasurementHandling(weightedAverage); // DAF makes no sense otherwise
65  kalman_->setMaxIterations(1);
66 
67  if (dynamic_cast<KalmanFitterRefTrack*>(kalman_.get()) != NULL) {
68  static_cast<KalmanFitterRefTrack*>(kalman_.get())->setRefitAll();
69  }
70 
71  setAnnealingScheme(100, 0.1, 5); // also sets maxIterations_
72  setProbCut(0.01);
73 }
74 
75 
76 void DAF::processTrackWithRep(Track* tr, const AbsTrackRep* rep, bool resortHits) {
77 
78  if (debugLvl_ > 0) {
79  debugOut<<"DAF::processTrack //////////////////////////////////////////////////////////////// \n";
80  }
81 
82  KalmanFitStatus* status = 0;
83  bool oneLastIter = false;
84 
85  double lastPval = -1;
86 
87  for(unsigned int iBeta = 0;; ++iBeta) {
88 
89  if (debugLvl_ > 0) {
90  debugOut<<"DAF::processTrack, trackRep " << rep << ", iteration " << iBeta+1 << ", beta = " << betas_.at(iBeta) << "\n";
91  }
92 
93  kalman_->processTrackWithRep(tr, rep, resortHits);
94 
95  status = static_cast<KalmanFitStatus*>(tr->getFitStatus(rep));
96  status->setIsFittedWithDaf();
97  status->setNumIterations(iBeta+1);
98 
99 
100  // check break conditions
101 
102  if (! status->isFitted()){
103  if (debugLvl_ > 0) {
104  debugOut << "DAF::Kalman could not fit!\n";
105  }
106  status->setIsFitted(false);
107  break;
108  }
109 
110  if( oneLastIter == true){
111  if (debugLvl_ > 0) {
112  debugOut << "DAF::break after one last iteration\n";
113  }
114  status->setIsFitConvergedFully(status->getNFailedPoints() == 0);
115  status->setIsFitConvergedPartially();
116  break;
117  }
118 
119  if(iBeta >= maxIterations_-1){
120  status->setIsFitConvergedFully(false);
121  status->setIsFitConvergedPartially(false);
122  if (debugLvl_ > 0) {
123  debugOut << "DAF::number of max iterations reached!\n";
124  }
125  break;
126  }
127 
128 
129  // get and update weights
130  bool converged(false);
131  try{
132  converged = calcWeights(tr, rep, betas_.at(iBeta));
133  if (!converged && iBeta >= minIterations_-1 &&
134  status->getBackwardPVal() != 0 && abs(lastPval - status->getBackwardPVal()) < this->deltaPval_) {
135  if (debugLvl_ > 0) {
136  debugOut << "converged by Pval = " << status->getBackwardPVal() << " even though weights changed at iBeta = " << iBeta << std::endl;
137  }
138  converged = true;
139  }
140  lastPval = status->getBackwardPVal();
141  } catch(Exception& e) {
142  errorOut<<e.what();
143  e.info();
144  //errorOut << "calc weights failed" << std::endl;
145  //mini_trk->getTrackRep(0)->setStatusFlag(1);
146  status->setIsFitted(false);
147  status->setIsFitConvergedFully(false);
148  status->setIsFitConvergedPartially(false);
149  break;
150  }
151 
152  // check if converged
153  if (iBeta >= minIterations_-1 && converged) {
154  if (debugLvl_ > 0) {
155  debugOut << "DAF::convergence reached in iteration " << iBeta+1 << " -> Do one last iteration with updated weights.\n";
156  }
157  oneLastIter = true;
158  status->setIsFitConvergedFully(status->getNFailedPoints() == 0);
159  status->setIsFitConvergedPartially();
160  }
161 
162  } // end loop over betas
163 
164 
165  if (status->getForwardPVal() == 0. &&
166  status->getBackwardPVal() == 0.) {
167  status->setIsFitConvergedFully(false);
168  status->setIsFitConvergedPartially(false);
169  }
170 
171 }
172 
173 
174 void DAF::setProbCut(const double prob_cut){
175  for ( int i = 1; i < 7; ++i){
176  addProbCut(prob_cut, i);
177  }
178 }
179 
180 void DAF::addProbCut(const double prob_cut, const int measDim){
181  if ( prob_cut > 1.0 || prob_cut < 0.0){
182  Exception exc("DAF::addProbCut prob_cut is not between 0 and 1",__LINE__,__FILE__);
183  exc.setFatal();
184  throw exc;
185  }
186  if ( measDim < 1 || measDim > 6 ){
187  Exception exc("DAF::addProbCut measDim must be in the range [1,6]",__LINE__,__FILE__);
188  exc.setFatal();
189  throw exc;
190  }
191  chi2Cuts_[measDim] = ROOT::Math::chisquared_quantile_c( prob_cut, measDim);
192 }
193 
194 
195 void DAF::setBetas(double b1,double b2,double b3,double b4,double b5,double b6,double b7,double b8,double b9,double b10){
196  betas_.clear();
197  assert(b1>0);betas_.push_back(b1);
198  if(b2>0){
199  assert(b2<=b1);betas_.push_back(b2);
200  if(b3>=0.) {
201  assert(b3<=b2);betas_.push_back(b3);
202  if(b4>=0.) {
203  assert(b4<=b3);betas_.push_back(b4);
204  if(b5>=0.) {
205  assert(b5<=b4);betas_.push_back(b5);
206  if(b6>=0.) {
207  assert(b6<=b5);betas_.push_back(b6);
208  if(b7>=0.) {
209  assert(b7<=b6);betas_.push_back(b7);
210  if(b8>=0.) {
211  assert(b8<=b7);betas_.push_back(b8);
212  if(b9>=0.) {
213  assert(b9<=b8);betas_.push_back(b9);
214  if(b10>=0.) {
215  assert(b10<=b9);betas_.push_back(b10);
216  }
217  }
218  }
219  }
220  }
221  }
222  }
223  }
224  }
225  minIterations_ = betas_.size();
226  maxIterations_ = betas_.size() + 4;
227  betas_.resize(maxIterations_,betas_.back()); //make sure main loop has a maximum of maxIterations_ and also make sure the last beta value is used for if more iterations are needed then the ones set by the user.
228 }
229 
230 
231 void DAF::setAnnealingScheme(double bStart, double bFinal, unsigned int nSteps) {
232  // The betas are calculated as a geometric series that takes nSteps
233  // steps to go from bStart to bFinal.
234  assert(bStart > bFinal);
235  assert(bFinal > 1.E-10);
236  assert(nSteps > 1);
237 
238  minIterations_ = nSteps;
239  maxIterations_ = nSteps + 4;
240 
241  betas_.clear();
242 
243  for (unsigned int i=0; i<nSteps; ++i) {
244  betas_.push_back(bStart * pow(bFinal / bStart, i / (nSteps - 1.)));
245  }
246 
247  betas_.resize(maxIterations_,betas_.back()); //make sure main loop has a maximum of 10 iterations and also make sure the last beta value is used for if more iterations are needed then the ones set by the user.
248 
249  /*for (unsigned int i=0; i<betas_.size(); ++i) {
250  debugOut<< betas_.at(i) << ", ";
251  }*/
252 }
253 
254 
255 bool DAF::calcWeights(Track* tr, const AbsTrackRep* rep, double beta) {
256 
257  if (debugLvl_ > 0) {
258  debugOut<<"DAF::calcWeights \n";
259  }
260 
261  bool converged(true);
262  double maxAbsChange(0);
263 
264  const std::vector< TrackPoint* >& trackPoints = tr->getPointsWithMeasurement();
265  for (std::vector< TrackPoint* >::const_iterator tp = trackPoints.begin(); tp != trackPoints.end(); ++tp) {
266  if (! (*tp)->hasFitterInfo(rep)) {
267  continue;
268  }
269  AbsFitterInfo* fi = (*tp)->getFitterInfo(rep);
270  if (dynamic_cast<KalmanFitterInfo*>(fi) == NULL){
271  Exception exc("DAF::getWeights ==> can only use KalmanFitterInfos",__LINE__,__FILE__);
272  throw exc;
273  }
274  KalmanFitterInfo* kfi = static_cast<KalmanFitterInfo*>(fi);
275 
276  if (kfi->areWeightsFixed()) {
277  if (debugLvl_ > 0) {
278  debugOut<<"weights are fixed, continue \n";
279  }
280  continue;
281  }
282 
283  unsigned int nMeas = kfi->getNumMeasurements();
284 
285  std::vector<double> phi(nMeas, 0.);
286  double phi_sum = 0;
287  double phi_cut = 0;
288  for(unsigned int j=0; j<nMeas; j++) {
289 
290  try{
291  const MeasurementOnPlane& residual = kfi->getResidual(j, true, true);
292  const TVectorD& resid(residual.getState());
293  TMatrixDSym Vinv(residual.getCov());
294  double detV;
295  tools::invertMatrix(Vinv, &detV); // can throw an Exception
296  int hitDim = resid.GetNrows();
297  // Needed for normalization, special cases for the two common cases,
298  // shouldn't matter, but the original code made some efforts to make
299  // this calculation faster, and it's not complex ...
300  double twoPiN = 2.*M_PI;
301  if (hitDim == 2)
302  twoPiN *= twoPiN;
303  else if (hitDim > 2)
304  twoPiN = pow(twoPiN, hitDim);
305 
306  double chi2 = Vinv.Similarity(resid);
307  if (debugLvl_ > 1) {
308  debugOut<<"chi2 = " << chi2 << "\n";
309  }
310 
311  // The common factor beta is eliminated.
312  double norm = 1./sqrt(twoPiN * detV);
313 
314  phi[j] = norm*exp(-0.5*chi2/beta);
315  phi_sum += phi[j];
316  //errorOut << "hitDim " << hitDim << " fchi2Cuts[hitDim] " << fchi2Cuts[hitDim] << std::endl;
317  double cutVal = chi2Cuts_[hitDim];
318  assert(cutVal>1.E-6);
319  //the following assumes that in the competing hits could have different V otherwise calculation could be simplified
320  phi_cut += norm*exp(-0.5*cutVal/beta);
321  }
322  catch(Exception& e) {
323  errorOut << e.what();
324  e.info();
325  }
326  }
327 
328  for(unsigned int j=0; j<nMeas; j++) {
329  double weight = phi[j]/(phi_sum+phi_cut);
330  //debugOut << phi_sum << " " << phi_cut << " " << weight << std::endl;
331 
332  // check convergence
333  double absChange(fabs(weight - kfi->getMeasurementOnPlane(j)->getWeight()));
334  if (converged && absChange > deltaWeight_) {
335  converged = false;
336  if (absChange > maxAbsChange)
337  maxAbsChange = absChange;
338  }
339 
340  if (debugLvl_ > 0) {
341  if (debugLvl_ > 1 || absChange > deltaWeight_) {
342  debugOut<<"\t old weight: " << kfi->getMeasurementOnPlane(j)->getWeight();
343  debugOut<<"\t new weight: " << weight;
344  }
345  }
346 
347  kfi->getMeasurementOnPlane(j)->setWeight(weight);
348  }
349  }
350 
351  if (debugLvl_ > 0) {
352  debugOut << "\t ";
353  debugOut << "max abs weight change = " << maxAbsChange << "\n";
354  }
355 
356  return converged;
357 }
358 
359 
360 // Customized from generated Streamer.
361 void DAF::Streamer(TBuffer &R__b)
362 {
363  // Stream an object of class genfit::DAF.
364 
365  //This works around a msvc bug and should be harmless on other platforms
366  typedef ::genfit::DAF thisClass;
367  UInt_t R__s, R__c;
368  if (R__b.IsReading()) {
369  Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v) { }
370  //This works around a msvc bug and should be harmless on other platforms
371  typedef genfit::AbsKalmanFitter baseClass0;
372  baseClass0::Streamer(R__b);
373  R__b >> deltaWeight_;
374  // weights_ are only of intermediate use -> not saved
375  {
376  std::vector<double> &R__stl = betas_;
377  R__stl.clear();
378  int R__i, R__n;
379  R__b >> R__n;
380  R__stl.reserve(R__n);
381  for (R__i = 0; R__i < R__n; R__i++) {
382  double R__t;
383  R__b >> R__t;
384  R__stl.push_back(R__t);
385  }
386  }
387  if (R__v == 1) {
388  // Old versions kept the chi2Cuts_ in a map.
389  // We ignore non-sensical dimensionalities when reading it again.
390  int R__i, R__n;
391  R__b >> R__n;
392  for (R__i = 0; R__i < R__n; R__i++) {
393  memset(chi2Cuts_, 0, sizeof(chi2Cuts_));
394  int R__t;
395  R__b >> R__t;
396  double R__t2;
397  R__b >> R__t2;
398  if (R__t >= 1 && R__t <= 6)
399  chi2Cuts_[R__t] = R__t2;
400  }
401  } else {
402  char n_chi2Cuts; // should be six
403  R__b >> n_chi2Cuts;
404  assert(n_chi2Cuts == 6); // Cannot be different as long as sanity prevails.
405  chi2Cuts_[0] = 0; // nonsensical.
406  R__b.ReadFastArray(&chi2Cuts_[1], n_chi2Cuts);
407  }
408  AbsKalmanFitter *p;
409  R__b >> p;
410  kalman_.reset(p);
411  R__b.CheckByteCount(R__s, R__c, thisClass::IsA());
412  } else {
413  R__c = R__b.WriteVersion(thisClass::IsA(), kTRUE);
414  //This works around a msvc bug and should be harmless on other platforms
415  typedef genfit::AbsKalmanFitter baseClass0;
416  baseClass0::Streamer(R__b);
417  R__b << deltaWeight_;
418  {
419  std::vector<double> &R__stl = betas_;
420  int R__n=int(R__stl.size());
421  R__b << R__n;
422  if(R__n) {
423  std::vector<double>::iterator R__k;
424  for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
425  R__b << (*R__k);
426  }
427  }
428  }
429  {
430  R__b << (char)6; // Number of chi2Cuts_
431  R__b.WriteFastArray(&chi2Cuts_[1], 6);
432  }
433  R__b << kalman_.get();
434  R__b.SetByteCount(R__c, kTRUE);
435  }
436 }
437 
438 } /* End of namespace genfit */
void processTrackWithRep(Track *tr, const AbsTrackRep *rep, bool resortHits=false)
Process a track using the DAF.
Definition: DAF.cc:76
DAF(const DAF &)
bool areWeightsFixed() const
Are the weights fixed?
unsigned int maxIterations_
Maximum number of iterations to attempt. Forward and backward are counted as one iteration.
double getBackwardPVal() const
void setProbCut(const double prob_cut)
Set the probability cut for the weight calculation for the hits.
Definition: DAF.cc:174
unsigned int debugLvl_
Definition: AbsFitter.h:55
const TVectorD & getState() const
Definition: StateOnPlane.h:61
std::vector< double > betas_
Definition: DAF.h:121
double chi2Cuts_[7]
Definition: DAF.h:122
unsigned int minIterations_
Minimum number of iterations to attempt. Forward and backward are counted as one iteration.
void setBetas(double b1, double b2=-1, 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.
Definition: DAF.cc:195
Kalman filter implementation with linearization around a reference track.
double deltaPval_
Convergence criterion.
bool isFitted() const
Has the track been fitted?
Definition: FitStatus.h:94
bool calcWeights(Track *trk, const AbsTrackRep *rep, double beta)
Calculate and set the weights for the next fitting pass. Return if convergence is met...
Definition: DAF.cc:255
double deltaWeight_
Definition: DAF.h:120
std::ostream debugOut
void setIsFittedWithDaf(bool fittedWithDaf=true)
void setNumIterations(unsigned int numIterations)
Simple Kalman filter implementation.
Definition: KalmanFitter.h:42
unsigned int getNumMeasurements() const
const std::vector< genfit::TrackPoint * > & getPointsWithMeasurement() const
Definition: Track.h:113
void setIsFitted(bool fitted=true)
Definition: FitStatus.h:128
virtual const char * what() const
Standard error message handling for exceptions. use like "std::cerr << e.what();".
Definition: Exception.cc:52
void setFatal(bool b=true)
Set fatal flag.
Definition: Exception.h:61
Abstract base class for a track representation.
Definition: AbsTrackRep.h:66
Collects information needed and produced by a AbsKalmanFitter implementations and is specific to one ...
void addProbCut(const double prob_cut, const int measDim)
Set the probability cut for the weight calculation for the hits for a specific measurement dimensiona...
Definition: DAF.cc:180
void info()
Print information in the exception object.
Definition: Exception.cc:56
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
Definition: Exception.h:48
void invertMatrix(const TMatrixDSym &mat, TMatrixDSym &inv, double *determinant=NULL)
Invert a matrix, throwing an Exception when inversion fails. Optional calculation of determinant...
Definition: Tools.cc:40
void setWeight(double weight)
double getForwardPVal() const
Abstract base class for Kalman fitter and derived fitting algorithms.
std::ostream errorOut
FitStatus for use with AbsKalmanFitter implementations.
int getNFailedPoints() const
Definition: FitStatus.h:110
MeasurementOnPlane getResidual(unsigned int iMeasurement=0, bool biased=false, bool onlyMeasurementErrors=true) const
Get unbiased (default) or biased residual from ith measurement.
Collection of TrackPoint objects, AbsTrackRep objects and FitStatus objects.
Definition: Track.h:71
Measured coordinates on a plane.
void setAnnealingScheme(double bStart, double bFinal, unsigned int nSteps)
Configure the annealing scheme.
Definition: DAF.cc:231
void setIsFitConvergedPartially(bool fitConverged=true)
Definition: FitStatus.h:130
const TMatrixDSym & getCov() const
void setIsFitConvergedFully(bool fitConverged=true)
Definition: FitStatus.h:129
AbsKalmanFitter(unsigned int maxIterations=4, double deltaPval=1e-3, double blowUpFactor=1e3)
boost::scoped_ptr< AbsKalmanFitter > kalman_
Definition: DAF.h:127
MeasurementOnPlane * getMeasurementOnPlane(int i=0) const
Defines for I/O streams used for error and debug printing.
This class collects all information needed and produced by a specific AbsFitter and is specific to on...
Definition: AbsFitterInfo.h:42
FitStatus * getFitStatus(const AbsTrackRep *rep=NULL) const
Get FitStatus for a AbsTrackRep. Per default, return FitStatus for cardinalRep.
Definition: Track.h:152