Context Navigation

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

Visit:

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

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

Download in other formats: