source: trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestCam.cc@ 7593

Last change on this file since 7593 was 6857, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 14.3 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 MSignalCam into
29// MHCalibrationPix 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// +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit())
52// In case this does not make the fit valid, the histogram means and RMS's are
53// taken directly (see MHCalibrationPix::BypassFit()) and the following flags are set:
54// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) and
55// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun )
56//
57// Outliers of more than MHCalibrationPix::fPickupLimit (default: 5) sigmas
58// from the mean are counted as Pickup events (stored in MHCalibrationPix::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
86#include "MHCalibrationPix.h"
87
88#include "MLog.h"
89#include "MLogManip.h"
90
91#include "MParList.h"
92
93#include "MCalibrationCam.h"
94#include "MCalibrationPix.h"
95
96#include "MSignalCam.h"
97#include "MSignalPix.h"
98
99#include "MGeomCam.h"
100#include "MGeomPix.h"
101
102#include "MBadPixelsCam.h"
103#include "MBadPixelsPix.h"
104
105#include <TOrdCollection.h>
106
107ClassImp(MHCalibrationTestCam);
108
109using namespace std;
110
111const Int_t MHCalibrationTestCam::fgNbins = 1000;
112const Axis_t MHCalibrationTestCam::fgFirst = -1.;
113const Axis_t MHCalibrationTestCam::fgLast = 1999.;
114const Float_t MHCalibrationTestCam::fgProbLimit = 0.00000001;
115const TString MHCalibrationTestCam::gsHistName = "Test";
116const TString MHCalibrationTestCam::gsHistTitle = "Calibrated Calibration Signals";
117const TString MHCalibrationTestCam::gsHistXTitle = "Nr. Photons";
118const TString MHCalibrationTestCam::gsHistYTitle = "Nr. events";
119// --------------------------------------------------------------------------
120//
121// Default Constructor.
122//
123// Sets:
124// - fNbins to fgNbins
125// - fFirst to fgFirst
126// - fLast to fgLast
127//
128// - fHistName to gsHistName
129// - fHistTitle to gsHistTitle
130// - fHistXTitle to gsHistXTitle
131// - fHistYTitle to gsHistYTitle
132//
133MHCalibrationTestCam::MHCalibrationTestCam(const char *name, const char *title)
134{
135
136 fName = name ? name : "MHCalibrationTestCam";
137 fTitle = title ? title : "Histogram class for testing the calibration";
138
139 SetNbins(fgNbins);
140 SetFirst(fgFirst);
141 SetLast (fgLast );
142
143 SetProbLimit(fgProbLimit);
144
145 SetHistName (gsHistName .Data());
146 SetHistTitle (gsHistTitle .Data());
147 SetHistXTitle(gsHistXTitle.Data());
148 SetHistYTitle(gsHistYTitle.Data());
149
150 SetLoGain(kFALSE);
151
152}
153
154// --------------------------------------------------------------------------
155//
156// Searches pointer to:
157// - MSignalCam
158//
159// Calls:
160// - MHCalibrationCam::InitHiGainArrays()
161//
162// Sets:
163// - SetLoGain(kFALSE);
164// - fMeanMeanPhotPerArea to nareas
165// - fRmsMeanPhotPerArea to nareas
166// - fMeanSigmaPhotPerArea to nareas
167// - fRmsSigmaPhotPerArea to nareas
168//
169Bool_t MHCalibrationTestCam::ReInitHists(MParList *pList)
170{
171
172 if (!InitCams(pList,"Test"))
173 return kFALSE;
174
175
176 MSignalCam *signal = (MSignalCam*)pList->FindObject("MSignalCam");
177 if (!signal)
178 {
179 *fLog << err << "MSignalCam not found... abort." << endl;
180 return kFALSE;
181 }
182
183
184 const Int_t npixels = fGeom->GetNumPixels();
185 const Int_t nsectors = fGeom->GetNumSectors();
186 const Int_t nareas = fGeom->GetNumAreas();
187
188 InitHiGainArrays(npixels,nareas,nsectors);
189 InitLoGainArrays(npixels,nareas,nsectors);
190
191 fMeanMeanPhotPerArea.Set(nareas);
192 fRmsMeanPhotPerArea .Set(nareas);
193 fMeanSigmaPhotPerArea.Set(nareas);
194 fRmsSigmaPhotPerArea.Set(nareas);
195
196 return kTRUE;
197}
198
199// -------------------------------------------------------------------------------
200//
201// Retrieves pointer to MSignalCam:
202//
203// Retrieves from MGeomCam:
204// - number of pixels
205// - number of pixel areas
206// - number of sectors
207//
208// Fills HiGain histograms (MHGausEvents::FillHistAndArray())
209// with:
210// - MSignalPix::GetNumPhotons(pixid);
211//
212Bool_t MHCalibrationTestCam::FillHists(const MParContainer *par, const Stat_t w)
213{
214
215 MSignalCam *calibration = (MSignalCam*)par;
216 if (!calibration)
217 {
218 gLog << err << "No argument in MHCalibrationTestCam::Fill... abort." << endl;
219 return kFALSE;
220 }
221
222 const Int_t npixels = fGeom->GetNumPixels();
223 const Int_t nareas = fGeom->GetNumAreas();
224 const Int_t nsectors = fGeom->GetNumSectors();
225
226 TArrayF sumareahi (nareas);
227 TArrayF sumsectorhi(nsectors);
228 TArrayI numareahi (nareas);
229 TArrayI numsectorhi(nsectors);
230
231 for (Int_t i=0; i<npixels; i++)
232 {
233
234 MHCalibrationPix &histhi = (*this)[i];
235
236 const MSignalPix *pix = calibration->GetPixById(i);
237 if (!pix)
238 continue;
239
240 const Float_t signal = pix->GetNumPhotons();
241
242 if (signal < 0.0001)
243 continue;
244
245 const Int_t aidx = (*fGeom)[i].GetAidx();
246 const Int_t sector = (*fGeom)[i].GetSector();
247
248 histhi.FillHistAndArray(signal) ;
249
250 sumareahi [aidx] += signal;
251 numareahi [aidx] ++;
252 sumsectorhi[sector] += signal;
253 numsectorhi[sector] ++;
254 }
255
256 for (Int_t j=0; j<nareas; j++)
257 {
258 MHCalibrationPix &histhi = GetAverageHiGainArea(j);
259 histhi.FillHistAndArray(numareahi[j] == 0 ? 0. : sumareahi[j]/numareahi[j]);
260 }
261
262 for (Int_t j=0; j<nsectors; j++)
263 {
264 MHCalibrationPix &histhi = GetAverageHiGainSector(j);
265 histhi.FillHistAndArray(numsectorhi[j] == 0 ? 0. : sumsectorhi[j]/numsectorhi[j]);
266
267 }
268
269 return kTRUE;
270}
271
272// --------------------------------------------------------------------------
273//
274// Calls:
275// - MHCalibrationCam::FitHiGainArrays() with flags:
276// MBadPixelsPix::kTestNotFitted and MBadPixelsPix::kTestOscillating
277//
278Bool_t MHCalibrationTestCam::FinalizeHists()
279{
280
281 *fLog << endl;
282
283 TArrayI numaidx;
284 numaidx.Set(fGeom->GetNumAreas());
285
286 for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
287 {
288
289 MHCalibrationPix &hist = (*this)[i];
290
291 if (hist.IsEmpty())
292 continue;
293
294 if (!hist.FitGaus())
295 if (!hist.RepeatFit())
296 hist.BypassFit();
297
298 hist.CreateFourierSpectrum();
299
300 const Float_t area = (*fGeom)[i].GetA();
301 const Int_t aidx = (*fGeom)[i].GetAidx();
302
303 fMeanMeanPhotPerArea[aidx] += hist.GetMean() / area;
304 fRmsMeanPhotPerArea [aidx] += hist.GetMean() / area * hist.GetMean() / area;
305 fMeanSigmaPhotPerArea[aidx] += hist.GetSigma()/ area;
306 fRmsSigmaPhotPerArea [aidx] += hist.GetSigma()/ area * hist.GetSigma() / area;
307 numaidx[aidx]++;
308 }
309
310
311 for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++)
312 {
313
314 MHCalibrationPix &hist = GetAverageHiGainArea(j);
315 if (hist.IsEmpty())
316 continue;
317
318 if (!hist.FitGaus())
319 if (!hist.RepeatFit())
320 hist.BypassFit();
321
322 hist.CreateFourierSpectrum();
323
324 fRmsMeanPhotPerArea [j] -= fMeanMeanPhotPerArea [j]*fMeanMeanPhotPerArea [j]/numaidx[j];
325 fRmsSigmaPhotPerArea[j] -= fMeanSigmaPhotPerArea[j]*fMeanSigmaPhotPerArea[j]/numaidx[j];
326
327 fMeanMeanPhotPerArea [j] /= numaidx[j];
328 fMeanSigmaPhotPerArea[j] /= numaidx[j];
329 fRmsMeanPhotPerArea [j] /= numaidx[j]-1.;
330 fRmsSigmaPhotPerArea [j] /= numaidx[j]-1.;
331
332 if (fRmsMeanPhotPerArea [j] > 0.)
333 fRmsMeanPhotPerArea [j] = TMath::Sqrt(fRmsMeanPhotPerArea [j]);
334 if (fRmsSigmaPhotPerArea [j] > 0.)
335 fRmsSigmaPhotPerArea [j] = TMath::Sqrt(fRmsSigmaPhotPerArea [j]);
336 }
337
338 for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
339 {
340
341 MHCalibrationPix &hist = GetAverageHiGainSector(j);
342 if (hist.IsEmpty())
343 continue;
344
345 if (!hist.FitGaus())
346 if (!hist.RepeatFit())
347 hist.BypassFit();
348
349 hist.CreateFourierSpectrum();
350 }
351
352 return kTRUE;
353}
354
355
356// --------------------------------------------------------------------------
357//
358// The types are as follows:
359//
360// Fitted values:
361// ==============
362//
363// 0: Fitted Mean Test Calibration (MHGausEvents::GetMean())
364// 1: Error Mean Test Calibration (MHGausEvents::GetMeanErr())
365// 2: Sigma fitted Test Calibration (MHGausEvents::GetSigma())
366// 3: Error Sigma Test Calibration (MHGausEvents::GetSigmaErr())
367//
368// Useful variables derived from the fit results:
369// =============================================
370//
371// 4: Returned probability of Gauss fit (calls: MHGausEvents::GetProb())
372//
373// Localized defects:
374// ==================
375//
376// 5: Gaus fit not OK (calls: MHGausEvents::IsGausFitOK())
377// 6: Fourier spectrum not OK (calls: MHGausEvents::IsFourierSpectrumOK())
378//
379// Converted values:
380// =================
381//
382// 7: Fitted Mean Test Calibration (MHGausEvents::GetMean()) by MGeomPix::GetA()
383// 8: Fitted Mean Error Calibration (MHGausEvents::GetMeanErr()) by MGeomPix::GetA()
384// 9: Fitted Sigma Test Calibration (MHGausEvents::GetSigma()) by MGeomPix::GetA()
385// 10: Fitted Sigma Error Calibration (MHGausEvents::GetSigmaErr()) by MGeomPix::GetA()
386//
387Bool_t MHCalibrationTestCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
388{
389
390 if (fHiGainArray->GetSize() <= idx)
391 return kFALSE;
392
393 const MHCalibrationPix &pix = (*this)[idx];
394
395 if (pix.IsEmpty())
396 return kFALSE;
397
398 switch (type)
399 {
400 case 0:
401 val = pix.GetMean();
402 break;
403 case 1:
404 val = pix.GetMeanErr();
405 break;
406 case 2:
407 val = pix.GetSigma();
408 break;
409 case 3:
410 val = pix.GetSigmaErr();
411 break;
412 case 4:
413 val = pix.GetProb();
414 break;
415 case 5:
416 if (!pix.IsGausFitOK())
417 val = 1.;
418 break;
419 case 6:
420 if (!pix.IsFourierSpectrumOK())
421 val = 1.;
422 break;
423 case 7:
424 val = pix.GetMean()/cam[idx].GetA();
425 break;
426 case 8:
427 val = pix.GetMeanErr()/cam[idx].GetA();
428 break;
429 case 9:
430 val = pix.GetSigma()/cam[idx].GetA();
431 break;
432 case 10:
433 val = pix.GetSigmaErr()/cam[idx].GetA();
434 break;
435 default:
436 return kFALSE;
437 }
438 return kTRUE;
439}
440
441// --------------------------------------------------------------------------
442//
443// Calls MHCalibrationPix::DrawClone() for pixel idx
444//
445void MHCalibrationTestCam::DrawPixelContent(Int_t idx) const
446{
447 (*this)[idx].DrawClone();
448}
449
450
451//------------------------------------------------------------
452//
453// For all averaged areas, the fitted sigma is multiplied with the square root of
454// the number involved pixels
455//
456void MHCalibrationTestCam::CalcAverageSigma()
457{
458
459 for (UInt_t j=0; j<fGeom->GetNumAreas(); j++)
460 {
461
462 MHCalibrationPix &hist = GetAverageHiGainArea(j);
463
464 const Float_t numsqr = TMath::Sqrt((Float_t)fAverageAreaNum[j]);
465 fAverageAreaSigma[j] = hist.GetSigma () * numsqr;
466 fAverageAreaSigmaVar[j] = hist.GetSigmaErr () * hist.GetSigmaErr() * numsqr;
467
468 fAverageAreaRelSigma [j] = fAverageAreaSigma[j] / hist.GetMean();
469 fAverageAreaRelSigmaVar[j] = fAverageAreaSigmaVar[j] / (fAverageAreaSigma[j]*fAverageAreaSigma[j]);
470 fAverageAreaRelSigmaVar[j] += hist.GetMeanErr()*hist.GetMeanErr()/hist.GetMean()/hist.GetMean();
471 fAverageAreaRelSigmaVar[j] *= fAverageAreaRelSigma[j];
472 }
473}
Note: See TracBrowser for help on using the repository browser.