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

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