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