]>
Commit | Line | Data |
---|---|---|
0d73200d | 1 | // |
2 | // | |
3 | // | |
4 | ||
5 | #include <TF1.h> | |
6 | #include <TH1.h> | |
7 | #include <TMath.h> | |
8 | ||
9 | #include "AliLog.h" | |
10 | #include <AliRsnFitResult.h> | |
11 | ||
12 | ClassImp(AliRsnFitResult) | |
13 | ||
14 | //__________________________________________________________________________________________________ | |
15 | AliRsnFitResult::AliRsnFitResult(const char *name, ESignalType sType, EBackgroundType bType) : | |
16 | TNamed(name, ""), | |
17 | fSignalType(sType), | |
18 | fBackgroundType(bType) | |
19 | { | |
20 | // | |
21 | // Default constructor | |
22 | // | |
23 | ||
24 | InitFunctions(); | |
25 | } | |
26 | ||
27 | //__________________________________________________________________________________________________ | |
28 | Double_t AliRsnFitResult::NormSquareRoot(Double_t *x, Double_t *par) | |
29 | { | |
30 | // | |
31 | // Computes a square-root like function normalized | |
32 | // to the integral in the background range specified | |
33 | // for this object. | |
34 | // This is obtained by dividing the function value by | |
35 | // its integral analytically computed in that range. | |
36 | // Returns 0 in all cases where a floating point exception | |
37 | // could be raised by computing square root of a negative number. | |
38 | // | |
39 | ||
40 | if (x[0] < par[1]) return 0.0; | |
41 | ||
42 | Double_t fcnVal = TMath::Sqrt(x[0] - par[1]); | |
43 | ||
44 | Double_t fcnMax = fPeakRange[1] - par[1]; | |
45 | Double_t fcnMin = fPeakRange[0] - par[1]; | |
46 | Double_t fcnInt = (2.0 / 3.0) * ((fcnMax*TMath::Sqrt(fcnMax) - fcnMin*TMath::Sqrt(fcnMin))); | |
47 | ||
48 | return par[0] * fcnVal / fcnInt; | |
49 | } | |
50 | ||
51 | //__________________________________________________________________________________________________ | |
52 | Double_t AliRsnFitResult::NormLinear(Double_t *x, Double_t *par) | |
53 | { | |
54 | // | |
55 | // Computes a 1st order polynomial function normalized | |
56 | // to the integral in the background range specified | |
57 | // for this object. | |
58 | // This is obtained by dividing the function value by | |
59 | // its integral analytically computed in that range. | |
60 | // Returns 0 in all cases where a floating point exception | |
61 | // could be raised by computing square root of a negative number. | |
62 | // | |
63 | ||
64 | Double_t fcnVal = 1.0 + x[0] * par[1]; | |
65 | Double_t fcnInt = (fPeakRange[1] - fPeakRange[0]) + 0.5*par[1]*(fPeakRange[1]*fPeakRange[1] - fPeakRange[0]*fPeakRange[0]); | |
66 | ||
67 | return par[0] * fcnVal / fcnInt; | |
68 | } | |
69 | ||
70 | //__________________________________________________________________________________________________ | |
71 | Double_t AliRsnFitResult::NormPoly2(Double_t *x, Double_t *par) | |
72 | { | |
73 | // | |
74 | // Computes a 2nd order polynomial function normalized | |
75 | // to the integral in the background range specified | |
76 | // for this object. | |
77 | // This is obtained by dividing the function value by | |
78 | // its integral analytically computed in that range. | |
79 | // Returns 0 in all cases where a floating point exception | |
80 | // could be raised by computing square root of a negative number. | |
81 | // | |
82 | ||
83 | Double_t fcnVal = 1.0 + x[0] * par[1] + x[0]*x[0] * par[2]; | |
84 | Double_t fcnInt = (fPeakRange[1] - fPeakRange[0]) + par[1]*(fPeakRange[1]*fPeakRange[1] - fPeakRange[0]*fPeakRange[0])/2.0 + par[2]*(fPeakRange[1]*fPeakRange[1]*fPeakRange[1] - fPeakRange[0]*fPeakRange[0]*fPeakRange[0])/3.0; | |
85 | ||
86 | return par[0] * fcnVal / fcnInt; | |
87 | } | |
88 | ||
89 | //__________________________________________________________________________________________________ | |
90 | Double_t AliRsnFitResult::NormBreitWigner(Double_t *x, Double_t *par) | |
91 | { | |
92 | // | |
93 | // Computes a Breit-Wigner function normalized | |
94 | // to the integral in the background range specified | |
95 | // for this object. | |
96 | // This is obtained by dividing the function value by | |
97 | // its integral analytically computed in that range. | |
98 | // Returns 0 in all cases where a floating point exception | |
99 | // could be raised by computing square root of a negative number. | |
100 | // | |
101 | ||
102 | Double_t fcnVal = 1.0 / ((x[0] - par[1])*(x[0] - par[1]) + 0.25*par[2]*par[2]); | |
103 | Double_t fcnInt = 2.0 / par[2] * (TMath::ATan((fPeakRange[1] - par[1]) / (0.5 * par[2])) - TMath::ATan((fPeakRange[0] - par[1]) / (0.5 * par[2]))); | |
104 | ||
105 | return par[0] * fcnVal / fcnInt; | |
106 | } | |
107 | ||
108 | //__________________________________________________________________________________________________ | |
109 | Double_t AliRsnFitResult::NormGaus(Double_t *x, Double_t *par) | |
110 | { | |
111 | // | |
112 | // Computes a Gaussian function normalized | |
113 | // to the integral in the background range specified | |
114 | // for this object. | |
115 | // This is obtained by dividing the function value by | |
116 | // its integral analytically computed in that range. | |
117 | // Returns 0 in all cases where a floating point exception | |
118 | // could be raised by computing square root of a negative number. | |
119 | // | |
120 | ||
121 | Double_t fcnVal = TMath::Gaus(x[0], par[1], par[2], kTRUE); | |
122 | Double_t fcnIntLeft = 0.5 * TMath::Erf(((par[1] - fPeakRange[0]) / par[2]) / TMath::Sqrt(2.0)); | |
123 | Double_t fcnIntRight = 0.5 * TMath::Erf(((fPeakRange[1] - par[1]) / par[2]) / TMath::Sqrt(2.0)); | |
124 | Double_t fcnInt = fcnIntLeft + fcnIntRight; | |
125 | ||
126 | return par[0] * fcnVal / fcnInt; | |
127 | } | |
128 | ||
129 | //__________________________________________________________________________________________________ | |
130 | Double_t AliRsnFitResult::Signal(Double_t *x, Double_t *par) | |
131 | { | |
132 | // | |
133 | // Computes the signal function according to chosen option | |
134 | // | |
135 | ||
136 | switch (fSignalType) | |
137 | { | |
138 | case kBreitWigner: | |
139 | return NormBreitWigner(x, par); | |
140 | case kGaus: | |
141 | return NormGaus(x, par); | |
142 | default: | |
143 | AliError("Invalid signal function"); | |
144 | return 0.0; | |
145 | } | |
146 | } | |
147 | ||
148 | //__________________________________________________________________________________________________ | |
149 | Double_t AliRsnFitResult::Background(Double_t *x, Double_t *par) | |
150 | { | |
151 | // | |
152 | // Computes the background function according to chosen option | |
153 | // | |
154 | ||
155 | switch (fBackgroundType) | |
156 | { | |
157 | case kSquareRoot: | |
158 | return NormSquareRoot(x, par); | |
159 | case kLinear: | |
160 | return NormLinear(x, par); | |
161 | case kPoly2: | |
162 | return NormPoly2(x, par); | |
163 | default: | |
164 | AliError("Invalid background function"); | |
165 | return 0.0; | |
166 | } | |
167 | } | |
168 | ||
169 | //__________________________________________________________________________________________________ | |
170 | Double_t AliRsnFitResult::Sum(Double_t *x, Double_t *par) | |
171 | { | |
172 | // | |
173 | // Computes the sum of the signal and the background chosen. | |
174 | // First 3 parameters are for the signal (they are always 3), | |
175 | // and other parameters, up to the necessary amount, are for | |
176 | // the background. | |
177 | // | |
178 | ||
179 | Double_t signal = Signal(x, par); | |
180 | Double_t background = Background(x, &par[3]); | |
181 | ||
182 | return signal + background; | |
183 | } | |
184 | ||
185 | //__________________________________________________________________________________________________ | |
186 | Int_t AliRsnFitResult::GetNParBackground() | |
187 | { | |
188 | // | |
189 | // Tells how many parameters has the background function | |
190 | // | |
191 | ||
192 | switch (fBackgroundType) | |
193 | { | |
194 | case kSquareRoot: return 2; | |
195 | case kLinear : return 2; | |
196 | case kPoly2 : return 3; | |
197 | default : return 0; | |
198 | } | |
199 | } | |
200 | ||
201 | //__________________________________________________________________________________________________ | |
202 | Bool_t AliRsnFitResult::InitFunctions() | |
203 | { | |
204 | // | |
205 | // Initialize functions | |
206 | // | |
207 | ||
208 | //if (!fSignalFcn) fSignalFcn = new TF1(Form("sg%s", GetName()), Signal, fFullRange[0], fFullRange[1], 3); | |
209 | //if (!fBackgroundFcn) fBackgroundFcn = new TF1(Form("bg%s", GetName()), Background, fFullRange[0], fFullRange[1], GetNParBackground()); | |
210 | //if (!fSumFcn) fSumFcn = new TF1(Form("sum%s", GetName()), Sum, fFullRange[0], fFullRange[1], 3+GetNParBackground()); | |
211 | } | |
212 | ||
213 | //__________________________________________________________________________________________________ | |
214 | void AliRsnFitResult::SetHistogram(TH1F *histogram) | |
215 | { | |
216 | // | |
217 | // Clones the passed histogram to proceed with fits | |
218 | // | |
219 | ||
220 | if (fHistogram) delete fHistogram; | |
221 | ||
b2b08ca2 | 222 | fHistogram = (TH1D*)histogram->Clone(); |
0d73200d | 223 | fHistogram->SetName(Form("h_%s_%s", GetName(), histogram->GetName())); |
224 | fHistogram->SetTitle(histogram->GetTitle()); | |
225 | } | |
226 | ||
227 | //__________________________________________________________________________________________________ | |
228 | Bool_t AliRsnFitResult::SingleFit(const char *opt, Double_t mass, Double_t width) | |
229 | { | |
230 | // | |
231 | // Executes a single fit of the function | |
232 | // | |
233 | /* | |
234 | // require histogram | |
235 | if (!fHistogram) | |
236 | { | |
237 | AliError("Require an initialized histogram!"); | |
238 | return kFALSE; | |
239 | } | |
240 | ||
241 | // if necessary, initialize functions | |
242 | if (!fSignalFcn || !fBackgroundFcn || !fSumFcn) InitFunctions(); | |
243 | ||
244 | // step #0: know how many parameter we have | |
245 | Int_t npar = GetNParBackground(); | |
246 | ||
247 | // step #1: | |
248 | // fit outside peak to roughly initialize the background | |
249 | for (Int_t i = 0; i < npar; i++) fBackgroundFcn->SetParameter(i, 0.5); | |
250 | status = (Int_t)h->Fit(fcnOut, opt); | |
251 | if (status) return kFALSE; | |
252 | ||
253 | // step #2: | |
254 | // estimate signal/background in the peak | |
255 | Int_t imax = h->GetMaximumBin(); | |
256 | Double_t xmax = h->GetBinCenter(imax); | |
257 | Double_t ymax = h->GetMaximum() - fBackgroundFcn->Eval(xmax); | |
258 | ||
259 | // step #3: | |
260 | // fit whole range with signal + background | |
261 | // which is initialized to the reference values for mass and width | |
262 | fSignalFcn->SetParameter(0, ymax); | |
263 | fSignalFcn->SetParameter(1, mass); | |
264 | fSignalFcn->SetParameter(2, width); | |
265 | for (Int_t i = 0; i < npar; i++) fcnSum->SetParameter(i + 3, fBackgroundFcn->GetParameter(i)); | |
266 | status = (Int_t)h->Fit(fcnSum, opt); | |
267 | if (status) return kFALSE; | |
268 | ||
269 | // get integral and its error here | |
270 | fValue[kSumIntegralError] = fSumFcn->IntegralError(fitPeakRangeMin, fitPeakRangeMax); | |
271 | fValue[kSumIntegral] = fSumFcn->Integral (fitPeakRangeMin, fitPeakRangeMax); | |
272 | fValue[kChi2] = fSumFcn->GetChisquare(); | |
273 | fValue[kNDF] = fSumFcn->GetNDF(); | |
274 | */ | |
275 | return kTRUE; | |
276 | } |