source: trunk/MagicSoft/Mars/mcalib/MCalibrationChargeBlindPix.cc@ 4203

Last change on this file since 4203 was 4203, checked in by gaug, 20 years ago
*** empty log message ***
File size: 13.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-2004
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MCalibrationChargeBlindPix
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// - gkBlindPixelArea: 100 mm^2
36// - Average QE of Blind Pixel:
37// gkBlindPixelQEGreen: 0.154
38// gkBlindPixelQEBlue : 0.226
39// gkBlindPixelQEUV : 0.247
40// gkBlindPixelQECT1 : 0.247
41// - Average QE Error of Blind Pixel:
42// gkBlindPixelQEGreenErr: 0.015;
43// gkBlindPixelQEBlueErr : 0.02;
44// gkBlindPixelQEUVErr : 0.02;
45// gkBlindPixelQECT1Err : 0.02;
46// - Attenuation factor Blind Pixel:
47// gkBlindPixelAttGreen : 1.97;
48// gkBlindPixelAttBlue : 1.96;
49// gkBlindPixelAttUV : 1.95;
50// gkBlindPixelAttCT1 : 1.95;
51//
52//
53/////////////////////////////////////////////////////////////////////////////
54#include "MCalibrationChargeBlindPix.h"
55#include "MCalibrationCam.h"
56
57#include <TH1.h>
58
59#include "MLog.h"
60#include "MLogManip.h"
61
62ClassImp(MCalibrationChargeBlindPix);
63
64using namespace std;
65const Float_t MCalibrationChargeBlindPix::gkBlindPixelArea = 100;
66const Float_t MCalibrationChargeBlindPix::gkBlindPixelAttGreen = 1.97;
67const Float_t MCalibrationChargeBlindPix::gkBlindPixelAttBlue = 1.96;
68const Float_t MCalibrationChargeBlindPix::gkBlindPixelAttUV = 1.95;
69const Float_t MCalibrationChargeBlindPix::gkBlindPixelAttCT1 = 1.95;
70const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEGreen = 0.154;
71const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEBlue = 0.226;
72const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUV = 0.247;
73const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECT1 = 0.247;
74const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEGreenErr = 0.015;
75const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEBlueErr = 0.02;
76const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUVErr = 0.02;
77const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECT1Err = 0.02;
78// --------------------------------------------------------------------------
79//
80// Default Constructor.
81//
82// Calls:
83// - Clear()
84//
85MCalibrationChargeBlindPix::MCalibrationChargeBlindPix(const char *name, const char *title)
86{
87
88 fName = name ? name : "MCalibrationChargeBlindPix";
89 fTitle = title ? title : "Container of the fit results of the blind pixel";
90
91 Clear();
92}
93
94
95// ------------------------------------------------------------------------
96//
97// Sets:
98// - all flags to kFALSE
99// - all variables to -1.
100// - the fColor to MCalibrationCam::kNONE
101//
102// Calls:
103// - MCalibrationChargePix::Clear()
104//
105void MCalibrationChargeBlindPix::Clear(Option_t *o)
106{
107
108 fFluxInsidePlexiglass = -1.;
109 fFluxInsidePlexiglassVar = -1.;
110 fLambda = -1.;
111 fLambdaCheck = -1.;
112 fLambdaVar = -1.;
113 fMu0 = -1.;
114 fMu0Err = -1.;
115 fMu1 = -1.;
116 fMu1Err = -1.;
117 fSigma0 = -1.;
118 fSigma0Err = -1.;
119 fSigma1 = -1.;
120 fSigma1Err = -1.;
121
122 SetOscillating ( kFALSE );
123 SetExcluded ( kFALSE );
124 SetChargeFitValid ( kFALSE );
125 SetPedestalFitOK ( kFALSE );
126 SetSinglePheFitOK ( kFALSE );
127 SetFluxInsidePlexiglassAvailable ( kFALSE );
128
129 SetColor(MCalibrationCam::kNONE);
130
131 MCalibrationChargePix::Clear();
132}
133
134void MCalibrationChargeBlindPix::SetFluxInsidePlexiglassAvailable( const Bool_t b)
135{
136 b ? SETBIT(fFlags,kFluxInsidePlexiglassAvailable) : CLRBIT(fFlags,kFluxInsidePlexiglassAvailable);
137}
138
139// --------------------------------------------------------------------------
140//
141// Set the Oscillating Bit from outside
142//
143void MCalibrationChargeBlindPix::SetOscillating( const Bool_t b)
144{
145 b ? SETBIT(fFlags,kOscillating) : CLRBIT(fFlags,kOscillating);
146}
147
148// --------------------------------------------------------------------------
149//
150// Set the ChargeFitValid Bit from outside
151//
152void MCalibrationChargeBlindPix::SetChargeFitValid( const Bool_t b)
153{
154 b ? SETBIT(fFlags,kChargeFitValid) : CLRBIT(fFlags,kChargeFitValid);
155}
156
157// --------------------------------------------------------------------------
158//
159// Set the PedestalFitValid Bit from outside
160//
161void MCalibrationChargeBlindPix::SetPedestalFitOK( const Bool_t b)
162{
163 b ? SETBIT(fFlags,kPedestalFitOK) : CLRBIT(fFlags,kPedestalFitOK);
164}
165
166// --------------------------------------------------------------------------
167//
168// Set the SinglePheFitValid Bit from outside
169//
170void MCalibrationChargeBlindPix::SetSinglePheFitOK( const Bool_t b)
171{
172 b ? SETBIT(fFlags,kSinglePheFitOK) : CLRBIT(fFlags,kSinglePheFitOK);
173}
174
175// --------------------------------------------------------------------------
176//
177// Return -1 if fFluxInsidePlexiglassVar is smaller than 0.
178// Return square root of fFluxInsidePlexiglassVar
179//
180Float_t MCalibrationChargeBlindPix::GetFluxInsidePlexiglassErr() const
181{
182 if (fFluxInsidePlexiglassVar < 0.)
183 return -1.;
184
185 return TMath::Sqrt(fFluxInsidePlexiglassVar);
186}
187
188// --------------------------------------------------------------------------
189//
190// Return -1 if fFluxInsidePlexiglassVar is smaller than 0.
191// Return -1 if fFluxInsidePlexiglass is 0.
192// Return fFluxInsidePlexiglassVar / fFluxInsidePlexiglass^2
193//
194Float_t MCalibrationChargeBlindPix::GetFluxInsidePlexiglassRelVar() const
195{
196 if (fFluxInsidePlexiglassVar < 0.)
197 return -1.;
198
199 if (fFluxInsidePlexiglass == 0.)
200 return -1.;
201
202 return fFluxInsidePlexiglassVar / (fFluxInsidePlexiglass * fFluxInsidePlexiglass) ;
203}
204
205// --------------------------------------------------------------------------
206//
207// Return -1 if fLambdaVar is smaller than 0.
208// Return square root of fLambdaVar
209//
210Float_t MCalibrationChargeBlindPix::GetLambdaErr() const
211{
212 if (fLambdaVar < 0.)
213 return -1.;
214
215 return TMath::Sqrt(fLambdaVar);
216}
217
218// --------------------------------------------------------------------------
219//
220// Return -1 if fLambdaVar is smaller than 0.
221// Return -1 if fLambda is 0.
222// Return fLambdaVar / (fLambda * fLambda )
223//
224Float_t MCalibrationChargeBlindPix::GetLambdaRelVar() const
225{
226 if (fLambdaVar < 0.)
227 return -1.;
228
229 if (fLambda == 0.)
230 return -1.;
231
232 return fLambdaVar / fLambda / fLambda ;
233}
234
235// --------------------------------------------------------------------------
236//
237// Return -1 if gkBlindPixelQEGreenErr is smaller than 0.
238// Return -1 if gkBlindPixelQEGreen is 0.
239// Return gkBlindPixelQEGreenErr^2 / (gkBlindPixelQEGreen^2 )
240//
241const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEGreenRelVar() const
242{
243 if (gkBlindPixelQEGreenErr < 0.)
244 return -1.;
245
246 if (gkBlindPixelQEGreen == 0.)
247 return -1.;
248
249 return gkBlindPixelQEGreenErr * gkBlindPixelQEGreenErr / gkBlindPixelQEGreen / gkBlindPixelQEGreen ;
250}
251
252// --------------------------------------------------------------------------
253//
254// Return -1 if gkBlindPixelQEBlueErr is smaller than 0.
255// Return -1 if gkBlindPixelQEBlue is 0.
256// Return gkBlindPixelQEBlueErr^2 / gkBlindPixelQEBlue^2
257//
258const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEBlueRelVar() const
259{
260 if (gkBlindPixelQEBlueErr < 0.)
261 return -1.;
262
263 if (gkBlindPixelQEBlue == 0.)
264 return -1.;
265
266 return gkBlindPixelQEBlueErr * gkBlindPixelQEBlueErr / gkBlindPixelQEBlue / gkBlindPixelQEBlue ;
267}
268
269// --------------------------------------------------------------------------
270//
271// Return -1 if gkBlindPixelQEUVErr is smaller than 0.
272// Return -1 if gkBlindPixelQEUV is 0.
273// Return gkBlindPixelQEUVErr ^2 / gkBlindPixelQEUV^2
274//
275const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEUVRelVar() const
276{
277 if (gkBlindPixelQEUVErr < 0.)
278 return -1.;
279
280 if (gkBlindPixelQEUV == 0.)
281 return -1.;
282
283 return gkBlindPixelQEUVErr * gkBlindPixelQEUVErr / gkBlindPixelQEUV / gkBlindPixelQEUV ;
284}
285
286// --------------------------------------------------------------------------
287//
288// Return -1 if gkBlindPixelQECT1Err is smaller than 0.
289// Return -1 if gkBlindPixelQECT1 is 0.
290// Return gkBlindPixelQECT1Err ^2 / gkBlindPixelQECT1^2
291//
292const Float_t MCalibrationChargeBlindPix::GetBlindPixelQECT1RelVar() const
293{
294 if (gkBlindPixelQECT1Err < 0.)
295 return -1.;
296
297 if (gkBlindPixelQECT1 == 0.)
298 return -1.;
299
300 return gkBlindPixelQECT1Err * gkBlindPixelQECT1Err / gkBlindPixelQECT1 / gkBlindPixelQECT1 ;
301}
302
303// --------------------------------------------------------------------------
304//
305// Test bit kChargeFitValid
306//
307Bool_t MCalibrationChargeBlindPix::IsChargeFitValid() const
308{
309 return TESTBIT(fFlags,kChargeFitValid);
310}
311
312// --------------------------------------------------------------------------
313//
314// Test bit kOscillating
315//
316Bool_t MCalibrationChargeBlindPix::IsOscillating() const
317{
318 return TESTBIT(fFlags,kOscillating);
319}
320
321// --------------------------------------------------------------------------
322//
323// Test bit kPedestalFitValid
324//
325Bool_t MCalibrationChargeBlindPix::IsPedestalFitOK() const
326{
327 return TESTBIT(fFlags,kPedestalFitOK);
328}
329
330// --------------------------------------------------------------------------
331//
332// Test bit kSinglePheFitValid
333//
334Bool_t MCalibrationChargeBlindPix::IsSinglePheFitOK() const
335{
336 return TESTBIT(fFlags,kSinglePheFitOK);
337}
338
339// --------------------------------------------------------------------------
340//
341// Test bit kFluxInsidePlexiglassAvailable
342//
343Bool_t MCalibrationChargeBlindPix::IsFluxInsidePlexiglassAvailable() const
344{
345 return TESTBIT(fFlags,kFluxInsidePlexiglassAvailable);
346}
347
348
349// --------------------------------------------------------------------------
350//
351// Return kFALSE if IsChargeFitValid() is kFALSE
352//
353// Calculate fFluxInsidePlexiglass with the formula:
354// - fFluxInsidePlexiglass = fLambda * gkBlindPixelArea / gkBlindPixelQE * 10**gkBlindPixelAtt
355// - fFluxInsidePlexiglassVar = sqrt( fLambdaVar / ( fLambda * fLambda )
356// + ( gkBlindPixelQEErr * gkBlindPixelQEErr / gkBlindPixelQE / gkBlindPixelQE )
357// ) * fFluxInsidePlexiglass * * fFluxInsidePlexiglass
358//
359// If the fFluxInsidePlexiglass is smaller than 0., return kFALSE
360// If the Variance is smaller than 0., return kFALSE
361//
362// SetFluxInsidePlexiglassAvailable() and return kTRUE
363//
364Bool_t MCalibrationChargeBlindPix::CalcFluxInsidePlexiglass()
365{
366
367 if (IsChargeFitValid())
368 return kFALSE;
369
370
371 //
372 // Start calculation of number of photons
373 // The blind pixel has exactly 100 mm^2 area (with negligible error),
374 //
375 switch (fColor)
376 {
377 case MCalibrationCam::kGREEN:
378 fFluxInsidePlexiglass = fLambda / gkBlindPixelQEGreen * TMath::Power(10,gkBlindPixelAttGreen);
379 // attenuation has negligible error
380 fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQEGreenRelVar();
381 break;
382 case MCalibrationCam::kBLUE:
383 fFluxInsidePlexiglass = fLambda / gkBlindPixelQEBlue * TMath::Power(10,gkBlindPixelAttBlue);
384 // attenuation has negligible error
385 fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQEBlueRelVar();
386 break;
387 case MCalibrationCam::kUV:
388 fFluxInsidePlexiglass = fLambda / gkBlindPixelQEUV * TMath::Power(10,gkBlindPixelAttUV);
389 // attenuation has negligible error
390 fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQEUVRelVar();
391 break;
392 case MCalibrationCam::kCT1:
393 default:
394 fFluxInsidePlexiglass = fLambda / gkBlindPixelQECT1 * TMath::Power(10,gkBlindPixelAttCT1);
395 // attenuation has negligible error
396 fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQECT1RelVar();
397 break;
398 }
399
400 fFluxInsidePlexiglass /= gkBlindPixelArea;
401
402 //
403 // Finish calculation of errors -> convert from relative variance to absolute variance
404 //
405 fFluxInsidePlexiglassVar *= fFluxInsidePlexiglass * fFluxInsidePlexiglass;
406
407 if (fFluxInsidePlexiglass < 0.)
408 return kFALSE;
409
410 if (fFluxInsidePlexiglassVar < 0.)
411 return kFALSE;
412
413 SetFluxInsidePlexiglassAvailable(kTRUE);
414
415 *fLog << inf << endl;
416 *fLog << inf << GetDescriptor()
417 << ": Photon flux [ph/mm^2] inside Plexiglass: "
418 << Form("%5.3f%s%5.3f",fFluxInsidePlexiglass," +- ",GetFluxInsidePlexiglassErr()) << endl;
419
420 return kTRUE;
421}
422
423
424
425
426
427
428
429
430
431
Note: See TracBrowser for help on using the repository browser.