source: trunk/MagicSoft/Mars/manalysis/MCalibrationFits.h@ 2539

Last change on this file since 2539 was 2525, checked in by gaug, 21 years ago
*** empty log message ***
File size: 7.7 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 gfKto5(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 Double_t mu5 = (5.*mu1)-(4.*mu0);
35
36 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
37 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
38 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
39 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
40
41 Double_t lambda2 = lambda*lambda;
42 Double_t lambda3 = lambda2*lambda;
43 Double_t lambda4 = lambda3*lambda;
44 Double_t lambda5 = lambda4*lambda;
45
46 // k=0:
47 arg = (x[0] - mu0)/sigma0;
48 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
49
50 // k=1:
51 arg = (x[0] - mu1)/sigma1;
52 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
53
54 // k=2:
55 arg = (x[0] - mu2)/sigma2;
56 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
57
58 // k=3:
59 arg = (x[0] - mu3)/sigma3;
60 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
61
62 // k=4:
63 arg = (x[0] - mu4)/sigma4;
64 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
65
66 // k=5:
67 arg = (x[0] - mu5)/sigma5;
68 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
69
70 return par[5]*sum;
71
72};
73
74inline Double_t gfKto6(Double_t *x, Double_t *par)
75{
76
77 Double_t lambda = par[0];
78
79 Double_t sum = 0.;
80 Double_t arg = 0.;
81
82 Double_t mu0 = par[1];
83 Double_t mu1 = par[2];
84
85 if (mu1 < mu0)
86 return GIMMEABREAK;
87
88 Double_t sigma0 = par[3];
89 Double_t sigma1 = par[4];
90
91 if (sigma1 < sigma0)
92 return GIMMEABREAK;
93
94
95 Double_t mu2 = (2.*mu1)-mu0;
96 Double_t mu3 = (3.*mu1)-(2.*mu0);
97 Double_t mu4 = (4.*mu1)-(3.*mu0);
98 Double_t mu5 = (5.*mu1)-(4.*mu0);
99 Double_t mu6 = (6.*mu1)-(5.*mu0);
100
101 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
102 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
103 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
104 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
105 Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
106
107 Double_t lambda2 = lambda*lambda;
108 Double_t lambda3 = lambda2*lambda;
109 Double_t lambda4 = lambda3*lambda;
110 Double_t lambda5 = lambda4*lambda;
111 Double_t lambda6 = lambda5*lambda;
112
113 // k=0:
114 arg = (x[0] - mu0)/sigma0;
115 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
116
117 // k=1:
118 arg = (x[0] - mu1)/sigma1;
119 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
120
121 // k=2:
122 arg = (x[0] - mu2)/sigma2;
123 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
124
125 // k=3:
126 arg = (x[0] - mu3)/sigma3;
127 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
128
129 // k=4:
130 arg = (x[0] - mu4)/sigma4;
131 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
132
133 // k=5:
134 arg = (x[0] - mu5)/sigma5;
135 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
136
137 // k=6:
138 arg = (x[0] - mu6)/sigma6;
139 sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
140
141 return par[5]*sum;
142
143};
144
145inline Double_t gfKto7(Double_t *x, Double_t *par)
146{
147
148 Double_t lambda = par[0];
149
150 Double_t sum = 0.;
151 Double_t arg = 0.;
152
153 Double_t mu0 = par[1];
154 Double_t mu1 = par[2];
155
156 if (mu1 < mu0)
157 return GIMMEABREAK;
158
159 Double_t sigma0 = par[3];
160 Double_t sigma1 = par[4];
161
162 if (sigma1 < sigma0)
163 return GIMMEABREAK;
164
165
166 Double_t mu2 = (2.*mu1)-mu0;
167 Double_t mu3 = (3.*mu1)-(2.*mu0);
168 Double_t mu4 = (4.*mu1)-(3.*mu0);
169 Double_t mu5 = (5.*mu1)-(4.*mu0);
170 Double_t mu6 = (6.*mu1)-(5.*mu0);
171 Double_t mu7 = (7.*mu1)-(6.*mu0);
172
173 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
174 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
175 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
176 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
177 Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
178 Double_t sigma7 = TMath::Sqrt((7.*sigma1*sigma1) - (6.*sigma0*sigma0));
179
180 Double_t lambda2 = lambda*lambda;
181 Double_t lambda3 = lambda2*lambda;
182 Double_t lambda4 = lambda3*lambda;
183 Double_t lambda5 = lambda4*lambda;
184 Double_t lambda6 = lambda5*lambda;
185 Double_t lambda7 = lambda6*lambda;
186
187 // k=0:
188 arg = (x[0] - mu0)/sigma0;
189 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
190
191 // k=1:
192 arg = (x[0] - mu1)/sigma1;
193 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
194
195 // k=2:
196 arg = (x[0] - mu2)/sigma2;
197 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
198
199 // k=3:
200 arg = (x[0] - mu3)/sigma3;
201 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
202
203 // k=4:
204 arg = (x[0] - mu4)/sigma4;
205 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
206
207 // k=5:
208 arg = (x[0] - mu5)/sigma5;
209 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
210
211 // k=6:
212 arg = (x[0] - mu6)/sigma6;
213 sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
214
215 // k=7:
216 arg = (x[0] - mu7)/sigma7;
217 sum += 0.000198412698413*lambda7*TMath::Exp(-0.5*arg*arg)/sigma7;
218
219 return par[5]*sum;
220
221};
222
223
224inline Double_t gfKto8(Double_t *x, Double_t *par)
225{
226
227 Double_t lambda = par[0];
228
229 Double_t sum = 0.;
230 Double_t arg = 0.;
231
232 Double_t mu0 = par[1];
233 Double_t mu1 = par[2];
234
235 if (mu1 < mu0)
236 return GIMMEABREAK;
237
238 Double_t sigma0 = par[3];
239 Double_t sigma1 = par[4];
240
241 if (sigma1 < sigma0)
242 return GIMMEABREAK;
243
244
245 Double_t mu2 = (2.*mu1)-mu0;
246 Double_t mu3 = (3.*mu1)-(2.*mu0);
247 Double_t mu4 = (4.*mu1)-(3.*mu0);
248 Double_t mu5 = (5.*mu1)-(4.*mu0);
249 Double_t mu6 = (6.*mu1)-(5.*mu0);
250 Double_t mu7 = (7.*mu1)-(6.*mu0);
251 Double_t mu8 = (8.*mu1)-(7.*mu0);
252
253 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
254 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
255 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
256 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
257 Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
258 Double_t sigma7 = TMath::Sqrt((7.*sigma1*sigma1) - (6.*sigma0*sigma0));
259 Double_t sigma8 = TMath::Sqrt((8.*sigma1*sigma1) - (7.*sigma0*sigma0));
260
261 Double_t lambda2 = lambda*lambda;
262 Double_t lambda3 = lambda2*lambda;
263 Double_t lambda4 = lambda3*lambda;
264 Double_t lambda5 = lambda4*lambda;
265 Double_t lambda6 = lambda5*lambda;
266 Double_t lambda7 = lambda6*lambda;
267 Double_t lambda8 = lambda7*lambda;
268
269 // k=0:
270 arg = (x[0] - mu0)/sigma0;
271 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
272
273 // k=1:
274 arg = (x[0] - mu1)/sigma1;
275 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
276
277 // k=2:
278 arg = (x[0] - mu2)/sigma2;
279 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
280
281 // k=3:
282 arg = (x[0] - mu3)/sigma3;
283 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
284
285 // k=4:
286 arg = (x[0] - mu4)/sigma4;
287 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
288
289 // k=5:
290 arg = (x[0] - mu5)/sigma5;
291 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
292
293 // k=6:
294 arg = (x[0] - mu6)/sigma6;
295 sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
296
297 // k=7:
298 arg = (x[0] - mu7)/sigma7;
299 sum += 0.000198412698413*lambda7*TMath::Exp(-0.5*arg*arg)/sigma7;
300
301 // k=8:
302 arg = (x[0] - mu8)/sigma8;
303 sum += 0.000024801587315*lambda8*TMath::Exp(-0.5*arg*arg)/sigma7;
304
305 return par[5]*sum;
306
307};
308
309#endif /* MARS_MCalibrationFits */
310
Note: See TracBrowser for help on using the repository browser.