Logo ROOT   6.12/07
Reference Guide
combinedFit.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_fit
3 /// \notebook
4 /// Combined (simultaneous) fit of two histogram with separate functions
5 /// and some common parameters
6 ///
7 /// See http://root.cern.ch/phpBB3//viewtopic.php?f=3&t=11740#p50908
8 /// for a modified version working with Fumili or GSLMultiFit
9 ///
10 /// N.B. this macro must be compiled with ACliC
11 ///
12 /// \macro_image
13 /// \macro_output
14 /// \macro_code
15 ///
16 /// \author Lorenzo Moneta
17 
18 #include "Fit/Fitter.h"
19 #include "Fit/BinData.h"
20 #include "Fit/Chi2FCN.h"
21 #include "TH1.h"
22 #include "TList.h"
23 #include "Math/WrappedMultiTF1.h"
24 #include "HFitInterface.h"
25 #include "TCanvas.h"
26 #include "TStyle.h"
27 
28 
29 // definition of shared parameter
30 // background function
31 int iparB[2] = { 0, // exp amplitude in B histo
32  2 // exp common parameter
33 };
34 
35 // signal + background function
36 int iparSB[5] = { 1, // exp amplitude in S+B histo
37  2, // exp common parameter
38  3, // gaussian amplitude
39  4, // gaussian mean
40  5 // gaussian sigma
41 };
42 
43 // Create the GlobalCHi2 structure
44 
45 struct GlobalChi2 {
46  GlobalChi2( ROOT::Math::IMultiGenFunction & f1,
48  fChi2_1(&f1), fChi2_2(&f2) {}
49 
50  // parameter vector is first background (in common 1 and 2)
51  // and then is signal (only in 2)
52  double operator() (const double *par) const {
53  double p1[2];
54  for (int i = 0; i < 2; ++i) p1[i] = par[iparB[i] ];
55 
56  double p2[5];
57  for (int i = 0; i < 5; ++i) p2[i] = par[iparSB[i] ];
58 
59  return (*fChi2_1)(p1) + (*fChi2_2)(p2);
60  }
61 
62  const ROOT::Math::IMultiGenFunction * fChi2_1;
63  const ROOT::Math::IMultiGenFunction * fChi2_2;
64 };
65 
66 void combinedFit() {
67 
68  TH1D * hB = new TH1D("hB","histo B",100,0,100);
69  TH1D * hSB = new TH1D("hSB","histo S+B",100, 0,100);
70 
71  TF1 * fB = new TF1("fB","expo",0,100);
72  fB->SetParameters(1,-0.05);
73  hB->FillRandom("fB");
74 
75  TF1 * fS = new TF1("fS","gaus",0,100);
76  fS->SetParameters(1,30,5);
77 
78  hSB->FillRandom("fB",2000);
79  hSB->FillRandom("fS",1000);
80 
81  // perform now global fit
82 
83  TF1 * fSB = new TF1("fSB","expo + gaus(2)",0,100);
84 
85  ROOT::Math::WrappedMultiTF1 wfB(*fB,1);
86  ROOT::Math::WrappedMultiTF1 wfSB(*fSB,1);
87 
89  ROOT::Fit::DataRange rangeB;
90  // set the data range
91  rangeB.SetRange(10,90);
92  ROOT::Fit::BinData dataB(opt,rangeB);
93  ROOT::Fit::FillData(dataB, hB);
94 
95  ROOT::Fit::DataRange rangeSB;
96  rangeSB.SetRange(10,50);
97  ROOT::Fit::BinData dataSB(opt,rangeSB);
98  ROOT::Fit::FillData(dataSB, hSB);
99 
100  ROOT::Fit::Chi2Function chi2_B(dataB, wfB);
101  ROOT::Fit::Chi2Function chi2_SB(dataSB, wfSB);
102 
103  GlobalChi2 globalChi2(chi2_B, chi2_SB);
104 
105  ROOT::Fit::Fitter fitter;
106 
107  const int Npar = 6;
108  double par0[Npar] = { 5,5,-0.1,100, 30,10};
109 
110  // create before the parameter settings in order to fix or set range on them
111  fitter.Config().SetParamsSettings(6,par0);
112  // fix 5-th parameter
113  fitter.Config().ParSettings(4).Fix();
114  // set limits on the third and 4-th parameter
115  fitter.Config().ParSettings(2).SetLimits(-10,-1.E-4);
116  fitter.Config().ParSettings(3).SetLimits(0,10000);
117  fitter.Config().ParSettings(3).SetStepSize(5);
118 
119  fitter.Config().MinimizerOptions().SetPrintLevel(0);
120  fitter.Config().SetMinimizer("Minuit2","Migrad");
121 
122  // fit FCN function directly
123  // (specify optionally data size and flag to indicate that is a chi2 fit)
124  fitter.FitFCN(6,globalChi2,0,dataB.Size()+dataSB.Size(),true);
125  ROOT::Fit::FitResult result = fitter.Result();
126  result.Print(std::cout);
127 
128  TCanvas * c1 = new TCanvas("Simfit","Simultaneous fit of two histograms",
129  10,10,700,700);
130  c1->Divide(1,2);
131  c1->cd(1);
132  gStyle->SetOptFit(1111);
133 
134  fB->SetFitResult( result, iparB);
135  fB->SetRange(rangeB().first, rangeB().second);
136  fB->SetLineColor(kBlue);
137  hB->GetListOfFunctions()->Add(fB);
138  hB->Draw();
139 
140  c1->cd(2);
141  fSB->SetFitResult( result, iparSB);
142  fSB->SetRange(rangeSB().first, rangeSB().second);
143  fSB->SetLineColor(kRed);
144  hSB->GetListOfFunctions()->Add(fSB);
145  hSB->Draw();
146 
147 
148 }
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:628
void SetPrintLevel(int level)
set print level
return c1
Definition: legend1.C:41
Definition: Rtypes.h:59
R__EXTERN TStyle * gStyle
Definition: TStyle.h:402
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:688
virtual void SetRange(Double_t xmin, Double_t xmax)
Initialize the upper and lower bounds to draw the function.
Definition: TF1.cxx:3406
bool FitFCN(unsigned int npar, Function &fcn, const double *params=0, unsigned int dataSize=0, bool chi2fit=false)
Fit using the a generic FCN function as a C++ callable object implementing double () (const double *)...
Definition: Fitter.h:590
ROOT::Math::MinimizerOptions & MinimizerOptions()
access to the minimizer control parameter (non const method)
Definition: FitConfig.h:166
TRObject operator()(const T1 &t1) const
const ParameterSettings & ParSettings(unsigned int i) const
get the parameter settings for the i-th parameter (const method)
Definition: FitConfig.h:75
const FitResult & Result() const
get fit result
Definition: Fitter.h:365
Class to Wrap a ROOT Function class (like TF1) in a IParamMultiFunction interface of multi-dimensions...
static double p2(double t, double a, double b, double c)
void Fix()
fix the parameter
const FitConfig & Config() const
access to the fit configuration (const method)
Definition: Fitter.h:393
static constexpr double second
Documentation for the abstract class IBaseFunctionMultiDim.
Definition: IFunction.h:62
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:40
Chi2FCN class for binnned fits using the least square methods.
Definition: Chi2FCN.h:49
void SetMinimizer(const char *type, const char *algo=0)
set minimizer type
Definition: FitConfig.h:180
void SetStepSize(double err)
set the step size
void FillData(BinData &dv, const TH1 *hist, TF1 *func=0)
fill the data vector from a TH1.
virtual void FillRandom(const char *fname, Int_t ntimes=5000)
Fill histogram following distribution in function fname.
Definition: TH1.cxx:3414
DataOptions : simple structure holding the options on how the data are filled.
Definition: DataOptions.h:28
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2969
void SetRange(unsigned int icoord, double xmin, double xmax)
set a range [xmin,xmax] for the new coordinate icoord If more range exists for other coordinates...
Definition: DataRange.cxx:124
Fitter class, entry point for performing all type of fits.
Definition: Fitter.h:77
void SetOptFit(Int_t fit=1)
The type of information about fit parameters printed in the histogram statistics box can be selected ...
Definition: TStyle.cxx:1218
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition: BinData.h:53
static double p1(double t, double a, double b)
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:610
constexpr Double_t E()
Definition: TMath.h:74
The Canvas class.
Definition: TCanvas.h:31
class containg the result of the fit and all the related information (fitted parameter values...
Definition: FitResult.h:48
class describing the range in the coordinates it supports multiple range in a coordinate.
Definition: DataRange.h:34
void Print(std::ostream &os, bool covmat=false) const
print the result and optionaly covariance matrix and correlations
Definition: FitResult.cxx:446
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1153
virtual void Add(TObject *obj)
Definition: TList.h:87
void SetParamsSettings(unsigned int npar, const double *params, const double *vstep=0)
set the parameter settings from number of parameters and a vector of values and optionally step value...
Definition: FitConfig.cxx:136
1-Dim function class
Definition: TF1.h:211
virtual void SetFitResult(const ROOT::Fit::FitResult &result, const Int_t *indpar=0)
Set the result from the fit parameter values, errors, chi2, etc...
Definition: TF1.cxx:3242
TF1 * f1
Definition: legend1.C:11
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
Definition: THist.hxx:284
void SetLimits(double low, double up)
set a double side limit, if low == up the parameter is fixed if low > up the limits are removed The c...
Definition: Rtypes.h:59
Definition: first.py:1
TList * GetListOfFunctions() const
Definition: TH1.h:238