ROOT  6.06/09
Reference Guide
TLinearFitter.cxx
Go to the documentation of this file.
1 // @(#)root/minuit:$Id$
2 // Author: Anna Kreshuk 04/03/2005
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2005, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 #include "TLinearFitter.h"
13 #include "TMath.h"
14 #include "TDecompChol.h"
15 #include "TGraph.h"
16 #include "TGraph2D.h"
17 #include "TMultiGraph.h"
18 #include "TRandom2.h"
19 #include "TObjString.h"
20 #include "TF2.h"
21 #include "TH1.h"
22 #include "TList.h"
23 #include "TClass.h"
24 #include "TROOT.h"
25 
27 
28 
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
96  as additive parts
97  --TMath functions can be used too
98  -Functions, used as additive parts, shouldn't have any parameters,
99  even if those parameters are set.
100  --for example, if normalizing a sum of a gaus(0, 1) and a
101  gaus(0, 2), don't use the built-in "gaus" of TFormula,
102  because it has parameters, take TMath::Gaus(x, 0, 1) instead.
103  -Polynomials can be used like "pol3", .."polN"
104  -If fitting a more than 3-dimensional formula, variables should
105  be numbered as follows:
106  -- x[0], x[1], x[2]... For example, to fit "1 ++ x[0] ++ x[1] ++ x[2] ++ x[3]*x[3]"
107 
108 #### 2.2 Setting the formula:
109 
110 ##### 2.2.1 If fitting a 1-2-3-dimensional formula, one can create a
111  TF123 based on a linear expression and pass this function
112  to the fitter:
113  --Example:
114 ~~~~
115  TLinearFitter *lf = new TLinearFitter();
116  TF2 *f2 = new TF2("f2", "x ++ y ++ x*x*y*y", -2, 2, -2, 2);
117  lf->SetFormula(f2);
118 ~~~~
119  --The results of the fit are then stored in the function,
120  just like when the TH1::Fit or TGraph::Fit is used
121  --A linear function of this kind is by no means different
122  from any other function, it can be drawn, evaluated, etc.
123 
124  --For multidimensional fitting, TFormulas of the form:
125  x[0]++...++x[n] can be used
126 ##### 2.2.2 There is no need to create the function if you don't want to,
127  the formula can be set by expression:
128  --Example:
129 ~~~~
130  // 2 is the number of dimensions
131  TLinearFitter *lf = new TLinearFitter(2);
132  lf->SetFormula("x ++ y ++ x*x*y*y");
133 ~~~~
134 
135 ##### 2.2.3 The fastest functions to compute are polynomials and hyperplanes.
136  --Polynomials are set the usual way: "pol1", "pol2",...
137  --Hyperplanes are set by expression "hyp3", "hyp4", ...
138  ---The "hypN" expressions only work when the linear fitter
139  is used directly, not through TH1::Fit or TGraph::Fit.
140  To fit a graph or a histogram with a hyperplane, define
141  the function as "1++x++y".
142  ---A constant term is assumed for a hyperplane, when using
143  the "hypN" expression, so "hyp3" is in fact fitting with
144  "1++x++y++z" function.
145  --Fitting hyperplanes is much faster than fitting other
146  expressions so if performance is vital, calculate the
147  function values beforehand and give them to the fitter
148  as variables
149  --Example:
150  You want to fit "sin(x)|cos(2*x)" very fast. Calculate
151  sin(x) and cos(2*x) beforehand and store them in array *data.
152  Then:
153  TLinearFitter *lf=new TLinearFitter(2, "hyp2");
154  lf->AssignData(npoint, 2, data, y);
155 
156 #### 2.3 Resetting the formula
157 
158 ##### 2.3.1 If the input data is stored (or added via AssignData() function),
159  the fitting formula can be reset without re-adding all the points.
160  --Example:
161 ~~~~
162  TLinearFitter *lf=new TLinearFitter("1++x++x*x");
163  lf->AssignData(n, 1, x, y, e);
164  lf->Eval()
165  //looking at the parameter significance, you see,
166  // that maybe the fit will improve, if you take out
167  // the constant term
168  lf->SetFormula("x++x*x");
169  lf->Eval();
170  ...
171 ~~~~
172 
173 ##### 2.3.2 If the input data is not stored, the fitter will have to be
174  cleared and the data will have to be added again to try a
175  different formula.
176 
177 ### 3.Accessing the fit results
178 
179 #### 3.1 There are methods in the fitter to access all relevant information:
180  --GetParameters, GetCovarianceMatrix, etc
181  --the t-values of parameters and their significance can be reached by
182  GetParTValue() and GetParSignificance() methods
183 
184 #### 3.2 If fitting with a pre-defined TF123, the fit results are also
185  written into this function.
186 
187 
188 ### 4.Robust fitting - Least Trimmed Squares regression (LTS)
189  Outliers are atypical(by definition), infrequant observations; data points
190  which do not appear to follow the characteristic distribution of the rest
191  of the data. These may reflect genuine properties of the underlying
192  phenomenon(variable), or be due to measurement errors or anomalies which
193  shouldn't be modelled. (StatSoft electronic textbook)
194 
195  Even a single gross outlier can greatly influence the results of least-
196  squares fitting procedure, and in this case use of robust(resistant) methods
197  is recommended.
198 
199  The method implemented here is based on the article and algorithm:
200  "Computing LTS Regression for Large Data Sets" by
201  P.J.Rousseeuw and Katrien Van Driessen
202  The idea of the method is to find the fitting coefficients for a subset
203  of h observations (out of n) with the smallest sum of squared residuals.
204  The size of the subset h should lie between (npoints + nparameters +1)/2
205  and n, and represents the minimal number of good points in the dataset.
206  The default value is set to (npoints + nparameters +1)/2, but of course
207  if you are sure that the data contains less outliers it's better to change
208  h according to your data.
209 
210  To perform a robust fit, call EvalRobust() function instead of Eval() after
211  adding the points and setting the fitting function.
212  Note, that standard errors on parameters are not computed!
213 
214 */
215 
216 
217 
218 ////////////////////////////////////////////////////////////////////////////////
219 ///default c-tor, input data is stored
220 ///If you don't want to store the input data,
221 ///run the function StoreData(kFALSE) after constructor
222 
225  fParams(),
226  fParCovar(),
227  fTValues(),
228  fDesign(),
229  fDesignTemp(),
230  fDesignTemp2(),
231  fDesignTemp3(),
232  fAtb(),
233  fAtbTemp(),
234  fAtbTemp2(),
235  fAtbTemp3(),
236  fFunctions(),
237  fY(),
238  fX(),
239  fE(),
240  fVal()
241 {
242  fChisquare =0;
243  fNpoints =0;
244  fNdim =0;
245  fY2 =0;
246  fY2Temp =0;
247  fNfixed =0;
248  fIsSet =kFALSE;
249  fFormula =0;
250  fFixedParams=0;
251  fSpecial =0;
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){
597  Error("AddData", "Those points are already added");
598  return;
599  }
600  Bool_t same=kFALSE;
601  if (fX.GetMatrixArray()==x && fY.GetMatrixArray()==y){
602  if (e){
603  if (fE.GetMatrixArray()==e)
604  same=kTRUE;
605  }
606  }
607 
608  fX.Use(npoints, xncols, x);
609  fY.Use(npoints, y);
610  if (e)
611  fE.Use(npoints, e);
612  else {
613  fE.ResizeTo(npoints);
614  fE=1;
615  }
616  Int_t xfirst;
617  if (!fFunctions.IsEmpty() || fInputFunction || fSpecial>200) {
618  if (same)
619  xfirst=fNpoints;
620 
621  else
622  xfirst=0;
623  for (Int_t i=xfirst; i<npoints; i++)
624  AddToDesign(TMatrixDRow(fX, i).GetPtr(), fY(i), fE(i));
625  }
626  fNpoints=npoints;
627 }
628 
629 ////////////////////////////////////////////////////////////////////////////////
630 ///Add a point to the AtA matrix and to the Atb vector.
631 
633 {
634 
635 
636  Int_t i, j, ii;
637  y/=e;
638 
639 // Double_t fVal[1000];
640 
641  if ((fSpecial>100)&&(fSpecial<200)){
642  //polynomial fitting
643  Int_t npar=fSpecial-100;
644  fVal[0]=1;
645  for (i=1; i<npar; i++)
646  fVal[i]=fVal[i-1]*x[0];
647  for (i=0; i<npar; i++)
648  fVal[i]/=e;
649  } else if (fSpecial>200){
650  //Hyperplane fitting. Constant term is added
651  Int_t npar=fSpecial-201;
652  fVal[0]=1./e;
653  for (i=0; i<npar; i++)
654  fVal[i+1]=x[i]/e;
655  } else {
656  //general case
657  for (ii=0; ii<fNfunctions; ii++){
658  if (!fFunctions.IsEmpty()){
659  // ffunctions can be TF1 or TFormula depending on how they are created
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)));
877  fParSign(i) = 2*(1-TMath::StudentI(TMath::Abs(fTValues(i)),fNpoints-fNfunctions+fNfixed));
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 
921  AddTempMatrices();
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());
984  ((TF1*)fInputFunction)->SetNDF(fNpoints-fNfunctions+fNfixed);
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())
1079  v = fAtb;
1080 }
1081 
1082 ////////////////////////////////////////////////////////////////////////////////
1083 /// Get the Chisquare.
1084 
1086 {
1087  if (fChisquare > 1e-16)
1088  return fChisquare;
1089  else {
1090  Chisquare();
1091  return fChisquare;
1092  }
1093 }
1094 
1095 ////////////////////////////////////////////////////////////////////////////////
1096 ///Computes point-by-point confidence intervals for the fitted function
1097 ///Parameters:
1098 ///n - number of points
1099 ///ndim - dimensions of points
1100 ///x - points, at which to compute the intervals, for ndim > 1
1101 /// should be in order: (x0,y0, x1, y1, ... xn, yn)
1102 ///ci - computed intervals are returned in this array
1103 ///cl - confidence level, default=0.95
1104 ///
1105 ///NOTE, that this method can only be used when the fitting function inherits from a TF1,
1106 ///so it's not possible when the fitting function was set as a string or as a pure TFormula
1107 
1109 {
1110  if (fInputFunction){
1111  Double_t *grad = new Double_t[fNfunctions];
1112  Double_t *sum_vector = new Double_t[fNfunctions];
1113  Double_t c=0;
1115  Double_t t = TMath::StudentQuantile(0.5 + cl/2, df);
1116  Double_t chidf = TMath::Sqrt(fChisquare/df);
1117 
1118  for (Int_t ipoint=0; ipoint<n; ipoint++){
1119  c=0;
1120  ((TF1*)(fInputFunction))->GradientPar(x+ndim*ipoint, grad);
1121  //multiply the covariance matrix by gradient
1122  for (Int_t irow=0; irow<fNfunctions; irow++){
1123  sum_vector[irow]=0;
1124  for (Int_t icol=0; icol<fNfunctions; icol++){
1125  sum_vector[irow]+=fParCovar(irow,icol)*grad[icol];
1126  }
1127  }
1128  for (Int_t i=0; i<fNfunctions; i++)
1129  c+=grad[i]*sum_vector[i];
1130  c=TMath::Sqrt(c);
1131  ci[ipoint]=c*t*chidf;
1132  }
1133 
1134  delete [] grad;
1135  delete [] sum_vector;
1136  }
1137 }
1138 
1139 ////////////////////////////////////////////////////////////////////////////////
1140 ///Computes confidence intervals at level cl. Default is 0.95
1141 ///The TObject parameter can be a TGraphErrors, a TGraph2DErrors or a TH123.
1142 ///For Graphs, confidence intervals are computed for each point,
1143 ///the value of the graph at that point is set to the function value at that
1144 ///point, and the graph y-errors (or z-errors) are set to the value of
1145 ///the confidence interval at that point
1146 ///For Histograms, confidence intervals are computed for each bin center
1147 ///The bin content of this bin is then set to the function value at the bin
1148 ///center, and the bin error is set to the confidence interval value.
1149 ///Allowed combinations:
1150 ///Fitted object Passed object
1151 ///TGraph TGraphErrors, TH1
1152 ///TGraphErrors, AsymmErrors TGraphErrors, TH1
1153 ///TH1 TGraphErrors, TH1
1154 ///TGraph2D TGraph2DErrors, TH2
1155 ///TGraph2DErrors TGraph2DErrors, TH2
1156 ///TH2 TGraph2DErrors, TH2
1157 ///TH3 TH3
1158 
1160 {
1161  if (!fInputFunction) {
1162  Error("GetConfidenceIntervals", "The case of fitting not with a TFormula is not yet implemented");
1163  return;
1164  }
1165 
1166  //TGraph//////////////////
1167 
1168  if (obj->InheritsFrom(TGraph::Class())) {
1169  TGraph *gr = (TGraph*)obj;
1170  if (!gr->GetEY()){
1171  Error("GetConfidenceIntervals", "A TGraphErrors should be passed instead of a graph");
1172  return;
1173  }
1175  Error("GetConfidenceIntervals", "A TGraph2DErrors should be passed instead of a graph");
1176  return;
1177  }
1179  if (((TH1*)(fObjectFit))->GetDimension()>1){
1180  Error("GetConfidenceIntervals", "A TGraph2DErrors or a TH23 should be passed instead of a graph");
1181  return;
1182  }
1183  }
1184 
1185  GetConfidenceIntervals(gr->GetN(),1,gr->GetX(), gr->GetEY(), cl);
1186  for (Int_t i=0; i<gr->GetN(); i++)
1187  gr->SetPoint(i, gr->GetX()[i], fInputFunction->Eval(gr->GetX()[i]));
1188  }
1189 
1190  //TGraph2D///////////////
1191  else if (obj->InheritsFrom(TGraph2D::Class())) {
1192  TGraph2D *gr2 = (TGraph2D*)obj;
1193  if (!gr2->GetEZ()){
1194  Error("GetConfidenceIntervals", "A TGraph2DErrors should be passed instead of a TGraph2D");
1195  return;
1196  }
1198  Error("GetConfidenceIntervals", "A TGraphErrors should be passed instead of a TGraph2D");
1199  return;
1200  }
1202  if (((TH1*)(fObjectFit))->GetDimension()==1){
1203  Error("GetConfidenceIntervals", "A TGraphErrors or a TH1 should be passed instead of a graph");
1204  return;
1205  }
1206  }
1207  Double_t xy[2];
1208  Int_t np = gr2->GetN();
1209  Double_t *grad = new Double_t[fNfunctions];
1210  Double_t *sum_vector = new Double_t[fNfunctions];
1211  Double_t *x = gr2->GetX();
1212  Double_t *y = gr2->GetY();
1213  Double_t t = TMath::StudentQuantile(0.5 + cl/2, ((TF1*)(fInputFunction))->GetNDF());
1214  Double_t chidf = TMath::Sqrt(fChisquare/((TF1*)(fInputFunction))->GetNDF());
1215  Double_t c = 0;
1216  for (Int_t ipoint=0; ipoint<np; ipoint++){
1217  c=0;
1218  xy[0]=x[ipoint];
1219  xy[1]=y[ipoint];
1220  ((TF1*)(fInputFunction))->GradientPar(xy, grad);
1221  for (Int_t irow=0; irow<fNfunctions; irow++){
1222  sum_vector[irow]=0;
1223  for (Int_t icol=0; icol<fNfunctions; icol++)
1224  sum_vector[irow]+=fParCovar(irow, icol)*grad[icol];
1225  }
1226  for (Int_t i=0; i<fNfunctions; i++)
1227  c+=grad[i]*sum_vector[i];
1228  c=TMath::Sqrt(c);
1229  gr2->SetPoint(ipoint, xy[0], xy[1], fInputFunction->EvalPar(xy));
1230  gr2->GetEZ()[ipoint]=c*t*chidf;
1231  }
1232  delete [] grad;
1233  delete [] sum_vector;
1234  }
1235 
1236  //TH1////////////////////////
1237  else if (obj->InheritsFrom(TH1::Class())) {
1239  if (((TH1*)obj)->GetDimension()>1){
1240  Error("GetConfidenceIntervals", "Fitted graph and passed histogram have different number of dimensions");
1241  return;
1242  }
1243  }
1245  if (((TH1*)obj)->GetDimension()!=2){
1246  Error("GetConfidenceIntervals", "Fitted graph and passed histogram have different number of dimensions");
1247  return;
1248  }
1249  }
1251  if (((TH1*)(fObjectFit))->GetDimension()!=((TH1*)(obj))->GetDimension()){
1252  Error("GetConfidenceIntervals", "Fitted and passed histograms have different number of dimensions");
1253  return;
1254  }
1255  }
1256 
1257 
1258  TH1 *hfit = (TH1*)obj;
1259  Double_t *grad = new Double_t[fNfunctions];
1260  Double_t *sum_vector = new Double_t[fNfunctions];
1261  Double_t x[3];
1262  Int_t hxfirst = hfit->GetXaxis()->GetFirst();
1263  Int_t hxlast = hfit->GetXaxis()->GetLast();
1264  Int_t hyfirst = hfit->GetYaxis()->GetFirst();
1265  Int_t hylast = hfit->GetYaxis()->GetLast();
1266  Int_t hzfirst = hfit->GetZaxis()->GetFirst();
1267  Int_t hzlast = hfit->GetZaxis()->GetLast();
1268 
1269  TAxis *xaxis = hfit->GetXaxis();
1270  TAxis *yaxis = hfit->GetYaxis();
1271  TAxis *zaxis = hfit->GetZaxis();
1272  Double_t t = TMath::StudentQuantile(0.5 + cl/2, ((TF1*)(fInputFunction))->GetNDF());
1273  Double_t chidf = TMath::Sqrt(fChisquare/((TF1*)(fInputFunction))->GetNDF());
1274  Double_t c=0;
1275  for (Int_t binz=hzfirst; binz<=hzlast; binz++){
1276  x[2]=zaxis->GetBinCenter(binz);
1277  for (Int_t biny=hyfirst; biny<=hylast; biny++) {
1278  x[1]=yaxis->GetBinCenter(biny);
1279  for (Int_t binx=hxfirst; binx<=hxlast; binx++) {
1280  x[0]=xaxis->GetBinCenter(binx);
1281  ((TF1*)(fInputFunction))->GradientPar(x, grad);
1282  c=0;
1283  for (Int_t irow=0; irow<fNfunctions; irow++){
1284  sum_vector[irow]=0;
1285  for (Int_t icol=0; icol<fNfunctions; icol++)
1286  sum_vector[irow]+=fParCovar(irow, icol)*grad[icol];
1287  }
1288  for (Int_t i=0; i<fNfunctions; i++)
1289  c+=grad[i]*sum_vector[i];
1290  c=TMath::Sqrt(c);
1291  hfit->SetBinContent(binx, biny, binz, fInputFunction->EvalPar(x));
1292  hfit->SetBinError(binx, biny, binz, c*t*chidf);
1293  }
1294  }
1295  }
1296  delete [] grad;
1297  delete [] sum_vector;
1298  }
1299  else {
1300  Error("GetConfidenceIntervals", "This object type is not supported");
1301  return;
1302  }
1303 }
1304 
1305 ////////////////////////////////////////////////////////////////////////////////
1306 ///Returns covariance matrix
1307 
1309 {
1310  Double_t *p = const_cast<Double_t*>(fParCovar.GetMatrixArray());
1311  return p;
1312 }
1313 
1314 ////////////////////////////////////////////////////////////////////////////////
1315 ///Returns covariance matrix
1316 
1318 {
1319  if (matr.GetNrows()!=fNfunctions || matr.GetNcols()!=fNfunctions){
1321  }
1322  matr = fParCovar;
1323 }
1324 
1325 ////////////////////////////////////////////////////////////////////////////////
1326 ///Returns the internal design matrix
1327 
1329 {
1330  if (matr.GetNrows()!=fNfunctions || matr.GetNcols()!=fNfunctions){
1332  }
1333  matr = fDesign;
1334 }
1335 
1336 ////////////////////////////////////////////////////////////////////////////////
1337 ///Returns parameter errors
1338 
1340 {
1341  if (vpar.GetNoElements()!=fNfunctions) {
1342  vpar.ResizeTo(fNfunctions);
1343  }
1344  for (Int_t i=0; i<fNfunctions; i++)
1345  vpar(i) = TMath::Sqrt(fParCovar(i, i));
1346 
1347 }
1348 
1349 ////////////////////////////////////////////////////////////////////////////////
1350 ///Returns parameter values
1351 
1353 {
1354  if (vpar.GetNoElements()!=fNfunctions) {
1355  vpar.ResizeTo(fNfunctions);
1356  }
1357  vpar=fParams;
1358 }
1359 
1360 ////////////////////////////////////////////////////////////////////////////////
1361 ///Returns the value and the name of the parameter #ipar
1362 ///NB: In the calling function the argument name must be set large enough
1363 
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())) {
1459  Error("Add","Attempt to add object of class: %s to a %s",lfit->ClassName(),this->ClassName());
1460  return -1;
1461  }
1462  Add(lfit);
1463  }
1464  return 0;
1465 }
1466 ////////////////////////////////////////////////////////////////////////////////
1467 ///set the basis functions in case the fitting function is not
1468 /// set directly
1469 /// The TLinearFitter will manage and delete the functions contained in the list
1470 
1472 {
1473  fFunctions = *(functions);
1475  int size = fFunctions.GetEntries();
1476 
1477  fNfunctions=size;
1478  //change the size of design matrix
1479  fDesign.ResizeTo(size, size);
1480  fAtb.ResizeTo(size);
1481  fDesignTemp.ResizeTo(size, size);
1482  fDesignTemp2.ResizeTo(size, size);
1483  fDesignTemp3.ResizeTo(size, size);
1484  fAtbTemp.ResizeTo(size);
1485  fAtbTemp2.ResizeTo(size);
1486  fAtbTemp3.ResizeTo(size);
1487  if (fFixedParams)
1488  delete [] fFixedParams;
1489  fFixedParams=new Bool_t[size];
1490  fDesign.Zero();
1491  fAtb.Zero();
1492  fDesignTemp.Zero();
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 = 0;
1539  fFormulaSize = strlen(formula);
1540  fFormula = new char[fFormulaSize+1];
1541  strlcpy(fFormula, formula,fFormulaSize+1);
1542  fSpecial = 0;
1543  //in case of a hyperplane:
1544  char *fstring;
1545  fstring = (char *)strstr(fFormula, "hyp");
1546  if (fstring!=0){
1547  // isHyper = kTRUE;
1548  fstring+=3;
1549  sscanf(fstring, "%d", &size);
1550  //+1 for the constant term
1551  size++;
1552  fSpecial=200+size;
1553  }
1554 
1555  if (fSpecial==0) {
1556  //in case it's not a hyperplane
1557  TString sstring(fFormula);
1558  sstring = sstring.ReplaceAll("++", 2, "|", 1);
1559  TString replaceformula;
1560 
1561  //count the number of functions
1562 
1563  TObjArray *oa = sstring.Tokenize("|");
1564 
1565  //change the size of functions array and clear it
1566  if (!fFunctions.IsEmpty())
1567  fFunctions.Clear();
1568 
1569  // do not own the functions in this case
1571 
1572  fNfunctions = oa->GetEntriesFast();
1574 
1575  //check if the old notation of xi is used somewhere instead of x[i]
1576  char pattern[5];
1577  char replacement[6];
1578  for (i=0; i<fNdim; i++){
1579  snprintf(pattern,5, "x%d", i);
1580  snprintf(replacement,6, "x[%d]", i);
1581  sstring = sstring.ReplaceAll(pattern, Int_t(i/10)+2, replacement, Int_t(i/10)+4);
1582  }
1583 
1584  //fill the array of functions
1585  oa = sstring.Tokenize("|");
1586  for (i=0; i<fNfunctions; i++) {
1587  replaceformula = ((TObjString *)oa->UncheckedAt(i))->GetString();
1588  // look first if exists in the map
1589  TFormula * f = nullptr;
1590  auto elem = fgFormulaMap.find(replaceformula );
1591  if (elem != fgFormulaMap.end() )
1592  f = elem->second;
1593  else {
1594  // create a new formula and add in the static map
1595  f = new TFormula("f", replaceformula.Data());
1596  {
1598  fgFormulaMap[replaceformula]=f;
1599  }
1600  }
1601  if (!f) {
1602  Error("TLinearFitter", "f_linear not allocated");
1603  return;
1604  }
1605  special=f->GetNumber();
1606  fFunctions.Add(f);
1607  }
1608 
1609  if ((fNfunctions==1)&&(special>299)&&(special<310)){
1610  //if fitting a polynomial
1611  size=special-299;
1612  fSpecial=100+size;
1613  } else
1614  size=fNfunctions;
1615  oa->Delete();
1616  delete oa;
1617  }
1618  fNfunctions=size;
1619  //change the size of design matrix
1620  fDesign.ResizeTo(size, size);
1621  fAtb.ResizeTo(size);
1622  fDesignTemp.ResizeTo(size, size);
1623  fDesignTemp2.ResizeTo(size, size);
1624  fDesignTemp3.ResizeTo(size, size);
1625  fAtbTemp.ResizeTo(size);
1626  fAtbTemp2.ResizeTo(size);
1627  fAtbTemp3.ResizeTo(size);
1628  if (fFixedParams)
1629  delete [] fFixedParams;
1630  fFixedParams=new Bool_t[size];
1631  fDesign.Zero();
1632  fAtb.Zero();
1633  fDesignTemp.Zero();
1634  fDesignTemp2.Zero();
1635  fDesignTemp3.Zero();
1636  fAtbTemp.Zero();
1637  fAtbTemp2.Zero();
1638  fAtbTemp3.Zero();
1639  fY2Temp=0;
1640  fY2=0;
1641  for (i=0; i<size; i++)
1642  fFixedParams[i]=0;
1643  fIsSet=kFALSE;
1644  fChisquare=0;
1645 
1646 }
1647 
1648 ////////////////////////////////////////////////////////////////////////////////
1649 ///Set the fitting function.
1650 
1652 {
1653  Int_t special, size;
1654  fInputFunction=function;
1656  fSpecial = 0;
1657  special=fInputFunction->GetNumber();
1658  if (!fFunctions.IsEmpty())
1659  fFunctions.Delete();
1660 
1661  if ((special>299)&&(special<310)){
1662  //if fitting a polynomial
1663  size=special-299;
1664  fSpecial=100+size;
1665  } else
1666  size=fNfunctions;
1667 
1668  fNfunctions=size;
1669  //change the size of design matrix
1670  fDesign.ResizeTo(size, size);
1671  fAtb.ResizeTo(size);
1672  fDesignTemp.ResizeTo(size, size);
1673  fAtbTemp.ResizeTo(size);
1674 
1675  fDesignTemp2.ResizeTo(size, size);
1676  fDesignTemp3.ResizeTo(size, size);
1677 
1678  fAtbTemp2.ResizeTo(size);
1679  fAtbTemp3.ResizeTo(size);
1680  //
1681  if (fFixedParams)
1682  delete [] fFixedParams;
1683  fFixedParams=new Bool_t[size];
1684  fDesign.Zero();
1685  fAtb.Zero();
1686  fDesignTemp.Zero();
1687  fAtbTemp.Zero();
1688 
1689  fDesignTemp2.Zero();
1690  fDesignTemp3.Zero();
1691 
1692  fAtbTemp2.Zero();
1693  fAtbTemp3.Zero();
1694  fY2Temp=0;
1695  fY2=0;
1696  for (Int_t i=0; i<size; i++)
1697  fFixedParams[i]=0;
1698  //check if any parameters are fixed (not for pure TFormula)
1699 
1700  if (function->InheritsFrom(TF1::Class())){
1701  Double_t al,bl;
1702  for (Int_t i=0;i<fNfunctions;i++) {
1703  ((TF1*)function)->GetParLimits(i,al,bl);
1704  if (al*bl !=0 && al >= bl) {
1705  FixParameter(i, function->GetParameter(i));
1706  }
1707  }
1708  }
1709 
1710  fIsSet=kFALSE;
1711  fChisquare=0;
1712 
1713 }
1714 
1715 ////////////////////////////////////////////////////////////////////////////////
1716 ///Update the design matrix after the formula has been changed.
1717 
1719 {
1720  if (fStoreData) {
1721  for (Int_t i=0; i<fNpoints; i++) {
1722  AddToDesign(TMatrixDRow(fX, i).GetPtr(), fY(i), fE(i));
1723  }
1724  return 1;
1725  } else
1726  return 0;
1727 
1728 }
1729 
1730 ////////////////////////////////////////////////////////////////////////////////
1731 ///To use in TGraph::Fit and TH1::Fit().
1732 
1733 Int_t TLinearFitter::ExecuteCommand(const char *command, Double_t *args, Int_t /*nargs*/)
1734 {
1735  if (!strcmp(command, "FitGraph")){
1736  if (args) return GraphLinearFitter(args[0]);
1737  else return GraphLinearFitter(0);
1738  }
1739  if (!strcmp(command, "FitGraph2D")){
1740  if (args) return Graph2DLinearFitter(args[0]);
1741  else return Graph2DLinearFitter(0);
1742  }
1743  if (!strcmp(command, "FitMultiGraph")){
1744  if (args) return MultiGraphLinearFitter(args[0]);
1745  else return MultiGraphLinearFitter(0);
1746  }
1747  if (!strcmp(command, "FitHist")) return HistLinearFitter();
1748 // if (!strcmp(command, "FitMultiGraph")) MultiGraphLinearFitter();
1749 
1750  return 0;
1751 }
1752 
1753 ////////////////////////////////////////////////////////////////////////////////
1754 /// Level = 3 (to be consistent with minuit) prints parameters and parameter
1755 /// errors.
1756 
1757 void TLinearFitter::PrintResults(Int_t level, Double_t /*amin*/) const
1758 {
1759  if (level==3){
1760  if (!fRobust){
1761  printf("Fitting results:\nParameters:\nNO.\t\tVALUE\t\tERROR\n");
1762  for (Int_t i=0; i<fNfunctions; i++){
1763  printf("%d\t%e\t%e\n", i, fParams(i), TMath::Sqrt(fParCovar(i, i)));
1764  }
1765  } else {
1766  printf("Fitting results:\nParameters:\nNO.\t\tVALUE\n");
1767  for (Int_t i=0; i<fNfunctions; i++){
1768  printf("%d\t%e\n", i, fParams(i));
1769  }
1770  }
1771  }
1772 }
1773 
1774 ////////////////////////////////////////////////////////////////////////////////
1775 ///Used in TGraph::Fit().
1776 
1778 {
1779  StoreData(kFALSE);
1780  TGraph *grr=(TGraph*)GetObjectFit();
1781  TF1 *f1=(TF1*)GetUserFunc();
1782  Foption_t fitOption=GetFitOption();
1783 
1784  //Int_t np=0;
1785  Double_t *x=grr->GetX();
1786  Double_t *y=grr->GetY();
1787  Double_t e;
1788 
1789  Int_t fitResult = 0;
1790  //set the fitting formula
1791  SetDim(1);
1792  SetFormula(f1->GetFormula());
1793 
1794  if (fitOption.Robust){
1795  fRobust=kTRUE;
1796  StoreData(kTRUE);
1797  }
1798  //put the points into the fitter
1799  Int_t n=grr->GetN();
1800  for (Int_t i=0; i<n; i++){
1801  if (!f1->IsInside(&x[i])) continue;
1802  e=grr->GetErrorY(i);
1803  if (e<0 || fitOption.W1)
1804  e=1;
1805  AddPoint(&x[i], y[i], e);
1806  }
1807 
1808  if (fitOption.Robust){
1809  return EvalRobust(h);
1810  }
1811 
1812  fitResult = Eval();
1813 
1814  //calculate the precise chisquare
1815  if (!fitResult){
1816  if (!fitOption.Nochisq){
1817  Double_t temp, temp2, sumtotal=0;
1818  for (Int_t i=0; i<n; i++){
1819  if (!f1->IsInside(&x[i])) continue;
1820  temp=f1->Eval(x[i]);
1821  temp2=(y[i]-temp)*(y[i]-temp);
1822  e=grr->GetErrorY(i);
1823  if (e<0 || fitOption.W1)
1824  e=1;
1825  temp2/=(e*e);
1826 
1827  sumtotal+=temp2;
1828  }
1829  fChisquare=sumtotal;
1830  f1->SetChisquare(fChisquare);
1831  }
1832  }
1833  return fitResult;
1834 }
1835 
1836 ////////////////////////////////////////////////////////////////////////////////
1837 ///Minimisation function for a TGraph2D
1838 
1840 {
1841  StoreData(kFALSE);
1842 
1844  TF2 *f2=(TF2*)GetUserFunc();
1845 
1846  Foption_t fitOption=GetFitOption();
1847  Int_t n = gr->GetN();
1848  Double_t *gx = gr->GetX();
1849  Double_t *gy = gr->GetY();
1850  Double_t *gz = gr->GetZ();
1851  Double_t x[2];
1852  Double_t z, e;
1853  Int_t fitResult=0;
1854  SetDim(2);
1855  SetFormula(f2->GetFormula());
1856 
1857  if (fitOption.Robust){
1858  fRobust=kTRUE;
1859  StoreData(kTRUE);
1860  }
1861 
1862  for (Int_t bin=0;bin<n;bin++) {
1863  x[0] = gx[bin];
1864  x[1] = gy[bin];
1865  if (!f2->IsInside(x)) {
1866  continue;
1867  }
1868  z = gz[bin];
1869  e=gr->GetErrorZ(bin);
1870  if (e<0 || fitOption.W1)
1871  e=1;
1872  AddPoint(x, z, e);
1873  }
1874 
1875  if (fitOption.Robust){
1876  return EvalRobust(h);
1877  }
1878 
1879  fitResult = Eval();
1880 
1881  if (!fitResult){
1882  if (!fitOption.Nochisq){
1883  //calculate the precise chisquare
1884  Double_t temp, temp2, sumtotal=0;
1885  for (Int_t bin=0; bin<n; bin++){
1886  x[0] = gx[bin];
1887  x[1] = gy[bin];
1888  if (!f2->IsInside(x)) {
1889  continue;
1890  }
1891  z = gz[bin];
1892 
1893  temp=f2->Eval(x[0], x[1]);
1894  temp2=(z-temp)*(z-temp);
1895  e=gr->GetErrorZ(bin);
1896  if (e<0 || fitOption.W1)
1897  e=1;
1898  temp2/=(e*e);
1899 
1900  sumtotal+=temp2;
1901  }
1902  fChisquare=sumtotal;
1903  f2->SetChisquare(fChisquare);
1904  }
1905  }
1906  return fitResult;
1907 }
1908 
1909 ////////////////////////////////////////////////////////////////////////////////
1910 ///Minimisation function for a TMultiGraph
1911 
1913 {
1914  Int_t n, i;
1915  Double_t *gx, *gy;
1916  Double_t e;
1918  TMultiGraph *mg = (TMultiGraph*)grFitter->GetObjectFit();
1919  TF1 *f1 = (TF1*)grFitter->GetUserFunc();
1920  Foption_t fitOption = grFitter->GetFitOption();
1921  Int_t fitResult=0;
1922  SetDim(1);
1923 
1924  if (fitOption.Robust){
1925  fRobust=kTRUE;
1926  StoreData(kTRUE);
1927  }
1928  SetFormula(f1->GetFormula());
1929 
1930  TGraph *gr;
1931  TIter next(mg->GetListOfGraphs());
1932  while ((gr = (TGraph*) next())) {
1933  n = gr->GetN();
1934  gx = gr->GetX();
1935  gy = gr->GetY();
1936  for (i=0; i<n; i++){
1937  if (!f1->IsInside(&gx[i])) continue;
1938  e=gr->GetErrorY(i);
1939  if (e<0 || fitOption.W1)
1940  e=1;
1941  AddPoint(&gx[i], gy[i], e);
1942  }
1943  }
1944 
1945  if (fitOption.Robust){
1946  return EvalRobust(h);
1947  }
1948 
1949  fitResult = Eval();
1950 
1951  //calculate the chisquare
1952  if (!fitResult){
1953  if (!fitOption.Nochisq){
1954  Double_t temp, temp2, sumtotal=0;
1955  next.Reset();
1956  while((gr = (TGraph*)next())) {
1957  n = gr->GetN();
1958  gx = gr->GetX();
1959  gy = gr->GetY();
1960  for (i=0; i<n; i++){
1961  if (!f1->IsInside(&gx[i])) continue;
1962  temp=f1->Eval(gx[i]);
1963  temp2=(gy[i]-temp)*(gy[i]-temp);
1964  e=gr->GetErrorY(i);
1965  if (e<0 || fitOption.W1)
1966  e=1;
1967  temp2/=(e*e);
1968 
1969  sumtotal+=temp2;
1970  }
1971 
1972  }
1973  fChisquare=sumtotal;
1974  f1->SetChisquare(fChisquare);
1975  }
1976  }
1977  return fitResult;
1978 }
1979 
1980 ////////////////////////////////////////////////////////////////////////////////
1981 /// Minimization function for H1s using a Chisquare method.
1982 
1984 {
1985  StoreData(kFALSE);
1986  Double_t cu,eu;
1987  // Double_t dersum[100], grad[100];
1988  Double_t x[3];
1989  Int_t bin,binx,biny,binz;
1990  // Axis_t binlow, binup, binsize;
1991 
1992  TH1 *hfit = (TH1*)GetObjectFit();
1993  TF1 *f1 = (TF1*)GetUserFunc();
1994  Int_t fitResult;
1995  Foption_t fitOption = GetFitOption();
1996  //SetDim(hfit->GetDimension());
1997  SetDim(f1->GetNdim());
1998  SetFormula(f1->GetFormula());
1999 
2000  Int_t hxfirst = GetXfirst();
2001  Int_t hxlast = GetXlast();
2002  Int_t hyfirst = GetYfirst();
2003  Int_t hylast = GetYlast();
2004  Int_t hzfirst = GetZfirst();
2005  Int_t hzlast = GetZlast();
2006  TAxis *xaxis = hfit->GetXaxis();
2007  TAxis *yaxis = hfit->GetYaxis();
2008  TAxis *zaxis = hfit->GetZaxis();
2009 
2010  for (binz=hzfirst;binz<=hzlast;binz++) {
2011  x[2] = zaxis->GetBinCenter(binz);
2012  for (biny=hyfirst;biny<=hylast;biny++) {
2013  x[1] = yaxis->GetBinCenter(biny);
2014  for (binx=hxfirst;binx<=hxlast;binx++) {
2015  x[0] = xaxis->GetBinCenter(binx);
2016  if (!f1->IsInside(x)) continue;
2017  bin = hfit->GetBin(binx,biny,binz);
2018  cu = hfit->GetBinContent(bin);
2019  if (f1->GetNdim() < hfit->GetDimension()) {
2020  if (hfit->GetDimension() > 2) cu = x[2];
2021  else cu = x[1];
2022  }
2023  if (fitOption.W1) {
2024  if (fitOption.W1==1 && cu == 0) continue;
2025  eu = 1;
2026  } else {
2027  eu = hfit->GetBinError(bin);
2028  if (eu <= 0) continue;
2029  }
2030  AddPoint(x, cu, eu);
2031  }
2032  }
2033  }
2034 
2035  fitResult = Eval();
2036 
2037  if (!fitResult){
2038  if (!fitOption.Nochisq){
2039  Double_t temp, temp2, sumtotal=0;
2040  for (binz=hzfirst;binz<=hzlast;binz++) {
2041  x[2] = zaxis->GetBinCenter(binz);
2042  for (biny=hyfirst;biny<=hylast;biny++) {
2043  x[1] = yaxis->GetBinCenter(biny);
2044  for (binx=hxfirst;binx<=hxlast;binx++) {
2045  x[0] = xaxis->GetBinCenter(binx);
2046  if (!f1->IsInside(x)) continue;
2047  bin = hfit->GetBin(binx,biny,binz);
2048  cu = hfit->GetBinContent(bin);
2049 
2050  if (fitOption.W1) {
2051  if (fitOption.W1==1 && cu == 0) continue;
2052  eu = 1;
2053  } else {
2054  eu = hfit->GetBinError(bin);
2055  if (eu <= 0) continue;
2056  }
2057  temp=f1->EvalPar(x);
2058  temp2=(cu-temp)*(cu-temp);
2059  temp2/=(eu*eu);
2060  sumtotal+=temp2;
2061  }
2062  }
2063  }
2064 
2065  fChisquare=sumtotal;
2066  f1->SetChisquare(fChisquare);
2067  }
2068  }
2069  return fitResult;
2070 }
2071 
2072 ////////////////////////////////////////////////////////////////////////////////
2073 
2074 void TLinearFitter::Streamer(TBuffer &R__b)
2075 {
2076  if (R__b.IsReading()) {
2077  Int_t old_matr_size = fNfunctions;
2079  if (old_matr_size < fNfunctions) {
2082 
2085 
2088  }
2089  } else {
2090  if (fAtb.NonZeros()==0) AddTempMatrices();
2092  }
2093 }
2094 
2095 ////////////////////////////////////////////////////////////////////////////////
2096 ///Finds the parameters of the fitted function in case data contains
2097 ///outliers.
2098 ///Parameter h stands for the minimal fraction of good points in the
2099 ///dataset (h < 1, i.e. for 70% of good points take h=0.7).
2100 ///The default value of h*Npoints is (Npoints + Nparameters+1)/2
2101 ///If the user provides a value of h smaller than above, default is taken
2102 ///See class description for the algorithm details
2103 
2105 {
2106  fRobust = kTRUE;
2107  Double_t kEps = 1e-13;
2108  Int_t nmini = 300;
2109  Int_t i, j, maxind=0, k, k1 = 500;
2110  Int_t nbest = 10;
2111  Double_t chi2 = -1;
2112 
2113  if (fFunctions.IsEmpty()&&(!fInputFunction)&&(fSpecial<=200)){
2114  Error("TLinearFitter::EvalRobust", "The formula hasn't been set");
2115  return 1;
2116  }
2117 
2118 
2119  Double_t *bestchi2 = new Double_t[nbest];
2120  for (i=0; i<nbest; i++)
2121  bestchi2[i]=1e30;
2122 
2123  Int_t hdef=Int_t((fNpoints+fNfunctions+1)/2);
2124 
2125  if (h>0.000001 && h<1 && fNpoints*h > hdef)
2126  fH = Int_t(fNpoints*h);
2127  else {
2128  // print message only when h is not zero
2129  if (h>0) Warning("Fitting:", "illegal value of H, default is taken, h = %3.2f",double(hdef)/fNpoints);
2130  fH=hdef;
2131  }
2132 
2136 
2137  Int_t *index = new Int_t[fNpoints];
2138  Double_t *residuals = new Double_t[fNpoints];
2139 
2140  if (fNpoints < 2*nmini) {
2141  //when number of cases is small
2142 
2143  //to store the best coefficients (columnwise)
2144  TMatrixD cstock(fNfunctions, nbest);
2145  for (k = 0; k < k1; k++) {
2146  CreateSubset(fNpoints, fH, index);
2147  chi2 = CStep(1, fH, residuals,index, index, -1, -1);
2148  chi2 = CStep(2, fH, residuals,index, index, -1, -1);
2149  maxind = TMath::LocMax(nbest, bestchi2);
2150  if (chi2 < bestchi2[maxind]) {
2151  bestchi2[maxind] = chi2;
2152  for (i=0; i<fNfunctions; i++)
2153  cstock(i, maxind) = fParams(i);
2154  }
2155  }
2156 
2157  //for the nbest best results, perform CSteps until convergence
2158  Int_t *bestindex = new Int_t[fH];
2159  Double_t currentbest;
2160  for (i=0; i<nbest; i++) {
2161  for (j=0; j<fNfunctions; j++)
2162  fParams(j) = cstock(j, i);
2163  chi2 = 1;
2164  while (chi2 > kEps) {
2165  chi2 = CStep(2, fH, residuals,index, index, -1, -1);
2166  if (TMath::Abs(chi2 - bestchi2[i]) < kEps)
2167  break;
2168  else
2169  bestchi2[i] = chi2;
2170  }
2171  currentbest = TMath::MinElement(nbest, bestchi2);
2172  if (chi2 <= currentbest + kEps) {
2173  for (j=0; j<fH; j++){
2174  bestindex[j]=index[j];
2175  }
2176  maxind = i;
2177  }
2178  for (j=0; j<fNfunctions; j++)
2179  cstock(j, i) = fParams(j);
2180  }
2181  //report the result with the lowest chisquare
2182  for (j=0; j<fNfunctions; j++)
2183  fParams(j) = cstock(j, maxind);
2185  for (j=0; j<fH; j++){
2186  //printf("bestindex[%d]=%d\n", j, bestindex[j]);
2187  fFitsample.SetBitNumber(bestindex[j]);
2188  }
2190  ((TF1*)fInputFunction)->SetChisquare(bestchi2[maxind]);
2191  ((TF1*)fInputFunction)->SetNumberFitPoints(fH);
2192  ((TF1*)fInputFunction)->SetNDF(fH-fNfunctions);
2193  }
2194  delete [] index;
2195  delete [] bestindex;
2196  delete [] residuals;
2197  delete [] bestchi2;
2198  return 0;
2199  }
2200  //if n is large, the dataset should be partitioned
2201  Int_t indsubdat[5];
2202  for (i=0; i<5; i++)
2203  indsubdat[i] = 0;
2204 
2205  Int_t nsub = Partition(nmini, indsubdat);
2206  Int_t hsub;
2207 
2208  Int_t sum = TMath::Min(nmini*5, fNpoints);
2209 
2210  Int_t *subdat = new Int_t[sum]; //to store the indices of selected cases
2211  RDraw(subdat, indsubdat);
2212 
2213  TMatrixD cstockbig(fNfunctions, nbest*5);
2214  Int_t *beststock = new Int_t[nbest];
2215  Int_t i_start = 0;
2216  Int_t i_end = indsubdat[0];
2217  Int_t k2 = Int_t(k1/nsub);
2218  for (Int_t kgroup = 0; kgroup < nsub; kgroup++) {
2219 
2220  hsub = Int_t(fH * indsubdat[kgroup]/fNpoints);
2221  for (i=0; i<nbest; i++)
2222  bestchi2[i] = 1e16;
2223  for (k=0; k<k2; k++) {
2224  CreateSubset(indsubdat[kgroup], hsub, index);
2225  chi2 = CStep(1, hsub, residuals, index, subdat, i_start, i_end);
2226  chi2 = CStep(2, hsub, residuals, index, subdat, i_start, i_end);
2227  maxind = TMath::LocMax(nbest, bestchi2);
2228  if (chi2 < bestchi2[maxind]){
2229  for (i=0; i<fNfunctions; i++)
2230  cstockbig(i, nbest*kgroup + maxind) = fParams(i);
2231  bestchi2[maxind] = chi2;
2232  }
2233  }
2234  if (kgroup != nsub - 1){
2235  i_start += indsubdat[kgroup];
2236  i_end += indsubdat[kgroup+1];
2237  }
2238  }
2239 
2240  for (i=0; i<nbest; i++)
2241  bestchi2[i] = 1e30;
2242  //on the pooled subset
2243  Int_t hsub2 = Int_t(fH*sum/fNpoints);
2244  for (k=0; k<nbest*5; k++) {
2245  for (i=0; i<fNfunctions; i++)
2246  fParams(i)=cstockbig(i, k);
2247  chi2 = CStep(1, hsub2, residuals, index, subdat, 0, sum);
2248  chi2 = CStep(2, hsub2, residuals, index, subdat, 0, sum);
2249  maxind = TMath::LocMax(nbest, bestchi2);
2250  if (chi2 < bestchi2[maxind]){
2251  beststock[maxind] = k;
2252  bestchi2[maxind] = chi2;
2253  }
2254  }
2255 
2256  //now the array beststock keeps indices of 10 best candidates in cstockbig matrix
2257  for (k=0; k<nbest; k++) {
2258  for (i=0; i<fNfunctions; i++)
2259  fParams(i) = cstockbig(i, beststock[k]);
2260  chi2 = CStep(1, fH, residuals, index, index, -1, -1);
2261  chi2 = CStep(2, fH, residuals, index, index, -1, -1);
2262  bestchi2[k] = chi2;
2263  }
2264 
2265  maxind = TMath::LocMin(nbest, bestchi2);
2266  for (i=0; i<fNfunctions; i++)
2267  fParams(i)=cstockbig(i, beststock[maxind]);
2268 
2269  chi2 = 1;
2270  while (chi2 > kEps) {
2271  chi2 = CStep(2, fH, residuals, index, index, -1, -1);
2272  if (TMath::Abs(chi2 - bestchi2[maxind]) < kEps)
2273  break;
2274  else
2275  bestchi2[maxind] = chi2;
2276  }
2277 
2279  for (j=0; j<fH; j++)
2280  fFitsample.SetBitNumber(index[j]);
2281  if (fInputFunction){
2282  ((TF1*)fInputFunction)->SetChisquare(bestchi2[maxind]);
2283  ((TF1*)fInputFunction)->SetNumberFitPoints(fH);
2284  ((TF1*)fInputFunction)->SetNDF(fH-fNfunctions);
2285  }
2286 
2287  delete [] subdat;
2288  delete [] beststock;
2289  delete [] bestchi2;
2290  delete [] residuals;
2291  delete [] index;
2292 
2293  return 0;
2294 }
2295 
2296 ////////////////////////////////////////////////////////////////////////////////
2297 ///Creates a p-subset to start
2298 ///ntotal - total number of points from which the subset is chosen
2299 
2301 {
2302  Int_t i, j;
2303  Bool_t repeat=kFALSE;
2304  Int_t nindex=0;
2305  Int_t num;
2306  for(i=0; i<ntotal; i++)
2307  index[i] = ntotal+1;
2308 
2309  TRandom2 r;
2310  //create a p-subset
2311  for (i=0; i<fNfunctions; i++) {
2312  num=Int_t(r.Uniform(0, 1)*(ntotal-1));
2313  if (i>0){
2314  for(j=0; j<=i-1; j++) {
2315  if(index[j]==num)
2316  repeat = kTRUE;
2317  }
2318  }
2319  if(repeat==kTRUE) {
2320  i--;
2321  repeat = kFALSE;
2322  } else {
2323  index[i] = num;
2324  nindex++;
2325  }
2326  }
2327 
2328  //compute the coefficients of a hyperplane through the p-subset
2329  fDesign.Zero();
2330  fAtb.Zero();
2331  for (i=0; i<fNfunctions; i++){
2332  AddToDesign(TMatrixDRow(fX, index[i]).GetPtr(), fY(index[i]), fE(index[i]));
2333  }
2334  Bool_t ok;
2335 
2336  ok = Linf();
2337 
2338  //if the chosen points don't define a hyperplane, add more
2339  while (!ok && (nindex < h)) {
2340  repeat=kFALSE;
2341  do{
2342  num=Int_t(r.Uniform(0,1)*(ntotal-1));
2343  repeat=kFALSE;
2344  for(i=0; i<nindex; i++) {
2345  if(index[i]==num) {
2346  repeat=kTRUE;
2347  break;
2348  }
2349  }
2350  } while(repeat==kTRUE);
2351 
2352  index[nindex] = num;
2353  nindex++;
2354  //check if the system is of full rank now
2355  AddToDesign(TMatrixDRow(fX, index[nindex-1]).GetPtr(), fY(index[nindex-1]), fE(index[nindex-1]));
2356  ok = Linf();
2357  }
2358 }
2359 
2360 ////////////////////////////////////////////////////////////////////////////////
2361 ///The CStep procedure, as described in the article
2362 
2363 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)
2364 {
2366 
2367  Int_t i, j, itemp, n;
2368  Double_t func;
2369  Double_t val[100];
2370  Int_t npar;
2371  if (start > -1) {
2372  n = end - start;
2373  for (i=0; i<n; i++) {
2374  func = 0;
2375  itemp = subdat[start+i];
2376  if (fInputFunction){
2378  func=fInputFunction->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
2379  } else {
2380  func=0;
2381  if ((fSpecial>100)&&(fSpecial<200)){
2382  npar = fSpecial-100;
2383  val[0] = 1;
2384  for (j=1; j<npar; j++)
2385  val[j] = val[j-1]*fX(itemp, 0);
2386  for (j=0; j<npar; j++)
2387  func += fParams(j)*val[j];
2388  } else if (fSpecial>200) {
2389  //hyperplane case
2390  npar = fSpecial-201;
2391  func+=fParams(0);
2392  for (j=0; j<npar; j++)
2393  func += fParams(j+1)*fX(itemp, j);
2394  } else {
2395  // general case
2396  for (j=0; j<fNfunctions; j++) {
2397  TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
2398  val[j] = f1->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
2399  func += fParams(j)*val[j];
2400  }
2401  }
2402  }
2403  residuals[i] = (fY(itemp) - func)*(fY(itemp) - func)/(fE(i)*fE(i));
2404  }
2405  } else {
2406  n=fNpoints;
2407  for (i=0; i<fNpoints; i++) {
2408  func = 0;
2409  if (fInputFunction){
2411  func=fInputFunction->EvalPar(TMatrixDRow(fX, i).GetPtr());
2412  } else {
2413  func=0;
2414  if ((fSpecial>100)&&(fSpecial<200)){
2415  npar = fSpecial-100;
2416  val[0] = 1;
2417  for (j=1; j<npar; j++)
2418  val[j] = val[j-1]*fX(i, 0);
2419  for (j=0; j<npar; j++)
2420  func += fParams(j)*val[j];
2421  } else if (fSpecial>200) {
2422  //hyperplane case
2423  npar = fSpecial-201;
2424  func+=fParams(0);
2425  for (j=0; j<npar; j++)
2426  func += fParams(j+1)*fX(i, j);
2427  } else {
2428  for (j=0; j<fNfunctions; j++) {
2429  TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
2430  val[j] = f1->EvalPar(TMatrixDRow(fX, i).GetPtr());
2431  func += fParams(j)*val[j];
2432  }
2433  }
2434  }
2435  residuals[i] = (fY(i) - func)*(fY(i) - func)/(fE(i)*fE(i));
2436  }
2437  }
2438  //take h with smallest residuals
2439  TMath::KOrdStat(n, residuals, h-1, index);
2440  //add them to the design matrix
2441  fDesign.Zero();
2442  fAtb.Zero();
2443  for (i=0; i<h; i++)
2444  AddToDesign(TMatrixDRow(fX, index[i]).GetPtr(), fY(index[i]), fE(index[i]));
2445 
2446  Linf();
2447 
2448  //don't calculate the chisquare at the 1st cstep
2449  if (step==1) return 0;
2450  Double_t sum=0;
2451 
2452 
2453  if (start > -1) {
2454  for (i=0; i<h; i++) {
2455  itemp = subdat[start+index[i]];
2456  if (fInputFunction){
2458  func=fInputFunction->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
2459  } else {
2460  func=0;
2461  if ((fSpecial>100)&&(fSpecial<200)){
2462  npar = fSpecial-100;
2463  val[0] = 1;
2464  for (j=1; j<npar; j++)
2465  val[j] = val[j-1]*fX(itemp, 0);
2466  for (j=0; j<npar; j++)
2467  func += fParams(j)*val[j];
2468  } else if (fSpecial>200) {
2469  //hyperplane case
2470  npar = fSpecial-201;
2471  func+=fParams(0);
2472  for (j=0; j<npar; j++)
2473  func += fParams(j+1)*fX(itemp, j);
2474  } else {
2475  for (j=0; j<fNfunctions; j++) {
2476  TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
2477  val[j] = f1->EvalPar(TMatrixDRow(fX, itemp).GetPtr());
2478  func += fParams(j)*val[j];
2479  }
2480  }
2481  }
2482  sum+=(fY(itemp)-func)*(fY(itemp)-func)/(fE(itemp)*fE(itemp));
2483  }
2484  } else {
2485  for (i=0; i<h; i++) {
2486  if (fInputFunction){
2488  func=fInputFunction->EvalPar(TMatrixDRow(fX, index[i]).GetPtr());
2489  } else {
2490  func=0;
2491  if ((fSpecial>100)&&(fSpecial<200)){
2492  npar = fSpecial-100;
2493  val[0] = 1;
2494  for (j=1; j<npar; j++)
2495  val[j] = val[j-1]*fX(index[i], 0);
2496  for (j=0; j<npar; j++)
2497  func += fParams(j)*val[j];
2498  } else {
2499  if (fSpecial>200) {
2500  //hyperplane case
2501  npar = fSpecial-201;
2502  func+=fParams(0);
2503  for (j=0; j<npar; j++)
2504  func += fParams(j+1)*fX(index[i], j);
2505  } else {
2506  for (j=0; j<fNfunctions; j++) {
2507  TF1 *f1 = (TF1*)(fFunctions.UncheckedAt(j));
2508  val[j] = f1->EvalPar(TMatrixDRow(fX, index[i]).GetPtr());
2509  func += fParams(j)*val[j];
2510  }
2511  }
2512  }
2513  }
2514 
2515  sum+=(fY(index[i])-func)*(fY(index[i])-func)/(fE(index[i])*fE(index[i]));
2516  }
2517  }
2518 
2519  return sum;
2520 }
2521 
2522 ////////////////////////////////////////////////////////////////////////////////
2523 
2525 {
2526  //currently without the intercept term
2530  fDesignTemp3.Zero();
2531  fDesignTemp2.Zero();
2532  fDesignTemp.Zero();
2535  fAtb+=fAtbTemp;
2536  fAtbTemp3.Zero();
2537  fAtbTemp2.Zero();
2538  fAtbTemp.Zero();
2539 
2540  fY2+=fY2Temp;
2541  fY2Temp=0;
2542 
2543 
2544  TDecompChol chol(fDesign);
2545  TVectorD temp(fNfunctions);
2546  Bool_t ok;
2547  temp = chol.Solve(fAtb, ok);
2548  if (!ok){
2549  Error("Linf","Matrix inversion failed");
2550  //fDesign.Print();
2551  fParams.Zero();
2552  return kFALSE;
2553  }
2554  fParams = temp;
2555  return ok;
2556 }
2557 
2558 ////////////////////////////////////////////////////////////////////////////////
2559 ///divides the elements into approximately equal subgroups
2560 ///number of elements in each subgroup is stored in indsubdat
2561 ///number of subgroups is returned
2562 
2564 {
2565  Int_t nsub;
2566 
2567  if ((fNpoints>=2*nmini) && (fNpoints<=(3*nmini-1))) {
2568  if (fNpoints%2==1){
2569  indsubdat[0]=Int_t(fNpoints*0.5);
2570  indsubdat[1]=Int_t(fNpoints*0.5)+1;
2571  } else
2572  indsubdat[0]=indsubdat[1]=Int_t(fNpoints/2);
2573  nsub=2;
2574  }
2575  else{
2576  if((fNpoints>=3*nmini) && (fNpoints<(4*nmini -1))) {
2577  if(fNpoints%3==0){
2578  indsubdat[0]=indsubdat[1]=indsubdat[2]=Int_t(fNpoints/3);
2579  } else {
2580  indsubdat[0]=Int_t(fNpoints/3);
2581  indsubdat[1]=Int_t(fNpoints/3)+1;
2582  if (fNpoints%3==1) indsubdat[2]=Int_t(fNpoints/3);
2583  else indsubdat[2]=Int_t(fNpoints/3)+1;
2584  }
2585  nsub=3;
2586  }
2587  else{
2588  if((fNpoints>=4*nmini)&&(fNpoints<=(5*nmini-1))){
2589  if (fNpoints%4==0) indsubdat[0]=indsubdat[1]=indsubdat[2]=indsubdat[3]=Int_t(fNpoints/4);
2590  else {
2591  indsubdat[0]=Int_t(fNpoints/4);
2592  indsubdat[1]=Int_t(fNpoints/4)+1;
2593  if(fNpoints%4==1) indsubdat[2]=indsubdat[3]=Int_t(fNpoints/4);
2594  if(fNpoints%4==2) {
2595  indsubdat[2]=Int_t(fNpoints/4)+1;
2596  indsubdat[3]=Int_t(fNpoints/4);
2597  }
2598  if(fNpoints%4==3) indsubdat[2]=indsubdat[3]=Int_t(fNpoints/4)+1;
2599  }
2600  nsub=4;
2601  } else {
2602  for(Int_t i=0; i<5; i++)
2603  indsubdat[i]=nmini;
2604  nsub=5;
2605  }
2606  }
2607  }
2608  return nsub;
2609 }
2610 
2611 ////////////////////////////////////////////////////////////////////////////////
2612 ///Draws ngroup nonoverlapping subdatasets out of a dataset of size n
2613 ///such that the selected case numbers are uniformly distributed from 1 to n
2614 
2615 void TLinearFitter::RDraw(Int_t *subdat, Int_t *indsubdat)
2616 {
2617  Int_t jndex = 0;
2618  Int_t nrand;
2619  Int_t i, k, m, j;
2620  Int_t ngroup=0;
2621  for (i=0; i<5; i++) {
2622  if (indsubdat[i]!=0)
2623  ngroup++;
2624  }
2625  TRandom r;
2626  for (k=1; k<=ngroup; k++) {
2627  for (m=1; m<=indsubdat[k-1]; m++) {
2628  nrand = Int_t(r.Uniform(0, 1) * (fNpoints-jndex)) + 1;
2629  jndex++;
2630  if (jndex==1) {
2631  subdat[0] = nrand;
2632  } else {
2633  subdat[jndex-1] = nrand + jndex - 2;
2634  for (i=1; i<=jndex-1; i++) {
2635  if(subdat[i-1] > nrand+i-2) {
2636  for(j=jndex; j>=i+1; j--) {
2637  subdat[j-1] = subdat[j-2];
2638  }
2639  subdat[i-1] = nrand+i-2;
2640  break; //breaking the loop for(i=1...
2641  }
2642  }
2643  }
2644  }
2645  }
2646 }
2647 
2648 
virtual void GetAtbVector(TVectorD &v)
Get the Atb vector - a vector, used for internal computations.
virtual void GetErrors(TVectorD &vpar)
Returns parameter errors.
Int_t GetFirst() const
Return first bin on the axis i.e.
Definition: TAxis.cxx:429
virtual void StoreData(Bool_t store)
static std::map< TString, TFormula * > fgFormulaMap
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...
An array of TObjects.
Definition: TObjArray.h:39
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Definition: TVectorT.cxx:291
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
virtual TFormula * GetFormula()
Definition: TF1.h:339
void Clear(Option_t *option="")
Clear the value.
Definition: TBits.cxx:80
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4629
Long64_t LocMax(Long64_t n, const T *a)
Definition: TMath.h:724
Double_t Eval(Double_t x) const
Definition: TFormula.cxx:2538
virtual void FixParameter(Int_t ipar)
Fixes paramter #ipar at its current value.
Bool_t IsReading() const
Definition: TBuffer.h:81
virtual Bool_t InheritsFrom(const char *classname) const
Returns kTRUE if object inherits from class "classname".
Definition: TObject.cxx:487
ClassImp(TSeqCollection) Int_t TSeqCollection TIter next(this)
Return index of object in collection.
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
TList * GetListOfGraphs() const
Definition: TMultiGraph.h:71
Bool_t IsEmpty() const
Definition: TObjArray.h:72
Collectable string class.
Definition: TObjString.h:32
virtual void Clear(Option_t *option="")
Remove all objects from the array.
Definition: TObjArray.cxx:297
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:328
The Linear Fitter - For fitting functions that are LINEAR IN PARAMETERS.
virtual Foption_t GetFitOption() const
Random number generator class based on the maximally quidistributed combined Tausworthe generator by ...
Definition: TRandom2.h:29
TString & ReplaceAll(const TString &s1, const TString &s2)
Definition: TString.h:635
virtual void AddTempMatrices()
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, at which the quantile is computed 2nd argument - the number of degrees of freedom of the Student distribution When the 3rd argument lower_tail is kTRUE (default)- the algorithm returns such x0, that P(x < x0)=p upper tail (lower_tail is kFALSE)- the algorithm returns such x0, that P(x > x0)=p the algorithm was taken from G.W.Hill, "Algorithm 396, Student's t-quantiles" "Communications of the ACM", 13(10), October 1970.
Definition: TMath.cxx:2553
virtual Int_t GetDimension() const
Definition: TH1.h:283
Double_t * GetX() const
Definition: TGraph2D.h:128
TObjArray fFunctions
map of basis functions and formula
void SetParameters(const Double_t *params)
Definition: TFormula.cxx:2401
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.
TH1 * h
Definition: legend2.C:5
A TMultiGraph is a collection of TGraph (or derived) objects.
Definition: TMultiGraph.h:37
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1101
Buffer base class used for serializing objects.
Definition: TBuffer.h:40
virtual void ClearPoints()
To be used when different sets of points are fitted with the same formula.
#define R__ASSERT(e)
Definition: TError.h:98
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:137
virtual TObject * GetUserFunc() const
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:170
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
R__EXTERN TVirtualMutex * gROOTMutex
Definition: TROOT.h:63
const Bool_t kFALSE
Definition: Rtypes.h:92
virtual void GetDesignMatrix(TMatrixD &matr)
Returns the internal design matrix.
Int_t Partition(Int_t nmini, Int_t *indsubdat)
divides the elements into approximately equal subgroups number of elements in each subgroup is stored...
STL namespace.
Int_t GetEntriesFast() const
Definition: TObjArray.h:66
virtual void Clear(Option_t *="")
Definition: TMatrixTSym.h:96
TFormula * fInputFunction
Int_t GetN() const
Definition: TGraph.h:132
virtual Int_t GetZlast() const
const TObject * GetLinearPart(Int_t i) const
Definition: TFormula.cxx:1968
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:1201
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
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:2529
virtual Double_t * GetEY() const
Definition: TGraph.h:142
virtual void PrintResults(Int_t level, Double_t amin=0) const
Level = 3 (to be consistent with minuit) prints parameters and parameter errors.
Double_t fY2Temp
virtual TObject * Clone(const char *newname="") const
Make a clone of an collection using the Streamer facility.
const char * Data() const
Definition: TString.h:349
Double_t * GetY() const
Definition: TGraph.h:140
TVectorT< Element > & Use(Int_t lwb, Int_t upb, Element *data)
Use the array data to fill the vector lwb..upb].
Definition: TVectorT.cxx:346
Double_t x[n]
Definition: legend1.C:17
Double_t EvalPar(const Double_t *x, const Double_t *params=0) const
Definition: TFormula.cxx:2509
void Class()
Definition: Class.C:29
int GetDimension(const TH1 *h1)
Definition: HFitImpl.cxx:49
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:29
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
UChar_t mod R__LOCKGUARD2(gSrvAuthenticateMutex)
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 void Clear(Option_t *option="")
Clears everything. Used in TH1::Fit and TGraph::Fit().
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:4535
static const std::string pattern("pattern")
Bool_t Invert(TMatrixDSym &inv)
For a symmetric matrix A(m,m), its inverse A_inv(m,m) is returned .
TMatrixDSym fDesign
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:918
void Clear(Option_t *="")
Definition: TVectorT.h:180
virtual void SetBinError(Int_t bin, Double_t error)
see convention for numbering bins in TH1::GetBin
Definition: TH1.cxx:8528
virtual Double_t * GetEZ() const
Definition: TGraph2D.h:133
virtual const Element * GetMatrixArray() const
Definition: TMatrixTSym.h:183
TMatrixTRow< Double_t > TMatrixDRow
virtual void Chisquare()
Calculates the chisquare.
Int_t GetNrows() const
Definition: TMatrixTBase.h:134
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:84
virtual Double_t GetParTValue(Int_t ipar)
Returns the t-value for parameter #ipar.
TObject * UncheckedAt(Int_t i) const
Definition: TObjArray.h:91
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 void SetChisquare(Double_t chi2)
Definition: TF1.h:419
TVectorD fTValues
virtual Int_t GetYlast() const
void function(const char *name_, T fun, const char *docstring=0)
Definition: RExports.h:159
Int_t NonZeros() const
Compute the number of elements != 0.0.
Definition: TVectorT.cxx:617
Double_t * GetX() const
Definition: TGraph.h:139
TVectorD fAtbTemp
Int_t HistLinearFitter()
Minimization function for H1s using a Chisquare method.
Int_t GetNpar() const
Definition: TFormula.h:175
ROOT::R::TRInterface & r
Definition: Object.C:4
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:36
virtual Double_t GetParSignificance(Int_t ipar)
Returns the significance of parameter #ipar.
SVector< double, 2 > v
Definition: Dict.h:5
XPoint xy[kMAXMK]
Definition: TGX11.cxx:122
Int_t MultiGraphLinearFitter(Double_t h)
Minimisation function for a TMultiGraph.
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:187
The F O R M U L A class.
Definition: TFormula.h:89
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:8543
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:48
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)
Double_t * GetZ() const
Definition: TGraph2D.h:130
TVectorT< Element > & Zero()
Set vector elements to zero.
Definition: TVectorT.cxx:450
TMarker * m
Definition: textangle.C:8
Int_t GraphLinearFitter(Double_t h)
Used in TGraph::Fit().
virtual void ReleaseParameter(Int_t ipar)
Releases parameter #ipar.
virtual Int_t GetNdim() const
Definition: TF1.h:350
TLinearFitter()
default c-tor, input data is stored If you don't want to store the input data, run the function Store...
static TVirtualFitter * GetFitter()
static: return the current Fitter
TMatrixD fX
temporary variable used for num.stability
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
virtual const Element * GetMatrixArray() const
Definition: TMatrixT.h:211
TAxis * GetYaxis()
Definition: TH1.h:320
A 2-Dim function with parameters.
Definition: TF2.h:33
virtual Double_t GetParameter(Int_t ipar) const
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
TObjArray * Tokenize(const TString &delim) const
This function is used to isolate sequential tokens in a TString.
Definition: TString.cxx:2240
int W1
Definition: Foption.h:36
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
virtual Int_t GetXlast() const
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
int Robust
Definition: Foption.h:48
#define ClassImp(name)
Definition: Rtypes.h:279
double f(double x)
virtual void Add(TLinearFitter *tlf)
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:1707
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:415
double Double_t
Definition: RtypesCore.h:55
Int_t GetNumber() const
Definition: TFormula.h:176
virtual void Clear(Option_t *="")
Definition: TMatrixT.h:121
TVectorD fParSign
TLinearFitter & operator=(const TLinearFitter &tlf)
Assignment operator.
virtual void GetParameters(TVectorD &vpar)
Returns parameter values.
ClassImp(TMCParticle) void TMCParticle printf(": p=(%7.3f,%7.3f,%9.3f) ;", fPx, fPy, fPz)
Int_t GetNcols() const
Definition: TMatrixTBase.h:137
Double_t y[n]
Definition: legend1.C:17
double func(double *x, double *p)
Definition: stressTF1.cxx:213
virtual void Expand(Int_t newSize)
Expand or shrink the array to newSize elements.
Definition: TObjArray.cxx:369
The TH1 histogram class.
Definition: TH1.h:80
virtual void AssignData(Int_t npoints, Int_t xncols, Double_t *x, Double_t *y, Double_t *e=0)
This function is to use when you already have all the data in arrays and don't want to copy them into...
virtual void SetFormula(const char *formula)
Additive parts should be separated by "++".
Int_t GetEntries() const
Return the number of objects in array (i.e.
Definition: TObjArray.cxx:493
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:606
TVectorD fParams
virtual void SetDim(Int_t n)
set the number of dimensions
static const double eu
Definition: Vavilov.cxx:52
TMatrixDSym fDesignTemp
void ComputeTValues()
Computes parameters' t-values and significance.
int Nochisq
Definition: Foption.h:45
#define name(a, b)
Definition: linkTestLib0.cpp:5
Abstract Base Class for Fitting.
virtual Int_t GetXfirst() const
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:1044
Int_t GetLast() const
Return last bin on the axis i.e.
Definition: TAxis.cxx:440
TAxis * GetZaxis()
Definition: TH1.h:321
virtual const char * GetParName(Int_t ipar) const
Returns name of parameter #ipar.
const char * GetParName(Int_t ipar) const
Definition: TFormula.cxx:2288
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:449
Mother of all ROOT objects.
Definition: TObject.h:58
virtual Double_t GetErrorZ(Int_t bin) const
This function is called by Graph2DFitChisquare.
Definition: TGraph2D.cxx:1034
TMatrixDSym fParCovar
Container of bits.
Definition: TBits.h:33
Int_t GetNoElements() const
Definition: TVectorT.h:82
virtual Double_t GetErrorY(Int_t bin) const
This function is called by GraphFitChisquare.
Definition: TGraph.cxx:1380
virtual void SetPoint(Int_t i, Double_t x, Double_t y)
Set x and y values for point number i.
Definition: TGraph.cxx:2127
Element KOrdStat(Size n, const Element *a, Size k, Size *work=0)
Definition: TMath.h:1167
double f2(const double *x)
Int_t GetN() const
Definition: TGraph2D.h:127
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:567
1-Dim function class
Definition: TF1.h:149
virtual TObject * GetObjectFit() const
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:53
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:1185
TF1 * f1
Definition: legend1.C:11
Double_t * GetY() const
Definition: TGraph2D.h:129
virtual Double_t GetParError(Int_t ipar) const
Returns the error of parameter #ipar.
void Add(TObject *obj)
Definition: TObjArray.h:75
virtual Bool_t UpdateMatrix()
Update the design matrix after the formula has been changed.
Bool_t TestBitNumber(UInt_t bitnumber) const
Definition: TBits.h:223
virtual Int_t GetYfirst() const
Int_t fNpoints
temporary
virtual Bool_t IsInside(const Double_t *x) const
Return kTRUE is the point is inside the function range.
Definition: TF2.cxx:640
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
TMatrixDSym fDesignTemp3
Graphics object made of three arrays X, Y and Z with the same number of points each.
Definition: TGraph2D.h:50
virtual Bool_t IsInside(const Double_t *x) const
return kTRUE if the point is inside the function range
Definition: TF1.h:411
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8395
Long64_t LocMin(Long64_t n, const T *a)
Definition: TMath.h:695
const Bool_t kTRUE
Definition: Rtypes.h:91
virtual Double_t EvalPar(const Double_t *x, const Double_t *params=0)
Evaluate function with given coordinates and parameters.
Definition: TF1.cxx:1215
TObject * obj
TVectorD fAtbTemp2
temporary vector, used for num.stability
float value
Definition: math.cpp:443
virtual Int_t GetZfirst() const
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:198
TAxis * GetXaxis()
Definition: TH1.h:319
virtual Double_t * GetCovarianceMatrix() const
Returns covariance matrix.
T MinElement(Long64_t n, const T *a)
Definition: TMath.h:681
virtual void AddPoint(Double_t *x, Double_t y, Double_t e=1)
Adds 1 point to the fitter.
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:904