Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TLinearFitter.cxx
Go to the documentation of this file.
1// @(#)root/minuit:$Id$
2// Author: Anna Kreshuk 04/03/2005
3
4/*************************************************************************
5 * Copyright (C) 1995-2005, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12#include "TLinearFitter.h"
13#include "TMath.h"
14#include "TDecompChol.h"
15#include "TGraph.h"
16#include "TGraph2D.h"
17#include "TMultiGraph.h"
18#include "TRandom2.h"
19#include "TObjString.h"
20#include "TF2.h"
21#include "TH1.h"
22#include "TList.h"
23#include "TClass.h"
24#include "TBuffer.h"
25#include "TROOT.h"
26#include "strlcpy.h"
27
28#include <string>
29
30
31
32std::map<TString,TFormula*> TLinearFitter::fgFormulaMap;
33
34
35
36/**
37
38\class TLinearFitter
39
40\note An alternative to this class is to use ROOT::Fit::Fitter, calling the LinearFit() method.
41
42\ingroup MinuitOld
43
44The Linear Fitter - For fitting functions that are LINEAR IN PARAMETERS
45
46## The Linear Fitter
47
48Linear fitter is used to fit a set of data points with a linear
49combination of specified functions. Note, that "linear" in the name
50stands only for the model dependency on parameters, the specified
51functions can be nonlinear.
52The general form of this kind of model is
53~~~~
54 y(x) = a[0] + a[1]*f[1](x)+...a[n]*f[n](x)
55~~~~
56
57Functions f are fixed functions of x. For example, fitting with a
58polynomial is linear fitting in this sense.
59
60### Introduction
61
62#### The fitting method
63
64The fit is performed using the Normal Equations method with Cholesky
65decomposition.
66
67#### Why should it be used?
68
69The linear fitter is considerably faster than general non-linear
70fitters and doesn't require to set the initial values of parameters.
71
72### Using the fitter:
73
74### 1.Adding the data points:
75
76#### 1.1 To store or not to store the input data?
77 - There are 2 options in the constructor - to store or not
78 store the input data. The advantages of storing the data
79 are that you'll be able to reset the fitting model without
80 adding all the points again, and that for very large sets
81 of points the chisquare is calculated more precisely.
82 The obvious disadvantage is the amount of memory used to
83 keep all the points.
84 - Before you start adding the points, you can change the
85 store/not store option by StoreData() method.
86
87#### 1.2 The data can be added:
88 - simply point by point - AddPoint() method
89 - an array of points at once:
90 If the data is already stored in some arrays, this data
91 can be assigned to the linear fitter without physically
92 coping bytes, thanks to the Use() method of
93 TVector and TMatrix classes - AssignData() method
94
95### 2.Setting the formula
96
97#### 2.1 The linear formula syntax:
98 -Additive parts are separated by 2 plus signs "++"
99 --for example "1 ++ x" - for fitting a straight line
100 -All standard functions, undrestood by TFormula, can be used
101 as additive parts
102 --TMath functions can be used too
103 -Functions, used as additive parts, shouldn't have any parameters,
104 even if those parameters are set.
105 --for example, if normalizing a sum of a gaus(0, 1) and a
106 gaus(0, 2), don't use the built-in "gaus" of TFormula,
107 because it has parameters, take TMath::Gaus(x, 0, 1) instead.
108 -Polynomials can be used like "pol3", .."polN"
109 -If fitting a more than 3-dimensional formula, variables should
110 be numbered as follows:
111 -- x[0], x[1], x[2]... For example, to fit "1 ++ x[0] ++ x[1] ++ x[2] ++ x[3]*x[3]"
112
113#### 2.2 Setting the formula:
114
115##### 2.2.1 If fitting a 1-2-3-dimensional formula, one can create a
116 TF123 based on a linear expression and pass this function
117 to the fitter:
118 --Example:
119~~~~
120 TLinearFitter *lf = new TLinearFitter();
121 TF2 *f2 = new TF2("f2", "x ++ y ++ x*x*y*y", -2, 2, -2, 2);
122 lf->SetFormula(f2);
123~~~~
124 --The results of the fit are then stored in the function,
125 just like when the TH1::Fit or TGraph::Fit is used
126 --A linear function of this kind is by no means different
127 from any other function, it can be drawn, evaluated, etc.
128
129 --For multidimensional fitting, TFormulas of the form:
130 x[0]++...++x[n] can be used
131##### 2.2.2 There is no need to create the function if you don't want to,
132 the formula can be set by expression:
133 --Example:
134~~~~
135 // 2 is the number of dimensions
136 TLinearFitter *lf = new TLinearFitter(2);
137 lf->SetFormula("x ++ y ++ x*x*y*y");
138~~~~
139
140##### 2.2.3 The fastest functions to compute are polynomials and hyperplanes.
141 --Polynomials are set the usual way: "pol1", "pol2",...
142 --Hyperplanes are set by expression "hyp3", "hyp4", ...
143 ---The "hypN" expressions only work when the linear fitter
144 is used directly, not through TH1::Fit or TGraph::Fit.
145 To fit a graph or a histogram with a hyperplane, define
146 the function as "1++x++y".
147 ---A constant term is assumed for a hyperplane, when using
148 the "hypN" expression, so "hyp3" is in fact fitting with
149 "1++x++y++z" function.
150 --Fitting hyperplanes is much faster than fitting other
151 expressions so if performance is vital, calculate the
152 function values beforehand and give them to the fitter
153 as variables
154 --Example:
155 You want to fit "sin(x)|cos(2*x)" very fast. Calculate
156 sin(x) and cos(2*x) beforehand and store them in array *data.
157 Then:
158 TLinearFitter *lf=new TLinearFitter(2, "hyp2");
159 lf->AssignData(npoint, 2, data, y);
160
161#### 2.3 Resetting the formula
162
163##### 2.3.1 If the input data is stored (or added via AssignData() function),
164 the fitting formula can be reset without re-adding all the points.
165 --Example:
166~~~~
167 TLinearFitter *lf=new TLinearFitter("1++x++x*x");
168 lf->AssignData(n, 1, x, y, e);
169 lf->Eval()
170 //looking at the parameter significance, you see,
171 // that maybe the fit will improve, if you take out
172 // the constant term
173 lf->SetFormula("x++x*x");
174 lf->Eval();
175 ...
176~~~~
177
178##### 2.3.2 If the input data is not stored, the fitter will have to be
179 cleared and the data will have to be added again to try a
180 different formula.
181
182### 3.Accessing the fit results
183
184#### 3.1 There are methods in the fitter to access all relevant information:
185 --GetParameters, GetCovarianceMatrix, etc
186 --the t-values of parameters and their significance can be reached by
187 GetParTValue() and GetParSignificance() methods
188
189#### 3.2 If fitting with a pre-defined TF123, the fit results are also
190 written into this function.
191
192
193### 4.Robust fitting - Least Trimmed Squares regression (LTS)
194 Outliers are atypical(by definition), infrequant observations; data points
195 which do not appear to follow the characteristic distribution of the rest
196 of the data. These may reflect genuine properties of the underlying
197 phenomenon(variable), or be due to measurement errors or anomalies which
198 shouldn't be modelled. (StatSoft electronic textbook)
199
200 Even a single gross outlier can greatly influence the results of least-
201 squares fitting procedure, and in this case use of robust(resistant) methods
202 is recommended.
203
204 The method implemented here is based on the article and algorithm:
205 "Computing LTS Regression for Large Data Sets" by
206 P.J.Rousseeuw and Katrien Van Driessen
207 The idea of the method is to find the fitting coefficients for a subset
208 of h observations (out of n) with the smallest sum of squared residuals.
209 The size of the subset h should lie between (npoints + nparameters +1)/2
210 and n, and represents the minimal number of good points in the dataset.
211 The default value is set to (npoints + nparameters +1)/2, but of course
212 if you are sure that the data contains less outliers it's better to change
213 h according to your data.
214
215 To perform a robust fit, call EvalRobust() function instead of Eval() after
216 adding the points and setting the fitting function.
217 Note, that standard errors on parameters are not computed!
218
219*/
220
221
222
223////////////////////////////////////////////////////////////////////////////////
224///default c-tor, input data is stored
225///If you don't want to store the input data,
226///run the function StoreData(kFALSE) after constructor
227
229 fVal()
230{
231 fChisquare =0;
232 fNpoints =0;
233 fNdim =0;
234 fY2 =0;
235 fY2Temp =0;
236 fNfixed =0;
237 fIsSet =kFALSE;
238 fFormula =nullptr;
239 fFixedParams=nullptr;
240 fSpecial =0;
241 fInputFunction=nullptr;
244 fNfunctions = 0;
245 fFormulaSize = 0;
246 fH = 0;
247}
248
249////////////////////////////////////////////////////////////////////////////////
250///The parameter stands for number of dimensions in the fitting formula
251///The input data is stored. If you don't want to store the input data,
252///run the function StoreData(kFALSE) after constructor
253
255 fVal()
256{
257 fNdim =ndim;
258 fNpoints =0;
259 fY2 =0;
260 fY2Temp =0;
261 fNfixed =0;
262 fFixedParams=nullptr;
263 fFormula =nullptr;
264 fIsSet =kFALSE;
265 fChisquare=0;
266 fSpecial =0;
267 fInputFunction=nullptr;
269 fRobust = kFALSE;
270 fNfunctions = 0;
271 fFormulaSize = 0;
272 fH = 0;
273}
274
275////////////////////////////////////////////////////////////////////////////////
276///First parameter stands for number of dimensions in the fitting formula
277///Second parameter is the fitting formula: see class description for formula syntax
278///Options:
279///The option is to store or not to store the data
280///If you don't want to store the data, choose "" for the option, or run
281///StoreData(kFalse) member function after the constructor
282
283TLinearFitter::TLinearFitter(Int_t ndim, const char *formula, Option_t *opt)
284{
285 fNdim=ndim;
286 fNpoints=0;
287 fChisquare=0;
288 fY2=0;
289 fNfixed=0;
290 fFixedParams=nullptr;
291 fSpecial=0;
292 fInputFunction=nullptr;
293 fFormula = nullptr;
294 TString option=opt;
295 option.ToUpper();
296 if (option.Contains("D"))
298 else
301 SetFormula(formula);
302}
303
304////////////////////////////////////////////////////////////////////////////////
305///This constructor uses a linear function. How to create it?
306///TFormula now accepts formulas of the following kind:
307///TFormula("f", "x++y++z++x*x") or
308///TFormula("f", "x[0]++x[1]++x[2]*x[2]");
309///Other than the look, it's in no
310///way different from the regular formula, it can be evaluated,
311///drawn, etc.
312///The option is to store or not to store the data
313///If you don't want to store the data, choose "" for the option, or run
314///StoreData(kFalse) member function after the constructor
315
317{
318 fNdim=function->GetNdim();
319 if (!function->IsLinear()){
320 Int_t number=function->GetNumber();
321 if (number<299 || number>310){
322 Error("TLinearFitter", "Trying to fit with a nonlinear function");
323 return;
324 }
325 }
326 fNpoints=0;
327 fChisquare=0;
328 fY2=0;
329 fNfixed=0;
330 fFixedParams=nullptr;
331 fSpecial=0;
332 fFormula = nullptr;
333 TString option=opt;
334 option.ToUpper();
335 if (option.Contains("D"))
337 else
341 fInputFunction=nullptr;
342
343 SetFormula(function);
344}
345
346
347////////////////////////////////////////////////////////////////////////////////
348/// Linear fitter cleanup.
349
351{
352 if (fFormula) {
353 delete [] fFormula;
354 fFormula = nullptr;
355 }
356 if (fFixedParams) {
357 delete [] fFixedParams;
358 fFixedParams = nullptr;
359 }
360 fInputFunction = nullptr;
361
362 //fFunctions.Delete();
364
365}
366
367////////////////////////////////////////////////////////////////////////////////
368///Add another linear fitter to this linear fitter. Points and Design matrices
369///are added, but the previous fitting results (if any) are deleted.
370///Fitters must have same formulas (this is not checked). Fixed parameters are not changed
371
373{
374 fParams.Zero();
375 fParCovar.Zero();
376 fTValues.Zero();
377 fParSign.Zero();
378
379 fDesign += tlf->fDesign;
380
381 fDesignTemp += tlf->fDesignTemp;
382 fDesignTemp2 += tlf->fDesignTemp2;
383 fDesignTemp3 += tlf->fDesignTemp3;
384 fAtb += tlf->fAtb;
385 fAtbTemp += tlf->fAtbTemp;
386 fAtbTemp2 += tlf->fAtbTemp2;
387 fAtbTemp3 += tlf->fAtbTemp3;
388
389 if (fStoreData){
391 Int_t newsize = fNpoints+tlf->fNpoints;
392 if (size<newsize){
396 }
397 for (Int_t i=fNpoints; i<newsize; i++){
398 fY(i)=tlf->fY(i-fNpoints);
399 fE(i)=tlf->fE(i-fNpoints);
400 for (Int_t j=0; j<fNdim; j++){
401 fX(i,j)=tlf->fX(i-fNpoints, j);
402 }
403 }
404 }
405 fY2 += tlf->fY2;
406 fY2Temp += tlf->fY2Temp;
407 fNpoints += tlf->fNpoints;
408 //fInputFunction=(TFormula*)tlf.fInputFunction->Clone();
409
410 fChisquare=0;
411 fH=0;
412 fRobust=false;
413}
414
415////////////////////////////////////////////////////////////////////////////////
416///Adds 1 point to the fitter.
417///First parameter stands for the coordinates of the point, where the function is measured
418///Second parameter - the value being fitted
419///Third parameter - weight(measurement error) of this point (=1 by default)
420
422{
423 Int_t size;
424 fNpoints++;
425 if (fStoreData){
427 if (size<fNpoints){
431 }
432
433 Int_t j=fNpoints-1;
434 fY(j)=y;
435 fE(j)=e;
436 for (Int_t i=0; i<fNdim; i++)
437 fX(j,i)=x[i];
438 }
439 //add the point to the design matrix, if the formula has been set
440 // (LM: why condition !fRobust ?? - remove it to fix Coverity 11602)
441 // when 100< fSpecial < 200 (Polynomial) fFunctions is not empty
442 // (see implementation of SetFormula method)
443 if (fFunctions.IsEmpty()&&(!fInputFunction)&&(fSpecial<=200)){
444 Error("AddPoint", "Point can't be added, because the formula hasn't been set");
445 return;
446 }
447 if (!fRobust) AddToDesign(x, y, e);
448}
449
450////////////////////////////////////////////////////////////////////////////////
451///This function is to use when you already have all the data in arrays
452///and don't want to copy them into the fitter. In this function, the Use() method
453///of TVectorD and TMatrixD is used, so no bytes are physically moved around.
454///First parameter - number of points to fit
455///Second parameter - number of variables in the model
456///Third parameter - the variables of the model, stored in the following way:
457///(x0(0), x1(0), x2(0), x3(0), x0(1), x1(1), x2(1), x3(1),...
458
460{
461 if (npoints<fNpoints){
462 Error("AddData", "Those points are already added");
463 return;
464 }
466 if (fX.GetMatrixArray()==x && fY.GetMatrixArray()==y){
467 if (e){
468 if (fE.GetMatrixArray()==e)
469 same=kTRUE;
470 }
471 }
472
473 fX.Use(npoints, xncols, x);
474 fY.Use(npoints, y);
475 if (e)
476 fE.Use(npoints, e);
477 else {
479 fE=1;
480 }
481 Int_t xfirst;
482 if (!fFunctions.IsEmpty() || fInputFunction || fSpecial>200) {
483 if (same)
484 xfirst=fNpoints;
485
486 else
487 xfirst=0;
488 for (Int_t i=xfirst; i<npoints; i++)
489 AddToDesign(TMatrixDRow(fX, i).GetPtr(), fY(i), fE(i));
490 }
492}
493
494////////////////////////////////////////////////////////////////////////////////
495///Add a point to the AtA matrix and to the Atb vector.
496
498{
499
500
501 Int_t i, j, ii;
502 y/=e;
503
504// Double_t fVal[1000];
505
506 if ((fSpecial>100)&&(fSpecial<200)){
507 //polynomial fitting
508 Int_t npar=fSpecial-100;
509 fVal[0]=1;
510 for (i=1; i<npar; i++)
511 fVal[i]=fVal[i-1]*x[0];
512 for (i=0; i<npar; i++)
513 fVal[i]/=e;
514 } else if (fSpecial>200){
515 //Hyperplane fitting. Constant term is added
516 Int_t npar=fSpecial-201;
517 fVal[0]=1./e;
518 for (i=0; i<npar; i++)
519 fVal[i+1]=x[i]/e;
520 } else {
521 //general case
522 for (ii=0; ii<fNfunctions; ii++){
523 if (!fFunctions.IsEmpty()){
524 // ffunctions can be TF1 or TFormula depending on how they are created
526 if (obj->IsA() == TFormula::Class() ) {
527 TFormula *f1 = (TFormula*)(obj);
528 fVal[ii]=f1->EvalPar(x)/e;
529 }
530 else if (obj->IsA() == TF1::Class() ) {
531 TF1 *f1 = (TF1*)(obj);
532 fVal[ii]=f1->EvalPar(x)/e;
533 }
534 else {
535 Error("AddToDesign","Basis Function %s is of an invalid type %s",obj->GetName(),obj->IsA()->GetName());
536 return;
537 }
538 } else {
540 if (!f){
541 Error("AddToDesign","Function %s has no linear parts - maybe missing a ++ in the formula expression",fInputFunction->GetName());
542 return;
543 }
544 fVal[ii]=f->EvalPar(x)/e;
545 }
546 }
547 }
548 //additional matrices for numerical stability
549 for (i=0; i<fNfunctions; i++){
550 for (j=0; j<i; j++)
551 fDesignTemp3(j, i)+=fVal[i]*fVal[j];
552 fDesignTemp3(i, i)+=fVal[i]*fVal[i];
553 fAtbTemp3(i)+=fVal[i]*y;
554
555 }
556 fY2Temp+=y*y;
558
559 if (fNpoints % 100 == 0 && fNpoints>100){
563 fAtbTemp3.Zero();
564 if (fNpoints % 10000 == 0 && fNpoints>10000){
568 fAtbTemp2.Zero();
569 fY2+=fY2Temp;
570 fY2Temp=0;
571 if (fNpoints % 1000000 == 0 && fNpoints>1000000){
574 fAtb+=fAtbTemp;
575 fAtbTemp.Zero();
576 }
577 }
578 }
579}
580
581////////////////////////////////////////////////////////////////////////////////
582
603
604////////////////////////////////////////////////////////////////////////////////
605///Clears everything. Used in TH1::Fit and TGraph::Fit().
606
608{
609 fParams.Clear();
611 fTValues.Clear();
612 fParSign.Clear();
613 fDesign.Clear();
617 fAtb.Clear();
618 fAtbTemp.Clear();
622 fInputFunction=nullptr;
623 fY.Clear();
624 fX.Clear();
625 fE.Clear();
626
627 fNpoints=0;
628 fNfunctions=0;
629 fFormulaSize=0;
630 fNdim=0;
631 if (fFormula) delete [] fFormula;
632 fFormula=nullptr;
633 fIsSet=false;
634 if (fFixedParams) delete [] fFixedParams;
635 fFixedParams=nullptr;
636
637 fChisquare=0;
638 fY2=0;
639 fSpecial=0;
642}
643
644////////////////////////////////////////////////////////////////////////////////
645///To be used when different sets of points are fitted with the same formula.
646
648{
649 fDesign.Zero();
650 fAtb.Zero();
654 fAtbTemp.Zero();
655 fAtbTemp2.Zero();
656 fAtbTemp3.Zero();
657
658 fParams.Zero();
659 fParCovar.Zero();
660 fTValues.Zero();
661 fParSign.Zero();
662
663 for (Int_t i=0; i<fNfunctions; i++)
664 fFixedParams[i]=false;
665 fChisquare=0;
666 fNpoints=0;
667
668}
669
670////////////////////////////////////////////////////////////////////////////////
671///Calculates the chisquare.
672
674{
675 Int_t i, j;
677 Double_t temp, temp2;
678
679 if (!fStoreData){
680 sumtotal2 = 0;
681 for (i=0; i<fNfunctions; i++){
682 for (j=0; j<i; j++){
683 sumtotal2 += 2*fParams(i)*fParams(j)*fDesign(j, i);
684 }
685 sumtotal2 += fParams(i)*fParams(i)*fDesign(i, i);
686 sumtotal2 -= 2*fParams(i)*fAtb(i);
687 }
688 sumtotal2 += fY2;
689 } else {
690 sumtotal2 = 0;
691 if (fInputFunction){
692 for (i=0; i<fNpoints; i++){
693 temp = fInputFunction->EvalPar(TMatrixDRow(fX, i).GetPtr());
694 temp2 = (fY(i)-temp)*(fY(i)-temp);
695 temp2 /= fE(i)*fE(i);
696 sumtotal2 += temp2;
697 }
698 } else {
699 sumtotal2 = 0;
700 Double_t val[100];
701 for (Int_t point=0; point<fNpoints; point++){
702 temp = 0;
703 if ((fSpecial>100)&&(fSpecial<200)){
704 Int_t npar = fSpecial-100;
705 val[0] = 1;
706 for (i=1; i<npar; i++)
707 val[i] = val[i-1]*fX(point, 0);
708 for (i=0; i<npar; i++)
709 temp += fParams(i)*val[i];
710 } else {
711 if (fSpecial>200) {
712 //hyperplane case
713 Int_t npar = fSpecial-201;
714 temp+=fParams(0);
715 for (i=0; i<npar; i++)
716 temp += fParams(i+1)*fX(point, i);
717 } else {
718 for (j=0; j<fNfunctions; j++) {
720 val[j] = f1->EvalPar(TMatrixDRow(fX, point).GetPtr());
721 temp += fParams(j)*val[j];
722 }
723 }
724 }
725 temp2 = (fY(point)-temp)*(fY(point)-temp);
726 temp2 /= fE(point)*fE(point);
727 sumtotal2 += temp2;
728 }
729 }
730 }
732
733}
734
735////////////////////////////////////////////////////////////////////////////////
736/// Computes parameters' t-values and significance
737
739{
740 for (Int_t i=0; i<fNfunctions; i++){
741 fTValues(i) = fParams(i)/(TMath::Sqrt(fParCovar(i, i)));
743 }
744}
745
746////////////////////////////////////////////////////////////////////////////////
747/// Perform the fit and evaluate the parameters
748/// Returns 0 if the fit is ok, 1 if there are errors
749
751{
752 Double_t e;
753 if (fFunctions.IsEmpty()&&(!fInputFunction)&&(fSpecial<=200)){
754 Error("TLinearFitter::Eval", "The formula hasn't been set");
755 return 1;
756 }
757 //
762
763 fChisquare=0;
764
765 if (!fIsSet){
767 if (!update){
768 //no points to fit
769 fParams.Zero();
770 fParCovar.Zero();
771 fTValues.Zero();
772 fParSign.Zero();
773 fChisquare=0;
774 if (fInputFunction){
776 for (Int_t i=0; i<fNfunctions; i++)
777 ((TF1*)fInputFunction)->SetParError(i, 0);
778 ((TF1*)fInputFunction)->SetChisquare(0);
779 ((TF1*)fInputFunction)->SetNDF(0);
780 ((TF1*)fInputFunction)->SetNumberFitPoints(0);
781 }
782 return 1;
783 }
784 }
785
787
788//fixing fixed parameters, if there are any
789 Int_t i, ii, j=0;
790 if (fNfixed>0){
791 for (ii=0; ii<fNfunctions; ii++)
793 for (i=0; i<fNfunctions; i++){
794 if (fFixedParams[i]){
795 for (ii=0; ii<i; ii++)
796 fDesignTemp(ii, j) = fDesign(ii, i);
797 for (ii=i; ii<fNfunctions; ii++)
798 fDesignTemp(ii, j) = fDesign(i, ii);
799 j++;
800 for (ii=0; ii<fNfunctions; ii++){
801 fAtb(ii)-=fParams(i)*(fDesignTemp(ii, j-1));
802 }
803 }
804 }
805 for (i=0; i<fNfunctions; i++){
806 if (fFixedParams[i]){
807 for (ii=0; ii<fNfunctions; ii++){
808 fDesign(ii, i) = 0;
809 fDesign(i, ii) = 0;
810 }
811 fDesign (i, i) = 1;
812 fAtb(i) = fParams(i);
813 }
814 }
815 }
816
818 Bool_t ok;
819 TVectorD coef(fNfunctions);
820 coef=chol.Solve(fAtb, ok);
821 if (!ok){
822 Error("Eval","Matrix inversion failed");
823 fParams.Zero();
824 fParCovar.Zero();
825 fTValues.Zero();
826 fParSign.Zero();
827 if (fInputFunction){
830 ((TF1*)fInputFunction)->SetChisquare(0);
832 ((TF1*)fInputFunction)->SetNumberFitPoints(fNpoints);
833 }
834 }
835 return 1;
836 }
837 fParams=coef;
838 fParCovar=chol.Invert();
839
840 if (fInputFunction){
843 for (i=0; i<fNfunctions; i++){
844 e = TMath::Sqrt(fParCovar(i, i));
845 ((TF1*)fInputFunction)->SetParError(i, e);
846 }
847 if (!fObjectFit)
848 ((TF1*)fInputFunction)->SetChisquare(GetChisquare());
850 ((TF1*)fInputFunction)->SetNumberFitPoints(fNpoints);
851 }
852 }
853
854 //if parameters were fixed, change the design matrix back as it was before fixing
855 j = 0;
856 if (fNfixed>0){
857 for (i=0; i<fNfunctions; i++){
858 if (fFixedParams[i]){
859 for (ii=0; ii<i; ii++){
860 fDesign(ii, i) = fDesignTemp(ii, j);
862 }
863 for (ii=i; ii<fNfunctions; ii++){
864 fDesign(i, ii) = fDesignTemp(ii, j);
866 }
867 j++;
868 }
869 }
870 }
871 return 0;
872}
873
874////////////////////////////////////////////////////////////////////////////////
875///Fixes paramter `#ipar` at its current value.
876
878{
879 if (fParams.NonZeros()<1){
880 Error("FixParameter", "no value available to fix the parameter");
881 return;
882 }
883 if (ipar>fNfunctions || ipar<0){
884 Error("FixParameter", "illegal parameter value");
885 return;
886 }
887 if (fNfixed==fNfunctions) {
888 Error("FixParameter", "no free parameters left");
889 return;
890 }
891 if (!fFixedParams)
893 fFixedParams[ipar] = true;
894 fNfixed++;
895}
896
897////////////////////////////////////////////////////////////////////////////////
898///Fixes parameter `#ipar` at value `parvalue`.
899
901{
902 if (ipar>fNfunctions || ipar<0){
903 Error("FixParameter", "illegal parameter value");
904 return;
905 }
906 if (fNfixed==fNfunctions) {
907 Error("FixParameter", "no free parameters left");
908 return;
909 }
910 if(!fFixedParams)
912 fFixedParams[ipar] = true;
915 fParams(ipar) = parvalue;
916 fNfixed++;
917}
918
919////////////////////////////////////////////////////////////////////////////////
920///Releases parameter `#ipar`.
921
923{
924 if (ipar>fNfunctions || ipar<0){
925 Error("ReleaseParameter", "illegal parameter value");
926 return;
927 }
928 if (!fFixedParams[ipar]){
929 Warning("ReleaseParameter","This parameter is not fixed\n");
930 return;
931 } else {
932 fFixedParams[ipar] = false;
933 fNfixed--;
934 }
935}
936
937////////////////////////////////////////////////////////////////////////////////
938///Get the Atb vector - a vector, used for internal computations
939
941{
942 if (v.GetNoElements()!=fAtb.GetNoElements())
943 v.ResizeTo(fAtb.GetNoElements());
944 v = fAtb;
945}
946
947////////////////////////////////////////////////////////////////////////////////
948/// Get the Chisquare.
949
951{
952 if (fChisquare > 1e-16)
953 return fChisquare;
954 else {
955 Chisquare();
956 return fChisquare;
957 }
958}
959
960////////////////////////////////////////////////////////////////////////////////
961///Computes point-by-point confidence intervals for the fitted function
962///Parameters:
963///n - number of points
964///ndim - dimensions of points
965///x - points, at which to compute the intervals, for ndim > 1
966/// should be in order: (x0,y0, x1, y1, ... xn, yn)
967///ci - computed intervals are returned in this array
968///cl - confidence level, default=0.95
969///
970///NOTE, that this method can only be used when the fitting function inherits from a TF1,
971///so it's not possible when the fitting function was set as a string or as a pure TFormula
972
974{
975 if (fInputFunction){
976 Double_t *grad = new Double_t[fNfunctions];
978 Double_t c=0;
980 Double_t t = TMath::StudentQuantile(0.5 + cl/2, df);
982
983 for (Int_t ipoint=0; ipoint<n; ipoint++){
984 c=0;
985 ((TF1*)(fInputFunction))->GradientPar(x+ndim*ipoint, grad);
986 //multiply the covariance matrix by gradient
987 for (Int_t irow=0; irow<fNfunctions; irow++){
988 sum_vector[irow]=0;
989 for (Int_t icol=0; icol<fNfunctions; icol++){
991 }
992 }
993 for (Int_t i=0; i<fNfunctions; i++)
994 c+=grad[i]*sum_vector[i];
995 c=TMath::Sqrt(c);
996 ci[ipoint]=c*t*chidf;
997 }
998
999 delete [] grad;
1000 delete [] sum_vector;
1001 }
1002}
1003
1004////////////////////////////////////////////////////////////////////////////////
1005///Computes confidence intervals at level cl. Default is 0.95
1006///The TObject parameter can be a TGraphErrors, a TGraph2DErrors or a TH123.
1007///For Graphs, confidence intervals are computed for each point,
1008///the value of the graph at that point is set to the function value at that
1009///point, and the graph y-errors (or z-errors) are set to the value of
1010///the confidence interval at that point
1011///For Histograms, confidence intervals are computed for each bin center
1012///The bin content of this bin is then set to the function value at the bin
1013///center, and the bin error is set to the confidence interval value.
1014///Allowed combinations:
1015///Fitted object Passed object
1016///TGraph TGraphErrors, TH1
1017///TGraphErrors, AsymmErrors TGraphErrors, TH1
1018///TH1 TGraphErrors, TH1
1019///TGraph2D TGraph2DErrors, TH2
1020///TGraph2DErrors TGraph2DErrors, TH2
1021///TH2 TGraph2DErrors, TH2
1022///TH3 TH3
1023
1025{
1026 if (!fInputFunction) {
1027 Error("GetConfidenceIntervals", "The case of fitting not with a TFormula is not yet implemented");
1028 return;
1029 }
1030
1031 //TGraph//////////////////
1032
1033 if (obj->InheritsFrom(TGraph::Class())) {
1034 TGraph *gr = (TGraph*)obj;
1035 if (!gr->GetEY()){
1036 Error("GetConfidenceIntervals", "A TGraphErrors should be passed instead of a graph");
1037 return;
1038 }
1040 Error("GetConfidenceIntervals", "A TGraph2DErrors should be passed instead of a graph");
1041 return;
1042 }
1044 if (((TH1*)(fObjectFit))->GetDimension()>1){
1045 Error("GetConfidenceIntervals", "A TGraph2DErrors or a TH23 should be passed instead of a graph");
1046 return;
1047 }
1048 }
1049
1050 GetConfidenceIntervals(gr->GetN(),1,gr->GetX(), gr->GetEY(), cl);
1051 for (Int_t i=0; i<gr->GetN(); i++)
1052 gr->SetPoint(i, gr->GetX()[i], fInputFunction->Eval(gr->GetX()[i]));
1053 }
1054
1055 //TGraph2D///////////////
1056 else if (obj->InheritsFrom(TGraph2D::Class())) {
1057 TGraph2D *gr2 = (TGraph2D*)obj;
1058 if (!gr2->GetEZ()){
1059 Error("GetConfidenceIntervals", "A TGraph2DErrors should be passed instead of a TGraph2D");
1060 return;
1061 }
1063 Error("GetConfidenceIntervals", "A TGraphErrors should be passed instead of a TGraph2D");
1064 return;
1065 }
1067 if (((TH1*)(fObjectFit))->GetDimension()==1){
1068 Error("GetConfidenceIntervals", "A TGraphErrors or a TH1 should be passed instead of a graph");
1069 return;
1070 }
1071 }
1072 Double_t xy[2];
1073 Int_t np = gr2->GetN();
1074 Double_t *grad = new Double_t[fNfunctions];
1076 Double_t *x = gr2->GetX();
1077 Double_t *y = gr2->GetY();
1078 Double_t t = TMath::StudentQuantile(0.5 + cl/2, ((TF1*)(fInputFunction))->GetNDF());
1080 Double_t c = 0;
1081 for (Int_t ipoint=0; ipoint<np; ipoint++){
1082 c=0;
1083 xy[0]=x[ipoint];
1084 xy[1]=y[ipoint];
1085 ((TF1*)(fInputFunction))->GradientPar(xy, grad);
1086 for (Int_t irow=0; irow<fNfunctions; irow++){
1087 sum_vector[irow]=0;
1088 for (Int_t icol=0; icol<fNfunctions; icol++)
1090 }
1091 for (Int_t i=0; i<fNfunctions; i++)
1092 c+=grad[i]*sum_vector[i];
1093 c=TMath::Sqrt(c);
1094 gr2->SetPoint(ipoint, xy[0], xy[1], fInputFunction->EvalPar(xy));
1095 gr2->GetEZ()[ipoint]=c*t*chidf;
1096 }
1097 delete [] grad;
1098 delete [] sum_vector;
1099 }
1100
1101 //TH1////////////////////////
1102 else if (obj->InheritsFrom(TH1::Class())) {
1104 if (((TH1*)obj)->GetDimension()>1){
1105 Error("GetConfidenceIntervals", "Fitted graph and passed histogram have different number of dimensions");
1106 return;
1107 }
1108 }
1110 if (((TH1*)obj)->GetDimension()!=2){
1111 Error("GetConfidenceIntervals", "Fitted graph and passed histogram have different number of dimensions");
1112 return;
1113 }
1114 }
1116 if (((TH1*)(fObjectFit))->GetDimension()!=((TH1*)(obj))->GetDimension()){
1117 Error("GetConfidenceIntervals", "Fitted and passed histograms have different number of dimensions");
1118 return;
1119 }
1120 }
1121
1122
1123 TH1 *hfit = (TH1*)obj;
1124 Double_t *grad = new Double_t[fNfunctions];
1126 Double_t x[3];
1127 Int_t hxfirst = hfit->GetXaxis()->GetFirst();
1128 Int_t hxlast = hfit->GetXaxis()->GetLast();
1129 Int_t hyfirst = hfit->GetYaxis()->GetFirst();
1130 Int_t hylast = hfit->GetYaxis()->GetLast();
1131 Int_t hzfirst = hfit->GetZaxis()->GetFirst();
1132 Int_t hzlast = hfit->GetZaxis()->GetLast();
1133
1134 TAxis *xaxis = hfit->GetXaxis();
1135 TAxis *yaxis = hfit->GetYaxis();
1136 TAxis *zaxis = hfit->GetZaxis();
1137 Double_t t = TMath::StudentQuantile(0.5 + cl/2, ((TF1*)(fInputFunction))->GetNDF());
1139 Double_t c=0;
1140 for (Int_t binz=hzfirst; binz<=hzlast; binz++){
1141 x[2]=zaxis->GetBinCenter(binz);
1142 for (Int_t biny=hyfirst; biny<=hylast; biny++) {
1143 x[1]=yaxis->GetBinCenter(biny);
1144 for (Int_t binx=hxfirst; binx<=hxlast; binx++) {
1145 x[0]=xaxis->GetBinCenter(binx);
1146 ((TF1*)(fInputFunction))->GradientPar(x, grad);
1147 c=0;
1148 for (Int_t irow=0; irow<fNfunctions; irow++){
1149 sum_vector[irow]=0;
1150 for (Int_t icol=0; icol<fNfunctions; icol++)
1152 }
1153 for (Int_t i=0; i<fNfunctions; i++)
1154 c+=grad[i]*sum_vector[i];
1155 c=TMath::Sqrt(c);
1156 hfit->SetBinContent(binx, biny, binz, fInputFunction->EvalPar(x));
1157 hfit->SetBinError(binx, biny, binz, c*t*chidf);
1158 }
1159 }
1160 }
1161 delete [] grad;
1162 delete [] sum_vector;
1163 }
1164 else {
1165 Error("GetConfidenceIntervals", "This object type is not supported");
1166 return;
1167 }
1168}
1169
1170////////////////////////////////////////////////////////////////////////////////
1171///Returns covariance matrix
1172
1174{
1175 Double_t *p = const_cast<Double_t*>(fParCovar.GetMatrixArray());
1176 return p;
1177}
1178
1179////////////////////////////////////////////////////////////////////////////////
1180///Returns covariance matrix
1181
1183{
1184 if (matr.GetNrows()!=fNfunctions || matr.GetNcols()!=fNfunctions){
1186 }
1187 matr = fParCovar;
1188}
1189
1190////////////////////////////////////////////////////////////////////////////////
1191///Returns the internal design matrix
1192
1194{
1195 if (matr.GetNrows()!=fNfunctions || matr.GetNcols()!=fNfunctions){
1197 }
1198 matr = fDesign;
1199}
1200
1201////////////////////////////////////////////////////////////////////////////////
1202///Returns parameter errors
1203
1205{
1206 if (vpar.GetNoElements()!=fNfunctions) {
1207 vpar.ResizeTo(fNfunctions);
1208 }
1209 for (Int_t i=0; i<fNfunctions; i++)
1210 vpar(i) = TMath::Sqrt(fParCovar(i, i));
1211
1212}
1213
1214////////////////////////////////////////////////////////////////////////////////
1215///Returns parameter values
1216
1218{
1219 if (vpar.GetNoElements()!=fNfunctions) {
1220 vpar.ResizeTo(fNfunctions);
1221 }
1222 vpar=fParams;
1223}
1224
1225////////////////////////////////////////////////////////////////////////////////
1226///Returns the value and the name of the parameter `#ipar`
1227///NB: In the calling function the argument name must be set large enough
1228
1229Int_t TLinearFitter::GetParameter(Int_t ipar,char* name,Double_t& value,Double_t& /*verr*/,Double_t& /*vlow*/, Double_t& /*vhigh*/) const
1230{
1232 Error("GetParError", "illegal value of parameter");
1233 return 0;
1234 }
1235 value = fParams(ipar);
1236 if (fInputFunction)
1238 else
1239 name = nullptr;
1240 return 1;
1241}
1242
1243
1244////////////////////////////////////////////////////////////////////////////////
1245///Returns the error of parameter `#ipar`
1246
1248{
1250 Error("GetParError", "illegal value of parameter");
1251 return 0;
1252 }
1253
1254 return TMath::Sqrt(fParCovar(ipar, ipar));
1255}
1256
1257
1258////////////////////////////////////////////////////////////////////////////////
1259///Returns name of parameter `#ipar`
1260
1261const char *TLinearFitter::GetParName(Int_t ipar) const
1262{
1264 Error("GetParError", "illegal value of parameter");
1265 return nullptr;
1266 }
1267 if (fInputFunction)
1268 return fInputFunction->GetParName(ipar);
1269 return "";
1270}
1271
1272////////////////////////////////////////////////////////////////////////////////
1273///Returns the t-value for parameter `#ipar`
1274
1276{
1278 Error("GetParTValue", "illegal value of parameter");
1279 return 0;
1280 }
1281 if (!fTValues.NonZeros())
1283 return fTValues(ipar);
1284}
1285
1286////////////////////////////////////////////////////////////////////////////////
1287///Returns the significance of parameter `#ipar`
1288
1290{
1292 Error("GetParSignificance", "illegal value of parameter");
1293 return 0;
1294 }
1295 if (!fParSign.NonZeros())
1297 return fParSign(ipar);
1298}
1299
1300////////////////////////////////////////////////////////////////////////////////
1301///For robust lts fitting, returns the sample, on which the best fit was based
1302
1304{
1305 if (!fRobust){
1306 Error("GetFitSample", "there is no fit sample in ordinary least-squares fit");
1307 return;
1308 }
1309 for (Int_t i=0; i<fNpoints; i++)
1311
1312}
1313
1314////////////////////////////////////////////////////////////////////////////////
1315///Merge objects in list
1316
1318{
1319 if (!list) return -1;
1320 TIter next(list);
1321 TLinearFitter *lfit = nullptr;
1322 while ((lfit = (TLinearFitter*)next())) {
1323 if (!lfit->InheritsFrom(TLinearFitter::Class())) {
1324 Error("Add","Attempt to add object of class: %s to a %s",lfit->ClassName(),this->ClassName());
1325 return -1;
1326 }
1327 Add(lfit);
1328 }
1329 return 0;
1330}
1331////////////////////////////////////////////////////////////////////////////////
1332///set the basis functions in case the fitting function is not
1333/// set directly
1334/// The TLinearFitter will manage and delete the functions contained in the list
1335
1337{
1338 fFunctions = *(functions);
1340 int size = fFunctions.GetEntries();
1341
1343 //change the size of design matrix
1352 if (fFixedParams)
1353 delete [] fFixedParams;
1355 fDesign.Zero();
1356 fAtb.Zero();
1357 fDesignTemp.Zero();
1360 fAtbTemp.Zero();
1361 fAtbTemp2.Zero();
1362 fAtbTemp3.Zero();
1363 fY2Temp=0;
1364 fY2=0;
1365 for (int i=0; i<size; i++)
1366 fFixedParams[i]=false;
1367 fIsSet=kFALSE;
1368 fChisquare=0;
1369
1370}
1371
1372
1373////////////////////////////////////////////////////////////////////////////////
1374///set the number of dimensions
1375
1377{
1378 fNdim=ndim;
1379 fY.ResizeTo(ndim+1);
1380 fX.ResizeTo(ndim+1, ndim);
1381 fE.ResizeTo(ndim+1);
1382
1383 fNpoints=0;
1384 fIsSet=kFALSE;
1385}
1386
1387////////////////////////////////////////////////////////////////////////////////
1388///Additive parts should be separated by "++".
1389///Examples (ai are parameters to fit):
1390///1.fitting function: a0*x0 + a1*x1 + a2*x2
1391/// input formula "x[0]++x[1]++x[2]"
1392///2.TMath functions can be used:
1393/// fitting function: a0*TMath::Gaus(x, 0, 1) + a1*y
1394/// input formula: "TMath::Gaus(x, 0, 1)++y"
1395///fills the array of functions
1396
1397void TLinearFitter::SetFormula(const char *formula)
1398{
1399 Int_t size = 0, special = 0;
1400 Int_t i;
1401 //Int_t len = strlen(formula);
1402 if (fInputFunction)
1403 fInputFunction = nullptr;
1404 if (fFormula)
1405 delete [] fFormula;
1406 fFormulaSize = strlen(formula);
1407 fFormula = new char[fFormulaSize+1];
1408 strlcpy(fFormula, formula,fFormulaSize+1);
1409 fSpecial = 0;
1410 //in case of a hyperplane:
1411 char *fstring;
1412 fstring = (char *)strstr(fFormula, "hyp");
1413 if (fstring!=nullptr){
1414 // isHyper = kTRUE;
1415 fstring+=3;
1416 sscanf(fstring, "%d", &size);
1417 //+1 for the constant term
1418 size++;
1419 fSpecial=200+size;
1420 }
1421
1422 if (fSpecial==0) {
1423 //in case it's not a hyperplane
1425 sstring = sstring.ReplaceAll("++", 2, "|", 1);
1427
1428 //count the number of functions
1429
1430 TObjArray *oa = sstring.Tokenize("|");
1431
1432 //change the size of functions array and clear it
1433 if (!fFunctions.IsEmpty())
1434 fFunctions.Clear();
1435
1436 // do not own the functions in this case
1438
1439 fNfunctions = oa->GetEntriesFast();
1441
1442 //check if the old notation of xi is used somewhere instead of x[i]
1443 for (i=0; i<fNdim; i++){
1444 std::string pattern = "x" + std::to_string(i);
1445 std::string replacement = "x[" + std::to_string(i) + "]";
1446 sstring = sstring.ReplaceAll(pattern.c_str(), replacement.c_str());
1447 }
1448
1449 //fill the array of functions
1450 oa = sstring.Tokenize("|");
1451 for (i=0; i<fNfunctions; i++) {
1452 replaceformula = ((TObjString *)oa->UncheckedAt(i))->GetString();
1453 // look first if exists in the map
1454 TFormula * f = nullptr;
1455 auto elem = fgFormulaMap.find(replaceformula );
1456 if (elem != fgFormulaMap.end() )
1457 f = elem->second;
1458 else {
1459 // create a new formula and add in the static map
1460 f = new TFormula("f", replaceformula.Data());
1461 {
1464 }
1465 }
1466 if (!f) {
1467 Error("TLinearFitter", "f_linear not allocated");
1468 return;
1469 }
1470 special=f->GetNumber();
1471 fFunctions.Add(f);
1472 }
1473
1474 if ((fNfunctions==1)&&(special>299)&&(special<310)){
1475 //if fitting a polynomial
1476 size=special-299;
1477 fSpecial=100+size;
1478 } else
1480 oa->Delete();
1481 delete oa;
1482 }
1484 //change the size of design matrix
1493 if (fFixedParams)
1494 delete [] fFixedParams;
1496 fDesign.Zero();
1497 fAtb.Zero();
1498 fDesignTemp.Zero();
1501 fAtbTemp.Zero();
1502 fAtbTemp2.Zero();
1503 fAtbTemp3.Zero();
1504 fY2Temp=0;
1505 fY2=0;
1506 for (i=0; i<size; i++)
1507 fFixedParams[i]=false;
1508 fIsSet=kFALSE;
1509 fChisquare=0;
1510
1511}
1512
1513////////////////////////////////////////////////////////////////////////////////
1514///Set the fitting function.
1515
1517{
1519 fInputFunction=function;
1521 fSpecial = 0;
1523 if (!fFunctions.IsEmpty())
1525
1526 if ((special>299)&&(special<310)){
1527 //if fitting a polynomial
1528 size=special-299;
1529 fSpecial=100+size;
1530 } else
1532
1534 //change the size of design matrix
1539
1542
1545 //
1546 if (fFixedParams)
1547 delete [] fFixedParams;
1549 fDesign.Zero();
1550 fAtb.Zero();
1551 fDesignTemp.Zero();
1552 fAtbTemp.Zero();
1553
1556
1557 fAtbTemp2.Zero();
1558 fAtbTemp3.Zero();
1559 fY2Temp=0;
1560 fY2=0;
1561 for (Int_t i=0; i<size; i++)
1562 fFixedParams[i]=false;
1563 //check if any parameters are fixed (not for pure TFormula)
1564
1565 if (function->InheritsFrom(TF1::Class())){
1566 Double_t al,bl;
1567 for (Int_t i=0;i<fNfunctions;i++) {
1568 ((TF1*)function)->GetParLimits(i,al,bl);
1569 if (al*bl !=0 && al >= bl) {
1570 FixParameter(i, function->GetParameter(i));
1571 }
1572 }
1573 }
1574
1575 fIsSet=kFALSE;
1576 fChisquare=0;
1577
1578}
1579
1580////////////////////////////////////////////////////////////////////////////////
1581///Update the design matrix after the formula has been changed.
1582
1584{
1585 if (fStoreData) {
1586 for (Int_t i=0; i<fNpoints; i++) {
1587 AddToDesign(TMatrixDRow(fX, i).GetPtr(), fY(i), fE(i));
1588 }
1589 return true;
1590 } else
1591 return false;
1592
1593}
1594
1595////////////////////////////////////////////////////////////////////////////////
1596///To use in TGraph::Fit and TH1::Fit().
1597
1599{
1600 if (!strcmp(command, "FitGraph")){
1601 if (args) return GraphLinearFitter(args[0]);
1602 else return GraphLinearFitter(0);
1603 }
1604 if (!strcmp(command, "FitGraph2D")){
1605 if (args) return Graph2DLinearFitter(args[0]);
1606 else return Graph2DLinearFitter(0);
1607 }
1608 if (!strcmp(command, "FitMultiGraph")){
1609 if (args) return MultiGraphLinearFitter(args[0]);
1610 else return MultiGraphLinearFitter(0);
1611 }
1612 if (!strcmp(command, "FitHist")) return HistLinearFitter();
1613// if (!strcmp(command, "FitMultiGraph")) MultiGraphLinearFitter();
1614
1615 return 0;
1616}
1617
1618////////////////////////////////////////////////////////////////////////////////
1619/// Level = 3 (to be consistent with minuit) prints parameters and parameter
1620/// errors.
1621
1622void TLinearFitter::PrintResults(Int_t level, Double_t /*amin*/) const
1623{
1624 if (level==3){
1625 if (!fRobust){
1626 printf("Fitting results:\nParameters:\nNO.\t\tVALUE\t\tERROR\n");
1627 for (Int_t i=0; i<fNfunctions; i++){
1628 printf("%d\t%e\t%e\n", i, fParams(i), TMath::Sqrt(fParCovar(i, i)));
1629 }
1630 } else {
1631 printf("Fitting results:\nParameters:\nNO.\t\tVALUE\n");
1632 for (Int_t i=0; i<fNfunctions; i++){
1633 printf("%d\t%e\n", i, fParams(i));
1634 }
1635 }
1636 }
1637}
1638
1639////////////////////////////////////////////////////////////////////////////////
1640///Used in TGraph::Fit().
1641
1643{
1646 TF1 *f1=(TF1*)GetUserFunc();
1648
1649 //Int_t np=0;
1650 Double_t *x=grr->GetX();
1651 Double_t *y=grr->GetY();
1652 Double_t e;
1653
1654 Int_t fitResult = 0;
1655 //set the fitting formula
1656 SetDim(1);
1658
1659 if (fitOption.Robust){
1660 fRobust=kTRUE;
1662 }
1663 //put the points into the fitter
1664 Int_t n=grr->GetN();
1665 for (Int_t i=0; i<n; i++){
1666 if (!f1->IsInside(&x[i])) continue;
1667 e=grr->GetErrorY(i);
1668 if (e<0 || fitOption.W1)
1669 e=1;
1670 AddPoint(&x[i], y[i], e);
1671 }
1672
1673 if (fitOption.Robust){
1674 return EvalRobust(h);
1675 }
1676
1677 fitResult = Eval();
1678
1679 //calculate the precise chisquare
1680 if (!fitResult){
1681 if (!fitOption.Nochisq){
1682 Double_t temp, temp2, sumtotal=0;
1683 for (Int_t i=0; i<n; i++){
1684 if (!f1->IsInside(&x[i])) continue;
1685 temp=f1->Eval(x[i]);
1686 temp2=(y[i]-temp)*(y[i]-temp);
1687 e=grr->GetErrorY(i);
1688 if (e<0 || fitOption.W1)
1689 e=1;
1690 temp2/=(e*e);
1691
1692 sumtotal+=temp2;
1693 }
1696 }
1697 }
1698 return fitResult;
1699}
1700
1701////////////////////////////////////////////////////////////////////////////////
1702///Minimisation function for a TGraph2D
1703
1705{
1707
1709 TF2 *f2=(TF2*)GetUserFunc();
1710
1712 Int_t n = gr->GetN();
1713 Double_t *gx = gr->GetX();
1714 Double_t *gy = gr->GetY();
1715 Double_t *gz = gr->GetZ();
1716 Double_t x[2];
1717 Double_t z, e;
1718 Int_t fitResult=0;
1719 SetDim(2);
1720 SetFormula(f2->GetFormula());
1721
1722 if (fitOption.Robust){
1723 fRobust=kTRUE;
1725 }
1726
1727 for (Int_t bin=0;bin<n;bin++) {
1728 x[0] = gx[bin];
1729 x[1] = gy[bin];
1730 if (!f2->IsInside(x)) {
1731 continue;
1732 }
1733 z = gz[bin];
1734 e=gr->GetErrorZ(bin);
1735 if (e<0 || fitOption.W1)
1736 e=1;
1737 AddPoint(x, z, e);
1738 }
1739
1740 if (fitOption.Robust){
1741 return EvalRobust(h);
1742 }
1743
1744 fitResult = Eval();
1745
1746 if (!fitResult){
1747 if (!fitOption.Nochisq){
1748 //calculate the precise chisquare
1749 Double_t temp, temp2, sumtotal=0;
1750 for (Int_t bin=0; bin<n; bin++){
1751 x[0] = gx[bin];
1752 x[1] = gy[bin];
1753 if (!f2->IsInside(x)) {
1754 continue;
1755 }
1756 z = gz[bin];
1757
1758 temp=f2->Eval(x[0], x[1]);
1759 temp2=(z-temp)*(z-temp);
1760 e=gr->GetErrorZ(bin);
1761 if (e<0 || fitOption.W1)
1762 e=1;
1763 temp2/=(e*e);
1764
1765 sumtotal+=temp2;
1766 }
1769 }
1770 }
1771 return fitResult;
1772}
1773
1774////////////////////////////////////////////////////////////////////////////////
1775///Minimisation function for a TMultiGraph
1776
1778{
1779 Int_t n, i;
1780 Double_t *gx, *gy;
1781 Double_t e;
1783 TMultiGraph *mg = (TMultiGraph*)grFitter->GetObjectFit();
1784 TF1 *f1 = (TF1*)grFitter->GetUserFunc();
1785 Foption_t fitOption = grFitter->GetFitOption();
1786 Int_t fitResult=0;
1787 SetDim(1);
1788
1789 if (fitOption.Robust){
1790 fRobust=kTRUE;
1792 }
1794
1795 TGraph *gr;
1796 TIter next(mg->GetListOfGraphs());
1797 while ((gr = (TGraph*) next())) {
1798 n = gr->GetN();
1799 gx = gr->GetX();
1800 gy = gr->GetY();
1801 for (i=0; i<n; i++){
1802 if (!f1->IsInside(&gx[i])) continue;
1803 e=gr->GetErrorY(i);
1804 if (e<0 || fitOption.W1)
1805 e=1;
1806 AddPoint(&gx[i], gy[i], e);
1807 }
1808 }
1809
1810 if (fitOption.Robust){
1811 return EvalRobust(h);
1812 }
1813
1814 fitResult = Eval();
1815
1816 //calculate the chisquare
1817 if (!fitResult){
1818 if (!fitOption.Nochisq){
1819 Double_t temp, temp2, sumtotal=0;
1820 next.Reset();
1821 while((gr = (TGraph*)next())) {
1822 n = gr->GetN();
1823 gx = gr->GetX();
1824 gy = gr->GetY();
1825 for (i=0; i<n; i++){
1826 if (!f1->IsInside(&gx[i])) continue;
1827 temp=f1->Eval(gx[i]);
1828 temp2=(gy[i]-temp)*(gy[i]-temp);
1829 e=gr->GetErrorY(i);
1830 if (e<0 || fitOption.W1)
1831 e=1;
1832 temp2/=(e*e);
1833
1834 sumtotal+=temp2;
1835 }
1836
1837 }
1840 }
1841 }
1842 return fitResult;
1843}
1844
1845////////////////////////////////////////////////////////////////////////////////
1846/// Minimization function for H1s using a Chisquare method.
1847
1849{
1851 Double_t cu,eu;
1852 // Double_t dersum[100], grad[100];
1853 Double_t x[3];
1854 Int_t bin,binx,biny,binz;
1855 // Axis_t binlow, binup, binsize;
1856
1857 TH1 *hfit = (TH1*)GetObjectFit();
1858 TF1 *f1 = (TF1*)GetUserFunc();
1859 Int_t fitResult;
1861 //SetDim(hfit->GetDimension());
1862 SetDim(f1->GetNdim());
1864
1866 Int_t hxlast = GetXlast();
1868 Int_t hylast = GetYlast();
1870 Int_t hzlast = GetZlast();
1871 TAxis *xaxis = hfit->GetXaxis();
1872 TAxis *yaxis = hfit->GetYaxis();
1873 TAxis *zaxis = hfit->GetZaxis();
1874
1875 for (binz=hzfirst;binz<=hzlast;binz++) {
1876 x[2] = zaxis->GetBinCenter(binz);
1877 for (biny=hyfirst;biny<=hylast;biny++) {
1878 x[1] = yaxis->GetBinCenter(biny);
1879 for (binx=hxfirst;binx<=hxlast;binx++) {
1880 x[0] = xaxis->GetBinCenter(binx);
1881 if (!f1->IsInside(x)) continue;
1882 bin = hfit->GetBin(binx,biny,binz);
1883 cu = hfit->GetBinContent(bin);
1884 if (f1->GetNdim() < hfit->GetDimension()) {
1885 if (hfit->GetDimension() > 2) cu = x[2];
1886 else cu = x[1];
1887 }
1888 if (fitOption.W1) {
1889 if (fitOption.W1==1 && cu == 0) continue;
1890 eu = 1;
1891 } else {
1892 eu = hfit->GetBinError(bin);
1893 if (eu <= 0) continue;
1894 }
1895 AddPoint(x, cu, eu);
1896 }
1897 }
1898 }
1899
1900 fitResult = Eval();
1901
1902 if (!fitResult){
1903 if (!fitOption.Nochisq){
1904 Double_t temp, temp2, sumtotal=0;
1905 for (binz=hzfirst;binz<=hzlast;binz++) {
1906 x[2] = zaxis->GetBinCenter(binz);
1907 for (biny=hyfirst;biny<=hylast;biny++) {
1908 x[1] = yaxis->GetBinCenter(biny);
1909 for (binx=hxfirst;binx<=hxlast;binx++) {
1910 x[0] = xaxis->GetBinCenter(binx);
1911 if (!f1->IsInside(x)) continue;
1912 bin = hfit->GetBin(binx,biny,binz);
1913 cu = hfit->GetBinContent(bin);
1914
1915 if (fitOption.W1) {
1916 if (fitOption.W1==1 && cu == 0) continue;
1917 eu = 1;
1918 } else {
1919 eu = hfit->GetBinError(bin);
1920 if (eu <= 0) continue;
1921 }
1922 temp=f1->EvalPar(x);
1923 temp2=(cu-temp)*(cu-temp);
1924 temp2/=(eu*eu);
1925 sumtotal+=temp2;
1926 }
1927 }
1928 }
1929
1932 }
1933 }
1934 return fitResult;
1935}
1936
1937////////////////////////////////////////////////////////////////////////////////
1938
1940{
1941 if (R__b.IsReading()) {
1943 R__b.ReadClassBuffer(TLinearFitter::Class(),this);
1944 if (old_matr_size < fNfunctions) {
1947
1950
1953 }
1954 } else {
1955 if (fAtb.NonZeros()==0) AddTempMatrices();
1956 R__b.WriteClassBuffer(TLinearFitter::Class(),this);
1957 }
1958}
1959
1960////////////////////////////////////////////////////////////////////////////////
1961///Finds the parameters of the fitted function in case data contains
1962///outliers.
1963///Parameter h stands for the minimal fraction of good points in the
1964///dataset (h < 1, i.e. for 70% of good points take h=0.7).
1965///The default value of h*Npoints is (Npoints + Nparameters+1)/2
1966///If the user provides a value of h smaller than above, default is taken
1967///See class description for the algorithm details
1968
1970{
1971 fRobust = kTRUE;
1972 Double_t kEps = 1e-13;
1973 Int_t nmini = 300;
1974 Int_t i, j, maxind=0, k, k1 = 500;
1975 Int_t nbest = 10;
1976 Double_t chi2 = -1;
1977
1978 if (fFunctions.IsEmpty()&&(!fInputFunction)&&(fSpecial<=200)){
1979 Error("TLinearFitter::EvalRobust", "The formula hasn't been set");
1980 return 1;
1981 }
1982
1983
1985 for (i=0; i<nbest; i++)
1986 bestchi2[i]=1e30;
1987
1989
1990 if (h>0.000001 && h<1 && fNpoints*h > hdef)
1991 fH = Int_t(fNpoints*h);
1992 else {
1993 // print message only when h is not zero
1994 if (h>0) Warning("Fitting:", "illegal value of H, default is taken, h = %3.2f",double(hdef)/fNpoints);
1995 fH=hdef;
1996 }
1997
2001
2002 Int_t *index = new Int_t[fNpoints];
2004
2005 if (fNpoints < 2*nmini) {
2006 //when number of cases is small
2007
2008 //to store the best coefficients (columnwise)
2010 for (k = 0; k < k1; k++) {
2012 chi2 = CStep(1, fH, residuals,index, index, -1, -1);
2013 chi2 = CStep(2, fH, residuals,index, index, -1, -1);
2015 if (chi2 < bestchi2[maxind]) {
2016 bestchi2[maxind] = chi2;
2017 for (i=0; i<fNfunctions; i++)
2018 cstock(i, maxind) = fParams(i);
2019 }
2020 }
2021
2022 //for the nbest best results, perform CSteps until convergence
2023 Int_t *bestindex = new Int_t[fH];
2025 for (i=0; i<nbest; i++) {
2026 for (j=0; j<fNfunctions; j++)
2027 fParams(j) = cstock(j, i);
2028 chi2 = 1;
2029 while (chi2 > kEps) {
2030 chi2 = CStep(2, fH, residuals,index, index, -1, -1);
2031 if (TMath::Abs(chi2 - bestchi2[i]) < kEps)
2032 break;
2033 else
2034 bestchi2[i] = chi2;
2035 }
2037 if (chi2 <= currentbest + kEps) {
2038 for (j=0; j<fH; j++){
2039 bestindex[j]=index[j];
2040 }
2041 maxind = i;
2042 }
2043 for (j=0; j<fNfunctions; j++)
2044 cstock(j, i) = fParams(j);
2045 }
2046 //report the result with the lowest chisquare
2047 for (j=0; j<fNfunctions; j++)
2048 fParams(j) = cstock(j, maxind);
2050 for (j=0; j<fH; j++){
2051 //printf("bestindex[%d]=%d\n", j, bestindex[j]);
2053 }
2055 ((TF1*)fInputFunction)->SetChisquare(bestchi2[maxind]);
2056 ((TF1*)fInputFunction)->SetNumberFitPoints(fH);
2057 ((TF1*)fInputFunction)->SetNDF(fH-fNfunctions);
2058 }
2059 delete [] index;
2060 delete [] bestindex;
2061 delete [] residuals;
2062 delete [] bestchi2;
2063 return 0;
2064 }
2065 //if n is large, the dataset should be partitioned
2066 Int_t indsubdat[5];
2067 for (i=0; i<5; i++)
2068 indsubdat[i] = 0;
2069
2071 Int_t hsub;
2072
2074
2075 Int_t *subdat = new Int_t[sum]; //to store the indices of selected cases
2077
2079 Int_t *beststock = new Int_t[nbest];
2080 Int_t i_start = 0;
2081 Int_t i_end = indsubdat[0];
2082 Int_t k2 = Int_t(k1/nsub);
2083 for (Int_t kgroup = 0; kgroup < nsub; kgroup++) {
2084
2086 for (i=0; i<nbest; i++)
2087 bestchi2[i] = 1e16;
2088 for (k=0; k<k2; k++) {
2093 if (chi2 < bestchi2[maxind]){
2094 for (i=0; i<fNfunctions; i++)
2095 cstockbig(i, nbest*kgroup + maxind) = fParams(i);
2096 bestchi2[maxind] = chi2;
2097 }
2098 }
2099 if (kgroup != nsub - 1){
2101 i_end += indsubdat[kgroup+1];
2102 }
2103 }
2104
2105 for (i=0; i<nbest; i++)
2106 bestchi2[i] = 1e30;
2107 //on the pooled subset
2109 for (k=0; k<nbest*5; k++) {
2110 for (i=0; i<fNfunctions; i++)
2111 fParams(i)=cstockbig(i, k);
2112 chi2 = CStep(1, hsub2, residuals, index, subdat, 0, sum);
2113 chi2 = CStep(2, hsub2, residuals, index, subdat, 0, sum);
2115 if (chi2 < bestchi2[maxind]){
2116 beststock[maxind] = k;
2117 bestchi2[maxind] = chi2;
2118 }
2119 }
2120
2121 //now the array beststock keeps indices of 10 best candidates in cstockbig matrix
2122 for (k=0; k<nbest; k++) {
2123 for (i=0; i<fNfunctions; i++)
2124 fParams(i) = cstockbig(i, beststock[k]);
2125 chi2 = CStep(1, fH, residuals, index, index, -1, -1);
2126 chi2 = CStep(2, fH, residuals, index, index, -1, -1);
2127 bestchi2[k] = chi2;
2128 }
2129
2131 for (i=0; i<fNfunctions; i++)
2133
2134 chi2 = 1;
2135 while (chi2 > kEps) {
2136 chi2 = CStep(2, fH, residuals, index, index, -1, -1);
2138 break;
2139 else
2140 bestchi2[maxind] = chi2;
2141 }
2142
2144 for (j=0; j<fH; j++)
2146 if (fInputFunction){
2147 ((TF1*)fInputFunction)->SetChisquare(bestchi2[maxind]);
2148 ((TF1*)fInputFunction)->SetNumberFitPoints(fH);
2149 ((TF1*)fInputFunction)->SetNDF(fH-fNfunctions);
2150 }
2151
2152 delete [] subdat;
2153 delete [] beststock;
2154 delete [] bestchi2;
2155 delete [] residuals;
2156 delete [] index;
2157
2158 return 0;
2159}
2160
2161////////////////////////////////////////////////////////////////////////////////
2162///Creates a p-subset to start
2163///ntotal - total number of points from which the subset is chosen
2164
2166{
2167 Int_t i, j;
2169 Int_t nindex=0;
2170 Int_t num;
2171 for(i=0; i<ntotal; i++)
2172 index[i] = ntotal+1;
2173
2174 TRandom2 r;
2175 //create a p-subset
2176 for (i=0; i<fNfunctions; i++) {
2177 num=Int_t(r.Uniform(0, 1)*(ntotal-1));
2178 if (i>0){
2179 for(j=0; j<=i-1; j++) {
2180 if(index[j]==num)
2181 repeat = kTRUE;
2182 }
2183 }
2184 if(repeat==kTRUE) {
2185 i--;
2186 repeat = kFALSE;
2187 } else {
2188 index[i] = num;
2189 nindex++;
2190 }
2191 }
2192
2193 //compute the coefficients of a hyperplane through the p-subset
2194 fDesign.Zero();
2195 fAtb.Zero();
2196 for (i=0; i<fNfunctions; i++){
2197 AddToDesign(TMatrixDRow(fX, index[i]).GetPtr(), fY(index[i]), fE(index[i]));
2198 }
2199 Bool_t ok;
2200
2201 ok = Linf();
2202
2203 //if the chosen points don't define a hyperplane, add more
2204 while (!ok && (nindex < h)) {
2205 repeat=kFALSE;
2206 do{
2207 num=Int_t(r.Uniform(0,1)*(ntotal-1));
2208 repeat=kFALSE;
2209 for(i=0; i<nindex; i++) {
2210 if(index[i]==num) {
2211 repeat=kTRUE;
2212 break;
2213 }
2214 }
2215 } while(repeat==kTRUE);
2216
2217 index[nindex] = num;
2218 nindex++;
2219 //check if the system is of full rank now
2220 AddToDesign(TMatrixDRow(fX, index[nindex-1]).GetPtr(), fY(index[nindex-1]), fE(index[nindex-1]));
2221 ok = Linf();
2222 }
2223}
2224
2225////////////////////////////////////////////////////////////////////////////////
2226///The CStep procedure, as described in the article
2227
2229{
2231
2232 Int_t i, j, itemp, n;
2233 Double_t func;
2234 Double_t val[100];
2235 Int_t npar;
2236 if (start > -1) {
2237 n = end - start;
2238 for (i=0; i<n; i++) {
2239 func = 0;
2240 itemp = subdat[start+i];
2241 if (fInputFunction){
2243 func=fInputFunction->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
2244 } else {
2245 func=0;
2246 if ((fSpecial>100)&&(fSpecial<200)){
2247 npar = fSpecial-100;
2248 val[0] = 1;
2249 for (j=1; j<npar; j++)
2250 val[j] = val[j-1]*fX(itemp, 0);
2251 for (j=0; j<npar; j++)
2252 func += fParams(j)*val[j];
2253 } else if (fSpecial>200) {
2254 //hyperplane case
2255 npar = fSpecial-201;
2256 func+=fParams(0);
2257 for (j=0; j<npar; j++)
2258 func += fParams(j+1)*fX(itemp, j);
2259 } else {
2260 // general case
2261 for (j=0; j<fNfunctions; j++) {
2262 TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
2263 val[j] = f1->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
2264 func += fParams(j)*val[j];
2265 }
2266 }
2267 }
2268 residuals[i] = (fY(itemp) - func)*(fY(itemp) - func)/(fE(i)*fE(i));
2269 }
2270 } else {
2271 n=fNpoints;
2272 for (i=0; i<fNpoints; i++) {
2273 func = 0;
2274 if (fInputFunction){
2276 func=fInputFunction->EvalPar(TMatrixDRow(fX, i).GetPtr());
2277 } else {
2278 func=0;
2279 if ((fSpecial>100)&&(fSpecial<200)){
2280 npar = fSpecial-100;
2281 val[0] = 1;
2282 for (j=1; j<npar; j++)
2283 val[j] = val[j-1]*fX(i, 0);
2284 for (j=0; j<npar; j++)
2285 func += fParams(j)*val[j];
2286 } else if (fSpecial>200) {
2287 //hyperplane case
2288 npar = fSpecial-201;
2289 func+=fParams(0);
2290 for (j=0; j<npar; j++)
2291 func += fParams(j+1)*fX(i, j);
2292 } else {
2293 for (j=0; j<fNfunctions; j++) {
2294 TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
2295 val[j] = f1->EvalPar(TMatrixDRow(fX, i).GetPtr());
2296 func += fParams(j)*val[j];
2297 }
2298 }
2299 }
2300 residuals[i] = (fY(i) - func)*(fY(i) - func)/(fE(i)*fE(i));
2301 }
2302 }
2303 //take h with smallest residuals
2305 //add them to the design matrix
2306 fDesign.Zero();
2307 fAtb.Zero();
2308 for (i=0; i<h; i++)
2309 AddToDesign(TMatrixDRow(fX, index[i]).GetPtr(), fY(index[i]), fE(index[i]));
2310
2311 Linf();
2312
2313 //don't calculate the chisquare at the 1st cstep
2314 if (step==1) return 0;
2315 Double_t sum=0;
2316
2317
2318 if (start > -1) {
2319 for (i=0; i<h; i++) {
2320 itemp = subdat[start+index[i]];
2321 if (fInputFunction){
2323 func=fInputFunction->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
2324 } else {
2325 func=0;
2326 if ((fSpecial>100)&&(fSpecial<200)){
2327 npar = fSpecial-100;
2328 val[0] = 1;
2329 for (j=1; j<npar; j++)
2330 val[j] = val[j-1]*fX(itemp, 0);
2331 for (j=0; j<npar; j++)
2332 func += fParams(j)*val[j];
2333 } else if (fSpecial>200) {
2334 //hyperplane case
2335 npar = fSpecial-201;
2336 func+=fParams(0);
2337 for (j=0; j<npar; j++)
2338 func += fParams(j+1)*fX(itemp, j);
2339 } else {
2340 for (j=0; j<fNfunctions; j++) {
2341 TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
2342 val[j] = f1->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
2343 func += fParams(j)*val[j];
2344 }
2345 }
2346 }
2347 sum+=(fY(itemp)-func)*(fY(itemp)-func)/(fE(itemp)*fE(itemp));
2348 }
2349 } else {
2350 for (i=0; i<h; i++) {
2351 if (fInputFunction){
2353 func=fInputFunction->EvalPar(TMatrixDRow(fX, index[i]).GetPtr());
2354 } else {
2355 func=0;
2356 if ((fSpecial>100)&&(fSpecial<200)){
2357 npar = fSpecial-100;
2358 val[0] = 1;
2359 for (j=1; j<npar; j++)
2360 val[j] = val[j-1]*fX(index[i], 0);
2361 for (j=0; j<npar; j++)
2362 func += fParams(j)*val[j];
2363 } else {
2364 if (fSpecial>200) {
2365 //hyperplane case
2366 npar = fSpecial-201;
2367 func+=fParams(0);
2368 for (j=0; j<npar; j++)
2369 func += fParams(j+1)*fX(index[i], j);
2370 } else {
2371 for (j=0; j<fNfunctions; j++) {
2372 TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
2373 val[j] = f1->EvalPar(TMatrixDRow(fX, index[i]).GetPtr());
2374 func += fParams(j)*val[j];
2375 }
2376 }
2377 }
2378 }
2379
2380 sum+=(fY(index[i])-func)*(fY(index[i])-func)/(fE(index[i])*fE(index[i]));
2381 }
2382 }
2383
2384 return sum;
2385}
2386
2387////////////////////////////////////////////////////////////////////////////////
2388
2390{
2391 //currently without the intercept term
2397 fDesignTemp.Zero();
2400 fAtb+=fAtbTemp;
2401 fAtbTemp3.Zero();
2402 fAtbTemp2.Zero();
2403 fAtbTemp.Zero();
2404
2405 fY2+=fY2Temp;
2406 fY2Temp=0;
2407
2408
2410 TVectorD temp(fNfunctions);
2411 Bool_t ok;
2412 temp = chol.Solve(fAtb, ok);
2413 if (!ok){
2414 Error("Linf","Matrix inversion failed");
2415 //fDesign.Print();
2416 fParams.Zero();
2417 return kFALSE;
2418 }
2419 fParams = temp;
2420 return ok;
2421}
2422
2423////////////////////////////////////////////////////////////////////////////////
2424///divides the elements into approximately equal subgroups
2425///number of elements in each subgroup is stored in indsubdat
2426///number of subgroups is returned
2427
2429{
2430 Int_t nsub;
2431
2432 if ((fNpoints>=2*nmini) && (fNpoints<=(3*nmini-1))) {
2433 if (fNpoints%2==1){
2434 indsubdat[0]=Int_t(fNpoints*0.5);
2435 indsubdat[1]=Int_t(fNpoints*0.5)+1;
2436 } else
2438 nsub=2;
2439 }
2440 else{
2441 if((fNpoints>=3*nmini) && (fNpoints<(4*nmini -1))) {
2442 if(fNpoints%3==0){
2444 } else {
2445 indsubdat[0]=Int_t(fNpoints/3);
2446 indsubdat[1]=Int_t(fNpoints/3)+1;
2447 if (fNpoints%3==1) indsubdat[2]=Int_t(fNpoints/3);
2448 else indsubdat[2]=Int_t(fNpoints/3)+1;
2449 }
2450 nsub=3;
2451 }
2452 else{
2453 if((fNpoints>=4*nmini)&&(fNpoints<=(5*nmini-1))){
2454 if (fNpoints%4==0) indsubdat[0]=indsubdat[1]=indsubdat[2]=indsubdat[3]=Int_t(fNpoints/4);
2455 else {
2456 indsubdat[0]=Int_t(fNpoints/4);
2457 indsubdat[1]=Int_t(fNpoints/4)+1;
2458 if(fNpoints%4==1) indsubdat[2]=indsubdat[3]=Int_t(fNpoints/4);
2459 if(fNpoints%4==2) {
2460 indsubdat[2]=Int_t(fNpoints/4)+1;
2461 indsubdat[3]=Int_t(fNpoints/4);
2462 }
2463 if(fNpoints%4==3) indsubdat[2]=indsubdat[3]=Int_t(fNpoints/4)+1;
2464 }
2465 nsub=4;
2466 } else {
2467 for(Int_t i=0; i<5; i++)
2468 indsubdat[i]=nmini;
2469 nsub=5;
2470 }
2471 }
2472 }
2473 return nsub;
2474}
2475
2476////////////////////////////////////////////////////////////////////////////////
2477///Draws ngroup nonoverlapping subdatasets out of a dataset of size n
2478///such that the selected case numbers are uniformly distributed from 1 to n
2479
2481{
2482 Int_t jndex = 0;
2483 Int_t nrand;
2484 Int_t i, k, m, j;
2485 Int_t ngroup=0;
2486 for (i=0; i<5; i++) {
2487 if (indsubdat[i]!=0)
2488 ngroup++;
2489 }
2490 TRandom r;
2491 for (k=1; k<=ngroup; k++) {
2492 for (m=1; m<=indsubdat[k-1]; m++) {
2493 nrand = Int_t(r.Uniform(0, 1) * (fNpoints-jndex)) + 1;
2494 jndex++;
2495 if (jndex==1) {
2496 subdat[0] = nrand;
2497 } else {
2498 subdat[jndex-1] = nrand + jndex - 2;
2499 for (i=1; i<=jndex-1; i++) {
2500 if(subdat[i-1] > nrand+i-2) {
2501 for(j=jndex; j>=i+1; j--) {
2502 subdat[j-1] = subdat[j-2];
2503 }
2504 subdat[i-1] = nrand+i-2;
2505 break; //breaking the loop for(i=1...
2506 }
2507 }
2508 }
2509 }
2510 }
2511}
2512
2513
#define f(i)
Definition RSha256.hxx:104
#define c(i)
Definition RSha256.hxx:101
#define h(i)
Definition RSha256.hxx:106
#define e(i)
Definition RSha256.hxx:103
static void update(gsl_integration_workspace *workspace, double a1, double b1, double area1, double error1, double a2, double b2, double area2, double error2)
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
int Int_t
Signed integer 4 bytes (int)
Definition RtypesCore.h:59
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
const char Option_t
Option string (const char)
Definition RtypesCore.h:80
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
#define R__ASSERT(e)
Checks condition e and reports a fatal error if it's false.
Definition TError.h:125
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t option
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t np
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
Option_t Option_t TPoint xy
char name[80]
Definition TGX11.cxx:110
TMatrixTRow< Double_t > TMatrixDRow
R__EXTERN TVirtualMutex * gROOTMutex
Definition TROOT.h:63
#define R__LOCKGUARD(mutex)
const_iterator end() const
Class to manage histogram axis.
Definition TAxis.h:32
Container of bits.
Definition TBits.h:26
void Clear(Option_t *option="") override
Clear the value.
Definition TBits.cxx:83
Bool_t TestBitNumber(UInt_t bitnumber) const
Definition TBits.h:222
void SetBitNumber(UInt_t bitnumber, Bool_t value=kTRUE)
Definition TBits.h:206
Buffer base class used for serializing objects.
Definition TBuffer.h:43
Collection abstract base class.
Definition TCollection.h:65
virtual void SetOwner(Bool_t enable=kTRUE)
Set whether this collection is the owner (enable==true) of its content.
Cholesky Decomposition class.
Definition TDecompChol.h:25
1-Dim function class
Definition TF1.h:182
static TClass * Class()
virtual void SetChisquare(Double_t chi2)
Definition TF1.h:596
virtual TFormula * GetFormula()
Definition TF1.h:433
virtual Double_t EvalPar(const Double_t *x, const Double_t *params=nullptr)
Evaluate function with given coordinates and parameters.
Definition TF1.cxx:1475
virtual Double_t Eval(Double_t x, Double_t y=0, Double_t z=0, Double_t t=0) const
Evaluate this function.
Definition TF1.cxx:1446
virtual Bool_t IsInside(const Double_t *x) const
return kTRUE if the point is inside the function range
Definition TF1.h:582
virtual Int_t GetNdim() const
Definition TF1.h:465
A 2-Dim function with parameters.
Definition TF2.h:29
Bool_t IsInside(const Double_t *x) const override
Return kTRUE is the point is inside the function range.
Definition TF2.cxx:694
The Formula class.
Definition TFormula.h:89
const TObject * GetLinearPart(Int_t i) const
Return linear part.
Double_t Eval(Args... args) const
Set first 1, 2, 3 or 4 variables (e.g.
Definition TFormula.h:325
static TClass * Class()
const char * GetParName(Int_t ipar) const
Return parameter name given by integer.
Int_t GetNumber() const
Definition TFormula.h:262
void SetParameters(const Double_t *params)
Set a vector of parameters value.
Int_t GetNpar() const
Definition TFormula.h:261
Double_t EvalPar(const Double_t *x, const Double_t *params=nullptr) const
Graphics object made of three arrays X, Y and Z with the same number of points each.
Definition TGraph2D.h:41
static TClass * Class()
Double_t GetErrorY(Int_t bin) const override
It returns the error along Y at point i.
Double_t * GetEY() const override
A TGraph is an object made of two arrays X and Y with npoints each.
Definition TGraph.h:41
static TClass * Class()
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition TGraph.cxx:2290
Double_t * GetY() const
Definition TGraph.h:139
Int_t GetN() const
Definition TGraph.h:131
Double_t * GetX() const
Definition TGraph.h:138
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:109
static TClass * Class()
void Reset()
virtual void AddTempMatrices()
TMatrixDSym fDesignTemp2
temporary matrix, used for num.stability
Int_t GraphLinearFitter(Double_t h)
Used in TGraph::Fit().
~TLinearFitter() override
Linear fitter cleanup.
TMatrixDSym fDesignTemp
virtual Double_t GetChisquare()
Get the Chisquare.
Int_t Partition(Int_t nmini, Int_t *indsubdat)
divides the elements into approximately equal subgroups number of elements in each subgroup is stored...
virtual void GetErrors(TVectorD &vpar)
Returns parameter errors.
Double_t CStep(Int_t step, Int_t h, Double_t *residuals, Int_t *index, Int_t *subdat, Int_t start, Int_t end)
The CStep procedure, as described in the article.
virtual Int_t Merge(TCollection *list)
Merge objects in list.
Double_t fChisquare
void ComputeTValues()
Computes parameters' t-values and significance.
Bool_t * fFixedParams
TVectorD fParSign
void Clear(Option_t *option="") override
Clears everything. Used in TH1::Fit and TGraph::Fit().
TLinearFitter()
default c-tor, input data is stored If you don't want to store the input data, run the function Store...
const char * GetParName(Int_t ipar) const override
Returns name of parameter #ipar
Int_t MultiGraphLinearFitter(Double_t h)
Minimisation function for a TMultiGraph.
virtual Double_t GetParSignificance(Int_t ipar)
Returns the significance of parameter #ipar
TMatrixDSym fDesign
virtual Int_t Eval()
Perform the fit and evaluate the parameters Returns 0 if the fit is ok, 1 if there are errors.
TVectorD fAtbTemp2
temporary vector, used for num.stability
void PrintResults(Int_t level, Double_t amin=0) const override
Level = 3 (to be consistent with minuit) prints parameters and parameter errors.
Int_t HistLinearFitter()
Minimization function for H1s using a Chisquare method.
Int_t Graph2DLinearFitter(Double_t h)
Minimisation function for a TGraph2D.
static TClass * Class()
TVectorD fAtbTemp
virtual void ClearPoints()
To be used when different sets of points are fitted with the same formula.
TMatrixDSym fDesignTemp3
TObjArray fFunctions
map of basis functions and formula
Int_t ExecuteCommand(const char *command, Double_t *args, Int_t nargs) override
To use in TGraph::Fit and TH1::Fit().
virtual void GetFitSample(TBits &bits)
For robust lts fitting, returns the sample, on which the best fit was based.
virtual void Add(TLinearFitter *tlf)
Add another linear fitter to this linear fitter.
virtual void GetDesignMatrix(TMatrixD &matr)
Returns the internal design matrix.
void Streamer(TBuffer &) override
Stream an object of class TObject.
virtual void AssignData(Int_t npoints, Int_t xncols, Double_t *x, Double_t *y, Double_t *e=nullptr)
This function is to use when you already have all the data in arrays and don't want to copy them into...
virtual void GetParameters(TVectorD &vpar)
Returns parameter values.
void RDraw(Int_t *subdat, Int_t *indsubdat)
Draws ngroup nonoverlapping subdatasets out of a dataset of size n such that the selected case number...
static std::map< TString, TFormula * > fgFormulaMap
virtual void SetDim(Int_t n)
set the number of dimensions
TFormula * fInputFunction
void GetConfidenceIntervals(Int_t n, Int_t ndim, const Double_t *x, Double_t *ci, Double_t cl=0.95) override
Computes point-by-point confidence intervals for the fitted function Parameters: n - number of points...
TMatrixD fX
temporary variable used for num.stability
virtual Bool_t UpdateMatrix()
Update the design matrix after the formula has been changed.
virtual void GetAtbVector(TVectorD &v)
Get the Atb vector - a vector, used for internal computations.
virtual void Chisquare()
Calculates the chisquare.
virtual void SetBasisFunctions(TObjArray *functions)
set the basis functions in case the fitting function is not set directly The TLinearFitter will manag...
Int_t fNpoints
temporary
Double_t fVal[1000]
virtual Int_t EvalRobust(Double_t h=-1)
Finds the parameters of the fitted function in case data contains outliers.
TVectorD fTValues
Double_t GetParError(Int_t ipar) const override
Returns the error of parameter #ipar
void AddToDesign(Double_t *x, Double_t y, Double_t e)
Add a point to the AtA matrix and to the Atb vector.
TVectorD fAtbTemp3
void FixParameter(Int_t ipar) override
Fixes paramter #ipar at its current value.
void CreateSubset(Int_t ntotal, Int_t h, Int_t *index)
Creates a p-subset to start ntotal - total number of points from which the subset is chosen.
virtual void SetFormula(const char *formula)
Additive parts should be separated by "++".
void ReleaseParameter(Int_t ipar) override
Releases parameter #ipar.
virtual void AddPoint(Double_t *x, Double_t y, Double_t e=1)
Adds 1 point to the fitter.
TMatrixDSym fParCovar
virtual Double_t GetParTValue(Int_t ipar)
Returns the t-value for parameter #ipar
Double_t * GetCovarianceMatrix() const override
Returns covariance matrix.
Double_t GetParameter(Int_t ipar) const override
virtual void StoreData(Bool_t store)
Int_t GetNrows() const
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
Int_t GetNcols() const
TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1) override
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
void Clear(Option_t *="") override
Definition TMatrixTSym.h:96
const Element * GetMatrixArray() const override
TMatrixT< Element > & Use(Int_t row_lwb, Int_t row_upb, Int_t col_lwb, Int_t col_upb, Element *data)
Use the array data to fill the matrix ([row_lwb..row_upb] x [col_lwb..col_upb])
void Clear(Option_t *="") override
Definition TMatrixT.h:123
const Element * GetMatrixArray() const override
Definition TMatrixT.h:228
TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1) override
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
A TMultiGraph is a collection of TGraph (or derived) objects.
Definition TMultiGraph.h:34
TList * GetListOfGraphs() const
Definition TMultiGraph.h:68
const char * GetName() const override
Returns name of object.
Definition TNamed.h:49
An array of TObjects.
Definition TObjArray.h:31
virtual void Expand(Int_t newSize)
Expand or shrink the array to newSize elements.
void Clear(Option_t *option="") override
Remove all objects from the array.
Int_t GetEntries() const override
Return the number of objects in array (i.e.
void Delete(Option_t *option="") override
Remove all objects from the array AND delete all heap based objects.
TObject * UncheckedAt(Int_t i) const
Definition TObjArray.h:84
Bool_t IsEmpty() const override
Definition TObjArray.h:65
void Add(TObject *obj) override
Definition TObjArray.h:68
Collectable string class.
Definition TObjString.h:28
Mother of all ROOT objects.
Definition TObject.h:41
virtual const char * GetName() const
Returns name of object.
Definition TObject.cxx:457
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:1057
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition TObject.cxx:543
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1071
virtual TClass * IsA() const
Definition TObject.h:246
Random number generator class based on the maximally quidistributed combined Tausworthe generator by ...
Definition TRandom2.h:27
This is the base class for the ROOT Random number generators.
Definition TRandom.h:27
Basic string class.
Definition TString.h:138
TVectorT< Element > & Zero()
Set vector elements to zero.
Definition TVectorT.cxx:452
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition TVectorT.cxx:293
TVectorT< Element > & Use(Int_t lwb, Int_t upb, Element *data)
Use the array data to fill the vector lwb..upb].
Definition TVectorT.cxx:348
void Clear(Option_t *="") override
Definition TVectorT.h:194
Int_t NonZeros() const
Compute the number of elements != 0.0.
Definition TVectorT.cxx:619
Int_t GetNoElements() const
Definition TVectorT.h:76
Element * GetMatrixArray()
Definition TVectorT.h:78
Abstract Base Class for Fitting.
virtual Int_t GetXlast() const
virtual Int_t GetYfirst() const
virtual TObject * GetObjectFit() const
virtual Foption_t GetFitOption() const
virtual Int_t GetZfirst() const
virtual Int_t GetZlast() const
virtual Int_t GetXfirst() const
TObject * fObjectFit
Pointer to object being fitted.
virtual Int_t GetYlast() const
static TVirtualFitter * GetFitter()
static: return the current Fitter
virtual TObject * GetUserFunc() const
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
TGraphErrors * gr
Definition legend1.C:25
TF1 * f1
Definition legend1.C:11
Element KOrdStat(Size n, const Element *a, Size k, Size *work=0)
Returns k_th order statistic of the array a of size n (k_th smallest element out of n elements).
Definition TMath.h:1568
Long64_t LocMin(Long64_t n, const T *a)
Returns index of array with the minimum element.
Definition TMath.h:993
T MinElement(Long64_t n, const T *a)
Returns minimum of array a of length n.
Definition TMath.h:971
Long64_t LocMax(Long64_t n, const T *a)
Returns index of array with the maximum element.
Definition TMath.h:1099
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673
Short_t Min(Short_t a, Short_t b)
Returns the smallest of a and b.
Definition TMathBase.h:199
Double_t StudentI(Double_t T, Double_t ndf)
Calculates the cumulative distribution function of Student's t-distribution second parameter stands f...
Definition TMath.cxx:2648
Double_t StudentQuantile(Double_t p, Double_t ndf, Bool_t lower_tail=kTRUE)
Computes quantiles of the Student's t-distribution 1st argument is the probability,...
Definition TMath.cxx:2676
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:124
TMarker m
Definition textangle.C:8
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2339