GENFIT
Rev:NoNumberAvailable
|
Go to the documentation of this file.
37 #include <Math/QuantFuncMathCore.h>
43 DAF::DAF(
bool useRefKalman,
double deltaPval,
double deltaWeight)
48 static_cast<KalmanFitterRefTrack*>(
kalman_.get())->setRefitAll();
67 if (dynamic_cast<KalmanFitterRefTrack*>(
kalman_.get()) != NULL) {
68 static_cast<KalmanFitterRefTrack*>(
kalman_.get())->setRefitAll();
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);
95 status = static_cast<KalmanFitStatus*>(tr->
getFitStatus(rep));
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 setWeight(double weight)
bool areWeightsFixed() const
Are the weights fixed?
Collection of TrackPoint objects, AbsTrackRep objects and FitStatus objects.
Collects information needed and produced by a AbsKalmanFitter implementations and is specific to one ...
void setAnnealingScheme(double bStart, double bFinal, unsigned int nSteps)
Configure the annealing scheme.
unsigned int getNumMeasurements() const
const TMatrixDSym & getCov() const
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.
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...
Exception class for error handling in GENFIT (provides storage for diagnostic information)
bool isFitted() const
Has the track been fitted?
Defines for I/O streams used for error and debug printing.
void setProbCut(const double prob_cut)
Set the probability cut for the weight calculation for the hits.
FitStatus * getFitStatus(const AbsTrackRep *rep=NULL) const
Get FitStatus for a AbsTrackRep. Per default, return FitStatus for cardinalRep.
void setIsFitted(bool fitted=true)
Abstract base class for a track representation.
void setFatal(bool b=true)
Set fatal flag.
virtual const char * what() const
Standard error message handling for exceptions. use like "std::cerr << e.what();".
Kalman filter implementation with linearization around a reference track.
AbsKalmanFitter(unsigned int maxIterations=4, double deltaPval=1e-3, double blowUpFactor=1e3)
This class collects all information needed and produced by a specific AbsFitter and is specific to on...
Simple Kalman filter implementation.
void info()
Print information in the exception object.
int getNFailedPoints() const
Abstract base class for Kalman fitter and derived fitting algorithms.
MeasurementOnPlane getResidual(unsigned int iMeasurement=0, bool biased=false, bool onlyMeasurementErrors=true) const
Get unbiased (default) or biased residual from ith measurement.
MeasurementOnPlane * getMeasurementOnPlane(int i=0) const
std::vector< double > betas_
bool calcWeights(Track *trk, const AbsTrackRep *rep, double beta)
Calculate and set the weights for the next fitting pass. Return if convergence is met....
unsigned int maxIterations_
Maximum number of iterations to attempt. Forward and backward are counted as one iteration.
Measured coordinates on a plane.
double getForwardPVal() const
void setIsFittedWithDaf(bool fittedWithDaf=true)
void setIsFitConvergedPartially(bool fitConverged=true)
void setIsFitConvergedFully(bool fitConverged=true)
const TVectorD & getState() const
void processTrackWithRep(Track *tr, const AbsTrackRep *rep, bool resortHits=false)
Process a track using the DAF.
void setNumIterations(unsigned int numIterations)
const std::vector< genfit::TrackPoint * > & getPointsWithMeasurement() const
unsigned int minIterations_
Minimum number of iterations to attempt. Forward and backward are counted as one iteration.
FitStatus for use with AbsKalmanFitter implementations.
double deltaPval_
Convergence criterion.
double getBackwardPVal() const
boost::scoped_ptr< AbsKalmanFitter > kalman_