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

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