Logo 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 /*************************************************************************
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
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 
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());
985  ((TF1*)fInputFunction)->SetNumberFitPoints(fNpoints);
986  }
987  }
988 
989  //if parameters were fixed, change the design matrix back as it was before fixing
990  j = 0;
991  if (fNfixed>0){
992  for (i=0; i<fNfunctions; i++){
993  if (fFixedParams[i]){
994  for (ii=0; ii<i; ii++){
995  fDesign(ii, i) = fDesignTemp(ii, j);
996  fAtb(ii) = fDesignTemp(ii, fNfixed);
997  }
998  for (ii=i; ii<fNfunctions; ii++){
999  fDesign(i, ii) = fDesignTemp(ii, j);
1000  fAtb(ii) = fDesignTemp(ii, fNfixed);
1001  }
1002  j++;
1003  }
1004  }
1005  }
1006  return 0;
1007 }
1008 
1009 ////////////////////////////////////////////////////////////////////////////////
1010 ///Fixes paramter #ipar at its current value.
1011 
1013 {
1014  if (fParams.NonZeros()<1){
1015  Error("FixParameter", "no value available to fix the parameter");
1016  return;
1017  }
1018  if (ipar>fNfunctions || ipar<0){
1019  Error("FixParameter", "illegal parameter value");
1020  return;
1021  }
1022  if (fNfixed==fNfunctions) {
1023  Error("FixParameter", "no free parameters left");
1024  return;
1025  }
1026  if (!fFixedParams)
1028  fFixedParams[ipar] = 1;
1029  fNfixed++;
1030 }
1031 
1032 ////////////////////////////////////////////////////////////////////////////////
1033 ///Fixes parameter #ipar at value parvalue.
1034 
1036 {
1037  if (ipar>fNfunctions || ipar<0){
1038  Error("FixParameter", "illegal parameter value");
1039  return;
1040  }
1041  if (fNfixed==fNfunctions) {
1042  Error("FixParameter", "no free parameters left");
1043  return;
1044  }
1045  if(!fFixedParams)
1047  fFixedParams[ipar] = 1;
1050  fParams(ipar) = parvalue;
1051  fNfixed++;
1052 }
1053 
1054 ////////////////////////////////////////////////////////////////////////////////
1055 ///Releases parameter #ipar.
1056 
1058 {
1059  if (ipar>fNfunctions || ipar<0){
1060  Error("ReleaseParameter", "illegal parameter value");
1061  return;
1062  }
1063  if (!fFixedParams[ipar]){
1064  Warning("ReleaseParameter","This parameter is not fixed\n");
1065  return;
1066  } else {
1067  fFixedParams[ipar] = 0;
1068  fNfixed--;
1069  }
1070 }
1071 
1072 ////////////////////////////////////////////////////////////////////////////////
1073 ///Get the Atb vector - a vector, used for internal computations
1074 
1076 {
1077  if (v.GetNoElements()!=fAtb.GetNoElements())
1078  v.ResizeTo(fAtb.GetNoElements());
1079  v = fAtb;
1080 }
1081 
1082 ////////////////////////////////////////////////////////////////////////////////
1083 /// Get the Chisquare.
1084 
1086 {
1087  if (fChisquare > 1e-16)
1088  return fChisquare;
1089  else {
1090  Chisquare();
1091  return fChisquare;
1092  }
1093 }
1094 
1095 ////////////////////////////////////////////////////////////////////////////////
1096 ///Computes point-by-point confidence intervals for the fitted function
1097 ///Parameters:
1098 ///n - number of points
1099 ///ndim - dimensions of points
1100 ///x - points, at which to compute the intervals, for ndim > 1
1101 /// should be in order: (x0,y0, x1, y1, ... xn, yn)
1102 ///ci - computed intervals are returned in this array
1103 ///cl - confidence level, default=0.95
1104 ///
1105 ///NOTE, that this method can only be used when the fitting function inherits from a TF1,
1106 ///so it's not possible when the fitting function was set as a string or as a pure TFormula
1107 
1109 {
1110  if (fInputFunction){
1111  Double_t *grad = new Double_t[fNfunctions];
1112  Double_t *sum_vector = new Double_t[fNfunctions];
1113  Double_t c=0;
1115  Double_t t = TMath::StudentQuantile(0.5 + cl/2, df);
1116  Double_t chidf = TMath::Sqrt(fChisquare/df);
1117 
1118  for (Int_t ipoint=0; ipoint<n; ipoint++){
1119  c=0;
1120  ((TF1*)(fInputFunction))->GradientPar(x+ndim*ipoint, grad);
1121  //multiply the covariance matrix by gradient
1122  for (Int_t irow=0; irow<fNfunctions; irow++){
1123  sum_vector[irow]=0;
1124  for (Int_t icol=0; icol<fNfunctions; icol++){
1125  sum_vector[irow]+=fParCovar(irow,icol)*grad[icol];
1126  }
1127  }
1128  for (Int_t i=0; i<fNfunctions; i++)
1129  c+=grad[i]*sum_vector[i];
1130  c=TMath::Sqrt(c);
1131  ci[ipoint]=c*t*chidf;
1132  }
1133 
1134  delete [] grad;
1135  delete [] sum_vector;
1136  }
1137 }
1138 
1139 ////////////////////////////////////////////////////////////////////////////////
1140 ///Computes confidence intervals at level cl. Default is 0.95
1141 ///The TObject parameter can be a TGraphErrors, a TGraph2DErrors or a TH123.
1142 ///For Graphs, confidence intervals are computed for each point,
1143 ///the value of the graph at that point is set to the function value at that
1144 ///point, and the graph y-errors (or z-errors) are set to the value of
1145 ///the confidence interval at that point
1146 ///For Histograms, confidence intervals are computed for each bin center
1147 ///The bin content of this bin is then set to the function value at the bin
1148 ///center, and the bin error is set to the confidence interval value.
1149 ///Allowed combinations:
1150 ///Fitted object Passed object
1151 ///TGraph TGraphErrors, TH1
1152 ///TGraphErrors, AsymmErrors TGraphErrors, TH1
1153 ///TH1 TGraphErrors, TH1
1154 ///TGraph2D TGraph2DErrors, TH2
1155 ///TGraph2DErrors TGraph2DErrors, TH2
1156 ///TH2 TGraph2DErrors, TH2
1157 ///TH3 TH3
1158 
1160 {
1161  if (!fInputFunction) {
1162  Error("GetConfidenceIntervals", "The case of fitting not with a TFormula is not yet implemented");
1163  return;
1164  }
1165 
1166  //TGraph//////////////////
1167 
1168  if (obj->InheritsFrom(TGraph::Class())) {
1169  TGraph *gr = (TGraph*)obj;
1170  if (!gr->GetEY()){
1171  Error("GetConfidenceIntervals", "A TGraphErrors should be passed instead of a graph");
1172  return;
1173  }
1175  Error("GetConfidenceIntervals", "A TGraph2DErrors should be passed instead of a graph");
1176  return;
1177  }
1179  if (((TH1*)(fObjectFit))->GetDimension()>1){
1180  Error("GetConfidenceIntervals", "A TGraph2DErrors or a TH23 should be passed instead of a graph");
1181  return;
1182  }
1183  }
1184 
1185  GetConfidenceIntervals(gr->GetN(),1,gr->GetX(), gr->GetEY(), cl);
1186  for (Int_t i=0; i<gr->GetN(); i++)
1187  gr->SetPoint(i, gr->GetX()[i], fInputFunction->Eval(gr->GetX()[i]));
1188  }
1189 
1190  //TGraph2D///////////////
1191  else if (obj->InheritsFrom(TGraph2D::Class())) {
1192  TGraph2D *gr2 = (TGraph2D*)obj;
1193  if (!gr2->GetEZ()){
1194  Error("GetConfidenceIntervals", "A TGraph2DErrors should be passed instead of a TGraph2D");
1195  return;
1196  }
1198  Error("GetConfidenceIntervals", "A TGraphErrors should be passed instead of a TGraph2D");
1199  return;
1200  }
1202  if (((TH1*)(fObjectFit))->GetDimension()==1){
1203  Error("GetConfidenceIntervals", "A TGraphErrors or a TH1 should be passed instead of a graph");
1204  return;
1205  }
1206  }
1207  Double_t xy[2];
1208  Int_t np = gr2->GetN();
1209  Double_t *grad = new Double_t[fNfunctions];
1210  Double_t *sum_vector = new Double_t[fNfunctions];
1211  Double_t *x = gr2->GetX();
1212  Double_t *y = gr2->GetY();
1213  Double_t t = TMath::StudentQuantile(0.5 + cl/2, ((TF1*)(fInputFunction))->GetNDF());
1214  Double_t chidf = TMath::Sqrt(fChisquare/((TF1*)(fInputFunction))->GetNDF());
1215  Double_t c = 0;
1216  for (Int_t ipoint=0; ipoint<np; ipoint++){
1217  c=0;
1218  xy[0]=x[ipoint];
1219  xy[1]=y[ipoint];
1220  ((TF1*)(fInputFunction))->GradientPar(xy, grad);
1221  for (Int_t irow=0; irow<fNfunctions; irow++){
1222  sum_vector[irow]=0;
1223  for (Int_t icol=0; icol<fNfunctions; icol++)
1224  sum_vector[irow]+=fParCovar(irow, icol)*grad[icol];
1225  }
1226  for (Int_t i=0; i<fNfunctions; i++)
1227  c+=grad[i]*sum_vector[i];
1228  c=TMath::Sqrt(c);
1229  gr2->SetPoint(ipoint, xy[0], xy[1], fInputFunction->EvalPar(xy));
1230  gr2->GetEZ()[ipoint]=c*t*chidf;
1231  }
1232  delete [] grad;
1233  delete [] sum_vector;
1234  }
1235 
1236  //TH1////////////////////////
1237  else if (obj->InheritsFrom(TH1::Class())) {
1239  if (((TH1*)obj)->GetDimension()>1){
1240  Error("GetConfidenceIntervals", "Fitted graph and passed histogram have different number of dimensions");
1241  return;
1242  }
1243  }
1245  if (((TH1*)obj)->GetDimension()!=2){
1246  Error("GetConfidenceIntervals", "Fitted graph and passed histogram have different number of dimensions");
1247  return;
1248  }
1249  }
1251  if (((TH1*)(fObjectFit))->GetDimension()!=((TH1*)(obj))->GetDimension()){
1252  Error("GetConfidenceIntervals", "Fitted and passed histograms have different number of dimensions");
1253  return;
1254  }
1255  }
1256 
1257 
1258  TH1 *hfit = (TH1*)obj;
1259  Double_t *grad = new Double_t[fNfunctions];
1260  Double_t *sum_vector = new Double_t[fNfunctions];
1261  Double_t x[3];
1262  Int_t hxfirst = hfit->GetXaxis()->GetFirst();
1263  Int_t hxlast = hfit->GetXaxis()->GetLast();
1264  Int_t hyfirst = hfit->GetYaxis()->GetFirst();
1265  Int_t hylast = hfit->GetYaxis()->GetLast();
1266  Int_t hzfirst = hfit->GetZaxis()->GetFirst();
1267  Int_t hzlast = hfit->GetZaxis()->GetLast();
1268 
1269  TAxis *xaxis = hfit->GetXaxis();
1270  TAxis *yaxis = hfit->GetYaxis();
1271  TAxis *zaxis = hfit->GetZaxis();
1272  Double_t t = TMath::StudentQuantile(0.5 + cl/2, ((TF1*)(fInputFunction))->GetNDF());
1273  Double_t chidf = TMath::Sqrt(fChisquare/((TF1*)(fInputFunction))->GetNDF());
1274  Double_t c=0;
1275  for (Int_t binz=hzfirst; binz<=hzlast; binz++){
1276  x[2]=zaxis->GetBinCenter(binz);
1277  for (Int_t biny=hyfirst; biny<=hylast; biny++) {
1278  x[1]=yaxis->GetBinCenter(biny);
1279  for (Int_t binx=hxfirst; binx<=hxlast; binx++) {
1280  x[0]=xaxis->GetBinCenter(binx);
1281  ((TF1*)(fInputFunction))->GradientPar(x, grad);
1282  c=0;
1283  for (Int_t irow=0; irow<fNfunctions; irow++){
1284  sum_vector[irow]=0;
1285  for (Int_t icol=0; icol<fNfunctions; icol++)
1286  sum_vector[irow]+=fParCovar(irow, icol)*grad[icol];
1287  }
1288  for (Int_t i=0; i<fNfunctions; i++)
1289  c+=grad[i]*sum_vector[i];
1290  c=TMath::Sqrt(c);
1291  hfit->SetBinContent(binx, biny, binz, fInputFunction->EvalPar(x));
1292  hfit->SetBinError(binx, biny, binz, c*t*chidf);
1293  }
1294  }
1295  }
1296  delete [] grad;
1297  delete [] sum_vector;
1298  }
1299  else {
1300  Error("GetConfidenceIntervals", "This object type is not supported");
1301  return;
1302  }
1303 }
1304 
1305 ////////////////////////////////////////////////////////////////////////////////
1306 ///Returns covariance matrix
1307 
1309 {
1310  Double_t *p = const_cast<Double_t*>(fParCovar.GetMatrixArray());
1311  return p;
1312 }
1313 
1314 ////////////////////////////////////////////////////////////////////////////////
1315 ///Returns covariance matrix
1316 
1318 {
1319  if (matr.GetNrows()!=fNfunctions || matr.GetNcols()!=fNfunctions){
1321  }
1322  matr = fParCovar;
1323 }
1324 
1325 ////////////////////////////////////////////////////////////////////////////////
1326 ///Returns the internal design matrix
1327 
1329 {
1330  if (matr.GetNrows()!=fNfunctions || matr.GetNcols()!=fNfunctions){
1332  }
1333  matr = fDesign;
1334 }
1335 
1336 ////////////////////////////////////////////////////////////////////////////////
1337 ///Returns parameter errors
1338 
1340 {
1341  if (vpar.GetNoElements()!=fNfunctions) {
1342  vpar.ResizeTo(fNfunctions);
1343  }
1344  for (Int_t i=0; i<fNfunctions; i++)
1345  vpar(i) = TMath::Sqrt(fParCovar(i, i));
1346 
1347 }
1348 
1349 ////////////////////////////////////////////////////////////////////////////////
1350 ///Returns parameter values
1351 
1353 {
1354  if (vpar.GetNoElements()!=fNfunctions) {
1355  vpar.ResizeTo(fNfunctions);
1356  }
1357  vpar=fParams;
1358 }
1359 
1360 ////////////////////////////////////////////////////////////////////////////////
1361 ///Returns the value and the name of the parameter #ipar
1362 ///NB: In the calling function the argument name must be set large enough
1363 
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 = 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();
1608  fFunctions.Add(f);
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++) {
1724  AddToDesign(TMatrixDRow(fX, i).GetPtr(), fY(i), fE(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;
1807  AddPoint(&x[i], y[i], e);
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;
1874  AddPoint(x, z, e);
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;
1943  AddPoint(&gx[i], gy[i], e);
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;
1989  // Double_t dersum[100], grad[100];
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  }
2032  AddPoint(x, cu, eu);
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 {
2078  if (R__b.IsReading()) {
2079  Int_t old_matr_size = fNfunctions;
2081  if (old_matr_size < fNfunctions) {
2084 
2087 
2090  }
2091  } else {
2092  if (fAtb.NonZeros()==0) AddTempMatrices();
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++){
2334  AddToDesign(TMatrixDRow(fX, index[i]).GetPtr(), fY(index[i]), fE(index[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
2357  AddToDesign(TMatrixDRow(fX, index[nindex-1]).GetPtr(), fY(index[nindex-1]), fE(index[nindex-1]));
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++)
2446  AddToDesign(TMatrixDRow(fX, index[i]).GetPtr(), fY(index[i]), fE(index[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
Bool_t IsReading() const
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
virtual void AddTempMatrices()
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
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: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
void Add(TObject *obj)
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()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:316
Double_t * GetX() const
Definition: TGraph2D.h:119
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8485
T MinElement(Long64_t n, const T *a)
Return minimum of array a of length n.
Definition: TMath.h:942
const char * Data() const
Definition: TString.h:364
virtual void AddPoint(Double_t *x, Double_t y, Double_t e=1)
Adds 1 point to the fitter.