source: trunk/MagicSoft/Mars/mcalib/MCalibrationBlindPix.cc@ 4991

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