source: trunk/MagicSoft/Mars/mhcalib/MHPedestalCam.cc@ 8761

Last change on this file since 8761 was 8504, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 26.2 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//
27// MHPedestalCam
28//
29// Fills the extracted pedestals of MExtractedSignalCam into
30// the MHCalibrationPix-classes MHPedestalPix for every:
31//
32// - Pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray
33// or MHCalibrationCam::fHiGainArray, respectively
34//
35// - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera),
36// stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas and
37// MHCalibrationCam::fAverageHiGainAreas
38//
39// - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera),
40// stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors
41// and MHCalibrationCam::fAverageHiGainSectors
42//
43// Every pedestal is directly taken from MExtractedSignalCam, filled by the
44// appropriate extractor.
45//
46// The pedestals are filled into a histogram and an array, in order to perform
47// a Fourier analysis (see MHGausEvents). The signals are moreover averaged on an
48// event-by-event basis and written into the corresponding average pixels.
49//
50// The histograms are fitted to a Gaussian, mean and sigma with its errors
51// and the fit probability are extracted. If none of these values are NaN's and
52// if the probability is bigger than MHGausEvents::fProbLimit (default: 0.5%),
53// the fit is declared valid.
54// Otherwise, the fit is repeated within ranges of the previous mean
55// +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit())
56// In case this does not make the fit valid, the histogram means and RMS's are
57// taken directly (see MHCalibrationPix::BypassFit()).
58//
59// Outliers of more than MHCalibrationPix::fPickupLimit (default: 5) sigmas
60// from the mean are counted as Pickup events (stored in MHCalibrationPix::fPickup)
61//
62// The number of summed FADC slices is taken to re-normalize
63// the result to a pedestal per pixel with the formulas (see MHPedestalPix::Renorm()):
64// - Mean Pedestal / slice = Mean Pedestal / Number slices
65// - Mean Pedestal Error / slice = Mean Pedestal Error / Number slices
66// - Sigma Pedestal / slice = Sigma Pedestal / Sqrt (Number slices)
67// - Sigma Pedestal Error / slice = Sigma Pedestal Error / Sqrt (Number slices)
68//
69// The class also fills arrays with the signal vs. event number, creates a fourier
70// spectrum (see MHGausEvents::CreateFourierSpectrum()) and investigates if the
71// projected fourier components follow an exponential distribution.
72//
73// This same procedure is performed for the average pixels.
74//
75// The following results are written into MPedestalCam:
76//
77// - MCalibrationPix::SetMean()
78// - MCalibrationPix::SetMeanErr()
79// - MCalibrationPix::SetSigma()
80// - MCalibrationPix::SetSigmaErr()
81// - MCalibrationPix::SetProb()
82// - MCalibrationPix::SetNumPickup()
83//
84// For all averaged areas, the fitted sigma is multiplied with the square root of
85// the number involved pixels in order to be able to compare it to the average of
86// sigmas in the camera.
87//
88/////////////////////////////////////////////////////////////////////////////
89#include "MHPedestalCam.h"
90#include "MHPedestalPix.h"
91
92#include "MLog.h"
93#include "MLogManip.h"
94
95#include "MParList.h"
96
97#include "MExtractedSignalCam.h"
98#include "MExtractedSignalPix.h"
99
100#include "MPedestalCam.h"
101#include "MPedestalPix.h"
102
103#include "MCalibrationPix.h"
104
105#include "MBadPixelsCam.h"
106
107#include "MGeomCam.h"
108#include "MGeomPix.h"
109
110#include <TH1F.h>
111#include <TOrdCollection.h>
112
113ClassImp(MHPedestalCam);
114
115using namespace std;
116
117const Int_t MHPedestalCam::fgNbins = 50;
118const Axis_t MHPedestalCam::fgFirst = -57.5;
119const Axis_t MHPedestalCam::fgLast = 192.5;
120const TString MHPedestalCam::gsHistName = "Pedestal";
121const TString MHPedestalCam::gsHistTitle = "Pedestal";
122const TString MHPedestalCam::gsHistXTitle = "Charge [FADC slices]";
123const TString MHPedestalCam::gsHistYTitle = "Nr. events";
124const TString MHPedestalCam::fgNamePedestalCam = "MPedestalCam";
125
126// --------------------------------------------------------------------------
127//
128// Default constructor.
129//
130// Initializes:
131// - fExtractHiGainSlices to 0.
132// - the event frequency to 1200 Hz.
133// - the fRenormflag to kFALSE
134//
135// - fNbins to fgNbins
136// - fFirst to fgFirst
137// - fLast to fgLast
138//
139// - fHistName to gsHistName
140// - fHistTitle to gsHistTitle
141// - fHistXTitle to gsHistXTitle
142// - fHistYTitle to gsHistYTitle
143//
144MHPedestalCam::MHPedestalCam(const char *name, const char *title)
145 : fNamePedestalCamOut(fgNamePedestalCam), fNumEvents(0),
146 fExtractHiGainSlices(0.), fPedestalsOut(NULL)
147{
148
149 fName = name ? name : "MHPedestalCam";
150 fTitle = title ? title : "Class to fill the pedestal histograms";
151
152 SetPulserFrequency(1200);
153 SetRenorm(kFALSE);
154 SetLoGain(kFALSE);
155
156 SetBinning(fgNbins, fgFirst, fgLast);
157
158 SetHistName (gsHistName .Data());
159 SetHistTitle (gsHistTitle .Data());
160 SetHistXTitle(gsHistXTitle.Data());
161 SetHistYTitle(gsHistYTitle.Data());
162
163 SetFitStart();
164}
165
166// --------------------------------------------------------------------------
167//
168// Calls MHCalibrationCam::ResetHists()
169//
170// Resets:
171// - fSum
172// - fSumSquare
173// - fAreaSum
174// - fAreaSumSquare
175// - fAreaNum
176// - fSectorSum
177// - fSectorSumSquare
178// - fSectorNum
179//
180void MHPedestalCam::ResetHists()
181{
182
183 MHCalibrationCam::ResetHists();
184
185 fNumEvents = 0;
186
187 // If the size is yet set, set the size
188 if (fSum.GetSize()>0)
189 {
190 // Reset contents of arrays.
191 fSum.Reset();
192 fSumSquare.Reset();
193
194 fAreaSum. Reset();
195 fAreaSumSquare.Reset();
196 fAreaNum.Reset();
197
198 fSectorSum. Reset();
199 fSectorSumSquare.Reset();
200 fSectorNum.Reset();
201 }
202}
203// --------------------------------------------------------------------------
204//
205// Creates new MHCalibrationChargeCam only with the averaged areas:
206// the rest has to be retrieved directly, e.g. via:
207// MHPedestalCam *cam = MParList::FindObject("MHPedestalCam");
208// - cam->GetAverageSector(5).DrawClone();
209// - (*cam)[100].DrawClone()
210//
211TObject *MHPedestalCam::Clone(const char *) const
212{
213
214 MHPedestalCam *cam = new MHPedestalCam();
215
216 //
217 // Copy the data members
218 //
219 cam->fRunNumbers = fRunNumbers;
220 cam->fPulserFrequency = fPulserFrequency;
221 cam->fFlags = fFlags;
222 cam->fNbins = fNbins;
223 cam->fFirst = fFirst;
224 cam->fLast = fLast;
225
226 if (!IsAverageing())
227 return cam;
228
229 const Int_t navhi = fAverageHiGainAreas->GetSize();
230
231 for (int i=0; i<navhi; i++)
232 cam->fAverageHiGainAreas->AddAt(GetAverageHiGainArea(i).Clone(),i);
233
234 return cam;
235}
236
237// --------------------------------------------------------------------------
238//
239// Searches pointer to:
240// - MPedestalCam
241// - MExtractedSignalCam
242//
243// Searches or creates:
244// - MCalibrationPedCam
245//
246// Retrieves from MExtractedSignalCam:
247// - number of used High Gain FADC slices (MExtractedSignalCam::GetNumUsedHiGainFADCSlices())
248//
249// Initializes, if empty to MGeomCam::GetNumPixels():
250// - MHCalibrationCam::fHiGainArray
251//
252// Initializes, if empty to MGeomCam::GetNumAreas() for:
253// - MHCalibrationCam::fAverageHiGainAreas
254//
255// Initializes, if empty to MGeomCam::GetNumSectors() for:
256// - MHCalibrationCam::fAverageHiGainSectors
257//
258// Calls MHCalibrationCam::InitPedHists() for every entry in:
259// - MHCalibrationCam::fHiGainArray
260// - MHCalibrationCam::fAverageHiGainAreas
261// - MHCalibrationCam::fAverageHiGainSectors
262//
263// Sets Titles and Names for the Histograms
264// - MHCalibrationCam::fAverageHiGainAreas
265//
266Bool_t MHPedestalCam::ReInitHists(MParList *pList)
267{
268
269 MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam");
270 if (!signal && fRenorm)
271 {
272 *fLog << err << "Cannot find MExtractedSignalCam, but re-normalization switched on..."
273 << " Cannot find number of used slices ...abort." << endl;
274 return kFALSE;
275 }
276
277
278 if (!fPedestalsOut)
279 fPedestalsOut = (MPedestalCam*)pList->FindObject(AddSerialNumber(fNamePedestalCamOut),"MPedestalCam");
280
281 if (!fPedestalsOut)
282 {
283 *fLog << err << "Cannot find nor create outgoing MPedestalCam ..." << endl;
284 return kFALSE;
285 }
286
287 const Int_t npixels = fGeom->GetNumPixels();
288 const Int_t nsectors = fGeom->GetNumSectors();
289 const Int_t nareas = fGeom->GetNumAreas();
290
291 fCam = (MCalibrationCam*)pList->FindObject("MCalibrationPedCam", "MCalibrationCam");
292 if (!fCam)
293 {
294 fCam = (MCalibrationCam*)pList->FindCreateObj("MCalibrationCam", AddSerialNumber("MCalibrationPedCam"));
295 if (!fCam)
296 return kFALSE;
297
298 fCam->Init(*fGeom);
299 }
300
301 InitHiGainArrays(npixels,nareas,nsectors);
302
303 if (fSum.GetSize()==0)
304 {
305 fSum. Set(npixels);
306 fSumSquare. Set(npixels);
307 fAreaSum. Set(nareas);
308 fAreaSumSquare. Set(nareas);
309 fAreaNum. Set(nareas);
310 fSectorSum. Set(nsectors);
311 fSectorSumSquare.Set(nsectors);
312 fSectorNum.Set(nsectors);
313 }
314
315 fNumEvents = 0;
316
317 if (!signal)
318 return kTRUE;
319
320
321 Float_t sliceshi = signal->GetNumUsedHiGainFADCSlices();
322
323 if (sliceshi == 0.)
324 {
325 *fLog << err << "Number of used signal slices in MExtractedSignalCam is zero ... abort."
326 << endl;
327 return kFALSE;
328 }
329
330 if (fExtractHiGainSlices != 0. && sliceshi != fExtractHiGainSlices )
331 {
332 *fLog << err << "Number of used High Gain signal slices changed in MExtractedSignalCam ... abort."
333 << endl;
334 return kFALSE;
335 }
336
337 fExtractHiGainSlices = sliceshi;
338 *fLog << endl;
339 *fLog << inf << GetDescriptor()
340 << ": Number of found High Gain signal slices: " << fExtractHiGainSlices << endl;
341
342 return kTRUE;
343}
344
345// --------------------------------------------------------------------------
346//
347// Initializes the High Gain Arrays:
348//
349// - For every entry in the expanded arrays:
350// * Initialize an MHCalibrationChargePix
351// * Set Binning from fNbins, fFirst and fLast
352// * Set Binning of Abs Times histogram from fAbsNbins, fAbsFirst and fAbsLast
353// * Set Histgram names and titles from fHistName and fHistTitle
354// * Set Abs Times Histgram names and titles from fAbsHistName and fAbsHistTitle
355// * Set X-axis and Y-axis titles from fHistXTitle and fHistYTitle
356// * Set X-axis and Y-axis titles of Abs Times Histogram from fAbsHistXTitle and fAbsHistYTitle
357// * Call InitHists
358//
359//
360void MHPedestalCam::InitHiGainArrays(const Int_t npixels, const Int_t nareas, const Int_t nsectors)
361{
362
363 if (fHiGainArray->GetSize()==0)
364 {
365 for (Int_t i=0; i<npixels; i++)
366 {
367 fHiGainArray->AddAt(new MHPedestalPix(Form("%sHiGainPix%04d",fHistName.Data(),i),
368 Form("%s High Gain Pixel%04d",fHistTitle.Data(),i)),i);
369
370 MHPedestalPix &pix = (MHPedestalPix&)(*this)[i];
371 pix.SetBinning(fNbins, fFirst, fLast);
372 pix.SetProbLimit(fProbLimit);
373
374 InitHists(pix, (*fBadPixels)[i], i);
375 }
376 }
377
378 if (!IsAverageing())
379 return;
380
381 if (fAverageHiGainAreas->GetSize()==0)
382 {
383 for (Int_t j=0; j<nareas; j++)
384 {
385 fAverageHiGainAreas->AddAt(new MHPedestalPix(Form("%sHiGainArea%d",fHistName.Data(),j),
386 Form("%s High Gain Area Idx %d",fHistTitle.Data(),j)),j);
387
388 MHPedestalPix &pix = (MHPedestalPix&)GetAverageHiGainArea(j);
389
390 pix.SetBinning(fNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas),
391 fFirst, fLast);
392
393 InitHists(pix, fCam->GetAverageBadArea(j), j);
394
395 }
396 }
397
398 if (fAverageHiGainSectors->GetSize()==0)
399 {
400 for (Int_t j=0; j<nsectors; j++)
401 {
402 fAverageHiGainSectors->AddAt(new MHPedestalPix(Form("%sHiGainSector%02d",fHistName.Data(),j),
403 Form("%s High Gain Sector %02d",fHistTitle.Data(),j)),j);
404
405 MHPedestalPix &pix = (MHPedestalPix&)GetAverageHiGainSector(j);
406
407 pix.SetBinning(fNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas),
408 fFirst, fLast);
409
410 InitHists(pix, fCam->GetAverageBadSector(j), j);
411
412 }
413 }
414}
415
416// -------------------------------------------------------------------------------
417//
418// Retrieves pointer to MExtractedSignalCam:
419//
420// Retrieves from MGeomCam:
421// - number of pixels
422// - number of pixel areas
423// - number of sectors
424//
425// Fills HiGain histograms (MHGausEvents::FillHistAndArray()), respectively
426// with the signals MExtractedSignalCam::GetExtractedSignalHiGain().
427//
428Bool_t MHPedestalCam::FillHists(const MParContainer *par, const Stat_t w)
429{
430
431 MPedestalCam *pedestal = (MPedestalCam*)par;
432 if (!pedestal)
433 {
434 *fLog << err << "No argument in MHPedestalCam::Fill... abort." << endl;
435 return kFALSE;
436 }
437
438 const UInt_t npixels = fGeom->GetNumPixels();
439 const UInt_t nareas = fGeom->GetNumAreas();
440 const UInt_t nsectors = fGeom->GetNumSectors();
441
442 MArrayF sumareahi(nareas);
443 MArrayF sumsectorhi(nsectors);
444 MArrayI numareahi(nareas);
445 MArrayI numsectorhi(nsectors);
446
447 for (UInt_t i=0; i<npixels; i++)
448 {
449
450 MHCalibrationPix &histhi = (*this)[i];
451
452 if (histhi.IsExcluded())
453 continue;
454
455 const MPedestalPix &ped = (*pedestal)[i];
456
457 const Float_t pedes = ped.GetPedestal();
458
459 const Int_t aidx = (*fGeom)[i].GetAidx();
460 const Int_t sector = (*fGeom)[i].GetSector();
461
462 histhi.FillHistAndArray(pedes) ;
463
464 fSum[i] += pedes;
465 fAreaSum[aidx] += pedes;
466 fSectorSum[sector] += pedes;
467
468 const Float_t sqrsum = pedes*pedes;
469 fSumSquare[i] += sqrsum;
470 fAreaSumSquare[aidx] += sqrsum;
471 fSectorSumSquare[sector] += sqrsum;
472
473 sumareahi [aidx] += pedes;
474 numareahi [aidx] ++;
475 sumsectorhi[sector] += pedes;
476 numsectorhi[sector] ++;
477
478 }
479
480 for (UInt_t j=0; j<nareas; j++)
481 {
482 MHCalibrationPix &histhi = GetAverageHiGainArea(j);
483 histhi.FillHistAndArray(numareahi[j] == 0 ? 0. : sumareahi[j]/numareahi[j]);
484 }
485
486 for (UInt_t j=0; j<nsectors; j++)
487 {
488 MHCalibrationPix &histhi = GetAverageHiGainSector(j);
489 histhi.FillHistAndArray(numsectorhi[j] == 0 ? 0. : sumsectorhi[j]/numsectorhi[j]);
490
491 }
492
493 fNumEvents++;
494
495 return kTRUE;
496}
497
498// --------------------------------------------------------------------------
499//
500// Calls:
501// - MHCalibrationCam::FitHiGainArrays() with Bad Pixels flags 0
502//
503Bool_t MHPedestalCam::FinalizeHists()
504{
505
506 if (fNumEvents <= 1)
507 {
508 *fLog << err << GetDescriptor()
509 << ": Number of processed events smaller or equal to 1" << endl;
510 return kFALSE;
511 }
512
513 FitHists();
514
515 if (fRenorm)
516 RenormResults();
517
518 return kTRUE;
519
520}
521
522void MHPedestalCam::FitHists()
523{
524
525 for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
526 {
527
528 MHPedestalPix &hist = (MHPedestalPix&)(*this)[i];
529 MCalibrationPix &pix = (*fCam)[i];
530
531 //
532 // 1) Store calculated means and variances in the low-gain slices
533 //
534 pix.SetLoGainMean ( fSum[i] / fNumEvents );
535 const Double_t diff = fSumSquare[i] - fSum[i]*fSum[i]/fNumEvents;
536 pix.SetLoGainSigma ( diff > 0. ? TMath::Sqrt( diff / (fNumEvents-1) ) : 0.);
537 pix.SetLoGainMeanVar ( pix.GetLoGainSigma() * pix.GetLoGainSigma() / fNumEvents );
538 pix.SetLoGainSigmaVar( pix.GetLoGainMeanVar() / 4. );
539
540 if (hist.IsEmpty() || hist.IsOnlyOverflow() || hist.IsOnlyUnderflow())
541 continue;
542
543 //
544 // 2) Fit the Hi Gain histograms with a Gaussian
545 //
546 // if (i%100)
547 // hist.FitTripleGaus();
548 // else
549 TH1F *gaush = hist.GetHGausHist();
550 hist.FitGaus("RQ0",fFitStart,gaush->GetBinCenter(gaush->GetXaxis()->GetLast()));
551 //
552 // 4) Check for oscillations
553 //
554 if (IsOscillations())
555 hist.CreateFourierSpectrum();
556 //
557 // 5) Retrieve the results and store them in this class
558 //
559 pix.SetHiGainMean ( hist.GetMean() );
560 pix.SetHiGainMeanVar ( hist.GetMeanErr() * hist.GetMeanErr() );
561 pix.SetHiGainRms ( hist.GetHistRms() );
562 pix.SetHiGainSigma ( hist.GetSigma() );
563 pix.SetHiGainSigmaVar ( hist.GetSigmaErr()* hist.GetSigmaErr() );
564 pix.SetHiGainProb ( hist.GetProb() );
565 pix.SetHiGainNumBlackout( hist.GetBlackout() );
566 pix.SetHiGainNumPickup ( hist.GetPickup() );
567 }
568
569 if (!IsAverageing())
570 return;
571
572 for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++)
573 {
574
575 MHCalibrationPix &hist = GetAverageHiGainArea(j);
576 MCalibrationPix &pix = fCam->GetAverageArea(j);
577
578 if (hist.IsEmpty() || hist.IsOnlyOverflow() || hist.IsOnlyUnderflow())
579 continue;
580
581 //
582 // 2) Fit the Hi Gain histograms with a Gaussian
583 //
584 TH1F *gaush = hist.GetHGausHist();
585 hist.FitGaus("RQ0",fFitStart,gaush->GetBinCenter(gaush->GetXaxis()->GetLast()));
586 //
587 // 4) Check for oscillations
588 //
589 if (IsOscillations())
590 hist.CreateFourierSpectrum();
591 //
592 // 5) Retrieve the results and store them in this class
593 //
594 pix.SetHiGainMean ( hist.GetMean() );
595 pix.SetHiGainMeanVar ( hist.GetMeanErr() * hist.GetMeanErr() );
596 pix.SetHiGainRms ( hist.GetHistRms() );
597 pix.SetHiGainSigma ( hist.GetSigma() );
598 pix.SetHiGainSigmaVar ( hist.GetSigmaErr()* hist.GetSigmaErr() );
599 pix.SetHiGainProb ( hist.GetProb() );
600 pix.SetHiGainNumBlackout( hist.GetBlackout() );
601 pix.SetHiGainNumPickup ( hist.GetPickup() );
602 //
603 // 6) Store calculated means and variances in the low-gain slices
604 //
605 const ULong_t aevts = fNumEvents * fAreaNum[j];
606 if (aevts <= 1)
607 continue;
608
609 pix.SetLoGainMean ( fAreaSum[j] / aevts );
610 const Double_t diff = fAreaSumSquare[j] - fAreaSum[j]*fAreaSum[j]/aevts ;
611 pix.SetLoGainSigma( diff > 0. ? TMath::Sqrt( diff / (aevts-1) ) : 0.);
612 }
613
614 for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
615 {
616 MHCalibrationPix &hist = GetAverageHiGainSector(j);
617 MCalibrationPix &pix = fCam->GetAverageSector(j);
618
619 if (hist.IsEmpty() || hist.IsOnlyOverflow() || hist.IsOnlyUnderflow())
620 continue;
621
622 //
623 // 2) Fit the Hi Gain histograms with a Gaussian
624 //
625 TH1F *gaush = hist.GetHGausHist();
626 hist.FitGaus("RQ0",fFitStart,gaush->GetBinCenter(gaush->GetXaxis()->GetLast()));
627 //
628 // 4) Check for oscillations
629 //
630 if (IsOscillations())
631 hist.CreateFourierSpectrum();
632 //
633 // 5) Retrieve the results and store them in this class
634 //
635 pix.SetHiGainMean ( hist.GetMean() );
636 pix.SetHiGainMeanVar ( hist.GetMeanErr() * hist.GetMeanErr() );
637 pix.SetHiGainRms ( hist.GetHistRms() );
638 pix.SetHiGainSigma ( hist.GetSigma() );
639 pix.SetHiGainSigmaVar ( hist.GetSigmaErr()* hist.GetSigmaErr() );
640 pix.SetHiGainProb ( hist.GetProb() );
641 pix.SetHiGainNumBlackout( hist.GetBlackout() );
642 pix.SetHiGainNumPickup ( hist.GetPickup() );
643 //
644 // 6) Store calculated means and variances in the low-gain slices
645 //
646 const ULong_t sevts = fNumEvents * fSectorNum[j];
647 if (sevts <= 1)
648 continue;
649
650 pix.SetLoGainMean ( fSectorSum[j] / sevts );
651 const Double_t diff = fSectorSumSquare[j] - fSectorSum[j]*fSectorSum[j]/sevts ;
652 pix.SetLoGainSigma( diff > 0. ? TMath::Sqrt( diff / (sevts-1) ) : 0.);
653 }
654}
655
656// --------------------------------------------------------------------------
657//
658// Renormalizes the pedestal fit results by the following formulae:
659//
660// - Mean Pedestal / slice -> Mean Pedestal / Number slices
661// - Mean Pedestal Error / slice -> Mean Pedestal Error / Number slices
662// - Sigma Pedestal / slice -> Sigma Pedestal / Sqrt (Number slices)
663// - Sigma Pedestal Error / slice -> Sigma Pedestal Error / Sqrt (Number slices)
664//
665void MHPedestalCam::RenormResults()
666{
667
668 //
669 // One never knows...
670 //
671 if (fExtractHiGainSlices <= 0)
672 return;
673
674 const Float_t sqslices = TMath::Sqrt(fExtractHiGainSlices);
675
676 for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
677 {
678
679 MCalibrationPix &pix = (*fCam)[i];
680 MPedestalPix &ped = (*fPedestalsOut)[i];
681 pix.SetHiGainMean ( pix.GetHiGainMean() / fExtractHiGainSlices );
682 pix.SetLoGainMean ( pix.GetLoGainMean() / fExtractHiGainSlices );
683
684 ped.SetPedestal(pix.GetHiGainMean());
685 //
686 // Mean error goes with PedestalRMS/Sqrt(entries) -> scale with sqrt(slices)
687 //
688 pix.SetHiGainMeanVar ( pix.GetHiGainMeanVar() / fExtractHiGainSlices );
689 pix.SetLoGainMeanVar ( pix.GetHiGainMeanVar() / fExtractHiGainSlices );
690
691 //
692 // Sigma goes like PedestalRMS -> scale with sqrt(slices)
693 //
694 pix.SetHiGainSigma ( pix.GetHiGainSigma() / sqslices );
695 pix.SetLoGainSigma ( pix.GetLoGainSigma() / sqslices );
696
697 ped.SetPedestalRms(pix.GetHiGainSigma());
698 //
699 // Sigma error goes like PedestalRMS/2.(entries) -> scale with sqrt(slices)
700 //
701 pix.SetHiGainSigmaVar ( pix.GetHiGainSigmaVar() / fExtractHiGainSlices );
702 pix.SetLoGainSigmaVar ( pix.GetLoGainSigmaVar() / fExtractHiGainSlices );
703 }
704
705 if (!IsAverageing())
706 return;
707
708 for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++)
709 {
710
711 MCalibrationPix &pix = fCam->GetAverageArea(j);
712 pix.SetHiGainMean ( pix.GetHiGainMean() / fExtractHiGainSlices );
713 pix.SetLoGainMean ( pix.GetLoGainMean() / fExtractHiGainSlices );
714 pix.SetHiGainMeanVar ( pix.GetHiGainMeanVar() / fExtractHiGainSlices );
715 pix.SetHiGainSigma ( pix.GetHiGainSigma() / sqslices );
716 pix.SetLoGainSigma ( pix.GetLoGainSigma() / sqslices );
717 pix.SetHiGainSigmaVar ( pix.GetHiGainSigmaVar() / fExtractHiGainSlices );
718 }
719
720 for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
721 {
722 MCalibrationPix &pix = fCam->GetAverageSector(j);
723 pix.SetHiGainMean ( pix.GetHiGainMean() / fExtractHiGainSlices );
724 pix.SetLoGainMean ( pix.GetLoGainMean() / fExtractHiGainSlices );
725 pix.SetHiGainMeanVar ( pix.GetHiGainMeanVar() / fExtractHiGainSlices );
726 pix.SetHiGainSigma ( pix.GetHiGainSigma() / sqslices );
727 pix.SetLoGainSigma ( pix.GetLoGainSigma() / sqslices );
728 pix.SetHiGainSigmaVar ( pix.GetHiGainSigmaVar() / fExtractHiGainSlices );
729 }
730
731}
732
733
734// ------------------------------------------------------------------
735//
736// The types are as follows:
737//
738// Fitted values:
739// ==============
740//
741// 0: Fitted Charge
742// 1: Error of fitted Charge
743// 2: Sigma of fitted Charge
744// 3: Error of Sigma of fitted Charge
745//
746//
747// Useful variables derived from the fit results:
748// =============================================
749//
750// 4: Returned probability of Gauss fit to Charge distribution
751// 5: Relative differenc of calculated pedestal (per slice) and fitted (per slice)
752// 6: Error of the Relative differenc of calculated pedestal (per slice) and fitted (per slice)
753// 7: Relative difference of the error of the mean pedestal (per slice) - calculated and fitted
754// 8: Relative differenc of calculated pedestal RMS (per slice) and fitted sigma (per slice)
755// 9: Error of Relative differenc of calculated pedestal RMS (per slice) and fitted sigma (per slice)
756// 10: Relative difference of the error of the pedestal RMS (per slice) - calculated and fitted
757//
758// Localized defects:
759// ==================
760//
761// 11: Gaus fit not OK
762// 12: Fourier spectrum not OK
763//
764Bool_t MHPedestalCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
765{
766
767 if (fHiGainArray->GetSize() <= idx)
768 return kFALSE;
769
770 if ((*this)[idx].IsExcluded())
771 return kFALSE;
772
773 const MPedestalPix &ppix = (*fPedestalsOut)[idx];
774
775 const Float_t ped = ppix.GetPedestal();
776 const Float_t rms = ppix.GetPedestalRms();
777
778 const Float_t entsqr = TMath::Sqrt((Float_t)fPedestalsOut->GetNumSlices()*ppix.GetNumEvents());
779
780 const Float_t pederr = rms/entsqr;
781 const Float_t rmserr = rms/entsqr/2.;
782
783 const Float_t mean = (*this)[idx].GetMean();
784 const Float_t meanerr = (*this)[idx].GetMeanErr();
785 const Float_t sigma = (*this)[idx].GetSigma() ;
786 const Float_t sigmaerr = (*this)[idx].GetSigmaErr();
787
788 switch (type)
789 {
790 case 0:
791 val = mean;
792 break;
793 case 1:
794 val = meanerr;
795 break;
796 case 2:
797 val = sigma;
798 break;
799 case 3:
800 val = sigmaerr;
801 break;
802 case 4:
803 val = (*this)[idx].GetProb();
804 break;
805 case 5:
806 val = 2.*(mean-ped)/(ped+mean);
807 break;
808 case 6:
809 val = TMath::Sqrt((pederr*pederr + meanerr*meanerr) * (ped*ped + mean*mean))
810 *2./(ped+mean)/(ped+mean);
811 break;
812 case 7:
813 val = 2.*(meanerr-pederr)/(pederr + meanerr);
814 break;
815 case 8:
816 val = 2.*(sigma-rms)/(sigma+rms);
817 break;
818 case 9:
819 val = TMath::Sqrt((rmserr*rmserr + sigmaerr*sigmaerr) * (rms*rms + sigma*sigma))
820 *2./(rms+sigma)/(rms+sigma);
821 break;
822 case 10:
823 val = 2.*(sigmaerr - rmserr)/(sigmaerr + rmserr);
824 break;
825 case 11:
826 if (!(*this)[idx].IsGausFitOK())
827 val = 1.;
828 break;
829 case 12:
830 if (!(*this)[idx].IsFourierSpectrumOK())
831 val = 1.;
832 break;
833 default:
834 return kFALSE;
835 }
836 return kTRUE;
837}
838
839// --------------------------------------------------------------------------
840//
841// Calls MHGausEvents::DrawClone() for pixel idx
842//
843void MHPedestalCam::DrawPixelContent(Int_t idx) const
844{
845 (*this)[idx].DrawClone();
846}
Note: See TracBrowser for help on using the repository browser.