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

Last change on this file since 4231 was 4231, checked in by gaug, 20 years ago
*** empty log message ***
File size: 19.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::gkBlindPixelQEUnCoatedGreen = 0.154;
71const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUnCoatedBlue = 0.226;
72const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUnCoatedUV = 0.247;
73const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUnCoatedCT1 = 0.247;
74const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUnCoatedGreenErr = 0.015;
75const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUnCoatedBlueErr = 0.02;
76const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUnCoatedUVErr = 0.02;
77const Float_t MCalibrationChargeBlindPix::gkBlindPixelQEUnCoatedCT1Err = 0.02;
78const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedGreen = 0.192;
79const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedBlue = 0.27;
80const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedUV = 0.285;
81const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedCT1 = 0.285;
82const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedGreenErr = 0.05;
83const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedBlueErr = 0.07;
84const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedUVErr = 0.07;
85const Float_t MCalibrationChargeBlindPix::gkBlindPixelQECoatedCT1Err = 0.07;
86const Float_t MCalibrationChargeBlindPix::gkBlindPixelCollectionEff = 0.90;
87const Float_t MCalibrationChargeBlindPix::gkBlindPixelCollectionEffErr = 0.02;
88// --------------------------------------------------------------------------
89//
90// Default Constructor.
91//
92// Calls:
93// - Clear()
94// - SetCoated()
95//
96MCalibrationChargeBlindPix::MCalibrationChargeBlindPix(const char *name, const char *title)
97{
98
99 fName = name ? name : "MCalibrationChargeBlindPix";
100 fTitle = title ? title : "Container of the fit results of the blind pixel";
101
102 SetCoated();
103 Clear();
104}
105
106
107// ------------------------------------------------------------------------
108//
109// Sets:
110// - all flags to kFALSE
111// - all variables to -1.
112// - the fColor to MCalibrationCam::kNONE
113//
114// Calls:
115// - MCalibrationChargePix::Clear()
116//
117void MCalibrationChargeBlindPix::Clear(Option_t *o)
118{
119
120 fFluxInsidePlexiglass = -1.;
121 fFluxInsidePlexiglassVar = -1.;
122 fLambda = -1.;
123 fLambdaCheck = -1.;
124 fLambdaVar = -1.;
125 fMu0 = -1.;
126 fMu0Err = -1.;
127 fMu1 = -1.;
128 fMu1Err = -1.;
129 fSigma0 = -1.;
130 fSigma0Err = -1.;
131 fSigma1 = -1.;
132 fSigma1Err = -1.;
133
134 SetOscillating ( kFALSE );
135 SetExcluded ( kFALSE );
136 SetChargeFitValid ( kFALSE );
137 SetPedestalFitOK ( kFALSE );
138 SetSinglePheFitOK ( kFALSE );
139 SetFluxInsidePlexiglassAvailable ( kFALSE );
140
141 SetColor(MCalibrationCam::kNONE);
142
143 MCalibrationChargePix::Clear();
144}
145
146void MCalibrationChargeBlindPix::SetFluxInsidePlexiglassAvailable( const Bool_t b)
147{
148 b ? SETBIT(fFlags,kFluxInsidePlexiglassAvailable) : CLRBIT(fFlags,kFluxInsidePlexiglassAvailable);
149}
150
151// --------------------------------------------------------------------------
152//
153// Set the Coated Bit from outside
154//
155void MCalibrationChargeBlindPix::SetCoated( const Bool_t b)
156{
157 b ? SETBIT(fFlags,kCoated) : CLRBIT(fFlags,kCoated);
158}
159
160// --------------------------------------------------------------------------
161//
162// Set the Oscillating Bit from outside
163//
164void MCalibrationChargeBlindPix::SetOscillating( const Bool_t b)
165{
166 b ? SETBIT(fFlags,kOscillating) : CLRBIT(fFlags,kOscillating);
167}
168
169// --------------------------------------------------------------------------
170//
171// Set the ChargeFitValid Bit from outside
172//
173void MCalibrationChargeBlindPix::SetChargeFitValid( const Bool_t b)
174{
175 b ? SETBIT(fFlags,kChargeFitValid) : CLRBIT(fFlags,kChargeFitValid);
176}
177
178// --------------------------------------------------------------------------
179//
180// Set the PedestalFitValid Bit from outside
181//
182void MCalibrationChargeBlindPix::SetPedestalFitOK( const Bool_t b)
183{
184 b ? SETBIT(fFlags,kPedestalFitOK) : CLRBIT(fFlags,kPedestalFitOK);
185}
186
187// --------------------------------------------------------------------------
188//
189// Set the SinglePheFitValid Bit from outside
190//
191void MCalibrationChargeBlindPix::SetSinglePheFitOK( const Bool_t b)
192{
193 b ? SETBIT(fFlags,kSinglePheFitOK) : CLRBIT(fFlags,kSinglePheFitOK);
194}
195
196// --------------------------------------------------------------------------
197//
198// Return -1 if fFluxInsidePlexiglassVar is smaller than 0.
199// Return square root of fFluxInsidePlexiglassVar
200//
201Float_t MCalibrationChargeBlindPix::GetFluxInsidePlexiglassErr() const
202{
203 if (fFluxInsidePlexiglassVar < 0.)
204 return -1.;
205
206 return TMath::Sqrt(fFluxInsidePlexiglassVar);
207}
208
209// --------------------------------------------------------------------------
210//
211// Return -1 if fFluxInsidePlexiglassVar is smaller than 0.
212// Return -1 if fFluxInsidePlexiglass is 0.
213// Return fFluxInsidePlexiglassVar / fFluxInsidePlexiglass^2
214//
215Float_t MCalibrationChargeBlindPix::GetFluxInsidePlexiglassRelVar() const
216{
217 if (fFluxInsidePlexiglassVar < 0.)
218 return -1.;
219
220 if (fFluxInsidePlexiglass == 0.)
221 return -1.;
222
223 return fFluxInsidePlexiglassVar / (fFluxInsidePlexiglass * fFluxInsidePlexiglass) ;
224}
225
226// --------------------------------------------------------------------------
227//
228// Return -1 if fLambdaVar is smaller than 0.
229// Return square root of fLambdaVar
230//
231Float_t MCalibrationChargeBlindPix::GetLambdaErr() const
232{
233 if (fLambdaVar < 0.)
234 return -1.;
235
236 return TMath::Sqrt(fLambdaVar);
237}
238
239// --------------------------------------------------------------------------
240//
241// Return -1 if fLambdaVar is smaller than 0.
242// Return -1 if fLambda is 0.
243// Return fLambdaVar / (fLambda * fLambda )
244//
245Float_t MCalibrationChargeBlindPix::GetLambdaRelVar() const
246{
247 if (fLambdaVar < 0.)
248 return -1.;
249
250 if (fLambda == 0.)
251 return -1.;
252
253 return fLambdaVar / fLambda / fLambda ;
254}
255
256// --------------------------------------------------------------------------
257//
258// Return -1 if gkBlindPixelQEGreenErr is smaller than 0.
259// Return -1 if gkBlindPixelQEGreen is 0.
260// Return gkBlindPixelQEGreenErr^2 / (gkBlindPixelQEGreen^2 )
261//
262const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEGreen() const
263{
264
265 if (IsCoated())
266 {
267 if (gkBlindPixelQECoatedGreen < 0.)
268 return -1.;
269
270 return gkBlindPixelQECoatedGreen;
271 }
272 else
273 {
274 if (gkBlindPixelQEUnCoatedGreen < 0.)
275 return -1.;
276
277 return gkBlindPixelQEUnCoatedGreen;
278 }
279
280}
281
282// --------------------------------------------------------------------------
283//
284// Return -1 if gkBlindPixelQEBlueErr is smaller than 0.
285// Return -1 if gkBlindPixelQEBlue is 0.
286// Return gkBlindPixelQEBlueErr^2 / gkBlindPixelQEBlue^2
287//
288const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEBlue() const
289{
290 if (IsCoated())
291 {
292 if (gkBlindPixelQECoatedBlue < 0.)
293 return -1.;
294
295 return gkBlindPixelQECoatedBlue;
296 }
297 else
298 {
299 if (gkBlindPixelQEUnCoatedBlue < 0.)
300 return -1.;
301
302 return gkBlindPixelQEUnCoatedBlue;
303 }
304}
305
306// --------------------------------------------------------------------------
307//
308// Return -1 if gkBlindPixelQEUVErr is smaller than 0.
309// Return -1 if gkBlindPixelQEUV is 0.
310// Return gkBlindPixelQEUVErr ^2 / gkBlindPixelQEUV^2
311//
312const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEUV() const
313{
314
315 if (IsCoated())
316 {
317 if (gkBlindPixelQECoatedUV < 0.)
318 return -1.;
319
320 return gkBlindPixelQECoatedUV;
321 }
322 else
323 {
324 if (gkBlindPixelQEUnCoatedUV < 0.)
325 return -1.;
326
327
328 return gkBlindPixelQEUnCoatedUV;
329 }
330}
331
332// --------------------------------------------------------------------------
333//
334// Return -1 if gkBlindPixelQECT1Err is smaller than 0.
335// Return -1 if gkBlindPixelQECT1 is 0.
336// Return gkBlindPixelQECT1Err ^2 / gkBlindPixelQECT1^2
337//
338const Float_t MCalibrationChargeBlindPix::GetBlindPixelQECT1() const
339{
340
341 if (IsCoated())
342 {
343 if (gkBlindPixelQECoatedCT1 < 0.)
344 return -1.;
345
346 return gkBlindPixelQECoatedCT1;
347 }
348 else
349 {
350 if (gkBlindPixelQEUnCoatedCT1 < 0.)
351 return -1.;
352
353 return gkBlindPixelQEUnCoatedCT1;
354 }
355}
356
357// --------------------------------------------------------------------------
358//
359// Return -1 if gkBlindPixelQEGreenErr is smaller than 0.
360// Return -1 if gkBlindPixelQEGreen is 0.
361// Return gkBlindPixelQEGreenErr^2 / (gkBlindPixelQEGreen^2 )
362//
363const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEGreenRelVar() const
364{
365
366 if (IsCoated())
367 {
368 if (gkBlindPixelQECoatedGreenErr < 0.)
369 return -1.;
370
371 if (gkBlindPixelQECoatedGreen == 0.)
372 return -1.;
373
374 return gkBlindPixelQECoatedGreenErr * gkBlindPixelQECoatedGreenErr
375 / gkBlindPixelQECoatedGreen / gkBlindPixelQECoatedGreen ;
376 }
377 else
378 {
379 if (gkBlindPixelQEUnCoatedGreenErr < 0.)
380 return -1.;
381
382 if (gkBlindPixelQEUnCoatedGreen == 0.)
383 return -1.;
384
385 return gkBlindPixelQEUnCoatedGreenErr * gkBlindPixelQEUnCoatedGreenErr
386 / gkBlindPixelQEUnCoatedGreen / gkBlindPixelQEUnCoatedGreen ;
387 }
388
389}
390
391// --------------------------------------------------------------------------
392//
393// Return -1 if gkBlindPixelQEBlueErr is smaller than 0.
394// Return -1 if gkBlindPixelQEBlue is 0.
395// Return gkBlindPixelQEBlueErr^2 / gkBlindPixelQEBlue^2
396//
397const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEBlueRelVar() const
398{
399 if (IsCoated())
400 {
401 if (gkBlindPixelQECoatedBlueErr < 0.)
402 return -1.;
403
404 if (gkBlindPixelQECoatedBlue == 0.)
405 return -1.;
406
407 return gkBlindPixelQECoatedBlueErr * gkBlindPixelQECoatedBlueErr
408 / gkBlindPixelQECoatedBlue / gkBlindPixelQECoatedBlue ;
409 }
410 else
411 {
412 if (gkBlindPixelQEUnCoatedBlueErr < 0.)
413 return -1.;
414
415 if (gkBlindPixelQEUnCoatedBlue == 0.)
416 return -1.;
417
418 return gkBlindPixelQEUnCoatedBlueErr * gkBlindPixelQEUnCoatedBlueErr
419 / gkBlindPixelQEUnCoatedBlue / gkBlindPixelQEUnCoatedBlue ;
420 }
421}
422
423// --------------------------------------------------------------------------
424//
425// Return -1 if gkBlindPixelQEUVErr is smaller than 0.
426// Return -1 if gkBlindPixelQEUV is 0.
427// Return gkBlindPixelQEUVErr ^2 / gkBlindPixelQEUV^2
428//
429const Float_t MCalibrationChargeBlindPix::GetBlindPixelQEUVRelVar() const
430{
431
432 if (IsCoated())
433 {
434 if (gkBlindPixelQECoatedUVErr < 0.)
435 return -1.;
436
437 if (gkBlindPixelQECoatedUV == 0.)
438 return -1.;
439
440 return gkBlindPixelQECoatedUVErr * gkBlindPixelQECoatedUVErr
441 / gkBlindPixelQECoatedUV / gkBlindPixelQECoatedUV ;
442 }
443 else
444 {
445 if (gkBlindPixelQEUnCoatedUVErr < 0.)
446 return -1.;
447
448 if (gkBlindPixelQEUnCoatedUV == 0.)
449 return -1.;
450
451 return gkBlindPixelQEUnCoatedUVErr * gkBlindPixelQEUnCoatedUVErr
452 / gkBlindPixelQEUnCoatedUV / gkBlindPixelQEUnCoatedUV ;
453 }
454}
455
456// --------------------------------------------------------------------------
457//
458// Return -1 if gkBlindPixelQECT1Err is smaller than 0.
459// Return -1 if gkBlindPixelQECT1 is 0.
460// Return gkBlindPixelQECT1Err ^2 / gkBlindPixelQECT1^2
461//
462const Float_t MCalibrationChargeBlindPix::GetBlindPixelQECT1RelVar() const
463{
464
465 if (IsCoated())
466 {
467 if (gkBlindPixelQECoatedCT1Err < 0.)
468 return -1.;
469
470 if (gkBlindPixelQECoatedCT1 == 0.)
471 return -1.;
472
473 return gkBlindPixelQECoatedCT1Err * gkBlindPixelQECoatedCT1Err
474 / gkBlindPixelQECoatedCT1 / gkBlindPixelQECoatedCT1 ;
475 }
476 else
477 {
478 if (gkBlindPixelQEUnCoatedCT1Err < 0.)
479 return -1.;
480
481 if (gkBlindPixelQEUnCoatedCT1 == 0.)
482 return -1.;
483
484 return gkBlindPixelQEUnCoatedCT1Err * gkBlindPixelQEUnCoatedCT1Err
485 / gkBlindPixelQEUnCoatedCT1 / gkBlindPixelQEUnCoatedCT1 ;
486 }
487}
488
489// --------------------------------------------------------------------------
490//
491// Return gkBlindPixelCollectionEffErr^2 / (gkBlindPixelCollectionEff^2 )
492//
493const Float_t MCalibrationChargeBlindPix::GetBlindPixelCollectionEffRelVar() const
494{
495
496 return gkBlindPixelCollectionEffErr * gkBlindPixelCollectionEffErr
497 / gkBlindPixelCollectionEff / gkBlindPixelCollectionEff ;
498}
499
500// --------------------------------------------------------------------------
501//
502// Test bit kChargeFitValid
503//
504Bool_t MCalibrationChargeBlindPix::IsChargeFitValid() const
505{
506 return TESTBIT(fFlags,kChargeFitValid);
507}
508
509// --------------------------------------------------------------------------
510//
511// Test bit kCoated
512//
513Bool_t MCalibrationChargeBlindPix::IsCoated() const
514{
515 return TESTBIT(fFlags,kCoated);
516}
517
518// --------------------------------------------------------------------------
519//
520// Test bit kOscillating
521//
522Bool_t MCalibrationChargeBlindPix::IsOscillating() const
523{
524 return TESTBIT(fFlags,kOscillating);
525}
526
527// --------------------------------------------------------------------------
528//
529// Test bit kPedestalFitValid
530//
531Bool_t MCalibrationChargeBlindPix::IsPedestalFitOK() const
532{
533 return TESTBIT(fFlags,kPedestalFitOK);
534}
535
536// --------------------------------------------------------------------------
537//
538// Test bit kSinglePheFitValid
539//
540Bool_t MCalibrationChargeBlindPix::IsSinglePheFitOK() const
541{
542 return TESTBIT(fFlags,kSinglePheFitOK);
543}
544
545// --------------------------------------------------------------------------
546//
547// Test bit kFluxInsidePlexiglassAvailable
548//
549Bool_t MCalibrationChargeBlindPix::IsFluxInsidePlexiglassAvailable() const
550{
551 return TESTBIT(fFlags,kFluxInsidePlexiglassAvailable);
552}
553
554
555// --------------------------------------------------------------------------
556//
557// Return kFALSE if IsChargeFitValid() is kFALSE
558//
559// Calculate fFluxInsidePlexiglass with the formula:
560// - fFluxInsidePlexiglass = fLambda
561// / GetBlindPixelCollectionEff()
562// / GetBlindPixelQE()
563// * 10**gkBlindPixelAtt[color]
564// / gkBlindPixelArea
565// - fFluxInsidePlexiglassVar = sqrt( fLambdaVar / ( fLambda * fLambda )
566// + GetBlindPixelQERelVar()
567// + GetBlindPixelCollectionEffRelVar()
568// ) * fFluxInsidePlexiglass * * fFluxInsidePlexiglass
569//
570// If the fFluxInsidePlexiglass is smaller than 0., return kFALSE
571// If the Variance is smaller than 0., return kFALSE
572//
573// SetFluxInsidePlexiglassAvailable() and return kTRUE
574//
575Bool_t MCalibrationChargeBlindPix::CalcFluxInsidePlexiglass()
576{
577
578 if (IsChargeFitValid())
579 return kFALSE;
580
581
582 //
583 // Start calculation of number of photons
584 // The blind pixel has exactly 100 mm^2 area (with negligible error),
585 //
586 switch (fColor)
587 {
588 case MCalibrationCam::kGREEN:
589 fFluxInsidePlexiglass = fLambda / GetBlindPixelQEGreen() * TMath::Power(10,gkBlindPixelAttGreen);
590 // attenuation has negligible error
591 fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQEGreenRelVar();
592 break;
593 case MCalibrationCam::kBLUE:
594 fFluxInsidePlexiglass = fLambda / GetBlindPixelQEBlue() * TMath::Power(10,gkBlindPixelAttBlue);
595 // attenuation has negligible error
596 fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQEBlueRelVar();
597 break;
598 case MCalibrationCam::kUV:
599 fFluxInsidePlexiglass = fLambda / GetBlindPixelQEUV() * TMath::Power(10,gkBlindPixelAttUV);
600 // attenuation has negligible error
601 fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQEUVRelVar();
602 break;
603 case MCalibrationCam::kCT1:
604 default:
605 fFluxInsidePlexiglass = fLambda / GetBlindPixelQECT1() * TMath::Power(10,gkBlindPixelAttCT1);
606 // attenuation has negligible error
607 fFluxInsidePlexiglassVar = GetLambdaRelVar() + GetBlindPixelQECT1RelVar();
608 break;
609 }
610
611 fFluxInsidePlexiglass /= gkBlindPixelArea;
612 fFluxInsidePlexiglass /= gkBlindPixelCollectionEff;
613 //
614 // Finish calculation of errors -> convert from relative variance to absolute variance
615 //
616 fFluxInsidePlexiglassVar += GetBlindPixelCollectionEffRelVar();
617 fFluxInsidePlexiglassVar *= fFluxInsidePlexiglass * fFluxInsidePlexiglass;
618
619 if (fFluxInsidePlexiglass < 0.)
620 return kFALSE;
621
622 if (fFluxInsidePlexiglassVar < 0.)
623 return kFALSE;
624
625 SetFluxInsidePlexiglassAvailable(kTRUE);
626
627 *fLog << inf << endl;
628 *fLog << inf << GetDescriptor()
629 << ": Photon flux [ph/mm^2] inside Plexiglass: "
630 << Form("%5.3f%s%5.3f",fFluxInsidePlexiglass," +- ",GetFluxInsidePlexiglassErr()) << endl;
631
632 return kTRUE;
633}
634
635
636
637
638
639
640
641
642
643
Note: See TracBrowser for help on using the repository browser.