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