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