source: trunk/MagicSoft/Mars/mcalib/MCalibrationPix.cc@ 2840

Last change on this file since 2840 was 2839, checked in by gaug, 21 years ago
*** empty log message ***
File size: 16.3 KB
Line 
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/////////////////////////////////////////////////////////////////////////////
33#include "MCalibrationPix.h"
34#include "MCalibrationConfig.h"
35
36#include "MLog.h"
37#include "MLogManip.h"
38
39ClassImp(MCalibrationPix);
40
41using namespace std;
42
43// --------------------------------------------------------------------------
44//
45// Default Constructor:
46//
47// The following values are initialized to meaningful values:
48//
49// - The Electronic Rms to 1.5 per FADC slice
50// - The uncertainty about the Electronic RMS to 0.3 per slice
51// - The F-Factor is assumed to have been measured in Munich to 1.13 - 1.17.
52// We use here the Square of the Munich definition, thus:
53// Mean F-Factor = 1.15*1.15 = 1.32
54// Error F-Factor = 2.*0.02 = 0.04
55//
56MCalibrationPix::MCalibrationPix(const char *name, const char *title)
57 : fPixId(-1),
58 fCharge(-1.),
59 fErrCharge(-1.),
60 fSigmaCharge(-1.),
61 fErrSigmaCharge(-1.),
62 fRSigmaSquare(-1.),
63 fChargeProb(-1.),
64 fPed(-1.),
65 fPedRms(-1.),
66 fErrPedRms(0.),
67 fElectronicPedRms(1.5),
68 fErrElectronicPedRms(0.3),
69 fTime(-1.),
70 fSigmaTime(-1.),
71 fTimeChiSquare(-1.),
72 fFactor(1.32),
73 fFactorError(0.04),
74 fPheFFactorMethod(-1.),
75 fPheFFactorMethodError(-1.),
76 fConversionFFactorMethod(-1.),
77 fConversionBlindPixelMethod(-1.),
78 fConversionPINDiodeMethod(-1.),
79 fConversionErrorFFactorMethod(-1.),
80 fConversionErrorBlindPixelMethod(-1.),
81 fConversionErrorPINDiodeMethod(-1.),
82 fConversionSigmaFFactorMethod(-1.),
83 fConversionSigmaBlindPixelMethod(-1.),
84 fConversionSigmaPINDiodeMethod(-1.),
85 fHiGainSaturation(kFALSE),
86 fExcluded(kFALSE),
87 fFitValid(kFALSE),
88 fFitted(kFALSE),
89 fBlindPixelMethodValid(kFALSE),
90 fFFactorMethodValid(kFALSE),
91 fPINDiodeMethodValid(kFALSE)
92{
93
94 fName = name ? name : "MCalibrationPixel";
95 fTitle = title ? title : "Container of the MHCalibrationPixels and the fit results";
96
97 //
98 // At the moment, we don't have a database, yet,
99 // so we get it from the configuration file
100 //
101 fConversionHiLo = gkConversionHiLo;
102 fConversionHiLoError = gkConversionHiLoError;
103
104 fHist = new MHCalibrationPixel("MHCalibrationPixel","Calibration Histograms Pixel ");
105
106 if (!fHist)
107 *fLog << warn << dbginf << " Could not create MHCalibrationPixel " << endl;
108
109}
110
111MCalibrationPix::~MCalibrationPix()
112{
113 delete fHist;
114}
115
116
117void MCalibrationPix::DefinePixId(Int_t i)
118{
119
120 fPixId = i;
121 fHist->ChangeHistId(i);
122
123}
124
125
126// ------------------------------------------------------------------------
127//
128// Invalidate values
129//
130void MCalibrationPix::Clear(Option_t *o)
131{
132 fHist->Reset();
133}
134
135
136// --------------------------------------------------------------------------
137//
138// 1) Return if the charge distribution is already succesfully fitted
139// or if the histogram is empty
140// 2) Set a lower Fit range according to 1.5 Pedestal RMS in order to avoid
141// possible remaining cosmics to spoil the fit.
142// 3) Decide if the LoGain Histogram is fitted or the HiGain Histogram
143// 4) Fit the histograms with a Gaussian
144// 5) In case of failure print out the fit results
145// 6) Retrieve the results and store them in this class
146// 7) Calculate the number of photo-electrons after the F-Factor method
147// 8) Calculate the errors of the F-Factor method
148//
149// The fits are declared valid (fFitValid = kTRUE), if:
150//
151// 1) Pixel has a fitted charge greater than 5*PedRMS
152// 2) Pixel has a fit error greater than 0.
153// 3) Pixel has a fit Probability greater than 0.0001
154// 4) Pixel has a charge sigma bigger than its Pedestal RMS
155// 5) If FitTimes is used,
156// the mean arrival time is at least 1.0 slices from the used edge slices
157// (this stage is only performed in the times fit)
158//
159// If the histogram is empty, all values are set to -1.
160//
161// The conversion factor after the F-Factor method is declared valid, if:
162//
163// 1) fFitValid is kTRUE
164// 2) Conversion Factor is bigger than 0.
165// 3) The error of the conversion factor is smaller than 10%
166//
167Bool_t MCalibrationPix::FitCharge()
168{
169
170 //
171 // 1) Return if the charge distribution is already succesfully fitted
172 // or if the histogram is empty
173 //
174 if (fHist->IsFitOK() || fHist->IsEmpty())
175 return kTRUE;
176
177 //
178 // 2) Set a lower Fit range according to 1.5 Pedestal RMS in order to avoid
179 // possible remaining cosmics to spoil the fit.
180 //
181 if (fPed && fPedRms)
182 fHist->SetLowerFitRange(1.5*fPedRms);
183 else
184 *fLog << warn << "WARNING: Cannot set lower fit range: Pedestals not available" << endl;
185
186 //
187 // 3) Decide if the LoGain Histogram is fitted or the HiGain Histogram
188 //
189 if (fHist->UseLoGain())
190 {
191
192 SetHiGainSaturation();
193
194 //
195 // 4) Fit the Lo Gain histograms with a Gaussian
196 //
197 if(!fHist->FitChargeLoGain())
198 {
199 *fLog << warn << "WARNING: Could not fit Lo Gain charges of pixel " << fPixId << endl;
200 //
201 // 5) In case of failure print out the fit results
202 //
203 fHist->PrintChargeFitResult();
204 }
205 }
206 else
207 {
208 //
209 // 4) Fit the Hi Gain histograms with a Gaussian
210 //
211 if(!fHist->FitChargeHiGain())
212 {
213 *fLog << warn << "WARNING: Could not fit Hi Gain charges of pixel " << fPixId << endl;
214 //
215 // 5) In case of failure print out the fit results
216 //
217 fHist->PrintChargeFitResult();
218 }
219 }
220
221
222 //
223 // 6) Retrieve the results and store them in this class
224 //
225 fCharge = fHist->GetChargeMean();
226 fErrCharge = fHist->GetChargeMeanErr();
227 fSigmaCharge = fHist->GetChargeSigma();
228 fErrSigmaCharge = fHist->GetChargeSigmaErr();
229 fChargeProb = fHist->GetChargeProb();
230
231 if (fCharge <= 0.)
232 {
233 *fLog << warn << "WARNING: Cannot apply calibration: Mean Fitted Charges are smaller than 0 in pixel "
234 << fPixId << endl;
235 return kFALSE;
236 }
237
238 if (fErrCharge > 0.)
239 fFitted = kTRUE;
240
241 if (CheckChargeFitValidity())
242 fFitValid = kTRUE;
243
244
245 //
246 // 7) Calculate the number of photo-electrons after the F-Factor method
247 // 8) Calculate the errors of the F-Factor method
248 //
249 if ((fPed > 0.) && (fPedRms > 0.))
250 {
251
252 //
253 // Square all variables in order to avoid applications of square root
254 //
255 // First the relative error squares
256 //
257 const Float_t chargeSquare = fCharge* fCharge;
258 const Float_t chargeSquareRelErrSquare = 4.*fErrCharge*fErrCharge / chargeSquare;
259
260 const Float_t fFactorRelErrSquare = fFactorError * fFactorError / (fFactor * fFactor);
261 //
262 // Now the absolute error squares
263 //
264 const Float_t sigmaSquare = fSigmaCharge* fSigmaCharge;
265 const Float_t sigmaSquareErrSquare = 4.*fErrSigmaCharge*fErrSigmaCharge * sigmaSquare;
266
267 const Float_t elecRmsSquare = fElectronicPedRms* fElectronicPedRms;
268 const Float_t elecRmsSquareErrSquare = 4.*fErrElectronicPedRms*fErrElectronicPedRms * elecRmsSquare;
269
270 Float_t pedRmsSquare = fPedRms* fPedRms;
271 Float_t pedRmsSquareErrSquare = 4.*fErrPedRms*fErrPedRms * pedRmsSquare;
272
273 if (fHiGainSaturation)
274 {
275
276 //
277 // We do not know the Lo Gain Pedestal RMS, so we have to retrieve it
278 // from the Hi Gain:
279 //
280 // We extract the pure NSB contribution:
281 //
282 Float_t nsbSquare = pedRmsSquare - elecRmsSquare;
283 Float_t nsbSquareRelErrSquare = (pedRmsSquareErrSquare + elecRmsSquareErrSquare)
284 / (nsbSquare * nsbSquare) ;
285
286 if (nsbSquare < 0.)
287 nsbSquare = 0.;
288
289 //
290 // Now, we divide the NSB by the conversion factor and
291 // add it quadratically to the electronic noise
292 //
293 const Float_t conversionSquare = fConversionHiLo *fConversionHiLo;
294 const Float_t conversionSquareRelErrSquare = 4.*fConversionHiLoError*fConversionHiLoError/conversionSquare;
295
296 //
297 // Calculate the new "Pedestal RMS"
298 //
299 const Float_t convertedNsbSquare = nsbSquare / conversionSquare;
300 const Float_t convertedNsbSquareErrSquare = (nsbSquareRelErrSquare + conversionSquareRelErrSquare)
301 * convertedNsbSquare * convertedNsbSquare;
302
303 pedRmsSquare = convertedNsbSquare + elecRmsSquare;
304 pedRmsSquareErrSquare = convertedNsbSquareErrSquare + elecRmsSquareErrSquare;
305
306 } /* if (fHiGainSaturation) */
307
308 //
309 // Calculate the reduced sigmas
310 //
311 fRSigmaSquare = sigmaSquare - pedRmsSquare;
312 if (fRSigmaSquare <= 0.)
313 {
314 *fLog << warn
315 << "WARNING: Cannot apply F-Factor calibration: Reduced Sigma smaller than 0 in pixel "
316 << fPixId << endl;
317 if (fHiGainSaturation)
318 ApplyLoGainConversion();
319 return kFALSE;
320 }
321
322 const Float_t rSigmaSquareRelErrSquare = (sigmaSquareErrSquare + pedRmsSquareErrSquare)
323 / (fRSigmaSquare * fRSigmaSquare) ;
324
325 //
326 // Calculate the number of phe's from the F-Factor method
327 // (independent on Hi Gain or Lo Gain)
328 //
329 fPheFFactorMethod = fFactor * chargeSquare / fRSigmaSquare;
330
331 const Float_t pheFFactorRelErrSquare = fFactorRelErrSquare
332 + chargeSquareRelErrSquare
333 + rSigmaSquareRelErrSquare ;
334
335 fPheFFactorMethodError = TMath::Sqrt(pheFFactorRelErrSquare) * fPheFFactorMethod;
336
337 //
338 // Calculate the conversion factors
339 //
340 if (fHiGainSaturation)
341 ApplyLoGainConversion();
342
343 const Float_t chargeRelErrSquare = fErrCharge*fErrCharge / (fCharge * fCharge);
344
345 fConversionFFactorMethod = fPheFFactorMethod / fCharge ;
346 fConversionErrorFFactorMethod = ( pheFFactorRelErrSquare + chargeRelErrSquare )
347 * fConversionFFactorMethod * fConversionFFactorMethod;
348
349 if ( IsFitValid() &&
350 (fConversionFFactorMethod > 0.) &&
351 (fConversionErrorFFactorMethod/fConversionFFactorMethod < 0.1) )
352 fFFactorMethodValid = kTRUE;
353
354
355 } /* if ((fPed > 0.) && (fPedRms > 0.)) */
356
357 return kTRUE;
358
359}
360
361//
362// The check return kTRUE if:
363//
364// 1) Pixel has a fitted charge greater than 5*PedRMS
365// 2) Pixel has a fit error greater than 0.
366// 3) Pixel has a fit Probability greater than 0.0001
367// 4) Pixel has a charge sigma bigger than its Pedestal RMS
368//
369Bool_t MCalibrationPix::CheckChargeFitValidity()
370{
371
372 Float_t equivpedestal = GetPedRms();
373
374 if (fHiGainSaturation)
375 equivpedestal /= fConversionHiLo;
376
377 if (fCharge < 5.*equivpedestal)
378 {
379 *fLog << warn << "WARNING: Fitted Charge is smaller than 5 Pedestal RMS in Pixel " << fPixId << endl;
380 return kFALSE;
381 }
382
383 if (fErrCharge < 0.)
384 {
385 *fLog << warn << "WARNING: Error of Fitted Charge is smaller than 0 in Pixel " << fPixId << endl;
386 return kFALSE;
387 }
388
389 if (!fHist->IsFitOK())
390 {
391 *fLog << warn << "WARNING: Probability of Fitted Charge too low in Pixel " << fPixId << endl;
392 return kFALSE;
393 }
394
395 if (fSigmaCharge < equivpedestal)
396 {
397 *fLog << warn << "WARNING: Sigma of Fitted Charge smaller than Pedestal RMS in Pixel " << fPixId << endl;
398 return kFALSE;
399 }
400 return kTRUE;
401}
402
403//
404// The check returns kTRUE if:
405//
406// The mean arrival time is at least 1.0 slices from the used edge slices
407//
408Bool_t MCalibrationPix::CheckTimeFitValidity()
409{
410
411 Float_t lowerrange;
412 Float_t upperrange;
413
414 if (fHiGainSaturation)
415 {
416 lowerrange = (Float_t)fHist->GetTimeLowerFitRangeLoGain()+1.;
417 upperrange = (Float_t)fHist->GetTimeUpperFitRangeLoGain()+1.;
418 }
419 else
420 {
421 lowerrange = (Float_t)fHist->GetTimeLowerFitRangeHiGain()+1.;
422 upperrange = (Float_t)fHist->GetTimeUpperFitRangeHiGain()+1.;
423 }
424
425
426 if (fTime < lowerrange)
427 {
428 *fLog << warn
429 << "WARNING: Mean Fitted Time inside or smaller than first used FADC slice in Pixel "
430 << fPixId << " time: " << fTime << " Range: " << lowerrange << endl;
431 return kFALSE;
432 }
433
434 if (fTime > upperrange)
435 {
436 *fLog << warn
437 << "WARNING: Mean Fitted Time inside or greater than last used FADC slice in Pixel "
438 << fPixId << " time: " << fTime << " Range: " << upperrange << endl;
439 return kFALSE;
440 }
441
442 return kTRUE;
443}
444
445
446//
447// The check returns kTRUE if:
448//
449//
450//
451Bool_t MCalibrationPix::CheckOscillations()
452{
453
454
455 return kTRUE;
456}
457
458
459
460void MCalibrationPix::ApplyLoGainConversion()
461{
462
463 const Float_t chargeRelErrSquare = fErrCharge*fErrCharge
464 /( fCharge * fCharge);
465 const Float_t sigmaRelErrSquare = fErrSigmaCharge*fErrSigmaCharge
466 /( fSigmaCharge * fSigmaCharge);
467 const Float_t conversionRelErrSquare = fConversionHiLoError*fConversionHiLoError
468 /(fConversionHiLo * fConversionHiLo);
469
470 fCharge *= fConversionHiLo;
471 fErrCharge = TMath::Sqrt(chargeRelErrSquare + conversionRelErrSquare) * fCharge;
472
473 fSigmaCharge *= fConversionHiLo;
474 fErrSigmaCharge = TMath::Sqrt(sigmaRelErrSquare + conversionRelErrSquare) * fSigmaCharge;
475
476}
477
478
479// --------------------------------------------------------------------------
480//
481// Set the pedestals from outside
482//
483void MCalibrationPix::SetPedestal(Float_t ped, Float_t pedrms)
484{
485
486 fPed = ped;
487 fPedRms = pedrms;
488
489}
490
491// --------------------------------------------------------------------------
492//
493// 1) Fit the arrival times
494// 2) Retrieve the results
495// 3) Note that because of the low number of bins, the NDf is sometimes 0, so
496// Root does not give a reasonable Probability, the Chisquare is more significant
497//
498// This fit has to be done AFTER the Charges fit,
499// otherwise only the Hi Gain will be fitted, even if there are no entries
500//
501//
502Bool_t MCalibrationPix::FitTime()
503{
504
505 //
506 // Fit the Low Gain
507 //
508 if (fHiGainSaturation)
509 {
510 if(!fHist->FitTimeLoGain())
511 {
512 *fLog << warn << "WARNING: Could not fit Lo Gain times of pixel " << fPixId << endl;
513 fHist->PrintTimeFitResult();
514 return kFALSE;
515 }
516 }
517
518 //
519 // Fit the High Gain
520 //
521 else
522 {
523 if(!fHist->FitTimeHiGain())
524 {
525 *fLog << warn << "WARNING: Could not fit Hi Gain times of pixel " << fPixId << endl;
526 fHist->PrintTimeFitResult();
527 return kFALSE;
528 }
529 }
530
531 fTime = fHist->GetTimeMean();
532 fSigmaTime = fHist->GetTimeSigma();
533 fTimeChiSquare = fHist->GetTimeChiSquare();
534 fTimeProb = fHist->GetTimeProb();
535
536 if (!CheckTimeFitValidity())
537 fFitValid = kFALSE;
538
539 return kTRUE;
540}
541
Note: See TracBrowser for help on using the repository browser.