source: trunk/MagicSoft/Mars/mcalib/MHCalibrationTestCam.cc@ 4613

Last change on this file since 4613 was 4572, checked in by gaug, 22 years ago
*** empty log message ***
File size: 17.9 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// MHCalibrationTestCam
27//
28// Fills the calibrated signal from an MCerPhotEvt into
29// MHCalibrationTestPix for every:
30//
31// - Pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray
32// or MHCalibrationCam::fHiGainArray, respectively.
33//
34// - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera),
35// stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas and
36// MHCalibrationCam::fAverageHiGainAreas
37//
38// - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera),
39// stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors
40// and MHCalibrationCam::fAverageHiGainSectors
41//
42// The signals are filled into a histogram and an array, in order to perform
43// a Fourier analysis (see MHGausEvents). The signals are moreover averaged on an
44// event-by-event basis and written into the corresponding average pixels.
45//
46// The histograms are fitted to a Gaussian, mean and sigma with its errors
47// and the fit probability are extracted. If none of these values are NaN's and
48// if the probability is bigger than MHGausEvents::fProbLimit (default: 0.5%),
49// the fit is declared valid.
50// Otherwise, the fit is repeated within ranges of the previous mean
51// +- MHGausEvents::fPickupLimit (default: 5) sigma (see MHGausEvents::RepeatFit())
52// In case this does not make the fit valid, the histogram means and RMS's are
53// taken directly (see MHGausEvents::BypassFit()) and the following flags are set:
54// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) and
55// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun )
56//
57// Outliers of more than MHGausEvents::fPickupLimit (default: 5) sigmas
58// from the mean are counted as Pickup events (stored in MHGausEvents::fPickup)
59//
60// The class also fills arrays with the signal vs. event number, creates a fourier
61// spectrum (see MHGausEvents::CreateFourierSpectrum()) and investigates if the
62// projected fourier components follow an exponential distribution.
63// In case that the probability of the exponential fit is less than
64// MHGausEvents::fProbLimit (default: 0.5%), the following flags are set:
65// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating ) and
66// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun )
67//
68// This same procedure is performed for the average pixels.
69//
70// The following results are written into an MCalibrationCam:
71//
72// - MCalibrationPix::SetMean()
73// - MCalibrationPix::SetMeanErr()
74// - MCalibrationPix::SetSigma()
75// - MCalibrationPix::SetSigmaErr()
76// - MCalibrationPix::SetProb()
77// - MCalibrationPix::SetNumPickup()
78//
79// For all averaged areas, the fitted sigma is multiplied with the square root of
80// the number involved pixels in order to be able to compare it to the average of
81// sigmas in the camera.
82//
83/////////////////////////////////////////////////////////////////////////////
84#include "MHCalibrationTestCam.h"
85#include "MHCalibrationTestPix.h"
86
87#include "MLog.h"
88#include "MLogManip.h"
89
90#include "MParList.h"
91
92#include "MCalibrationCam.h"
93#include "MCalibrationPix.h"
94
95#include "MCerPhotEvt.h"
96#include "MCerPhotPix.h"
97
98#include "MGeomCam.h"
99#include "MGeomPix.h"
100
101#include "MBadPixelsCam.h"
102#include "MBadPixelsPix.h"
103
104ClassImp(MHCalibrationTestCam);
105
106using namespace std;
107// --------------------------------------------------------------------------
108//
109// Default Constructor.
110//
111MHCalibrationTestCam::MHCalibrationTestCam(const char *name, const char *title)
112{
113
114 fName = name ? name : "MHCalibrationTestCam";
115 fTitle = title ? title : "Histogram class for testing the calibration";
116
117 SetAverageNbins(5000);
118
119 fNotInterpolateablePixels.Set(0);
120}
121
122// --------------------------------------------------------------------------
123//
124// Gets or creates the pointers to:
125// - MCalibrationTestCam
126//
127// Searches pointer to:
128// - MCerPhotEvt
129//
130// Initializes, if empty to MGeomCam::GetNumAreas() for:
131// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
132//
133// Initializes, if empty to MGeomCam::GetNumSectors() for:
134// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
135//
136// Calls MHCalibrationCam::InitHists() for every entry in:
137// - MHCalibrationCam::fHiGainArray
138// - MHCalibrationCam::fAverageHiGainAreas
139// - MHCalibrationCam::fAverageHiGainSectors
140//
141// Sets Titles and Names for the Histograms
142// - MHCalibrationCam::fAverageHiGainAreas
143// - MHCalibrationCam::fAverageHiGainSectors
144//
145// Sets number of bins to MHCalibrationCam::fAverageNbins for:
146// - MHCalibrationCam::fAverageHiGainAreas
147// - MHCalibrationCam::fAverageHiGainSectors
148//
149Bool_t MHCalibrationTestCam::ReInitHists(MParList *pList)
150{
151
152 MCerPhotEvt *signal = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
153 if (!signal)
154 {
155 *fLog << err << "MCerPhotEvt not found... abort." << endl;
156 return kFALSE;
157 }
158
159 const Int_t npixels = fGeom->GetNumPixels();
160 const Int_t nsectors = fGeom->GetNumSectors();
161 const Int_t nareas = fGeom->GetNumAreas();
162
163 if (fHiGainArray->GetEntries()==0)
164 {
165 fHiGainArray->Expand(npixels);
166 for (Int_t i=0; i<npixels; i++)
167 {
168 (*fHiGainArray)[i] = new MHCalibrationTestPix("Calibrated Events",
169 "Test Calibration Pixel");
170 InitHists((*this)[i],(*fBadPixels)[i],i);
171 }
172 }
173
174 if (fLoGainArray->GetEntries()==0)
175 {
176 fLoGainArray->Expand(npixels);
177 for (Int_t i=0; i<npixels; i++)
178 {
179 (*fLoGainArray)[i] = new MHCalibrationTestPix("Calibrated Events",
180 "Test Calibration Pixel");
181 InitHists((*this)(i),(*fBadPixels)[i],i);
182 }
183 }
184
185
186 if (fAverageHiGainAreas->GetEntries()==0)
187 {
188 fAverageHiGainAreas->Expand(nareas);
189
190 for (Int_t j=0; j<nareas; j++)
191 {
192 (*fAverageHiGainAreas)[j] =
193 new MHCalibrationTestPix("MHCalibrationTestAverageArea",
194 "Average Test Calibrations Area Idx ");
195
196 GetAverageHiGainArea(j).GetHGausHist()->SetTitle("Test Calibrations Area Idx ");
197 GetAverageHiGainArea(j).SetNbins(fAverageNbins);
198 GetAverageHiGainArea(j).InitBins();
199 GetAverageHiGainArea(j).ChangeHistId(j);
200 GetAverageHiGainArea(j).SetEventFrequency(fPulserFrequency);
201
202 TH1F *h = GetAverageHiGainArea(j).GetHGausHist();
203 h->SetTitle( Form("%s%s", h->GetTitle()," Runs: "));
204 }
205 }
206
207 if (fAverageLoGainAreas->GetEntries()==0)
208 {
209 fAverageLoGainAreas->Expand(nareas);
210
211 for (Int_t j=0; j<nareas; j++)
212 {
213 (*fAverageLoGainAreas)[j] =
214 new MHCalibrationTestPix("MHCalibrationTestAverageArea",
215 "Average Test Calibrations Area Idx ");
216
217 GetAverageLoGainArea(j).GetHGausHist()->SetTitle("Test Calibrations Area Idx ");
218 GetAverageLoGainArea(j).SetNbins(fAverageNbins);
219 GetAverageLoGainArea(j).InitBins();
220 GetAverageLoGainArea(j).ChangeHistId(j);
221 GetAverageLoGainArea(j).SetEventFrequency(fPulserFrequency);
222
223 }
224 }
225
226
227 if (fAverageHiGainSectors->GetEntries()==0)
228 {
229 fAverageHiGainSectors->Expand(nsectors);
230
231 for (Int_t j=0; j<nsectors; j++)
232 {
233 (*fAverageHiGainSectors)[j] =
234 new MHCalibrationTestPix("MHCalibrationTestAverageSector",
235 "Average Test Calibrations Sector ");
236
237 GetAverageHiGainSector(j).GetHGausHist()->SetTitle("Test Calibrations Sector ");
238 GetAverageHiGainSector(j).SetNbins(fAverageNbins);
239 GetAverageHiGainSector(j).InitBins();
240 GetAverageHiGainSector(j).ChangeHistId(j);
241 GetAverageHiGainSector(j).SetEventFrequency(fPulserFrequency);
242
243 TH1F *h = GetAverageHiGainSector(j).GetHGausHist();
244 h->SetTitle( Form("%s%s", h->GetTitle()," Runs: "));
245
246 }
247 }
248
249
250 if (fAverageLoGainSectors->GetEntries()==0)
251 {
252 fAverageLoGainSectors->Expand(nsectors);
253
254 for (Int_t j=0; j<nsectors; j++)
255 {
256 (*fAverageLoGainSectors)[j] =
257 new MHCalibrationTestPix("MHCalibrationTestAverageSector",
258 "Average Test Calibrations Sector ");
259
260 GetAverageLoGainSector(j).GetHGausHist()->SetTitle("Test Calibrations Sector ");
261 GetAverageLoGainSector(j).SetNbins(fAverageNbins);
262 GetAverageLoGainSector(j).InitBins();
263 GetAverageLoGainSector(j).ChangeHistId(j);
264 GetAverageLoGainSector(j).SetEventFrequency(fPulserFrequency);
265 }
266 }
267
268
269 fMeanMeanPhotPerArea.Set(nareas);
270 fRmsMeanPhotPerArea .Set(nareas);
271 fMeanSigmaPhotPerArea.Set(nareas);
272 fRmsSigmaPhotPerArea.Set(nareas);
273
274 return kTRUE;
275}
276
277
278// -------------------------------------------------------------------------------
279//
280// Retrieves pointer to MCerPhotEvt:
281//
282// Retrieves from MGeomCam:
283// - number of pixels
284// - number of pixel areas
285// - number of sectors
286//
287// Fills HiGain or LoGain histograms (MHGausEvents::FillHistAndArray()), respectively
288// depending on MCerPhotPix::IsLoGainUsed(), with:
289// - MCerPhotPix::GetArrivalTime(pixid) - MCerPhotPix::GetArrivalTime(1);
290// (i.e. the time difference between pixel i and pixel 1 (hardware number: 2) )
291//
292Bool_t MHCalibrationTestCam::FillHists(const MParContainer *par, const Stat_t w)
293{
294
295 MCerPhotEvt *calibration = (MCerPhotEvt*)par;
296 if (!calibration)
297 {
298 gLog << err << "No argument in MHCalibrationRelTimeCam::Fill... abort." << endl;
299 return kFALSE;
300 }
301
302 const Int_t npixels = fGeom->GetNumPixels();
303 const Int_t nareas = fGeom->GetNumAreas();
304 const Int_t nsectors = fGeom->GetNumSectors();
305
306 Float_t sumareahi [nareas];
307 Float_t sumsectorhi[nsectors];
308 Int_t numareahi [nareas];
309 Int_t numsectorhi[nsectors];
310
311 memset(sumareahi, 0, nareas * sizeof(Float_t));
312 memset(sumsectorhi, 0, nsectors*sizeof(Float_t));
313 memset(numareahi, 0, nareas * sizeof(Int_t));
314 memset(numsectorhi, 0, nsectors*sizeof(Int_t));
315
316 for (Int_t i=0; i<npixels; i++)
317 {
318
319 MHGausEvents &histhi = (*this)[i];
320
321 const MCerPhotPix *pix = calibration->GetPixById(i);
322 if (!pix)
323 continue;
324
325
326 const Float_t signal = pix->GetNumPhotons();
327
328 if (signal < 0.0001)
329 continue;
330
331 const Int_t aidx = (*fGeom)[i].GetAidx();
332 const Int_t sector = (*fGeom)[i].GetSector();
333
334 histhi.FillHistAndArray(signal) ;
335
336 sumareahi [aidx] += signal;
337 numareahi [aidx] ++;
338 sumsectorhi[sector] += signal;
339 numsectorhi[sector] ++;
340 }
341
342 for (Int_t j=0; j<nareas; j++)
343 {
344 MHGausEvents &histhi = GetAverageHiGainArea(j);
345 histhi.FillHistAndArray(numareahi[j] == 0 ? 0. : sumareahi[j]/numareahi[j]);
346 }
347
348 for (Int_t j=0; j<nsectors; j++)
349 {
350 MHGausEvents &histhi = GetAverageHiGainSector(j);
351 histhi.FillHistAndArray(numsectorhi[j] == 0 ? 0. : sumsectorhi[j]/numsectorhi[j]);
352
353 }
354
355 return kTRUE;
356}
357
358// --------------------------------------------------------------------------
359//
360// Calls:
361// - MHCalibrationCam::FitHiGainArrays() with flags:
362// MBadPixelsPix::kTestNotFitted and MBadPixelsPix::kTestOscillating
363// - MHCalibrationCam::FitLoGainArrays() with flags:
364// MBadPixelsPix::kTestNotFitted and MBadPixelsPix::kTestOscillating
365//
366Bool_t MHCalibrationTestCam::FinalizeHists()
367{
368
369 *fLog << endl;
370
371 TArrayI numaidx;
372 numaidx.Set(fGeom->GetNumAreas());
373
374 for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
375 {
376
377 MHGausEvents &hist = (*this)[i];
378
379 if (hist.IsEmpty())
380 {
381 const Int_t size = fNotInterpolateablePixels.GetSize();
382 fNotInterpolateablePixels.Set(size+1);
383 fNotInterpolateablePixels[size] = i;
384 continue;
385 }
386
387 if (!hist.FitGaus())
388 if (!hist.RepeatFit())
389 hist.BypassFit();
390
391 hist.CreateFourierSpectrum();
392
393 const Float_t area = (*fGeom)[i].GetA();
394 const Int_t aidx = (*fGeom)[i].GetAidx();
395
396 fMeanMeanPhotPerArea[aidx] += hist.GetMean() / area;
397 fRmsMeanPhotPerArea [aidx] += hist.GetMean() / area * hist.GetMean() / area;
398 fMeanSigmaPhotPerArea[aidx] += hist.GetSigma()/ area;
399 fRmsSigmaPhotPerArea [aidx] += hist.GetSigma()/ area * hist.GetSigma() / area;
400 numaidx[aidx]++;
401 }
402
403
404 for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++)
405 {
406
407 MHGausEvents &hist = GetAverageHiGainArea(j);
408 if (hist.IsEmpty())
409 continue;
410
411 if (!hist.FitGaus())
412 if (!hist.RepeatFit())
413 hist.BypassFit();
414
415 hist.CreateFourierSpectrum();
416
417 fRmsMeanPhotPerArea [j] -= fMeanMeanPhotPerArea [j]*fMeanMeanPhotPerArea [j]/numaidx[j];
418 fRmsSigmaPhotPerArea[j] -= fMeanSigmaPhotPerArea[j]*fMeanSigmaPhotPerArea[j]/numaidx[j];
419
420 fMeanMeanPhotPerArea [j] /= numaidx[j];
421 fMeanSigmaPhotPerArea[j] /= numaidx[j];
422 fRmsMeanPhotPerArea [j] /= numaidx[j]-1.;
423 fRmsSigmaPhotPerArea [j] /= numaidx[j]-1.;
424
425 if (fRmsMeanPhotPerArea [j] > 0.)
426 fRmsMeanPhotPerArea [j] = TMath::Sqrt(fRmsMeanPhotPerArea [j]);
427 if (fRmsSigmaPhotPerArea [j] > 0.)
428 fRmsSigmaPhotPerArea [j] = TMath::Sqrt(fRmsSigmaPhotPerArea [j]);
429 }
430
431 for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
432 {
433
434 MHGausEvents &hist = GetAverageHiGainSector(j);
435 if (hist.IsEmpty())
436 continue;
437
438 if (!hist.FitGaus())
439 if (!hist.RepeatFit())
440 hist.BypassFit();
441
442 hist.CreateFourierSpectrum();
443 }
444
445 return kTRUE;
446}
447
448
449// --------------------------------------------------------------------------
450//
451// The types are as follows:
452//
453// Fitted values:
454// ==============
455//
456// 0: Fitted Mean Test Calibration (MHGausEvents::GetMean())
457// 1: Error Mean Test Calibration (MHGausEvents::GetMeanErr())
458// 2: Sigma fitted Test Calibration (MHGausEvents::GetSigma())
459// 3: Error Sigma Test Calibration (MHGausEvents::GetSigmaErr())
460//
461// Useful variables derived from the fit results:
462// =============================================
463//
464// 4: Returned probability of Gauss fit (calls: MHGausEvents::GetProb())
465//
466// Localized defects:
467// ==================
468//
469// 5: Gaus fit not OK (calls: MHGausEvents::IsGausFitOK())
470// 6: Fourier spectrum not OK (calls: MHGausEvents::IsFourierSpectrumOK())
471//
472// Converted values:
473// =================
474//
475// 7: Fitted Mean Test Calibration (MHGausEvents::GetMean()) by MGeomPix::GetA()
476// 8: Fitted Mean Error Calibration (MHGausEvents::GetMeanErr()) by MGeomPix::GetA()
477// 9: Fitted Sigma Test Calibration (MHGausEvents::GetSigma()) by MGeomPix::GetA()
478// 10: Fitted Sigma Error Calibration (MHGausEvents::GetSigmaErr()) by MGeomPix::GetA()
479//
480Bool_t MHCalibrationTestCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
481{
482
483 if (fHiGainArray->GetSize() <= idx)
484 return kFALSE;
485
486 const MHGausEvents &pix = (*this)[idx];
487
488 if (pix.IsEmpty())
489 return kFALSE;
490
491 switch (type)
492 {
493 case 0:
494 val = pix.GetMean();
495 break;
496 case 1:
497 val = pix.GetMeanErr();
498 break;
499 case 2:
500 val = pix.GetSigma();
501 break;
502 case 3:
503 val = pix.GetSigmaErr();
504 break;
505 case 4:
506 val = pix.GetProb();
507 break;
508 case 5:
509 if (!pix.IsGausFitOK())
510 val = 1.;
511 break;
512 case 6:
513 if (!pix.IsFourierSpectrumOK())
514 val = 1.;
515 break;
516 case 7:
517 val = pix.GetMean()/cam[idx].GetA();
518 break;
519 case 8:
520 val = pix.GetMeanErr()/cam[idx].GetA();
521 break;
522 case 9:
523 val = pix.GetSigma()/cam[idx].GetA();
524 break;
525 case 10:
526 val = pix.GetSigmaErr()/cam[idx].GetA();
527 break;
528 default:
529 return kFALSE;
530 }
531 return kTRUE;
532}
533
534// --------------------------------------------------------------------------
535//
536// Calls MHGausEvents::DrawClone() for pixel idx
537//
538void MHCalibrationTestCam::DrawPixelContent(Int_t idx) const
539{
540 (*this)[idx].DrawClone();
541}
542
543
544//------------------------------------------------------------
545//
546// For all averaged areas, the fitted sigma is multiplied with the square root of
547// the number involved pixels
548//
549void MHCalibrationTestCam::CalcAverageSigma()
550{
551
552 for (UInt_t j=0; j<fGeom->GetNumAreas(); j++)
553 {
554
555 MHGausEvents &hist = GetAverageHiGainArea(j);
556
557 const Float_t numsqr = TMath::Sqrt((Float_t)fAverageAreaNum[j]);
558 fAverageAreaSigma[j] = hist.GetSigma () * numsqr;
559 fAverageAreaSigmaVar[j] = hist.GetSigmaErr () * hist.GetSigmaErr() * numsqr;
560
561 fAverageAreaRelSigma [j] = fAverageAreaSigma[j] / hist.GetMean();
562 fAverageAreaRelSigmaVar[j] = fAverageAreaSigmaVar[j] / (fAverageAreaSigma[j]*fAverageAreaSigma[j]);
563 fAverageAreaRelSigmaVar[j] += hist.GetMeanErr()*hist.GetMeanErr()/hist.GetMean()/hist.GetMean();
564 fAverageAreaRelSigmaVar[j] *= fAverageAreaRelSigma[j];
565 }
566}
Note: See TracBrowser for help on using the repository browser.