source: trunk/MagicSoft/Mars/mcalib/MCalibrationFits.h@ 2812

Last change on this file since 2812 was 2793, checked in by gaug, 22 years ago
*** empty log message ***
File size: 9.0 KB
Line 
1#ifndef MARS_MCalibrationFits
2#define MARS_MCalibrationFits
3
4#ifndef ROOT_TMath
5#include <TMath.h>
6#endif
7
8#define GIMMEABREAK 10000000000.0
9
10inline Double_t gfKto4(Double_t *x, Double_t *par)
11{
12
13 Double_t lambda = par[0];
14
15 Double_t sum = 0.;
16 Double_t arg = 0.;
17
18 Double_t mu0 = par[1];
19 Double_t mu1 = par[2];
20
21 if (mu1 < mu0)
22 return GIMMEABREAK;
23
24 Double_t sigma0 = par[3];
25 Double_t sigma1 = par[4];
26
27 if (sigma1 < sigma0)
28 return GIMMEABREAK;
29
30
31 Double_t mu2 = (2.*mu1)-mu0;
32 Double_t mu3 = (3.*mu1)-(2.*mu0);
33 Double_t mu4 = (4.*mu1)-(3.*mu0);
34
35 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
36 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
37 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
38
39 Double_t lambda2 = lambda*lambda;
40 Double_t lambda3 = lambda2*lambda;
41 Double_t lambda4 = lambda3*lambda;
42
43 // k=0:
44 arg = (x[0] - mu0)/sigma0;
45 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
46
47 // k=1:
48 arg = (x[0] - mu1)/sigma1;
49 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
50
51 // k=2:
52 arg = (x[0] - mu2)/sigma2;
53 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
54
55 // k=3:
56 arg = (x[0] - mu3)/sigma3;
57 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
58
59 // k=4:
60 arg = (x[0] - mu4)/sigma4;
61 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
62
63 return TMath::Exp(-1.*lambda)*par[5]*sum;
64
65};
66
67inline Double_t gfKto5(Double_t *x, Double_t *par)
68{
69
70 Double_t lambda = par[0];
71
72 Double_t sum = 0.;
73 Double_t arg = 0.;
74
75 Double_t mu0 = par[1];
76 Double_t mu1 = par[2];
77
78 if (mu1 < mu0)
79 return GIMMEABREAK;
80
81 Double_t sigma0 = par[3];
82 Double_t sigma1 = par[4];
83
84 if (sigma1 < sigma0)
85 return GIMMEABREAK;
86
87
88 Double_t mu2 = (2.*mu1)-mu0;
89 Double_t mu3 = (3.*mu1)-(2.*mu0);
90 Double_t mu4 = (4.*mu1)-(3.*mu0);
91 Double_t mu5 = (5.*mu1)-(4.*mu0);
92
93 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
94 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
95 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
96 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
97
98 Double_t lambda2 = lambda*lambda;
99 Double_t lambda3 = lambda2*lambda;
100 Double_t lambda4 = lambda3*lambda;
101 Double_t lambda5 = lambda4*lambda;
102
103 // k=0:
104 arg = (x[0] - mu0)/sigma0;
105 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
106
107 // k=1:
108 arg = (x[0] - mu1)/sigma1;
109 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
110
111 // k=2:
112 arg = (x[0] - mu2)/sigma2;
113 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
114
115 // k=3:
116 arg = (x[0] - mu3)/sigma3;
117 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
118
119 // k=4:
120 arg = (x[0] - mu4)/sigma4;
121 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
122
123 // k=5:
124 arg = (x[0] - mu5)/sigma5;
125 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
126
127 return TMath::Exp(-1.*lambda)*par[5]*sum;
128
129};
130
131
132inline Double_t gfKto6(Double_t *x, Double_t *par)
133{
134
135 Double_t lambda = par[0];
136
137 Double_t sum = 0.;
138 Double_t arg = 0.;
139
140 Double_t mu0 = par[1];
141 Double_t mu1 = par[2];
142
143 if (mu1 < mu0)
144 return GIMMEABREAK;
145
146 Double_t sigma0 = par[3];
147 Double_t sigma1 = par[4];
148
149 if (sigma1 < sigma0)
150 return GIMMEABREAK;
151
152
153 Double_t mu2 = (2.*mu1)-mu0;
154 Double_t mu3 = (3.*mu1)-(2.*mu0);
155 Double_t mu4 = (4.*mu1)-(3.*mu0);
156 Double_t mu5 = (5.*mu1)-(4.*mu0);
157 Double_t mu6 = (6.*mu1)-(5.*mu0);
158
159 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
160 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
161 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
162 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
163 Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
164
165 Double_t lambda2 = lambda*lambda;
166 Double_t lambda3 = lambda2*lambda;
167 Double_t lambda4 = lambda3*lambda;
168 Double_t lambda5 = lambda4*lambda;
169 Double_t lambda6 = lambda5*lambda;
170
171 // k=0:
172 arg = (x[0] - mu0)/sigma0;
173 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
174
175 // k=1:
176 arg = (x[0] - mu1)/sigma1;
177 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
178
179 // k=2:
180 arg = (x[0] - mu2)/sigma2;
181 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
182
183 // k=3:
184 arg = (x[0] - mu3)/sigma3;
185 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
186
187 // k=4:
188 arg = (x[0] - mu4)/sigma4;
189 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
190
191 // k=5:
192 arg = (x[0] - mu5)/sigma5;
193 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
194
195 // k=6:
196 arg = (x[0] - mu6)/sigma6;
197 sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
198
199 return TMath::Exp(-1.*lambda)*par[5]*sum;
200
201};
202
203inline Double_t gfKto7(Double_t *x, Double_t *par)
204{
205
206 Double_t lambda = par[0];
207
208 Double_t sum = 0.;
209 Double_t arg = 0.;
210
211 Double_t mu0 = par[1];
212 Double_t mu1 = par[2];
213
214 if (mu1 < mu0)
215 return GIMMEABREAK;
216
217 Double_t sigma0 = par[3];
218 Double_t sigma1 = par[4];
219
220 if (sigma1 < sigma0)
221 return GIMMEABREAK;
222
223
224 Double_t mu2 = (2.*mu1)-mu0;
225 Double_t mu3 = (3.*mu1)-(2.*mu0);
226 Double_t mu4 = (4.*mu1)-(3.*mu0);
227 Double_t mu5 = (5.*mu1)-(4.*mu0);
228 Double_t mu6 = (6.*mu1)-(5.*mu0);
229 Double_t mu7 = (7.*mu1)-(6.*mu0);
230
231 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
232 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
233 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
234 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
235 Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
236 Double_t sigma7 = TMath::Sqrt((7.*sigma1*sigma1) - (6.*sigma0*sigma0));
237
238 Double_t lambda2 = lambda*lambda;
239 Double_t lambda3 = lambda2*lambda;
240 Double_t lambda4 = lambda3*lambda;
241 Double_t lambda5 = lambda4*lambda;
242 Double_t lambda6 = lambda5*lambda;
243 Double_t lambda7 = lambda6*lambda;
244
245 // k=0:
246 arg = (x[0] - mu0)/sigma0;
247 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
248
249 // k=1:
250 arg = (x[0] - mu1)/sigma1;
251 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
252
253 // k=2:
254 arg = (x[0] - mu2)/sigma2;
255 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
256
257 // k=3:
258 arg = (x[0] - mu3)/sigma3;
259 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
260
261 // k=4:
262 arg = (x[0] - mu4)/sigma4;
263 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
264
265 // k=5:
266 arg = (x[0] - mu5)/sigma5;
267 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
268
269 // k=6:
270 arg = (x[0] - mu6)/sigma6;
271 sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
272
273 // k=7:
274 arg = (x[0] - mu7)/sigma7;
275 sum += 0.000198412698413*lambda7*TMath::Exp(-0.5*arg*arg)/sigma7;
276
277 return TMath::Exp(-1.*lambda)*par[5]*sum;
278
279};
280
281
282inline Double_t gfKto8(Double_t *x, Double_t *par)
283{
284
285 Double_t lambda = par[0];
286
287 Double_t sum = 0.;
288 Double_t arg = 0.;
289
290 Double_t mu0 = par[1];
291 Double_t mu1 = par[2];
292
293 if (mu1 < mu0)
294 return GIMMEABREAK;
295
296 Double_t sigma0 = par[3];
297 Double_t sigma1 = par[4];
298
299 if (sigma1 < sigma0)
300 return GIMMEABREAK;
301
302
303 Double_t mu2 = (2.*mu1)-mu0;
304 Double_t mu3 = (3.*mu1)-(2.*mu0);
305 Double_t mu4 = (4.*mu1)-(3.*mu0);
306 Double_t mu5 = (5.*mu1)-(4.*mu0);
307 Double_t mu6 = (6.*mu1)-(5.*mu0);
308 Double_t mu7 = (7.*mu1)-(6.*mu0);
309 Double_t mu8 = (8.*mu1)-(7.*mu0);
310
311 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
312 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
313 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
314 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
315 Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
316 Double_t sigma7 = TMath::Sqrt((7.*sigma1*sigma1) - (6.*sigma0*sigma0));
317 Double_t sigma8 = TMath::Sqrt((8.*sigma1*sigma1) - (7.*sigma0*sigma0));
318
319 Double_t lambda2 = lambda*lambda;
320 Double_t lambda3 = lambda2*lambda;
321 Double_t lambda4 = lambda3*lambda;
322 Double_t lambda5 = lambda4*lambda;
323 Double_t lambda6 = lambda5*lambda;
324 Double_t lambda7 = lambda6*lambda;
325 Double_t lambda8 = lambda7*lambda;
326
327 // k=0:
328 arg = (x[0] - mu0)/sigma0;
329 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
330
331 // k=1:
332 arg = (x[0] - mu1)/sigma1;
333 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
334
335 // k=2:
336 arg = (x[0] - mu2)/sigma2;
337 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
338
339 // k=3:
340 arg = (x[0] - mu3)/sigma3;
341 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
342
343 // k=4:
344 arg = (x[0] - mu4)/sigma4;
345 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
346
347 // k=5:
348 arg = (x[0] - mu5)/sigma5;
349 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
350
351 // k=6:
352 arg = (x[0] - mu6)/sigma6;
353 sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
354
355 // k=7:
356 arg = (x[0] - mu7)/sigma7;
357 sum += 0.000198412698413*lambda7*TMath::Exp(-0.5*arg*arg)/sigma7;
358
359 // k=8:
360 arg = (x[0] - mu8)/sigma8;
361 sum += 0.000024801587315*lambda8*TMath::Exp(-0.5*arg*arg)/sigma7;
362
363 return TMath::Exp(-1.*lambda)*par[5]*sum;
364
365};
366
367#endif /* MARS_MCalibrationFits */
368
Note: See TracBrowser for help on using the repository browser.