source: trunk/MagicSoft/Mars/manalysis/MHPedestalCam.cc@ 5977

Last change on this file since 5977 was 4964, checked in by gaug, 20 years ago
*** empty log message ***
File size: 17.4 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
115ClassImp(MHPedestalCam);
116
117using namespace std;
118// --------------------------------------------------------------------------
119//
120// Default constructor.
121//
122// Initializes:
123// - fExtractHiGainSlices to 0.
124// - fExtractLoGainSlices to 0.
125// - the event frequency to 1200 Hz.
126// - the fRenorm flag to kTRUE
127//
128MHPedestalCam::MHPedestalCam(const char *name, const char *title)
129 : fExtractHiGainSlices(0.), fExtractLoGainSlices(0.)
130{
131
132 fName = name ? name : "MHPedestalCam";
133 fTitle = title ? title : "";
134
135 SetPulserFrequency(1200);
136 SetRenorm();
137
138}
139
140
141
142// --------------------------------------------------------------------------
143//
144// Searches pointer to:
145// - MPedestalCam
146// - MExtractedSignalCam
147//
148// Searches or creates:
149// - MCalibrationPedCam
150//
151// Retrieves from MExtractedSignalCam:
152// - number of used High Gain FADC slices (MExtractedSignalCam::GetNumUsedHiGainFADCSlices())
153// - number of used Low Gain FADC slices (MExtractedSignalCam::GetNumUsedLoGainFADCSlices())
154//
155// Initializes, if empty to MGeomCam::GetNumPixels():
156// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
157//
158// Initializes, if empty to MGeomCam::GetNumAreas() for:
159// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
160//
161// Initializes, if empty to MGeomCam::GetNumSectors() for:
162// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
163//
164// Calls MHCalibrationCam::InitPedHists() for every entry in:
165// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray
166// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas
167// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors
168//
169// Sets Titles and Names for the Histograms
170// - MHCalibrationCam::fAverageHiGainAreas
171// - MHCalibrationCam::fAverageHiGainSectors
172//
173Bool_t MHPedestalCam::ReInitHists(MParList *pList)
174{
175
176 MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
177 if (!signal)
178 {
179 gLog << err << "Cannot find MExtractedSignalCam... abort." << endl;
180 return kFALSE;
181 }
182
183
184 fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam");
185
186 if (!fPedestals)
187 {
188 gLog << err << "Cannot find MPedestalCam ... abort." << endl;
189 return kFALSE;
190 }
191
192 const Int_t npixels = fGeom->GetNumPixels();
193 const Int_t nsectors = fGeom->GetNumSectors();
194 const Int_t nareas = fGeom->GetNumAreas();
195
196 fCam = (MCalibrationCam*)pList->FindObject("MCalibrationPedCam");
197 if (!fCam)
198 {
199 fCam = (MCalibrationCam*)pList->FindCreateObj(AddSerialNumber("MCalibrationPedCam"));
200 if (!fCam)
201 return kFALSE;
202
203 fCam->Init(*fGeom);
204 }
205
206 Float_t sliceshi = signal->GetNumUsedHiGainFADCSlices();
207 Float_t sliceslo = signal->GetNumUsedLoGainFADCSlices();
208
209 if (sliceshi == 0.)
210 {
211 gLog << err << "Number of used signal slices in MExtractedSignalCam is zero ... abort."
212 << endl;
213 return kFALSE;
214 }
215
216 if (fExtractHiGainSlices != 0. && sliceshi != fExtractHiGainSlices )
217 {
218 gLog << err << "Number of used High Gain signal slices changed in MExtractedSignalCam ... abort."
219 << endl;
220 return kFALSE;
221 }
222
223 if (fExtractLoGainSlices != 0. && sliceslo != fExtractLoGainSlices )
224 {
225 gLog << err << "Number of used Low Gain signal slices changed in MExtractedSignalCam ... abort."
226 << endl;
227 return kFALSE;
228 }
229
230 fExtractHiGainSlices = sliceshi;
231 fExtractLoGainSlices = sliceslo;
232
233 if (fHiGainArray->GetEntries()==0)
234 {
235 fHiGainArray->Expand(npixels);
236 for (Int_t i=0; i<npixels; i++)
237 {
238 (*fHiGainArray)[i] = new MHPedestalPix;
239 InitPedHists((MHPedestalPix&)(*this)[i],i,fExtractHiGainSlices);
240
241 if ((*fBadPixels)[i].IsBad())
242 (*this)[i].SetExcluded();
243
244 (*fPedestals)[i].InitUseHists();
245 }
246
247 }
248
249 if (fLoGainArray->GetEntries()==0)
250 {
251 fLoGainArray->Expand(npixels);
252 for (Int_t i=0; i<npixels; i++)
253 {
254 (*fLoGainArray)[i] = new MHPedestalPix;
255 InitPedHists((MHPedestalPix&)(*this)(i),i,fExtractLoGainSlices);
256
257 if ((*fBadPixels)[i].IsBad())
258 (*this)(i).SetExcluded();
259 }
260 }
261
262 if (fAverageHiGainAreas->GetEntries()==0)
263 {
264 fAverageHiGainAreas->Expand(nareas);
265
266 for (Int_t j=0; j<nareas; j++)
267 {
268 (*fAverageHiGainAreas)[j] =
269 new MHPedestalPix("AverageHiGainArea",
270 "Average Pedestals area idx ");
271
272 GetAverageHiGainArea(j).GetHGausHist()->SetTitle("Pedestals average Area Idx ");
273
274 InitPedHists((MHPedestalPix&)GetAverageHiGainArea(j),j,fExtractHiGainSlices);
275
276 }
277 }
278
279 if (fAverageLoGainAreas->GetEntries()==0)
280 {
281 fAverageLoGainAreas->Expand(nareas);
282
283 for (Int_t j=0; j<nareas; j++)
284 {
285 (*fAverageLoGainAreas)[j] =
286 new MHPedestalPix("AverageLoGainArea",
287 "Pedestals average Area idx ");
288
289 GetAverageLoGainArea(j).GetHGausHist()->SetTitle("Pedestals average Area Idx ");
290
291 InitPedHists((MHPedestalPix&)GetAverageLoGainArea(j),j,fExtractLoGainSlices);
292
293 }
294 }
295
296 if (fAverageHiGainSectors->GetEntries()==0)
297 {
298 fAverageHiGainSectors->Expand(nsectors);
299
300 for (Int_t j=0; j<nsectors; j++)
301 {
302 (*fAverageHiGainSectors)[j] =
303 new MHPedestalPix("AverageHiGainSector",
304 "Pedestals average sector ");
305
306 GetAverageHiGainSector(j).GetHGausHist()->SetTitle("Pedestals average Sector ");
307
308 InitPedHists((MHPedestalPix&)GetAverageHiGainSector(j),j,fExtractHiGainSlices);
309
310 }
311 }
312
313 if (fAverageLoGainSectors->GetEntries()==0)
314 {
315 fAverageLoGainSectors->Expand(nsectors);
316
317 for (Int_t j=0; j<nsectors; j++)
318 {
319 (*fAverageLoGainSectors)[j] =
320 new MHPedestalPix("AverageLoGainSector",
321 "Pedestals average sector ");
322
323 GetAverageLoGainSector(j).GetHGausHist()->SetTitle("Pedestals average Sector ");
324
325 InitPedHists((MHPedestalPix&)GetAverageLoGainSector(j),j,fExtractLoGainSlices);
326
327 }
328 }
329
330 return kTRUE;
331}
332
333
334// -------------------------------------------------------------
335//
336// If MBadPixelsPix::IsBad():
337// - calls MHCalibrationPix::SetExcluded()
338//
339// Calls:
340// - MHGausEvents::InitBins()
341// - MHCalibrationPix::ChangeHistId(i)
342// - MHCalibrationPix::SetEventFrequency(fPulserFrequency)
343// - MHPedestalPix::SetNSlices(nslices)
344//
345void MHPedestalCam::InitPedHists(MHPedestalPix &hist, const Int_t i, const Float_t nslices)
346{
347
348 hist.InitBins();
349 hist.ChangeHistId(i);
350 hist.SetEventFrequency(fPulserFrequency);
351
352 if (fRenorm)
353 hist.SetNSlices(nslices);
354
355 hist.SetProbLimit(0.);
356
357 TH1F *h = hist.GetHGausHist();
358 h->SetTitle( Form("%s%s", h->GetTitle()," Runs: "));
359}
360// -------------------------------------------------------------------------------
361//
362// Retrieves pointer to MExtractedSignalCam:
363//
364// Retrieves from MGeomCam:
365// - number of pixels
366// - number of pixel areas
367// - number of sectors
368//
369// Fills HiGain or LoGain histograms (MHGausEvents::FillHistAndArray()), respectively
370// with the signals MExtractedSignalCam::GetExtractedSignalHiGain() and
371// MExtractedSignalCam::GetExtractedSignalLoGain(), respectively.
372//
373Bool_t MHPedestalCam::FillHists(const MParContainer *par, const Stat_t w)
374{
375
376 MExtractedSignalCam *signal = (MExtractedSignalCam*)par;
377 if (!signal)
378 {
379 gLog << err << "No argument in MExtractedSignalCam::Fill... abort." << endl;
380 return kFALSE;
381 }
382
383
384
385 const UInt_t npixels = fGeom->GetNumPixels();
386 const UInt_t nareas = fGeom->GetNumAreas();
387 const UInt_t nsectors = fGeom->GetNumSectors();
388
389 TArrayF sumareahi(nareas);
390 TArrayF sumarealo(nareas);
391 TArrayF sumsectorhi(nsectors);
392 TArrayD sumsectorlo(nsectors);
393 TArrayI numareahi(nareas);
394 TArrayI numarealo(nareas);
395 TArrayI numsectorhi(nsectors);
396 TArrayI numsectorlo(nsectors);
397
398 for (UInt_t i=0; i<npixels; i++)
399 {
400 MHCalibrationPix &histhi = (*this)[i];
401 MHCalibrationPix &histlo = (*this)(i);
402
403 if (histhi.IsExcluded())
404 continue;
405
406 const MExtractedSignalPix &pix = (*signal)[i];
407
408 const Float_t pedhi = pix.GetExtractedSignalHiGain();
409 const Float_t pedlo = pix.GetExtractedSignalLoGain();
410
411 const Int_t aidx = (*fGeom)[i].GetAidx();
412 const Int_t sector = (*fGeom)[i].GetSector();
413
414 histhi.FillHistAndArray(pedhi) ;
415 sumareahi [aidx] += pedhi;
416 numareahi [aidx] ++;
417 sumsectorhi[sector] += pedhi;
418 numsectorhi[sector] ++;
419
420 histlo.FillHistAndArray(pedlo);
421 sumarealo [aidx] += pedlo;
422 numarealo [aidx] ++;
423 sumsectorlo[sector] += pedlo;
424 numsectorlo[sector] ++;
425
426 }
427
428 for (UInt_t j=0; j<nareas; j++)
429 {
430 MHCalibrationPix &histhi = GetAverageHiGainArea(j);
431 histhi.FillHistAndArray(numareahi[j] == 0 ? 0. : sumareahi[j]/numareahi[j]);
432
433 MHCalibrationPix &histlo = GetAverageLoGainArea(j);
434 histlo.FillHistAndArray(numarealo[j] == 0 ? 0. : sumarealo[j]/numarealo[j]);
435 }
436
437 for (UInt_t j=0; j<nsectors; j++)
438 {
439 MHCalibrationPix &histhi = GetAverageHiGainSector(j);
440 histhi.FillHistAndArray(numsectorhi[j] == 0 ? 0. : sumsectorhi[j]/numsectorhi[j]);
441
442 MHCalibrationPix &histlo = GetAverageLoGainSector(j);
443 histlo.FillHistAndArray(numsectorlo[j] == 0 ? 0. : sumsectorlo[j]/numsectorlo[j]);
444 }
445
446 return kTRUE;
447}
448
449// --------------------------------------------------------------------------
450//
451// Calls:
452// - MHCalibrationCam::FitHiGainArrays() with Bad Pixels flags 0
453// - MHCalibrationCam::FitLoGainArrays() with Bad Pixels flags 0
454//
455Bool_t MHPedestalCam::FinalizeHists()
456{
457
458 FitHiGainArrays((*fCam),*fBadPixels,
459 MBadPixelsPix::kHiGainNotFitted,
460 MBadPixelsPix::kHiGainOscillating);
461 FitLoGainArrays((*fCam),*fBadPixels,
462 MBadPixelsPix::kLoGainNotFitted,
463 MBadPixelsPix::kLoGainOscillating);
464
465
466 return kTRUE;
467
468
469}
470
471// ------------------------------------------------------------------
472//
473// The types are as follows:
474//
475// Fitted values:
476// ==============
477//
478// 0: Fitted Charge
479// 1: Error of fitted Charge
480// 2: Sigma of fitted Charge
481// 3: Error of Sigma of fitted Charge
482//
483//
484// Useful variables derived from the fit results:
485// =============================================
486//
487// 4: Returned probability of Gauss fit to Charge distribution
488// 5: Relative differenc of calculated pedestal (per slice) and fitted (per slice)
489// 6: Error of the Relative differenc of calculated pedestal (per slice) and fitted (per slice)
490// 7: Relative difference of the error of the mean pedestal (per slice) - calculated and fitted
491// 8: Relative differenc of calculated pedestal RMS (per slice) and fitted sigma (per slice)
492// 9: Error of Relative differenc of calculated pedestal RMS (per slice) and fitted sigma (per slice)
493// 10: Relative difference of the error of the pedestal RMS (per slice) - calculated and fitted
494//
495// Localized defects:
496// ==================
497//
498// 11: Gaus fit not OK
499// 12: Fourier spectrum not OK
500//
501Bool_t MHPedestalCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
502{
503
504 if (fHiGainArray->GetSize() <= idx)
505 return kFALSE;
506
507 if ((*this)[idx].IsExcluded())
508 return kFALSE;
509
510 const Float_t ped = (*fPedestals)[idx].GetPedestal();
511 const Float_t rms = (*fPedestals)[idx].GetPedestalRms();
512
513 const Float_t entsqr = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries());
514
515 const Float_t pederr = rms/entsqr;
516 const Float_t rmserr = rms/entsqr/2.;
517
518 const Float_t mean = (*this)[idx].GetMean();
519 const Float_t meanerr = (*this)[idx].GetMeanErr();
520 const Float_t sigma = (*this)[idx].GetSigma() ;
521 const Float_t sigmaerr = (*this)[idx].GetSigmaErr();
522
523 switch (type)
524 {
525 case 0:
526 val = mean;
527 break;
528 case 1:
529 val = meanerr;
530 break;
531 case 2:
532 val = sigma;
533 break;
534 case 3:
535 val = sigmaerr;
536 break;
537 case 4:
538 val = (*this)[idx].GetProb();
539 break;
540 case 5:
541 val = 2.*(mean-ped)/(ped+mean);
542 break;
543 case 6:
544 val = TMath::Sqrt((pederr*pederr + meanerr*meanerr) * (ped*ped + mean*mean))
545 *2./(ped+mean)/(ped+mean);
546 break;
547 case 7:
548 val = 2.*(meanerr-pederr)/(pederr + meanerr);
549 break;
550 case 8:
551 val = 2.*(sigma-rms)/(sigma+rms);
552 break;
553 case 9:
554 val = TMath::Sqrt((rmserr*rmserr + sigmaerr*sigmaerr) * (rms*rms + sigma*sigma))
555 *2./(rms+sigma)/(rms+sigma);
556 break;
557 case 10:
558 val = 2.*(sigmaerr - rmserr)/(sigmaerr + rmserr);
559 break;
560 case 11:
561 if (!(*this)[idx].IsGausFitOK())
562 val = 1.;
563 break;
564 case 12:
565 if (!(*this)[idx].IsFourierSpectrumOK())
566 val = 1.;
567 break;
568 default:
569 return kFALSE;
570 }
571 return kTRUE;
572}
573
574// --------------------------------------------------------------------------
575//
576// Calls MHGausEvents::DrawClone() for pixel idx
577//
578void MHPedestalCam::DrawPixelContent(Int_t idx) const
579{
580 (*this)[idx].DrawClone();
581}
Note: See TracBrowser for help on using the repository browser.