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

Last change on this file since 3243 was 3242, checked in by gaug, 21 years ago
*** empty log message ***
File size: 28.0 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// 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#include "MCalibrationConfig.h"
89
90#include "MLog.h"
91#include "MLogManip.h"
92
93ClassImp(MCalibrationPix);
94
95using namespace std;
96
97const Float_t MCalibrationPix::gkElectronicPedRms = 1.5;
98const Float_t MCalibrationPix::gkElectronicPedRmsErr = 0.3;
99const Float_t MCalibrationPix::gkFFactor = 1.15;
100const Float_t MCalibrationPix::gkFFactorErr = 0.02;
101const Float_t MCalibrationPix::gkChargeLimit = 3.;
102const Float_t MCalibrationPix::gkChargeErrLimit = 0.;
103const Float_t MCalibrationPix::gkChargeRelErrLimit = 1.;
104const Float_t MCalibrationPix::gkTimeLimit = 1.5;
105const Float_t MCalibrationPix::gkTimeErrLimit = 3.;
106const Float_t MCalibrationPix::gkConvFFactorRelErrLimit = 0.1;
107
108const Float_t MCalibrationPix::gkAverageQE = 0.18;
109const Float_t MCalibrationPix::gkAverageQEErr = 0.02;
110const Float_t MCalibrationPix::gkConversionHiLo = 10.;
111const Float_t MCalibrationPix::gkConversionHiLoErr = 2.5;
112// --------------------------------------------------------------------------
113//
114// Default Constructor:
115//
116MCalibrationPix::MCalibrationPix(const char *name, const char *title)
117 : fPixId(-1),
118 fFlags(0)
119{
120
121 fName = name ? name : "MCalibrationPixel";
122 fTitle = title ? title : "Container of the MHCalibrationPixels and the fit results";
123
124 fHist = new MHCalibrationPixel("MHCalibrationPixel","Calibration Histograms Pixel ");
125
126 if (!fHist)
127 *fLog << warn << dbginf << " Could not create MHCalibrationPixel " << endl;
128
129 Clear();
130
131 //
132 // At the moment, we don't have a database, yet,
133 // so we get it from the configuration file
134 //
135 SetConversionHiLo();
136 SetConversionHiLoErr();
137
138 SetAverageQE();
139}
140
141MCalibrationPix::~MCalibrationPix()
142{
143 delete fHist;
144}
145
146
147// ------------------------------------------------------------------------
148//
149// Invalidate values
150//
151void MCalibrationPix::Clear(Option_t *o)
152{
153
154 fHist->Reset();
155
156 CLRBIT(fFlags, kHiGainSaturation);
157 CLRBIT(fFlags, kExcluded);
158 CLRBIT(fFlags, kExcludeQualityCheck);
159 CLRBIT(fFlags, kChargeValid);
160 CLRBIT(fFlags, kFitted);
161 CLRBIT(fFlags, kOscillating);
162 CLRBIT(fFlags, kBlindPixelMethodValid);
163 CLRBIT(fFlags, kFFactorMethodValid);
164 CLRBIT(fFlags, kPINDiodeMethodValid);
165 CLRBIT(fFlags, kCombinedMethodValid);
166
167 fCharge = -1.;
168 fChargeErr = -1.;
169 fSigmaCharge = -1.;
170 fSigmaChargeErr = -1.;
171 fRSigmaCharge = -1.;
172 fRSigmaChargeErr = -1.;
173
174 fChargeProb = -1.;
175 fPed = -1.;
176 fPedRms = -1.;
177 fPedRmsErr = -1.;
178
179 fNumHiGainSamples = -1.;
180 fNumLoGainSamples = -1.;
181
182 fTimeFirstHiGain = 0 ;
183 fTimeLastHiGain = 0 ;
184 fTimeFirstLoGain = 0 ;
185 fTimeLastLoGain = 0 ;
186
187 fAbsTimeMean = -1.;
188 fAbsTimeRms = -1.;
189
190 fPheFFactorMethod = -1.;
191 fPheFFactorMethodErr = -1.;
192
193 fMeanConversionFFactorMethod = -1.;
194 fMeanConversionBlindPixelMethod = -1.;
195 fMeanConversionPINDiodeMethod = -1.;
196 fMeanConversionCombinedMethod = -1.;
197
198 fConversionFFactorMethodErr = -1.;
199 fConversionBlindPixelMethodErr = -1.;
200 fConversionPINDiodeMethodErr = -1.;
201 fConversionCombinedMethodErr = -1.;
202
203 fSigmaConversionFFactorMethod = -1.;
204 fSigmaConversionBlindPixelMethod = -1.;
205 fSigmaConversionPINDiodeMethod = -1.;
206 fSigmaConversionCombinedMethod = -1.;
207
208 fTotalFFactorFFactorMethod = -1.;
209 fTotalFFactorBlindPixelMethod = -1.;
210 fTotalFFactorPINDiodeMethod = -1.;
211 fTotalFFactorCombinedMethod = -1.;
212
213}
214
215
216void MCalibrationPix::DefinePixId(Int_t i)
217{
218
219 fPixId = i;
220 fHist->ChangeHistId(i);
221
222}
223
224
225// --------------------------------------------------------------------------
226//
227// Set the pedestals from outside
228//
229void MCalibrationPix::SetPedestal(const Float_t ped, const Float_t pedrms,
230 const Float_t higainsamp, const Float_t logainsamp )
231{
232
233 fPed = ped;
234 fPedRms = pedrms;
235
236 fNumHiGainSamples = higainsamp;
237 fNumLoGainSamples = logainsamp;
238
239}
240
241// --------------------------------------------------------------------------
242//
243// Set the conversion factors from outside (only for MC)
244//
245void MCalibrationPix::SetConversionFFactorMethod(Float_t c, Float_t err, Float_t sig)
246{
247 fMeanConversionFFactorMethod = c;
248 fConversionFFactorMethodErr = err;
249 fSigmaConversionFFactorMethod = sig;
250}
251
252// --------------------------------------------------------------------------
253//
254// Set the conversion factors from outside (only for MC)
255//
256void MCalibrationPix::SetConversionCombinedMethod(Float_t c, Float_t err, Float_t sig)
257{
258 fMeanConversionCombinedMethod = c;
259 fConversionCombinedMethodErr = err;
260 fSigmaConversionCombinedMethod = sig;
261}
262
263
264// --------------------------------------------------------------------------
265//
266// Set the conversion factors from outside (only for MC)
267//
268void MCalibrationPix::SetConversionBlindPixelMethod(Float_t c, Float_t err, Float_t sig)
269{
270 fMeanConversionBlindPixelMethod = c;
271 fConversionBlindPixelMethodErr = err;
272 fSigmaConversionBlindPixelMethod = sig;
273}
274
275// --------------------------------------------------------------------------
276//
277// Set the conversion factors from outside (only for MC)
278//
279void MCalibrationPix::SetConversionPINDiodeMethod(Float_t c, Float_t err, Float_t sig)
280{
281 fMeanConversionPINDiodeMethod = c ;
282 fConversionPINDiodeMethodErr = err;
283 fSigmaConversionPINDiodeMethod = sig;
284}
285
286// --------------------------------------------------------------------------
287//
288// Set the Hi Gain Saturation Bit from outside (only for MC)
289//
290void MCalibrationPix::SetHiGainSaturation(Bool_t b)
291{
292
293 if (b)
294 {
295 SETBIT(fFlags, kHiGainSaturation);
296 fHist->SetUseLoGain(1);
297 }
298 else
299 {
300 CLRBIT(fFlags, kHiGainSaturation);
301 fHist->SetUseLoGain(0);
302 }
303}
304
305// --------------------------------------------------------------------------
306//
307// Set the Excluded Bit from outside
308//
309void MCalibrationPix::SetExcluded(Bool_t b )
310{
311 b ? SETBIT(fFlags, kExcluded) : CLRBIT(fFlags, kExcluded);
312}
313
314
315// --------------------------------------------------------------------------
316//
317// Set the Excluded Bit from outside
318//
319void MCalibrationPix::SetExcludeQualityCheck(Bool_t b )
320{
321 b ? SETBIT(fFlags, kExcludeQualityCheck) : CLRBIT(fFlags, kExcludeQualityCheck);
322}
323
324// --------------------------------------------------------------------------
325//
326// Set the Excluded Bit from outside
327//
328void MCalibrationPix::SetChargeValid(Bool_t b )
329{
330 b ? SETBIT(fFlags, kChargeValid) : CLRBIT(fFlags, kChargeValid);
331}
332
333// --------------------------------------------------------------------------
334//
335// Set the Excluded Bit from outside
336//
337void MCalibrationPix::SetFitted(Bool_t b )
338{
339 b ? SETBIT(fFlags, kFitted) : CLRBIT(fFlags, kFitted);
340}
341
342// --------------------------------------------------------------------------
343//
344// Set the Excluded Bit from outside
345//
346void MCalibrationPix::SetOscillating(Bool_t b )
347{
348 b ? SETBIT(fFlags, kOscillating) : CLRBIT(fFlags, kOscillating);
349}
350
351// --------------------------------------------------------------------------
352//
353// Set the Excluded Bit from outside
354//
355void MCalibrationPix::SetBlindPixelMethodValid(Bool_t b )
356{
357 b ? SETBIT(fFlags, kBlindPixelMethodValid) : CLRBIT(fFlags, kBlindPixelMethodValid);
358}
359
360// --------------------------------------------------------------------------
361//
362// Set the Excluded Bit from outside
363//
364void MCalibrationPix::SetFFactorMethodValid(Bool_t b )
365{
366 b ? SETBIT(fFlags, kFFactorMethodValid) : CLRBIT(fFlags, kFFactorMethodValid);
367}
368
369// --------------------------------------------------------------------------
370//
371// Set the Excluded Bit from outside
372//
373void MCalibrationPix::SetPINDiodeMethodValid(Bool_t b )
374{
375 b ? SETBIT(fFlags, kPINDiodeMethodValid) : CLRBIT(fFlags, kPINDiodeMethodValid);
376}
377
378void MCalibrationPix::SetAbsTimeBordersHiGain(Byte_t f, Byte_t l)
379{
380
381 fTimeFirstHiGain = f;
382 fTimeLastHiGain = l;
383
384}
385
386void MCalibrationPix::SetAbsTimeBordersLoGain(Byte_t f, Byte_t l)
387{
388
389 fTimeFirstLoGain = f;
390 fTimeLastLoGain = l;
391
392}
393
394Bool_t MCalibrationPix::IsExcluded() const
395{
396 return TESTBIT(fFlags,kExcluded);
397}
398
399Bool_t MCalibrationPix::IsExcludeQualityCheck() const
400{
401 return TESTBIT(fFlags,kExcludeQualityCheck);
402}
403
404Bool_t MCalibrationPix::IsHiGainSaturation() const
405{
406 return TESTBIT(fFlags,kHiGainSaturation);
407}
408
409Bool_t MCalibrationPix::IsChargeValid() const
410{
411 return TESTBIT(fFlags, kChargeValid);
412}
413
414Bool_t MCalibrationPix::IsFitted() const
415{
416 return TESTBIT(fFlags, kFitted);
417}
418
419Bool_t MCalibrationPix::IsOscillating() const
420{
421 return TESTBIT(fFlags, kOscillating);
422}
423
424Bool_t MCalibrationPix::IsBlindPixelMethodValid() const
425{
426 return TESTBIT(fFlags, kBlindPixelMethodValid);
427}
428
429Bool_t MCalibrationPix::IsFFactorMethodValid() const
430{
431 return TESTBIT(fFlags, kFFactorMethodValid);
432}
433
434Bool_t MCalibrationPix::IsPINDiodeMethodValid() const
435{
436 return TESTBIT(fFlags, kPINDiodeMethodValid);
437}
438
439Bool_t MCalibrationPix::IsCombinedMethodValid() const
440{
441 return TESTBIT(fFlags, kCombinedMethodValid);
442}
443
444
445// --------------------------------------------------------------------------
446//
447// 1) Return if the charge distribution is already succesfully fitted
448// or if the histogram is empty
449// 2) Set a lower Fit range according to 1.5 Pedestal RMS in order to avoid
450// possible remaining cosmics to spoil the fit.
451// 3) Decide if the LoGain Histogram is fitted or the HiGain Histogram
452// 4) Fit the histograms with a Gaussian
453// 5) In case of failure set the bit kFitted to false
454// 6) Retrieve the results and store them in this class
455// 7) Calculate the number of photo-electrons after the F-Factor method
456// 8) Calculate the errors of the F-Factor method
457//
458// The fits are declared valid (fFitValid = kTRUE), if:
459//
460// 1) Pixel has a fitted charge greater than 3*PedRMS
461// 2) Pixel has a fit error greater than 0.
462// 3) Pixel has a fit Probability greater than 0.0001
463// 4) Pixel has a charge sigma bigger than its Pedestal RMS
464// 5) If FitTimes is used,
465// the mean arrival time is at least 1.0 slices from the used edge slices
466// (this stage is only performed in the times fit)
467//
468// If the histogram is empty, all values are set to -1.
469//
470// The conversion factor after the F-Factor method is declared valid, if:
471//
472// 1) fFitValid is kTRUE
473// 2) Conversion Factor is bigger than 0.
474// 3) The error of the conversion factor is smaller than 10%
475//
476Bool_t MCalibrationPix::FitCharge()
477{
478
479 //
480 // 1) Return if the charge distribution is already succesfully fitted
481 // or if the histogram is empty
482 //
483 if (fHist->IsChargeFitOK() || fHist->IsEmpty())
484 return kTRUE;
485
486 //
487 // 2) Set a lower Fit range according to 1.5 Pedestal RMS in order to avoid
488 // possible remaining cosmics to spoil the fit.
489 //
490 // if (fPed && fPedRms)
491 // fHist->SetLowerFitRange(1.5*fPedRms);
492 // else
493 // *fLog << warn << "WARNING: Cannot set lower fit range: Pedestals not available" << endl;
494
495 //
496 // 3) Decide if the LoGain Histogram is fitted or the HiGain Histogram
497 //
498 if (fHist->UseLoGain())
499 SetHiGainSaturation();
500
501 //
502 // 4) Fit the Lo Gain histograms with a Gaussian
503 //
504 if (fHist->FitCharge())
505 SETBIT(fFlags,kFitted);
506 else
507 {
508 *fLog << warn << "WARNING: Could not fit charges of pixel " << fPixId << endl;
509 //
510 // 5) In case of failure set the bit kFitted to false
511 //
512 CLRBIT(fFlags,kFitted);
513 }
514
515 //
516 // 6) Retrieve the results and store them in this class
517 // If fFitted is false, we get the eans and RMS of the histogram!!
518 //
519 fCharge = fHist->GetChargeMean();
520 fChargeErr = fHist->GetChargeMeanErr();
521 fSigmaCharge = fHist->GetChargeSigma();
522 fSigmaChargeErr = fHist->GetChargeSigmaErr();
523 fChargeProb = fHist->GetChargeProb();
524
525
526 fAbsTimeMean = fHist->GetAbsTimeMean();
527 fAbsTimeMeanErr = fHist->GetAbsTimeMeanErr();
528 fAbsTimeRms = fHist->GetAbsTimeRms();
529
530 if (CheckTimeFitValidity())
531 SETBIT(fFlags,kTimeFitValid);
532 else
533 CLRBIT(fFlags,kTimeFitValid);
534
535 //
536 //Calculate the conversion factors
537 //
538 if (IsHiGainSaturation())
539 ApplyLoGainConversion();
540
541 if (CheckChargeValidity())
542 SETBIT(fFlags,kChargeValid);
543 else
544 {
545 CLRBIT(fFlags,kChargeValid);
546 return kFALSE;
547 }
548
549 return kTRUE;
550
551}
552
553//
554// Calculate the number of photo-electrons after the F-Factor method
555// Calculate the errors of the F-Factor method
556//
557Bool_t MCalibrationPix::CalcFFactorMethod()
558{
559
560 if ( (fCharge == -1.)
561 || (fChargeErr < 0.)
562 || (fSigmaCharge < 0.)
563 || (fPedRms < 0.) )
564 {
565 *fLog << warn << GetDescriptor() << "Cannot calculate the FFactor Method! "
566 << "Some of the needed parameters are not available ";
567 CLRBIT(fFlags,kFFactorMethodValid);
568 return kFALSE;
569 }
570
571 //
572 // Square all variables in order to avoid applications of square root
573 //
574 // First the relative error squares
575 //
576 const Float_t chargeSquare = fCharge * fCharge;
577 const Float_t chargeSquareRelErrSquare = 4.* fChargeErr * fChargeErr / chargeSquare;
578
579 const Float_t chargeRelErrSquare = fChargeErr * fChargeErr
580 / (fCharge * fCharge);
581
582 const Float_t ffactorsquare = gkFFactor * gkFFactor;
583 const Float_t ffactorsquareRelErrSquare = 4.*gkFFactorErr * gkFFactorErr / ffactorsquare;
584
585 const Float_t avQERelErrSquare = fAverageQEErr * fAverageQEErr / fAverageQE / fAverageQE;
586
587 const Float_t avQEFFactor = TMath::Sqrt( ( 1. - fAverageQE ) / fAverageQE );
588 const Float_t avQEFFactorErr = 1./ ( 2. * avQEFFactor ) * fAverageQEErr
589 / ( fAverageQE * fAverageQE );
590 const Float_t avQEFFactorRelErrSquare = avQEFFactorErr * avQEFFactorErr
591 / ( avQEFFactor * avQEFFactor) ;
592 //
593 // Now the absolute error squares
594 //
595 const Float_t sigmaSquare = fSigmaCharge * fSigmaCharge;
596 const Float_t sigmaSquareErrSquare = 4.*fSigmaChargeErr* fSigmaChargeErr * sigmaSquare;
597
598 Float_t pedRmsSquare = fPedRms * fPedRms;
599 Float_t pedRmsSquareErrSquare = 4.*fPedRmsErr * fPedRmsErr * pedRmsSquare;
600
601 if (!IsHiGainSaturation())
602 { /* HiGain */
603
604 pedRmsSquare *= fNumHiGainSamples;
605 pedRmsSquareErrSquare *= fNumHiGainSamples*fNumHiGainSamples;
606 }
607 else
608 { /* LoGain */
609
610 //
611 // We do not know the Lo Gain Pedestal RMS, so we have to retrieve it
612 // from the HI GAIN (all calculation per slice up to now):
613 //
614 // We extract the pure NSB contribution:
615 //
616 const Float_t elecRmsSquare = gkElectronicPedRms * gkElectronicPedRms;
617 const Float_t elecRmsSquareErrSquare = 4.*gkElectronicPedRmsErr * gkElectronicPedRmsErr * elecRmsSquare;
618
619 Float_t nsbSquare = pedRmsSquare - elecRmsSquare;
620 Float_t nsbSquareRelErrSquare = (pedRmsSquareErrSquare + elecRmsSquareErrSquare)
621 / (nsbSquare * nsbSquare) ;
622
623 if (nsbSquare < 0.)
624 nsbSquare = 0.;
625
626 //
627 // Now, we divide the NSB by the conversion factor and
628 // add it quadratically to the electronic noise
629 //
630 const Float_t conversionSquare = fConversionHiLo * fConversionHiLo;
631 const Float_t conversionSquareRelErrSquare = 4.*fConversionHiLoErr * fConversionHiLoErr / conversionSquare;
632
633 const Float_t convertedNsbSquare = nsbSquare / conversionSquare;
634 const Float_t convertedNsbSquareErrSquare = (nsbSquareRelErrSquare + conversionSquareRelErrSquare)
635 * convertedNsbSquare * convertedNsbSquare;
636
637 pedRmsSquare = convertedNsbSquare + elecRmsSquare;
638 pedRmsSquareErrSquare = convertedNsbSquareErrSquare + elecRmsSquareErrSquare;
639
640 //
641 // Now, correct for the number of used FADC slices in the LoGain:
642 //
643 pedRmsSquare *= fNumLoGainSamples;
644 pedRmsSquareErrSquare *= fNumLoGainSamples*fNumLoGainSamples;
645 //
646 // Correct also for the conversion to Hi-Gain:
647 //
648 pedRmsSquare *= fConversionHiLo*fConversionHiLo;
649 pedRmsSquareErrSquare *= fConversionHiLo*fConversionHiLo*fConversionHiLo*fConversionHiLo;
650
651 } /* if (HiGainSaturation) */
652
653 //
654 // Calculate the reduced sigmas
655 //
656 const Float_t rsigmachargesquare = sigmaSquare - pedRmsSquare;
657 if (rsigmachargesquare <= 0.)
658 {
659 *fLog << warn
660 << "WARNING: Cannot apply F-Factor calibration: Reduced Sigma smaller than 0 in pixel "
661 << fPixId << endl;
662 CLRBIT(fFlags,kFFactorMethodValid);
663 return kFALSE;
664 }
665
666 const Float_t rSigmaSquareRelErrSquare = (sigmaSquareErrSquare + pedRmsSquareErrSquare)
667 / (rsigmachargesquare * rsigmachargesquare) ;
668
669 fRSigmaCharge = TMath::Sqrt(rsigmachargesquare);
670 fRSigmaChargeErr = TMath::Sqrt(sigmaSquareErrSquare + pedRmsSquareErrSquare);
671
672
673 //
674 // Calculate the number of phe's from the F-Factor method
675 // (independent on Hi Gain or Lo Gain)
676 //
677 fPheFFactorMethod = ffactorsquare * chargeSquare / rsigmachargesquare;
678 //
679 // Calculate the number of photons from the F-Factor method
680 // FIXME: This is a preliminary solution, the qe shall be
681 // calibrated itself!
682 //
683 fPheFFactorMethod /= fAverageQE;
684
685 const Float_t pheFFactorRelErrSquare = ffactorsquareRelErrSquare
686 + chargeSquareRelErrSquare
687 + rSigmaSquareRelErrSquare
688 + avQERelErrSquare;
689
690 fPheFFactorMethodErr = TMath::Sqrt(pheFFactorRelErrSquare) * fPheFFactorMethod;
691
692 fMeanConversionFFactorMethod = fPheFFactorMethod / fCharge ;
693 fConversionFFactorMethodErr = ( pheFFactorRelErrSquare + chargeRelErrSquare )
694 * fMeanConversionFFactorMethod * fMeanConversionFFactorMethod;
695
696 const Float_t convrelerror = fConversionFFactorMethodErr
697 / fMeanConversionFFactorMethod;
698
699 if ( (fMeanConversionFFactorMethod > 0.) && (convrelerror < gkConvFFactorRelErrLimit))
700 SETBIT(fFlags,kFFactorMethodValid);
701
702 fSigmaConversionFFactorMethod = GetTotalFFactorFFactorMethod()*TMath::Sqrt(fMeanConversionFFactorMethod);
703
704 //
705 // Calculate the Total F-Factor of the camera ( in photons )
706 //
707 if (fPheFFactorMethod > 0)
708 {
709 fTotalFFactorFFactorMethod = (fRSigmaCharge/fCharge)*TMath::Sqrt(fPheFFactorMethod);
710 fTotalFFactorFFactorMethod *= avQEFFactor;
711 }
712
713 //
714 // Calculate the error of the Total F-Factor of the camera ( in photons )
715 //
716 const Float_t rSigmaChargeRelErrSquare = fRSigmaChargeErr * fRSigmaChargeErr
717 / (fRSigmaCharge * fRSigmaCharge) ;
718
719 fTotalFFactorErrFFactorMethod = TMath::Sqrt( rSigmaChargeRelErrSquare
720 + chargeRelErrSquare
721 + pheFFactorRelErrSquare
722 + avQEFFactorRelErrSquare );
723
724 fTotalFFactorErrFFactorMethod *= fTotalFFactorFFactorMethod;
725
726 return kTRUE;
727}
728
729
730//
731// The check returns kTRUE if:
732//
733// 0) Pixel has BIT fitted set:
734// This means:
735// a) No result is a nan
736// b) The NDF is not smaller than fNDFLimit (5)
737// c) The Probability is greater than gkProbLimit (default 0.001 == 99.9%)
738// 1) Pixel has a fitted charge greater than 3*PedRMS
739// 2) Pixel has a fit error greater than 0.
740// 3) Pixel has a fitted charge greater its charge error
741// 4) Pixel has a fit Probability greater than 0.0001
742// 5) Pixel has a charge sigma bigger than its Pedestal RMS
743//
744Bool_t MCalibrationPix::CheckChargeValidity()
745{
746
747 if (!IsFitted())
748 return kFALSE;
749
750 if (IsExcludeQualityCheck())
751 return kTRUE;
752
753 Float_t pedestal;
754
755 if (!IsHiGainSaturation()) /* higain */
756 pedestal = GetPedRms()*TMath::Sqrt(fNumHiGainSamples);
757 else /* logain */
758 pedestal = GetPedRms()*TMath::Sqrt(fNumLoGainSamples);
759
760
761 if (fCharge < gkChargeLimit*pedestal)
762 {
763 *fLog << warn << "WARNING: Fitted Charge is smaller than "
764 << gkChargeLimit << " Pedestal RMS in Pixel " << fPixId << endl;
765 return kFALSE;
766 }
767
768 if (fChargeErr < gkChargeErrLimit)
769 {
770 *fLog << warn << "WARNING: Err of Fitted Charge is smaller than "
771 << gkChargeErrLimit << " in Pixel " << fPixId << endl;
772 return kFALSE;
773 }
774
775 if (fCharge < gkChargeRelErrLimit*fChargeErr)
776 {
777 *fLog << warn << "WARNING: Fitted Charge is smaller than "
778 << gkChargeRelErrLimit << "* its error in Pixel " << fPixId << endl;
779 return kFALSE;
780 }
781
782 if (!fHist->IsChargeFitOK())
783 {
784 *fLog << warn << "WARNING: Probability of Fitted Charge too low in Pixel "
785 << fPixId << endl;
786 return kFALSE;
787 }
788
789 if (fSigmaCharge < pedestal)
790 {
791 *fLog << warn << "WARNING: Sigma of Fitted Charge smaller than Pedestal RMS in Pixel "
792 << fPixId << endl;
793 return kFALSE;
794 }
795 return kTRUE;
796}
797
798//
799// The check return kTRUE if:
800//
801// 0) No value is nan
802// 1) Pixel has a fitted rel. time smaller than 3*FADC slices
803// 2) Pixel has a fit error greater than 0.
804// 4) Pixel has a fit Probability greater than 0.001
805// 5) The absolute arrival time is at least 1.0 slices from the used edge slices
806//
807Bool_t MCalibrationPix::CheckTimeFitValidity()
808{
809
810
811 if (IsExcludeQualityCheck())
812 return kTRUE;
813
814 if (IsHiGainSaturation())
815 {
816
817 if (fAbsTimeMean < (Float_t)fTimeFirstLoGain+1)
818 {
819 *fLog << warn
820 << "WARNING: Some absolute times smaller than limit in Pixel "
821 << fPixId << " time: " << fAbsTimeMean
822 << " Limit: " << (Float_t)fTimeFirstLoGain+1. << endl;
823 return kFALSE;
824 }
825
826 if (fAbsTimeMean > (Float_t)fTimeLastLoGain-1)
827 {
828 *fLog << warn
829 << "WARNING: Some absolute times bigger than limit in Pixel "
830 << fPixId << " time: " << fAbsTimeMean
831 << " Limit: " << (Float_t)fTimeLastLoGain-1. << endl;
832 return kFALSE;
833 }
834
835 }
836 else
837 {
838
839 if (fAbsTimeMean < (Float_t)fTimeFirstHiGain+1.)
840 {
841 *fLog << warn
842 << "WARNING: Some absolute times smaller than limit in Pixel "
843 << fPixId << " time: " << fAbsTimeMean
844 << " Limit: " << (Float_t)fTimeFirstHiGain+1. << endl;
845 // return kFALSE;
846 }
847
848 if (fAbsTimeMean > (Float_t)fTimeLastHiGain-1.)
849 {
850 *fLog << warn
851 << "WARNING: Some absolute times bigger than limit in Pixel "
852 << fPixId << " time: " << fAbsTimeMean
853 << " Limit: " << (Float_t)fTimeLastHiGain-1. << endl;
854 // return kFALSE;
855 }
856
857 }
858
859
860
861 return kTRUE;
862}
863
864
865void MCalibrationPix::CheckOscillations()
866{
867 fHist->CheckOscillations();
868}
869
870void MCalibrationPix::ApplyLoGainConversion()
871{
872
873 const Float_t chargeRelErrSquare = fChargeErr * fChargeErr
874 /( fCharge * fCharge );
875 const Float_t sigmaRelErrSquare = fSigmaChargeErr * fSigmaChargeErr
876 /( fSigmaCharge * fSigmaCharge );
877 const Float_t conversionRelErrSquare = fConversionHiLoErr * fConversionHiLoErr
878 /( fConversionHiLo * fConversionHiLo );
879
880 fCharge *= fConversionHiLo;
881 fChargeErr = TMath::Sqrt(chargeRelErrSquare + conversionRelErrSquare) * fCharge;
882
883 fSigmaCharge *= fConversionHiLo;
884 fSigmaChargeErr = TMath::Sqrt(sigmaRelErrSquare + conversionRelErrSquare) * fSigmaCharge;
885
886}
Note: See TracBrowser for help on using the repository browser.