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

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