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

Last change on this file since 2943 was 2943, checked in by gaug, 21 years ago
*** empty log message ***
File size: 18.7 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// 4)
38//
39/////////////////////////////////////////////////////////////////////////////
40#include "MCalibrationCam.h"
41
42#include <TH2.h>
43#include <TCanvas.h>
44#include <TClonesArray.h>
45
46#include "MLog.h"
47#include "MLogManip.h"
48
49#include "MGeomCam.h"
50
51#include "MCalibrationPix.h"
52#include "MCalibrationConfig.h"
53#include "MCalibrationBlindPix.h"
54#include "MCalibrationPINDiode.h"
55
56#include "MHCalibrationPixel.h"
57
58ClassImp(MCalibrationCam);
59
60using namespace std;
61// --------------------------------------------------------------------------
62//
63// Default constructor.
64//
65// Creates a TClonesArray of MCalibrationPix containers, initialized to 1 entry
66// Later, a call to MCalibrationCam::InitSize(Int_t size) has to be performed
67//
68// Creates an MCalibrationBlindPix container
69// Creates an MCalibrationPINDiode container
70//
71MCalibrationCam::MCalibrationCam(const char *name, const char *title)
72 : fOffsets(NULL),
73 fSlopes(NULL),
74 fOffvsSlope(NULL)
75{
76 fName = name ? name : "MCalibrationCam";
77 fTitle = title ? title : "Storage container for the Calibration Information in the camera";
78
79 fPixels = new TClonesArray("MCalibrationPix",1);
80 fBlindPixel = new MCalibrationBlindPix();
81 fPINDiode = new MCalibrationPINDiode();
82
83 Clear();
84}
85
86// --------------------------------------------------------------------------
87//
88// Delete the TClonesArray of MCalibrationPix containers
89// Delete the MCalibrationPINDiode and the MCalibrationBlindPix
90//
91// Delete the histograms if they exist
92//
93MCalibrationCam::~MCalibrationCam()
94{
95
96 //
97 // delete fPixels should delete all Objects stored inside
98 //
99 delete fPixels;
100 delete fBlindPixel;
101 delete fPINDiode;
102
103 if (fOffsets)
104 delete fOffsets;
105 if (fSlopes)
106 delete fSlopes;
107 if (fOffvsSlope)
108 delete fOffvsSlope;
109
110}
111
112// -------------------------------------------------------------------
113//
114// This function simply allocates memory via the ROOT command:
115// (TObject**) TStorage::ReAlloc(fCont, newSize * sizeof(TObject*),
116// fSize * sizeof(TObject*));
117// newSize corresponds to size in our case
118// fSize is the old size (in most cases: 1)
119//
120void MCalibrationCam::InitSize(const UInt_t i)
121{
122
123 //
124 // check if we have already initialized to size
125 //
126 if (CheckBounds(i))
127 return;
128
129 fPixels->ExpandCreate(i);
130
131}
132
133// --------------------------------------------------------------------------
134//
135// This function returns the current size of the TClonesArray
136// independently if the MCalibrationPix is filled with values or not.
137//
138// It is the size of the array fPixels.
139//
140Int_t MCalibrationCam::GetSize() const
141{
142 return fPixels->GetEntriesFast();
143}
144
145// --------------------------------------------------------------------------
146//
147// Check if position i is inside the current bounds of the TClonesArray
148//
149Bool_t MCalibrationCam::CheckBounds(Int_t i) const
150{
151 return i < GetSize();
152}
153
154
155// --------------------------------------------------------------------------
156//
157// Get i-th pixel (pixel number)
158//
159MCalibrationPix &MCalibrationCam::operator[](Int_t i)
160{
161 return *static_cast<MCalibrationPix*>(fPixels->UncheckedAt(i));
162}
163
164// --------------------------------------------------------------------------
165//
166// Get i-th pixel (pixel number)
167//
168MCalibrationPix &MCalibrationCam::operator[](Int_t i) const
169{
170 return *static_cast<MCalibrationPix*>(fPixels->UncheckedAt(i));
171}
172
173
174// --------------------------------------
175//
176void MCalibrationCam::Clear(Option_t *o)
177{
178
179 fPixels->ForEach(TObject, Clear)();
180 fBlindPixel->Clear();
181 fPINDiode->Clear();
182
183 fNumPhotInsidePlexiglassAvailable = kFALSE;
184 fMeanPhotInsidePlexiglass = -1.;
185 fMeanPhotErrInsidePlexiglass = -1.;
186 fNumPhotOutsidePlexiglassAvailable = kFALSE;
187 fMeanPhotOutsidePlexiglass = -1.;
188 fMeanPhotErrOutsidePlexiglass = -1.;
189
190 fNumExcludedPixels = 0;
191
192 return;
193}
194
195void MCalibrationCam::SetBlindPixelMethodValid(const Bool_t b)
196{
197
198 if (b)
199 SETBIT(fFlags, kBlindPixelMethodValid);
200 else
201 CLRBIT(fFlags, kBlindPixelMethodValid);
202
203}
204
205void MCalibrationCam::SetPINDiodeMethodValid(const Bool_t b)
206{
207
208 if (b)
209 SETBIT(fFlags, kPINDiodeMethodValid);
210 else
211 CLRBIT(fFlags, kPINDiodeMethodValid);
212
213
214}
215
216Bool_t MCalibrationCam::IsBlindPixelMethodValid() const
217{
218 return TESTBIT(fFlags,kBlindPixelMethodValid);
219}
220
221Bool_t MCalibrationCam::IsPINDiodeMethodValid() const
222{
223 return TESTBIT(fFlags,kPINDiodeMethodValid);
224}
225
226
227// --------------------------------------------------------------------------
228//
229// Print first the well fitted pixels
230// and then the ones which are not FitValid
231//
232void MCalibrationCam::Print(Option_t *o) const
233{
234
235 *fLog << all << GetDescriptor() << ":" << endl;
236 int id = 0;
237
238 *fLog << all << "Succesfully calibrated pixels:" << endl;
239 *fLog << all << endl;
240
241 TIter Next(fPixels);
242 MCalibrationPix *pix;
243 while ((pix=(MCalibrationPix*)Next()))
244 {
245
246 if (pix->IsFitValid() && !pix->IsExcluded())
247 {
248
249 Float_t rsigma = pix->GetRSigmaSquare();
250 if (rsigma > 0.)
251 rsigma = TMath::Sqrt(rsigma);
252
253 *fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- "
254 << pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- "
255 << pix->GetSigmaCharge() << " Reduced Sigma: " << rsigma
256 << " Nr Phe's: " << pix->GetPheFFactorMethod() << endl;
257 id++;
258 }
259 }
260
261 *fLog << all << id << " succesful pixels :-))" << endl;
262 id = 0;
263
264 *fLog << all << endl;
265 *fLog << all << "Pixels with errors:" << endl;
266 *fLog << all << endl;
267
268 TIter Next2(fPixels);
269 while ((pix=(MCalibrationPix*)Next2()))
270 {
271
272 if (!pix->IsFitValid() && !pix->IsExcluded())
273 {
274
275 Float_t rsigma = pix->GetRSigmaSquare();
276 if (rsigma > 0.)
277 rsigma = TMath::Sqrt(rsigma);
278
279 *fLog << all << pix->GetPixId() << " Pedestals: " << pix->GetPed() << " +- "
280 << pix->GetPedRms() << " Reduced Charge: " << pix->GetCharge() << " +- "
281 << pix->GetSigmaCharge() << " Reduced Sigma: " << rsigma << endl;
282 id++;
283 }
284 }
285 *fLog << all << id << " pixels with errors :-((" << endl;
286
287 *fLog << all << endl;
288 *fLog << all << "Excluded pixels:" << endl;
289 *fLog << all << endl;
290
291 TIter Next3(fPixels);
292 while ((pix=(MCalibrationPix*)Next3()))
293 if (pix->IsExcluded())
294 *fLog << all << pix->GetPixId() << endl;
295
296 *fLog << all << fNumExcludedPixels << " excluded pixels " << endl;
297}
298
299// --------------------------------------------------------------------------
300//
301// Return true if pixel is inside bounds of the TClonesArray fPixels
302//
303Bool_t MCalibrationCam::IsPixelUsed(Int_t idx) const
304{
305 if (!CheckBounds(idx))
306 return kFALSE;
307
308 return kTRUE;
309}
310
311// --------------------------------------------------------------------------
312//
313// Return true if pixel has already been fitted once (independent of the result)
314//
315Bool_t MCalibrationCam::IsPixelFitted(Int_t idx) const
316{
317
318 if (!CheckBounds(idx))
319 return kFALSE;
320
321 return (*this)[idx].IsFitted();
322}
323
324// --------------------------------------------------------------------------
325//
326// Sets the user ranges of all histograms such that
327// empty bins at the edges are not used. Additionally, it rebins the
328// histograms such that in total, 50 bins are used.
329//
330void MCalibrationCam::CutEdges()
331{
332
333 fBlindPixel->GetHist()->CutAllEdges();
334 fPINDiode->GetHist()->CutAllEdges();
335
336 TIter Next(fPixels);
337 MCalibrationPix *pix;
338 while ((pix=(MCalibrationPix*)Next()))
339 {
340 pix->GetHist()->CutAllEdges();
341 }
342
343 return;
344}
345
346
347// The types are as follows:
348//
349// 0: Fitted Charge
350// 1: Error of fitted Charge
351// 2: Sigma of fitted Charge
352// 3: Error of Sigma of fitted Charge
353// 4: Returned probability of Gauss fit to Charge distribution
354// 5: Mean arrival time
355// 6: Sigma of the arrival time
356// 7: Chi-square of the Gauss fit to the arrival times
357// 8: Pedestal
358// 9: Pedestal RMS
359// 10: Reduced Sigma Square
360// 11: Number of Photo-electrons after the F-Factor method
361// 12: Error on the Number of Photo-electrons after the F-Factor method
362// 13: Mean conversion factor after the F-Factor method
363// 14: Error on the conversion factor after the F-Factor method
364// 15: Number of Photons after the Blind Pixel method
365// 16: Mean conversion factor after the Blind Pixel method
366//
367Bool_t MCalibrationCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
368{
369
370 if (idx > GetSize())
371 return kFALSE;
372
373 if ( (!(*this)[idx].IsFitValid()) || (*this)[idx].IsExcluded())
374 return kFALSE;
375
376 if (idx == gkCalibrationBlindPixelId)
377 return kFALSE;
378
379 switch (type)
380 {
381 case 0:
382 val = (*this)[idx].GetCharge();
383 break;
384 case 1:
385 val = (*this)[idx].GetErrCharge();
386 break;
387 case 2:
388 val = (*this)[idx].GetSigmaCharge();
389 break;
390 case 3:
391 val = (*this)[idx].GetErrSigmaCharge();
392 break;
393 case 4:
394 val = (*this)[idx].GetChargeProb();
395 break;
396 case 5:
397 val = (*this)[idx].GetTime();
398 break;
399 case 6:
400 val = (*this)[idx].GetSigmaTime();
401 break;
402 case 7:
403 val = (*this)[idx].GetTimeChiSquare();
404 break;
405 case 8:
406 val = (*this)[idx].GetPed();
407 break;
408 case 9:
409 val = (*this)[idx].GetPedRms();
410 break;
411 case 10:
412 if ((*this)[idx].GetRSigmaSquare() > 0.)
413 val = TMath::Sqrt((*this)[idx].GetRSigmaSquare());
414 else
415 val = -1.;
416 break;
417 case 11:
418 val = (*this)[idx].GetPheFFactorMethod();
419 break;
420 case 12:
421 val = (*this)[idx].GetPheFFactorMethodError();
422 break;
423 case 13:
424 val = (*this)[idx].GetMeanConversionFFactorMethod();
425 break;
426 case 14:
427 val = (*this)[idx].GetErrorConversionFFactorMethod();
428 break;
429 case 15:
430 if (idx < 397)
431 val = (double)fMeanPhotInsidePlexiglass;
432 else
433 val = (double)fMeanPhotInsidePlexiglass*gkCalibrationOutervsInnerPixelArea;
434 break;
435 case 16:
436 if (idx < 397)
437 val = (*this)[idx].GetMeanConversionBlindPixelMethod();
438 else
439 val = (*this)[idx].GetMeanConversionBlindPixelMethod()*gkCalibrationOutervsInnerPixelArea;
440 break;
441 case 17:
442 if ( (*this)[idx].GetRSigmaSquare() > 0. && (*this)[idx].GetCharge() > 0. )
443 val = TMath::Sqrt((*this)[idx].GetRSigmaSquare()) / (*this)[idx].GetCharge();
444 else
445 val = -1.;
446 break;
447 default:
448 return kFALSE;
449 }
450 return val!=-1.;
451}
452
453// --------------------------------------------------------------------------
454//
455// What MHCamera needs in order to draw an individual pixel in the camera
456//
457void MCalibrationCam::DrawPixelContent(Int_t idx) const
458{
459 (*this)[idx].Draw();
460}
461
462
463// --------------------------------------------------------------------------
464//
465//
466//
467Bool_t MCalibrationCam::CalcNumPhotInsidePlexiglass()
468{
469
470 if (!fBlindPixel->IsFitOK())
471 return kFALSE;
472
473 const Float_t mean = fBlindPixel->GetLambda();
474 const Float_t merr = fBlindPixel->GetErrLambda();
475
476 switch (fColor)
477 {
478 case kECGreen:
479 fMeanPhotInsidePlexiglass = (mean / gkCalibrationBlindPixelQEGreen) // real photons
480 *TMath::Power(10,gkCalibrationBlindPixelAttGreen) // correct for absorption
481 * gkCalibrationInnerPixelArea; // correct for area
482
483
484 break;
485 case kECBlue:
486 fMeanPhotInsidePlexiglass = (mean / gkCalibrationBlindPixelQEBlue )
487 *TMath::Power(10,gkCalibrationBlindPixelAttBlue)
488 * gkCalibrationInnerPixelArea;
489 break;
490 case kECUV:
491 fMeanPhotInsidePlexiglass = (mean / gkCalibrationBlindPixelQEUV )
492 *TMath::Power(10,gkCalibrationBlindPixelAttUV)
493 * gkCalibrationInnerPixelArea;
494 break;
495 case kECCT1:
496 default:
497 fMeanPhotInsidePlexiglass = (mean / gkCalibrationBlindPixelQECT1 )
498 *TMath::Power(10,gkCalibrationBlindPixelAttCT1)
499 * gkCalibrationInnerPixelArea;
500 break;
501 }
502
503 fNumPhotInsidePlexiglassAvailable = kTRUE;
504
505 *fLog << inf << endl;
506 *fLog << inf << "Mean number of Photons for an Inner Pixel (inside Plexiglass): "
507 << fMeanPhotInsidePlexiglass << endl;
508
509 TIter Next(fPixels);
510 MCalibrationPix *pix;
511 while ((pix=(MCalibrationPix*)Next()))
512 {
513 if((pix->GetCharge() > 0.) && (fMeanPhotInsidePlexiglass > 0.))
514 {
515
516 Float_t conversion = fMeanPhotInsidePlexiglass/pix->GetCharge();
517 Float_t conversionerr = 0.;
518 Float_t conversionsigma = 0.;
519 pix->SetConversionBlindPixelMethod(conversion, conversionerr, conversionsigma);
520
521 if (conversionerr/conversion < 0.1)
522 pix->SetBlindPixelMethodValid();
523 }
524 }
525 return kTRUE;
526}
527
528
529Bool_t MCalibrationCam::CalcNumPhotOutsidePlexiglass()
530{
531
532 if (!fPINDiode->IsFitValid())
533 return kFALSE;
534
535 const Float_t mean = fPINDiode->GetCharge();
536 const Float_t merr = fPINDiode->GetErrCharge();
537
538 switch (fColor)
539 {
540 case kECGreen:
541 fMeanPhotOutsidePlexiglass = (mean / gkCalibrationPINDiodeQEGreen) // real photons
542 * gkCalibrationInnerPixelvsPINDiodeArea; // correct for area
543 break;
544 case kECBlue:
545 fMeanPhotOutsidePlexiglass = (mean / gkCalibrationPINDiodeQEBlue )
546 * gkCalibrationInnerPixelvsPINDiodeArea;
547 break;
548 case kECUV:
549 fMeanPhotOutsidePlexiglass = (mean / gkCalibrationPINDiodeQEUV )
550 * gkCalibrationInnerPixelvsPINDiodeArea;
551 break;
552 case kECCT1:
553 default:
554 fMeanPhotOutsidePlexiglass = (mean / gkCalibrationPINDiodeQECT1 )
555 * gkCalibrationInnerPixelvsPINDiodeArea;
556 break;
557 }
558
559 fNumPhotOutsidePlexiglassAvailable = kTRUE;
560
561 *fLog << inf << endl;
562 *fLog << inf << mean << " Mean number of Photons for an Inner Pixel (outside Plexiglass): "
563 << fMeanPhotOutsidePlexiglass << endl;
564 *fLog << inf << endl;
565
566 TIter Next(fPixels);
567 MCalibrationPix *pix;
568 while ((pix=(MCalibrationPix*)Next()))
569 {
570
571 if((pix->GetCharge() > 0.) && (fMeanPhotInsidePlexiglass > 0.))
572 pix->SetConversionPINDiodeMethod(fMeanPhotOutsidePlexiglass/pix->GetCharge(), 0., 0.);
573 }
574 return kTRUE;
575}
576
577
578
579Bool_t MCalibrationCam::GetConversionFactorBlindPixel(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
580{
581
582 if (ipx < 0 || !IsPixelFitted(ipx))
583 return kFALSE;
584
585 if (!fNumPhotInsidePlexiglassAvailable)
586 if (!CalcNumPhotInsidePlexiglass())
587 return kFALSE;
588
589 mean = (*this)[ipx].GetMeanConversionBlindPixelMethod();
590 err = (*this)[ipx].GetErrorConversionBlindPixelMethod();
591 sigma = (*this)[ipx].GetSigmaConversionBlindPixelMethod();
592
593 return kTRUE;
594}
595
596
597Bool_t MCalibrationCam::GetConversionFactorFFactor(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
598{
599
600 if (ipx < 0 || !IsPixelFitted(ipx))
601 return kFALSE;
602
603 Float_t conv = (*this)[ipx].GetMeanConversionFFactorMethod();
604
605 if (conv < 0.)
606 return kFALSE;
607
608 mean = conv;
609 err = (*this)[ipx].GetErrorConversionFFactorMethod();
610 sigma = (*this)[ipx].GetSigmaConversionFFactorMethod();
611
612 return kTRUE;
613}
614
615
616//-----------------------------------------------------------------------------------
617//
618// Calculates the conversion factor between the integral of FADCs slices
619// (as defined in the signal extractor MExtractSignal.cc)
620// and the number of photons reaching the plexiglass for one Inner Pixel
621//
622// FIXME: The PINDiode is still not working and so is the code
623//
624Bool_t MCalibrationCam::GetConversionFactorPINDiode(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
625{
626
627 if (ipx < 0 || !IsPixelFitted(ipx))
628 return kFALSE;
629
630 return kFALSE;
631
632}
633
634//-----------------------------------------------------------------------------------
635//
636// Calculates the best combination of the three used methods possible
637// between the integral of FADCs slices
638// (as defined in the signal extractor MExtractSignal.cc)
639// and the number of photons reaching one Inner Pixel.
640// The procedure is not yet defined.
641//
642// FIXME: The PINDiode is still not working and so is the code
643//
644Bool_t MCalibrationCam::GetConversionFactorCombined(Int_t ipx, Float_t &mean, Float_t &err, Float_t &sigma)
645{
646
647 if (ipx < 0 || !IsPixelFitted(ipx))
648 return kFALSE;
649
650 return kFALSE;
651
652}
653
654
655void MCalibrationCam::DrawHiLoFits()
656{
657
658 if (!fOffsets)
659 fOffsets = new TH1D("pp","Offsets of the HiGain LoGain Fit",100,-600.,400.);
660 if (!fSlopes)
661 fSlopes = new TH1D("mm","Slopes of the HiGain LoGain Fit",100,-2.,2.);
662 if (!fOffvsSlope)
663 fOffvsSlope = new TH2D("aa","Slopes vs Offsets of the HiGain LoGain Fit",100,-600.,400.,100,-2.,2.);
664
665 TIter Next(fPixels);
666 MCalibrationPix *pix;
667 MHCalibrationPixel *hist;
668 while ((pix=(MCalibrationPix*)Next()))
669 {
670 hist = pix->GetHist();
671 hist->FitHiGainvsLoGain();
672 fOffsets->Fill(hist->GetOffset(),1.);
673 fSlopes->Fill(hist->GetSlope(),1.);
674 fOffvsSlope->Fill(hist->GetOffset(),hist->GetSlope(),1.);
675 }
676
677 TCanvas *c1 = new TCanvas();
678
679 c1->Divide(1,3);
680 c1->cd(1);
681 fOffsets->Draw();
682 gPad->Modified();
683 gPad->Update();
684
685 c1->cd(2);
686 fSlopes->Draw();
687 gPad->Modified();
688 gPad->Update();
689
690 c1->cd(3);
691 fOffvsSlope->Draw("col1");
692 gPad->Modified();
693 gPad->Update();
694}
695
Note: See TracBrowser for help on using the repository browser.