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

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