1 | /* ======================================================================== *\
|
---|
2 | !
|
---|
3 | ! *
|
---|
4 | ! * This file is part of MARS, the MAGIC Analysis and Reconstruction
|
---|
5 | ! * Software. It is distributed to you in the hope that it can be a useful
|
---|
6 | ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
|
---|
7 | ! * It is distributed WITHOUT ANY WARRANTY.
|
---|
8 | ! *
|
---|
9 | ! * Permission to use, copy, modify and distribute this software and its
|
---|
10 | ! * documentation for any purpose is hereby granted without fee,
|
---|
11 | ! * provided that the above copyright notice appear in all copies and
|
---|
12 | ! * that both that copyright notice and this permission notice appear
|
---|
13 | ! * in supporting documentation. It is provided "as is" without express
|
---|
14 | ! * or implied warranty.
|
---|
15 | ! *
|
---|
16 | !
|
---|
17 | !
|
---|
18 | ! Author(s): Markus Gaug 11/2003 <mailto:markus@ifae.es>
|
---|
19 | !
|
---|
20 | ! Copyright: MAGIC Software Development, 2000-2001
|
---|
21 | !
|
---|
22 | !
|
---|
23 | \* ======================================================================== */
|
---|
24 |
|
---|
25 | /////////////////////////////////////////////////////////////////////////////
|
---|
26 | // //
|
---|
27 | // MCalibrationPix //
|
---|
28 | // //
|
---|
29 | // This is the storage container to hold informations about the pedestal //
|
---|
30 | // (offset) value of one Pixel (PMT). //
|
---|
31 | // //
|
---|
32 | // The following values are initialized to meaningful values:
|
---|
33 | //
|
---|
34 | // - The Electronic Rms to 1.5 per FADC slice
|
---|
35 | // - The uncertainty about the Electronic RMS to 0.3 per slice
|
---|
36 | // - The F-Factor is assumed to have been measured in Munich to 1.13 - 1.17.
|
---|
37 | // with the Munich definition of the F-Factor, thus:
|
---|
38 | // F = Sigma(Out)/Mean(Out) * Mean(In)/Sigma(In)
|
---|
39 | // Mean F-Factor = 1.15
|
---|
40 | // Error F-Factor = 0.02
|
---|
41 | //
|
---|
42 | // - Average QE: (email David Paneque, 14.2.04):
|
---|
43 | //
|
---|
44 | // The conversion factor that comes purely from QE folded to a Cherenkov
|
---|
45 | // spectrum has to be multiplied by:
|
---|
46 | // * Plexiglass window -->> 0.96 X 0.96
|
---|
47 | // * PMT photoelectron collection efficiency -->> 0.9
|
---|
48 | // * Light guides efficiency -->> 0.94
|
---|
49 | //
|
---|
50 | // Concerning the light guides effiency estimation... Daniel Ferenc
|
---|
51 | // is preparing some work (simulations) to estimate it. Yet so far, he has
|
---|
52 | // been busy with other stuff, and this work is still UNfinished.
|
---|
53 | //
|
---|
54 | // The estimation I did comes from:
|
---|
55 | // 1) Reflectivity of light guide walls is 85 % (aluminum)
|
---|
56 | // 2) At ZERO degree light incidence, 37% of the light hits such walls
|
---|
57 | // (0.15X37%= 5.6% of light lost)
|
---|
58 | // 3) When increasing the light incidence angle, more and more light hits
|
---|
59 | // the walls.
|
---|
60 | //
|
---|
61 | // However, the loses due to larger amount of photons hitting the walls is more
|
---|
62 | // or less counteracted by the fact that more and more photon trajectories cross
|
---|
63 | // the PMT photocathode twice, increasing the effective sensitivity of the PMT.
|
---|
64 | //
|
---|
65 | // Jurgen Gebauer did some quick measurements about this issue. I attach a
|
---|
66 | // plot. You can see that the angular dependence is (more or less) in agreement
|
---|
67 | // with a CosTheta function (below 20-25 degrees),
|
---|
68 | // which is the variation of teh entrance window cross section. So, in
|
---|
69 | // first approximation, no loses when increasing light incidence angle;
|
---|
70 | // and therefore, the factor 0.94.
|
---|
71 | //
|
---|
72 | // So, summarizing... I would propose the following conversion factors
|
---|
73 | // (while working with CT1 cal box) in order to get the final number of photons
|
---|
74 | // from the detected measured size in ADC counts.
|
---|
75 | //
|
---|
76 | // Nph = ADC * FmethodConversionFactor * ConvPhe-PhFactor
|
---|
77 | //
|
---|
78 | // FmethodConversionFactor ; measured for individual pmts
|
---|
79 | //
|
---|
80 | // ConvPhe-PhFactor = 0.98 * 0.23 * 0.90 * 0.94 * 0.96 * 0.96 = 0.18
|
---|
81 | //
|
---|
82 | // I would not apply any smearing of this factor (which we have in nature),
|
---|
83 | // since we might be applying it to PMTs in the totally wrong direction.
|
---|
84 | //
|
---|
85 | //
|
---|
86 | /////////////////////////////////////////////////////////////////////////////
|
---|
87 | #include "MCalibrationPix.h"
|
---|
88 |
|
---|
89 | #include "MLog.h"
|
---|
90 | #include "MLogManip.h"
|
---|
91 |
|
---|
92 | ClassImp(MCalibrationPix);
|
---|
93 |
|
---|
94 | using namespace std;
|
---|
95 |
|
---|
96 | const Float_t MCalibrationPix::gkElectronicPedRms = 1.5;
|
---|
97 | const Float_t MCalibrationPix::gkElectronicPedRmsErr = 0.3;
|
---|
98 | const Float_t MCalibrationPix::gkFFactor = 1.15;
|
---|
99 | const Float_t MCalibrationPix::gkFFactorErr = 0.02;
|
---|
100 | const Float_t MCalibrationPix::gkChargeLimit = 3.;
|
---|
101 | const Float_t MCalibrationPix::gkChargeErrLimit = 0.;
|
---|
102 | const Float_t MCalibrationPix::gkChargeRelErrLimit = 1.;
|
---|
103 | const Float_t MCalibrationPix::gkTimeLimit = 1.5;
|
---|
104 | const Float_t MCalibrationPix::gkTimeErrLimit = 3.;
|
---|
105 | const Float_t MCalibrationPix::gkConvFFactorRelErrLimit = 0.1;
|
---|
106 |
|
---|
107 | const Float_t MCalibrationPix::gkAverageQE = 0.18;
|
---|
108 | const Float_t MCalibrationPix::gkAverageQEErr = 0.02;
|
---|
109 | const Float_t MCalibrationPix::gkConversionHiLo = 10.;
|
---|
110 | const Float_t MCalibrationPix::gkConversionHiLoErr = 2.5;
|
---|
111 | // --------------------------------------------------------------------------
|
---|
112 | //
|
---|
113 | // Default Constructor:
|
---|
114 | //
|
---|
115 | MCalibrationPix::MCalibrationPix(const char *name, const char *title)
|
---|
116 | : fPixId(-1),
|
---|
117 | fFlags(0)
|
---|
118 | {
|
---|
119 |
|
---|
120 | fName = name ? name : "MCalibrationPixel";
|
---|
121 | fTitle = title ? title : "Container of the MHCalibrationPixels and the fit results";
|
---|
122 |
|
---|
123 | fHist = new MHCalibrationPixel("MHCalibrationPixel","Calibration Histograms Pixel ");
|
---|
124 |
|
---|
125 | if (!fHist)
|
---|
126 | *fLog << warn << dbginf << " Could not create MHCalibrationPixel " << endl;
|
---|
127 |
|
---|
128 | Clear();
|
---|
129 |
|
---|
130 | //
|
---|
131 | // At the moment, we don't have a database, yet,
|
---|
132 | // so we get it from the configuration file
|
---|
133 | //
|
---|
134 | SetConversionHiLo();
|
---|
135 | SetConversionHiLoErr();
|
---|
136 |
|
---|
137 | SetAverageQE();
|
---|
138 | }
|
---|
139 |
|
---|
140 | MCalibrationPix::~MCalibrationPix()
|
---|
141 | {
|
---|
142 | delete fHist;
|
---|
143 | }
|
---|
144 |
|
---|
145 |
|
---|
146 | // ------------------------------------------------------------------------
|
---|
147 | //
|
---|
148 | // Invalidate values
|
---|
149 | //
|
---|
150 | void MCalibrationPix::Clear(Option_t *o)
|
---|
151 | {
|
---|
152 |
|
---|
153 | fHist->Reset();
|
---|
154 |
|
---|
155 | CLRBIT(fFlags, kHiGainSaturation);
|
---|
156 | CLRBIT(fFlags, kExcluded);
|
---|
157 | CLRBIT(fFlags, kExcludeQualityCheck);
|
---|
158 | CLRBIT(fFlags, kChargeValid);
|
---|
159 | CLRBIT(fFlags, kFitted);
|
---|
160 | CLRBIT(fFlags, kOscillating);
|
---|
161 | CLRBIT(fFlags, kBlindPixelMethodValid);
|
---|
162 | CLRBIT(fFlags, kFFactorMethodValid);
|
---|
163 | CLRBIT(fFlags, kPINDiodeMethodValid);
|
---|
164 | CLRBIT(fFlags, kCombinedMethodValid);
|
---|
165 |
|
---|
166 | fCharge = -1.;
|
---|
167 | fChargeErr = -1.;
|
---|
168 | fSigmaCharge = -1.;
|
---|
169 | fSigmaChargeErr = -1.;
|
---|
170 | fRSigmaCharge = -1.;
|
---|
171 | fRSigmaChargeErr = -1.;
|
---|
172 |
|
---|
173 | fChargeProb = -1.;
|
---|
174 | fPed = -1.;
|
---|
175 | fPedRms = -1.;
|
---|
176 | fPedRmsErr = -1.;
|
---|
177 |
|
---|
178 | fNumHiGainSamples = -1.;
|
---|
179 | fNumLoGainSamples = -1.;
|
---|
180 |
|
---|
181 | fTimeFirstHiGain = 0 ;
|
---|
182 | fTimeLastHiGain = 0 ;
|
---|
183 | fTimeFirstLoGain = 0 ;
|
---|
184 | fTimeLastLoGain = 0 ;
|
---|
185 |
|
---|
186 | fAbsTimeMean = -1.;
|
---|
187 | fAbsTimeRms = -1.;
|
---|
188 |
|
---|
189 | fPheFFactorMethod = -1.;
|
---|
190 | fPheFFactorMethodErr = -1.;
|
---|
191 |
|
---|
192 | fMeanConversionFFactorMethod = -1.;
|
---|
193 | fMeanConversionBlindPixelMethod = -1.;
|
---|
194 | fMeanConversionPINDiodeMethod = -1.;
|
---|
195 | fMeanConversionCombinedMethod = -1.;
|
---|
196 |
|
---|
197 | fConversionFFactorMethodErr = -1.;
|
---|
198 | fConversionBlindPixelMethodErr = -1.;
|
---|
199 | fConversionPINDiodeMethodErr = -1.;
|
---|
200 | fConversionCombinedMethodErr = -1.;
|
---|
201 |
|
---|
202 | fSigmaConversionFFactorMethod = -1.;
|
---|
203 | fSigmaConversionBlindPixelMethod = -1.;
|
---|
204 | fSigmaConversionPINDiodeMethod = -1.;
|
---|
205 | fSigmaConversionCombinedMethod = -1.;
|
---|
206 |
|
---|
207 | fTotalFFactorFFactorMethod = -1.;
|
---|
208 | fTotalFFactorBlindPixelMethod = -1.;
|
---|
209 | fTotalFFactorPINDiodeMethod = -1.;
|
---|
210 | fTotalFFactorCombinedMethod = -1.;
|
---|
211 |
|
---|
212 | }
|
---|
213 |
|
---|
214 |
|
---|
215 | void MCalibrationPix::DefinePixId(Int_t i)
|
---|
216 | {
|
---|
217 |
|
---|
218 | fPixId = i;
|
---|
219 | fHist->ChangeHistId(i);
|
---|
220 |
|
---|
221 | }
|
---|
222 |
|
---|
223 |
|
---|
224 | // --------------------------------------------------------------------------
|
---|
225 | //
|
---|
226 | // Set the pedestals from outside
|
---|
227 | //
|
---|
228 | void MCalibrationPix::SetPedestal(const Float_t ped, const Float_t pedrms,
|
---|
229 | const Float_t higainsamp, const Float_t logainsamp )
|
---|
230 | {
|
---|
231 |
|
---|
232 | fPed = ped;
|
---|
233 | fPedRms = pedrms;
|
---|
234 |
|
---|
235 | fNumHiGainSamples = higainsamp;
|
---|
236 | fNumLoGainSamples = logainsamp;
|
---|
237 |
|
---|
238 | }
|
---|
239 |
|
---|
240 | // --------------------------------------------------------------------------
|
---|
241 | //
|
---|
242 | // Set the conversion factors from outside (only for MC)
|
---|
243 | //
|
---|
244 | void MCalibrationPix::SetConversionFFactorMethod(Float_t c, Float_t err, Float_t sig)
|
---|
245 | {
|
---|
246 | fMeanConversionFFactorMethod = c;
|
---|
247 | fConversionFFactorMethodErr = err;
|
---|
248 | fSigmaConversionFFactorMethod = sig;
|
---|
249 | }
|
---|
250 |
|
---|
251 | // --------------------------------------------------------------------------
|
---|
252 | //
|
---|
253 | // Set the conversion factors from outside (only for MC)
|
---|
254 | //
|
---|
255 | void MCalibrationPix::SetConversionCombinedMethod(Float_t c, Float_t err, Float_t sig)
|
---|
256 | {
|
---|
257 | fMeanConversionCombinedMethod = c;
|
---|
258 | fConversionCombinedMethodErr = err;
|
---|
259 | fSigmaConversionCombinedMethod = sig;
|
---|
260 | }
|
---|
261 |
|
---|
262 |
|
---|
263 | // --------------------------------------------------------------------------
|
---|
264 | //
|
---|
265 | // Set the conversion factors from outside (only for MC)
|
---|
266 | //
|
---|
267 | void MCalibrationPix::SetConversionBlindPixelMethod(Float_t c, Float_t err, Float_t sig)
|
---|
268 | {
|
---|
269 | fMeanConversionBlindPixelMethod = c;
|
---|
270 | fConversionBlindPixelMethodErr = err;
|
---|
271 | fSigmaConversionBlindPixelMethod = sig;
|
---|
272 | }
|
---|
273 |
|
---|
274 | // --------------------------------------------------------------------------
|
---|
275 | //
|
---|
276 | // Set the conversion factors from outside (only for MC)
|
---|
277 | //
|
---|
278 | void MCalibrationPix::SetConversionPINDiodeMethod(Float_t c, Float_t err, Float_t sig)
|
---|
279 | {
|
---|
280 | fMeanConversionPINDiodeMethod = c ;
|
---|
281 | fConversionPINDiodeMethodErr = err;
|
---|
282 | fSigmaConversionPINDiodeMethod = sig;
|
---|
283 | }
|
---|
284 |
|
---|
285 | // --------------------------------------------------------------------------
|
---|
286 | //
|
---|
287 | // Set the Hi Gain Saturation Bit from outside (only for MC)
|
---|
288 | //
|
---|
289 | void MCalibrationPix::SetHiGainSaturation(Bool_t b)
|
---|
290 | {
|
---|
291 |
|
---|
292 | if (b)
|
---|
293 | {
|
---|
294 | SETBIT(fFlags, kHiGainSaturation);
|
---|
295 | fHist->SetUseLoGain(1);
|
---|
296 | }
|
---|
297 | else
|
---|
298 | {
|
---|
299 | CLRBIT(fFlags, kHiGainSaturation);
|
---|
300 | fHist->SetUseLoGain(0);
|
---|
301 | }
|
---|
302 | }
|
---|
303 |
|
---|
304 | // --------------------------------------------------------------------------
|
---|
305 | //
|
---|
306 | // Set the Excluded Bit from outside
|
---|
307 | //
|
---|
308 | void MCalibrationPix::SetExcluded(Bool_t b )
|
---|
309 | {
|
---|
310 | b ? SETBIT(fFlags, kExcluded) : CLRBIT(fFlags, kExcluded);
|
---|
311 | }
|
---|
312 |
|
---|
313 |
|
---|
314 | // --------------------------------------------------------------------------
|
---|
315 | //
|
---|
316 | // Set the Excluded Bit from outside
|
---|
317 | //
|
---|
318 | void MCalibrationPix::SetExcludeQualityCheck(Bool_t b )
|
---|
319 | {
|
---|
320 | b ? SETBIT(fFlags, kExcludeQualityCheck) : CLRBIT(fFlags, kExcludeQualityCheck);
|
---|
321 | }
|
---|
322 |
|
---|
323 | // --------------------------------------------------------------------------
|
---|
324 | //
|
---|
325 | // Set the Excluded Bit from outside
|
---|
326 | //
|
---|
327 | void MCalibrationPix::SetChargeValid(Bool_t b )
|
---|
328 | {
|
---|
329 | b ? SETBIT(fFlags, kChargeValid) : CLRBIT(fFlags, kChargeValid);
|
---|
330 | }
|
---|
331 |
|
---|
332 | // --------------------------------------------------------------------------
|
---|
333 | //
|
---|
334 | // Set the Excluded Bit from outside
|
---|
335 | //
|
---|
336 | void MCalibrationPix::SetFitted(Bool_t b )
|
---|
337 | {
|
---|
338 | b ? SETBIT(fFlags, kFitted) : CLRBIT(fFlags, kFitted);
|
---|
339 | }
|
---|
340 |
|
---|
341 | // --------------------------------------------------------------------------
|
---|
342 | //
|
---|
343 | // Set the Excluded Bit from outside
|
---|
344 | //
|
---|
345 | void MCalibrationPix::SetOscillating(Bool_t b )
|
---|
346 | {
|
---|
347 | b ? SETBIT(fFlags, kOscillating) : CLRBIT(fFlags, kOscillating);
|
---|
348 | }
|
---|
349 |
|
---|
350 | // --------------------------------------------------------------------------
|
---|
351 | //
|
---|
352 | // Set the Excluded Bit from outside
|
---|
353 | //
|
---|
354 | void MCalibrationPix::SetBlindPixelMethodValid(Bool_t b )
|
---|
355 | {
|
---|
356 | b ? SETBIT(fFlags, kBlindPixelMethodValid) : CLRBIT(fFlags, kBlindPixelMethodValid);
|
---|
357 | }
|
---|
358 |
|
---|
359 | // --------------------------------------------------------------------------
|
---|
360 | //
|
---|
361 | // Set the Excluded Bit from outside
|
---|
362 | //
|
---|
363 | void MCalibrationPix::SetFFactorMethodValid(Bool_t b )
|
---|
364 | {
|
---|
365 | b ? SETBIT(fFlags, kFFactorMethodValid) : CLRBIT(fFlags, kFFactorMethodValid);
|
---|
366 | }
|
---|
367 |
|
---|
368 | // --------------------------------------------------------------------------
|
---|
369 | //
|
---|
370 | // Set the Excluded Bit from outside
|
---|
371 | //
|
---|
372 | void MCalibrationPix::SetPINDiodeMethodValid(Bool_t b )
|
---|
373 | {
|
---|
374 | b ? SETBIT(fFlags, kPINDiodeMethodValid) : CLRBIT(fFlags, kPINDiodeMethodValid);
|
---|
375 | }
|
---|
376 |
|
---|
377 | void MCalibrationPix::SetAbsTimeBordersHiGain(Byte_t f, Byte_t l)
|
---|
378 | {
|
---|
379 |
|
---|
380 | fTimeFirstHiGain = f;
|
---|
381 | fTimeLastHiGain = l;
|
---|
382 |
|
---|
383 | }
|
---|
384 |
|
---|
385 | void MCalibrationPix::SetAbsTimeBordersLoGain(Byte_t f, Byte_t l)
|
---|
386 | {
|
---|
387 |
|
---|
388 | fTimeFirstLoGain = f;
|
---|
389 | fTimeLastLoGain = l;
|
---|
390 |
|
---|
391 | }
|
---|
392 |
|
---|
393 | Bool_t MCalibrationPix::IsExcluded() const
|
---|
394 | {
|
---|
395 | return TESTBIT(fFlags,kExcluded);
|
---|
396 | }
|
---|
397 |
|
---|
398 | Bool_t MCalibrationPix::IsExcludeQualityCheck() const
|
---|
399 | {
|
---|
400 | return TESTBIT(fFlags,kExcludeQualityCheck);
|
---|
401 | }
|
---|
402 |
|
---|
403 | Bool_t MCalibrationPix::IsHiGainSaturation() const
|
---|
404 | {
|
---|
405 | return TESTBIT(fFlags,kHiGainSaturation);
|
---|
406 | }
|
---|
407 |
|
---|
408 | Bool_t MCalibrationPix::IsChargeValid() const
|
---|
409 | {
|
---|
410 | return TESTBIT(fFlags, kChargeValid);
|
---|
411 | }
|
---|
412 |
|
---|
413 | Bool_t MCalibrationPix::IsFitted() const
|
---|
414 | {
|
---|
415 | return TESTBIT(fFlags, kFitted);
|
---|
416 | }
|
---|
417 |
|
---|
418 | Bool_t MCalibrationPix::IsOscillating() const
|
---|
419 | {
|
---|
420 | return TESTBIT(fFlags, kOscillating);
|
---|
421 | }
|
---|
422 |
|
---|
423 | Bool_t MCalibrationPix::IsBlindPixelMethodValid() const
|
---|
424 | {
|
---|
425 | return TESTBIT(fFlags, kBlindPixelMethodValid);
|
---|
426 | }
|
---|
427 |
|
---|
428 | Bool_t MCalibrationPix::IsFFactorMethodValid() const
|
---|
429 | {
|
---|
430 | return TESTBIT(fFlags, kFFactorMethodValid);
|
---|
431 | }
|
---|
432 |
|
---|
433 | Bool_t MCalibrationPix::IsPINDiodeMethodValid() const
|
---|
434 | {
|
---|
435 | return TESTBIT(fFlags, kPINDiodeMethodValid);
|
---|
436 | }
|
---|
437 |
|
---|
438 | Bool_t MCalibrationPix::IsCombinedMethodValid() const
|
---|
439 | {
|
---|
440 | return TESTBIT(fFlags, kCombinedMethodValid);
|
---|
441 | }
|
---|
442 |
|
---|
443 |
|
---|
444 | // --------------------------------------------------------------------------
|
---|
445 | //
|
---|
446 | // 1) Return if the charge distribution is already succesfully fitted
|
---|
447 | // or if the histogram is empty
|
---|
448 | // 2) Set a lower Fit range according to 1.5 Pedestal RMS in order to avoid
|
---|
449 | // possible remaining cosmics to spoil the fit.
|
---|
450 | // 3) Decide if the LoGain Histogram is fitted or the HiGain Histogram
|
---|
451 | // 4) Fit the histograms with a Gaussian
|
---|
452 | // 5) In case of failure set the bit kFitted to false
|
---|
453 | // 6) Retrieve the results and store them in this class
|
---|
454 | // 7) Calculate the number of photo-electrons after the F-Factor method
|
---|
455 | // 8) Calculate the errors of the F-Factor method
|
---|
456 | //
|
---|
457 | // The fits are declared valid (fFitValid = kTRUE), if:
|
---|
458 | //
|
---|
459 | // 1) Pixel has a fitted charge greater than 3*PedRMS
|
---|
460 | // 2) Pixel has a fit error greater than 0.
|
---|
461 | // 3) Pixel has a fit Probability greater than 0.0001
|
---|
462 | // 4) Pixel has a charge sigma bigger than its Pedestal RMS
|
---|
463 | // 5) If FitTimes is used,
|
---|
464 | // the mean arrival time is at least 1.0 slices from the used edge slices
|
---|
465 | // (this stage is only performed in the times fit)
|
---|
466 | //
|
---|
467 | // If the histogram is empty, all values are set to -1.
|
---|
468 | //
|
---|
469 | // The conversion factor after the F-Factor method is declared valid, if:
|
---|
470 | //
|
---|
471 | // 1) fFitValid is kTRUE
|
---|
472 | // 2) Conversion Factor is bigger than 0.
|
---|
473 | // 3) The error of the conversion factor is smaller than 10%
|
---|
474 | //
|
---|
475 | Bool_t MCalibrationPix::FitCharge()
|
---|
476 | {
|
---|
477 |
|
---|
478 | //
|
---|
479 | // 1) Return if the charge distribution is already succesfully fitted
|
---|
480 | // or if the histogram is empty
|
---|
481 | //
|
---|
482 | if (fHist->IsChargeFitOK() || fHist->IsEmpty())
|
---|
483 | return kTRUE;
|
---|
484 |
|
---|
485 | //
|
---|
486 | // 2) Set a lower Fit range according to 1.5 Pedestal RMS in order to avoid
|
---|
487 | // possible remaining cosmics to spoil the fit.
|
---|
488 | //
|
---|
489 | // if (fPed && fPedRms)
|
---|
490 | // fHist->SetLowerFitRange(1.5*fPedRms);
|
---|
491 | // else
|
---|
492 | // *fLog << warn << "WARNING: Cannot set lower fit range: Pedestals not available" << endl;
|
---|
493 |
|
---|
494 | //
|
---|
495 | // 3) Decide if the LoGain Histogram is fitted or the HiGain Histogram
|
---|
496 | //
|
---|
497 | if (fHist->UseLoGain())
|
---|
498 | SetHiGainSaturation();
|
---|
499 |
|
---|
500 | //
|
---|
501 | // 4) Fit the Lo Gain histograms with a Gaussian
|
---|
502 | //
|
---|
503 | if (fHist->FitCharge())
|
---|
504 | SETBIT(fFlags,kFitted);
|
---|
505 | else
|
---|
506 | {
|
---|
507 | *fLog << warn << "WARNING: Could not fit charges of pixel " << fPixId << endl;
|
---|
508 | //
|
---|
509 | // 5) In case of failure set the bit kFitted to false
|
---|
510 | //
|
---|
511 | CLRBIT(fFlags,kFitted);
|
---|
512 | }
|
---|
513 |
|
---|
514 | //
|
---|
515 | // 6) Retrieve the results and store them in this class
|
---|
516 | // If fFitted is false, we get the eans and RMS of the histogram!!
|
---|
517 | //
|
---|
518 | fCharge = fHist->GetChargeMean();
|
---|
519 | fChargeErr = fHist->GetChargeMeanErr();
|
---|
520 | fSigmaCharge = fHist->GetChargeSigma();
|
---|
521 | fSigmaChargeErr = fHist->GetChargeSigmaErr();
|
---|
522 | fChargeProb = fHist->GetChargeProb();
|
---|
523 |
|
---|
524 |
|
---|
525 | fAbsTimeMean = fHist->GetAbsTimeMean();
|
---|
526 | fAbsTimeMeanErr = fHist->GetAbsTimeMeanErr();
|
---|
527 | fAbsTimeRms = fHist->GetAbsTimeRms();
|
---|
528 |
|
---|
529 | if (CheckTimeFitValidity())
|
---|
530 | SETBIT(fFlags,kTimeFitValid);
|
---|
531 | else
|
---|
532 | CLRBIT(fFlags,kTimeFitValid);
|
---|
533 |
|
---|
534 | //
|
---|
535 | //Calculate the conversion factors
|
---|
536 | //
|
---|
537 | if (IsHiGainSaturation())
|
---|
538 | ApplyLoGainConversion();
|
---|
539 |
|
---|
540 | if (CheckChargeValidity())
|
---|
541 | SETBIT(fFlags,kChargeValid);
|
---|
542 | else
|
---|
543 | {
|
---|
544 | CLRBIT(fFlags,kChargeValid);
|
---|
545 | return kFALSE;
|
---|
546 | }
|
---|
547 |
|
---|
548 | return kTRUE;
|
---|
549 |
|
---|
550 | }
|
---|
551 |
|
---|
552 | //
|
---|
553 | // Calculate the number of photo-electrons after the F-Factor method
|
---|
554 | // Calculate the errors of the F-Factor method
|
---|
555 | //
|
---|
556 | Bool_t MCalibrationPix::CalcFFactorMethod()
|
---|
557 | {
|
---|
558 |
|
---|
559 | if ( (fCharge == -1.)
|
---|
560 | || (fChargeErr < 0.)
|
---|
561 | || (fSigmaCharge < 0.)
|
---|
562 | || (fPedRms < 0.) )
|
---|
563 | {
|
---|
564 | *fLog << warn << GetDescriptor() << "Cannot calculate the FFactor Method! "
|
---|
565 | << "Some of the needed parameters are not available ";
|
---|
566 | CLRBIT(fFlags,kFFactorMethodValid);
|
---|
567 | return kFALSE;
|
---|
568 | }
|
---|
569 |
|
---|
570 | //
|
---|
571 | // Square all variables in order to avoid applications of square root
|
---|
572 | //
|
---|
573 | // First the relative error squares
|
---|
574 | //
|
---|
575 | const Float_t chargeSquare = fCharge * fCharge;
|
---|
576 | const Float_t chargeSquareRelErrSquare = 4.* fChargeErr * fChargeErr / chargeSquare;
|
---|
577 |
|
---|
578 | const Float_t chargeRelErrSquare = fChargeErr * fChargeErr
|
---|
579 | / (fCharge * fCharge);
|
---|
580 |
|
---|
581 | const Float_t ffactorsquare = gkFFactor * gkFFactor;
|
---|
582 | const Float_t ffactorsquareRelErrSquare = 4.*gkFFactorErr * gkFFactorErr / ffactorsquare;
|
---|
583 |
|
---|
584 | const Float_t avQERelErrSquare = fAverageQEErr * fAverageQEErr / fAverageQE / fAverageQE;
|
---|
585 |
|
---|
586 | const Float_t avQEFFactor = TMath::Sqrt( ( 1. - fAverageQE ) / fAverageQE );
|
---|
587 | const Float_t avQEFFactorErr = 1./ ( 2. * avQEFFactor ) * fAverageQEErr
|
---|
588 | / ( fAverageQE * fAverageQE );
|
---|
589 | const Float_t avQEFFactorRelErrSquare = avQEFFactorErr * avQEFFactorErr
|
---|
590 | / ( avQEFFactor * avQEFFactor) ;
|
---|
591 | //
|
---|
592 | // Now the absolute error squares
|
---|
593 | //
|
---|
594 | const Float_t sigmaSquare = fSigmaCharge * fSigmaCharge;
|
---|
595 | const Float_t sigmaSquareErrSquare = 4.*fSigmaChargeErr* fSigmaChargeErr * sigmaSquare;
|
---|
596 |
|
---|
597 | Float_t pedRmsSquare = fPedRms * fPedRms;
|
---|
598 | Float_t pedRmsSquareErrSquare = 4.*fPedRmsErr * fPedRmsErr * pedRmsSquare;
|
---|
599 |
|
---|
600 | if (!IsHiGainSaturation())
|
---|
601 | { /* HiGain */
|
---|
602 |
|
---|
603 | pedRmsSquare *= fNumHiGainSamples;
|
---|
604 | pedRmsSquareErrSquare *= fNumHiGainSamples*fNumHiGainSamples;
|
---|
605 | }
|
---|
606 | else
|
---|
607 | { /* LoGain */
|
---|
608 |
|
---|
609 | //
|
---|
610 | // We do not know the Lo Gain Pedestal RMS, so we have to retrieve it
|
---|
611 | // from the HI GAIN (all calculation per slice up to now):
|
---|
612 | //
|
---|
613 | // We extract the pure NSB contribution:
|
---|
614 | //
|
---|
615 | const Float_t elecRmsSquare = gkElectronicPedRms * gkElectronicPedRms;
|
---|
616 | const Float_t elecRmsSquareErrSquare = 4.*gkElectronicPedRmsErr * gkElectronicPedRmsErr * elecRmsSquare;
|
---|
617 |
|
---|
618 | Float_t nsbSquare = pedRmsSquare - elecRmsSquare;
|
---|
619 | Float_t nsbSquareRelErrSquare = (pedRmsSquareErrSquare + elecRmsSquareErrSquare)
|
---|
620 | / (nsbSquare * nsbSquare) ;
|
---|
621 |
|
---|
622 | if (nsbSquare < 0.)
|
---|
623 | nsbSquare = 0.;
|
---|
624 |
|
---|
625 | //
|
---|
626 | // Now, we divide the NSB by the conversion factor and
|
---|
627 | // add it quadratically to the electronic noise
|
---|
628 | //
|
---|
629 | const Float_t conversionSquare = fConversionHiLo * fConversionHiLo;
|
---|
630 | const Float_t conversionSquareRelErrSquare = 4.*fConversionHiLoErr * fConversionHiLoErr / conversionSquare;
|
---|
631 |
|
---|
632 | const Float_t convertedNsbSquare = nsbSquare / conversionSquare;
|
---|
633 | const Float_t convertedNsbSquareErrSquare = (nsbSquareRelErrSquare + conversionSquareRelErrSquare)
|
---|
634 | * convertedNsbSquare * convertedNsbSquare;
|
---|
635 |
|
---|
636 | pedRmsSquare = convertedNsbSquare + elecRmsSquare;
|
---|
637 | pedRmsSquareErrSquare = convertedNsbSquareErrSquare + elecRmsSquareErrSquare;
|
---|
638 |
|
---|
639 | //
|
---|
640 | // Now, correct for the number of used FADC slices in the LoGain:
|
---|
641 | //
|
---|
642 | pedRmsSquare *= fNumLoGainSamples;
|
---|
643 | pedRmsSquareErrSquare *= fNumLoGainSamples*fNumLoGainSamples;
|
---|
644 | //
|
---|
645 | // Correct also for the conversion to Hi-Gain:
|
---|
646 | //
|
---|
647 | pedRmsSquare *= fConversionHiLo*fConversionHiLo;
|
---|
648 | pedRmsSquareErrSquare *= fConversionHiLo*fConversionHiLo*fConversionHiLo*fConversionHiLo;
|
---|
649 |
|
---|
650 | } /* if (HiGainSaturation) */
|
---|
651 |
|
---|
652 | //
|
---|
653 | // Calculate the reduced sigmas
|
---|
654 | //
|
---|
655 | const Float_t rsigmachargesquare = sigmaSquare - pedRmsSquare;
|
---|
656 | if (rsigmachargesquare <= 0.)
|
---|
657 | {
|
---|
658 | *fLog << warn
|
---|
659 | << "WARNING: Cannot apply F-Factor calibration: Reduced Sigma smaller than 0 in pixel "
|
---|
660 | << fPixId << endl;
|
---|
661 | CLRBIT(fFlags,kFFactorMethodValid);
|
---|
662 | return kFALSE;
|
---|
663 | }
|
---|
664 |
|
---|
665 | const Float_t rSigmaSquareRelErrSquare = (sigmaSquareErrSquare + pedRmsSquareErrSquare)
|
---|
666 | / (rsigmachargesquare * rsigmachargesquare) ;
|
---|
667 |
|
---|
668 | fRSigmaCharge = TMath::Sqrt(rsigmachargesquare);
|
---|
669 | fRSigmaChargeErr = TMath::Sqrt(sigmaSquareErrSquare + pedRmsSquareErrSquare);
|
---|
670 |
|
---|
671 |
|
---|
672 | //
|
---|
673 | // Calculate the number of phe's from the F-Factor method
|
---|
674 | // (independent on Hi Gain or Lo Gain)
|
---|
675 | //
|
---|
676 | fPheFFactorMethod = ffactorsquare * chargeSquare / rsigmachargesquare;
|
---|
677 | //
|
---|
678 | // Calculate the number of photons from the F-Factor method
|
---|
679 | // FIXME: This is a preliminary solution, the qe shall be
|
---|
680 | // calibrated itself!
|
---|
681 | //
|
---|
682 | fPheFFactorMethod /= fAverageQE;
|
---|
683 |
|
---|
684 | const Float_t pheFFactorRelErrSquare = ffactorsquareRelErrSquare
|
---|
685 | + chargeSquareRelErrSquare
|
---|
686 | + rSigmaSquareRelErrSquare
|
---|
687 | + avQERelErrSquare;
|
---|
688 |
|
---|
689 | fPheFFactorMethodErr = TMath::Sqrt(pheFFactorRelErrSquare) * fPheFFactorMethod;
|
---|
690 |
|
---|
691 | fMeanConversionFFactorMethod = fPheFFactorMethod / fCharge ;
|
---|
692 | fConversionFFactorMethodErr = ( pheFFactorRelErrSquare + chargeRelErrSquare )
|
---|
693 | * fMeanConversionFFactorMethod * fMeanConversionFFactorMethod;
|
---|
694 |
|
---|
695 | const Float_t convrelerror = fConversionFFactorMethodErr
|
---|
696 | / fMeanConversionFFactorMethod;
|
---|
697 |
|
---|
698 | if ( (fMeanConversionFFactorMethod > 0.) && (convrelerror < gkConvFFactorRelErrLimit))
|
---|
699 | SETBIT(fFlags,kFFactorMethodValid);
|
---|
700 |
|
---|
701 | fSigmaConversionFFactorMethod = GetTotalFFactorFFactorMethod()*TMath::Sqrt(fMeanConversionFFactorMethod);
|
---|
702 |
|
---|
703 | //
|
---|
704 | // Calculate the Total F-Factor of the camera ( in photons )
|
---|
705 | //
|
---|
706 | if (fPheFFactorMethod > 0)
|
---|
707 | {
|
---|
708 | fTotalFFactorFFactorMethod = (fRSigmaCharge/fCharge)*TMath::Sqrt(fPheFFactorMethod);
|
---|
709 | fTotalFFactorFFactorMethod *= avQEFFactor;
|
---|
710 | }
|
---|
711 |
|
---|
712 | //
|
---|
713 | // Calculate the error of the Total F-Factor of the camera ( in photons )
|
---|
714 | //
|
---|
715 | const Float_t rSigmaChargeRelErrSquare = fRSigmaChargeErr * fRSigmaChargeErr
|
---|
716 | / (fRSigmaCharge * fRSigmaCharge) ;
|
---|
717 |
|
---|
718 | fTotalFFactorErrFFactorMethod = TMath::Sqrt( rSigmaChargeRelErrSquare
|
---|
719 | + chargeRelErrSquare
|
---|
720 | + pheFFactorRelErrSquare
|
---|
721 | + avQEFFactorRelErrSquare );
|
---|
722 |
|
---|
723 | fTotalFFactorErrFFactorMethod *= fTotalFFactorFFactorMethod;
|
---|
724 |
|
---|
725 | return kTRUE;
|
---|
726 | }
|
---|
727 |
|
---|
728 |
|
---|
729 | //
|
---|
730 | // The check returns kTRUE if:
|
---|
731 | //
|
---|
732 | // 0) Pixel has BIT fitted set:
|
---|
733 | // This means:
|
---|
734 | // a) No result is a nan
|
---|
735 | // b) The NDF is not smaller than fNDFLimit (5)
|
---|
736 | // c) The Probability is greater than gkProbLimit (default 0.001 == 99.9%)
|
---|
737 | // 1) Pixel has a fitted charge greater than 3*PedRMS
|
---|
738 | // 2) Pixel has a fit error greater than 0.
|
---|
739 | // 3) Pixel has a fitted charge greater its charge error
|
---|
740 | // 4) Pixel has a fit Probability greater than 0.0001
|
---|
741 | // 5) Pixel has a charge sigma bigger than its Pedestal RMS
|
---|
742 | //
|
---|
743 | Bool_t MCalibrationPix::CheckChargeValidity()
|
---|
744 | {
|
---|
745 |
|
---|
746 | if (!IsFitted())
|
---|
747 | return kFALSE;
|
---|
748 |
|
---|
749 | if (IsExcludeQualityCheck())
|
---|
750 | return kTRUE;
|
---|
751 |
|
---|
752 | Float_t pedestal;
|
---|
753 |
|
---|
754 | if (!IsHiGainSaturation()) /* higain */
|
---|
755 | pedestal = GetPedRms()*TMath::Sqrt(fNumHiGainSamples);
|
---|
756 | else /* logain */
|
---|
757 | pedestal = GetPedRms()*TMath::Sqrt(fNumLoGainSamples);
|
---|
758 |
|
---|
759 |
|
---|
760 | if (fCharge < gkChargeLimit*pedestal)
|
---|
761 | {
|
---|
762 | *fLog << warn << "WARNING: Fitted Charge is smaller than "
|
---|
763 | << gkChargeLimit << " Pedestal RMS in Pixel " << fPixId << endl;
|
---|
764 | return kFALSE;
|
---|
765 | }
|
---|
766 |
|
---|
767 | if (fChargeErr < gkChargeErrLimit)
|
---|
768 | {
|
---|
769 | *fLog << warn << "WARNING: Err of Fitted Charge is smaller than "
|
---|
770 | << gkChargeErrLimit << " in Pixel " << fPixId << endl;
|
---|
771 | return kFALSE;
|
---|
772 | }
|
---|
773 |
|
---|
774 | if (fCharge < gkChargeRelErrLimit*fChargeErr)
|
---|
775 | {
|
---|
776 | *fLog << warn << "WARNING: Fitted Charge is smaller than "
|
---|
777 | << gkChargeRelErrLimit << "* its error in Pixel " << fPixId << endl;
|
---|
778 | return kFALSE;
|
---|
779 | }
|
---|
780 |
|
---|
781 | if (!fHist->IsChargeFitOK())
|
---|
782 | {
|
---|
783 | *fLog << warn << "WARNING: Probability of Fitted Charge too low in Pixel "
|
---|
784 | << fPixId << endl;
|
---|
785 | return kFALSE;
|
---|
786 | }
|
---|
787 |
|
---|
788 | if (fSigmaCharge < pedestal)
|
---|
789 | {
|
---|
790 | *fLog << warn << "WARNING: Sigma of Fitted Charge smaller than Pedestal RMS in Pixel "
|
---|
791 | << fPixId << endl;
|
---|
792 | return kFALSE;
|
---|
793 | }
|
---|
794 | return kTRUE;
|
---|
795 | }
|
---|
796 |
|
---|
797 | //
|
---|
798 | // The check return kTRUE if:
|
---|
799 | //
|
---|
800 | // 0) No value is nan
|
---|
801 | // 1) Pixel has a fitted rel. time smaller than 3*FADC slices
|
---|
802 | // 2) Pixel has a fit error greater than 0.
|
---|
803 | // 4) Pixel has a fit Probability greater than 0.001
|
---|
804 | // 5) The absolute arrival time is at least 1.0 slices from the used edge slices
|
---|
805 | //
|
---|
806 | Bool_t MCalibrationPix::CheckTimeFitValidity()
|
---|
807 | {
|
---|
808 |
|
---|
809 |
|
---|
810 | if (IsExcludeQualityCheck())
|
---|
811 | return kTRUE;
|
---|
812 |
|
---|
813 | if (IsHiGainSaturation())
|
---|
814 | {
|
---|
815 |
|
---|
816 | if (fAbsTimeMean < (Float_t)fTimeFirstLoGain+1)
|
---|
817 | {
|
---|
818 | *fLog << warn
|
---|
819 | << "WARNING: Some absolute times smaller than limit in Pixel "
|
---|
820 | << fPixId << " time: " << fAbsTimeMean
|
---|
821 | << " Limit: " << (Float_t)fTimeFirstLoGain+1. << endl;
|
---|
822 | return kFALSE;
|
---|
823 | }
|
---|
824 |
|
---|
825 | if (fAbsTimeMean > (Float_t)fTimeLastLoGain-1)
|
---|
826 | {
|
---|
827 | *fLog << warn
|
---|
828 | << "WARNING: Some absolute times bigger than limit in Pixel "
|
---|
829 | << fPixId << " time: " << fAbsTimeMean
|
---|
830 | << " Limit: " << (Float_t)fTimeLastLoGain-1. << endl;
|
---|
831 | return kFALSE;
|
---|
832 | }
|
---|
833 |
|
---|
834 | }
|
---|
835 | else
|
---|
836 | {
|
---|
837 |
|
---|
838 | if (fAbsTimeMean < (Float_t)fTimeFirstHiGain+1.)
|
---|
839 | {
|
---|
840 | *fLog << warn
|
---|
841 | << "WARNING: Some absolute times smaller than limit in Pixel "
|
---|
842 | << fPixId << " time: " << fAbsTimeMean
|
---|
843 | << " Limit: " << (Float_t)fTimeFirstHiGain+1. << endl;
|
---|
844 | // return kFALSE;
|
---|
845 | }
|
---|
846 |
|
---|
847 | if (fAbsTimeMean > (Float_t)fTimeLastHiGain-1.)
|
---|
848 | {
|
---|
849 | *fLog << warn
|
---|
850 | << "WARNING: Some absolute times bigger than limit in Pixel "
|
---|
851 | << fPixId << " time: " << fAbsTimeMean
|
---|
852 | << " Limit: " << (Float_t)fTimeLastHiGain-1. << endl;
|
---|
853 | // return kFALSE;
|
---|
854 | }
|
---|
855 |
|
---|
856 | }
|
---|
857 |
|
---|
858 |
|
---|
859 |
|
---|
860 | return kTRUE;
|
---|
861 | }
|
---|
862 |
|
---|
863 |
|
---|
864 | void MCalibrationPix::CheckOscillations()
|
---|
865 | {
|
---|
866 | fHist->CheckOscillations();
|
---|
867 | }
|
---|
868 |
|
---|
869 | void MCalibrationPix::ApplyLoGainConversion()
|
---|
870 | {
|
---|
871 |
|
---|
872 | const Float_t chargeRelErrSquare = fChargeErr * fChargeErr
|
---|
873 | /( fCharge * fCharge );
|
---|
874 | const Float_t sigmaRelErrSquare = fSigmaChargeErr * fSigmaChargeErr
|
---|
875 | /( fSigmaCharge * fSigmaCharge );
|
---|
876 | const Float_t conversionRelErrSquare = fConversionHiLoErr * fConversionHiLoErr
|
---|
877 | /( fConversionHiLo * fConversionHiLo );
|
---|
878 |
|
---|
879 | fCharge *= fConversionHiLo;
|
---|
880 | fChargeErr = TMath::Sqrt(chargeRelErrSquare + conversionRelErrSquare) * fCharge;
|
---|
881 |
|
---|
882 | fSigmaCharge *= fConversionHiLo;
|
---|
883 | fSigmaChargeErr = TMath::Sqrt(sigmaRelErrSquare + conversionRelErrSquare) * fSigmaCharge;
|
---|
884 |
|
---|
885 | }
|
---|