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

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