Context Navigation

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

Visit:

Last change on this file since 8067 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
117	ClassImp(MHPedestalCam);
118
119	using namespace std;
120
121	const Int_t MHPedestalCam::fgNbins = 50;
122	const Axis_t MHPedestalCam::fgFirst = -57.5;
123	const Axis_t MHPedestalCam::fgLast = 192.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	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	//
185	void 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	//
216	TObject 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	//
271	Bool_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	//
365	void 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	//
440	Bool_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	//
515	Bool_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
534	void 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	//
677	void 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	//
776	Bool_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((pederrpederr + meanerrmeanerr) * (pedped + meanmean))
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((rmserrrmserr + sigmaerrsigmaerr) * (rmsrms + sigmasigma))
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	//
853	void MHPedestalCam::DrawPixelContent(Int_t idx) const
854	{
855	(*this)[idx].DrawClone();
856	}

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

Download in other formats: