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

Last change on this file since 9212 was 8917, checked in by tbretz, 17 years ago
*** empty log message ***
File size: 38.1 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-2008
21!
22\* ======================================================================== */
23
24/////////////////////////////////////////////////////////////////////////////
25//
26// MCalibrationChargePix
27//
28// Storage container of the calibrated Charge of one pixel.
29//
30// The following values are initialized to meaningful values:
31//
32// - The Electronic Rms to 1.5 per FADC slice
33// - The uncertainty about the Electronic RMS to 0.3 per slice
34// - The F-Factor is assumed to have been measured in Munich to 1.13 - 1.17.
35// with the Munich definition of the F-Factor, thus:
36// F = Sigma(Out)/Mean(Out) * Mean(In)/Sigma(In)
37// Mean F-Factor (gkFFactor) = 1.15
38// Error F-Factor (gkFFactorErr) = 0.02
39//
40// The following variables are calculated inside this class:
41// - fLoGainPedRmsSquare and fLoGainPedRmsSquareVar (see CalcLoGainPedestal())
42// - fRSigmaSquare and fRSigmaSquareVar (see CalcReducedSigma() )
43// - fPheFFactorMethod and fPheFFactorMethodVar (see CalcFFactor() )
44// - fMeanConvFADC2Phe and fMeanConvFADC2PheVar (see CalcConvFFactor() )
45//
46// The following variables are set by MHCalibrationChargeCam:
47// - fAbsTimeMean and fAbsTimeRms
48// - all variables in MCalibrationPix
49//
50// The following variables are set by MCalibrationChargeCalc:
51// - fPed, fPedVar and fPedRms
52// - fMeanConvFADC2Phe
53// - fConvFADC2PheVar
54// - fSigmaConvFADC2Phe
55// - fTotalFFactorFFactorMethod
56// - fTotalFFactorFFactorMethodVar
57//
58// The following variables are not yet implemented:
59// - fConversionHiLo and fConversionHiLoVar (now set fixed to 10. +- 2.5)
60//
61// Error of all variables are calculated by error-propagation. Note that internally,
62// all error variables contain Variances in order to save the CPU-intensive square rooting
63//
64// Low-Gain variables are stored internally unconverted, i.e. directly from the summed
65// FADC slices extraction results, but can be retrieved converted to High-Gain amplifications
66// by calls to: GetConvertedMean() or GetConvertedSigma()
67//
68// Inline Functions:
69// -----------------
70//
71// GetPedRms(): Get the pedestals RMS: Test bit kHiGainSaturation:
72// If yes, return square root of fLoGainPedRmsSquare (if greater than 0,
73// otherwise -1.), If no, return fPedRms
74//
75// GetConvertedMean(): Get the Low Gain Mean Charge converted to High Gain
76// amplification: Returns fLoGainMean multiplied with fConversionHiLo if
77// IsHiGainSaturation(), else return fHiGainMean
78//
79// GetConvertedSigma(): Get the Low Gain Sigma converted to High Gain
80// amplification: Returns fLoGainSigma multiplied with fConversionHiLo
81// if IsHiGainSaturation() else return fHiGainSigma
82//
83// GetConvertedRSigmaSquare(): Get the reduced Sigma Square:
84// If fRSigmaSquare is smaller than 0 (i.e. has not yet been set), return -1.
85// Test bit kHiGainSaturation:
86// If yes, return fRSigmaSquare, multiplied with fConversionHiLo^2,
87// If no , return fRSigmaSquare
88//
89// GetPheFFactorMethodRelVar(): Get the relative variance on the number of
90// photo-electrons (F-Factor Method): If fPheFFactorMethodVar is smaller
91// than 0 (i.e. has not yet been set), return -1. If fPheFFactorMethod
92// is 0, return -1. Else returns fPheFFactorMethodVar / fPheFFactorMethod^2
93//
94// GetMeanConvFADC2PheErr(): Get the error on the mean conversion factor
95// (FFactor  Method): If fMeanConvFADC2PheVar is smaller than 0 (i.e. has
96// not yet been set), return -1. Else returns the square root of
97// fMeanConvFADC2PheVar
98//
99//
100// See also: MCalibrationChargeCam, MCalibrationChargeCalc,
101// MHCalibrationChargeCam, MHCalibrationChargePix
102//
103//
104// ClassVersion 3:
105// + fNumSaturated
106//
107// Class Version 4:
108// + Float_t fConversionHiLoSigma; // Sigma of conversion factor betw. Hi and Lo Gain
109// - Float_t fMeanConvFADC2PheVar; // Variance conversion factor (F-factor method)
110// + Float_t fMeanConvFADC2PheStatVar; // Variance conversion factor, only stat. error
111// + Float_t fMeanConvFADC2PheSystVar; // Variance conversion factor, only syst. error
112// - Float_t fMeanFFactorFADC2PhotVar; // Variance mean F-Factor photons (F-factor method)
113// + Float_t fMeanFFactorFADC2PhotVar; // Variance mean F-Factor photons, only stat. error
114// - Float_t fPheFFactorMethod; // Number Phe's calculated (F-factor method)
115// - Float_t fPheFFactorMethodVar; // Variance number of Phe's (F-factor method)
116// + Float_t fPheFFactorMethod; // Number Phe's calculated with F-factor method)
117// + Float_t fPheFFactorMethodStatVar; // Variance number of Phe's, only stat. error
118// + Float_t fPheFFactorMethodSystVar; // Variance number of Phe's, only syst. error
119//
120/////////////////////////////////////////////////////////////////////////////
121#include "MCalibrationChargePix.h"
122
123#include <TMath.h>
124
125#include "MLog.h"
126#include "MLogManip.h"
127
128#include "MBadPixelsPix.h"
129
130ClassImp(MCalibrationChargePix);
131
132using namespace std;
133
134const Float_t MCalibrationChargePix::gkElectronicPedRms = 1.0;
135const Float_t MCalibrationChargePix::gkElectronicPedRmsErr = 0.35;
136const Float_t MCalibrationChargePix::gkFFactor = 1.15;
137const Float_t MCalibrationChargePix::gkFFactorErr = 0.02;
138
139const Float_t MCalibrationChargePix::fgConversionHiLo = 10.;
140const Float_t MCalibrationChargePix::fgConversionHiLoErr = 0.05;
141const Float_t MCalibrationChargePix::fgConversionHiLoSigma = 2.5;
142const Float_t MCalibrationChargePix::fgPheFFactorMethodLimit = 1.;
143const Float_t MCalibrationChargePix::fgConvFFactorRelErrLimit = 0.85;
144
145// --------------------------------------------------------------------------
146//
147// Default Constructor:
148//
149// Sets:
150// - fCalibFlags to 0
151// - fConversionHiLo to fgConversionHiLo
152// - fConversionHiLoVar to square of fgConversionHiLoErr
153// - fConvFFactorelErrLimit to fgConvFFactorRelErrLimit*fgConvFFactorelErrLimit
154// - fPheFFactorLimit to fgPheFFactorLimit
155//
156// Calls:
157// - Clear()
158//
159MCalibrationChargePix::MCalibrationChargePix(const char *name, const char *title)
160 : fCalibFlags(0), fNumSaturated(0)
161{
162
163 fName = name ? name : "MCalibrationChargePix";
164 fTitle = title ? title : "Container of the fit results of MHCalibrationChargePixs ";
165
166 //
167 // At the moment, we don't have a database, yet,
168 // so we get it from the configuration file
169 //
170 SetConversionHiLo();
171 SetConversionHiLoErr();
172
173 SetPheFFactorMethodLimit();
174 SetConvFFactorRelErrLimit();
175
176 Clear();
177}
178
179// ------------------------------------------------------------------------
180//
181// Sets:
182// - all flags to kFALSE
183// - all variables to -1.
184//
185// Calls:
186// - MCalibrationPix::Clear()
187//
188void MCalibrationChargePix::Clear(Option_t *o)
189{
190
191 SetFFactorMethodValid ( kFALSE );
192
193 fRSigmaSquare = -1.;
194 fRSigmaSquareVar = -1.;
195
196 fPed = -1.;
197 fPedRms = -1.;
198 fPedVar = -1.;
199
200 fLoGainPedRmsSquare = -1.;
201 fLoGainPedRmsSquareVar = -1.;
202
203 fAbsTimeMean = -1.;
204 fAbsTimeRms = -1.;
205
206 fPheFFactorMethod = -1.;
207 fPheFFactorMethodStatVar = -1.;
208 fPheFFactorMethodSystVar = -1.;
209
210 fMeanConvFADC2Phe = -1.;
211 fMeanConvFADC2PheStatVar = -1.;
212 fMeanConvFADC2PheSystVar = -1.;
213 fMeanFFactorFADC2Phot = -1.;
214 fMeanFFactorFADC2PhotVar = -1.;
215
216 fNumSaturated = 0;
217
218 MCalibrationPix::Clear();
219}
220
221// -----------------------------------------------------
222//
223void MCalibrationChargePix::Copy(TObject& object) const
224{
225 MCalibrationChargePix &pix = (MCalibrationChargePix&)object;
226
227 MCalibrationPix::Copy(pix);
228
229 //
230 // Copy the data members
231 //
232 pix.fAbsTimeMean = fAbsTimeMean ;
233 pix.fAbsTimeRms = fAbsTimeRms ;
234 pix.fCalibFlags = fCalibFlags ;
235 pix.fConversionHiLo = fConversionHiLo ;
236 pix.fConversionHiLoVar = fConversionHiLoVar ;
237 pix.fConversionHiLoSigma = fConversionHiLoSigma ;
238 pix.fConvFFactorRelVarLimit = fConvFFactorRelVarLimit ;
239 pix.fLoGainPedRmsSquare = fLoGainPedRmsSquare ;
240 pix.fLoGainPedRmsSquareVar = fLoGainPedRmsSquareVar ;
241 pix.fMeanConvFADC2Phe = fMeanConvFADC2Phe ;
242 pix.fMeanConvFADC2PheStatVar = fMeanConvFADC2PheStatVar ;
243 pix.fMeanConvFADC2PheSystVar = fMeanConvFADC2PheSystVar ;
244 pix.fMeanFFactorFADC2Phot = fMeanFFactorFADC2Phot ;
245 pix.fMeanFFactorFADC2PhotVar = fMeanFFactorFADC2PhotVar ;
246 pix.fPed = fPed ;
247 pix.fPedVar = fPedVar ;
248 pix.fPedRms = fPedRms ;
249 pix.fPedRmsVar = fPedRmsVar ;
250 pix.fPheFFactorMethod = fPheFFactorMethod ;
251 pix.fPheFFactorMethodStatVar = fPheFFactorMethodStatVar ;
252 pix.fPheFFactorMethodSystVar = fPheFFactorMethodSystVar ;
253 pix.fPheFFactorMethodLimit = fPheFFactorMethodLimit ;
254 pix.fRSigmaSquare = fRSigmaSquare ;
255 pix.fRSigmaSquareVar = fRSigmaSquareVar ;
256 pix.fNumSaturated = fNumSaturated ;
257}
258
259
260// --------------------------------------------------------------------------
261//
262// Set F-Factor Method Validity Bit from outside
263//
264void MCalibrationChargePix::SetFFactorMethodValid(const Bool_t b )
265{
266 b ? SETBIT(fCalibFlags, kFFactorMethodValid) : CLRBIT(fCalibFlags, kFFactorMethodValid);
267}
268
269// --------------------------------------------------------------------------
270//
271// Set pedestals from outside (done by MCalibrationChargeCalc)
272//
273void MCalibrationChargePix::SetPedestal(const Float_t ped, const Float_t pedrms, const Float_t pederr)
274{
275 fPed = ped;
276 fPedRms = pedrms;
277 fPedVar = pederr*pederr;
278}
279
280// --------------------------------------------------------------------------
281//
282// Set pedestals from outside (done by MCalibrationChargeCalc)
283//
284void MCalibrationChargePix::SetPed(const Float_t ped, const Float_t pederr)
285{
286 fPed = ped;
287 fPedVar = pederr*pederr;
288}
289
290// --------------------------------------------------------------------------
291//
292// Set pedestals RMS from outside (done by MHCalibrationChargeCam)
293//
294void MCalibrationChargePix::SetPedRMS( const Float_t pedrms, const Float_t pedrmserr)
295{
296 fPedRms = pedrms;
297 fPedRmsVar = pedrmserr*pedrmserr;
298}
299
300// --------------------------------------------------------------------------
301//
302// Get the conversion Error Hi-Gain to Low-Gain:
303// If fConversionHiLoVar is smaller than 0 (i.e. has not yet
304// been set), return -1.
305//
306Float_t MCalibrationChargePix::GetConversionHiLoErr() const
307{
308 return fConversionHiLoVar<0 ? -1 : TMath::Sqrt(fConversionHiLoVar);
309}
310
311// --------------------------------------------------------------------------
312//
313// Get the Error of the converted Low Gain Mean: Returns -1 if the
314// variable fLoGainMean or fLoGainMeanVar are smaller than 0. Returns the
315// square root of the quadratic sum of the relative variances of the
316// fLoGainMean and fConversionHiLo, mulitplied with GetConvertedMean() in
317// case of HiGain Saturation, else return GetMeanErr()
318//
319Float_t MCalibrationChargePix::GetConvertedMeanErr() const
320{
321 if (!IsHiGainSaturation())
322 return GetMeanErr();
323
324 const Float_t logainrelvar = GetLoGainMeanRelVar();
325 return logainrelvar<0 ? -1 : TMath::Sqrt(logainrelvar + GetConversionHiLoRelVar()) * GetConvertedMean();
326}
327
328// --------------------------------------------------------------------------
329//
330// Get the Error of the converted Sigma: Returns -1 if the variable
331// fLoGainSigma or fLoGainSigmaVar are smaller than 0. if
332// IsHiGainSaturatio() returns the square root of the quadratic sum of the
333// relative variances of the fLoGainSigma and fConversionHiLo, mulitplied
334// with GetConvertedSigma() else returns GetSigmaErr()
335//
336Float_t MCalibrationChargePix::GetConvertedSigmaErr() const
337{
338 if (!IsHiGainSaturation())
339 return GetSigmaErr();
340
341 if (fLoGainSigmaVar<0 || fLoGainSigma<0)
342 return -1.;
343
344 const Float_t sigmaRelVar = fLoGainSigmaVar/(fLoGainSigma*fLoGainSigma);
345 return TMath::Sqrt(sigmaRelVar+GetConversionHiLoRelVar()) * GetConvertedSigma();
346}
347
348// --------------------------------------------------------------------------
349//
350// Get the converted reduced Sigma: If fRSigmaSquare is smaller than 0
351// (i.e. has not yet been set), return -1. Test bit kHiGainSaturation: If
352// yes, return square root of fRSigmaSquare, multiplied with
353// fConversionHiLo, If no , return square root of fRSigmaSquare
354//
355Float_t MCalibrationChargePix::GetConvertedRSigma() const
356{
357 if (fRSigmaSquare < 0)
358 return -1;
359
360 const Float_t rsigma = TMath::Sqrt(fRSigmaSquare);
361 return IsHiGainSaturation() ? rsigma*fConversionHiLo : rsigma ;
362}
363
364// --------------------------------------------------------------------------
365//
366Float_t MCalibrationChargePix::GetMeanConvFADC2PheErr() const
367{
368 return fMeanConvFADC2PheStatVar<0 ? -1 : TMath::Sqrt(fMeanConvFADC2PheStatVar);
369}
370
371// --------------------------------------------------------------------------
372//
373Float_t MCalibrationChargePix::GetMeanConvFADC2PheSystErr() const
374{
375 return fMeanConvFADC2PheSystVar<0 ? -1 : TMath::Sqrt(fMeanConvFADC2PheSystVar);
376}
377
378// --------------------------------------------------------------------------
379//
380Float_t MCalibrationChargePix::GetMeanConvFADC2PheTotErr() const
381{
382 if (fMeanConvFADC2PheSystVar<0 || fMeanConvFADC2PheStatVar<0)
383 return -1.;
384
385 return TMath::Sqrt(fMeanConvFADC2PheSystVar+fMeanConvFADC2PheStatVar);
386}
387
388// --------------------------------------------------------------------------
389//
390// Get the error on the mean total F-Factor of the signal readout
391// (F-Factor Method): If fMeanFFactorFADC2PhotVar is smaller than 0
392// (i.e. has not yet been set), return -1. Else returns the square
393// root of fMeanFFactorFADC2PhotVar
394//
395Float_t MCalibrationChargePix::GetMeanFFactorFADC2PhotErr() const
396{
397 return fMeanFFactorFADC2PhotVar<0 ? -1. : TMath::Sqrt(fMeanFFactorFADC2PhotVar);
398}
399
400// --------------------------------------------------------------------------
401//
402// Get the relative variance of the conversion factor between higain and logain:
403// - If fConversionHiLo is 0, return -1.
404// - If fConversionHiLoVar is smaller than 0, return -1.
405// - Else returns: fConversionHiLoVar / fConversionHiLo^2
406//
407const Float_t MCalibrationChargePix::GetConversionHiLoRelVar() const
408{
409
410 if (fConversionHiLoVar < 0.)
411 return -1.;
412
413 if (fConversionHiLo == 0.)
414 return -1.;
415
416 return fConversionHiLoVar / (fConversionHiLo * fConversionHiLo);
417}
418
419// --------------------------------------------------------------------------
420//
421// Get the relative variance of the electronics pedestal RMS
422// - returns rel. variance of gkElectronicPedRms
423//
424const Float_t MCalibrationChargePix::GetElectronicPedRmsRelVar() const
425{
426 return gkElectronicPedRmsErr * gkElectronicPedRmsErr / gkElectronicPedRms / gkElectronicPedRms;
427}
428
429
430// --------------------------------------------------------------------------
431//
432// Get the relative variance of the conversion factor between higain and logain:
433// - If gkFFactor is 0, return -1.
434// - If gkFFactorErr is smaller than 0, return -1.
435// - Else returns: gkFFactorErr^2 / gkFFactor*^2
436//
437const Float_t MCalibrationChargePix::GetFFactorRelVar() const
438{
439
440 if (gkFFactorErr < 0.)
441 return -1.;
442
443 if (gkFFactor == 0.)
444 return -1.;
445
446 return gkFFactorErr * gkFFactorErr / (gkFFactor * gkFFactor);
447}
448
449// --------------------------------------------------------------------------
450//
451Float_t MCalibrationChargePix::GetPedErr() const
452{
453 return TMath::Sqrt(fPedVar);
454}
455
456// --------------------------------------------------------------------------
457//
458Float_t MCalibrationChargePix::GetPedRms() const
459{
460 if (!IsHiGainSaturation())
461 return fPedRms;
462
463 return fLoGainPedRmsSquare<0 ? -1 : TMath::Sqrt(fLoGainPedRmsSquare);
464}
465
466// --------------------------------------------------------------------------
467//
468// Get the Error of the pedestals RMS:
469// - Test bit kHiGainSaturation:
470// If yes, return square root of (0.25*fLoGainPedRmsSquareVar/ fLoGainPedRmsSquare) (if greater than 0, otherwise -1.)
471// If no , return square root of (fPedVar) (if greater than 0, otherwise -1.), divided by 2.
472//
473Float_t MCalibrationChargePix::GetPedRmsErr() const
474{
475 if (IsHiGainSaturation())
476 return fLoGainPedRmsSquareVar < 0. ? -1. :TMath::Sqrt(0.25*fLoGainPedRmsSquareVar/fLoGainPedRmsSquare);
477 else
478 return fPedVar < 0. ? -1. : TMath::Sqrt(fPedVar)/2.;
479}
480
481// --------------------------------------------------------------------------
482//
483// GetPheFFactorMethodErr(): Get the error on the number of photo-electrons
484// (F-Factor Method): If fPheFFactorMethodVar is smaller than 0 (i.e. has
485// not yet been set), return -1. Else returns the square root of
486// fPheFFactorMethodVar
487//
488Float_t MCalibrationChargePix::GetPheFFactorMethodErr() const
489{
490 return fPheFFactorMethodStatVar<0 ? -1 : TMath::Sqrt(fPheFFactorMethodStatVar);
491}
492
493// --------------------------------------------------------------------------
494//
495Float_t MCalibrationChargePix::GetPheFFactorMethodSystErr() const
496{
497 return fPheFFactorMethodSystVar<0 ? -1 : TMath::Sqrt(fPheFFactorMethodSystVar);
498}
499
500// --------------------------------------------------------------------------
501//
502Float_t MCalibrationChargePix::GetPheFFactorMethodTotErr() const
503{
504 if (fPheFFactorMethodStatVar<0 || fPheFFactorMethodSystVar<0)
505 return -1.;
506
507 return TMath::Sqrt(fPheFFactorMethodStatVar+fPheFFactorMethodSystVar);
508}
509
510// --------------------------------------------------------------------------
511//
512// Get the error of the converted reduced Sigma:
513// - If fRSigmaSquareVar is smaller than 0 (i.e. has not yet been set), return -1.
514// - Calculate the absolute variance of the reduced sigma with the formula:
515// reduced sigma variance = 0.25 * fRSigmaSquareVar / fRSigmaSquare
516// - Test bit kHiGainSaturation:
517// If yes, returns the square root of the quadratic sum of the relative variances of the
518// reduced sigma and fConversionHiLo, mulitplied with GetRSigma()
519// Else returns the square root of rel. (0.25*fRSigmaSquareVar / fRSigmaSquare)
520//
521Float_t MCalibrationChargePix::GetConvertedRSigmaErr() const
522{
523
524 if (fRSigmaSquareVar < 0)
525 return -1;
526
527 //
528 // SigmaSquareVar = 4. * Sigma * Sigma * Var(sigma)
529 // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma)
530 //
531 const Float_t rsigmaVar = 0.25 * fRSigmaSquareVar / fRSigmaSquare;
532
533 if (IsHiGainSaturation())
534 return TMath::Sqrt(rsigmaVar/fRSigmaSquare + GetConversionHiLoRelVar()) * GetRSigma();
535 else
536 return TMath::Sqrt(rsigmaVar);
537
538}
539
540// --------------------------------------------------------------------------
541//
542// Get the reduced Sigma:
543// If fRSigmaSquare is smaller than 0 (i.e. has not yet been set),
544// return -1.
545//
546Float_t MCalibrationChargePix::GetRSigma() const
547{
548 return fRSigmaSquare<0 ? -1 : TMath::Sqrt(fRSigmaSquare);
549}
550
551// --------------------------------------------------------------------------
552//
553// Get the error of the reduced Sigma:
554// - If fRSigmaSquareVar is smaller than 0 (i.e. has not yet been set), return -1.
555// - Calculate the absolute variance of the reduced sigma with the formula:
556// reduced sigma variance = 0.25 * fRSigmaSquareVar / fRSigmaSquare
557//
558Float_t MCalibrationChargePix::GetRSigmaErr() const
559{
560
561 if (fRSigmaSquareVar < 0)
562 return -1;
563
564 //
565 // SigmaSquareVar = 4. * Sigma * Sigma * Var(sigma)
566 // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma)
567 //
568 return TMath::Sqrt(0.25 * fRSigmaSquareVar / fRSigmaSquare);
569
570}
571
572// --------------------------------------------------------------------------
573//
574// Get the reduced Sigma per Charge:
575// - If GetRSigma() is smaller or equal 0. (i.e. has not yet been set), return -1.
576// - If GetMean() is 0. or -1. (i.e. has not yet been set), return -1.
577// - Return GetRSigma() / GetMean()
578//
579Float_t MCalibrationChargePix::GetRSigmaPerCharge() const
580{
581
582 const Float_t rsigma = GetRSigma();
583
584 if (rsigma <= 0)
585 return -1.;
586
587
588 const Float_t mean = GetMean();
589
590 if (mean == 0. || mean == -1.)
591 return -1.;
592
593 return rsigma / mean;
594}
595
596
597// --------------------------------------------------------------------------
598//
599// Get the error of the reduced Sigma per Charge:
600// - If GetRSigmaRelVar() is smaller or equal 0. (i.e. has not yet been set), return -1.
601// - If GetMeanRelVar() is smaller or equal 0. (i.e. has not yet been set), return -1.
602// - Return the propagated error of GetRSigmaPerCharge()
603//
604Float_t MCalibrationChargePix::GetRSigmaPerChargeErr() const
605{
606
607 const Float_t rsigmarelvar = GetRSigmaRelVar();
608
609 if (rsigmarelvar <= 0)
610 return -1.;
611
612
613 const Float_t meanrelvar = GetMeanRelVar();
614
615 if (meanrelvar <= 0.)
616 return -1.;
617
618 return TMath::Sqrt(rsigmarelvar + meanrelvar) * GetRSigmaPerCharge();
619}
620
621// --------------------------------------------------------------------------
622//
623// Get the relative variance of the reduced Sigma:
624// - If fRSigmaSquareVar is smaller than 0 (i.e. has not yet been set), return -1.
625// - Calculate the relative variance of the reduced sigma squares with the formula:
626// reduced sigma rel. variance = 0.25 * fRSigmaSquareVar / fRSigmaSquare / fRSigmaSquare
627// - Test bit kHiGainSaturation:
628// If yes, returns the sum of the relative variances of the reduced sigma and fConversionHiLo
629// Else returns the relative variance of the reduced sigma
630//
631Float_t MCalibrationChargePix::GetRSigmaRelVar() const
632{
633
634 if (fRSigmaSquareVar < 0)
635 return -1;
636
637 //
638 // SigmaSquareVar = 4. * Sigma * Sigma * Var(sigma)
639 // ==> Var(sigma) = 0.25 * SigmaSquareVar / (Sigma * Sigma)
640 //
641 return 0.25 * fRSigmaSquareVar / ( fRSigmaSquare * fRSigmaSquare );
642
643}
644
645
646
647// --------------------------------------------------------------------------
648//
649// Test bit kFFactorMethodValid
650//
651Bool_t MCalibrationChargePix::IsFFactorMethodValid() const
652{
653 return TESTBIT(fCalibFlags, kFFactorMethodValid);
654}
655
656
657// ----------------------------------------------------------------------------
658//
659// - If fSigma is smaller than 0 (i.e. has not yet been set), return kFALSE
660// - If fPedRms is smaller than 0 (i.e. has not yet been set), return kFALSE
661//
662// Calculate the reduced sigma of the low-Gain FADC slices:
663// - Test bit IsHiGainSaturation() for the Sigma:
664// If yes, take fLoGainSigma and fLoGainSigmaVar
665// If no , take fHiGainSigma and fHiGainSigmaVar
666//
667// - Test bit IsHiGainSaturation() for the pedRMS:
668// If yes, take fLoGainPedRmsSquare and fLoGainPedRmsSquareVar
669// If no , take fPedRms and fPedVar
670//
671// - Calculate the reduced sigma with the formula:
672// fRSigmaSquare = Sigma*Sigma - pedRMS*pedRMS
673//
674// - If fRSigmaSquare is smaller than 0, give a warning and return kFALSE
675//
676// - Calculate the variance of the reduced sigma with the formula:
677// fRSigmaSquareVar = 4.* (sigmaVar*Sigma*Sigma + pedRmsVar*pedRMS*pedRMS)
678//
679// A back-transformation to the corr. amplification factor of the High-Gain is done
680// in GetRSigma() and GetRSigmaErr()
681//
682Bool_t MCalibrationChargePix::CalcReducedSigma(const Float_t extractorres)
683{
684
685 if (GetSigma() < 0.)
686 return kFALSE;
687
688 if (GetPedRms() < 0.)
689 return kFALSE;
690
691 const Float_t sigma = IsHiGainSaturation() ? fLoGainSigma : fHiGainSigma ;
692 const Float_t sigmavar = IsHiGainSaturation() ? fLoGainSigmaVar : fHiGainSigmaVar;
693 const Float_t pedRmsSquare = IsHiGainSaturation() ? fLoGainPedRmsSquare : fPedRms*fPedRms;
694 const Float_t pedRmsSquareVar = IsHiGainSaturation() ? fLoGainPedRmsSquareVar : 0.25*fPedVar*pedRmsSquare;
695 const Float_t extractorresVar = extractorres * extractorres;
696
697 if (IsDebug())
698 {
699 *fLog << dbginf << "ID: " << GetPixId()
700 << " HiGainSaturation: " << IsHiGainSaturation()
701 << " Sigma: " << sigma
702 << " Var.Sigma: " << sigmavar
703 << " PedRmsSquare: " << pedRmsSquare
704 << " pedRmsSquareVar: " << pedRmsSquareVar
705 << " extractorresVar: " << extractorresVar
706 << endl;
707 }
708
709 const Float_t sigmaSquare = sigma * sigma;
710 const Float_t sigmaSquareVar = 4. * sigmavar * sigmaSquare;
711
712 //
713 // Calculate the reduced sigmas
714 //
715 fRSigmaSquare = sigmaSquare - pedRmsSquare - extractorresVar;
716
717 if (IsDebug())
718 {
719 *fLog << dbginf << "ID: " << GetPixId()
720 << " Red.Sigma Square: " << fRSigmaSquare
721 << endl;
722 }
723
724 if (fRSigmaSquare <= 0.)
725 {
726 if (IsDebug())
727 *fLog << warn
728 << "WARNING: Cannot calculate the reduced sigma: smaller than 0 in pixel "
729 << fPixId << endl;
730 return kFALSE;
731 }
732
733
734 fRSigmaSquareVar = 4. * (sigmaSquareVar + pedRmsSquareVar);
735
736 if (IsDebug())
737 {
738 *fLog << dbginf << "ID: " << GetPixId()
739 << " Var.Red.Sigma Square: " << fRSigmaSquareVar
740 << endl;
741 }
742
743 return kTRUE;
744}
745
746// ------------------------------------------------------------------
747//
748// If fRSigmaSquare is smaller than 0 (i.e. has not yet been set),
749// return kFALSE
750//
751// Calculate the number of photo-electrons with the F-Factor method:
752// - Test bit IsHiGainSaturation() for the Mean Sum of FADC slices:
753// If yes, take fLoGainMean and fLoGainMeanVar
754// If no , take fHiGainMean and fHiGainMeanVar
755//
756// - Test bit IsHiGainSaturation() for the pedRMS:
757// If yes, take fLoGainPedRmsSquare and fLoGainPedRmsSquareVar
758// If no , take fPedRms and fPedVar
759//
760// - Calculate the number of photo-electrons with the formula:
761// fPheFFactorMethod = gkFFactor*gkFFactor * Mean * Mean / fRSigmaSquare
762//
763// - Calculate the Variance on the photo-electrons with the formula:
764// fPheFFactorMethodVar = ( 4. * gkFFactorErr * gkFFactorErr / ( gkFFactor * gkFFactor )
765// + 4. * Mean Var. / ( Mean * Mean )
766// + fRSigmaSquareVar / fRSigmaSquare
767// ) * fPheFFactorMethod * fPheFFactorMethod
768//
769// - If fPheFFactorMethod is less than fPheFFactorMethodLimit,
770// set kFFactorMethodValid to kFALSE and return kFALSE
771//
772Bool_t MCalibrationChargePix::CalcFFactor()
773{
774 SetFFactorMethodValid(kFALSE);
775
776 if (fRSigmaSquare < 0.)
777 return kFALSE;
778
779 //
780 // Square all variables in order to avoid applications of square root
781 //
782 const Float_t meanSquare = GetMean() * GetMean();
783 const Float_t meanSquareRelVar = 4.* GetMeanRelVar();
784
785 const Float_t ffactorsquare = gkFFactor * gkFFactor;
786 const Float_t ffactorsquareRelVar = 4.* GetFFactorRelVar();
787
788 const Float_t rsigmaSquareRelVar = fRSigmaSquareVar / fRSigmaSquare / fRSigmaSquare;
789 //
790 // Calculate the number of phe's from the F-Factor method
791 // (independent on Hi Gain or Lo Gain)
792 //
793 fPheFFactorMethod = ffactorsquare * meanSquare / fRSigmaSquare;
794
795 if (IsDebug())
796 {
797 *fLog << dbginf << "ID: " << GetPixId()
798 << " F-Factor Square: " << ffactorsquare
799 << " Mean Square: " << meanSquare
800 << " Red.Sigma Square: " << fRSigmaSquare
801 << " Photo-electrons: " << fPheFFactorMethod
802 << endl;
803 }
804
805 if (fPheFFactorMethod < fPheFFactorMethodLimit)
806 return kFALSE;
807
808 //
809 // Calculate the Error of Nphe
810 //
811 const Float_t pheRelVar = meanSquareRelVar + rsigmaSquareRelVar;
812 fPheFFactorMethodStatVar = pheRelVar * fPheFFactorMethod * fPheFFactorMethod;
813 fPheFFactorMethodSystVar = ffactorsquareRelVar * fPheFFactorMethod * fPheFFactorMethod;
814
815 if (IsDebug())
816 {
817 *fLog << dbginf << "ID: " << GetPixId()
818 << " Rel.Var.F-Factor Square: " << ffactorsquareRelVar
819 << " Rel.Var. Mean Square: " << meanSquareRelVar
820 << " Rel.Var. Red.Sigma Square: " << rsigmaSquareRelVar
821 << " Rel.Var. Photo-electrons: " << pheRelVar
822 << endl;
823 }
824
825 if (fPheFFactorMethodStatVar < 0. )
826 return kFALSE;
827
828 return kTRUE;
829}
830
831// ------------------------------------------------------------------
832//
833// If fPheFFactorMethod is smaller than 0 (i.e. has not yet been set),
834// return kFALSE
835//
836// If GetCovertedMean() is smaller than 0 (i.e. has not yet been set),
837// return kFALSE
838//
839// Calculate fMeanConvFADC2Phe with the following formula:
840//
841// fMeanConvFADC2Phe = fPheFFactorMethod / GetConvMean();
842//
843// Calculate the rel. variance of fMeanConvFADC2Phe, taking into account that
844// in the calculation of the number of phe's one mean square has already been used.
845// Now, dividing by another mean, one mean calcels out, one cannot directly propagate
846// the errors, but instead havs to take into account this cancellation:
847//
848// convrelvar = ffactorsquareRelVar + GetMeanRelVar() + rsigmaSquareRelVar;
849//
850// If confrelvar is smaller than 0. or greater than fConvFFactorRelVarLimit,
851// return kFALSE
852//
853// Calculate the variance of fMeanConvFADC2Phe with the formula:
854//
855// fMeanConvFADC2PheVar = convrelvar * fMeanConvFADC2Phe * fMeanConvFADC2Phe;
856//
857// Set kFFactorMethodValid to kTRUE and
858// return kTRUE
859//
860Bool_t MCalibrationChargePix::CalcConvFFactor()
861{
862
863 if (fPheFFactorMethod <= 0.)
864 return kFALSE;
865
866 const Float_t convmean = GetConvertedMean();
867
868 if (convmean <= 0.)
869 return kFALSE;
870
871 fMeanConvFADC2Phe = fPheFFactorMethod / convmean;
872
873 if (IsDebug())
874 {
875 *fLog << dbginf << "ID: " << GetPixId()
876 << " Converted Mean: " << convmean
877 << " Conversion FADC2Phe: " << fMeanConvFADC2Phe
878 << endl;
879 }
880
881 const Float_t ffactorsquareRelVar = 4.* GetFFactorRelVar();
882 const Float_t rsigmaSquareRelVar = fRSigmaSquareVar / fRSigmaSquare / fRSigmaSquare;
883 //
884 // In the calculation of the number of phe's one mean square has already been used.
885 // Now, we divide by another mean, so one mean calcels out, we cannot directly propagate
886 // the errors, but have to take account of this cancellation:
887 //
888 Float_t convrelvar = GetMeanRelVar() + rsigmaSquareRelVar;
889 if (IsHiGainSaturation())
890 convrelvar += GetConversionHiLoRelVar();
891
892 const Float_t limit = IsHiGainSaturation() ? fConvFFactorRelVarLimit * 4. : fConvFFactorRelVarLimit;
893
894 //
895 // Also have to take into account the pixels labelled MBadPixelsPix::kChargeSigmaNotValid which do not
896 // have a fRSigmaSquareVar, calculate their error directly!
897 //
898 if (fRSigmaSquareVar < 0.)
899 convrelvar = GetMeanRelVar() + GetPheFFactorMethodRelVar();
900
901 if (IsDebug())
902 {
903 *fLog << dbginf << "ID: " << GetPixId()
904 << " Rel.Var.Red.Sigma: " << rsigmaSquareRelVar
905 << " Rel.Var.Mean: " << GetMeanRelVar()
906 << " Rel.Var.F-Factor: " << ffactorsquareRelVar
907 << " Rel.Var.Conversion FADC2Phe: " << convrelvar
908 << endl;
909 }
910
911 if (convrelvar > limit || convrelvar < 0.)
912 {
913 *fLog << warn << "pixel " << setw(4) << fPixId << ": Conv. F-Factor Method Rel. Var.: "
914 << Form("%4.3f out of limits: [0,%3.2f]",convrelvar,limit) << endl;
915 return kFALSE;
916 }
917
918 fMeanConvFADC2PheStatVar = convrelvar * fMeanConvFADC2Phe * fMeanConvFADC2Phe;
919 fMeanConvFADC2PheSystVar = ffactorsquareRelVar * fMeanConvFADC2Phe * fMeanConvFADC2Phe;
920
921 SetFFactorMethodValid(kTRUE);
922 return kTRUE;
923}
924
925// ----------------------------------------------------------------------------------
926//
927// If photflux is smaller or equal 0, return kFALSE
928//
929// Calculate the total F-Factor with the formula:
930// fMeanFFactorFADC2Phot = Sqrt ( fRSigmaSquare ) / GetMean() * sqrt(nphotons)
931//
932// Calculate the error of the total F-Factor
933//
934Bool_t MCalibrationChargePix::CalcMeanFFactor( const Float_t nphotons, const Float_t nphotonsrelvar )
935{
936
937
938 if (IsDebug())
939 {
940 *fLog << dbginf << "ID: " << GetPixId()
941 << " Number photons: " << nphotons
942 << " Rel.Var.Number photons: " << nphotonsrelvar
943 << " Red.Sigma Square: " << fRSigmaSquare
944 << " Mean: " << GetMean()
945 << endl;
946 }
947
948
949 if (nphotons <= 0.)
950 {
951 *fLog << warn << GetDescriptor() << ": Assumed photon flux is smaller or equal 0." << endl;
952 return kFALSE;
953 }
954
955 if (nphotonsrelvar < 0.)
956 {
957 *fLog << warn << GetDescriptor() << ": Assumed photon flux variance is smaller than 0." << endl;
958 return kFALSE;
959 }
960
961 fMeanFFactorFADC2Phot = TMath::Sqrt(fRSigmaSquare * nphotons) / GetMean() ;
962
963 if (IsDebug())
964 {
965 *fLog << dbginf << "ID: " << GetPixId()
966 << " F-Factor FADC2Phot: " << fMeanFFactorFADC2Phot
967 << endl;
968 }
969
970 if (fMeanFFactorFADC2Phot < 0.)
971 {
972 *fLog << warn << GetDescriptor() << ": F-Factor photons to FADC counts smaller than 0." << endl;
973 return kFALSE;
974 }
975
976 const Float_t ffactorrelvar = 0.25 * fRSigmaSquareVar / ( fRSigmaSquare * fRSigmaSquare)
977 + GetMeanRelVar()
978 + 0.25 * nphotonsrelvar;
979
980 fMeanFFactorFADC2PhotVar = ffactorrelvar * fMeanFFactorFADC2Phot * fMeanFFactorFADC2Phot;
981
982 if (IsDebug())
983 {
984 *fLog << dbginf << "ID: " << GetPixId()
985 << " Rel.Var.Red.Sigma: " << 0.25 * fRSigmaSquareVar / ( fRSigmaSquare * fRSigmaSquare)
986 << " Rel.Var.Mean: " << GetMeanRelVar()
987 << " Rel.Var.photons: " << 0.25 * nphotonsrelvar
988 << " Rel.Var.F-Factor FADC2Phot: " << ffactorrelvar
989 << endl;
990 }
991
992 return kTRUE;
993}
994
995
996// ----------------------------------------------------------------------------
997//
998// - If fPed is smaller than 0 (i.e. has not yet been set), return.
999// - If fPedVar is smaller than 0 (i.e. has not yet been set), return.
1000//
1001// Calculate the electronic pedestal RMS with the formula:
1002// - elec. pedestal = gkElectronicPedRms * sqrt(logainsamples)
1003//
1004// Calculate the night sky background ped. RMS contribution ("NSB") in the high-gain
1005// from the high gain Pedestal RMS with the formula:
1006// - HiGain NSB square = fPedRms * fPedRms - elec.ped.* elec.ped.
1007// - Var(HiGain NSB square) = fPedVar * fPedRms * fPedRms + 4.*elecPedRmsVar * elec.ped.* elec.ped.
1008//
1009// If HiGain NSB square is smaller than 0., set it to zero. (but not the error!)
1010//
1011// Convert the NSB ped. RMS contribution to the low-gain with the formula:
1012// - LoGain NSB square = HiGain NSB square / (fConversionHiLo*fConversionHiLo)
1013// - Var(LoGain NSB square) = ( Var(HiGain NSB square) / (HiGain NSB square * HiGain NSB square)
1014// + GetConversionHiLoRelVar()
1015// ) * LoGain NSB square * LoGain NSB square
1016//
1017// - Low Gain Ped RMS Square = LoGain NSB square + elec.ped. square
1018// Var (Low Gain Ped RMS Square) = Var(LoGain NSB square) + Var(elec.ped. square)
1019//
1020void MCalibrationChargePix::CalcLoGainPedestal(Float_t logainsamples)
1021{
1022
1023 if (fPedRms < 0.)
1024 return;
1025
1026 if (fPedVar < 0.)
1027 return;
1028
1029 const Float_t elecPedRms = gkElectronicPedRms * TMath::Sqrt(logainsamples);
1030 const Float_t elecPedRmsVar = GetElectronicPedRmsRelVar() * elecPedRms * elecPedRms;
1031
1032 Float_t pedRmsSquare = fPedRms * fPedRms;
1033 Float_t pedRmsSquareVar = fPedVar * pedRmsSquare; // fPedRmsErr = fPedErr/2.
1034
1035 //
1036 // We do not know the Low Gain Pedestal RMS, so we have to retrieve it
1037 // from the High Gain:
1038 //
1039 // We extract the pure NSB contribution:
1040 //
1041 const Float_t elecRmsSquare = elecPedRms * elecPedRms;
1042 const Float_t elecRmsSquareVar = 4.*elecPedRmsVar * elecRmsSquare;
1043
1044 if (IsDebug())
1045 {
1046 *fLog << dbginf << "ID: " << GetPixId()
1047 << " Ped.Rms Square: " << pedRmsSquare
1048 << " Elec.Rms Square: " << elecRmsSquare
1049 << " Ped.Rms.Square Var.: " << pedRmsSquareVar
1050 << " Elec.Rms Square Var.: " << elecRmsSquareVar
1051 << endl;
1052 }
1053
1054
1055 Float_t higainNsbSquare = pedRmsSquare - elecRmsSquare;
1056 Float_t higainNsbSquareVar = (pedRmsSquareVar + elecRmsSquareVar);
1057
1058 if (higainNsbSquare < 0.001)
1059 higainNsbSquare = 0.;
1060
1061 if (IsDebug())
1062 {
1063 *fLog << dbginf << "ID: " << GetPixId()
1064 << " HiGain NSB Square: " << higainNsbSquare
1065 << " Var.HiGain NSB Square: " << higainNsbSquareVar
1066 << endl;
1067 }
1068
1069 //
1070 // Now, we divide the NSB by the conversion factor and
1071 // add it quadratically to the electronic noise
1072 //
1073 const Float_t conversionSquare = fConversionHiLo * fConversionHiLo;
1074 const Float_t conversionSquareRelVar = 4.* GetConversionHiLoRelVar();
1075
1076 const Float_t logainNsbSquare = higainNsbSquare / conversionSquare;
1077 //
1078 // Calculation of variance of: c = a/b
1079 // Delta(c)^2 = ( Delta(a)^2 + a^2/b^2*(Delta(b)^2 ) / b^2
1080 //
1081 const Float_t logainNsbSquareVar = ( higainNsbSquareVar
1082 + conversionSquareRelVar * higainNsbSquare * higainNsbSquare )
1083 / conversionSquare / conversionSquare;
1084
1085 fLoGainPedRmsSquare = logainNsbSquare + elecRmsSquare;
1086 fLoGainPedRmsSquareVar = logainNsbSquareVar + elecRmsSquareVar;
1087
1088 if (IsDebug())
1089 {
1090 *fLog << dbginf << "ID: " << GetPixId()
1091 << " LoGain Ped Rms Square: " << fLoGainPedRmsSquare
1092 << " Var.Ped Rms Square: " << fLoGainPedRmsSquareVar
1093 << endl;
1094 }
1095
1096
1097}
1098
Note: See TracBrowser for help on using the repository browser.