37 #include <Math/QuantFuncMathCore.h>
43 DAF::DAF(
bool useRefKalman,
double deltaPval,
double deltaWeight)
67 if (dynamic_cast<KalmanFitterRefTrack*>(
kalman_.get()) !=
nullptr) {
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);
206 for (
unsigned int i=0;
i<nSteps; ++
i) {
207 betas_.push_back(bStart * pow(bFinal / bStart,
i / (nSteps - 1.)));
224 bool converged(
true);
225 double maxAbsChange(0);
228 for (std::vector< TrackPoint* >::const_iterator
tp = trackPoints.begin();
tp != trackPoints.end(); ++
tp) {
229 if (! (*tp)->hasFitterInfo(rep)) {
233 if (dynamic_cast<KalmanFitterInfo*>(fi) ==
nullptr){
234 Exception exc(
"DAF::getWeights ==> can only use KalmanFitterInfos",__LINE__,__FILE__);
241 debugOut<<
"weights are fixed, continue \n";
248 std::vector<double>
phi(nMeas, 0.);
251 for(
unsigned int j=0;
j<nMeas;
j++) {
255 const TVectorD& resid(residual.
getState());
256 TMatrixDSym Vinv(residual.
getCov());
259 int hitDim = resid.GetNrows();
263 double twoPiN = 2.*M_PI;
267 twoPiN = pow(twoPiN, hitDim);
269 double chi2 = Vinv.Similarity(resid);
271 debugOut<<
"chi2 = " << chi2 <<
"\n";
275 double norm = 1./sqrt(twoPiN * detV);
277 phi[
j] = norm*exp(-0.5*chi2/beta);
283 phi_cut += norm*exp(-0.5*cutVal/beta);
291 for(
unsigned int j=0;
j<nMeas;
j++) {
292 double weight = phi[
j]/(phi_sum+phi_cut);
299 if (absChange > maxAbsChange)
300 maxAbsChange = absChange;
306 debugOut<<
"\t new weight: " << weight;
316 debugOut <<
"max abs weight change = " << maxAbsChange <<
"\n";
324 void DAF::Streamer(TBuffer &R__b)
329 typedef ::genfit::DAF thisClass;
331 if (R__b.IsReading()) {
332 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
if (R__v) { }
335 baseClass0::Streamer(R__b);
339 std::vector<double> &R__stl =
betas_;
343 R__stl.reserve(R__n);
344 for (R__i = 0; R__i < R__n; R__i++) {
347 R__stl.push_back(R__t);
355 for (R__i = 0; R__i < R__n; R__i++) {
361 if (R__t >= 1 && R__t <= 6)
369 R__b.ReadFastArray(&
chi2Cuts_[1], n_chi2Cuts);
374 R__b.CheckByteCount(R__s, R__c, thisClass::IsA());
376 R__c = R__b.WriteVersion(thisClass::IsA(), kTRUE);
379 baseClass0::Streamer(R__b);
382 std::vector<double> &R__stl =
betas_;
383 int R__n=int(R__stl.size());
386 std::vector<double>::iterator R__k;
387 for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
397 R__b.SetByteCount(R__c, kTRUE);