Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DAF.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DAF.cc
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()) != nullptr) {
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 
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 && fabs(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 void DAF::setAnnealingScheme(double bStart, double bFinal, unsigned int nSteps) {
195  // The betas are calculated as a geometric series that takes nSteps
196  // steps to go from bStart to bFinal.
197  assert(bStart > bFinal);
198  assert(bFinal > 1.E-10);
199  assert(nSteps > 1);
200 
201  minIterations_ = nSteps;
202  maxIterations_ = nSteps + 4;
203 
204  betas_.clear();
205 
206  for (unsigned int i=0; i<nSteps; ++i) {
207  betas_.push_back(bStart * pow(bFinal / bStart, i / (nSteps - 1.)));
208  }
209 
210  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.
211 
212  /*for (unsigned int i=0; i<betas_.size(); ++i) {
213  debugOut<< betas_.at(i) << ", ";
214  }*/
215 }
216 
217 
218 bool DAF::calcWeights(Track* tr, const AbsTrackRep* rep, double beta) {
219 
220  if (debugLvl_ > 0) {
221  debugOut<<"DAF::calcWeights \n";
222  }
223 
224  bool converged(true);
225  double maxAbsChange(0);
226 
227  const std::vector< TrackPoint* >& trackPoints = tr->getPointsWithMeasurement();
228  for (std::vector< TrackPoint* >::const_iterator tp = trackPoints.begin(); tp != trackPoints.end(); ++tp) {
229  if (! (*tp)->hasFitterInfo(rep)) {
230  continue;
231  }
232  AbsFitterInfo* fi = (*tp)->getFitterInfo(rep);
233  if (dynamic_cast<KalmanFitterInfo*>(fi) == nullptr){
234  Exception exc("DAF::getWeights ==> can only use KalmanFitterInfos",__LINE__,__FILE__);
235  throw exc;
236  }
237  KalmanFitterInfo* kfi = static_cast<KalmanFitterInfo*>(fi);
238 
239  if (kfi->areWeightsFixed()) {
240  if (debugLvl_ > 0) {
241  debugOut<<"weights are fixed, continue \n";
242  }
243  continue;
244  }
245 
246  unsigned int nMeas = kfi->getNumMeasurements();
247 
248  std::vector<double> phi(nMeas, 0.);
249  double phi_sum = 0;
250  double phi_cut = 0;
251  for(unsigned int j=0; j<nMeas; j++) {
252 
253  try{
254  const MeasurementOnPlane& residual = kfi->getResidual(j, true, true);
255  const TVectorD& resid(residual.getState());
256  TMatrixDSym Vinv(residual.getCov());
257  double detV;
258  tools::invertMatrix(Vinv, &detV); // can throw an Exception
259  int hitDim = resid.GetNrows();
260  // Needed for normalization, special cases for the two common cases,
261  // shouldn't matter, but the original code made some efforts to make
262  // this calculation faster, and it's not complex ...
263  double twoPiN = 2.*M_PI;
264  if (hitDim == 2)
265  twoPiN *= twoPiN;
266  else if (hitDim > 2)
267  twoPiN = pow(twoPiN, hitDim);
268 
269  double chi2 = Vinv.Similarity(resid);
270  if (debugLvl_ > 1) {
271  debugOut<<"chi2 = " << chi2 << "\n";
272  }
273 
274  // The common factor beta is eliminated.
275  double norm = 1./sqrt(twoPiN * detV);
276 
277  phi[j] = norm*exp(-0.5*chi2/beta);
278  phi_sum += phi[j];
279  //errorOut << "hitDim " << hitDim << " fchi2Cuts[hitDim] " << fchi2Cuts[hitDim] << std::endl;
280  double cutVal = chi2Cuts_[hitDim];
281  assert(cutVal>1.E-6);
282  //the following assumes that in the competing hits could have different V otherwise calculation could be simplified
283  phi_cut += norm*exp(-0.5*cutVal/beta);
284  }
285  catch(Exception& e) {
286  errorOut << e.what();
287  e.info();
288  }
289  }
290 
291  for(unsigned int j=0; j<nMeas; j++) {
292  double weight = phi[j]/(phi_sum+phi_cut);
293  //debugOut << phi_sum << " " << phi_cut << " " << weight << std::endl;
294 
295  // check convergence
296  double absChange(fabs(weight - kfi->getMeasurementOnPlane(j)->getWeight()));
297  if (converged && absChange > deltaWeight_) {
298  converged = false;
299  if (absChange > maxAbsChange)
300  maxAbsChange = absChange;
301  }
302 
303  if (debugLvl_ > 0) {
304  if (debugLvl_ > 1 || absChange > deltaWeight_) {
305  debugOut<<"\t old weight: " << kfi->getMeasurementOnPlane(j)->getWeight();
306  debugOut<<"\t new weight: " << weight;
307  }
308  }
309 
310  kfi->getMeasurementOnPlane(j)->setWeight(weight);
311  }
312  }
313 
314  if (debugLvl_ > 0) {
315  debugOut << "\t ";
316  debugOut << "max abs weight change = " << maxAbsChange << "\n";
317  }
318 
319  return converged;
320 }
321 
322 
323 // Customized from generated Streamer.
324 void DAF::Streamer(TBuffer &R__b)
325 {
326  // Stream an object of class genfit::DAF.
327 
328  //This works around a msvc bug and should be harmless on other platforms
329  typedef ::genfit::DAF thisClass;
330  UInt_t R__s, R__c;
331  if (R__b.IsReading()) {
332  Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v) { }
333  //This works around a msvc bug and should be harmless on other platforms
334  typedef genfit::AbsKalmanFitter baseClass0;
335  baseClass0::Streamer(R__b);
336  R__b >> deltaWeight_;
337  // weights_ are only of intermediate use -> not saved
338  {
339  std::vector<double> &R__stl = betas_;
340  R__stl.clear();
341  int R__i, R__n;
342  R__b >> R__n;
343  R__stl.reserve(R__n);
344  for (R__i = 0; R__i < R__n; R__i++) {
345  double R__t;
346  R__b >> R__t;
347  R__stl.push_back(R__t);
348  }
349  }
350  if (R__v == 1) {
351  // Old versions kept the chi2Cuts_ in a map.
352  // We ignore non-sensical dimensionalities when reading it again.
353  int R__i, R__n;
354  R__b >> R__n;
355  for (R__i = 0; R__i < R__n; R__i++) {
356  memset(chi2Cuts_, 0, sizeof(chi2Cuts_));
357  int R__t;
358  R__b >> R__t;
359  double R__t2;
360  R__b >> R__t2;
361  if (R__t >= 1 && R__t <= 6)
362  chi2Cuts_[R__t] = R__t2;
363  }
364  } else {
365  char n_chi2Cuts; // should be six
366  R__b >> n_chi2Cuts;
367  assert(n_chi2Cuts == 6); // Cannot be different as long as sanity prevails.
368  chi2Cuts_[0] = 0; // nonsensical.
369  R__b.ReadFastArray(&chi2Cuts_[1], n_chi2Cuts);
370  }
372  R__b >> p;
373  kalman_.reset(p);
374  R__b.CheckByteCount(R__s, R__c, thisClass::IsA());
375  } else {
376  R__c = R__b.WriteVersion(thisClass::IsA(), kTRUE);
377  //This works around a msvc bug and should be harmless on other platforms
378  typedef genfit::AbsKalmanFitter baseClass0;
379  baseClass0::Streamer(R__b);
380  R__b << deltaWeight_;
381  {
382  std::vector<double> &R__stl = betas_;
383  int R__n=int(R__stl.size());
384  R__b << R__n;
385  if(R__n) {
386  std::vector<double>::iterator R__k;
387  for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
388  R__b << (*R__k);
389  }
390  }
391  }
392  {
393  R__b << (char)6; // Number of chi2Cuts_
394  R__b.WriteFastArray(&chi2Cuts_[1], 6);
395  }
396  R__b << kalman_.get();
397  R__b.SetByteCount(R__c, kTRUE);
398  }
399 }
400 
401 } /* End of namespace genfit */