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