source: trunk/WuerzburgSoft/Thomas/mphys/MFit.cc@ 10160

Last change on this file since 10160 was 1476, checked in by tbretz, 22 years ago
*** empty log message ***
File size: 3.2 KB
Line 
1#include "MFit.h"
2
3#include <iomanip.h>
4
5#include <TF1.h>
6#include <TH1.h>
7
8#include <TMinuit.h>
9
10ClassImp(MFit);
11
12TF1 *MFit::fgFunc=NULL;
13TH1 *MFit::fgHist=NULL;
14
15Int_t MFit::fgBinMin;
16Int_t MFit::fgBinMax;
17
18void MFit::Func(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t flag)
19{
20 Double_t chisq = 0;
21
22 for (int i=fgBinMin+1; i<fgBinMax; i++)
23 {
24 Double_t x = fgHist->GetBinCenter(i);
25
26 Double_t y = fgFunc->EvalPar(&x, par);
27 Double_t d = (y - fgHist->GetBinContent(i))/fgHist->GetBinError(i);
28
29 chisq += d*d;
30 }
31
32 f = chisq / fgFunc->GetNDF();
33 fgFunc->SetChisquare(f);
34}
35
36void MFit::FuncLog(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t flag)
37{
38 Double_t chisq = 0;
39
40 for (int i=fgBinMin+1; i<fgBinMax; i++)
41 {
42 Double_t x1 = fgHist->GetBinLowEdge(i+1);
43 Double_t x0 = fgHist->GetBinLowEdge(i);
44
45 Double_t x = sqrt(x1*x0);
46
47 Double_t y = fgFunc->EvalPar(&x, par);
48 Double_t d = (y - fgHist->GetBinContent(i))/fgHist->GetBinError(i);
49
50 chisq += d*d;
51 }
52
53 f = chisq / fgFunc->GetNDF();
54 fgFunc->SetChisquare(f);
55}
56
57void MFit::DeleteF()
58{
59 if (fFunc)
60 if (fFunc->TestBit(kIsOwner))
61 delete fFunc;
62
63 if (fMinuit)
64 delete fMinuit;
65}
66
67void MFit::DeleteH()
68{
69 if (fHist)
70 if (fHist->TestBit(kIsOwner))
71 delete fHist;
72}
73
74void MFit::SetFunc(TF1 *f)
75{
76 DeleteF();
77
78 fFunc = f;
79
80 fMinuit = new TMinuit(fFunc->GetNpar());
81 fMinuit->SetPrintLevel(-1);
82}
83
84void MFit::SetFunc(TString formula, Axis_t min, Axis_t max)
85{
86 TF1 *f = new TF1("", formula, min, max);
87 f->SetBit(kIsOwner);
88
89 SetFunc(f);
90}
91
92void MFit::SetRange(Axis_t min, Axis_t max)
93{
94 if (fFunc)
95 fFunc->SetRange(min, max);
96}
97
98void MFit::SetHist(TH1 *h, Bool_t candelete)
99{
100 DeleteH();
101
102 fHist = h;
103
104 if (candelete)
105 fHist->SetBit(kIsOwner);
106}
107
108void MFit::Fit(Bool_t log)
109{
110 fgBinMin = fHist->FindBin(fFunc->GetXmin()); // excluded
111 fgBinMax = fHist->FindBin(fFunc->GetXmax()); // excluded
112
113 fFunc->SetNDF(fgBinMax-fgBinMin-1);
114
115 fgHist = fHist;
116 fgFunc = fFunc;
117
118 fMinuit->SetFCN(log ? FuncLog : Func);
119 fMinuit->Migrad();
120
121 for (int i=0; i<fFunc->GetNpar(); i++)
122 {
123 Double_t par, err;
124 fMinuit->GetParameter(i, par, err);
125 fFunc->SetParameter(i, par);
126 fFunc->SetParError(i, err);
127 }
128}
129
130void MFit::Print(Option_t *opt) const
131{
132 for (int i=0; i<fFunc->GetNpar(); i++)
133 cout << "Par["<<i<<"]: " << GetParameter(i) << " +- " << GetParError(i) << endl;
134 cout << "ChiSq: " << fFunc->GetChisquare() << endl;
135 cout << "NDF: " << fFunc->GetNDF() << endl;
136}
137
138void MFit::SetParameter(Int_t n, TString name, Axis_t start, Axis_t min, Axis_t max)
139{
140 if (n>=fFunc->GetNpar())
141 return;
142
143 fMinuit->DefineParameter(n, name, start, min/10, min, max);
144}
145
146Double_t MFit::GetParameter(Int_t n) const
147{
148 return fFunc->GetParameter(n);
149}
150
151Double_t MFit::GetParError(Int_t n) const
152{
153 return fFunc->GetParError(n);
154}
155
156Double_t MFit::operator[](Int_t n) const
157{
158 return fFunc->GetParameter(n);
159}
160
161TF1 *MFit::DrawCopy(Option_t *o) const
162{
163 return fFunc->DrawCopy(o);
164}
165
Note: See TracBrowser for help on using the repository browser.