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

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