source: trunk/MagicSoft/Mars/mpedestal/MHPedestalCam.cc@ 5817

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