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

Last change on this file since 2934 was 2934, checked in by gaug, 22 years ago
*** empty log message ***
File size: 20.7 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 fElectronicPedRms(1.5),
59 fErrElectronicPedRms(0.3),
60 fFactor(1.32),
61 fFactorError(0.04),
62 fChargeLimit(3.),
63 fChargeErrLimit(0.),
64 fChargeRelErrLimit(1.),
65 fFlags(0)
66{
67
68 fName = name ? name : "MCalibrationPixel";
69 fTitle = title ? title : "Container of the MHCalibrationPixels and the fit results";
70
71 //
72 // At the moment, we don't have a database, yet,
73 // so we get it from the configuration file
74 //
75 fConversionHiLo = gkConversionHiLo;
76 fConversionHiLoError = gkConversionHiLoError;
77
78 fHist = new MHCalibrationPixel("MHCalibrationPixel","Calibration Histograms Pixel ");
79
80 if (!fHist)
81 *fLog << warn << dbginf << " Could not create MHCalibrationPixel " << endl;
82
83 Clear();
84}
85
86MCalibrationPix::~MCalibrationPix()
87{
88 delete fHist;
89}
90
91
92// ------------------------------------------------------------------------
93//
94// Invalidate values
95//
96void MCalibrationPix::Clear(Option_t *o)
97{
98
99 fHist->Reset();
100
101 CLRBIT(fFlags, kHiGainSaturation);
102 CLRBIT(fFlags, kExcluded);
103 CLRBIT(fFlags, kFitValid);
104 CLRBIT(fFlags, kFitted);
105 CLRBIT(fFlags, kBlindPixelMethodValid);
106 CLRBIT(fFlags, kFFactorMethodValid);
107 CLRBIT(fFlags, kPINDiodeMethodValid);
108
109 fCharge = -1.;
110 fErrCharge = -1.;
111 fSigmaCharge = -1.;
112 fErrSigmaCharge = -1.;
113 fRSigmaSquare = -1.;
114 fChargeProb = -1.;
115 fPed = -1.;
116 fPedRms = -1.;
117 fErrPedRms = 0.;
118 fTime = -1.;
119 fSigmaTime = -1.;
120 fTimeChiSquare = -1.;
121 fPheFFactorMethod = -1.;
122 fPheFFactorMethodError = -1.;
123 fConversionFFactorMethod = -1.;
124 fConversionBlindPixelMethod = -1.;
125 fConversionPINDiodeMethod = -1.;
126 fConversionErrorFFactorMethod = -1.;
127 fConversionErrorBlindPixelMethod = -1.;
128 fConversionErrorPINDiodeMethod = -1.;
129 fConversionSigmaFFactorMethod = -1.;
130 fConversionSigmaBlindPixelMethod = -1.;
131 fConversionSigmaPINDiodeMethod = -1.;
132
133}
134
135
136void MCalibrationPix::DefinePixId(Int_t i)
137{
138
139 fPixId = i;
140 fHist->ChangeHistId(i);
141
142}
143
144
145// --------------------------------------------------------------------------
146//
147// Set the pedestals from outside
148//
149void MCalibrationPix::SetPedestal(Float_t ped, Float_t pedrms)
150{
151
152 fPed = ped;
153 fPedRms = pedrms;
154
155}
156
157// --------------------------------------------------------------------------
158//
159// Set the conversion factors from outside (only for MC)
160//
161void MCalibrationPix::SetConversionFFactorMethod(Float_t c, Float_t err, Float_t sig)
162{
163 fConversionFFactorMethod = c;
164 fConversionErrorFFactorMethod = err;
165 fConversionSigmaFFactorMethod = sig;
166}
167
168
169// --------------------------------------------------------------------------
170//
171// Set the conversion factors from outside (only for MC)
172//
173void MCalibrationPix::SetConversionBlindPixelMethod(Float_t c, Float_t err, Float_t sig)
174{
175 fConversionBlindPixelMethod = c;
176 fConversionErrorBlindPixelMethod = err;
177 fConversionSigmaBlindPixelMethod = sig;
178}
179
180// --------------------------------------------------------------------------
181//
182// Set the conversion factors from outside (only for MC)
183//
184void MCalibrationPix::SetConversionPINDiodeMethod(Float_t c, Float_t err, Float_t sig)
185{
186 fConversionPINDiodeMethod = c ;
187 fConversionErrorPINDiodeMethod = err;
188 fConversionSigmaPINDiodeMethod = sig;
189}
190
191// --------------------------------------------------------------------------
192//
193// Set the Hi Gain Saturation Bit from outside (only for MC)
194//
195void MCalibrationPix::SetHiGainSaturation(Bool_t b)
196{
197
198 if (b)
199 {
200 SETBIT(fFlags, kHiGainSaturation);
201 fHist->SetUseLoGain(1);
202 }
203 else
204 {
205 CLRBIT(fFlags, kHiGainSaturation);
206 fHist->SetUseLoGain(0);
207 }
208}
209
210// --------------------------------------------------------------------------
211//
212// Set the Excluded Bit from outside
213//
214void MCalibrationPix::SetExcluded(Bool_t b )
215{
216 b ? SETBIT(fFlags, kExcluded) : CLRBIT(fFlags, kExcluded);
217}
218
219
220// --------------------------------------------------------------------------
221//
222// Set the Excluded Bit from outside
223//
224void MCalibrationPix::SetExcludeQualityCheck(Bool_t b )
225{
226 b ? SETBIT(fFlags, kExcludeQualityCheck) : CLRBIT(fFlags, kExcludeQualityCheck);
227}
228
229// --------------------------------------------------------------------------
230//
231// Set the Excluded Bit from outside
232//
233void MCalibrationPix::SetFitValid(Bool_t b )
234{
235 b ? SETBIT(fFlags, kFitValid) : CLRBIT(fFlags, kFitValid);
236}
237
238// --------------------------------------------------------------------------
239//
240// Set the Excluded Bit from outside
241//
242void MCalibrationPix::SetFitted(Bool_t b )
243{
244 b ? SETBIT(fFlags, kFitted) : CLRBIT(fFlags, kFitted);
245}
246
247// --------------------------------------------------------------------------
248//
249// Set the Excluded Bit from outside
250//
251void MCalibrationPix::SetBlindPixelMethodValid(Bool_t b )
252{
253 b ? SETBIT(fFlags, kBlindPixelMethodValid) : CLRBIT(fFlags, kBlindPixelMethodValid);
254}
255
256// --------------------------------------------------------------------------
257//
258// Set the Excluded Bit from outside
259//
260void MCalibrationPix::SetFFactorMethodValid(Bool_t b )
261{
262 b ? SETBIT(fFlags, kFFactorMethodValid) : CLRBIT(fFlags, kFFactorMethodValid);
263}
264
265// --------------------------------------------------------------------------
266//
267// Set the Excluded Bit from outside
268//
269void MCalibrationPix::SetPINDiodeMethodValid(Bool_t b )
270{
271 b ? SETBIT(fFlags, kPINDiodeMethodValid) : CLRBIT(fFlags, kPINDiodeMethodValid);
272}
273
274
275Bool_t MCalibrationPix::IsExcluded() const
276 {
277 return TESTBIT(fFlags,kExcluded);
278 }
279
280Bool_t MCalibrationPix::IsFitValid() const
281{
282 return TESTBIT(fFlags, kFitValid);
283}
284
285Bool_t MCalibrationPix::IsFitted() const
286{
287 return TESTBIT(fFlags, kFitted);
288}
289
290Bool_t MCalibrationPix::IsBlindPixelMethodValid() const
291{
292 return TESTBIT(fFlags, kBlindPixelMethodValid);
293}
294
295Bool_t MCalibrationPix::IsFFactorMethodValid() const
296{
297 return TESTBIT(fFlags, kFFactorMethodValid);
298}
299
300Bool_t MCalibrationPix::IsPINDiodeMethodValid() const
301{
302 return TESTBIT(fFlags, kPINDiodeMethodValid);
303}
304
305
306// --------------------------------------------------------------------------
307//
308// 1) Return if the charge distribution is already succesfully fitted
309// or if the histogram is empty
310// 2) Set a lower Fit range according to 1.5 Pedestal RMS in order to avoid
311// possible remaining cosmics to spoil the fit.
312// 3) Decide if the LoGain Histogram is fitted or the HiGain Histogram
313// 4) Fit the histograms with a Gaussian
314// 5) In case of failure print out the fit results
315// 6) Retrieve the results and store them in this class
316// 7) Calculate the number of photo-electrons after the F-Factor method
317// 8) Calculate the errors of the F-Factor method
318//
319// The fits are declared valid (fFitValid = kTRUE), if:
320//
321// 1) Pixel has a fitted charge greater than 5*PedRMS
322// 2) Pixel has a fit error greater than 0.
323// 3) Pixel has a fit Probability greater than 0.0001
324// 4) Pixel has a charge sigma bigger than its Pedestal RMS
325// 5) If FitTimes is used,
326// the mean arrival time is at least 1.0 slices from the used edge slices
327// (this stage is only performed in the times fit)
328//
329// If the histogram is empty, all values are set to -1.
330//
331// The conversion factor after the F-Factor method is declared valid, if:
332//
333// 1) fFitValid is kTRUE
334// 2) Conversion Factor is bigger than 0.
335// 3) The error of the conversion factor is smaller than 10%
336//
337Bool_t MCalibrationPix::FitCharge()
338{
339
340 //
341 // 1) Return if the charge distribution is already succesfully fitted
342 // or if the histogram is empty
343 //
344 if (fHist->IsFitOK() || fHist->IsEmpty())
345 return kTRUE;
346
347 //
348 // 2) Set a lower Fit range according to 1.5 Pedestal RMS in order to avoid
349 // possible remaining cosmics to spoil the fit.
350 //
351 // if (fPed && fPedRms)
352 // fHist->SetLowerFitRange(1.5*fPedRms);
353 // else
354 // *fLog << warn << "WARNING: Cannot set lower fit range: Pedestals not available" << endl;
355
356 //
357 // 3) Decide if the LoGain Histogram is fitted or the HiGain Histogram
358 //
359 if (fHist->UseLoGain())
360 SetHiGainSaturation();
361
362 //
363 // 4) Fit the Lo Gain histograms with a Gaussian
364 //
365 if(fHist->FitCharge())
366 {
367 SETBIT(fFlags,kFitted);
368 }
369 else
370 {
371 *fLog << warn << "WARNING: Could not fit charges of pixel " << fPixId << endl;
372 //
373 // 5) In case of failure print out the fit results
374 //
375 // fHist->PrintChargeFitResult();
376 CLRBIT(fFlags,kFitted);
377 }
378
379 //
380 // 6) Retrieve the results and store them in this class
381 //
382 fCharge = fHist->GetChargeMean();
383 fErrCharge = fHist->GetChargeMeanErr();
384 fSigmaCharge = fHist->GetChargeSigma();
385 fErrSigmaCharge = fHist->GetChargeSigmaErr();
386 fChargeProb = fHist->GetChargeProb();
387
388 if (CheckChargeFitValidity())
389 SETBIT(fFlags,kFitValid);
390 else
391 {
392 CLRBIT(fFlags,kFitValid);
393 return kFALSE;
394 }
395
396 //
397 // 7) Calculate the number of photo-electrons after the F-Factor method
398 // 8) Calculate the errors of the F-Factor method
399 //
400 if ((fPed > 0.) && (fPedRms > 0.))
401 {
402
403 //
404 // Square all variables in order to avoid applications of square root
405 //
406 // First the relative error squares
407 //
408 const Float_t chargeSquare = fCharge* fCharge;
409 const Float_t chargeSquareRelErrSquare = 4.*fErrCharge*fErrCharge / chargeSquare;
410
411 const Float_t fFactorRelErrSquare = fFactorError * fFactorError / (fFactor * fFactor);
412 //
413 // Now the absolute error squares
414 //
415 const Float_t sigmaSquare = fSigmaCharge* fSigmaCharge;
416 const Float_t sigmaSquareErrSquare = 4.*fErrSigmaCharge*fErrSigmaCharge * sigmaSquare;
417
418 const Float_t elecRmsSquare = fElectronicPedRms* fElectronicPedRms;
419 const Float_t elecRmsSquareErrSquare = 4.*fErrElectronicPedRms*fErrElectronicPedRms * elecRmsSquare;
420
421 Float_t pedRmsSquare = fPedRms* fPedRms;
422 Float_t pedRmsSquareErrSquare = 4.*fErrPedRms*fErrPedRms * pedRmsSquare;
423
424 if (TESTBIT(fFlags,kHiGainSaturation))
425 {
426
427 //
428 // We do not know the Lo Gain Pedestal RMS, so we have to retrieve it
429 // from the Hi Gain:
430 //
431 // We extract the pure NSB contribution:
432 //
433 Float_t nsbSquare = pedRmsSquare - elecRmsSquare;
434 Float_t nsbSquareRelErrSquare = (pedRmsSquareErrSquare + elecRmsSquareErrSquare)
435 / (nsbSquare * nsbSquare) ;
436
437 if (nsbSquare < 0.)
438 nsbSquare = 0.;
439
440 //
441 // Now, we divide the NSB by the conversion factor and
442 // add it quadratically to the electronic noise
443 //
444 const Float_t conversionSquare = fConversionHiLo *fConversionHiLo;
445 const Float_t conversionSquareRelErrSquare = 4.*fConversionHiLoError*fConversionHiLoError/conversionSquare;
446
447 //
448 // Calculate the new "Pedestal RMS"
449 //
450 const Float_t convertedNsbSquare = nsbSquare / conversionSquare;
451 const Float_t convertedNsbSquareErrSquare = (nsbSquareRelErrSquare + conversionSquareRelErrSquare)
452 * convertedNsbSquare * convertedNsbSquare;
453
454 pedRmsSquare = convertedNsbSquare + elecRmsSquare;
455 pedRmsSquareErrSquare = convertedNsbSquareErrSquare + elecRmsSquareErrSquare;
456
457 } /* if (kHiGainSaturation) */
458
459 //
460 // Calculate the reduced sigmas
461 //
462 fRSigmaSquare = sigmaSquare - pedRmsSquare;
463 if (fRSigmaSquare <= 0.)
464 {
465 *fLog << warn
466 << "WARNING: Cannot apply F-Factor calibration: Reduced Sigma smaller than 0 in pixel "
467 << fPixId << endl;
468 if (TESTBIT(fFlags,kHiGainSaturation))
469 ApplyLoGainConversion();
470 return kFALSE;
471 }
472
473 const Float_t rSigmaSquareRelErrSquare = (sigmaSquareErrSquare + pedRmsSquareErrSquare)
474 / (fRSigmaSquare * fRSigmaSquare) ;
475
476 //
477 // Calculate the number of phe's from the F-Factor method
478 // (independent on Hi Gain or Lo Gain)
479 //
480 fPheFFactorMethod = fFactor * chargeSquare / fRSigmaSquare;
481
482 const Float_t pheFFactorRelErrSquare = fFactorRelErrSquare
483 + chargeSquareRelErrSquare
484 + rSigmaSquareRelErrSquare ;
485
486 fPheFFactorMethodError = TMath::Sqrt(pheFFactorRelErrSquare) * fPheFFactorMethod;
487
488 //
489 // Calculate the conversion factors
490 //
491 if (TESTBIT(fFlags,kHiGainSaturation))
492 ApplyLoGainConversion();
493
494 const Float_t chargeRelErrSquare = fErrCharge*fErrCharge / (fCharge * fCharge);
495
496 fConversionFFactorMethod = fPheFFactorMethod / fCharge ;
497 fConversionErrorFFactorMethod = ( pheFFactorRelErrSquare + chargeRelErrSquare )
498 * fConversionFFactorMethod * fConversionFFactorMethod;
499
500 if ( IsFitValid() &&
501 (fConversionFFactorMethod > 0.) &&
502 (fConversionErrorFFactorMethod/fConversionFFactorMethod < 0.1) )
503 SETBIT(fFlags,kFFactorMethodValid);
504 else
505 CLRBIT(fFlags,kFFactorMethodValid);
506
507 } /* if ((fPed > 0.) && (fPedRms > 0.)) */
508
509 return kTRUE;
510
511}
512
513//
514// The check return kTRUE if:
515//
516// 1) Pixel has a fitted charge greater than 5*PedRMS
517// 2) Pixel has a fit error greater than 0.
518// 3) Pixel has a fitted charge greater its charge error
519// 4) Pixel has a fit Probability greater than 0.0001
520// 5) Pixel has a charge sigma bigger than its Pedestal RMS
521//
522Bool_t MCalibrationPix::CheckChargeFitValidity()
523{
524
525 if (TESTBIT(fFlags,kExcludeQualityCheck))
526 return kTRUE;
527
528 Float_t equivpedestal = GetPedRms();
529
530 if (TESTBIT(fFlags,kHiGainSaturation))
531 equivpedestal /= fConversionHiLo;
532
533 if (fCharge < fChargeLimit*equivpedestal)
534 {
535 *fLog << warn << "WARNING: Fitted Charge is smaller than "
536 << fChargeLimit << " Pedestal RMS in Pixel " << fPixId << endl;
537 return kFALSE;
538 }
539
540 if (fErrCharge < fChargeErrLimit)
541 {
542 *fLog << warn << "WARNING: Error of Fitted Charge is smaller than "
543 << fChargeErrLimit << " in Pixel " << fPixId << endl;
544 return kFALSE;
545 }
546
547 if (fCharge < fChargeRelErrLimit*fErrCharge)
548 {
549 *fLog << warn << "WARNING: Fitted Charge is smaller than "
550 << fChargeRelErrLimit << "* its error in Pixel " << fPixId << endl;
551 return kFALSE;
552 }
553
554 if (!fHist->IsFitOK())
555 {
556 *fLog << warn << "WARNING: Probability of Fitted Charge too low in Pixel " << fPixId << endl;
557 return kFALSE;
558 }
559
560 if (fSigmaCharge < equivpedestal)
561 {
562 *fLog << warn << "WARNING: Sigma of Fitted Charge smaller than Pedestal RMS in Pixel " << fPixId << endl;
563 return kFALSE;
564 }
565 return kTRUE;
566}
567
568//
569// The check returns kTRUE if:
570//
571// The mean arrival time is at least 1.0 slices from the used edge slices
572//
573Bool_t MCalibrationPix::CheckTimeFitValidity()
574{
575
576 if (TESTBIT(fFlags,kExcludeQualityCheck))
577 return kTRUE;
578
579 Float_t lowerrange;
580 Float_t upperrange;
581
582 if (TESTBIT(fFlags,kHiGainSaturation))
583 {
584 lowerrange = (Float_t)fHist->GetTimeLowerFitRangeLoGain()+1.;
585 upperrange = (Float_t)fHist->GetTimeUpperFitRangeLoGain()+1.;
586 }
587 else
588 {
589 lowerrange = (Float_t)fHist->GetTimeLowerFitRangeHiGain()+1.;
590 upperrange = (Float_t)fHist->GetTimeUpperFitRangeHiGain()+1.;
591 }
592
593
594 if (fTime < lowerrange)
595 {
596 *fLog << warn
597 << "WARNING: Mean Fitted Time inside or smaller than first used FADC slice in Pixel "
598 << fPixId << " time: " << fTime << " Range: " << lowerrange << endl;
599 return kFALSE;
600 }
601
602 if (fTime > upperrange)
603 {
604 *fLog << warn
605 << "WARNING: Mean Fitted Time inside or greater than last used FADC slice in Pixel "
606 << fPixId << " time: " << fTime << " Range: " << upperrange << endl;
607 return kFALSE;
608 }
609
610 return kTRUE;
611}
612
613
614//
615// The check returns kTRUE if:
616//
617//
618//
619Bool_t MCalibrationPix::CheckOscillations()
620{
621
622
623 return kTRUE;
624}
625
626
627
628void MCalibrationPix::ApplyLoGainConversion()
629{
630
631 const Float_t chargeRelErrSquare = fErrCharge*fErrCharge
632 /( fCharge * fCharge);
633 const Float_t sigmaRelErrSquare = fErrSigmaCharge*fErrSigmaCharge
634 /( fSigmaCharge * fSigmaCharge);
635 const Float_t conversionRelErrSquare = fConversionHiLoError*fConversionHiLoError
636 /(fConversionHiLo * fConversionHiLo);
637
638 fCharge *= fConversionHiLo;
639 fErrCharge = TMath::Sqrt(chargeRelErrSquare + conversionRelErrSquare) * fCharge;
640
641 fSigmaCharge *= fConversionHiLo;
642 fErrSigmaCharge = TMath::Sqrt(sigmaRelErrSquare + conversionRelErrSquare) * fSigmaCharge;
643
644}
645
646
647
648// --------------------------------------------------------------------------
649//
650// 1) Fit the arrival times
651// 2) Retrieve the results
652// 3) Note that because of the low number of bins, the NDf is sometimes 0, so
653// Root does not give a reasonable Probability, the Chisquare is more significant
654//
655// This fit has to be done AFTER the Charges fit,
656// otherwise only the Hi Gain will be fitted, even if there are no entries
657//
658//
659Bool_t MCalibrationPix::FitTime()
660{
661
662 //
663 // Fit the Low Gain
664 //
665 if (TESTBIT(fFlags,kHiGainSaturation))
666 {
667 if(!fHist->FitTimeLoGain())
668 {
669 *fLog << warn << "WARNING: Could not fit Lo Gain times of pixel " << fPixId << endl;
670 // fHist->PrintTimeFitResult();
671 return kFALSE;
672 }
673 }
674
675 //
676 // Fit the High Gain
677 //
678 else
679 {
680 if(!fHist->FitTimeHiGain())
681 {
682 *fLog << warn << "WARNING: Could not fit Hi Gain times of pixel " << fPixId << endl;
683 // fHist->PrintTimeFitResult();
684 return kFALSE;
685 }
686 }
687
688 fTime = fHist->GetTimeMean();
689 fSigmaTime = fHist->GetTimeSigma();
690 fTimeChiSquare = fHist->GetTimeChiSquare();
691 fTimeProb = fHist->GetTimeProb();
692
693 if (CheckTimeFitValidity())
694 SETBIT(fFlags,kFitValid);
695 else
696 CLRBIT(fFlags,kFitValid);
697
698 return kTRUE;
699}
700
Note: See TracBrowser for help on using the repository browser.