37 #include <Math/QuantFuncMathCore.h>
43 DAF::DAF(
bool useRefKalman,
double deltaPval,
double deltaWeight)
67 if (dynamic_cast<KalmanFitterRefTrack*>(
kalman_.get()) != NULL) {
79 debugOut<<
"DAF::processTrack //////////////////////////////////////////////////////////////// \n";
83 bool oneLastIter =
false;
87 for(
unsigned int iBeta = 0;; ++iBeta) {
90 debugOut<<
"DAF::processTrack, trackRep " << rep <<
", iteration " << iBeta+1 <<
", beta = " <<
betas_.at(iBeta) <<
"\n";
93 kalman_->processTrackWithRep(tr, rep, resortHits);
104 debugOut <<
"DAF::Kalman could not fit!\n";
110 if( oneLastIter ==
true){
112 debugOut <<
"DAF::break after one last iteration\n";
123 debugOut <<
"DAF::number of max iterations reached!\n";
130 bool converged(
false);
136 debugOut <<
"converged by Pval = " << status->
getBackwardPVal() <<
" even though weights changed at iBeta = " << iBeta << std::endl;
155 debugOut <<
"DAF::convergence reached in iteration " << iBeta+1 <<
" -> Do one last iteration with updated weights.\n";
175 for (
int i = 1; i < 7; ++i){
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__);
186 if ( measDim < 1 || measDim > 6 ){
187 Exception exc(
"DAF::addProbCut measDim must be in the range [1,6]",__LINE__,__FILE__);
191 chi2Cuts_[measDim] = ROOT::Math::chisquared_quantile_c( prob_cut, measDim);
195 void DAF::setBetas(
double b1,
double b2,
double b3,
double b4,
double b5,
double b6,
double b7,
double b8,
double b9,
double b10){
197 assert(b1>0);
betas_.push_back(b1);
199 assert(b2<=b1);
betas_.push_back(b2);
201 assert(b3<=b2);
betas_.push_back(b3);
203 assert(b4<=b3);
betas_.push_back(b4);
205 assert(b5<=b4);
betas_.push_back(b5);
207 assert(b6<=b5);
betas_.push_back(b6);
209 assert(b7<=b6);
betas_.push_back(b7);
211 assert(b8<=b7);
betas_.push_back(b8);
213 assert(b9<=b8);
betas_.push_back(b9);
215 assert(b10<=b9);
betas_.push_back(b10);
234 assert(bStart > bFinal);
235 assert(bFinal > 1.E-10);
243 for (
unsigned int i=0; i<nSteps; ++i) {
244 betas_.push_back(bStart * pow(bFinal / bStart, i / (nSteps - 1.)));
261 bool converged(
true);
262 double maxAbsChange(0);
265 for (std::vector< TrackPoint* >::const_iterator tp = trackPoints.begin(); tp != trackPoints.end(); ++tp) {
266 if (! (*tp)->hasFitterInfo(rep)) {
270 if (dynamic_cast<KalmanFitterInfo*>(fi) == NULL){
271 Exception exc(
"DAF::getWeights ==> can only use KalmanFitterInfos",__LINE__,__FILE__);
278 debugOut<<
"weights are fixed, continue \n";
285 std::vector<double> phi(nMeas, 0.);
288 for(
unsigned int j=0; j<nMeas; j++) {
292 const TVectorD& resid(residual.
getState());
293 TMatrixDSym Vinv(residual.
getCov());
296 int hitDim = resid.GetNrows();
300 double twoPiN = 2.*M_PI;
304 twoPiN = pow(twoPiN, hitDim);
306 double chi2 = Vinv.Similarity(resid);
308 debugOut<<
"chi2 = " << chi2 <<
"\n";
312 double norm = 1./sqrt(twoPiN * detV);
314 phi[j] = norm*exp(-0.5*chi2/beta);
318 assert(cutVal>1.E-6);
320 phi_cut += norm*exp(-0.5*cutVal/beta);
328 for(
unsigned int j=0; j<nMeas; j++) {
329 double weight = phi[j]/(phi_sum+phi_cut);
336 if (absChange > maxAbsChange)
337 maxAbsChange = absChange;
343 debugOut<<
"\t new weight: " << weight;
353 debugOut <<
"max abs weight change = " << maxAbsChange <<
"\n";
361 void DAF::Streamer(TBuffer &R__b)
366 typedef ::genfit::DAF thisClass;
368 if (R__b.IsReading()) {
369 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
if (R__v) { }
372 baseClass0::Streamer(R__b);
376 std::vector<double> &R__stl =
betas_;
380 R__stl.reserve(R__n);
381 for (R__i = 0; R__i < R__n; R__i++) {
384 R__stl.push_back(R__t);
392 for (R__i = 0; R__i < R__n; R__i++) {
398 if (R__t >= 1 && R__t <= 6)
404 assert(n_chi2Cuts == 6);
406 R__b.ReadFastArray(&
chi2Cuts_[1], n_chi2Cuts);
411 R__b.CheckByteCount(R__s, R__c, thisClass::IsA());
413 R__c = R__b.WriteVersion(thisClass::IsA(), kTRUE);
416 baseClass0::Streamer(R__b);
419 std::vector<double> &R__stl =
betas_;
420 int R__n=int(R__stl.size());
423 std::vector<double>::iterator R__k;
424 for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
434 R__b.SetByteCount(R__c, kTRUE);
void processTrackWithRep(Track *tr, const AbsTrackRep *rep, bool resortHits=false)
Process a track using the 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.
const TVectorD & getState() const
std::vector< double > betas_
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.
Kalman filter implementation with linearization around a reference track.
double deltaPval_
Convergence criterion.
bool isFitted() const
Has the track been fitted?
bool calcWeights(Track *trk, const AbsTrackRep *rep, double beta)
Calculate and set the weights for the next fitting pass. Return if convergence is met...
void setIsFittedWithDaf(bool fittedWithDaf=true)
void setNumIterations(unsigned int numIterations)
Simple Kalman filter implementation.
unsigned int getNumMeasurements() const
const std::vector< genfit::TrackPoint * > & getPointsWithMeasurement() const
void setIsFitted(bool fitted=true)
virtual const char * what() const
Standard error message handling for exceptions. use like "std::cerr << e.what();".
void setFatal(bool b=true)
Set fatal flag.
Abstract base class for a track representation.
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...
void info()
Print information in the exception object.
Exception class for error handling in GENFIT (provides storage for diagnostic information) ...
void setWeight(double weight)
double getForwardPVal() const
Abstract base class for Kalman fitter and derived fitting algorithms.
FitStatus for use with AbsKalmanFitter implementations.
int getNFailedPoints() const
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.
Measured coordinates on a plane.
void setAnnealingScheme(double bStart, double bFinal, unsigned int nSteps)
Configure the annealing scheme.
void setIsFitConvergedPartially(bool fitConverged=true)
const TMatrixDSym & getCov() const
void setIsFitConvergedFully(bool fitConverged=true)
AbsKalmanFitter(unsigned int maxIterations=4, double deltaPval=1e-3, double blowUpFactor=1e3)
boost::scoped_ptr< AbsKalmanFitter > kalman_
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...
FitStatus * getFitStatus(const AbsTrackRep *rep=NULL) const
Get FitStatus for a AbsTrackRep. Per default, return FitStatus for cardinalRep.