source: trunk/MagicSoft/Mars/mcalib/MCalibrationChargePix.cc@ 3337

Last change on this file since 3337 was 3288, checked in by gaug, 21 years ago
*** empty log message ***
File size: 26.4 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 02/2004 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26// //
27// MCalibrationChargePix //
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 "MCalibrationChargePix.h"
88
89#include "MLog.h"
90#include "MLogManip.h"
91
92ClassImp(MCalibrationChargePix);
93
94using namespace std;
95
96const Float_t MCalibrationChargePix::gkElectronicPedRms = 1.5;
97const Float_t MCalibrationChargePix::gkElectronicPedRmsErr = 0.3;
98const Float_t MCalibrationChargePix::gkFFactor = 1.15;
99const Float_t MCalibrationChargePix::gkFFactorErr = 0.02;
100
101const Float_t MCalibrationChargePix::gkAverageQE = 0.18;
102const Float_t MCalibrationChargePix::gkAverageQEErr = 0.02;
103const Float_t MCalibrationChargePix::gkConversionHiLo = 10.;
104const Float_t MCalibrationChargePix::gkConversionHiLoErr = 2.5;
105
106const Float_t MCalibrationChargePix::fgChargeLimit = 3.;
107const Float_t MCalibrationChargePix::fgChargeErrLimit = 0.;
108const Float_t MCalibrationChargePix::fgChargeRelErrLimit = 1.;
109const Float_t MCalibrationChargePix::fgTimeLimit = 1.5;
110const Float_t MCalibrationChargePix::fgTimeErrLimit = 3.;
111const Float_t MCalibrationChargePix::fgConvFFactorRelErrLimit = 0.25;
112// --------------------------------------------------------------------------
113//
114// Default Constructor:
115//
116MCalibrationChargePix::MCalibrationChargePix(const char *name, const char *title)
117 : fPixId(-1),
118 fFlags(0)
119{
120
121 fName = name ? name : "MCalibrationChargePix";
122 fTitle = title ? title : "Container of the fit results of MHCalibrationChargePixs ";
123
124 Clear();
125
126 //
127 // At the moment, we don't have a database, yet,
128 // so we get it from the configuration file
129 //
130 SetConversionHiLo();
131 SetConversionHiLoErr();
132
133 SetAverageQE();
134 SetChargeLimit();
135 SetChargeErrLimit();
136
137 SetChargeRelErrLimit();
138 SetTimeLimit();
139 SetTimeErrLimit();
140 SetConvFFactorRelErrLimit();
141}
142
143// ------------------------------------------------------------------------
144//
145// Invalidate values
146//
147void MCalibrationChargePix::Clear(Option_t *o)
148{
149
150 SetHiGainSaturation ( kFALSE );
151 SetLoGainSaturation ( kFALSE );
152 SetHiGainFitted ( kFALSE );
153 SetLoGainFitted ( kFALSE );
154 SetHiGainOscillating ( kFALSE );
155 SetLoGainOscillating ( kFALSE );
156 SetExcluded ( kFALSE );
157 SetChargeValid ( kFALSE );
158 SetBlindPixelMethodValid ( kFALSE );
159 SetFFactorMethodValid ( kFALSE );
160 SetPINDiodeMethodValid ( kFALSE );
161 SetCombinedMethodValid ( kFALSE );
162
163 SetMeanTimeInFirstBin ( kFALSE );
164 SetMeanTimeInLastBin ( kFALSE );
165
166 fHiGainMeanCharge = -1.;
167 fHiGainMeanChargeErr = -1.;
168 fHiGainSigmaCharge = -1.;
169 fHiGainSigmaChargeErr = -1.;
170 fHiGainChargeProb = -1.;
171
172 fLoGainMeanCharge = -1.;
173 fLoGainMeanChargeErr = -1.;
174 fLoGainSigmaCharge = -1.;
175 fLoGainSigmaChargeErr = -1.;
176 fLoGainChargeProb = -1.;
177
178 fRSigmaCharge = -1.;
179 fRSigmaChargeErr = -1.;
180
181 fHiGainNumPickup = -1;
182 fLoGainNumPickup = -1;
183
184 fNumLoGainSamples = -1.;
185
186 fPed = -1.;
187 fPedRms = -1.;
188 fPedErr = -1.;
189
190 fLoGainPedRms = -1.;
191 fLoGainPedRmsErr = -1.;
192
193 fTimeFirstHiGain = 0 ;
194 fTimeLastHiGain = 0 ;
195 fTimeFirstLoGain = 0 ;
196 fTimeLastLoGain = 0 ;
197
198 fAbsTimeMean = -1.;
199 fAbsTimeRms = -1.;
200
201 fPheFFactorMethod = -1.;
202 fPheFFactorMethodErr = -1.;
203
204 fMeanConversionFFactorMethod = -1.;
205 fMeanConversionBlindPixelMethod = -1.;
206 fMeanConversionPINDiodeMethod = -1.;
207 fMeanConversionCombinedMethod = -1.;
208
209 fConversionFFactorMethodErr = -1.;
210 fConversionBlindPixelMethodErr = -1.;
211 fConversionPINDiodeMethodErr = -1.;
212 fConversionCombinedMethodErr = -1.;
213
214 fSigmaConversionFFactorMethod = -1.;
215 fSigmaConversionBlindPixelMethod = -1.;
216 fSigmaConversionPINDiodeMethod = -1.;
217 fSigmaConversionCombinedMethod = -1.;
218
219 fTotalFFactorFFactorMethod = -1.;
220 fTotalFFactorBlindPixelMethod = -1.;
221 fTotalFFactorPINDiodeMethod = -1.;
222 fTotalFFactorCombinedMethod = -1.;
223
224}
225
226
227void MCalibrationChargePix::DefinePixId(Int_t i)
228{
229 fPixId = i;
230}
231
232
233// --------------------------------------------------------------------------
234//
235// Set the pedestals from outside
236//
237void MCalibrationChargePix::SetPedestal(const Float_t ped, const Float_t pedrms, const Float_t pederr)
238{
239
240 fPed = ped;
241 fPedRms = pedrms;
242 fPedErr = pederr;
243}
244
245
246// --------------------------------------------------------------------------
247//
248// Set the conversion factors from outside (only for MC)
249//
250void MCalibrationChargePix::SetConversionFFactorMethod(Float_t c, Float_t err, Float_t sig)
251{
252 fMeanConversionFFactorMethod = c;
253 fConversionFFactorMethodErr = err;
254 fSigmaConversionFFactorMethod = sig;
255}
256
257// --------------------------------------------------------------------------
258//
259// Set the conversion factors from outside (only for MC)
260//
261void MCalibrationChargePix::SetConversionCombinedMethod(Float_t c, Float_t err, Float_t sig)
262{
263 fMeanConversionCombinedMethod = c;
264 fConversionCombinedMethodErr = err;
265 fSigmaConversionCombinedMethod = sig;
266}
267
268
269// --------------------------------------------------------------------------
270//
271// Set the conversion factors from outside (only for MC)
272//
273void MCalibrationChargePix::SetConversionBlindPixelMethod(Float_t c, Float_t err, Float_t sig)
274{
275 fMeanConversionBlindPixelMethod = c;
276 fConversionBlindPixelMethodErr = err;
277 fSigmaConversionBlindPixelMethod = sig;
278}
279
280// --------------------------------------------------------------------------
281//
282// Set the conversion factors from outside (only for MC)
283//
284void MCalibrationChargePix::SetConversionPINDiodeMethod(Float_t c, Float_t err, Float_t sig)
285{
286 fMeanConversionPINDiodeMethod = c ;
287 fConversionPINDiodeMethodErr = err;
288 fSigmaConversionPINDiodeMethod = sig;
289}
290
291// --------------------------------------------------------------------------
292//
293// Set the Hi Gain Saturation Bit from outside
294//
295void MCalibrationChargePix::SetHiGainSaturation(Bool_t b)
296{
297 b ? SETBIT(fFlags, kHiGainSaturation) : CLRBIT(fFlags, kHiGainSaturation);
298}
299
300// --------------------------------------------------------------------------
301//
302// Set the Lo Gain Saturation Bit from outside
303//
304void MCalibrationChargePix::SetLoGainSaturation(Bool_t b)
305{
306 b ? SETBIT(fFlags, kLoGainSaturation) : CLRBIT(fFlags, kLoGainSaturation);
307}
308
309// --------------------------------------------------------------------------
310//
311// Set the Excluded Bit from outside
312//
313void MCalibrationChargePix::SetExcluded(Bool_t b )
314{
315 b ? SETBIT(fFlags, kExcluded) : CLRBIT(fFlags, kExcluded);
316}
317
318// --------------------------------------------------------------------------
319//
320// Set the Excluded Bit from outside
321//
322void MCalibrationChargePix::SetChargeValid(Bool_t b )
323{
324 b ? SETBIT(fFlags, kChargeValid) : CLRBIT(fFlags, kChargeValid);
325}
326
327// --------------------------------------------------------------------------
328//
329// Set the Fitted Bit from outside
330//
331void MCalibrationChargePix::SetHiGainFitted(Bool_t b )
332{
333 b ? SETBIT(fFlags, kHiGainFitted) : CLRBIT(fFlags, kHiGainFitted);
334}
335
336// --------------------------------------------------------------------------
337//
338// Set the Fitted Bit from outside
339//
340void MCalibrationChargePix::SetLoGainFitted(const Bool_t b )
341{
342 b ? SETBIT(fFlags, kLoGainFitted) : CLRBIT(fFlags, kLoGainFitted);
343}
344
345// --------------------------------------------------------------------------
346//
347// Set the Ocillating Bit from outside
348//
349void MCalibrationChargePix::SetHiGainOscillating(const Bool_t b )
350{
351 b ? SETBIT(fFlags, kHiGainOscillating) : CLRBIT(fFlags, kHiGainOscillating);
352}
353
354// --------------------------------------------------------------------------
355//
356// Set the Ocillating Bit from outside
357//
358void MCalibrationChargePix::SetLoGainOscillating(const Bool_t b )
359{
360 b ? SETBIT(fFlags, kLoGainOscillating) : CLRBIT(fFlags, kLoGainOscillating);
361}
362
363// --------------------------------------------------------------------------
364//
365// Set the Excluded Bit from outside
366//
367void MCalibrationChargePix::SetBlindPixelMethodValid(const Bool_t b )
368{
369 b ? SETBIT(fFlags, kBlindPixelMethodValid) : CLRBIT(fFlags, kBlindPixelMethodValid);
370}
371
372// --------------------------------------------------------------------------
373//
374// Set the Excluded Bit from outside
375//
376void MCalibrationChargePix::SetFFactorMethodValid(const Bool_t b )
377{
378 b ? SETBIT(fFlags, kFFactorMethodValid) : CLRBIT(fFlags, kFFactorMethodValid);
379}
380
381// --------------------------------------------------------------------------
382//
383// Set the Excluded Bit from outside
384//
385void MCalibrationChargePix::SetPINDiodeMethodValid(const Bool_t b )
386{
387 b ? SETBIT(fFlags, kPINDiodeMethodValid) : CLRBIT(fFlags, kPINDiodeMethodValid);
388}
389
390// --------------------------------------------------------------------------
391//
392// Set the Excluded Bit from outside
393//
394void MCalibrationChargePix::SetCombinedMethodValid(const Bool_t b )
395{
396 b ? SETBIT(fFlags, kCombinedMethodValid) : CLRBIT(fFlags, kCombinedMethodValid);
397}
398
399void MCalibrationChargePix::SetAbsTimeBordersHiGain(const Byte_t f, const Byte_t l)
400{
401
402 fTimeFirstHiGain = f;
403 fTimeLastHiGain = l;
404
405}
406
407void MCalibrationChargePix::SetAbsTimeBordersLoGain(const Byte_t f, const Byte_t l)
408{
409
410 fTimeFirstLoGain = f;
411 fTimeLastLoGain = l;
412
413}
414
415void MCalibrationChargePix::SetMeanTimeInFirstBin(const Bool_t b)
416{
417 b ? SETBIT(fFlags,kMeanTimeInFirstBin) : CLRBIT(fFlags,kMeanTimeInFirstBin);
418}
419
420void MCalibrationChargePix::SetMeanTimeInLastBin(const Bool_t b)
421{
422 b ? SETBIT(fFlags,kMeanTimeInLastBin) : CLRBIT(fFlags,kMeanTimeInLastBin);
423}
424
425Float_t MCalibrationChargePix::GetMeanCharge() const
426{
427 return IsHiGainSaturation() ? fLoGainMeanCharge : fHiGainMeanCharge ;
428}
429
430Float_t MCalibrationChargePix::GetMeanChargeErr() const
431{
432 return IsHiGainSaturation() ? fLoGainMeanChargeErr : fHiGainMeanChargeErr ;
433}
434
435Float_t MCalibrationChargePix::GetChargeProb() const
436{
437 return IsHiGainSaturation() ? fLoGainChargeProb : fHiGainChargeProb ;
438}
439
440Float_t MCalibrationChargePix::GetSigmaCharge() const
441{
442 return IsHiGainSaturation() ? fLoGainSigmaCharge : fHiGainSigmaCharge ;
443}
444
445Float_t MCalibrationChargePix::GetSigmaChargeErr() const
446{
447 return IsHiGainSaturation() ? fLoGainSigmaChargeErr : fHiGainSigmaChargeErr ;
448}
449
450Bool_t MCalibrationChargePix::IsFitted() const
451{
452 return IsHiGainSaturation() ? IsLoGainFitted() : IsHiGainFitted();
453}
454
455Bool_t MCalibrationChargePix::IsOscillating() const
456{
457 return IsHiGainSaturation() ? IsLoGainOscillating() : IsHiGainOscillating();
458}
459
460Bool_t MCalibrationChargePix::IsMeanTimeInFirstBin() const
461{
462 return TESTBIT(fFlags,kMeanTimeInFirstBin);
463}
464
465Bool_t MCalibrationChargePix::IsMeanTimeInLastBin() const
466{
467 return TESTBIT(fFlags,kMeanTimeInLastBin);
468}
469
470
471
472Bool_t MCalibrationChargePix::IsExcluded() const
473{
474 return TESTBIT(fFlags,kExcluded);
475}
476
477Bool_t MCalibrationChargePix::IsHiGainSaturation() const
478{
479 return TESTBIT(fFlags,kHiGainSaturation);
480}
481
482Bool_t MCalibrationChargePix::IsLoGainSaturation() const
483{
484 return TESTBIT(fFlags,kLoGainSaturation);
485}
486
487Bool_t MCalibrationChargePix::IsChargeValid() const
488{
489 return TESTBIT(fFlags, kChargeValid);
490}
491
492Bool_t MCalibrationChargePix::IsHiGainFitted() const
493{
494 return TESTBIT(fFlags, kHiGainFitted);
495}
496
497Bool_t MCalibrationChargePix::IsLoGainFitted() const
498{
499 return TESTBIT(fFlags, kLoGainFitted);
500}
501
502Bool_t MCalibrationChargePix::IsHiGainOscillating() const
503{
504 return TESTBIT(fFlags, kHiGainOscillating);
505}
506
507Bool_t MCalibrationChargePix::IsLoGainOscillating() const
508{
509 return TESTBIT(fFlags, kLoGainOscillating);
510}
511
512Bool_t MCalibrationChargePix::IsBlindPixelMethodValid() const
513{
514 return TESTBIT(fFlags, kBlindPixelMethodValid);
515}
516
517Bool_t MCalibrationChargePix::IsFFactorMethodValid() const
518{
519 return TESTBIT(fFlags, kFFactorMethodValid);
520}
521
522Bool_t MCalibrationChargePix::IsPINDiodeMethodValid() const
523{
524 return TESTBIT(fFlags, kPINDiodeMethodValid);
525}
526
527Bool_t MCalibrationChargePix::IsCombinedMethodValid() const
528{
529 return TESTBIT(fFlags, kCombinedMethodValid);
530}
531
532
533//
534// The check return kTRUE if:
535//
536// 1) Pixel has a fitted charge greater than fChargeLimit*PedRMS
537// 2) Pixel has a fit error greater than fChargeErrLimit
538// 3) Pixel has a fitted charge greater its fChargeRelErrLimit times its charge error
539// 4) Pixel has a charge sigma bigger than its Pedestal RMS
540//
541Bool_t MCalibrationChargePix::CheckChargeValidity()
542{
543
544 if (!IsFitted())
545 return kFALSE;
546
547 if (GetMeanCharge() < fChargeLimit*GetPedRms())
548 {
549 *fLog << warn << "WARNING: Fitted Charge is smaller than "
550 << fChargeLimit << " Pedestal RMS in Pixel " << fPixId << endl;
551 return kFALSE;
552 }
553
554 if (GetMeanChargeErr() < fChargeErrLimit)
555 {
556 *fLog << warn << "WARNING: Error of Fitted Charge is smaller than "
557 << fChargeErrLimit << " in Pixel " << fPixId << endl;
558 return kFALSE;
559 }
560
561 if (GetMeanCharge() < fChargeRelErrLimit*GetMeanChargeErr())
562 {
563 *fLog << warn << "WARNING: Fitted Charge is smaller than "
564 << fChargeRelErrLimit << "* its error in Pixel " << fPixId << endl;
565 return kFALSE;
566 }
567
568 if (GetSigmaCharge() < GetPedRms())
569 {
570 *fLog << warn << "WARNING: Sigma of Fitted Charge smaller than Pedestal RMS in Pixel " << fPixId << endl;
571 return kFALSE;
572 }
573
574 SetChargeValid();
575 return kTRUE;
576}
577
578//
579// The check returns kTRUE if:
580//
581// The mean arrival time is at least 1.0 slices from the used edge slices
582//
583Bool_t MCalibrationChargePix::CheckTimeValidity()
584{
585
586 const Float_t loweredge = IsHiGainSaturation() ? fTimeFirstHiGain : fTimeFirstLoGain;
587 const Float_t upperedge = IsHiGainSaturation() ? fTimeLastHiGain : fTimeLastLoGain;
588
589 if ( fAbsTimeMean < loweredge+1.)
590 {
591 *fLog << warn << "WARNING: Mean ArrivalTime in first extraction bin of the Pixel " << fPixId << endl;
592 SetMeanTimeInFirstBin();
593 return kFALSE;
594 }
595
596 if ( fAbsTimeMean > upperedge-1.)
597 {
598 *fLog << warn << "WARNING: Mean ArrivalTime in last extraction bin of the Pixel " << fPixId << endl;
599 SetMeanTimeInLastBin();
600 return kFALSE;
601 }
602
603 return kTRUE;
604}
605
606void MCalibrationChargePix::CalcLoGainPed()
607{
608
609 Float_t pedRmsSquare = fPedRms * fPedRms;
610 Float_t pedRmsSquareErrSquare = fPedErr * fPedErr * pedRmsSquare; // fPedRmsErr = fPedErr/2.
611 Float_t pedRmsSquareErr;
612
613 //
614 // We do not know the Lo Gain Pedestal RMS, so we have to retrieve it
615 // from the HI GAIN (all calculation per slice up to now):
616 //
617 // We extract the pure NSB contribution:
618 //
619 const Float_t elecRmsSquare = fElectronicPedRms * fElectronicPedRms;
620 const Float_t elecRmsSquareErrSquare = 4.*fElectronicPedRmsErr * fElectronicPedRmsErr * elecRmsSquare;
621
622 Float_t nsbSquare = pedRmsSquare - elecRmsSquare;
623 Float_t nsbSquareRelErrSquare = (pedRmsSquareErrSquare + elecRmsSquareErrSquare)
624 / (nsbSquare * nsbSquare) ;
625
626 if (nsbSquare < 0.)
627 nsbSquare = 0.;
628
629 //
630 // Now, we divide the NSB by the conversion factor and
631 // add it quadratically to the electronic noise
632 //
633 const Float_t conversionSquare = fConversionHiLo * fConversionHiLo;
634 const Float_t conversionSquareRelErrSquare = 4.*fConversionHiLoErr * fConversionHiLoErr / conversionSquare;
635
636 const Float_t convertedNsbSquare = nsbSquare / conversionSquare;
637 const Float_t convertedNsbSquareErrSquare = (nsbSquareRelErrSquare + conversionSquareRelErrSquare)
638 * convertedNsbSquare * convertedNsbSquare;
639
640 pedRmsSquare = convertedNsbSquare + elecRmsSquare;
641 pedRmsSquareErr = TMath::Sqrt( convertedNsbSquareErrSquare + elecRmsSquareErrSquare );
642
643 //
644 // Correct also for the conversion to Hi-Gain:
645 //
646 fLoGainPedRms = TMath::Sqrt(pedRmsSquare) * fConversionHiLo;
647 fLoGainPedRmsErr = 0.5 * pedRmsSquareErr / fLoGainPedRms * fConversionHiLo;
648}
649
650Bool_t MCalibrationChargePix::CalcReducedSigma()
651{
652
653 const Float_t sigmaSquare = GetSigmaCharge() * GetSigmaCharge();
654 const Float_t sigmaSquareErrSquare = 4.*GetSigmaChargeErr()* GetSigmaChargeErr() * sigmaSquare;
655
656 Float_t pedRmsSquare;
657 Float_t pedRmsSquareErrSquare;
658
659 if (IsHiGainSaturation())
660 {
661 CalcLoGainPed();
662
663 pedRmsSquare = fLoGainPedRms * fLoGainPedRms;
664 pedRmsSquareErrSquare = 4.* fLoGainPedRmsErr * fLoGainPedRmsErr * pedRmsSquare;
665 } /* if (HiGainSaturation) */
666 else
667 {
668
669 pedRmsSquare = fPedRms * fPedRms;
670 pedRmsSquareErrSquare = fPedErr * fPedErr * pedRmsSquare; // fPedRmsErr = fPedErr/2.
671 }
672 //
673 // Calculate the reduced sigmas
674 //
675 const Float_t rsigmachargesquare = sigmaSquare - pedRmsSquare;
676
677 if (rsigmachargesquare <= 0.)
678 {
679 *fLog << warn
680 << "WARNING: Cannot calculate the reduced sigma: smaller than 0 in pixel "
681 << fPixId << endl;
682 return kFALSE;
683 }
684
685 fRSigmaCharge = TMath::Sqrt(rsigmachargesquare);
686 fRSigmaChargeErr = TMath::Sqrt(sigmaSquareErrSquare + pedRmsSquareErrSquare) / 2. / fRSigmaCharge;
687
688 return kTRUE;
689}
690
691//
692// Calculate the number of photo-electrons after the F-Factor method
693// Calculate the errors of the F-Factor method
694//
695Bool_t MCalibrationChargePix::CalcFFactorMethod()
696{
697
698 if (!IsChargeValid())
699 return kFALSE;
700
701 if (fRSigmaCharge < 0.)
702 return kFALSE;
703
704 //
705 // Square all variables in order to avoid applications of square root
706 //
707 // First the relative error squares
708 //
709 const Float_t chargeSquare = GetMeanCharge() * GetMeanCharge();
710 const Float_t chargeSquareRelErrSquare = 4.* GetMeanChargeErr() * GetMeanChargeErr() / chargeSquare;
711
712 const Float_t chargeRelErrSquare = GetMeanChargeErr() * GetMeanChargeErr()
713 / (GetMeanCharge() * GetMeanCharge());
714
715 const Float_t ffactorsquare = gkFFactor * gkFFactor;
716 const Float_t ffactorsquareRelErrSquare = 4.*gkFFactorErr * gkFFactorErr / ffactorsquare;
717
718 const Float_t avQERelErrSquare = fAverageQEErr * fAverageQEErr / fAverageQE / fAverageQE;
719
720 const Float_t avQEFFactor = TMath::Sqrt( ( 1. - fAverageQE ) / fAverageQE );
721 const Float_t avQEFFactorErr = 1./ ( 2. * avQEFFactor ) * fAverageQEErr
722 / ( fAverageQE * fAverageQE );
723
724 const Float_t avQEFFactorRelErrSquare = avQEFFactorErr * avQEFFactorErr
725 / ( avQEFFactor * avQEFFactor) ;
726
727 const Float_t rsigmaSquare = fRSigmaCharge * fRSigmaCharge;
728 const Float_t rsigmaSquareRelErrSquare = 4.* fRSigmaChargeErr * fRSigmaChargeErr / rsigmaSquare;
729
730
731 //
732 // Calculate the number of phe's from the F-Factor method
733 // (independent on Hi Gain or Lo Gain)
734 //
735 fPheFFactorMethod = ffactorsquare * chargeSquare / rsigmaSquare;
736
737 //
738 // Calculate the Error of Nphe
739 //
740 const Float_t pheFFactorRelErrSquare = ffactorsquareRelErrSquare
741 + chargeSquareRelErrSquare
742 + rsigmaSquareRelErrSquare;
743 fPheFFactorMethodErr = TMath::Sqrt(pheFFactorRelErrSquare) * fPheFFactorMethod;
744
745 //
746 // Calculate the conversion factor between PHOTONS and FADC counts
747 //
748 // Nphot = Nphe / avQE
749 // conv = Nphot / FADC counts
750 //
751 fMeanConversionFFactorMethod = fPheFFactorMethod / fAverageQE / GetMeanCharge();
752
753 if (fMeanConversionFFactorMethod < 0.)
754 return kFALSE;
755
756 //
757 // Calculate the error of the mean conversion factor between PHOTONS and FADC counts
758 //
759 const Float_t convRelErrSquare = ( pheFFactorRelErrSquare + chargeRelErrSquare + avQERelErrSquare);
760
761 if (convRelErrSquare < 0.)
762 return kFALSE;
763
764
765 const Float_t convRelErr = TMath::Sqrt(convRelErrSquare);
766 fConversionFFactorMethodErr = convRelErr * fMeanConversionFFactorMethod;
767
768 if (convRelErr < fConvFFactorRelErrLimit)
769 SetFFactorMethodValid();
770
771 //
772 // Calculate the Total F-Factor of the camera (in photons)
773 //
774 fTotalFFactorFFactorMethod = (fRSigmaCharge/GetMeanCharge())*TMath::Sqrt(fPheFFactorMethod);
775 fTotalFFactorFFactorMethod *= avQEFFactor;
776
777 //
778 // Calculate the error of the Total F-Factor of the camera ( in photons )
779 //
780 const Float_t rsigmaChargeRelErrSquare = fRSigmaChargeErr * fRSigmaChargeErr
781 / (fRSigmaCharge * fRSigmaCharge) ;
782
783 fTotalFFactorErrFFactorMethod = TMath::Sqrt( rsigmaChargeRelErrSquare
784 + chargeRelErrSquare
785 + pheFFactorRelErrSquare
786 + avQEFFactorRelErrSquare );
787 fTotalFFactorErrFFactorMethod *= fTotalFFactorFFactorMethod;
788
789 //
790 // Calculate the sigma of the conversion from FADC counts to photons
791 //
792 fSigmaConversionFFactorMethod = GetTotalFFactorFFactorMethod()*TMath::Sqrt(fMeanConversionFFactorMethod);
793
794 return kTRUE;
795}
796
797void MCalibrationChargePix::ApplyLoGainConversion()
798{
799
800 const Float_t chargeRelErrSquare = fLoGainMeanChargeErr * fLoGainMeanChargeErr
801 /( fLoGainMeanCharge * fLoGainMeanCharge );
802 const Float_t sigmaRelErrSquare = fLoGainSigmaChargeErr * fLoGainSigmaChargeErr
803 /( fLoGainSigmaCharge * fLoGainSigmaCharge );
804 const Float_t conversionRelErrSquare = fConversionHiLoErr * fConversionHiLoErr
805 /( fConversionHiLo * fConversionHiLo );
806
807 fLoGainMeanCharge *= fConversionHiLo;
808 fLoGainMeanChargeErr = TMath::Sqrt(chargeRelErrSquare + conversionRelErrSquare) * fLoGainMeanCharge;
809
810 fLoGainSigmaCharge *= fConversionHiLo;
811 fLoGainSigmaChargeErr = TMath::Sqrt(sigmaRelErrSquare + conversionRelErrSquare) * fLoGainSigmaCharge;
812
813 fElectronicPedRms = gkElectronicPedRms * TMath::Sqrt(fNumLoGainSamples);
814 fElectronicPedRmsErr = gkElectronicPedRmsErr * TMath::Sqrt(fNumLoGainSamples);
815
816}
817
Note: See TracBrowser for help on using the repository browser.