source: branches/Mars_McMismatchStudy/mcalib/MCalibrationBlindPix.cc@ 18558

Last change on this file since 18558 was 10166, checked in by tbretz, 14 years ago
Removed the old obsolete cvs header line.
File size: 16.2 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// MCalibrationBlindPix
28//
29// Storage container of the fit results of the Blind Pixel signal
30// (from MHCalibrationChargeBlindPix).
31//
32// The Flux is calculated in photons per mm^2 in the camera plane.
33//
34// Currently, the following numbers are implemented:
35// - fArea: 100 mm^2
36// - Average QE of Blind Pixel:
37// fQEGreen: 0.154
38// fQEBlue : 0.226
39// fQEUV : 0.247
40// fQECT1 : 0.247
41// - Average QE Error of Blind Pixel:
42// fQEGreenErr: 0.015;
43// fQEBlueErr : 0.02;
44// fQEUVErr : 0.02;
45// fQECT1Err : 0.02;
46// - Attenuation factor Blind Pixel:
47// fAttGreen : 1.97;
48// fAttBlue : 1.96;
49// fAttUV : 1.95;
50// fAttCT1 : 1.95;
51//
52//
53/////////////////////////////////////////////////////////////////////////////
54#include "MCalibrationBlindPix.h"
55
56#include <TMath.h>
57#include <TH1.h>
58
59#include "MLog.h"
60#include "MLogManip.h"
61
62ClassImp(MCalibrationBlindPix);
63
64using namespace std;
65
66const Float_t MCalibrationBlindPix::fgArea = 100;
67const Float_t MCalibrationBlindPix::fgAttGreen = 1.97;
68const Float_t MCalibrationBlindPix::fgAttBlue = 1.96;
69const Float_t MCalibrationBlindPix::fgAttUV = 1.95;
70const Float_t MCalibrationBlindPix::fgAttCT1 = 1.95;
71const Float_t MCalibrationBlindPix::fgAttErr = 0.01;
72const Float_t MCalibrationBlindPix::fgQEGreen = 0.154;
73const Float_t MCalibrationBlindPix::fgQEBlue = 0.226;
74const Float_t MCalibrationBlindPix::fgQEUV = 0.247;
75const Float_t MCalibrationBlindPix::fgQECT1 = 0.247;
76const Float_t MCalibrationBlindPix::fgQEErrGreen = 0.005;
77const Float_t MCalibrationBlindPix::fgQEErrBlue = 0.007;
78const Float_t MCalibrationBlindPix::fgQEErrUV = 0.01;
79const Float_t MCalibrationBlindPix::fgQEErrCT1 = 0.01;
80const Float_t MCalibrationBlindPix::fgCollEffGreen = 0.99;
81const Float_t MCalibrationBlindPix::fgCollEffBlue = 0.93;
82const Float_t MCalibrationBlindPix::fgCollEffUV = 0.90;
83const Float_t MCalibrationBlindPix::fgCollEffCT1 = 0.90;
84const Float_t MCalibrationBlindPix::fgCollEffErr = 0.05;
85
86// --------------------------------------------------------------------------
87//
88// Default Constructor.
89//
90// Calls:
91// - Clear()
92//
93// For backward-compatibility reasons, quantum eff., coll. eff. and att.
94// are intialized from the static members. This should, however, be
95// overwritten by a class deriving from MCalibrationBlindCam.
96//
97MCalibrationBlindPix::MCalibrationBlindPix(const char *name, const char *title)
98{
99
100 fName = name ? name : "MCalibrationBlindPix";
101 fTitle = title ? title : "Container of the fit results of the blind pixel";
102
103 Clear();
104
105 fArea = fgArea;
106 fAreaErr = 0.;
107
108 Float_t att[MCalibrationCam::gkNumPulserColors];
109
110 att [ MCalibrationCam::kGREEN ] = fgAttGreen;
111 att [ MCalibrationCam::kBLUE ] = fgAttBlue;
112 att [ MCalibrationCam::kUV ] = fgAttUV;
113 att [ MCalibrationCam::kCT1 ] = fgAttCT1;
114
115 SetAtt(MCalibrationCam::gkNumPulserColors,att);
116
117 Float_t atterr[MCalibrationCam::gkNumPulserColors];
118
119 atterr [ MCalibrationCam::kGREEN ] = fgAttErr;
120 atterr [ MCalibrationCam::kBLUE ] = fgAttErr;
121 atterr [ MCalibrationCam::kUV ] = fgAttErr;
122 atterr [ MCalibrationCam::kCT1 ] = fgAttErr;
123
124 SetAttErr(MCalibrationCam::gkNumPulserColors,atterr);
125
126 Float_t qe[MCalibrationCam::gkNumPulserColors];
127
128 qe [ MCalibrationCam::kGREEN ] = fgQEGreen;
129 qe [ MCalibrationCam::kBLUE ] = fgQEBlue;
130 qe [ MCalibrationCam::kUV ] = fgQEUV;
131 qe [ MCalibrationCam::kCT1 ] = fgQECT1;
132
133 SetQE(MCalibrationCam::gkNumPulserColors,qe);
134
135 Float_t qeerr[MCalibrationCam::gkNumPulserColors];
136
137 qeerr [ MCalibrationCam::kGREEN ] = fgQEErrGreen;
138 qeerr [ MCalibrationCam::kBLUE ] = fgQEErrBlue;
139 qeerr [ MCalibrationCam::kUV ] = fgQEErrUV;
140 qeerr [ MCalibrationCam::kCT1 ] = fgQEErrCT1;
141
142 SetQEErr(MCalibrationCam::gkNumPulserColors,qeerr);
143
144 Float_t colleff[MCalibrationCam::gkNumPulserColors];
145
146 colleff [ MCalibrationCam::kGREEN ] = fgCollEffGreen;
147 colleff [ MCalibrationCam::kBLUE ] = fgCollEffBlue;
148 colleff [ MCalibrationCam::kUV ] = fgCollEffUV;
149 colleff [ MCalibrationCam::kCT1 ] = fgCollEffCT1;
150
151 SetCollEff(MCalibrationCam::gkNumPulserColors,colleff);
152
153 Float_t collefferr[MCalibrationCam::gkNumPulserColors];
154
155 collefferr [ MCalibrationCam::kGREEN ] = fgCollEffErr;
156 collefferr [ MCalibrationCam::kBLUE ] = fgCollEffErr;
157 collefferr [ MCalibrationCam::kUV ] = fgCollEffErr;
158 collefferr [ MCalibrationCam::kCT1 ] = fgCollEffErr;
159
160 SetCollEffErr(MCalibrationCam::gkNumPulserColors,collefferr);
161}
162
163// ------------------------------------------------------------------------
164//
165// Sets:
166// - all flags to kFALSE
167// - all variables to -1.
168// - the fColor to MCalibrationCam::kNONE
169//
170// Calls:
171// - MCalibrationPix::Clear()
172//
173void MCalibrationBlindPix::Clear(Option_t *o)
174{
175
176 fFluxInsidePlexiglass = -1.;
177 fFluxInsidePlexiglassVar = -1.;
178 fLambda = -1.;
179 fLambdaCheck = -1.;
180 fLambdaVar = -1.;
181 fMu0 = -1.;
182 fMu0Err = -1.;
183 fMu1 = -1.;
184 fMu1Err = -1.;
185 fSigma0 = -1.;
186 fSigma0Err = -1.;
187 fSigma1 = -1.;
188 fSigma1Err = -1.;
189
190 SetOscillating ( kFALSE );
191 SetExcluded ( kFALSE );
192 SetChargeFitValid ( kFALSE );
193 SetPedestalFitOK ( kFALSE );
194 SetSinglePheFitOK ( kFALSE );
195 SetFluxInsidePlexiglassAvailable ( kFALSE );
196
197 SetColor(MCalibrationCam::kNONE);
198
199 MCalibrationPix::Clear();
200}
201
202void MCalibrationBlindPix::SetFluxInsidePlexiglassAvailable( const Bool_t b)
203{
204 b ? SETBIT(fFlags,kFluxInsidePlexiglassAvailable) : CLRBIT(fFlags,kFluxInsidePlexiglassAvailable);
205}
206
207
208// --------------------------------------------------------------------------
209//
210// Set the Oscillating Bit from outside
211//
212void MCalibrationBlindPix::SetOscillating( const Bool_t b)
213{
214 b ? SETBIT(fFlags,kOscillating) : CLRBIT(fFlags,kOscillating);
215}
216
217// -----------------------------------------------------
218//
219// copy 'constructor'
220//
221void MCalibrationBlindPix::Copy(TObject& object) const
222{
223
224 MCalibrationBlindPix &pix = (MCalibrationBlindPix&)object;
225
226 //
227 // Copy the data members
228 //
229 pix.fPixId = fPixId;
230 pix.fFlags = fFlags;
231 pix.fArea = fArea;
232 pix.fAreaErr = fAreaErr;
233 pix.fColor = fColor;
234
235 //
236 // Copy arrays
237 //
238 pix.fAtt = fAtt;
239 pix.fAttErr = fAttErr;
240 pix.fQE = fQE;
241 pix.fQEErr = fQEErr;
242 pix.fCollEff = fCollEff;
243 pix.fCollEffErr = fCollEffErr;
244
245 pix.fLambda = fLambda;
246 pix.fLambdaCheck = fLambdaCheck;
247 pix.fLambdaCheckErr = fLambdaCheckErr;
248 pix.fLambdaVar = fLambdaVar;
249 pix.fFluxInsidePlexiglass = fFluxInsidePlexiglass;
250 pix.fFluxInsidePlexiglassVar = fFluxInsidePlexiglassVar;
251 pix.fMu0 = fMu0;
252 pix.fMu0Err = fMu0Err;
253 pix.fMu1 = fMu1;
254 pix.fMu1Err = fMu1Err;
255 pix.fSigma0 = fSigma0;
256 pix.fSigma0Err = fSigma0Err;
257 pix.fSigma1 = fSigma1;
258 pix.fSigma1Err = fSigma1Err;
259
260}
261
262
263// --------------------------------------------------------------------------
264//
265// Set the ChargeFitValid Bit from outside
266//
267void MCalibrationBlindPix::SetChargeFitValid( const Bool_t b)
268{
269 b ? SETBIT(fFlags,kChargeFitValid) : CLRBIT(fFlags,kChargeFitValid);
270}
271
272// --------------------------------------------------------------------------
273//
274// Set the PedestalFitValid Bit from outside
275//
276void MCalibrationBlindPix::SetPedestalFitOK( const Bool_t b)
277{
278 b ? SETBIT(fFlags,kPedestalFitOK) : CLRBIT(fFlags,kPedestalFitOK);
279}
280
281// --------------------------------------------------------------------------
282//
283// Set the SinglePheFitValid Bit from outside
284//
285void MCalibrationBlindPix::SetSinglePheFitOK( const Bool_t b)
286{
287 b ? SETBIT(fFlags,kSinglePheFitOK) : CLRBIT(fFlags,kSinglePheFitOK);
288}
289
290// --------------------------------------------------------------------------
291//
292// Return -1 if fFluxInsidePlexiglassVar is smaller than 0.
293// Return square root of fFluxInsidePlexiglassVar
294//
295const Float_t MCalibrationBlindPix::GetFluxInsidePlexiglassErr() const
296{
297 if (fFluxInsidePlexiglassVar < 0.)
298 return -1.;
299
300 return TMath::Sqrt(fFluxInsidePlexiglassVar);
301}
302
303// --------------------------------------------------------------------------
304//
305// Return -1 if fFluxInsidePlexiglassVar is smaller than 0.
306// Return -1 if fFluxInsidePlexiglass is 0.
307// Return fFluxInsidePlexiglassVar / fFluxInsidePlexiglass^2
308//
309const Float_t MCalibrationBlindPix::GetFluxInsidePlexiglassRelVar() const
310{
311 if (fFluxInsidePlexiglassVar < 0.)
312 return -1.;
313
314 if (fFluxInsidePlexiglass == 0.)
315 return -1.;
316
317 return fFluxInsidePlexiglassVar / (fFluxInsidePlexiglass * fFluxInsidePlexiglass) ;
318}
319
320// --------------------------------------------------------------------------
321//
322// Return -1 if fLambdaVar is smaller than 0.
323// Return square root of fLambdaVar
324//
325const Float_t MCalibrationBlindPix::GetLambdaErr() const
326{
327 if (fLambdaVar < 0.)
328 return -1.;
329
330 return TMath::Sqrt(fLambdaVar);
331}
332
333// --------------------------------------------------------------------------
334//
335// Return -1 if fLambdaVar is smaller than 0.
336// Return -1 if fLambda is 0.
337// Return fLambdaVar / (fLambda * fLambda )
338//
339const Float_t MCalibrationBlindPix::GetLambdaRelVar() const
340{
341 if (fLambdaVar < 0.)
342 return -1.;
343
344 if (fLambda == 0.)
345 return -1.;
346
347 return fLambdaVar / fLambda / fLambda ;
348}
349
350// --------------------------------------------------------------------------
351//
352// Return fAreaErr^2 / (fArea^2 )
353//
354const Float_t MCalibrationBlindPix::GetAreaRelVar() const
355{
356 return fAreaErr*fAreaErr/fArea/fArea;
357}
358
359
360// --------------------------------------------------------------------------
361//
362// Return TMath::Power(10,fAtt[fColor])
363//
364const Float_t MCalibrationBlindPix::GetAtt() const
365{
366 return TMath::Power(10,fAtt[fColor]);
367}
368
369// --------------------------------------------------------------------------
370//
371// Return -1 if fAttErr[fColor] is smaller than 0.
372// Error of TMath::Power(10,fAtt[fColor]) = TMath::Power(10,fAtt[fColor])*ln(10.)*fAttErr[fColor]
373// Return fAttErr^2 / (fAtt^2 )
374//
375const Float_t MCalibrationBlindPix::GetAttRelVar() const
376{
377
378 const Float_t ferr = fAttErr[fColor];
379
380 if (ferr < 0.)
381 return -1.;
382
383 return ferr*ferr*2.3;
384}
385
386// --------------------------------------------------------------------------
387//
388// Return fQE[fColor]
389//
390const Float_t MCalibrationBlindPix::GetQE() const
391{
392 return fQE[fColor];
393}
394
395// --------------------------------------------------------------------------
396//
397// Return -1 if fQEErr[fColor] is smaller than 0.
398// Return fQEErr^2 / (fQE^2 )
399//
400const Float_t MCalibrationBlindPix::GetQERelVar() const
401{
402
403 if (fQEErr[fColor] < 0.)
404 return -1.;
405
406 return fQEErr[fColor]* fQEErr[fColor] / GetQE() / GetQE();
407}
408
409// --------------------------------------------------------------------------
410//
411// Return fCollEff[fColor]
412//
413const Float_t MCalibrationBlindPix::GetCollEff() const
414{
415 return fCollEff[fColor];
416}
417
418// --------------------------------------------------------------------------
419//
420// Return -1 if fCollEffErr[fColor] is smaller than 0.
421// Return fCollEffErr^2 / (fCollEff^2 )
422//
423const Float_t MCalibrationBlindPix::GetCollEffRelVar() const
424{
425
426 if (fCollEffErr[fColor] < 0.)
427 return -1.;
428
429 return fCollEffErr[fColor]* fCollEffErr[fColor] / GetCollEff() / GetCollEff();
430}
431
432// --------------------------------------------------------------------------
433//
434// Test bit kChargeFitValid
435//
436const Bool_t MCalibrationBlindPix::IsChargeFitValid() const
437{
438 return TESTBIT(fFlags,kChargeFitValid);
439}
440
441// --------------------------------------------------------------------------
442//
443// Test bit kOscillating
444//
445const Bool_t MCalibrationBlindPix::IsOscillating() const
446{
447 return TESTBIT(fFlags,kOscillating);
448}
449
450// --------------------------------------------------------------------------
451//
452// Test bit kPedestalFitValid
453//
454const Bool_t MCalibrationBlindPix::IsPedestalFitOK() const
455{
456 return TESTBIT(fFlags,kPedestalFitOK);
457}
458
459// --------------------------------------------------------------------------
460//
461// Test bit kSinglePheFitValid
462//
463const Bool_t MCalibrationBlindPix::IsSinglePheFitOK() const
464{
465 return TESTBIT(fFlags,kSinglePheFitOK);
466}
467
468// --------------------------------------------------------------------------
469//
470// Test bit kFluxInsidePlexiglassAvailable
471//
472const Bool_t MCalibrationBlindPix::IsFluxInsidePlexiglassAvailable() const
473{
474 return TESTBIT(fFlags,kFluxInsidePlexiglassAvailable);
475}
476
477
478// --------------------------------------------------------------------------
479//
480// Return kFALSE if IsChargeFitValid() is kFALSE
481//
482// Calculate fFluxInsidePlexiglass with the formula:
483// - fFluxInsidePlexiglass = fLambda
484// / GetCollEff()
485// / GetQE()
486// * GetAtt()
487// / fArea
488// - fFluxInsidePlexiglassVar = sqrt( fLambdaVar / ( fLambda * fLambda )
489// + GetQERelVar()
490// + GetCollEffRelVar()
491// + GetAttRelVar()
492// ) * fFluxInsidePlexiglass * * fFluxInsidePlexiglass
493//
494// If the fFluxInsidePlexiglass is smaller than 0., return kFALSE
495// If the Variance is smaller than 0., return kFALSE
496//
497// SetFluxInsidePlexiglassAvailable() and return kTRUE
498//
499Bool_t MCalibrationBlindPix::CalcFluxInsidePlexiglass()
500{
501
502 if (IsChargeFitValid())
503 return kFALSE;
504
505
506 //
507 // Start calculation of number of photons
508 // The blind pixel has exactly 100 mm^2 area (with negligible error),
509 //
510 fFluxInsidePlexiglass = fLambda / GetQE() * GetAtt() / GetCollEff() / fArea;
511
512 if (fFluxInsidePlexiglass < 0.)
513 return kFALSE;
514
515 fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetQERelVar() + GetAttRelVar() + GetCollEffRelVar() + GetAreaRelVar();
516
517 //
518 // Finish calculation of errors -> convert from relative variance to absolute variance
519 //
520 fFluxInsidePlexiglassVar *= fFluxInsidePlexiglass * fFluxInsidePlexiglass;
521
522 if (fFluxInsidePlexiglassVar < 0.)
523 return kFALSE;
524
525 SetFluxInsidePlexiglassAvailable(kTRUE);
526
527 *fLog << inf << GetDescriptor()
528 << ": BlindPix #" << fPixId << ": Flux [ph/mm^2] inside Plexiglass: "
529 << Form("%5.3f +- %5.3f",fFluxInsidePlexiglass, GetFluxInsidePlexiglassErr()) << endl;
530
531 return kTRUE;
532}
533
534void MCalibrationBlindPix::Print(Option_t *opt) const
535{
536
537 *fLog << all << GetDescriptor()
538 << Form("BlindPixel: %3i" ,GetPixId())
539 << Form(" Lambda: %4.2f +- %4.2f",GetLambda(),GetLambdaErr())
540 << Form(" Mu0: %4.2f +- %4.2f" ,GetMu0(), GetMu0Err())
541 << Form(" Mu1: %4.2f +- %4.2f" ,GetMu1(), GetMu1Err())
542 << Form(" Sigma0: %4.2f +- %4.2f",GetSigma0(),GetSigma0Err())
543 << Form(" Sigma1: %4.2f +- %4.2f",GetSigma1(),GetSigma1Err())
544 << endl;
545 *fLog << all
546 << " Pedestal Fit OK? :" << IsPedestalFitOK()
547 << Form(" Lambda (Check): %4.2f +- %4.2f",GetLambdaCheck(),GetLambdaCheckErr())
548 << endl;
549 *fLog << all
550 << " Flux available? :" << IsFluxInsidePlexiglassAvailable()
551 << Form(" Flux: %4.2f +- %4.2f",GetFluxInsidePlexiglass(),GetFluxInsidePlexiglassErr())
552 << endl;
553}
Note: See TracBrowser for help on using the repository browser.