AliRoot Core  3dc7879 (3dc7879)
AliTMinuitToolkitTest.C
Go to the documentation of this file.
1 
4 
15 #include "AliTMinuitToolkit.h"
16 #include "TTree.h"
17 #include "TTreeStream.h"
18 #include "TF1.h"
19 #include "TVectorD.h"
20 #include "TString.h"
21 #include "TRandom.h"
22 #include "TH1.h"
23 #include "TCanvas.h"
24 #include "TLegend.h"
25 #include "TStyle.h"
26 #include "AliSysInfo.h"
27 #include "TGraphErrors.h"
28 #include "TStatToolkit.h"
29 #include "AliDrawStyle.h"
30 
31 Int_t fitEntries=400;
32 const Int_t kNDim=10;
33 const Int_t kLineColors[5]={1,2,4,3,6};
34 const Int_t kLineStyle[5]={1,7,10,1,2};
35 TF1 likeGausCachy("likeGausCachy", AliTMinuitToolkit::GaussCachyLogLike,-10,10,2);
36 TF1 likeAbs("likeAbs", "abs(x)",-10,10);
37 
38 TTree *inputTree=0;
39 void GenerateInput();
40 void Test1D(Int_t bootStrapIter=400);
41 void AliTMinuitToolkit_TestHistogram(Int_t nIter); // can be used
42 
46 void AliTMinuitToolkitTest(Int_t nIter, Int_t nPoints=1000){
48  AliDrawStyle::ApplyStyle("figTemplate");
49  fitEntries=nPoints;
51  GenerateInput();
52  Test1D(nIter);
53 }
54 
55 
58  //
59  // 0. Generate input data n dimensional array 1,1
60  TTreeSRedirector *pcstream = new TTreeSRedirector("AliTMinuiToolkitTestInput.root","recreate");
61  TVectorD xxx(kNDim);
62  for (Int_t i=0; i<fitEntries; i++){
63  for (Int_t iDim=0; iDim<kNDim; iDim++){
64  xxx[iDim]=2.*(gRandom->Rndm()-0.5);
65  }
66  Double_t noiseG=gRandom->Gaus();
67  Double_t noiseL=gRandom->Landau();
68  if (inputTree==NULL) {
69  (*pcstream)<<"data"<<
70  "x.="<<&xxx<<
71  "noiseG="<<noiseG<<
72  "noiseL="<<noiseL<<
73  "\n";
74  inputTree=((*pcstream)<<"data").GetTree();
75  }else{
76  inputTree->Fill();
77  }
78  }
79  delete pcstream;
80  TFile *f = TFile::Open("AliTMinuiToolkitTestInput.root");
81  inputTree= (TTree*)f->Get("data");
82 }
83 
84 
104  //
105  gStyle->SetOptStat(0);
106  TCanvas * canvas = new TCanvas("AliTMinuitToolkit_TestHistogram","AliTMinuitToolkit_TestHistogram",1200,1000);
107  canvas->Divide(3,3,0,0);
108  TMatrixD initParam(2, 4); // param,error,min,max
109  initParam(0, 0) = 20000;
110  initParam(0, 1) = 100;
111  initParam(0, 2) = 0;
112  initParam(0, 3) = 100000;
113  initParam(1, 0) = 1;
114  initParam(1, 1) = 1;
115  initParam(1, 2) = 0;
116  initParam(1, 3) = 10;
117  //
118  TTreeSRedirector *pcstream=new TTreeSRedirector("AliTMinuitToolkit_TestHistogram.root","recreate");
119  std::map<string,TF1*> fitMap;
120  TF1 *likeGausCachy = new TF1("likeGausCachy", AliTMinuitToolkit::GaussCachyLogLike, -10, 10, 2);
121  likeGausCachy->SetParameters(0.9, 1);
122  TF1 *likePseudoHuber = new TF1("likePseudoHuber", AliTMinuitToolkit::PseudoHuberLogLike, -10, 10, 2);
123  likePseudoHuber->SetParameter(0,3);
124 
125  for (Int_t iter=0; iter<nIter; iter++) {
126  Double_t slope=1+gRandom->Gaus(0,0.1);
127  if (iter <= 9) canvas->cd(iter+1);
128  // 0.1 provide some example histogram
129  TH1F *hist = new TH1F("test", "AliTMinuitToolkit: Test histogram fit (e^{-x}) with outliers", 50, 0, 4);
130  hist->SetMarkerStyle(25);
131  TRandom *random = new TRandom();
132  for (Int_t i = 0; i < 20000; i++) {
133  hist->Fill(random->Exp(slope));
134  }
135  for (Int_t iOutlier = 0; iOutlier < 10; iOutlier++) {
136  Double_t position = gRandom->Rndm() * 4;
137  Double_t value = gRandom->Rndm() * hist->GetEntries();
138  hist->Fill(position, value);
139  }
140  hist->SetMinimum(0);hist->SetMaximum(2500);
141  hist->Draw("error");
142  hist->Draw("same LHist");
143 
144  // 0.2 declare fit functions
145  TF1 *finput = new TF1("funFit1", "[0]*TMath::Exp(-[1]*x)", 0, 6);
146  TF1 *funFit1 = new TF1("funFit1", "[0]*TMath::Exp(-[1]*x)", 0, 6);
147  TF1 *funFit2 = new TF1("funFit2", "[0]*TMath::Exp(-[1]*x)", 0, 6);
148  TF1 *funFit3 = new TF1("funFit3", "[0]*TMath::Exp(-[1]*x)", 0, 6);
149  TVectorD oParam(2);
150  // 1.1) example fit without robust option
151  AliTMinuitToolkit *fitter = new AliTMinuitToolkit();
152  TF1 *aFormExp = new TF1("formExp", "[0]*TMath::Exp(-[1]*x)");
153  fitter->SetFitFunction(aFormExp, 0);
154  fitter->SetInitialParam(&initParam);
155  fitter->FitHistogram(hist);
156  funFit1->SetLineColor(1);
157  funFit1->SetParameters(fitter->GetParameters()->GetMatrixArray());
158  funFit1->Draw("same");
159  fitMap["chi2"]=(TF1*)fitter->GetFormula()->Clone();
160  // 1.2) Use "robust" custom user defined log likelihood function e.g gauss+background cachy
161  fitter->SetLogLikelihoodFunction(likeGausCachy);
162  fitter->SetInitialParam(&initParam);
163  fitter->FitHistogram(hist);
164  funFit2->SetLineColor(2);
165  funFit2->SetParameters(fitter->GetParameters()->GetMatrixArray());
166  funFit2->Draw("same");
167  // 1.3) Use predefined huber cost function
168  fitter->SetInitialParam(&initParam);
169  fitter->SetLogLikelihoodFunction(NULL);
170  fitter->EnableRobust(true);
171  fitter->FitHistogram(hist);
172  funFit3->SetLineColor(4);
173  funFit3->SetParameters(fitter->GetParameters()->GetMatrixArray());
174  funFit3->Draw("same");
175  fitMap["huber"]=(TF1*)fitter->GetFormula()->Clone();
176  //
177  // 2.) Test different fit strategies
178  fitter->SetLogLikelihoodFunction(likePseudoHuber);
179  AliTMinuitToolkit::SetPredefinedFitter("ExpFit", fitter);
180  AliTMinuitToolkit::Fit(hist, "ExpFit", "", NULL, "funOption(2,2,1)");
181  fitMap["pseudoHuber"]=(TF1*)fitter->GetFormula()->Clone();
182  AliTMinuitToolkit::Fit(hist, "ExpFit", "misac(10,50)", NULL, "funOption(2,2,1)");
183  fitMap["misacH(10,50)"]=(TF1*)fitter->GetFormula()->Clone();
184  AliTMinuitToolkit::Fit(hist, "ExpFit", "misac(10,100)", NULL, "funOption(4,2,2)");
185  fitMap["misacH(10,100)"]=(TF1*)fitter->GetFormula()->Clone();
186  AliTMinuitToolkit::Fit(hist, "ExpFit", "bootstrap50", NULL, "funOption(6,2,3)");
187  fitMap["bootstrapH50"]=(TF1*)fitter->GetFormula()->Clone();
188  //
189  fitter->SetLogLikelihoodFunction(likeGausCachy);
190  AliTMinuitToolkit::Fit(hist, "ExpFit", "misac(10,50)", NULL, "funOption(2,2,1)");
191  fitMap["misacGC(10,50)"]=(TF1*)fitter->GetFormula()->Clone();
192  AliTMinuitToolkit::Fit(hist, "ExpFit", "misac(10,100)", NULL, "funOption(4,2,2)");
193  fitMap["misacGC(10,100)"]=(TF1*)fitter->GetFormula()->Clone();
194  AliTMinuitToolkit::Fit(hist, "ExpFit", "bootstrap50", NULL, "funOption(6,2,3)");
195  fitMap["bootstrapGC50"]=(TF1*)fitter->GetFormula()->Clone();
196  //
197  //
198  // 2. Draw and save results
199  TLegend *legend = new TLegend(0.4, 0.7, 0.89, 0.89, "Test AliTMinuitToolkit - Exponential fit with outliers");
200  legend->SetBorderSize(0);
201  legend->AddEntry(hist, "Histogram");
202  legend->AddEntry(funFit1, "Default chi2 minimization");
203  legend->AddEntry(funFit2, "User defined likelihood (0.95*gaus+0.05*cachy)");
204  legend->AddEntry(funFit3, "Huber likelihood");
205  legend->Draw();
206  //
207  (*pcstream)<<"test"<<
208  "slope="<<slope<<
209  "his.="<<hist<<
210  "chi2.="<<fitMap["chi2"]<<
211  "huber.="<<fitMap["huber"]<<
212  "pseudoHuber.="<<fitMap["pseudoHuber"]<<
213  "misacH1050.="<<fitMap["misacH(10,50)"]<<
214  "misacH10100.="<<fitMap["misacH(10,100)"]<<
215  "bootstrapH50.="<<fitMap["bootstrapH50"]<<
216  "misacGC1050.="<<fitMap["misacGC(10,50)"]<<
217  "misacGC10100.="<<fitMap["misacGC(10,100)"]<<
218  "bootstrapGC50.="<<fitMap["bootstrapGC50"]<<
219 
220  "\n";
221  }
222  canvas->SaveAs("AliTMinuitToolkit_TestHistogram.png");
223  delete pcstream;
224 }
225 
226 
227 
228 
229 void Test1D(Int_t bootStrapIter){
230  //
231  // 1D fit example
232  // chi2, huber norm, abs(norm) and gaus+cachy log likelihood function
233  // Input:
234  // y=p[0]+p[1]*x; // p[0]=0; p[1]=2
235  //
236  TVectorD oParam(2); oParam[0]=0; oParam[1]=2;
237  TMatrixD initParam(2,4); // param,error,min,max
238  initParam(0,0)=1; initParam(0,1)=1; initParam(0,2)=0; initParam(0,3)=100000;
239  initParam(1,0)=1; initParam(1,1)=1; initParam(1,2)=0; initParam(1,3)=20;
240 
241  TF1 *fitFunctions[5];
242  TH1 *resHistograms[5];
243  TF1 formula1D("formula1D","[0]+[1]*x[0]",-1,1);
244  inputTree->SetAlias("test1D","2*x.fElements[0]");
245  inputTree->SetAlias("noise","(rndm<0.7)?noiseG:4*noiseL");
246  inputTree->SetAlias("X0","x.fElements[0]");
247 
248  TString selection="1";
249  AliTMinuitToolkit * tool1D = new AliTMinuitToolkit("AliTMinuitToolkitTest1D.root");
250  tool1D->SetVerbose(0x1);
251  tool1D->SetFitFunction(&formula1D,kTRUE);
252  tool1D->SetInitialParam(&initParam);
253  tool1D->FillFitter(inputTree,"test1D+noise:1/sqrt(12.+0)","X0", "", 0,fitEntries);
254  //
255  formula1D.SetParameters(oParam.GetMatrixArray());
256  fitFunctions[0]= (TF1*)formula1D.DrawClone("same");
257 
258  // 1.1) Standard fit
259  tool1D->EnableRobust(kFALSE);
260  tool1D->SetInitialParam(&initParam);
261  tool1D->Fit("");
262  inputTree->SetAlias("fitChi2Norm",tool1D->GetFitFunctionAsAlias().Data());
263  formula1D.SetParameters(tool1D->GetParameters()->GetMatrixArray());
264  fitFunctions[1]= (TF1*)formula1D.DrawClone("same");
265  tool1D->Bootstrap(bootStrapIter,"bootstrapChi2Norm");
266  tool1D->TwoFoldCrossValidation(bootStrapIter,"cvChi2Norm");
267  // 1.2) Huber cost function
268  tool1D->EnableRobust(kTRUE);
269  tool1D->SetInitialParam(&initParam);
270  tool1D->Fit("");
271  inputTree->SetAlias("fitHuberNorm",tool1D->GetFitFunctionAsAlias().Data());
272  formula1D.SetParameters(tool1D->GetParameters()->GetMatrixArray());
273  fitFunctions[2]= (TF1*)formula1D.DrawClone("same");
274  tool1D->Bootstrap(bootStrapIter,"bootstrapHuberNorm");
275  tool1D->TwoFoldCrossValidation(bootStrapIter,"cvHuberNorm");
276  // 1.3) User cost function
278  tool1D->Fit("");
279  inputTree->SetAlias("fitLikeAbs",tool1D->GetFitFunctionAsAlias().Data());
280  formula1D.SetParameters(tool1D->GetParameters()->GetMatrixArray());
281  fitFunctions[3]= (TF1*)formula1D.DrawClone("same");
282  tool1D->Bootstrap(bootStrapIter,"bootstrapLikeAbs");
283  tool1D->TwoFoldCrossValidation(bootStrapIter,"cvLikeAbs");
284  // 1.3) User cost function
285  likeGausCachy.SetParameters(0.8,1);
287  tool1D->Fit("");
288  inputTree->SetAlias("fitUserGausAndCachy",tool1D->GetFitFunctionAsAlias().Data());
289  formula1D.SetParameters(tool1D->GetParameters()->GetMatrixArray());
290  fitFunctions[4]= (TF1*)formula1D.DrawClone("same");
291  tool1D->Bootstrap(bootStrapIter,"bootstrapUserGausAndCachy");
292  tool1D->TwoFoldCrossValidation(bootStrapIter,"cvUserGausAndCachy");
293  delete tool1D;
294  //
295  // Draw Results
296  //
297  gStyle->SetOptStat(0);
298  gStyle->SetOptTitle(0);
299  TCanvas *canvasTest1D = new TCanvas("Test1D","Test1D",1400,1000);
300  canvasTest1D->Divide(1,3);
301  TLatex latex;
302  latex.SetTextSize(0.055);
303 
304  canvasTest1D->cd(1);
305  TGraphErrors *gr0 = TStatToolkit::MakeGraphErrors(inputTree,"test1D+noise:X0:1","",25,1,0.5,0,TMath::Min(fitEntries,100));
306  gr0->SetMaximum(30); gr0->SetMinimum(-30);
307  TGraphErrors *gr1 = TStatToolkit::MakeGraphErrors(inputTree,"test1D+noise:X0:1","",25,1,0.5,0,TMath::Min(fitEntries,100));
308  gr1->SetMaximum(5); gr1->SetMinimum(-5);
309 
310  gr0->Draw("ap");
311  for (Int_t iFit=0; iFit<5; iFit++){
312  fitFunctions[iFit]->SetLineColor(kLineColors[iFit]);
313  fitFunctions[iFit]->SetLineStyle(kLineStyle[iFit]);
314  fitFunctions[iFit]->SetLineWidth(3);
315  fitFunctions[iFit]->Draw("same");
316  }
317  latex.SetTextSize(0.07);
318  latex.DrawLatexNDC(0.11,0.8,"Input y=0+2*x+#epsilon ");
319  latex.DrawLatexNDC(0.11,0.7," #epsilon=(0.8 Gaus +0.2 Landau)");
320  latex.SetTextSize(0.055);
321  canvasTest1D->cd(2);
322 
323  gr1->Draw("ap");
324  for (Int_t iFit=0; iFit<5; iFit++){
325  fitFunctions[iFit]->SetLineColor(kLineColors[iFit]);
326  fitFunctions[iFit]->SetLineStyle(kLineStyle[iFit]);
327  fitFunctions[iFit]->SetLineWidth(3);
328  fitFunctions[iFit]->Draw("same");
329  }
330  TLegend * legend = new TLegend(0.11,0.7,0.6,0.89,"Unbinned1D fit for different cost function (ZOOM)");
331  legend->SetBorderSize(0); legend->SetNColumns(2);
332  legend->AddEntry(gr0, "Input data+noise","pl");
333  legend->AddEntry(fitFunctions[0],"MC input","l");
334  legend->AddEntry(fitFunctions[1],"Fit: Chi2 cost","l");
335  legend->AddEntry(fitFunctions[2],"Fit: Huber norm","l");
336  legend->AddEntry(fitFunctions[3],"Fit: LogLikeAbs","l");
337  legend->AddEntry(fitFunctions[4],"Fit: Log(80%Gaus+20%Cauchy)","l");
338  legend->Draw();
339  //
340  //for (Int_t iFit=0; iFit<3; iFit++) resHistograms[iFit]->Fit("gaus");
341  // TTre
342  TVirtualPad *pad =canvasTest1D->cd(3);
343  gStyle->SetOptTitle();
344  TFile *fReport = TFile::Open("AliTMinutiTookitTest1D.root");
345  TTree *tBootstrap = (TTree*)fReport->Get("bootstrap");
346  TTree *tCrossValidation = (TTree*)fReport->Get("crossValidation");
347  //
348  pad->Divide(2,1);
349  TH2 *h2D=0;
350  pad->cd(1);
351  tCrossValidation->SetAlias("dP0","(param0.fElements[0]-param1.fElements[0])");
352  tCrossValidation->SetAlias("dP1","(param0.fElements[1]-param1.fElements[1])");
353  tCrossValidation->Draw("dP0:reportName.GetName()>>hCrosValidDPar0(4,0,4,200,-1,1)","","colz");
354  h2D=(TH2*)(tCrossValidation->GetHistogram());
355  TGraph *grRMS0=TStatToolkit::MakeStat1D(h2D, 0, 1., 1, 21, 2);
356  grRMS0->Draw("p");
357  pad->cd(2);
358  tCrossValidation->Draw("dP1:reportName.GetName()>>hCrosValidDPar1(4,0,4,200,-1,1)","","colz");
359  h2D=(TH2*)(tCrossValidation->GetHistogram());
360  TGraph *grRMS1=TStatToolkit::MakeStat1D(h2D, 0, 1., 1, 21, 2);
361  grRMS1->Draw("p");
362  pad->cd(1);
363  latex.DrawLatexNDC(0.3,0.85,"Two fold cross validation. Param0");
364  pad->cd(2);
365  latex.DrawLatexNDC(0.3,0.85,"Two fold cross validation. Param1");
366  for (Int_t iFit=0; iFit<4; iFit++){
367  pad->cd(1);
368  latex.SetTextColor(kLineColors[iFit+1]);
369  latex.DrawLatexNDC(0.35,0.85-0.05*(iFit+1),TString::Format("%s: RMS %.2f ",h2D->GetXaxis()->GetBinLabel(1+iFit), grRMS0->GetY()[iFit]).Data());
370  pad->cd(2);
371  latex.SetTextColor(kLineColors[iFit+1]);
372  latex.DrawLatexNDC(0.35,0.85-0.05*(iFit+1),TString::Format("%s: RMS %.2f ",h2D->GetXaxis()->GetBinLabel(1+iFit), grRMS1->GetY()[iFit]).Data());
373  }
374  canvasTest1D->SaveAs("AliTMinuitToolkitTest.Test1D.png");
375  canvasTest1D->SaveAs("AliTMinuitToolkitTest.Test1D.pdf");
376 
377 }
378 
void FitHistogram(const TH1 *his, Option_t *option="")
Fit a one dimensional histogram.
static void ApplyStyle(const char *styleName)
void EnableRobust(Bool_t b)
void SetInitialParam(const TMatrixD *initialParam)
Set initial parameters matrix.
TFile * Open(const char *filename, Long64_t &nevents)
Int_t fitEntries
test of AliTMinuiToolkit class - log likelihood fits
TStyle * gStyle
void Fit(Option_t *option="")
Internal function that calls the fitter.
void Test1D(Int_t bootStrapIter=400)
TTreeSRedirector * pcstream
static Double_t PseudoHuberLogLike(const Double_t *x, const Double_t *p)
const Int_t kLineColors[5]
TGraphErrors * MakeGraphErrors(TTree *tree, const char *expr="Entry", const char *cut="1", Int_t mstyle=25, Int_t mcolor=1, Float_t msize=-1, Float_t offset=0.0, Int_t entries=10000000, Int_t firstEntry=0)
TF1 likeAbs("likeAbs","abs(x)",-10, 10)
const Int_t kLineStyle[5]
const TVectorD * GetParameters() const
void AliTMinuitToolkit_TestHistogram(Int_t nIter)
This test function shows the basic working principles of AliTMinuitToolkit.
TString GetFitFunctionAsAlias(Option_t *option="", TTree *tree=nullptr)
Format alias string for the for tree alias or for fit title.
static void SetDefaults()
SetDefault call RegisterDefaultLatexSymbols(), RegisterDefaultStyle(), RegisterDefaultMarkers();.
static void SetPredefinedFitter(const std::string &fitterName, AliTMinuitToolkit *fitter)
void SetVerbose(Int_t verbose)
void Bootstrap(ULong_t nIter, const char *reportName, Option_t *option=nullptr)
Bootstrap minimization.
TGraphErrors * MakeStat1D(TH2 *his, Int_t deltaBin, Double_t fraction, Int_t returnType, Int_t markerStyle, Int_t markerColor)
TF1 likeGausCachy("likeGausCachy", AliTMinuitToolkit::GaussCachyLogLike,-10, 10, 2)
TF1 * GetFormula() const
static Double_t GaussCachyLogLike(const Double_t *x, const Double_t *p)
void AliTMinuitToolkitTest(Int_t nIter, Int_t nPoints=1000)
TLatex latex
TF1 * f
Definition: interpolTest.C:21
The AliTMinuitToolkit serves as an easy to use interface for the TMinuit.
const Int_t kNDim
void GenerateInput()
Generate input data for the test.
void SetLogLikelihoodFunction(TF1 *logLikelihood)
TTree * inputTree
Long64_t FillFitter(TTree *inputTree, TString values, TString variables, TString selection, Int_t firstEntry, Int_t nEntries, Bool_t doReset=kTRUE)
Fill input matrices of the fitter.
void TwoFoldCrossValidation(UInt_t nIter, const char *reportName, Option_t *option=nullptr)
void SetFitFunction(TF1 *formula, Bool_t doReset)
Set formula of the fitted function and reset parameter. See doReset.
TTree * GetTree(Int_t ievent)