source: trunk/MagicSoft/Mars/mcalib/MCalibrationCam.cc@ 3065

Last change on this file since 3065 was 3065, checked in by gaug, 21 years ago
*** empty log message ***
File size: 32.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!
19! Author(s): Markus Gaug 11/2003 <mailto:markus@ifae.es>
20!
21! Copyright: MAGIC Software Development, 2000-2001
22!
23!
24\* ======================================================================== */
25
26/////////////////////////////////////////////////////////////////////////////
27//
28// MCalibrationCam
29//
30// Hold the whole Calibration results of the camera:
31//
32// 1) MCalibrationCam initializes a TClonesArray whose elements are
33// pointers to MCalibrationPix Containers
34// 2) It initializes a pointer to an MCalibrationBlindPix container
35// 3) It initializes a pointer to an MCalibrationPINDiode container
36//
37/////////////////////////////////////////////////////////////////////////////
38#include "MCalibrationCam.h"
39
40#include <TH2.h>
41#include <TCanvas.h>
42#include <TClonesArray.h>
43
44#include "MLog.h"
45#include "MLogManip.h"
46
47#include "MGeomCam.h"
48#include "MGeomPix.h"
49
50#include "MCalibrationPix.h"
51#include "MCalibrationConfig.h"
52#include "MCalibrationBlindPix.h"
53#include "MCalibrationPINDiode.h"
54
55#include "MHCalibrationPixel.h"
56
57ClassImp(MCalibrationCam);
58
59using namespace std;
60
61const Int_t MCalibrationCam::gkBlindPixelId = 559;
62const Int_t MCalibrationCam::gkPINDiodeId = 9999;
63
64// --------------------------------------------------------------------------
65//
66// Default constructor.
67//
68// Creates a TClonesArray of MCalibrationPix containers, initialized to 1 entry
69// Later, a call to MCalibrationCam::InitSize(Int_t size) has to be performed
70//
71// Creates an MCalibrationBlindPix container
72// Creates an MCalibrationPINDiode container
73//
74MCalibrationCam::MCalibrationCam(const char *name, const char *title)
75 : fOffsets(NULL),
76 fSlopes(NULL),
77 fOffvsSlope(NULL)
78{
79 fName = name ? name : "MCalibrationCam";
80 fTitle = title ? title : "Storage container for the Calibration Information in the camera";
81
82 fPixels = new TClonesArray("MCalibrationPix",1);
83 fBlindPixel = new MCalibrationBlindPix();
84 fPINDiode = new MCalibrationPINDiode();
85
86 Clear();
87}
88
89// --------------------------------------------------------------------------
90//
91// Delete the TClonesArray of MCalibrationPix containers
92// Delete the MCalibrationPINDiode and the MCalibrationBlindPix
93//
94// Delete the histograms if they exist
95//
96MCalibrationCam::~MCalibrationCam()
97{
98
99 //
100 // delete fPixels should delete all Objects stored inside
101 //
102 delete fPixels;
103 delete fBlindPixel;
104 delete fPINDiode;
105
106 if (fOffsets)
107 delete fOffsets;
108 if (fSlopes)
109 delete fSlopes;
110 if (fOffvsSlope)
111 delete fOffvsSlope;
112
113}
114
115// -------------------------------------------------------------------
116//
117// This function simply allocates memory via the ROOT command:
118// (TObject**) TStorage::ReAlloc(fCont, newSize * sizeof(TObject*),
119// fSize * sizeof(TObject*));
120// newSize corresponds to size in our case
121// fSize is the old size (in most cases: 1)
122//
123void MCalibrationCam::InitSize(const UInt_t i)
124{
125
126 //
127 // check if we have already initialized to size
128 //
129 if (CheckBounds(i))
130 return;
131
132 fPixels->ExpandCreate(i);
133
134}
135
136// --------------------------------------------------------------------------
137//
138// This function returns the current size of the TClonesArray
139// independently if the MCalibrationPix is filled with values or not.
140//
141// It is the size of the array fPixels.
142//
143Int_t MCalibrationCam::GetSize() const
144{
145 return fPixels->GetEntriesFast();
146}
147
148// --------------------------------------------------------------------------
149//
150// Check if position i is inside the current bounds of the TClonesArray
151//
152Bool_t MCalibrationCam::CheckBounds(Int_t i) const
153{
154 return i < GetSize();
155}
156
157
158// --------------------------------------------------------------------------
159//
160// Get i-th pixel (pixel number)
161//
162MCalibrationPix &MCalibrationCam::operator[](Int_t i)
163{
164 return *static_cast<MCalibrationPix*>(fPixels->UncheckedAt(i));
165}
166
167// --------------------------------------------------------------------------
168//
169// Get i-th pixel (pixel number)
170//
171MCalibrationPix &MCalibrationCam::operator[](Int_t i) const
172{
173 return *static_cast<MCalibrationPix*>(fPixels->UncheckedAt(i));
174}
175
176
177// --------------------------------------
178//
179void MCalibrationCam::Clear(Option_t *o)
180{
181
182 fPixels->ForEach(TObject, Clear)();
183 fBlindPixel->Clear();
184 fPINDiode->Clear();
185
186 fMeanPhotInsidePlexiglass = -1.;
187 fMeanPhotErrInsidePlexiglass = -1.;
188 fMeanPhotOutsidePlexiglass = -1.;
189 fMeanPhotErrOutsidePlexiglass = -1.;
190
191 fNumExcludedPixels = 0;
192
193 CLRBIT(fFlags,kBlindPixelMethodValid);
194 CLRBIT(fFlags,kPINDiodeMethodValid);
195 CLRBIT(fFlags,kNumPhotInsidePlexiglassAvailable);
196 CLRBIT(fFlags,kNumPhotOutsidePlexiglassAvailable);
197
198 return;
199}
200
201void MCalibrationCam::SetBlindPixelMethodValid(const Bool_t b)
202{
203
204 if (b)
205 SETBIT(fFlags, kBlindPixelMethodValid);
206 else
207 CLRBIT(fFlags, kBlindPixelMethodValid);
208
209}
210
211void MCalibrationCam::SetPINDiodeMethodValid(const Bool_t b)
212{
213
214 if (b)
215 SETBIT(fFlags, kPINDiodeMethodValid);
216 else
217 CLRBIT(fFlags, kPINDiodeMethodValid);
218
219
220}
221
222Bool_t MCalibrationCam::IsBlindPixelMethodValid() const
223{
224 return TESTBIT(fFlags,kBlindPixelMethodValid);
225}
226
227Bool_t MCalibrationCam::IsPINDiodeMethodValid() const
228{
229 return TESTBIT(fFlags,kPINDiodeMethodValid);
230}
231
232
233Bool_t MCalibrationCam::IsNumPhotInsidePlexiglassAvailable() const
234{
235 return TESTBIT(fFlags,kNumPhotInsidePlexiglassAvailable);
236}
237
238Bool_t MCalibrationCam::IsNumPhotOutsidePlexiglassAvailable() const
239{
240 return TESTBIT(fFlags,kNumPhotOutsidePlexiglassAvailable);
241}
242
243
244
245// --------------------------------------------------------------------------
246//
247// Print first the well fitted pixels
248// and then the ones which are not FitValid
249//
250void MCalibrationCam::Print(Option_t *o) const
251{
252
253 *fLog << all << GetDescriptor() << ":" << endl;
254 int id = 0;
255
256 *fLog << all << "Succesfully calibrated pixels:" << endl;
257 *fLog << all << endl;
258
259 TIter Next(fPixels);
260 MCalibrationPix *pix;
261 while ((pix=(MCalibrationPix*)Next()))
262 {
263
264 if (pix->IsChargeValid() && !pix->IsExcluded())
265 {
266
267 *fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- "
268 << pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- "
269 << pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigmaCharge()
270 << " Nr Phe's: " << pix->GetPheFFactorMethod() << endl;
271 id++;
272 }
273 }
274
275 *fLog << all << id << " succesful pixels :-))" << endl;
276 id = 0;
277
278 *fLog << all << endl;
279 *fLog << all << "Pixels with errors:" << endl;
280 *fLog << all << endl;
281
282 TIter Next2(fPixels);
283 while ((pix=(MCalibrationPix*)Next2()))
284 {
285
286 if (!pix->IsChargeValid() && !pix->IsExcluded())
287 {
288
289 *fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- "
290 << pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- "
291 << pix->GetSigmaCharge() << " Reduced Sigma: " << pix->GetRSigmaCharge() << endl;
292 id++;
293 }
294 }
295 *fLog << all << id << " pixels with errors :-((" << endl;
296
297 *fLog << all << endl;
298 *fLog << all << "Excluded pixels:" << endl;
299 *fLog << all << endl;
300
301 TIter Next3(fPixels);
302 while ((pix=(MCalibrationPix*)Next3()))
303 if (pix->IsExcluded())
304 *fLog << all << pix->GetPixId() << endl;
305
306 *fLog << all << fNumExcludedPixels << " excluded pixels " << endl;
307}
308
309// --------------------------------------------------------------------------
310//
311// Return true if pixel is inside bounds of the TClonesArray fPixels
312//
313Bool_t MCalibrationCam::IsPixelUsed(Int_t idx) const
314{
315 if (!CheckBounds(idx))
316 return kFALSE;
317
318 return kTRUE;
319}
320
321// --------------------------------------------------------------------------
322//
323// Return true if pixel has already been fitted once (independent of the result)
324//
325Bool_t MCalibrationCam::IsPixelFitted(Int_t idx) const
326{
327
328 if (!CheckBounds(idx))
329 return kFALSE;
330
331 return (*this)[idx].IsFitted();
332}
333
334// --------------------------------------------------------------------------
335//
336// Sets the user ranges of all histograms such that
337// empty bins at the edges are not used. Additionally, it rebins the
338// histograms such that in total, 50 bins are used.
339//
340void MCalibrationCam::CutEdges()
341{
342
343 fBlindPixel->GetHist()->CutAllEdges();
344 fPINDiode->GetHist()->CutAllEdges();
345
346 TIter Next(fPixels);
347 MCalibrationPix *pix;
348 while ((pix=(MCalibrationPix*)Next()))
349 {
350 pix->GetHist()->CutAllEdges();
351 }
352
353 return;
354}
355
356
357// The types are as follows:
358//
359// 0: Fitted Charge
360// 1: Error of fitted Charge
361// 2: Sigma of fitted Charge
362// 3: Error of Sigma of fitted Charge
363// 4: Returned probability of Gauss fit to Charge distribution
364// 5: Mean arrival time
365// 6: Sigma of the arrival time
366// 7: Chi-square of the Gauss fit to the arrival times
367// 8: Pedestal
368// 9: Pedestal RMS
369// 10: Reduced Sigma Square
370// 11: Number of Photo-electrons after the F-Factor method
371// 12: Error on the Number of Photo-electrons after the F-Factor method
372// 13: Mean conversion factor after the F-Factor method
373// 14: Error on the conversion factor after the F-Factor method
374// 15: Number of Photons after the Blind Pixel method
375// 16: Mean conversion factor after the Blind Pixel method
376//
377Bool_t MCalibrationCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
378{
379
380 if (idx > GetSize())
381 return kFALSE;
382
383 Float_t arearatio = cam.GetPixRatio(idx);
384
385 if (arearatio == 0)
386 return kFALSE;
387
388 switch (type)
389 {
390 case 0:
391 if ((*this)[idx].IsExcluded())
392 return kFALSE;
393 val = (*this)[idx].GetCharge();
394 break;
395 case 1:
396 if ((*this)[idx].IsExcluded())
397 return kFALSE;
398 val = (*this)[idx].GetErrCharge();
399 break;
400 case 2:
401 if ((*this)[idx].IsExcluded())
402 return kFALSE;
403 val = (*this)[idx].GetSigmaCharge();
404 break;
405 case 3:
406 if ((*this)[idx].IsExcluded())
407 return kFALSE;
408 val = (*this)[idx].GetErrSigmaCharge();
409 break;
410 case 4:
411 if ((*this)[idx].IsExcluded())
412 return kFALSE;
413 val = (*this)[idx].GetChargeProb();
414 break;
415 case 5:
416 if ((*this)[idx].IsExcluded())
417 return kFALSE;
418 val = (*this)[idx].GetRSigmaCharge();
419 break;
420 case 6:
421 if ((*this)[idx].IsExcluded())
422 return kFALSE;
423 val = (*this)[idx].GetErrRSigmaCharge();
424 break;
425 case 7:
426 if ((*this)[idx].IsExcluded())
427 return kFALSE;
428 val = (*this)[idx].GetRSigmaCharge() / (*this)[idx].GetCharge();
429 break;
430 case 8:
431 if ((*this)[idx].IsExcluded())
432 return kFALSE;
433 // relative error RsigmaCharge square
434 val = (*this)[idx].GetErrRSigmaCharge()* (*this)[idx].GetErrRSigmaCharge()
435 / ((*this)[idx].GetRSigmaCharge() * (*this)[idx].GetRSigmaCharge() );
436 // relative error Charge square
437 val += (*this)[idx].GetErrCharge() * (*this)[idx].GetErrCharge()
438 / ((*this)[idx].GetCharge() * (*this)[idx].GetCharge() );
439 // calculate relative error out of squares
440 val = TMath::Sqrt(val) ;
441 // multiply with value to get absolute error
442 val *= (*this)[idx].GetRSigmaCharge() / (*this)[idx].GetCharge();
443 break;
444 case 9:
445 if ((*this)[idx].IsExcluded())
446 return kFALSE;
447 val = (*this)[idx].GetPheFFactorMethod();
448 break;
449 case 10:
450 if ((*this)[idx].IsExcluded())
451 return kFALSE;
452 val = (*this)[idx].GetPheFFactorMethodError();
453 break;
454 case 11:
455 if ((*this)[idx].IsExcluded())
456 return kFALSE;
457 val = (*this)[idx].GetMeanConversionFFactorMethod();
458 break;
459 case 12:
460 if ((*this)[idx].IsExcluded())
461 return kFALSE;
462 val = (*this)[idx].GetErrorConversionFFactorMethod();
463 break;
464 case 13:
465 if ((*this)[idx].IsExcluded())
466 return kFALSE;
467 val = (*this)[idx].GetTotalFFactorFFactorMethod();
468 break;
469 case 14:
470 if ((*this)[idx].IsExcluded())
471 return kFALSE;
472 val = (*this)[idx].GetTotalFFactorErrorFFactorMethod();
473 break;
474 case 15:
475 if ((*this)[idx].IsExcluded())
476 return kFALSE;
477 if (arearatio == 1.)
478 val = GetMeanPhotInsidePlexiglass();
479 else
480 val = GetMeanPhotInsidePlexiglass()*gkCalibrationOutervsInnerPixelArea;
481 break;
482 case 16:
483 if ((*this)[idx].IsExcluded())
484 return kFALSE;
485 if (arearatio == 1.)
486 val = (double)fMeanPhotInsidePlexiglass;
487 else
488 val = (double)fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelArea;
489 break;
490 case 17:
491 if ((*this)[idx].IsExcluded())
492 return kFALSE;
493 if (arearatio == 1.)
494 val = (*this)[idx].GetMeanConversionBlindPixelMethod();
495 else
496 val = (*this)[idx].GetMeanConversionBlindPixelMethod()*gkCalibrationOutervsInnerPixelArea;
497 break;
498 case 18:
499 if ((*this)[idx].IsExcluded())
500 return kFALSE;
501 if (arearatio == 1.)
502 val = (*this)[idx].GetErrorConversionBlindPixelMethod();
503 else
504 {
505 val = (*this)[idx].GetErrorConversionBlindPixelMethod()*(*this)[idx].GetErrorConversionBlindPixelMethod()
506 * gkCalibrationOutervsInnerPixelArea * gkCalibrationOutervsInnerPixelArea;
507 val += gkCalibrationOutervsInnerPixelAreaError * gkCalibrationOutervsInnerPixelAreaError
508 * (*this)[idx].GetMeanConversionBlindPixelMethod() *(*this)[idx].GetMeanConversionBlindPixelMethod();
509 val = TMath::Sqrt(val);
510 }
511 break;
512 case 19:
513 if ((*this)[idx].IsExcluded())
514 return kFALSE;
515 val = (*this)[idx].GetTotalFFactorBlindPixelMethod();
516 break;
517 case 20:
518 if ((*this)[idx].IsExcluded())
519 return kFALSE;
520 val = (*this)[idx].GetTotalFFactorErrorBlindPixelMethod();
521 break;
522 case 21:
523 if ((*this)[idx].IsExcluded())
524 return kFALSE;
525 if (arearatio == 1.)
526 val = GetMeanPhotOutsidePlexiglass();
527 else
528 val = GetMeanPhotOutsidePlexiglass()*gkCalibrationOutervsInnerPixelArea;
529 break;
530 case 22:
531 if ((*this)[idx].IsExcluded())
532 return kFALSE;
533 if (arearatio == 1.)
534 val = (double)fMeanPhotOutsidePlexiglass;
535 else
536 val = (double)fMeanPhotOutsidePlexiglass*gkCalibrationOutervsInnerPixelArea;
537 break;
538 case 23:
539 if ((*this)[idx].IsExcluded())
540 return kFALSE;
541 if (arearatio == 1.)
542 val = (*this)[idx].GetMeanConversionPINDiodeMethod();
543 else
544 val = (*this)[idx].GetMeanConversionPINDiodeMethod()*gkCalibrationOutervsInnerPixelArea;
545 break;
546 case 24:
547 if ((*this)[idx].IsExcluded())
548 return kFALSE;
549 if (arearatio == 1.)
550 val = (*this)[idx].GetErrorConversionPINDiodeMethod();
551 else
552 {
553 val = (*this)[idx].GetErrorConversionPINDiodeMethod()*(*this)[idx].GetErrorConversionPINDiodeMethod()
554 * gkCalibrationOutervsInnerPixelArea * gkCalibrationOutervsInnerPixelArea;
555 val += gkCalibrationOutervsInnerPixelAreaError * gkCalibrationOutervsInnerPixelAreaError
556 * (*this)[idx].GetMeanConversionPINDiodeMethod() *(*this)[idx].GetMeanConversionPINDiodeMethod();
557 val = TMath::Sqrt(val);
558 }
559 break;
560 case 25:
561 if ((*this)[idx].IsExcluded())
562 return kFALSE;
563 val = (*this)[idx].GetTotalFFactorBlindPixelMethod();
564 break;
565 case 26:
566 if ((*this)[idx].IsExcluded())
567 return kFALSE;
568 val = (*this)[idx].GetTotalFFactorErrorBlindPixelMethod();
569 break;
570 case 27:
571 if ((*this)[idx].IsExcluded())
572 val = 1.;
573 else
574 return kFALSE;
575 break;
576 case 28:
577 if ((*this)[idx].IsExcluded())
578 return kFALSE;
579 if (!(*this)[idx].IsFitted())
580 val = 1;
581 else
582 return kFALSE;
583 break;
584 case 29:
585 if ((*this)[idx].IsExcluded())
586 return kFALSE;
587 if (!(*this)[idx].IsFitted())
588 return kFALSE;
589 if (!(*this)[idx].IsChargeValid())
590 val = 1;
591 else
592 return kFALSE;
593 break;
594 case 30:
595 if ((*this)[idx].IsExcluded())
596 return kFALSE;
597 if ((*this)[idx].IsOscillating())
598 val = 1;
599 else
600 return kFALSE;
601 break;
602 case 31:
603 if ((*this)[idx].IsExcluded())
604 return kFALSE;
605 if ((*this)[idx].IsHiGainSaturation())
606 val = 1;
607 else
608 return kFALSE;
609 break;
610 case 32:
611 if ((*this)[idx].IsExcluded())
612 return kFALSE;
613 if ((*this)[idx].IsFFactorMethodValid())
614 val = 1;
615 else
616 return kFALSE;
617 break;
618 case 33:
619 if ((*this)[idx].IsExcluded())
620 return kFALSE;
621 if ((*this)[idx].IsBlindPixelMethodValid())
622 val = 1;
623 else
624 return kFALSE;
625 break;
626 case 34:
627 if ((*this)[idx].IsExcluded())
628 return kFALSE;
629 if ((*this)[idx].IsPINDiodeMethodValid())
630 val = 1;
631 else
632 return kFALSE;
633 break;
634 case 35:
635 if ((*this)[idx].IsExcluded())
636 return kFALSE;
637 val = (*this)[idx].GetPed();
638 break;
639 case 36:
640 if ((*this)[idx].IsExcluded())
641 return kFALSE;
642 val = 1.;
643 // val = (*this)[idx].GetPedError();
644 break;
645 case 37:
646 if ((*this)[idx].IsExcluded())
647 return kFALSE;
648 val = (*this)[idx].GetPedRms();
649 break;
650 case 38:
651 if ((*this)[idx].IsExcluded())
652 return kFALSE;
653 val = 1.;
654 // val = (*this)[idx].GetPedRmsError();
655 break;
656 case 39:
657 if ((*this)[idx].IsExcluded())
658 return kFALSE;
659 val = (*this)[idx].GetAbsTimeMean();
660 break;
661 case 40:
662 if ((*this)[idx].IsExcluded())
663 return kFALSE;
664 val = (*this)[idx].GetAbsTimeMeanErr();
665 break;
666 case 41:
667 if ((*this)[idx].IsExcluded())
668 return kFALSE;
669 val = (*this)[idx].GetAbsTimeRms();
670 break;
671 case 42:
672 if ((*this)[idx].IsExcluded())
673 return kFALSE;
674 val = (*this)[idx].GetAbsTimeMeanErr()/TMath::Sqrt(2.);
675 break;
676 default:
677 return kFALSE;
678 }
679 return val!=-1.;
680}
681
682// --------------------------------------------------------------------------
683//
684// What MHCamera needs in order to draw an individual pixel in the camera
685//
686void MCalibrationCam::DrawPixelContent(Int_t idx) const
687{
688 (*this)[idx].Draw();
689}
690
691
692// --------------------------------------------------------------------------
693//
694//
695//
696Bool_t MCalibrationCam::CalcNumPhotInsidePlexiglass()
697{
698
699 if (!fBlindPixel->IsFitOK())
700 return kFALSE;
701
702 const Float_t mean = fBlindPixel->GetLambda();
703 const Float_t merr = fBlindPixel->GetErrLambda();
704
705 //
706 // Start calculation of number of photons
707 //
708 // The blind pixel has exactly one cm^2 area (with negligible error),
709 // thus we omi division by 1.
710 //
711 fMeanPhotInsidePlexiglass = mean * gkCalibrationInnerPixelArea;
712
713 // Start calculation of number of photons relative Variance (!!)
714 fMeanPhotErrInsidePlexiglass = merr*merr/mean/mean;
715 fMeanPhotErrInsidePlexiglass += gkCalibrationInnerPixelAreaError*gkCalibrationInnerPixelAreaError
716 / gkCalibrationInnerPixelArea/gkCalibrationInnerPixelArea;
717
718 switch (fColor)
719 {
720 case kECGreen:
721 fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQEGreen;
722 fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQEGreenError*gkCalibrationBlindPixelQEGreenError
723 / gkCalibrationBlindPixelQEGreen / gkCalibrationBlindPixelQEGreen;
724
725 fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttGreen); // correct for absorption
726 // attenuation has negligible error
727 break;
728 case kECBlue:
729 fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQEBlue;
730 fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQEBlueError*gkCalibrationBlindPixelQEBlueError
731 / gkCalibrationBlindPixelQEBlue / gkCalibrationBlindPixelQEBlue;
732
733 fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttBlue); // correct for absorption
734 // attenuation has negligible error
735 break;
736 case kECUV:
737 fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQEUV;
738 fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQEUVError*gkCalibrationBlindPixelQEUVError
739 / gkCalibrationBlindPixelQEUV / gkCalibrationBlindPixelQEUV;
740
741 fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttUV); // correct for absorption
742 // attenuation has negligible error
743 break;
744 case kECCT1:
745 default:
746 fMeanPhotInsidePlexiglass /= gkCalibrationBlindPixelQECT1;
747 fMeanPhotErrInsidePlexiglass += gkCalibrationBlindPixelQECT1Error*gkCalibrationBlindPixelQECT1Error
748 / gkCalibrationBlindPixelQECT1 / gkCalibrationBlindPixelQECT1;
749
750 fMeanPhotInsidePlexiglass *= TMath::Power(10,gkCalibrationBlindPixelAttCT1); // correct for absorption
751 // attenuation has negligible error
752 break;
753 }
754
755 *fLog << inf << endl;
756 *fLog << inf << " Mean number of Photons for an Inner Pixel (inside Plexiglass): "
757 << fMeanPhotInsidePlexiglass << endl;
758
759 if (fMeanPhotInsidePlexiglass > 0.)
760 SETBIT(fFlags,kNumPhotInsidePlexiglassAvailable);
761 else
762 {
763 CLRBIT(fFlags,kNumPhotInsidePlexiglassAvailable);
764 return kFALSE;
765 }
766
767 if (fMeanPhotErrInsidePlexiglass < 0.)
768 {
769 *fLog << warn << " Relative Variance on number of Photons for an Inner Pixel (inside Plexiglass): "
770 << fMeanPhotErrInsidePlexiglass << endl;
771 CLRBIT(fFlags,kNumPhotInsidePlexiglassAvailable);
772 return kFALSE;
773 }
774
775 // Finish calculation of errors -> convert from relative variance to absolute error
776 fMeanPhotErrInsidePlexiglass = TMath::Sqrt(fMeanPhotErrInsidePlexiglass);
777 fMeanPhotErrInsidePlexiglass *= fMeanPhotInsidePlexiglass;
778
779 *fLog << inf << " Error on number of Photons for an Inner Pixel (inside Plexiglass): "
780 << fMeanPhotErrInsidePlexiglass << endl;
781 *fLog << inf << endl;
782
783
784
785 TIter Next(fPixels);
786 MCalibrationPix *pix;
787 while ((pix=(MCalibrationPix*)Next()))
788 {
789
790 if(pix->IsChargeValid())
791 {
792
793 const Float_t charge = pix->GetCharge();
794 const Float_t ratio = fGeomCam->GetPixRatio(pix->GetPixId());
795 const Float_t chargeerr = pix->GetErrCharge();
796
797 Float_t nphot;
798 Float_t nphoterr;
799
800 if (ratio == 1.)
801 {
802 nphot = fMeanPhotInsidePlexiglass;
803 nphoterr = fMeanPhotErrInsidePlexiglass;
804 }
805 else
806 {
807 nphot = fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelArea;
808 nphoterr = fMeanPhotErrInsidePlexiglass*gkCalibrationOutervsInnerPixelArea
809 *fMeanPhotErrInsidePlexiglass*gkCalibrationOutervsInnerPixelArea;
810 nphoterr += fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelAreaError
811 *fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelAreaError;
812 nphoterr = TMath::Sqrt(nphoterr);
813 }
814
815
816
817 const Float_t conversion = nphot/charge;
818 Float_t conversionerr;
819
820 conversionerr = nphoterr/charge
821 * nphoterr/charge ;
822 conversionerr += chargeerr/charge
823 * chargeerr/charge
824 * conversion*conversion;
825 conversionerr = TMath::Sqrt(conversionerr);
826
827 const Float_t conversionsigma = 0.;
828
829 pix->SetConversionBlindPixelMethod(conversion, conversionerr, conversionsigma);
830
831 if (conversionerr/conversion < 0.1)
832 pix->SetBlindPixelMethodValid();
833 }
834 }
835 return kTRUE;
836}
837
838
839Bool_t MCalibrationCam::CalcNumPhotOutsidePlexiglass()
840{
841
842 if (!fPINDiode->IsChargeFitValid())
843 return kFALSE;
844
845 const Float_t mean = fPINDiode->GetCharge();
846 const Float_t merr = fPINDiode->GetErrCharge();
847
848 // Start calculation of number of photons
849 fMeanPhotOutsidePlexiglass = mean * gkCalibrationInnerPixelvsPINDiodeArea;
850
851 // Start calculation of number of photons relative Variance (!!)
852 fMeanPhotErrOutsidePlexiglass = merr*merr/mean/mean;
853 fMeanPhotErrOutsidePlexiglass += gkCalibrationInnerPixelvsPINDiodeAreaError*gkCalibrationInnerPixelvsPINDiodeAreaError
854 / gkCalibrationInnerPixelvsPINDiodeArea/gkCalibrationInnerPixelvsPINDiodeArea;
855
856 switch (fColor)
857 {
858 case kECGreen:
859 fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQEGreen;
860 fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQEGreenError*gkCalibrationPINDiodeQEGreenError
861 / gkCalibrationPINDiodeQEGreen/gkCalibrationPINDiodeQEGreen;
862 break;
863 case kECBlue:
864 fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQEBlue;
865 fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQEBlueError*gkCalibrationPINDiodeQEBlueError
866 / gkCalibrationPINDiodeQEBlue/gkCalibrationPINDiodeQEBlue;
867 break;
868 case kECUV:
869 fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQEUV;
870 fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQEUVError*gkCalibrationPINDiodeQEUVError
871 / gkCalibrationPINDiodeQEUV/gkCalibrationPINDiodeQEUV;
872 break;
873 case kECCT1:
874 default:
875 fMeanPhotOutsidePlexiglass /= gkCalibrationPINDiodeQECT1;
876 fMeanPhotErrOutsidePlexiglass += gkCalibrationPINDiodeQECT1Error*gkCalibrationPINDiodeQECT1Error
877 / gkCalibrationPINDiodeQECT1/gkCalibrationPINDiodeQECT1;
878 break;
879 }
880
881
882 *fLog << inf << endl;
883 *fLog << inf << " Mean number of Photons for an Inner Pixel (outside Plexiglass): "
884 << fMeanPhotOutsidePlexiglass << endl;
885
886 if (fMeanPhotOutsidePlexiglass > 0.)
887 SETBIT(fFlags,kNumPhotOutsidePlexiglassAvailable);
888 else
889 {
890 CLRBIT(fFlags,kNumPhotOutsidePlexiglassAvailable);
891 return kFALSE;
892 }
893
894 if (fMeanPhotErrOutsidePlexiglass < 0.)
895 {
896 *fLog << warn << "Relative Variance on number of Photons for an Inner Pixel (outside Plexiglass): "
897 << fMeanPhotErrOutsidePlexiglass << endl;
898 CLRBIT(fFlags,kNumPhotOutsidePlexiglassAvailable);
899 return kFALSE;
900 }
901
902 // Finish calculation of errors -> convert from relative variance to absolute error
903 fMeanPhotErrOutsidePlexiglass = TMath::Sqrt(fMeanPhotErrOutsidePlexiglass);
904 fMeanPhotErrOutsidePlexiglass *= fMeanPhotOutsidePlexiglass;
905
906 *fLog << inf << " Error on number of Photons for an Inner Pixel (outside Plexiglass): "
907 << fMeanPhotErrOutsidePlexiglass << endl;
908 *fLog << inf << endl;
909
910 TIter Next(fPixels);
911 MCalibrationPix *pix;
912 while ((pix=(MCalibrationPix*)Next()))
913 {
914
915 if (pix->IsChargeValid())
916 {
917
918 const Float_t charge = pix->GetCharge();
919 const Float_t ratio = fGeomCam->GetPixRatio(pix->GetPixId());
920 const Float_t chargeerr = pix->GetErrCharge();
921
922 Float_t nphot;
923 Float_t nphoterr;
924
925 if (ratio == 1.)
926 {
927 nphot = fMeanPhotInsidePlexiglass;
928 nphoterr = fMeanPhotErrInsidePlexiglass;
929 }
930 else
931 {
932 nphot = fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelArea;
933 nphoterr = fMeanPhotErrInsidePlexiglass*gkCalibrationOutervsInnerPixelArea
934 *fMeanPhotErrInsidePlexiglass*gkCalibrationOutervsInnerPixelArea;
935 nphoterr += fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelAreaError
936 *fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelAreaError;
937 nphoterr = TMath::Sqrt(nphoterr);
938 }
939
940 const Float_t conversion = nphot/charge;
941 Float_t conversionerr;
942
943 conversionerr = nphoterr/charge
944 * nphoterr/charge ;
945 conversionerr += chargeerr/charge
946 * chargeerr/charge
947 * conversion*conversion;
948 if (conversionerr > 0.)
949 conversionerr = TMath::Sqrt(conversionerr);
950
951 const Float_t conversionsigma = 0.;
952
953 pix->SetConversionBlindPixelMethod(conversion, conversionerr, conversionsigma);
954
955 if (conversionerr/conversion < 0.1)
956 pix->SetBlindPixelMethodValid();
957
958 }
959 }
960 return kTRUE;
961}
962
963
964
965Bool_t MCalibrationCam::GetConversionFactorBlindPixel(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
966{
967
968 if (ipx < 0 || !IsPixelFitted(ipx))
969 return kFALSE;
970
971 if (!IsNumPhotInsidePlexiglassAvailable())
972 if (!CalcNumPhotInsidePlexiglass())
973 return kFALSE;
974
975 mean = (*this)[ipx].GetMeanConversionBlindPixelMethod();
976 err = (*this)[ipx].GetErrorConversionBlindPixelMethod();
977 sigma = (*this)[ipx].GetSigmaConversionBlindPixelMethod();
978
979 return kTRUE;
980}
981
982
983Bool_t MCalibrationCam::GetConversionFactorFFactor(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
984{
985
986 if (ipx < 0 || !IsPixelFitted(ipx))
987 return kFALSE;
988
989 Float_t conv = (*this)[ipx].GetMeanConversionFFactorMethod();
990
991 if (conv < 0.)
992 return kFALSE;
993
994 mean = conv;
995 err = (*this)[ipx].GetErrorConversionFFactorMethod();
996 sigma = (*this)[ipx].GetSigmaConversionFFactorMethod();
997
998 return kTRUE;
999}
1000
1001
1002//-----------------------------------------------------------------------------------
1003//
1004// Calculates the conversion factor between the integral of FADCs slices
1005// (as defined in the signal extractor MExtractSignal.cc)
1006// and the number of photons reaching the plexiglass for one Inner Pixel
1007//
1008// FIXME: The PINDiode is still not working and so is the code
1009//
1010Bool_t MCalibrationCam::GetConversionFactorPINDiode(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
1011{
1012
1013 if (ipx < 0 || !IsPixelFitted(ipx))
1014 return kFALSE;
1015
1016 if (!IsNumPhotOutsidePlexiglassAvailable())
1017 if (!CalcNumPhotOutsidePlexiglass())
1018 return kFALSE;
1019
1020 mean = (*this)[ipx].GetMeanConversionPINDiodeMethod();
1021 err = (*this)[ipx].GetErrorConversionPINDiodeMethod();
1022 sigma = (*this)[ipx].GetSigmaConversionPINDiodeMethod();
1023
1024 return kFALSE;
1025
1026}
1027
1028//-----------------------------------------------------------------------------------
1029//
1030// Calculates the best combination of the three used methods possible
1031// between the integral of FADCs slices
1032// (as defined in the signal extractor MExtractSignal.cc)
1033// and the number of photons reaching one Inner Pixel.
1034// The procedure is not yet defined.
1035//
1036// FIXME: The PINDiode is still not working and so is the code
1037//
1038Bool_t MCalibrationCam::GetConversionFactorCombined(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
1039{
1040
1041 if (ipx < 0 || !IsPixelFitted(ipx))
1042 return kFALSE;
1043
1044 return kFALSE;
1045
1046}
1047
1048
1049void MCalibrationCam::DrawHiLoFits()
1050{
1051
1052 if (!fOffsets)
1053 fOffsets = new TH1D("pp","Offsets of the HiGain LoGain Fit",100,-600.,400.);
1054 if (!fSlopes)
1055 fSlopes = new TH1D("mm","Slopes of the HiGain LoGain Fit",100,-2.,2.);
1056 if (!fOffvsSlope)
1057 fOffvsSlope = new TH2D("aa","Slopes vs Offsets of the HiGain LoGain Fit",100,-600.,400.,100,-2.,2.);
1058
1059 TIter Next(fPixels);
1060 MCalibrationPix *pix;
1061 MHCalibrationPixel *hist;
1062 while ((pix=(MCalibrationPix*)Next()))
1063 {
1064 hist = pix->GetHist();
1065 hist->FitHiGainvsLoGain();
1066 fOffsets->Fill(hist->GetOffset(),1.);
1067 fSlopes->Fill(hist->GetSlope(),1.);
1068 fOffvsSlope->Fill(hist->GetOffset(),hist->GetSlope(),1.);
1069 }
1070
1071 TCanvas *c1 = new TCanvas();
1072
1073 c1->Divide(1,3);
1074 c1->cd(1);
1075 fOffsets->Draw();
1076 gPad->Modified();
1077 gPad->Update();
1078
1079 c1->cd(2);
1080 fSlopes->Draw();
1081 gPad->Modified();
1082 gPad->Update();
1083
1084 c1->cd(3);
1085 fOffvsSlope->Draw("col1");
1086 gPad->Modified();
1087 gPad->Update();
1088}
1089
Note: See TracBrowser for help on using the repository browser.