Context Navigation

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

Visit:

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: