StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
TF1Fitter.cxx
1 #include <assert.h>
2 #include "TF1Fitter.h"
3 
4 //static Color_t gMyColors[] = {kRed,kBlue,kGreen,kMagenta,kCyan};
5 //_____________________________________________________________________________
6 
7 
8 TF1Fitter::TF1Fitter(const char *name,int nPars):TF1(name, static_cast<DummyFuncPtr_t>(nullptr), 0,0,nPars)
9 {
10 fMean=0; fSigm=0;fTH1 = 0;
11 }
12 //_____________________________________________________________________________
13 void TF1Fitter::Copy(TObject& f1) const
14 {
15  TF1Fitter &obj = (TF1Fitter&)f1;
16  obj.fTH1 = fTH1;
17  TF1::Copy(obj);
18 }
19 //_____________________________________________________________________________
20 double TF1Fitter::operator()(const double* x, const double* params)
21 {
22  if (!params) params = GetParameters();
23  return EvalPar(x,params);
24 }
25 //_____________________________________________________________________________
26 double TF1Fitter::EvalPar(const double *x,const double *par)
27 {
28  assert(0);
29  return 0;
30 }
31 //_____________________________________________________________________________
32 void TF1Fitter::SetHist(TH1 *hist)
33 {
34 
35  fTH1= hist;
36  if(!fTH1) return;
37  if (*GetName()==0) SetName(fTH1->GetName());
38  fMean = fTH1->GetMean();
39  fSigm = fTH1->GetRMS();
40  if (fXmin >= fXmax) {
41  int nBins = fTH1->GetNbinsX();
42  double Xmin = fTH1->GetBinLowEdge(1);
43  double Xmax = fTH1->GetBinLowEdge(nBins)+fTH1->GetBinWidth(nBins);
44  SetRange(Xmin,Xmax );
45  }
46 }
47 
48 //_____________________________________________________________________________
49 void TF1Fitter::Draw(const char *opt)
50 {
51  SetLineColor(kRed);
52  TF1::Draw(opt);
53 }
54 
55 
56 //_____________________________________________________________________________
57 double TF1GausFitter::EvalPar(const double *x,const double *par)
58 {
59  double p1 = 1./(par[1]*par[1]);
60  return par[2]*exp(-0.5*p1*(x[0]-par[0])*(x[0]-par[0]));
61 }
62 //_____________________________________________________________________________
63 void TF1GausFitter::Init()
64 {
65 static const double SQR2 = sqrt(2.);
66  SetParName(0,"Mean");
67  SetParName(1,"Sigm");
68  SetParName(2,"Norm");
69 
70 
71  double gral = fTH1->Integral();
72  double erf = 0.5*TMath::Erf((fXmax-fMean)/(SQR2*fSigm));
73  SetParameter(0,fMean);
74  SetParameter(1,fSigm);
75  SetParameter(2,gral/erf);
76  double xLow = GetParameters()[0]-9*GetParameters()[1];
77  double xUpp = GetParameters()[0]+9*GetParameters()[1];
78  SetParLimits(0,xLow,xUpp);
79  SetParLimits(1,GetParameters()[1]*0.1,GetParameters()[1]*10);
80  SetParLimits(2, 0,GetParameters()[2]*10);
81 }
82 
83 
84 //_____________________________________________________________________________
85 double TF1TwoGausFitter::EvalPar(const double *x,const double *par)
86 {
87  double sig1 = par[1];
88  double sig2 = par[1]+par[4];
89  double w1 = 1./(sig1*sig1);
90  double w2 = 1./(sig2*sig2);
91  double c1 = par[2]*(1-par[5]);
92  double c2 = par[2]*( par[5]);
93 
94  return c1*exp(-0.5*w1*(x[0]-par[0])*(x[0]-par[0]))
95  +c2*exp(-0.5*w2*(x[0]-par[3])*(x[0]-par[3]));
96 }
97 
98 //_____________________________________________________________________________
99 void TF1TwoGausFitter::Init()
100 {
101 static const double SQR2 = sqrt(2.);
102  SetParName(0,"Mean");
103  SetParName(1,"Sigm");
104  SetParName(2,"Norm");
105  SetParName(3,"MeanBak");
106  SetParName(4,"AddSigmBak");
107  SetParName(5,"ContriBak");
108 
109 
110  double gral = fTH1->Integral();
111  double erf = 0.5*TMath::Erf((fXmax-fMean)/(SQR2*fSigm));
112  SetParameter(0,fMean);
113  SetParameter(1,fSigm);
114  SetParameter(2,gral/erf);
115  double xLow = GetParameters()[0]-9*GetParameters()[1];
116  double xUpp = GetParameters()[0]+9*GetParameters()[1];
117  SetParLimits(0,xLow,xUpp);
118  SetParLimits(1,GetParameters()[1]*0.1,GetParameters()[1]*10);
119  SetParLimits(2, 0,GetParameters()[2]*10);
120 
121  SetParameter(3,fMean);
122  SetParameter(4,fSigm);
123  SetParameter(5,0.1);
124 
125  SetParLimits(3,xLow,xUpp);
126  SetParLimits(4,0,GetParameters()[1]*10);
127  SetParLimits(5,0,1);
128 
129 }
130 //_____________________________________________________________________________
131 double TF1GausPol2Fitter::EvalPar(const double *x,const double *par)
132 {
133  double sig1 = par[1];
134  double w1 = 1./(sig1*sig1);
135  double c1 = par[2];
136  double p0 = par[3];
137  double p1 = par[4];
138  double p2 = par[5];
139 
140  double bak = p0+x[0]*(p1 + x[0]*p2);
141  if (bak<0) bak = 0;
142  return c1*exp(-0.5*w1*(x[0]-par[0])*(x[0]-par[0])) + bak;
143 }
144 //_____________________________________________________________________________
145 void TF1GausPol2Fitter::Init()
146 {
147 static const double SQR2 = sqrt(2.);
148  SetParName(0,"Mean");
149  SetParName(1,"Sigm");
150  SetParName(2,"Norm");
151  SetParName(3,"b0");
152  SetParName(4,"b1");
153  SetParName(5,"b2");
154 
155 
156  double gral = fTH1->Integral();
157  double erf = 0.5*TMath::Erf((fXmax-fMean)/(SQR2*fSigm));
158  SetParameter(0,fMean);
159  SetParameter(1,fSigm);
160  SetParameter(2,gral/erf);
161  double xLow = GetParameters()[0]-9*GetParameters()[1];
162  double xUpp = GetParameters()[0]+9*GetParameters()[1];
163  SetParLimits(0,xLow,xUpp);
164  SetParLimits(1,GetParameters()[1]*0.1,GetParameters()[1]*10);
165 // SetParLimits(2, 0,GetParameters()[2]*10);
166 // SetParLimits(3, 0,999999);
167 // SetParLimits(5,-999999, 0);
168 
169  SetParameter(3,0);
170  SetParameter(4,0);
171  SetParameter(5,0);
172 }
173 
174 
175 
176 
177 
178 
179 
180 
181 
182 
183 
184 
185 
186 
187 
188 
189 
190 #if 0
191 //_____________________________________________________________________________
192 TwoGausFitter::TwoGausFitter(TH1 *th,const char *name):TH1Fitter(th,name)
193 {
194 }
195 //_____________________________________________________________________________
196 void TwoGausFitter::Prep()
197 {
198 static const double SQR2 = sqrt(2.);
199  if (mTF1) return;
200  TF1 *gaus1 = new TF1Gaus("Signal"); Add(gaus1);
201  TF1 *gaus2 = new TF1Gaus("BackGround"); Add(gaus2);
202 
203  double mean = mTH1->GetMean();
204  printf("mean=%g\n",mean);
205  double sigm = mTH1->GetRMS();
206  printf("sigm=%g\n",sigm);
207  double gral = mTH1->Integral();
208  printf("gral=%g\n",gral);
209  double erf = 0.5*TMath::Erf((mXmax-mean)/(SQR2*sigm));
210  printf("erf=%g\n",erf);
211  erf+= 0.5*TMath::Erf((mean-mXmin)/(SQR2*sigm));
212  printf("erf=%g\n",erf);
213  mPars[0] = mean;
214  mPars[1] = sigm;
215  mPars[2] = gral/erf;
216  GetTF1();
217  mTF1->SetParLimits(0,mPars[0]-3*sigm,mPars[1]+10*sigm);
218  mTF1->SetParLimits(1,mPars[1]*0.1,mPars[1]*10);
219  mTF1->SetParLimits(2, 0,mPars[2]*1000);
220 
221 
222  mPars[3] = mean;
223  mPars[4] = mPars[1]*0.1;
224  mPars[5] = mPars[2]*0.1;
225 
226  mTF1->SetParLimits(3,mPars[0]-10*sigm,mPars[1]+10*sigm);
227  mTF1->SetParLimits(4,mPars[4]*0.1,mPars[4]*10);
228  mTF1->SetParLimits(5,0. ,mPars[5]*1000);
229 
230 }
231 //_____________________________________________________________________________
232 double TwoGausFitter::Fit(const char *opt)
233 {
234  if (!mTF1) Prep();
235  return TH1Fitter::Fit(opt);
236 }
237 //_____________________________________________________________________________
238 TH1Fitter::TH1Fitter(TH1 *th,const char *name):TNamed(name,"")
239 {
240  memset(mBeg,0,mEnd-mBeg);
241  Set(th);
242 }
243 
244 //_____________________________________________________________________________
245 void TH1Fitter::Set(TH1 *th)
246 {
247  mTH1 = th;
248  if (!mTH1) return;
249  if (*GetName()==0) SetName(mTH1->GetName());
250  if (mXmin < mXmax) return;
251  int nBins = mTH1->GetNbinsX();
252  mXmin = mTH1->GetBinLowEdge(1);
253  mXmax = mTH1->GetBinLowEdge(nBins)+mTH1->GetBinWidth(nBins);
254 }
255 //_____________________________________________________________________________
256 void TH1Fitter::Add(TF1 *fcn)
257 {
258  assert(!mTF1);
259  mTF1s[mNTF1s++] = fcn;
260  mNPars += fcn->GetNpar();
261  assert(mNTF1s<10);
262 }
263 //_____________________________________________________________________________
264 TF1 *TH1Fitter::GetTF1()
265 {
266  if (!mTF1){//Create gloabal TF1
267 // create TF1 main
268  mTF1 = new TF1AuxFitter(mNPars, mNTF1s, mTF1s);
269  mTF1->SetName(GetName());
270  }
271  mTF1->SetRange(mXmin,mXmax);
272  return mTF1;
273 }
274 //_____________________________________________________________________________
275 double TH1Fitter::Fit(const char *opt)
276 {
277  GetTF1();
278  mTF1->SetParameters(mPars);
279  mTH1->Fit(mTF1,opt);
280  mTF1->GetParameters(mPars);
281  return mTF1->GetChisquare();
282 }
283 //_____________________________________________________________________________
284 void TH1Fitter::Draw(const char *)
285 {
286  double par[100],*ipar=par;
287  mTF1->GetParameters(par);
288 
289  for (int i = -1;i<mNTF1s;i++) {
290  TF1 *tf = mTF1s[i];
291  tf->SetLineColor(gMyColors[(i+1)%5]);
292  tf->SetParameters(ipar);
293  tf->Draw("same");
294  if (i>-1) ipar += tf->GetNpar();
295  }
296 }
297 #endif
Generated by   doxygen 1.8.5