Logo ROOT  
Reference Guide
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#include "snprintf.h"
28
30
31
32std::map<TString,TFormula*> TLinearFitter::fgFormulaMap;
33
34
35
36/**
37
38\class TLinearFitter
39
40\ingroup MinuitOld
41
42The Linear Fitter - For fitting functions that are LINEAR IN PARAMETERS
43
44## The Linear Fitter
45
46Linear fitter is used to fit a set of data points with a linear
47combination of specified functions. Note, that "linear" in the name
48stands only for the model dependency on parameters, the specified
49functions can be nonlinear.
50The general form of this kind of model is
51~~~~
52 y(x) = a[0] + a[1]*f[1](x)+...a[n]*f[n](x)
53~~~~
54
55Functions f are fixed functions of x. For example, fitting with a
56polynomial is linear fitting in this sense.
57
58### Introduction
59
60#### The fitting method
61
62The fit is performed using the Normal Equations method with Cholesky
63decomposition.
64
65#### Why should it be used?
66
67The linear fitter is considerably faster than general non-linear
68fitters and doesn't require to set the initial values of parameters.
69
70### Using the fitter:
71
72### 1.Adding the data points:
73
74#### 1.1 To store or not to store the input data?
75 - There are 2 options in the constructor - to store or not
76 store the input data. The advantages of storing the data
77 are that you'll be able to reset the fitting model without
78 adding all the points again, and that for very large sets
79 of points the chisquare is calculated more precisely.
80 The obvious disadvantage is the amount of memory used to
81 keep all the points.
82 - Before you start adding the points, you can change the
83 store/not store option by StoreData() method.
84
85#### 1.2 The data can be added:
86 - simply point by point - AddPoint() method
87 - an array of points at once:
88 If the data is already stored in some arrays, this data
89 can be assigned to the linear fitter without physically
90 coping bytes, thanks to the Use() method of
91 TVector and TMatrix classes - AssignData() method
92
93### 2.Setting the formula
94
95#### 2.1 The linear formula syntax:
96 -Additive parts are separated by 2 plus signes "++"
97 --for example "1 ++ x" - for fitting a straight line
98 -All standard functions, undrestood by TFormula, can be used
99 as additive parts
100 --TMath functions can be used too
101 -Functions, used as additive parts, shouldn't have any parameters,
102 even if those parameters are set.
103 --for example, if normalizing a sum of a gaus(0, 1) and a
104 gaus(0, 2), don't use the built-in "gaus" of TFormula,
105 because it has parameters, take TMath::Gaus(x, 0, 1) instead.
106 -Polynomials can be used like "pol3", .."polN"
107 -If fitting a more than 3-dimensional formula, variables should
108 be numbered as follows:
109 -- x[0], x[1], x[2]... For example, to fit "1 ++ x[0] ++ x[1] ++ x[2] ++ x[3]*x[3]"
110
111#### 2.2 Setting the formula:
112
113##### 2.2.1 If fitting a 1-2-3-dimensional formula, one can create a
114 TF123 based on a linear expression and pass this function
115 to the fitter:
116 --Example:
117~~~~
118 TLinearFitter *lf = new TLinearFitter();
119 TF2 *f2 = new TF2("f2", "x ++ y ++ x*x*y*y", -2, 2, -2, 2);
120 lf->SetFormula(f2);
121~~~~
122 --The results of the fit are then stored in the function,
123 just like when the TH1::Fit or TGraph::Fit is used
124 --A linear function of this kind is by no means different
125 from any other function, it can be drawn, evaluated, etc.
126
127 --For multidimensional fitting, TFormulas of the form:
128 x[0]++...++x[n] can be used
129##### 2.2.2 There is no need to create the function if you don't want to,
130 the formula can be set by expression:
131 --Example:
132~~~~
133 // 2 is the number of dimensions
134 TLinearFitter *lf = new TLinearFitter(2);
135 lf->SetFormula("x ++ y ++ x*x*y*y");
136~~~~
137
138##### 2.2.3 The fastest functions to compute are polynomials and hyperplanes.
139 --Polynomials are set the usual way: "pol1", "pol2",...
140 --Hyperplanes are set by expression "hyp3", "hyp4", ...
141 ---The "hypN" expressions only work when the linear fitter
142 is used directly, not through TH1::Fit or TGraph::Fit.
143 To fit a graph or a histogram with a hyperplane, define
144 the function as "1++x++y".
145 ---A constant term is assumed for a hyperplane, when using
146 the "hypN" expression, so "hyp3" is in fact fitting with
147 "1++x++y++z" function.
148 --Fitting hyperplanes is much faster than fitting other
149 expressions so if performance is vital, calculate the
150 function values beforehand and give them to the fitter
151 as variables
152 --Example:
153 You want to fit "sin(x)|cos(2*x)" very fast. Calculate
154 sin(x) and cos(2*x) beforehand and store them in array *data.
155 Then:
156 TLinearFitter *lf=new TLinearFitter(2, "hyp2");
157 lf->AssignData(npoint, 2, data, y);
158
159#### 2.3 Resetting the formula
160
161##### 2.3.1 If the input data is stored (or added via AssignData() function),
162 the fitting formula can be reset without re-adding all the points.
163 --Example:
164~~~~
165 TLinearFitter *lf=new TLinearFitter("1++x++x*x");
166 lf->AssignData(n, 1, x, y, e);
167 lf->Eval()
168 //looking at the parameter significance, you see,
169 // that maybe the fit will improve, if you take out
170 // the constant term
171 lf->SetFormula("x++x*x");
172 lf->Eval();
173 ...
174~~~~
175
176##### 2.3.2 If the input data is not stored, the fitter will have to be
177 cleared and the data will have to be added again to try a
178 different formula.
179
180### 3.Accessing the fit results
181
182#### 3.1 There are methods in the fitter to access all relevant information:
183 --GetParameters, GetCovarianceMatrix, etc
184 --the t-values of parameters and their significance can be reached by
185 GetParTValue() and GetParSignificance() methods
186
187#### 3.2 If fitting with a pre-defined TF123, the fit results are also
188 written into this function.
189
190
191### 4.Robust fitting - Least Trimmed Squares regression (LTS)
192 Outliers are atypical(by definition), infrequant observations; data points
193 which do not appear to follow the characteristic distribution of the rest
194 of the data. These may reflect genuine properties of the underlying
195 phenomenon(variable), or be due to measurement errors or anomalies which
196 shouldn't be modelled. (StatSoft electronic textbook)
197
198 Even a single gross outlier can greatly influence the results of least-
199 squares fitting procedure, and in this case use of robust(resistant) methods
200 is recommended.
201
202 The method implemented here is based on the article and algorithm:
203 "Computing LTS Regression for Large Data Sets" by
204 P.J.Rousseeuw and Katrien Van Driessen
205 The idea of the method is to find the fitting coefficients for a subset
206 of h observations (out of n) with the smallest sum of squared residuals.
207 The size of the subset h should lie between (npoints + nparameters +1)/2
208 and n, and represents the minimal number of good points in the dataset.
209 The default value is set to (npoints + nparameters +1)/2, but of course
210 if you are sure that the data contains less outliers it's better to change
211 h according to your data.
212
213 To perform a robust fit, call EvalRobust() function instead of Eval() after
214 adding the points and setting the fitting function.
215 Note, that standard errors on parameters are not computed!
216
217*/
218
219
220
221////////////////////////////////////////////////////////////////////////////////
222///default c-tor, input data is stored
223///If you don't want to store the input data,
224///run the function StoreData(kFALSE) after constructor
225
228 fParams(),
229 fParCovar(),
230 fTValues(),
231 fDesign(),
232 fDesignTemp(),
233 fDesignTemp2(),
234 fDesignTemp3(),
235 fAtb(),
236 fAtbTemp(),
237 fAtbTemp2(),
238 fAtbTemp3(),
239 fFunctions(),
240 fY(),
241 fX(),
242 fE(),
243 fVal()
244{
245 fChisquare =0;
246 fNpoints =0;
247 fNdim =0;
248 fY2 =0;
249 fY2Temp =0;
250 fNfixed =0;
251 fIsSet =kFALSE;
252 fFormula =0;
253 fFixedParams=0;
254 fSpecial =0;
258 fNfunctions = 0;
259 fFormulaSize = 0;
260 fH = 0;
261}
262
263////////////////////////////////////////////////////////////////////////////////
264///The parameter stands for number of dimensions in the fitting formula
265///The input data is stored. If you don't want to store the input data,
266///run the function StoreData(kFALSE) after constructor
267
269 fVal()
270{
271 fNdim =ndim;
272 fNpoints =0;
273 fY2 =0;
274 fY2Temp =0;
275 fNfixed =0;
276 fFixedParams=0;
277 fFormula =0;
278 fIsSet =kFALSE;
279 fChisquare=0;
280 fSpecial =0;
283 fRobust = kFALSE;
284 fNfunctions = 0;
285 fFormulaSize = 0;
286 fH = 0;
287}
288
289////////////////////////////////////////////////////////////////////////////////
290///First parameter stands for number of dimensions in the fitting formula
291///Second parameter is the fitting formula: see class description for formula syntax
292///Options:
293///The option is to store or not to store the data
294///If you don't want to store the data, choose "" for the option, or run
295///StoreData(kFalse) member function after the constructor
296
297TLinearFitter::TLinearFitter(Int_t ndim, const char *formula, Option_t *opt)
298{
299 fNdim=ndim;
300 fNpoints=0;
301 fChisquare=0;
302 fY2=0;
303 fNfixed=0;
304 fFixedParams=0;
305 fSpecial=0;
307 fFormula = 0;
308 TString option=opt;
309 option.ToUpper();
310 if (option.Contains("D"))
312 else
315 SetFormula(formula);
316}
317
318////////////////////////////////////////////////////////////////////////////////
319///This constructor uses a linear function. How to create it?
320///TFormula now accepts formulas of the following kind:
321///TFormula("f", "x++y++z++x*x") or
322///TFormula("f", "x[0]++x[1]++x[2]*x[2]");
323///Other than the look, it's in no
324///way different from the regular formula, it can be evaluated,
325///drawn, etc.
326///The option is to store or not to store the data
327///If you don't want to store the data, choose "" for the option, or run
328///StoreData(kFalse) member function after the constructor
329
331{
332 fNdim=function->GetNdim();
333 if (!function->IsLinear()){
334 Int_t number=function->GetNumber();
335 if (number<299 || number>310){
336 Error("TLinearFitter", "Trying to fit with a nonlinear function");
337 return;
338 }
339 }
340 fNpoints=0;
341 fChisquare=0;
342 fY2=0;
343 fNfixed=0;
344 fFixedParams=0;
345 fSpecial=0;
346 fFormula = 0;
347 TString option=opt;
348 option.ToUpper();
349 if (option.Contains("D"))
351 else
356
358}
359
360////////////////////////////////////////////////////////////////////////////////
361/// Copy ctor
362
364 TVirtualFitter(tlf),
365 fParams(tlf.fParams),
366 fParCovar(tlf.fParCovar),
367 fTValues(tlf.fTValues),
368 fParSign(tlf.fParSign),
369 fDesign(tlf.fDesign),
370 fDesignTemp(tlf.fDesignTemp),
371 fDesignTemp2(tlf.fDesignTemp2),
372 fDesignTemp3(tlf.fDesignTemp3),
373 fAtb(tlf.fAtb),
374 fAtbTemp(tlf.fAtbTemp),
375 fAtbTemp2(tlf.fAtbTemp2),
376 fAtbTemp3(tlf.fAtbTemp3),
377 fFunctions( * (TObjArray *)tlf.fFunctions.Clone()),
378 fY(tlf.fY),
379 fY2(tlf.fY2),
380 fY2Temp(tlf.fY2Temp),
381 fX(tlf.fX),
382 fE(tlf.fE),
383 fInputFunction(tlf.fInputFunction),
384 fVal(),
385 fNpoints(tlf.fNpoints),
386 fNfunctions(tlf.fNfunctions),
387 fFormulaSize(tlf.fFormulaSize),
388 fNdim(tlf.fNdim),
389 fNfixed(tlf.fNfixed),
390 fSpecial(tlf.fSpecial),
391 fFormula(0),
392 fIsSet(tlf.fIsSet),
393 fStoreData(tlf.fStoreData),
394 fChisquare(tlf.fChisquare),
395 fH(tlf.fH),
396 fRobust(tlf.fRobust),
397 fFitsample(tlf.fFitsample),
398 fFixedParams(0)
399{
400 // make a deep copy of managed objects
401 // fFormula, fFixedParams and fFunctions
402
403 if ( tlf.fFixedParams && fNfixed > 0 ) {
405 for(Int_t i=0; i<fNfixed; ++i)
406 fFixedParams[i]=tlf.fFixedParams[i];
407 }
408 if (tlf.fFormula) {
409 fFormula = new char[fFormulaSize+1];
410 strlcpy(fFormula,tlf.fFormula,fFormulaSize+1);
411 }
412
413}
414
415
416////////////////////////////////////////////////////////////////////////////////
417/// Linear fitter cleanup.
418
420{
421 if (fFormula) {
422 delete [] fFormula;
423 fFormula = 0;
424 }
425 if (fFixedParams) {
426 delete [] fFixedParams;
427 fFixedParams = 0;
428 }
429 fInputFunction = 0;
430
431 //fFunctions.Delete();
433
434}
435
436////////////////////////////////////////////////////////////////////////////////
437/// Assignment operator
438
440{
441 if(this!=&tlf) {
442
452
453 fAtb.ResizeTo(tlf.fAtb); fAtb=tlf.fAtb;
457
458 // use clear instead of delete
461
462 fY.ResizeTo(tlf.fY); fY = tlf.fY;
463 fX.ResizeTo(tlf.fX); fX = tlf.fX;
464 fE.ResizeTo(tlf.fE); fE = tlf.fE;
465
466 fY2 = tlf.fY2;
467 fY2Temp = tlf.fY2Temp;
468 for(Int_t i = 0; i < 1000; i++) fVal[i] = tlf.fVal[i];
469
470 if(fInputFunction) { delete fInputFunction; fInputFunction = 0; }
472
473 fNpoints=tlf.fNpoints;
476 fNdim=tlf.fNdim;
477 fNfixed=tlf.fNfixed;
478 fSpecial=tlf.fSpecial;
479
480 if(fFormula) { delete [] fFormula; fFormula = 0; }
481 if (tlf.fFormula) {
482 fFormula = new char[fFormulaSize+1];
483 strlcpy(fFormula,tlf.fFormula,fFormulaSize+1);
484 }
485
486 fIsSet=tlf.fIsSet;
489
490 fH=tlf.fH;
491 fRobust=tlf.fRobust;
493
494 if(fFixedParams) { delete [] fFixedParams; fFixedParams = 0; }
495 if ( tlf.fFixedParams && fNfixed > 0 ) {
497 for(Int_t i=0; i< fNfixed; ++i)
498 fFixedParams[i]=tlf.fFixedParams[i];
499 }
500
501 }
502 return *this;
503}
504
505////////////////////////////////////////////////////////////////////////////////
506///Add another linear fitter to this linear fitter. Points and Design matrices
507///are added, but the previos fitting results (if any) are deleted.
508///Fitters must have same formulas (this is not checked). Fixed parameters are not changed
509
511{
512 fParams.Zero();
513 fParCovar.Zero();
514 fTValues.Zero();
515 fParSign.Zero();
516
517 fDesign += tlf->fDesign;
518
519 fDesignTemp += tlf->fDesignTemp;
522 fAtb += tlf->fAtb;
523 fAtbTemp += tlf->fAtbTemp;
524 fAtbTemp2 += tlf->fAtbTemp2;
525 fAtbTemp3 += tlf->fAtbTemp3;
526
527 if (fStoreData){
529 Int_t newsize = fNpoints+tlf->fNpoints;
530 if (size<newsize){
531 fY.ResizeTo(newsize);
532 fE.ResizeTo(newsize);
533 fX.ResizeTo(newsize, fNdim);
534 }
535 for (Int_t i=fNpoints; i<newsize; i++){
536 fY(i)=tlf->fY(i-fNpoints);
537 fE(i)=tlf->fE(i-fNpoints);
538 for (Int_t j=0; j<fNdim; j++){
539 fX(i,j)=tlf->fX(i-fNpoints, j);
540 }
541 }
542 }
543 fY2 += tlf->fY2;
544 fY2Temp += tlf->fY2Temp;
545 fNpoints += tlf->fNpoints;
546 //fInputFunction=(TFormula*)tlf.fInputFunction->Clone();
547
548 fChisquare=0;
549 fH=0;
550 fRobust=0;
551}
552
553////////////////////////////////////////////////////////////////////////////////
554///Adds 1 point to the fitter.
555///First parameter stands for the coordinates of the point, where the function is measured
556///Second parameter - the value being fitted
557///Third parameter - weight(measurement error) of this point (=1 by default)
558
560{
561 Int_t size;
562 fNpoints++;
563 if (fStoreData){
565 if (size<fNpoints){
569 }
570
571 Int_t j=fNpoints-1;
572 fY(j)=y;
573 fE(j)=e;
574 for (Int_t i=0; i<fNdim; i++)
575 fX(j,i)=x[i];
576 }
577 //add the point to the design matrix, if the formula has been set
578 // (LM: why condition !fRobust ?? - remove it to fix Coverity 11602)
579 // when 100< fSpecial < 200 (Polynomial) fFunctions is not empty
580 // (see implementation of SetFormula method)
581 if (fFunctions.IsEmpty()&&(!fInputFunction)&&(fSpecial<=200)){
582 Error("AddPoint", "Point can't be added, because the formula hasn't been set");
583 return;
584 }
585 if (!fRobust) AddToDesign(x, y, e);
586}
587
588////////////////////////////////////////////////////////////////////////////////
589///This function is to use when you already have all the data in arrays
590///and don't want to copy them into the fitter. In this function, the Use() method
591///of TVectorD and TMatrixD is used, so no bytes are physically moved around.
592///First parameter - number of points to fit
593///Second parameter - number of variables in the model
594///Third parameter - the variables of the model, stored in the following way:
595///(x0(0), x1(0), x2(0), x3(0), x0(1), x1(1), x2(1), x3(1),...
596
598{
599 if (npoints<fNpoints){
600 Error("AddData", "Those points are already added");
601 return;
602 }
603 Bool_t same=kFALSE;
604 if (fX.GetMatrixArray()==x && fY.GetMatrixArray()==y){
605 if (e){
606 if (fE.GetMatrixArray()==e)
607 same=kTRUE;
608 }
609 }
610
611 fX.Use(npoints, xncols, x);
612 fY.Use(npoints, y);
613 if (e)
614 fE.Use(npoints, e);
615 else {
616 fE.ResizeTo(npoints);
617 fE=1;
618 }
619 Int_t xfirst;
620 if (!fFunctions.IsEmpty() || fInputFunction || fSpecial>200) {
621 if (same)
622 xfirst=fNpoints;
623
624 else
625 xfirst=0;
626 for (Int_t i=xfirst; i<npoints; i++)
627 AddToDesign(TMatrixDRow(fX, i).GetPtr(), fY(i), fE(i));
628 }
629 fNpoints=npoints;
630}
631
632////////////////////////////////////////////////////////////////////////////////
633///Add a point to the AtA matrix and to the Atb vector.
634
636{
637
638
639 Int_t i, j, ii;
640 y/=e;
641
642// Double_t fVal[1000];
643
644 if ((fSpecial>100)&&(fSpecial<200)){
645 //polynomial fitting
646 Int_t npar=fSpecial-100;
647 fVal[0]=1;
648 for (i=1; i<npar; i++)
649 fVal[i]=fVal[i-1]*x[0];
650 for (i=0; i<npar; i++)
651 fVal[i]/=e;
652 } else if (fSpecial>200){
653 //Hyperplane fitting. Constant term is added
654 Int_t npar=fSpecial-201;
655 fVal[0]=1./e;
656 for (i=0; i<npar; i++)
657 fVal[i+1]=x[i]/e;
658 } else {
659 //general case
660 for (ii=0; ii<fNfunctions; ii++){
661 if (!fFunctions.IsEmpty()){
662 // ffunctions can be TF1 or TFormula depending on how they are created
663 TObject * obj = fFunctions.UncheckedAt(ii);
664 if (obj->IsA() == TFormula::Class() ) {
665 TFormula *f1 = (TFormula*)(obj);
666 fVal[ii]=f1->EvalPar(x)/e;
667 }
668 else if (obj->IsA() == TF1::Class() ) {
669 TF1 *f1 = (TF1*)(obj);
670 fVal[ii]=f1->EvalPar(x)/e;
671 }
672 else {
673 Error("AddToDesign","Basis Function %s is of an invalid type %s",obj->GetName(),obj->IsA()->GetName());
674 return;
675 }
676 } else {
678 if (!f){
679 Error("AddToDesign","Function %s has no linear parts - maybe missing a ++ in the formula expression",fInputFunction->GetName());
680 return;
681 }
682 fVal[ii]=f->EvalPar(x)/e;
683 }
684 }
685 }
686 //additional matrices for numerical stability
687 for (i=0; i<fNfunctions; i++){
688 for (j=0; j<i; j++)
689 fDesignTemp3(j, i)+=fVal[i]*fVal[j];
690 fDesignTemp3(i, i)+=fVal[i]*fVal[i];
691 fAtbTemp3(i)+=fVal[i]*y;
692
693 }
694 fY2Temp+=y*y;
696
697 if (fNpoints % 100 == 0 && fNpoints>100){
701 fAtbTemp3.Zero();
702 if (fNpoints % 10000 == 0 && fNpoints>10000){
706 fAtbTemp2.Zero();
707 fY2+=fY2Temp;
708 fY2Temp=0;
709 if (fNpoints % 1000000 == 0 && fNpoints>1000000){
712 fAtb+=fAtbTemp;
713 fAtbTemp.Zero();
714 }
715 }
716 }
717}
718
719////////////////////////////////////////////////////////////////////////////////
720
722{
723 if (fDesignTemp3.GetNrows()>0){
732 fAtb+=fAtbTemp;
733 fAtbTemp3.Zero();
734 fAtbTemp2.Zero();
735 fAtbTemp.Zero();
736
737 fY2+=fY2Temp;
738 fY2Temp=0;
739 }
740}
741
742////////////////////////////////////////////////////////////////////////////////
743///Clears everything. Used in TH1::Fit and TGraph::Fit().
744
746{
747 fParams.Clear();
749 fTValues.Clear();
750 fParSign.Clear();
751 fDesign.Clear();
755 fAtb.Clear();
756 fAtbTemp.Clear();
761 fY.Clear();
762 fX.Clear();
763 fE.Clear();
764
765 fNpoints=0;
766 fNfunctions=0;
767 fFormulaSize=0;
768 fNdim=0;
769 if (fFormula) delete [] fFormula;
770 fFormula=0;
771 fIsSet=0;
772 if (fFixedParams) delete [] fFixedParams;
773 fFixedParams=0;
774
775 fChisquare=0;
776 fY2=0;
777 fSpecial=0;
780}
781
782////////////////////////////////////////////////////////////////////////////////
783///To be used when different sets of points are fitted with the same formula.
784
786{
787 fDesign.Zero();
788 fAtb.Zero();
792 fAtbTemp.Zero();
793 fAtbTemp2.Zero();
794 fAtbTemp3.Zero();
795
796 fParams.Zero();
797 fParCovar.Zero();
798 fTValues.Zero();
799 fParSign.Zero();
800
801 for (Int_t i=0; i<fNfunctions; i++)
802 fFixedParams[i]=0;
803 fChisquare=0;
804 fNpoints=0;
805
806}
807
808////////////////////////////////////////////////////////////////////////////////
809///Calculates the chisquare.
810
812{
813 Int_t i, j;
814 Double_t sumtotal2;
815 Double_t temp, temp2;
816
817 if (!fStoreData){
818 sumtotal2 = 0;
819 for (i=0; i<fNfunctions; i++){
820 for (j=0; j<i; j++){
821 sumtotal2 += 2*fParams(i)*fParams(j)*fDesign(j, i);
822 }
823 sumtotal2 += fParams(i)*fParams(i)*fDesign(i, i);
824 sumtotal2 -= 2*fParams(i)*fAtb(i);
825 }
826 sumtotal2 += fY2;
827 } else {
828 sumtotal2 = 0;
829 if (fInputFunction){
830 for (i=0; i<fNpoints; i++){
831 temp = fInputFunction->EvalPar(TMatrixDRow(fX, i).GetPtr());
832 temp2 = (fY(i)-temp)*(fY(i)-temp);
833 temp2 /= fE(i)*fE(i);
834 sumtotal2 += temp2;
835 }
836 } else {
837 sumtotal2 = 0;
838 Double_t val[100];
839 for (Int_t point=0; point<fNpoints; point++){
840 temp = 0;
841 if ((fSpecial>100)&&(fSpecial<200)){
842 Int_t npar = fSpecial-100;
843 val[0] = 1;
844 for (i=1; i<npar; i++)
845 val[i] = val[i-1]*fX(point, 0);
846 for (i=0; i<npar; i++)
847 temp += fParams(i)*val[i];
848 } else {
849 if (fSpecial>200) {
850 //hyperplane case
851 Int_t npar = fSpecial-201;
852 temp+=fParams(0);
853 for (i=0; i<npar; i++)
854 temp += fParams(i+1)*fX(point, i);
855 } else {
856 for (j=0; j<fNfunctions; j++) {
858 val[j] = f1->EvalPar(TMatrixDRow(fX, point).GetPtr());
859 temp += fParams(j)*val[j];
860 }
861 }
862 }
863 temp2 = (fY(point)-temp)*(fY(point)-temp);
864 temp2 /= fE(point)*fE(point);
865 sumtotal2 += temp2;
866 }
867 }
868 }
869 fChisquare = sumtotal2;
870
871}
872
873////////////////////////////////////////////////////////////////////////////////
874/// Computes parameters' t-values and significance
875
877{
878 for (Int_t i=0; i<fNfunctions; i++){
879 fTValues(i) = fParams(i)/(TMath::Sqrt(fParCovar(i, i)));
881 }
882}
883
884////////////////////////////////////////////////////////////////////////////////
885/// Perform the fit and evaluate the parameters
886/// Returns 0 if the fit is ok, 1 if there are errors
887
889{
890 Double_t e;
891 if (fFunctions.IsEmpty()&&(!fInputFunction)&&(fSpecial<=200)){
892 Error("TLinearFitter::Eval", "The formula hasn't been set");
893 return 1;
894 }
895 //
900
901 fChisquare=0;
902
903 if (!fIsSet){
905 if (!update){
906 //no points to fit
907 fParams.Zero();
908 fParCovar.Zero();
909 fTValues.Zero();
910 fParSign.Zero();
911 fChisquare=0;
912 if (fInputFunction){
914 for (Int_t i=0; i<fNfunctions; i++)
915 ((TF1*)fInputFunction)->SetParError(i, 0);
916 ((TF1*)fInputFunction)->SetChisquare(0);
917 ((TF1*)fInputFunction)->SetNDF(0);
918 ((TF1*)fInputFunction)->SetNumberFitPoints(0);
919 }
920 return 1;
921 }
922 }
923
925
926//fixing fixed parameters, if there are any
927 Int_t i, ii, j=0;
928 if (fNfixed>0){
929 for (ii=0; ii<fNfunctions; ii++)
930 fDesignTemp(ii, fNfixed) = fAtb(ii);
931 for (i=0; i<fNfunctions; i++){
932 if (fFixedParams[i]){
933 for (ii=0; ii<i; ii++)
934 fDesignTemp(ii, j) = fDesign(ii, i);
935 for (ii=i; ii<fNfunctions; ii++)
936 fDesignTemp(ii, j) = fDesign(i, ii);
937 j++;
938 for (ii=0; ii<fNfunctions; ii++){
939 fAtb(ii)-=fParams(i)*(fDesignTemp(ii, j-1));
940 }
941 }
942 }
943 for (i=0; i<fNfunctions; i++){
944 if (fFixedParams[i]){
945 for (ii=0; ii<fNfunctions; ii++){
946 fDesign(ii, i) = 0;
947 fDesign(i, ii) = 0;
948 }
949 fDesign (i, i) = 1;
950 fAtb(i) = fParams(i);
951 }
952 }
953 }
954
955 TDecompChol chol(fDesign);
956 Bool_t ok;
957 TVectorD coef(fNfunctions);
958 coef=chol.Solve(fAtb, ok);
959 if (!ok){
960 Error("Eval","Matrix inversion failed");
961 fParams.Zero();
962 fParCovar.Zero();
963 fTValues.Zero();
964 fParSign.Zero();
965 if (fInputFunction){
968 ((TF1*)fInputFunction)->SetChisquare(0);
970 ((TF1*)fInputFunction)->SetNumberFitPoints(fNpoints);
971 }
972 }
973 return 1;
974 }
975 fParams=coef;
976 fParCovar=chol.Invert();
977
978 if (fInputFunction){
981 for (i=0; i<fNfunctions; i++){
982 e = TMath::Sqrt(fParCovar(i, i));
983 ((TF1*)fInputFunction)->SetParError(i, e);
984 }
985 if (!fObjectFit)
986 ((TF1*)fInputFunction)->SetChisquare(GetChisquare());
988 ((TF1*)fInputFunction)->SetNumberFitPoints(fNpoints);
989 }
990 }
991
992 //if parameters were fixed, change the design matrix back as it was before fixing
993 j = 0;
994 if (fNfixed>0){
995 for (i=0; i<fNfunctions; i++){
996 if (fFixedParams[i]){
997 for (ii=0; ii<i; ii++){
998 fDesign(ii, i) = fDesignTemp(ii, j);
999 fAtb(ii) = fDesignTemp(ii, fNfixed);
1000 }
1001 for (ii=i; ii<fNfunctions; ii++){
1002 fDesign(i, ii) = fDesignTemp(ii, j);
1003 fAtb(ii) = fDesignTemp(ii, fNfixed);
1004 }
1005 j++;
1006 }
1007 }
1008 }
1009 return 0;
1010}
1011
1012////////////////////////////////////////////////////////////////////////////////
1013///Fixes paramter `#ipar` at its current value.
1014
1016{
1017 if (fParams.NonZeros()<1){
1018 Error("FixParameter", "no value available to fix the parameter");
1019 return;
1020 }
1021 if (ipar>fNfunctions || ipar<0){
1022 Error("FixParameter", "illegal parameter value");
1023 return;
1024 }
1025 if (fNfixed==fNfunctions) {
1026 Error("FixParameter", "no free parameters left");
1027 return;
1028 }
1029 if (!fFixedParams)
1031 fFixedParams[ipar] = 1;
1032 fNfixed++;
1033}
1034
1035////////////////////////////////////////////////////////////////////////////////
1036///Fixes parameter `#ipar` at value `parvalue`.
1037
1039{
1040 if (ipar>fNfunctions || ipar<0){
1041 Error("FixParameter", "illegal parameter value");
1042 return;
1043 }
1044 if (fNfixed==fNfunctions) {
1045 Error("FixParameter", "no free parameters left");
1046 return;
1047 }
1048 if(!fFixedParams)
1050 fFixedParams[ipar] = 1;
1053 fParams(ipar) = parvalue;
1054 fNfixed++;
1055}
1056
1057////////////////////////////////////////////////////////////////////////////////
1058///Releases parameter `#ipar`.
1059
1061{
1062 if (ipar>fNfunctions || ipar<0){
1063 Error("ReleaseParameter", "illegal parameter value");
1064 return;
1065 }
1066 if (!fFixedParams[ipar]){
1067 Warning("ReleaseParameter","This parameter is not fixed\n");
1068 return;
1069 } else {
1070 fFixedParams[ipar] = 0;
1071 fNfixed--;
1072 }
1073}
1074
1075////////////////////////////////////////////////////////////////////////////////
1076///Get the Atb vector - a vector, used for internal computations
1077
1079{
1080 if (v.GetNoElements()!=fAtb.GetNoElements())
1081 v.ResizeTo(fAtb.GetNoElements());
1082 v = fAtb;
1083}
1084
1085////////////////////////////////////////////////////////////////////////////////
1086/// Get the Chisquare.
1087
1089{
1090 if (fChisquare > 1e-16)
1091 return fChisquare;
1092 else {
1093 Chisquare();
1094 return fChisquare;
1095 }
1096}
1097
1098////////////////////////////////////////////////////////////////////////////////
1099///Computes point-by-point confidence intervals for the fitted function
1100///Parameters:
1101///n - number of points
1102///ndim - dimensions of points
1103///x - points, at which to compute the intervals, for ndim > 1
1104/// should be in order: (x0,y0, x1, y1, ... xn, yn)
1105///ci - computed intervals are returned in this array
1106///cl - confidence level, default=0.95
1107///
1108///NOTE, that this method can only be used when the fitting function inherits from a TF1,
1109///so it's not possible when the fitting function was set as a string or as a pure TFormula
1110
1112{
1113 if (fInputFunction){
1114 Double_t *grad = new Double_t[fNfunctions];
1115 Double_t *sum_vector = new Double_t[fNfunctions];
1116 Double_t c=0;
1118 Double_t t = TMath::StudentQuantile(0.5 + cl/2, df);
1119 Double_t chidf = TMath::Sqrt(fChisquare/df);
1120
1121 for (Int_t ipoint=0; ipoint<n; ipoint++){
1122 c=0;
1123 ((TF1*)(fInputFunction))->GradientPar(x+ndim*ipoint, grad);
1124 //multiply the covariance matrix by gradient
1125 for (Int_t irow=0; irow<fNfunctions; irow++){
1126 sum_vector[irow]=0;
1127 for (Int_t icol=0; icol<fNfunctions; icol++){
1128 sum_vector[irow]+=fParCovar(irow,icol)*grad[icol];
1129 }
1130 }
1131 for (Int_t i=0; i<fNfunctions; i++)
1132 c+=grad[i]*sum_vector[i];
1133 c=TMath::Sqrt(c);
1134 ci[ipoint]=c*t*chidf;
1135 }
1136
1137 delete [] grad;
1138 delete [] sum_vector;
1139 }
1140}
1141
1142////////////////////////////////////////////////////////////////////////////////
1143///Computes confidence intervals at level cl. Default is 0.95
1144///The TObject parameter can be a TGraphErrors, a TGraph2DErrors or a TH123.
1145///For Graphs, confidence intervals are computed for each point,
1146///the value of the graph at that point is set to the function value at that
1147///point, and the graph y-errors (or z-errors) are set to the value of
1148///the confidence interval at that point
1149///For Histograms, confidence intervals are computed for each bin center
1150///The bin content of this bin is then set to the function value at the bin
1151///center, and the bin error is set to the confidence interval value.
1152///Allowed combinations:
1153///Fitted object Passed object
1154///TGraph TGraphErrors, TH1
1155///TGraphErrors, AsymmErrors TGraphErrors, TH1
1156///TH1 TGraphErrors, TH1
1157///TGraph2D TGraph2DErrors, TH2
1158///TGraph2DErrors TGraph2DErrors, TH2
1159///TH2 TGraph2DErrors, TH2
1160///TH3 TH3
1161
1163{
1164 if (!fInputFunction) {
1165 Error("GetConfidenceIntervals", "The case of fitting not with a TFormula is not yet implemented");
1166 return;
1167 }
1168
1169 //TGraph//////////////////
1170
1171 if (obj->InheritsFrom(TGraph::Class())) {
1172 TGraph *gr = (TGraph*)obj;
1173 if (!gr->GetEY()){
1174 Error("GetConfidenceIntervals", "A TGraphErrors should be passed instead of a graph");
1175 return;
1176 }
1178 Error("GetConfidenceIntervals", "A TGraph2DErrors should be passed instead of a graph");
1179 return;
1180 }
1182 if (((TH1*)(fObjectFit))->GetDimension()>1){
1183 Error("GetConfidenceIntervals", "A TGraph2DErrors or a TH23 should be passed instead of a graph");
1184 return;
1185 }
1186 }
1187
1188 GetConfidenceIntervals(gr->GetN(),1,gr->GetX(), gr->GetEY(), cl);
1189 for (Int_t i=0; i<gr->GetN(); i++)
1190 gr->SetPoint(i, gr->GetX()[i], fInputFunction->Eval(gr->GetX()[i]));
1191 }
1192
1193 //TGraph2D///////////////
1194 else if (obj->InheritsFrom(TGraph2D::Class())) {
1195 TGraph2D *gr2 = (TGraph2D*)obj;
1196 if (!gr2->GetEZ()){
1197 Error("GetConfidenceIntervals", "A TGraph2DErrors should be passed instead of a TGraph2D");
1198 return;
1199 }
1201 Error("GetConfidenceIntervals", "A TGraphErrors should be passed instead of a TGraph2D");
1202 return;
1203 }
1205 if (((TH1*)(fObjectFit))->GetDimension()==1){
1206 Error("GetConfidenceIntervals", "A TGraphErrors or a TH1 should be passed instead of a graph");
1207 return;
1208 }
1209 }
1210 Double_t xy[2];
1211 Int_t np = gr2->GetN();
1212 Double_t *grad = new Double_t[fNfunctions];
1213 Double_t *sum_vector = new Double_t[fNfunctions];
1214 Double_t *x = gr2->GetX();
1215 Double_t *y = gr2->GetY();
1216 Double_t t = TMath::StudentQuantile(0.5 + cl/2, ((TF1*)(fInputFunction))->GetNDF());
1217 Double_t chidf = TMath::Sqrt(fChisquare/((TF1*)(fInputFunction))->GetNDF());
1218 Double_t c = 0;
1219 for (Int_t ipoint=0; ipoint<np; ipoint++){
1220 c=0;
1221 xy[0]=x[ipoint];
1222 xy[1]=y[ipoint];
1223 ((TF1*)(fInputFunction))->GradientPar(xy, grad);
1224 for (Int_t irow=0; irow<fNfunctions; irow++){
1225 sum_vector[irow]=0;
1226 for (Int_t icol=0; icol<fNfunctions; icol++)
1227 sum_vector[irow]+=fParCovar(irow, icol)*grad[icol];
1228 }
1229 for (Int_t i=0; i<fNfunctions; i++)
1230 c+=grad[i]*sum_vector[i];
1231 c=TMath::Sqrt(c);
1232 gr2->SetPoint(ipoint, xy[0], xy[1], fInputFunction->EvalPar(xy));
1233 gr2->GetEZ()[ipoint]=c*t*chidf;
1234 }
1235 delete [] grad;
1236 delete [] sum_vector;
1237 }
1238
1239 //TH1////////////////////////
1240 else if (obj->InheritsFrom(TH1::Class())) {
1242 if (((TH1*)obj)->GetDimension()>1){
1243 Error("GetConfidenceIntervals", "Fitted graph and passed histogram have different number of dimensions");
1244 return;
1245 }
1246 }
1248 if (((TH1*)obj)->GetDimension()!=2){
1249 Error("GetConfidenceIntervals", "Fitted graph and passed histogram have different number of dimensions");
1250 return;
1251 }
1252 }
1254 if (((TH1*)(fObjectFit))->GetDimension()!=((TH1*)(obj))->GetDimension()){
1255 Error("GetConfidenceIntervals", "Fitted and passed histograms have different number of dimensions");
1256 return;
1257 }
1258 }
1259
1260
1261 TH1 *hfit = (TH1*)obj;
1262 Double_t *grad = new Double_t[fNfunctions];
1263 Double_t *sum_vector = new Double_t[fNfunctions];
1264 Double_t x[3];
1265 Int_t hxfirst = hfit->GetXaxis()->GetFirst();
1266 Int_t hxlast = hfit->GetXaxis()->GetLast();
1267 Int_t hyfirst = hfit->GetYaxis()->GetFirst();
1268 Int_t hylast = hfit->GetYaxis()->GetLast();
1269 Int_t hzfirst = hfit->GetZaxis()->GetFirst();
1270 Int_t hzlast = hfit->GetZaxis()->GetLast();
1271
1272 TAxis *xaxis = hfit->GetXaxis();
1273 TAxis *yaxis = hfit->GetYaxis();
1274 TAxis *zaxis = hfit->GetZaxis();
1275 Double_t t = TMath::StudentQuantile(0.5 + cl/2, ((TF1*)(fInputFunction))->GetNDF());
1276 Double_t chidf = TMath::Sqrt(fChisquare/((TF1*)(fInputFunction))->GetNDF());
1277 Double_t c=0;
1278 for (Int_t binz=hzfirst; binz<=hzlast; binz++){
1279 x[2]=zaxis->GetBinCenter(binz);
1280 for (Int_t biny=hyfirst; biny<=hylast; biny++) {
1281 x[1]=yaxis->GetBinCenter(biny);
1282 for (Int_t binx=hxfirst; binx<=hxlast; binx++) {
1283 x[0]=xaxis->GetBinCenter(binx);
1284 ((TF1*)(fInputFunction))->GradientPar(x, grad);
1285 c=0;
1286 for (Int_t irow=0; irow<fNfunctions; irow++){
1287 sum_vector[irow]=0;
1288 for (Int_t icol=0; icol<fNfunctions; icol++)
1289 sum_vector[irow]+=fParCovar(irow, icol)*grad[icol];
1290 }
1291 for (Int_t i=0; i<fNfunctions; i++)
1292 c+=grad[i]*sum_vector[i];
1293 c=TMath::Sqrt(c);
1294 hfit->SetBinContent(binx, biny, binz, fInputFunction->EvalPar(x));
1295 hfit->SetBinError(binx, biny, binz, c*t*chidf);
1296 }
1297 }
1298 }
1299 delete [] grad;
1300 delete [] sum_vector;
1301 }
1302 else {
1303 Error("GetConfidenceIntervals", "This object type is not supported");
1304 return;
1305 }
1306}
1307
1308////////////////////////////////////////////////////////////////////////////////
1309///Returns covariance matrix
1310
1312{
1313 Double_t *p = const_cast<Double_t*>(fParCovar.GetMatrixArray());
1314 return p;
1315}
1316
1317////////////////////////////////////////////////////////////////////////////////
1318///Returns covariance matrix
1319
1321{
1322 if (matr.GetNrows()!=fNfunctions || matr.GetNcols()!=fNfunctions){
1324 }
1325 matr = fParCovar;
1326}
1327
1328////////////////////////////////////////////////////////////////////////////////
1329///Returns the internal design matrix
1330
1332{
1333 if (matr.GetNrows()!=fNfunctions || matr.GetNcols()!=fNfunctions){
1335 }
1336 matr = fDesign;
1337}
1338
1339////////////////////////////////////////////////////////////////////////////////
1340///Returns parameter errors
1341
1343{
1344 if (vpar.GetNoElements()!=fNfunctions) {
1345 vpar.ResizeTo(fNfunctions);
1346 }
1347 for (Int_t i=0; i<fNfunctions; i++)
1348 vpar(i) = TMath::Sqrt(fParCovar(i, i));
1349
1350}
1351
1352////////////////////////////////////////////////////////////////////////////////
1353///Returns parameter values
1354
1356{
1357 if (vpar.GetNoElements()!=fNfunctions) {
1358 vpar.ResizeTo(fNfunctions);
1359 }
1360 vpar=fParams;
1361}
1362
1363////////////////////////////////////////////////////////////////////////////////
1364///Returns the value and the name of the parameter `#ipar`
1365///NB: In the calling function the argument name must be set large enough
1366
1367Int_t TLinearFitter::GetParameter(Int_t ipar,char* name,Double_t& value,Double_t& /*verr*/,Double_t& /*vlow*/, Double_t& /*vhigh*/) const
1368{
1369 if (ipar<0 || ipar>fNfunctions) {
1370 Error("GetParError", "illegal value of parameter");
1371 return 0;
1372 }
1373 value = fParams(ipar);
1374 if (fInputFunction)
1375 strcpy(name, fInputFunction->GetParName(ipar));
1376 else
1377 name = 0;
1378 return 1;
1379}
1380
1381
1382////////////////////////////////////////////////////////////////////////////////
1383///Returns the error of parameter `#ipar`
1384
1386{
1387 if (ipar<0 || ipar>fNfunctions) {
1388 Error("GetParError", "illegal value of parameter");
1389 return 0;
1390 }
1391
1392 return TMath::Sqrt(fParCovar(ipar, ipar));
1393}
1394
1395
1396////////////////////////////////////////////////////////////////////////////////
1397///Returns name of parameter `#ipar`
1398
1399const char *TLinearFitter::GetParName(Int_t ipar) const
1400{
1401 if (ipar<0 || ipar>fNfunctions) {
1402 Error("GetParError", "illegal value of parameter");
1403 return 0;
1404 }
1405 if (fInputFunction)
1406 return fInputFunction->GetParName(ipar);
1407 return "";
1408}
1409
1410////////////////////////////////////////////////////////////////////////////////
1411///Returns the t-value for parameter `#ipar`
1412
1414{
1415 if (ipar<0 || ipar>fNfunctions) {
1416 Error("GetParTValue", "illegal value of parameter");
1417 return 0;
1418 }
1419 if (!fTValues.NonZeros())
1421 return fTValues(ipar);
1422}
1423
1424////////////////////////////////////////////////////////////////////////////////
1425///Returns the significance of parameter `#ipar`
1426
1428{
1429 if (ipar<0 || ipar>fNfunctions) {
1430 Error("GetParSignificance", "illegal value of parameter");
1431 return 0;
1432 }
1433 if (!fParSign.NonZeros())
1435 return fParSign(ipar);
1436}
1437
1438////////////////////////////////////////////////////////////////////////////////
1439///For robust lts fitting, returns the sample, on which the best fit was based
1440
1442{
1443 if (!fRobust){
1444 Error("GetFitSample", "there is no fit sample in ordinary least-squares fit");
1445 return;
1446 }
1447 for (Int_t i=0; i<fNpoints; i++)
1449
1450}
1451
1452////////////////////////////////////////////////////////////////////////////////
1453///Merge objects in list
1454
1456{
1457 if (!list) return -1;
1458 TIter next(list);
1459 TLinearFitter *lfit = 0;
1460 while ((lfit = (TLinearFitter*)next())) {
1461 if (!lfit->InheritsFrom(TLinearFitter::Class())) {
1462 Error("Add","Attempt to add object of class: %s to a %s",lfit->ClassName(),this->ClassName());
1463 return -1;
1464 }
1465 Add(lfit);
1466 }
1467 return 0;
1468}
1469////////////////////////////////////////////////////////////////////////////////
1470///set the basis functions in case the fitting function is not
1471/// set directly
1472/// The TLinearFitter will manage and delete the functions contained in the list
1473
1475{
1476 fFunctions = *(functions);
1478 int size = fFunctions.GetEntries();
1479
1481 //change the size of design matrix
1490 if (fFixedParams)
1491 delete [] fFixedParams;
1493 fDesign.Zero();
1494 fAtb.Zero();
1495 fDesignTemp.Zero();
1498 fAtbTemp.Zero();
1499 fAtbTemp2.Zero();
1500 fAtbTemp3.Zero();
1501 fY2Temp=0;
1502 fY2=0;
1503 for (int i=0; i<size; i++)
1504 fFixedParams[i]=0;
1505 fIsSet=kFALSE;
1506 fChisquare=0;
1507
1508}
1509
1510
1511////////////////////////////////////////////////////////////////////////////////
1512///set the number of dimensions
1513
1515{
1516 fNdim=ndim;
1517 fY.ResizeTo(ndim+1);
1518 fX.ResizeTo(ndim+1, ndim);
1519 fE.ResizeTo(ndim+1);
1520
1521 fNpoints=0;
1522 fIsSet=kFALSE;
1523}
1524
1525////////////////////////////////////////////////////////////////////////////////
1526///Additive parts should be separated by "++".
1527///Examples (ai are parameters to fit):
1528///1.fitting function: a0*x0 + a1*x1 + a2*x2
1529/// input formula "x[0]++x[1]++x[2]"
1530///2.TMath functions can be used:
1531/// fitting function: a0*TMath::Gaus(x, 0, 1) + a1*y
1532/// input formula: "TMath::Gaus(x, 0, 1)++y"
1533///fills the array of functions
1534
1535void TLinearFitter::SetFormula(const char *formula)
1536{
1537 Int_t size = 0, special = 0;
1538 Int_t i;
1539 //Int_t len = strlen(formula);
1540 if (fInputFunction)
1541 fInputFunction = nullptr;
1542 if (fFormula)
1543 delete [] fFormula;
1544 fFormulaSize = strlen(formula);
1545 fFormula = new char[fFormulaSize+1];
1546 strlcpy(fFormula, formula,fFormulaSize+1);
1547 fSpecial = 0;
1548 //in case of a hyperplane:
1549 char *fstring;
1550 fstring = (char *)strstr(fFormula, "hyp");
1551 if (fstring!=0){
1552 // isHyper = kTRUE;
1553 fstring+=3;
1554 sscanf(fstring, "%d", &size);
1555 //+1 for the constant term
1556 size++;
1557 fSpecial=200+size;
1558 }
1559
1560 if (fSpecial==0) {
1561 //in case it's not a hyperplane
1562 TString sstring(fFormula);
1563 sstring = sstring.ReplaceAll("++", 2, "|", 1);
1564 TString replaceformula;
1565
1566 //count the number of functions
1567
1568 TObjArray *oa = sstring.Tokenize("|");
1569
1570 //change the size of functions array and clear it
1571 if (!fFunctions.IsEmpty())
1572 fFunctions.Clear();
1573
1574 // do not own the functions in this case
1576
1579
1580 //check if the old notation of xi is used somewhere instead of x[i]
1581 char pattern[12];
1582 char replacement[14];
1583 for (i=0; i<fNdim; i++){
1584 snprintf(pattern,sizeof(pattern), "x%d", i);
1585 snprintf(replacement,sizeof(replacement), "x[%d]", i);
1586 sstring = sstring.ReplaceAll(pattern, Int_t(i/10)+2, replacement, Int_t(i/10)+4);
1587 }
1588
1589 //fill the array of functions
1590 oa = sstring.Tokenize("|");
1591 for (i=0; i<fNfunctions; i++) {
1592 replaceformula = ((TObjString *)oa->UncheckedAt(i))->GetString();
1593 // look first if exists in the map
1594 TFormula * f = nullptr;
1595 auto elem = fgFormulaMap.find(replaceformula );
1596 if (elem != fgFormulaMap.end() )
1597 f = elem->second;
1598 else {
1599 // create a new formula and add in the static map
1600 f = new TFormula("f", replaceformula.Data());
1601 {
1603 fgFormulaMap[replaceformula]=f;
1604 }
1605 }
1606 if (!f) {
1607 Error("TLinearFitter", "f_linear not allocated");
1608 return;
1609 }
1610 special=f->GetNumber();
1611 fFunctions.Add(f);
1612 }
1613
1614 if ((fNfunctions==1)&&(special>299)&&(special<310)){
1615 //if fitting a polynomial
1616 size=special-299;
1617 fSpecial=100+size;
1618 } else
1620 oa->Delete();
1621 delete oa;
1622 }
1623 fNfunctions=size;
1624 //change the size of design matrix
1625 fDesign.ResizeTo(size, size);
1626 fAtb.ResizeTo(size);
1627 fDesignTemp.ResizeTo(size, size);
1628 fDesignTemp2.ResizeTo(size, size);
1629 fDesignTemp3.ResizeTo(size, size);
1630 fAtbTemp.ResizeTo(size);
1631 fAtbTemp2.ResizeTo(size);
1632 fAtbTemp3.ResizeTo(size);
1633 if (fFixedParams)
1634 delete [] fFixedParams;
1635 fFixedParams=new Bool_t[size];
1636 fDesign.Zero();
1637 fAtb.Zero();
1638 fDesignTemp.Zero();
1639 fDesignTemp2.Zero();
1640 fDesignTemp3.Zero();
1641 fAtbTemp.Zero();
1642 fAtbTemp2.Zero();
1643 fAtbTemp3.Zero();
1644 fY2Temp=0;
1645 fY2=0;
1646 for (i=0; i<size; i++)
1647 fFixedParams[i]=0;
1648 fIsSet=kFALSE;
1649 fChisquare=0;
1650
1651}
1652
1653////////////////////////////////////////////////////////////////////////////////
1654///Set the fitting function.
1655
1657{
1658 Int_t special, size;
1661 fSpecial = 0;
1662 special=fInputFunction->GetNumber();
1663 if (!fFunctions.IsEmpty())
1665
1666 if ((special>299)&&(special<310)){
1667 //if fitting a polynomial
1668 size=special-299;
1669 fSpecial=100+size;
1670 } else
1672
1674 //change the size of design matrix
1679
1682
1685 //
1686 if (fFixedParams)
1687 delete [] fFixedParams;
1689 fDesign.Zero();
1690 fAtb.Zero();
1691 fDesignTemp.Zero();
1692 fAtbTemp.Zero();
1693
1696
1697 fAtbTemp2.Zero();
1698 fAtbTemp3.Zero();
1699 fY2Temp=0;
1700 fY2=0;
1701 for (Int_t i=0; i<size; i++)
1702 fFixedParams[i]=0;
1703 //check if any parameters are fixed (not for pure TFormula)
1704
1705 if (function->InheritsFrom(TF1::Class())){
1706 Double_t al,bl;
1707 for (Int_t i=0;i<fNfunctions;i++) {
1708 ((TF1*)function)->GetParLimits(i,al,bl);
1709 if (al*bl !=0 && al >= bl) {
1710 FixParameter(i, function->GetParameter(i));
1711 }
1712 }
1713 }
1714
1715 fIsSet=kFALSE;
1716 fChisquare=0;
1717
1718}
1719
1720////////////////////////////////////////////////////////////////////////////////
1721///Update the design matrix after the formula has been changed.
1722
1724{
1725 if (fStoreData) {
1726 for (Int_t i=0; i<fNpoints; i++) {
1727 AddToDesign(TMatrixDRow(fX, i).GetPtr(), fY(i), fE(i));
1728 }
1729 return 1;
1730 } else
1731 return 0;
1732
1733}
1734
1735////////////////////////////////////////////////////////////////////////////////
1736///To use in TGraph::Fit and TH1::Fit().
1737
1738Int_t TLinearFitter::ExecuteCommand(const char *command, Double_t *args, Int_t /*nargs*/)
1739{
1740 if (!strcmp(command, "FitGraph")){
1741 if (args) return GraphLinearFitter(args[0]);
1742 else return GraphLinearFitter(0);
1743 }
1744 if (!strcmp(command, "FitGraph2D")){
1745 if (args) return Graph2DLinearFitter(args[0]);
1746 else return Graph2DLinearFitter(0);
1747 }
1748 if (!strcmp(command, "FitMultiGraph")){
1749 if (args) return MultiGraphLinearFitter(args[0]);
1750 else return MultiGraphLinearFitter(0);
1751 }
1752 if (!strcmp(command, "FitHist")) return HistLinearFitter();
1753// if (!strcmp(command, "FitMultiGraph")) MultiGraphLinearFitter();
1754
1755 return 0;
1756}
1757
1758////////////////////////////////////////////////////////////////////////////////
1759/// Level = 3 (to be consistent with minuit) prints parameters and parameter
1760/// errors.
1761
1762void TLinearFitter::PrintResults(Int_t level, Double_t /*amin*/) const
1763{
1764 if (level==3){
1765 if (!fRobust){
1766 printf("Fitting results:\nParameters:\nNO.\t\tVALUE\t\tERROR\n");
1767 for (Int_t i=0; i<fNfunctions; i++){
1768 printf("%d\t%e\t%e\n", i, fParams(i), TMath::Sqrt(fParCovar(i, i)));
1769 }
1770 } else {
1771 printf("Fitting results:\nParameters:\nNO.\t\tVALUE\n");
1772 for (Int_t i=0; i<fNfunctions; i++){
1773 printf("%d\t%e\n", i, fParams(i));
1774 }
1775 }
1776 }
1777}
1778
1779////////////////////////////////////////////////////////////////////////////////
1780///Used in TGraph::Fit().
1781
1783{
1785 TGraph *grr=(TGraph*)GetObjectFit();
1786 TF1 *f1=(TF1*)GetUserFunc();
1787 Foption_t fitOption=GetFitOption();
1788
1789 //Int_t np=0;
1790 Double_t *x=grr->GetX();
1791 Double_t *y=grr->GetY();
1792 Double_t e;
1793
1794 Int_t fitResult = 0;
1795 //set the fitting formula
1796 SetDim(1);
1798
1799 if (fitOption.Robust){
1800 fRobust=kTRUE;
1802 }
1803 //put the points into the fitter
1804 Int_t n=grr->GetN();
1805 for (Int_t i=0; i<n; i++){
1806 if (!f1->IsInside(&x[i])) continue;
1807 e=grr->GetErrorY(i);
1808 if (e<0 || fitOption.W1)
1809 e=1;
1810 AddPoint(&x[i], y[i], e);
1811 }
1812
1813 if (fitOption.Robust){
1814 return EvalRobust(h);
1815 }
1816
1817 fitResult = Eval();
1818
1819 //calculate the precise chisquare
1820 if (!fitResult){
1821 if (!fitOption.Nochisq){
1822 Double_t temp, temp2, sumtotal=0;
1823 for (Int_t i=0; i<n; i++){
1824 if (!f1->IsInside(&x[i])) continue;
1825 temp=f1->Eval(x[i]);
1826 temp2=(y[i]-temp)*(y[i]-temp);
1827 e=grr->GetErrorY(i);
1828 if (e<0 || fitOption.W1)
1829 e=1;
1830 temp2/=(e*e);
1831
1832 sumtotal+=temp2;
1833 }
1834 fChisquare=sumtotal;
1836 }
1837 }
1838 return fitResult;
1839}
1840
1841////////////////////////////////////////////////////////////////////////////////
1842///Minimisation function for a TGraph2D
1843
1845{
1847
1849 TF2 *f2=(TF2*)GetUserFunc();
1850
1851 Foption_t fitOption=GetFitOption();
1852 Int_t n = gr->GetN();
1853 Double_t *gx = gr->GetX();
1854 Double_t *gy = gr->GetY();
1855 Double_t *gz = gr->GetZ();
1856 Double_t x[2];
1857 Double_t z, e;
1858 Int_t fitResult=0;
1859 SetDim(2);
1860 SetFormula(f2->GetFormula());
1861
1862 if (fitOption.Robust){
1863 fRobust=kTRUE;
1865 }
1866
1867 for (Int_t bin=0;bin<n;bin++) {
1868 x[0] = gx[bin];
1869 x[1] = gy[bin];
1870 if (!f2->IsInside(x)) {
1871 continue;
1872 }
1873 z = gz[bin];
1874 e=gr->GetErrorZ(bin);
1875 if (e<0 || fitOption.W1)
1876 e=1;
1877 AddPoint(x, z, e);
1878 }
1879
1880 if (fitOption.Robust){
1881 return EvalRobust(h);
1882 }
1883
1884 fitResult = Eval();
1885
1886 if (!fitResult){
1887 if (!fitOption.Nochisq){
1888 //calculate the precise chisquare
1889 Double_t temp, temp2, sumtotal=0;
1890 for (Int_t bin=0; bin<n; bin++){
1891 x[0] = gx[bin];
1892 x[1] = gy[bin];
1893 if (!f2->IsInside(x)) {
1894 continue;
1895 }
1896 z = gz[bin];
1897
1898 temp=f2->Eval(x[0], x[1]);
1899 temp2=(z-temp)*(z-temp);
1900 e=gr->GetErrorZ(bin);
1901 if (e<0 || fitOption.W1)
1902 e=1;
1903 temp2/=(e*e);
1904
1905 sumtotal+=temp2;
1906 }
1907 fChisquare=sumtotal;
1909 }
1910 }
1911 return fitResult;
1912}
1913
1914////////////////////////////////////////////////////////////////////////////////
1915///Minimisation function for a TMultiGraph
1916
1918{
1919 Int_t n, i;
1920 Double_t *gx, *gy;
1921 Double_t e;
1923 TMultiGraph *mg = (TMultiGraph*)grFitter->GetObjectFit();
1924 TF1 *f1 = (TF1*)grFitter->GetUserFunc();
1925 Foption_t fitOption = grFitter->GetFitOption();
1926 Int_t fitResult=0;
1927 SetDim(1);
1928
1929 if (fitOption.Robust){
1930 fRobust=kTRUE;
1932 }
1934
1935 TGraph *gr;
1936 TIter next(mg->GetListOfGraphs());
1937 while ((gr = (TGraph*) next())) {
1938 n = gr->GetN();
1939 gx = gr->GetX();
1940 gy = gr->GetY();
1941 for (i=0; i<n; i++){
1942 if (!f1->IsInside(&gx[i])) continue;
1943 e=gr->GetErrorY(i);
1944 if (e<0 || fitOption.W1)
1945 e=1;
1946 AddPoint(&gx[i], gy[i], e);
1947 }
1948 }
1949
1950 if (fitOption.Robust){
1951 return EvalRobust(h);
1952 }
1953
1954 fitResult = Eval();
1955
1956 //calculate the chisquare
1957 if (!fitResult){
1958 if (!fitOption.Nochisq){
1959 Double_t temp, temp2, sumtotal=0;
1960 next.Reset();
1961 while((gr = (TGraph*)next())) {
1962 n = gr->GetN();
1963 gx = gr->GetX();
1964 gy = gr->GetY();
1965 for (i=0; i<n; i++){
1966 if (!f1->IsInside(&gx[i])) continue;
1967 temp=f1->Eval(gx[i]);
1968 temp2=(gy[i]-temp)*(gy[i]-temp);
1969 e=gr->GetErrorY(i);
1970 if (e<0 || fitOption.W1)
1971 e=1;
1972 temp2/=(e*e);
1973
1974 sumtotal+=temp2;
1975 }
1976
1977 }
1978 fChisquare=sumtotal;
1980 }
1981 }
1982 return fitResult;
1983}
1984
1985////////////////////////////////////////////////////////////////////////////////
1986/// Minimization function for H1s using a Chisquare method.
1987
1989{
1991 Double_t cu,eu;
1992 // Double_t dersum[100], grad[100];
1993 Double_t x[3];
1994 Int_t bin,binx,biny,binz;
1995 // Axis_t binlow, binup, binsize;
1996
1997 TH1 *hfit = (TH1*)GetObjectFit();
1998 TF1 *f1 = (TF1*)GetUserFunc();
1999 Int_t fitResult;
2000 Foption_t fitOption = GetFitOption();
2001 //SetDim(hfit->GetDimension());
2002 SetDim(f1->GetNdim());
2004
2005 Int_t hxfirst = GetXfirst();
2006 Int_t hxlast = GetXlast();
2007 Int_t hyfirst = GetYfirst();
2008 Int_t hylast = GetYlast();
2009 Int_t hzfirst = GetZfirst();
2010 Int_t hzlast = GetZlast();
2011 TAxis *xaxis = hfit->GetXaxis();
2012 TAxis *yaxis = hfit->GetYaxis();
2013 TAxis *zaxis = hfit->GetZaxis();
2014
2015 for (binz=hzfirst;binz<=hzlast;binz++) {
2016 x[2] = zaxis->GetBinCenter(binz);
2017 for (biny=hyfirst;biny<=hylast;biny++) {
2018 x[1] = yaxis->GetBinCenter(biny);
2019 for (binx=hxfirst;binx<=hxlast;binx++) {
2020 x[0] = xaxis->GetBinCenter(binx);
2021 if (!f1->IsInside(x)) continue;
2022 bin = hfit->GetBin(binx,biny,binz);
2023 cu = hfit->GetBinContent(bin);
2024 if (f1->GetNdim() < hfit->GetDimension()) {
2025 if (hfit->GetDimension() > 2) cu = x[2];
2026 else cu = x[1];
2027 }
2028 if (fitOption.W1) {
2029 if (fitOption.W1==1 && cu == 0) continue;
2030 eu = 1;
2031 } else {
2032 eu = hfit->GetBinError(bin);
2033 if (eu <= 0) continue;
2034 }
2035 AddPoint(x, cu, eu);
2036 }
2037 }
2038 }
2039
2040 fitResult = Eval();
2041
2042 if (!fitResult){
2043 if (!fitOption.Nochisq){
2044 Double_t temp, temp2, sumtotal=0;
2045 for (binz=hzfirst;binz<=hzlast;binz++) {
2046 x[2] = zaxis->GetBinCenter(binz);
2047 for (biny=hyfirst;biny<=hylast;biny++) {
2048 x[1] = yaxis->GetBinCenter(biny);
2049 for (binx=hxfirst;binx<=hxlast;binx++) {
2050 x[0] = xaxis->GetBinCenter(binx);
2051 if (!f1->IsInside(x)) continue;
2052 bin = hfit->GetBin(binx,biny,binz);
2053 cu = hfit->GetBinContent(bin);
2054
2055 if (fitOption.W1) {
2056 if (fitOption.W1==1 && cu == 0) continue;
2057 eu = 1;
2058 } else {
2059 eu = hfit->GetBinError(bin);
2060 if (eu <= 0) continue;
2061 }
2062 temp=f1->EvalPar(x);
2063 temp2=(cu-temp)*(cu-temp);
2064 temp2/=(eu*eu);
2065 sumtotal+=temp2;
2066 }
2067 }
2068 }
2069
2070 fChisquare=sumtotal;
2072 }
2073 }
2074 return fitResult;
2075}
2076
2077////////////////////////////////////////////////////////////////////////////////
2078
2080{
2081 if (R__b.IsReading()) {
2082 Int_t old_matr_size = fNfunctions;
2084 if (old_matr_size < fNfunctions) {
2087
2090
2093 }
2094 } else {
2095 if (fAtb.NonZeros()==0) AddTempMatrices();
2097 }
2098}
2099
2100////////////////////////////////////////////////////////////////////////////////
2101///Finds the parameters of the fitted function in case data contains
2102///outliers.
2103///Parameter h stands for the minimal fraction of good points in the
2104///dataset (h < 1, i.e. for 70% of good points take h=0.7).
2105///The default value of h*Npoints is (Npoints + Nparameters+1)/2
2106///If the user provides a value of h smaller than above, default is taken
2107///See class description for the algorithm details
2108
2110{
2111 fRobust = kTRUE;
2112 Double_t kEps = 1e-13;
2113 Int_t nmini = 300;
2114 Int_t i, j, maxind=0, k, k1 = 500;
2115 Int_t nbest = 10;
2116 Double_t chi2 = -1;
2117
2118 if (fFunctions.IsEmpty()&&(!fInputFunction)&&(fSpecial<=200)){
2119 Error("TLinearFitter::EvalRobust", "The formula hasn't been set");
2120 return 1;
2121 }
2122
2123
2124 Double_t *bestchi2 = new Double_t[nbest];
2125 for (i=0; i<nbest; i++)
2126 bestchi2[i]=1e30;
2127
2128 Int_t hdef=Int_t((fNpoints+fNfunctions+1)/2);
2129
2130 if (h>0.000001 && h<1 && fNpoints*h > hdef)
2131 fH = Int_t(fNpoints*h);
2132 else {
2133 // print message only when h is not zero
2134 if (h>0) Warning("Fitting:", "illegal value of H, default is taken, h = %3.2f",double(hdef)/fNpoints);
2135 fH=hdef;
2136 }
2137
2141
2142 Int_t *index = new Int_t[fNpoints];
2143 Double_t *residuals = new Double_t[fNpoints];
2144
2145 if (fNpoints < 2*nmini) {
2146 //when number of cases is small
2147
2148 //to store the best coefficients (columnwise)
2149 TMatrixD cstock(fNfunctions, nbest);
2150 for (k = 0; k < k1; k++) {
2152 chi2 = CStep(1, fH, residuals,index, index, -1, -1);
2153 chi2 = CStep(2, fH, residuals,index, index, -1, -1);
2154 maxind = TMath::LocMax(nbest, bestchi2);
2155 if (chi2 < bestchi2[maxind]) {
2156 bestchi2[maxind] = chi2;
2157 for (i=0; i<fNfunctions; i++)
2158 cstock(i, maxind) = fParams(i);
2159 }
2160 }
2161
2162 //for the nbest best results, perform CSteps until convergence
2163 Int_t *bestindex = new Int_t[fH];
2164 Double_t currentbest;
2165 for (i=0; i<nbest; i++) {
2166 for (j=0; j<fNfunctions; j++)
2167 fParams(j) = cstock(j, i);
2168 chi2 = 1;
2169 while (chi2 > kEps) {
2170 chi2 = CStep(2, fH, residuals,index, index, -1, -1);
2171 if (TMath::Abs(chi2 - bestchi2[i]) < kEps)
2172 break;
2173 else
2174 bestchi2[i] = chi2;
2175 }
2176 currentbest = TMath::MinElement(nbest, bestchi2);
2177 if (chi2 <= currentbest + kEps) {
2178 for (j=0; j<fH; j++){
2179 bestindex[j]=index[j];
2180 }
2181 maxind = i;
2182 }
2183 for (j=0; j<fNfunctions; j++)
2184 cstock(j, i) = fParams(j);
2185 }
2186 //report the result with the lowest chisquare
2187 for (j=0; j<fNfunctions; j++)
2188 fParams(j) = cstock(j, maxind);
2190 for (j=0; j<fH; j++){
2191 //printf("bestindex[%d]=%d\n", j, bestindex[j]);
2192 fFitsample.SetBitNumber(bestindex[j]);
2193 }
2195 ((TF1*)fInputFunction)->SetChisquare(bestchi2[maxind]);
2196 ((TF1*)fInputFunction)->SetNumberFitPoints(fH);
2197 ((TF1*)fInputFunction)->SetNDF(fH-fNfunctions);
2198 }
2199 delete [] index;
2200 delete [] bestindex;
2201 delete [] residuals;
2202 delete [] bestchi2;
2203 return 0;
2204 }
2205 //if n is large, the dataset should be partitioned
2206 Int_t indsubdat[5];
2207 for (i=0; i<5; i++)
2208 indsubdat[i] = 0;
2209
2210 Int_t nsub = Partition(nmini, indsubdat);
2211 Int_t hsub;
2212
2213 Int_t sum = TMath::Min(nmini*5, fNpoints);
2214
2215 Int_t *subdat = new Int_t[sum]; //to store the indices of selected cases
2216 RDraw(subdat, indsubdat);
2217
2218 TMatrixD cstockbig(fNfunctions, nbest*5);
2219 Int_t *beststock = new Int_t[nbest];
2220 Int_t i_start = 0;
2221 Int_t i_end = indsubdat[0];
2222 Int_t k2 = Int_t(k1/nsub);
2223 for (Int_t kgroup = 0; kgroup < nsub; kgroup++) {
2224
2225 hsub = Int_t(fH * indsubdat[kgroup]/fNpoints);
2226 for (i=0; i<nbest; i++)
2227 bestchi2[i] = 1e16;
2228 for (k=0; k<k2; k++) {
2229 CreateSubset(indsubdat[kgroup], hsub, index);
2230 chi2 = CStep(1, hsub, residuals, index, subdat, i_start, i_end);
2231 chi2 = CStep(2, hsub, residuals, index, subdat, i_start, i_end);
2232 maxind = TMath::LocMax(nbest, bestchi2);
2233 if (chi2 < bestchi2[maxind]){
2234 for (i=0; i<fNfunctions; i++)
2235 cstockbig(i, nbest*kgroup + maxind) = fParams(i);
2236 bestchi2[maxind] = chi2;
2237 }
2238 }
2239 if (kgroup != nsub - 1){
2240 i_start += indsubdat[kgroup];
2241 i_end += indsubdat[kgroup+1];
2242 }
2243 }
2244
2245 for (i=0; i<nbest; i++)
2246 bestchi2[i] = 1e30;
2247 //on the pooled subset
2248 Int_t hsub2 = Int_t(fH*sum/fNpoints);
2249 for (k=0; k<nbest*5; k++) {
2250 for (i=0; i<fNfunctions; i++)
2251 fParams(i)=cstockbig(i, k);
2252 chi2 = CStep(1, hsub2, residuals, index, subdat, 0, sum);
2253 chi2 = CStep(2, hsub2, residuals, index, subdat, 0, sum);
2254 maxind = TMath::LocMax(nbest, bestchi2);
2255 if (chi2 < bestchi2[maxind]){
2256 beststock[maxind] = k;
2257 bestchi2[maxind] = chi2;
2258 }
2259 }
2260
2261 //now the array beststock keeps indices of 10 best candidates in cstockbig matrix
2262 for (k=0; k<nbest; k++) {
2263 for (i=0; i<fNfunctions; i++)
2264 fParams(i) = cstockbig(i, beststock[k]);
2265 chi2 = CStep(1, fH, residuals, index, index, -1, -1);
2266 chi2 = CStep(2, fH, residuals, index, index, -1, -1);
2267 bestchi2[k] = chi2;
2268 }
2269
2270 maxind = TMath::LocMin(nbest, bestchi2);
2271 for (i=0; i<fNfunctions; i++)
2272 fParams(i)=cstockbig(i, beststock[maxind]);
2273
2274 chi2 = 1;
2275 while (chi2 > kEps) {
2276 chi2 = CStep(2, fH, residuals, index, index, -1, -1);
2277 if (TMath::Abs(chi2 - bestchi2[maxind]) < kEps)
2278 break;
2279 else
2280 bestchi2[maxind] = chi2;
2281 }
2282
2284 for (j=0; j<fH; j++)
2286 if (fInputFunction){
2287 ((TF1*)fInputFunction)->SetChisquare(bestchi2[maxind]);
2288 ((TF1*)fInputFunction)->SetNumberFitPoints(fH);
2289 ((TF1*)fInputFunction)->SetNDF(fH-fNfunctions);
2290 }
2291
2292 delete [] subdat;
2293 delete [] beststock;
2294 delete [] bestchi2;
2295 delete [] residuals;
2296 delete [] index;
2297
2298 return 0;
2299}
2300
2301////////////////////////////////////////////////////////////////////////////////
2302///Creates a p-subset to start
2303///ntotal - total number of points from which the subset is chosen
2304
2306{
2307 Int_t i, j;
2308 Bool_t repeat=kFALSE;
2309 Int_t nindex=0;
2310 Int_t num;
2311 for(i=0; i<ntotal; i++)
2312 index[i] = ntotal+1;
2313
2314 TRandom2 r;
2315 //create a p-subset
2316 for (i=0; i<fNfunctions; i++) {
2317 num=Int_t(r.Uniform(0, 1)*(ntotal-1));
2318 if (i>0){
2319 for(j=0; j<=i-1; j++) {
2320 if(index[j]==num)
2321 repeat = kTRUE;
2322 }
2323 }
2324 if(repeat==kTRUE) {
2325 i--;
2326 repeat = kFALSE;
2327 } else {
2328 index[i] = num;
2329 nindex++;
2330 }
2331 }
2332
2333 //compute the coefficients of a hyperplane through the p-subset
2334 fDesign.Zero();
2335 fAtb.Zero();
2336 for (i=0; i<fNfunctions; i++){
2337 AddToDesign(TMatrixDRow(fX, index[i]).GetPtr(), fY(index[i]), fE(index[i]));
2338 }
2339 Bool_t ok;
2340
2341 ok = Linf();
2342
2343 //if the chosen points don't define a hyperplane, add more
2344 while (!ok && (nindex < h)) {
2345 repeat=kFALSE;
2346 do{
2347 num=Int_t(r.Uniform(0,1)*(ntotal-1));
2348 repeat=kFALSE;
2349 for(i=0; i<nindex; i++) {
2350 if(index[i]==num) {
2351 repeat=kTRUE;
2352 break;
2353 }
2354 }
2355 } while(repeat==kTRUE);
2356
2357 index[nindex] = num;
2358 nindex++;
2359 //check if the system is of full rank now
2360 AddToDesign(TMatrixDRow(fX, index[nindex-1]).GetPtr(), fY(index[nindex-1]), fE(index[nindex-1]));
2361 ok = Linf();
2362 }
2363}
2364
2365////////////////////////////////////////////////////////////////////////////////
2366///The CStep procedure, as described in the article
2367
2368Double_t TLinearFitter::CStep(Int_t step, Int_t h, Double_t *residuals, Int_t *index, Int_t *subdat, Int_t start, Int_t end)
2369{
2371
2372 Int_t i, j, itemp, n;
2373 Double_t func;
2374 Double_t val[100];
2375 Int_t npar;
2376 if (start > -1) {
2377 n = end - start;
2378 for (i=0; i<n; i++) {
2379 func = 0;
2380 itemp = subdat[start+i];
2381 if (fInputFunction){
2383 func=fInputFunction->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
2384 } else {
2385 func=0;
2386 if ((fSpecial>100)&&(fSpecial<200)){
2387 npar = fSpecial-100;
2388 val[0] = 1;
2389 for (j=1; j<npar; j++)
2390 val[j] = val[j-1]*fX(itemp, 0);
2391 for (j=0; j<npar; j++)
2392 func += fParams(j)*val[j];
2393 } else if (fSpecial>200) {
2394 //hyperplane case
2395 npar = fSpecial-201;
2396 func+=fParams(0);
2397 for (j=0; j<npar; j++)
2398 func += fParams(j+1)*fX(itemp, j);
2399 } else {
2400 // general case
2401 for (j=0; j<fNfunctions; j++) {
2402 TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
2403 val[j] = f1->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
2404 func += fParams(j)*val[j];
2405 }
2406 }
2407 }
2408 residuals[i] = (fY(itemp) - func)*(fY(itemp) - func)/(fE(i)*fE(i));
2409 }
2410 } else {
2411 n=fNpoints;
2412 for (i=0; i<fNpoints; i++) {
2413 func = 0;
2414 if (fInputFunction){
2416 func=fInputFunction->EvalPar(TMatrixDRow(fX, i).GetPtr());
2417 } else {
2418 func=0;
2419 if ((fSpecial>100)&&(fSpecial<200)){
2420 npar = fSpecial-100;
2421 val[0] = 1;
2422 for (j=1; j<npar; j++)
2423 val[j] = val[j-1]*fX(i, 0);
2424 for (j=0; j<npar; j++)
2425 func += fParams(j)*val[j];
2426 } else if (fSpecial>200) {
2427 //hyperplane case
2428 npar = fSpecial-201;
2429 func+=fParams(0);
2430 for (j=0; j<npar; j++)
2431 func += fParams(j+1)*fX(i, j);
2432 } else {
2433 for (j=0; j<fNfunctions; j++) {
2434 TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
2435 val[j] = f1->EvalPar(TMatrixDRow(fX, i).GetPtr());
2436 func += fParams(j)*val[j];
2437 }
2438 }
2439 }
2440 residuals[i] = (fY(i) - func)*(fY(i) - func)/(fE(i)*fE(i));
2441 }
2442 }
2443 //take h with smallest residuals
2444 TMath::KOrdStat(n, residuals, h-1, index);
2445 //add them to the design matrix
2446 fDesign.Zero();
2447 fAtb.Zero();
2448 for (i=0; i<h; i++)
2449 AddToDesign(TMatrixDRow(fX, index[i]).GetPtr(), fY(index[i]), fE(index[i]));
2450
2451 Linf();
2452
2453 //don't calculate the chisquare at the 1st cstep
2454 if (step==1) return 0;
2455 Double_t sum=0;
2456
2457
2458 if (start > -1) {
2459 for (i=0; i<h; i++) {
2460 itemp = subdat[start+index[i]];
2461 if (fInputFunction){
2463 func=fInputFunction->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
2464 } else {
2465 func=0;
2466 if ((fSpecial>100)&&(fSpecial<200)){
2467 npar = fSpecial-100;
2468 val[0] = 1;
2469 for (j=1; j<npar; j++)
2470 val[j] = val[j-1]*fX(itemp, 0);
2471 for (j=0; j<npar; j++)
2472 func += fParams(j)*val[j];
2473 } else if (fSpecial>200) {
2474 //hyperplane case
2475 npar = fSpecial-201;
2476 func+=fParams(0);
2477 for (j=0; j<npar; j++)
2478 func += fParams(j+1)*fX(itemp, j);
2479 } else {
2480 for (j=0; j<fNfunctions; j++) {
2481 TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
2482 val[j] = f1->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
2483 func += fParams(j)*val[j];
2484 }
2485 }
2486 }
2487 sum+=(fY(itemp)-func)*(fY(itemp)-func)/(fE(itemp)*fE(itemp));
2488 }
2489 } else {
2490 for (i=0; i<h; i++) {
2491 if (fInputFunction){
2493 func=fInputFunction->EvalPar(TMatrixDRow(fX, index[i]).GetPtr());
2494 } else {
2495 func=0;
2496 if ((fSpecial>100)&&(fSpecial<200)){
2497 npar = fSpecial-100;
2498 val[0] = 1;
2499 for (j=1; j<npar; j++)
2500 val[j] = val[j-1]*fX(index[i], 0);
2501 for (j=0; j<npar; j++)
2502 func += fParams(j)*val[j];
2503 } else {
2504 if (fSpecial>200) {
2505 //hyperplane case
2506 npar = fSpecial-201;
2507 func+=fParams(0);
2508 for (j=0; j<npar; j++)
2509 func += fParams(j+1)*fX(index[i], j);
2510 } else {
2511 for (j=0; j<fNfunctions; j++) {
2512 TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
2513 val[j] = f1->EvalPar(TMatrixDRow(fX, index[i]).GetPtr());
2514 func += fParams(j)*val[j];
2515 }
2516 }
2517 }
2518 }
2519
2520 sum+=(fY(index[i])-func)*(fY(index[i])-func)/(fE(index[i])*fE(index[i]));
2521 }
2522 }
2523
2524 return sum;
2525}
2526
2527////////////////////////////////////////////////////////////////////////////////
2528
2530{
2531 //currently without the intercept term
2537 fDesignTemp.Zero();
2540 fAtb+=fAtbTemp;
2541 fAtbTemp3.Zero();
2542 fAtbTemp2.Zero();
2543 fAtbTemp.Zero();
2544
2545 fY2+=fY2Temp;
2546 fY2Temp=0;
2547
2548
2549 TDecompChol chol(fDesign);
2550 TVectorD temp(fNfunctions);
2551 Bool_t ok;
2552 temp = chol.Solve(fAtb, ok);
2553 if (!ok){
2554 Error("Linf","Matrix inversion failed");
2555 //fDesign.Print();
2556 fParams.Zero();
2557 return kFALSE;
2558 }
2559 fParams = temp;
2560 return ok;
2561}
2562
2563////////////////////////////////////////////////////////////////////////////////
2564///divides the elements into approximately equal subgroups
2565///number of elements in each subgroup is stored in indsubdat
2566///number of subgroups is returned
2567
2569{
2570 Int_t nsub;
2571
2572 if ((fNpoints>=2*nmini) && (fNpoints<=(3*nmini-1))) {
2573 if (fNpoints%2==1){
2574 indsubdat[0]=Int_t(fNpoints*0.5);
2575 indsubdat[1]=Int_t(fNpoints*0.5)+1;
2576 } else
2577 indsubdat[0]=indsubdat[1]=Int_t(fNpoints/2);
2578 nsub=2;
2579 }
2580 else{
2581 if((fNpoints>=3*nmini) && (fNpoints<(4*nmini -1))) {
2582 if(fNpoints%3==0){
2583 indsubdat[0]=indsubdat[1]=indsubdat[2]=Int_t(fNpoints/3);
2584 } else {
2585 indsubdat[0]=Int_t(fNpoints/3);
2586 indsubdat[1]=Int_t(fNpoints/3)+1;
2587 if (fNpoints%3==1) indsubdat[2]=Int_t(fNpoints/3);
2588 else indsubdat[2]=Int_t(fNpoints/3)+1;
2589 }
2590 nsub=3;
2591 }
2592 else{
2593 if((fNpoints>=4*nmini)&&(fNpoints<=(5*nmini-1))){
2594 if (fNpoints%4==0) indsubdat[0]=indsubdat[1]=indsubdat[2]=indsubdat[3]=Int_t(fNpoints/4);
2595 else {
2596 indsubdat[0]=Int_t(fNpoints/4);
2597 indsubdat[1]=Int_t(fNpoints/4)+1;
2598 if(fNpoints%4==1) indsubdat[2]=indsubdat[3]=Int_t(fNpoints/4);
2599 if(fNpoints%4==2) {
2600 indsubdat[2]=Int_t(fNpoints/4)+1;
2601 indsubdat[3]=Int_t(fNpoints/4);
2602 }
2603 if(fNpoints%4==3) indsubdat[2]=indsubdat[3]=Int_t(fNpoints/4)+1;
2604 }
2605 nsub=4;
2606 } else {
2607 for(Int_t i=0; i<5; i++)
2608 indsubdat[i]=nmini;
2609 nsub=5;
2610 }
2611 }
2612 }
2613 return nsub;
2614}
2615
2616////////////////////////////////////////////////////////////////////////////////
2617///Draws ngroup nonoverlapping subdatasets out of a dataset of size n
2618///such that the selected case numbers are uniformly distributed from 1 to n
2619
2620void TLinearFitter::RDraw(Int_t *subdat, Int_t *indsubdat)
2621{
2622 Int_t jndex = 0;
2623 Int_t nrand;
2624 Int_t i, k, m, j;
2625 Int_t ngroup=0;
2626 for (i=0; i<5; i++) {
2627 if (indsubdat[i]!=0)
2628 ngroup++;
2629 }
2630 TRandom r;
2631 for (k=1; k<=ngroup; k++) {
2632 for (m=1; m<=indsubdat[k-1]; m++) {
2633 nrand = Int_t(r.Uniform(0, 1) * (fNpoints-jndex)) + 1;
2634 jndex++;
2635 if (jndex==1) {
2636 subdat[0] = nrand;
2637 } else {
2638 subdat[jndex-1] = nrand + jndex - 2;
2639 for (i=1; i<=jndex-1; i++) {
2640 if(subdat[i-1] > nrand+i-2) {
2641 for(j=jndex; j>=i+1; j--) {
2642 subdat[j-1] = subdat[j-2];
2643 }
2644 subdat[i-1] = nrand+i-2;
2645 break; //breaking the loop for(i=1...
2646 }
2647 }
2648 }
2649 }
2650 }
2651}
2652
2653
#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
bool Bool_t
Definition: RtypesCore.h:63
int Int_t
Definition: RtypesCore.h:45
const Bool_t kFALSE
Definition: RtypesCore.h:101
const Bool_t kTRUE
Definition: RtypesCore.h:100
const char Option_t
Definition: RtypesCore.h:66
#define ClassImp(name)
Definition: Rtypes.h:375
#define R__ASSERT(e)
Definition: TError.h:118
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)
#define snprintf
Definition: civetweb.c:1540
Class to manage histogram axis.
Definition: TAxis.h:30
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:479
Int_t GetLast() const
Return last bin on the axis i.e.
Definition: TAxis.cxx:470
Int_t GetFirst() const
Return first bin on the axis i.e.
Definition: TAxis.cxx:459
Container of bits.
Definition: TBits.h:26
void Clear(Option_t *option="") override
Clear the value.
Definition: TBits.cxx:84
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
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
Bool_t IsReading() const
Definition: TBuffer.h:86
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
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.
TObject * Clone(const char *newname="") const override
Make a clone of an collection using the Streamer facility.
Cholesky Decomposition class.
Definition: TDecompChol.h:25
Bool_t Solve(TVectorD &b) override
Solve equations Ax=b assuming A has been factored by Cholesky.
Bool_t Invert(TMatrixDSym &inv)
For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
1-Dim function class
Definition: TF1.h:213
static TClass * Class()
virtual void SetChisquare(Double_t chi2)
Definition: TF1.h:617
virtual Double_t EvalPar(const Double_t *x, const Double_t *params=0)
Evaluate function with given coordinates and parameters.
Definition: TF1.cxx:1483
virtual TFormula * GetFormula()
Definition: TF1.h:458
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:1454
virtual Bool_t IsInside(const Double_t *x) const
return kTRUE if the point is inside the function range
Definition: TF1.h:603
virtual Int_t GetNdim() const
Definition: TF1.h:490
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:688
The Formula class.
Definition: TFormula.h:87
const TObject * GetLinearPart(Int_t i) const
Return linear part.
Definition: TFormula.cxx:2550
Double_t Eval(Double_t x) const
Sets first variable (e.g. x) and evaluate formula.
Definition: TFormula.cxx:3436
Double_t EvalPar(const Double_t *x, const Double_t *params=0) const
Definition: TFormula.cxx:3125
static TClass * Class()
const char * GetParName(Int_t ipar) const
Return parameter name given by integer.
Definition: TFormula.cxx:2856
Int_t GetNumber() const
Definition: TFormula.h:258
void SetParameters(const Double_t *params)
Set a vector of parameters value.
Definition: TFormula.cxx:2967
Int_t GetNpar() const
Definition: TFormula.h:257
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 * GetY() const
Definition: TGraph2D.h:122
Double_t * GetX() const
Definition: TGraph2D.h:121
virtual Double_t * GetEZ() const
Definition: TGraph2D.h:126
Int_t GetN() const
Definition: TGraph2D.h:120
virtual void SetPoint(Int_t point, Double_t x, Double_t y, Double_t z)
Sets point number n.
Definition: TGraph2D.cxx:1680
Double_t GetErrorY(Int_t bin) const override
It returns the error along Y at point i.
Double_t * GetEY() const override
Definition: TGraphErrors.h:68
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:2282
Double_t * GetY() const
Definition: TGraph.h:134
Int_t GetN() const
Definition: TGraph.h:126
Double_t * GetX() const
Definition: TGraph.h:133
virtual Double_t GetErrorY(Int_t bin) const
It always returns a negative value. Real implementation in TGraphErrors.
Definition: TGraph.cxx:1344
TH1 is the base class of all histogram classes in ROOT.
Definition: TH1.h:58
TAxis * GetZaxis()
Definition: TH1.h:321
static TClass * Class()
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8938
virtual Int_t GetDimension() const
Definition: TH1.h:281
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:319
virtual Int_t GetBin(Int_t binx, Int_t biny=0, Int_t binz=0) const
Return Global bin number corresponding to binx,y,z.
Definition: TH1.cxx:4935
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
Definition: TH1.cxx:9081
TAxis * GetYaxis()
Definition: TH1.h:320
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:9097
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:5035
void Reset()
Definition: TCollection.h:254
The Linear Fitter - For fitting functions that are LINEAR IN PARAMETERS.
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 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...
Double_t fY2Temp
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.
TLinearFitter & operator=(const TLinearFitter &tlf)
Assignment operator.
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.
TVectorD fParams
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 void AssignData(Int_t npoints, Int_t xncols, Double_t *x, Double_t *y, Double_t *e=0)
This function is to use when you already have all the data in arrays and don't want to copy them into...
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
Definition: TMatrixTBase.h:123
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
Int_t GetNcols() const
Definition: TMatrixTBase.h:126
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:92
const Element * GetMatrixArray() const override
Definition: TMatrixTSym.h:189
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])
Definition: TMatrixT.cxx:1054
void Clear(Option_t *="") override
Definition: TMatrixT.h:120
const Element * GetMatrixArray() const override
Definition: TMatrixT.h:222
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...
Definition: TMatrixT.cxx:1211
A TMultiGraph is a collection of TGraph (or derived) objects.
Definition: TMultiGraph.h:34
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
An array of TObjects.
Definition: TObjArray.h:31
Int_t GetEntriesFast() const
Definition: TObjArray.h:58
virtual void Expand(Int_t newSize)
Expand or shrink the array to newSize elements.
Definition: TObjArray.cxx:387
void Clear(Option_t *option="") override
Remove all objects from the array.
Definition: TObjArray.cxx:321
Int_t GetEntries() const override
Return the number of objects in array (i.e.
Definition: TObjArray.cxx:523
void Delete(Option_t *option="") override
Remove all objects from the array AND delete all heap based objects.
Definition: TObjArray.cxx:356
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:37
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:359
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:130
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:879
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition: TObject.cxx:445
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:893
virtual TClass * IsA() const
Definition: TObject.h:230
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:136
const char * Data() const
Definition: TString.h:369
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:692
TObjArray * Tokenize(const TString &delim) const
This function is used to isolate sequential tokens in a TString.
Definition: TString.cxx:2222
TVectorT< Element > & Zero()
Set vector elements to zero.
Definition: TVectorT.cxx:453
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:294
TVectorT< Element > & Use(Int_t lwb, Int_t upb, Element *data)
Use the array data to fill the vector lwb..upb].
Definition: TVectorT.cxx:349
void Clear(Option_t *="") override
Definition: TVectorT.h:174
Int_t NonZeros() const
Compute the number of elements != 0.0.
Definition: TVectorT.cxx:620
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
TVirtualFitter & operator=(const TVirtualFitter &tvf)
assignment operator
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
int GetDimension(const TH1 *h1)
Definition: HFitImpl.cxx:50
static const double eu
Definition: Vavilov.cxx:44
void function(const Char_t *name_, T fun, const Char_t *docstring=0)
Definition: RExports.h:167
static const std::string pattern("pattern")
static constexpr double mg
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:1283
Long64_t LocMin(Long64_t n, const T *a)
Return index of array with the minimum element.
Definition: TMath.h:922
T MinElement(Long64_t n, const T *a)
Return minimum of array a of length n.
Definition: TMath.h:902
Long64_t LocMax(Long64_t n, const T *a)
Return index of array with the maximum element.
Definition: TMath.h:950
Double_t Sqrt(Double_t x)
Definition: TMath.h:641
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:176
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:2609
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:2637
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
int W1
Definition: Foption.h:36
int Nochisq
Definition: Foption.h:45
int Robust
Definition: Foption.h:48
auto * m
Definition: textangle.C:8
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345