source: trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h@ 2852

Last change on this file since 2852 was 2852, checked in by gaug, 21 years ago
*** empty log message ***
File size: 10.8 KB
Line 
1#ifndef MARS_MHCalibrationBlindPixel
2#define MARS_MHCalibrationBlindPixel
3
4#ifndef MARS_MH
5#include "MH.h"
6#endif
7
8#ifndef MARS_MHCalibrationConfig
9#include "MHCalibrationConfig.h"
10#endif
11
12#ifndef ROOT_TH1
13#include "TH1.h"
14#endif
15
16#ifndef ROOT_TH1F
17#include "TH1F.h"
18#endif
19
20#ifndef ROOT_TF1
21#include "TF1.h"
22#endif
23
24#ifndef ROOT_TPaveText
25#include "TPaveText.h"
26#endif
27
28const Double_t NoWay = 10000000000.0;
29
30class TMath;
31class MParList;
32class MHCalibrationBlindPixel : public MH
33{
34private:
35
36 TH1F* fHBlindPixelCharge; //-> Histogram with the single Phe spectrum
37 TH1I* fHBlindPixelTime; //-> Variance of summed FADC slices
38 TH1I* fHBlindPixelChargevsN; //-> Summed Charge vs. Event Nr.
39
40 TF1 fSinglePheFit;
41 TF1 fTimeGausFit;
42 TF1 fSinglePhePedFit;
43
44 Axis_t fBlindPixelChargefirst;
45 Axis_t fBlindPixelChargelast;
46 Int_t fBlindPixelChargenbins;
47
48 void ResetBin(Int_t i);
49 void DrawLegend();
50
51 TPaveText *fFitLegend; //!
52 Bool_t fFitOK;
53
54 Double_t fLambda;
55 Double_t fMu0;
56 Double_t fMu1;
57 Double_t fSigma0;
58 Double_t fSigma1;
59
60 Double_t fLambdaErr;
61 Double_t fMu0Err;
62 Double_t fMu1Err;
63 Double_t fSigma0Err;
64 Double_t fSigma1Err;
65
66 Double_t fChisquare;
67 Double_t fProb;
68 Int_t fNdf;
69
70 Double_t fMeanTime;
71 Double_t fMeanTimeErr;
72 Double_t fSigmaTime;
73 Double_t fSigmaTimeErr;
74
75 Double_t fLambdaCheck;
76 Double_t fLambdaCheckErr;
77
78public:
79
80 MHCalibrationBlindPixel(const char *name=NULL, const char *title=NULL);
81 ~MHCalibrationBlindPixel();
82
83 Bool_t FillBlindPixelCharge(Float_t q) { return fHBlindPixelCharge->Fill(q) > -1; }
84 Bool_t FillBlindPixelTime(Int_t t) { return fHBlindPixelTime->Fill(t) > -1; }
85 Bool_t FillBlindPixelChargevsN(Stat_t rq, Int_t t) { return fHBlindPixelChargevsN->Fill(t,rq) > -1; }
86
87
88 //Getters
89 const Double_t GetLambda() const { return fLambda; }
90 const Double_t GetLambdaCheck() const { return fLambdaCheck; }
91 const Double_t GetMu0() const { return fMu0; }
92 const Double_t GetMu1() const { return fMu1; }
93 const Double_t GetSigma0() const { return fSigma0; }
94 const Double_t GetSigma1() const { return fSigma1; }
95
96 const Double_t GetLambdaErr() const { return fLambdaErr; }
97 const Double_t GetLambdaCheckErr() const { return fLambdaCheckErr; }
98 const Double_t GetMu0Err() const { return fMu0Err; }
99 const Double_t GetMu1Err() const { return fMu1Err; }
100 const Double_t GetSigma0Err() const { return fSigma0Err; }
101 const Double_t GetSigma1Err() const { return fSigma1Err; }
102
103 const Double_t GetChiSquare() const { return fChisquare; }
104 const Double_t GetProb() const { return fProb; }
105 const Int_t GetNdf() const { return fNdf; }
106
107 const Double_t GetMeanTime() const { return fMeanTime; }
108 const Double_t GetMeanTimeErr() const { return fMeanTimeErr; }
109 const Double_t GetSigmaTime() const { return fSigmaTime; }
110 const Double_t GetSigmaTimeErr() const { return fSigmaTimeErr; }
111
112 const Bool_t IsFitOK() { return fFitOK; }
113
114 // Draws
115 TObject *DrawClone(Option_t *option="") const;
116 void Draw(Option_t *option="");
117
118 // Fits
119 enum FitFunc_t { kEPoisson4, kEPoisson5, kEPoisson6, kEPoisson7, kEPolya, kEMichele };
120
121private:
122 FitFunc_t fFitFunc;
123
124public:
125 Bool_t FitSinglePhe(Axis_t rmin=0, Axis_t rmax=0, Option_t *opt="RL0+Q");
126 Bool_t FitTime(Axis_t rmin=0., Axis_t rmax=0.,Option_t *opt="R0+Q");
127 void ChangeFitFunc(FitFunc_t func) { fFitFunc = func; }
128
129 // Simulation
130 Bool_t SimulateSinglePhe(Double_t lambda,
131 Double_t mu0,Double_t mu1,
132 Double_t sigma0,Double_t sigma1);
133
134 // Others
135 void CutAllEdges();
136
137private:
138
139 void InitFit(TF1& f, Axis_t min, Axis_t max);
140 void ExitFit(TF1& f);
141
142 inline static Double_t fFitFuncMichele(Double_t *x, Double_t *par)
143 {
144
145 Double_t lambda = par[0];
146
147 Double_t sum = 0.;
148 Double_t arg = 0.;
149
150 Double_t mu0 = par[1];
151 Double_t mu1 = par[2];
152
153 if (mu1 < mu0)
154 return NoWay;
155
156 Double_t sigma0 = par[3];
157 Double_t sigma1 = par[4];
158
159 if (sigma1 < sigma0)
160 return NoWay;
161
162 Double_t mu2 = (2.*mu1)-mu0;
163 Double_t mu3 = (3.*mu1)-(2.*mu0);
164 Double_t mu4 = (4.*mu1)-(3.*mu0);
165
166 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
167 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
168 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
169
170 Double_t lambda2 = lambda*lambda;
171 Double_t lambda3 = lambda2*lambda;
172 Double_t lambda4 = lambda3*lambda;
173
174 // k=0:
175 arg = (x[0] - mu0)/sigma0;
176 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
177
178 // k=1:
179 arg = (x[0] - mu1)/sigma1;
180 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
181
182 // k=2:
183 arg = (x[0] - mu2)/sigma2;
184 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
185
186 // k=3:
187 arg = (x[0] - mu3)/sigma3;
188 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
189
190 // k=4:
191 arg = (x[0] - mu4)/sigma4;
192 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
193
194 return TMath::Exp(-1.*lambda)*par[5]*sum;
195
196 }
197
198 inline static Double_t fPoissonKto4(Double_t *x, Double_t *par)
199 {
200
201 Double_t lambda = par[0];
202
203 Double_t sum = 0.;
204 Double_t arg = 0.;
205
206 Double_t mu0 = par[1];
207 Double_t mu1 = par[2];
208
209 if (mu1 < mu0)
210 return NoWay;
211
212 Double_t sigma0 = par[3];
213 Double_t sigma1 = par[4];
214
215 if (sigma1 < sigma0)
216 return NoWay;
217
218 Double_t mu2 = (2.*mu1)-mu0;
219 Double_t mu3 = (3.*mu1)-(2.*mu0);
220 Double_t mu4 = (4.*mu1)-(3.*mu0);
221
222 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
223 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
224 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
225
226 Double_t lambda2 = lambda*lambda;
227 Double_t lambda3 = lambda2*lambda;
228 Double_t lambda4 = lambda3*lambda;
229
230 // k=0:
231 arg = (x[0] - mu0)/sigma0;
232 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
233
234 // k=1:
235 arg = (x[0] - mu1)/sigma1;
236 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
237
238 // k=2:
239 arg = (x[0] - mu2)/sigma2;
240 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
241
242 // k=3:
243 arg = (x[0] - mu3)/sigma3;
244 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
245
246 // k=4:
247 arg = (x[0] - mu4)/sigma4;
248 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
249
250 return TMath::Exp(-1.*lambda)*par[5]*sum;
251
252 }
253
254
255 inline static Double_t fPoissonKto5(Double_t *x, Double_t *par)
256 {
257
258 Double_t lambda = par[0];
259
260 Double_t sum = 0.;
261 Double_t arg = 0.;
262
263 Double_t mu0 = par[1];
264 Double_t mu1 = par[2];
265
266 if (mu1 < mu0)
267 return NoWay;
268
269 Double_t sigma0 = par[3];
270 Double_t sigma1 = par[4];
271
272 if (sigma1 < sigma0)
273 return NoWay;
274
275
276 Double_t mu2 = (2.*mu1)-mu0;
277 Double_t mu3 = (3.*mu1)-(2.*mu0);
278 Double_t mu4 = (4.*mu1)-(3.*mu0);
279 Double_t mu5 = (5.*mu1)-(4.*mu0);
280
281 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
282 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
283 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
284 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
285
286 Double_t lambda2 = lambda*lambda;
287 Double_t lambda3 = lambda2*lambda;
288 Double_t lambda4 = lambda3*lambda;
289 Double_t lambda5 = lambda4*lambda;
290
291 // k=0:
292 arg = (x[0] - mu0)/sigma0;
293 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
294
295 // k=1:
296 arg = (x[0] - mu1)/sigma1;
297 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
298
299 // k=2:
300 arg = (x[0] - mu2)/sigma2;
301 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
302
303 // k=3:
304 arg = (x[0] - mu3)/sigma3;
305 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
306
307 // k=4:
308 arg = (x[0] - mu4)/sigma4;
309 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
310
311 // k=5:
312 arg = (x[0] - mu5)/sigma5;
313 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
314
315 return TMath::Exp(-1.*lambda)*par[5]*sum;
316
317 }
318
319
320 inline static Double_t fPoissonKto6(Double_t *x, Double_t *par)
321 {
322
323 Double_t lambda = par[0];
324
325 Double_t sum = 0.;
326 Double_t arg = 0.;
327
328 Double_t mu0 = par[1];
329 Double_t mu1 = par[2];
330
331 if (mu1 < mu0)
332 return NoWay;
333
334 Double_t sigma0 = par[3];
335 Double_t sigma1 = par[4];
336
337 if (sigma1 < sigma0)
338 return NoWay;
339
340
341 Double_t mu2 = (2.*mu1)-mu0;
342 Double_t mu3 = (3.*mu1)-(2.*mu0);
343 Double_t mu4 = (4.*mu1)-(3.*mu0);
344 Double_t mu5 = (5.*mu1)-(4.*mu0);
345 Double_t mu6 = (6.*mu1)-(5.*mu0);
346
347 Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0));
348 Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0));
349 Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0));
350 Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0));
351 Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0));
352
353 Double_t lambda2 = lambda*lambda;
354 Double_t lambda3 = lambda2*lambda;
355 Double_t lambda4 = lambda3*lambda;
356 Double_t lambda5 = lambda4*lambda;
357 Double_t lambda6 = lambda5*lambda;
358
359 // k=0:
360 arg = (x[0] - mu0)/sigma0;
361 sum = TMath::Exp(-0.5*arg*arg)/sigma0;
362
363 // k=1:
364 arg = (x[0] - mu1)/sigma1;
365 sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1;
366
367 // k=2:
368 arg = (x[0] - mu2)/sigma2;
369 sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2;
370
371 // k=3:
372 arg = (x[0] - mu3)/sigma3;
373 sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3;
374
375 // k=4:
376 arg = (x[0] - mu4)/sigma4;
377 sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4;
378
379 // k=5:
380 arg = (x[0] - mu5)/sigma5;
381 sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5;
382
383 // k=6:
384 arg = (x[0] - mu6)/sigma6;
385 sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6;
386
387 return TMath::Exp(-1.*lambda)*par[5]*sum;
388
389 }
390
391 ClassDef(MHCalibrationBlindPixel, 1) // Histograms from the Calibration Blind Pixel
392};
393
394#endif /* MARS_MHCalibrationBlindPixel */
Note: See TracBrowser for help on using the repository browser.