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

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