Context Navigation

source: trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc@ 5365

Visit:

Last change on this file since 5365 was 5357, checked in by gaug, 21 years ago
* empty log message *
File size: 22.8 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): Josep Flix 04/2001 <mailto:jflix@ifae.es>
19	! Author(s): Thomas Bretz 05/2001 <mailto:tbretz@astro.uni-wuerzburg.de>
20	! Author(s): Sebastian Commichau 12/2003
21	! Author(s): Javier Rico 01/2004 <mailto:jrico@ifae.es>
22	! Author(s): Markus Gaug 01/2004 <mailto:markus@ifae.es>
23	!
24	! Copyright: MAGIC Software Development, 2000-2004
25	!
26	!
27	\* ======================================================================== */
28	/////////////////////////////////////////////////////////////////////////////
29	//
30	// MPedCalcPedRun
31	//
32	// This task takes a pedestal run file and fills MPedestalCam during
33	// the Process() with the pedestal and rms computed in an event basis.
34	// In the PostProcess() MPedestalCam is finally filled with the pedestal
35	// mean and rms computed in a run basis.
36	// More than one run (file) can be merged
37	//
38	// MPedCalcPedRun applies the following formula (1):
39	//
40	// Pedestal per slice = sum(x_i) / n / slices
41	// PedRMS per slice = Sqrt( ( sum(x_i^2) - sum(x_i)^2/n ) / n-1 / slices )
42	//
43	// where x_i is the sum of "slices" FADC slices and sum means the sum over all
44	// events. "n" is the number of events, "slices" is the number of summed FADC samples.
45	//
46	// Note that the slice-to-slice fluctuations are not Gaussian, but Poissonian, thus
47	// asymmetric and they are correlated.
48	//
49	// It is important to know that the Pedestal per slice and PedRMS per slice depend
50	// on the number of used FADC slices, as seen in the following plots:
51	//
52	//Begin_Html
53	/*
54	<img src="images/PedestalStudyInner.gif">
55	*/
56	//End_Html
57	//
58	//Begin_Html
59	/*
60	<img src="images/PedestalStudyOuter.gif">
61	*/
62	//
63	// The plots show the inner and outer pixels, respectivly and have the following meaning:
64	//
65	// 1) The calculated mean pedestal per slice (from MPedCalcPedRun)
66	// 2) The fitted mean pedestal per slice (from MHPedestalCam)
67	// 3) The calculated pedestal RMS per slice (from MPedCalcPedRun)
68	// 4) The fitted sigma of the pedestal distribution per slice
69	// (from MHPedestalCam)
70	// 5) The relative difference between calculation and histogram fit
71	// for the mean
72	// 6) The relative difference between calculation and histogram fit
73	// for the sigma or RMS, respectively.
74	//
75	// The calculated means do not change significantly except for the case of 2 slices,
76	// however the RMS changes from 5.7 per slice in the case of 2 extracted slices
77	// to 8.3 per slice in the case of 26 extracted slices. This change is very significant.
78	//
79	// The plots have been produced on run 20123. You can reproduce them using
80	// the macro pedestalstudies.C
81	//
82	// Usage of this class:
83	// ====================
84	//
85	// Call: SetRange(higainfirst, higainlast, logainfirst, logainlast)
86	// to modify the ranges in which the window is allowed to move.
87	// Defaults are:
88	//
89	// fHiGainFirst = fgHiGainFirst = 0
90	// fHiGainLast = fgHiGainLast = 29
91	// fLoGainFirst = fgLoGainFirst = 0
92	// fLoGainLast = fgLoGainLast = 14
93	//
94	// Call: SetWindowSize(windowhigain, windowlogain)
95	// to modify the sliding window widths. Windows have to be an even number.
96	// In case of odd numbers, the window will be modified.
97	//
98	// Defaults are:
99	//
100	// fHiGainWindowSize = fgHiGainWindowSize = 14
101	// fLoGainWindowSize = fgLoGainWindowSize = 0
102	//
103	//
104	// ToDo:
105	// - Take a setup file (ReadEnv-implementation) as input
106	//
107	//
108	// Input Containers:
109	// MRawEvtData
110	// MRawRunHeader
111	// MGeomCam
112	//
113	// Output Containers:
114	// MPedestalCam
115	//
116	// See also: MPedestalCam, MPedestalPix, MHPedestalCam, MExtractor
117	//
118	/////////////////////////////////////////////////////////////////////////////
119	#include "MPedCalcPedRun.h"
120	#include "MExtractor.h"
121
122	#include "MParList.h"
123
124	#include "MLog.h"
125	#include "MLogManip.h"
126
127	#include "MRawRunHeader.h"
128	#include "MRawEvtPixelIter.h"
129	#include "MRawEvtData.h"
130
131	#include "MPedestalPix.h"
132	#include "MPedestalCam.h"
133
134	#include "MGeomPix.h"
135	#include "MGeomCam.h"
136
137	ClassImp(MPedCalcPedRun);
138
139	using namespace std;
140
141	const Byte_t MPedCalcPedRun::fgHiGainFirst = 0;
142	const Byte_t MPedCalcPedRun::fgHiGainLast = 29;
143	const Byte_t MPedCalcPedRun::fgLoGainFirst = 0;
144	const Byte_t MPedCalcPedRun::fgLoGainLast = 14;
145	const Byte_t MPedCalcPedRun::fgHiGainWindowSize = 14;
146	const Byte_t MPedCalcPedRun::fgLoGainWindowSize = 0;
147	// --------------------------------------------------------------------------
148	//
149	// Default constructor:
150	//
151	// Sets:
152	// - all pointers to NULL
153	// - fWindowSizeHiGain to fgHiGainWindowSize
154	// - fWindowSizeLoGain to fgLoGainWindowSize
155	//
156	// Calls:
157	// - AddToBranchList("fHiGainPixId");
158	// - AddToBranchList("fHiGainFadcSamples");
159	// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
160	// - Clear()
161	//
162	MPedCalcPedRun::MPedCalcPedRun(const char name, const char title)
163	: fWindowSizeHiGain(fgHiGainWindowSize),
164	fWindowSizeLoGain(fgLoGainWindowSize),
165	fGeom(NULL)
166	{
167	fName = name ? name : "MPedCalcPedRun";
168	fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data";
169
170	AddToBranchList("fHiGainPixId");
171	AddToBranchList("fHiGainFadcSamples");
172
173	SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
174	Clear();
175	}
176
177	// --------------------------------------------------------------------------
178	//
179	// Sets:
180	// - fNumSamplesTot to 0
181	// - fRawEvt to NULL
182	// - fRunHeader to NULL
183	// - fPedestals to NULL
184	//
185	void MPedCalcPedRun::Clear(const Option_t *o)
186	{
187	fNumSamplesTot = 0;
188
189	fRawEvt = NULL;
190	fRunHeader = NULL;
191	fPedestals = NULL;
192
193	// If the size is yet set, set the size
194	if (fSumx.GetSize()>0)
195	{
196	// Reset contents of arrays.
197	fSumx.Reset();
198	fSumx2.Reset();
199	fSumAB0.Reset();
200	fSumAB1.Reset();
201
202	fAreaSumx. Reset();
203	fAreaSumx2.Reset();
204	fAreaSumAB0.Reset();
205	fAreaSumAB1.Reset();
206	fAreaValid.Reset();
207
208	fSectorSumx. Reset();
209	fSectorSumx2.Reset();
210	fSectorSumAB0.Reset();
211	fSectorSumAB1.Reset();
212	fSectorValid.Reset();
213	}
214	}
215
216	// --------------------------------------------------------------------------
217	//
218	// SetRange:
219	//
220	// Calls:
221	// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
222	// - SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
223	//
224	void MPedCalcPedRun::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
225	{
226	MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
227
228	//
229	// Redo the checks if the window is still inside the ranges
230	//
231	SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
232	}
233
234	// --------------------------------------------------------------------------
235	//
236	// Checks:
237	// - if a window is odd, subtract one
238	// - if a window is bigger than the one defined by the ranges, set it to the available range
239	// - if a window is smaller than 2, set it to 2
240	//
241	// Sets:
242	// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain
243	// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain
244	// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
245	// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)
246	//
247	void MPedCalcPedRun::SetWindowSize(Byte_t windowh, Byte_t windowl)
248	{
249
250	fWindowSizeHiGain = windowh & ~1;
251	fWindowSizeLoGain = windowl & ~1;
252
253	if (fWindowSizeHiGain != windowh)
254	*fLog << warn << GetDescriptor() << ": Hi Gain window size has to be even, set to: "
255	<< int(fWindowSizeHiGain) << " samples " << endl;
256
257	if (fWindowSizeLoGain != windowl)
258	*fLog << warn << GetDescriptor() << ": Lo Gain window size has to be even, set to: "
259	<< int(fWindowSizeLoGain) << " samples " << endl;
260
261	if (fWindowSizeHiGain == 0)
262	{
263	*fLog << warn;
264	*fLog << GetDescriptor() << ": HiGain window currently set to 0, will set it to 2 samples " << endl;
265	fWindowSizeHiGain = 2;
266	}
267
268	const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1;
269	const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1;
270
271	if (fWindowSizeHiGain > availhirange)
272	{
273	*fLog << warn << GetDescriptor()
274	<< Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",(int)fWindowSizeHiGain,
275	" is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl;
276	*fLog << warn << GetDescriptor()
277	<< ": Will set window size to: " << (int)availhirange << endl;
278	fWindowSizeHiGain = availhirange;
279	}
280
281	if (fWindowSizeLoGain > availlorange)
282	{
283	*fLog << warn << GetDescriptor()
284	<< Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",(int)fWindowSizeLoGain,
285	" is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl;
286	*fLog << warn << GetDescriptor()
287	<< ": Will set window size to: " << (int)availlorange << endl;
288	fWindowSizeLoGain = availlorange;
289	}
290
291
292	fNumHiGainSamples = (Float_t)fWindowSizeHiGain;
293	fNumLoGainSamples = (Float_t)fWindowSizeLoGain;
294
295	fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
296	fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
297	}
298
299	// --------------------------------------------------------------------------
300	//
301	// Look for the following input containers:
302	//
303	// - MRawEvtData
304	// - MRawRunHeader
305	// - MGeomCam
306	//
307	// The following output containers are also searched and created if
308	// they were not found:
309	//
310	// - MPedestalCam
311	//
312	Int_t MPedCalcPedRun::PreProcess( MParList *pList )
313	{
314	Clear();
315
316	fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
317	if (!fRawEvt)
318	{
319	*fLog << err << "MRawEvtData not found... aborting." << endl;
320	return kFALSE;
321	}
322
323	fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
324	if (!fRunHeader)
325	{
326	*fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl;
327	return kFALSE;
328	}
329
330	fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
331	if (!fGeom)
332	{
333	*fLog << err << "MGeomCam not found... aborting." << endl;
334	return kFALSE;
335	}
336
337	fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam");
338	if (!fPedestals)
339	return kFALSE;
340
341	return kTRUE;
342	}
343
344	// --------------------------------------------------------------------------
345	//
346	// The ReInit searches for:
347	// - MRawRunHeader::GetNumSamplesHiGain()
348	// - MRawRunHeader::GetNumSamplesLoGain()
349	//
350	// In case that the variables fHiGainLast and fLoGainLast are smaller than
351	// the even part of the number of samples obtained from the run header, a
352	// warning is given an the range is set back accordingly. A call to:
353	// - SetRange(fHiGainFirst, fHiGainLast-diff, fLoGainFirst, fLoGainLast) or
354	// - SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff)
355	// is performed in that case. The variable diff means here the difference
356	// between the requested range (fHiGainLast) and the available one. Note that
357	// the functions SetRange() are mostly overloaded and perform more checks,
358	// modifying the ranges again, if necessary.
359	//
360	Bool_t MPedCalcPedRun::ReInit(MParList *pList)
361	{
362
363	Int_t lastdesired = (Int_t)fLoGainLast;
364	Int_t lastavailable = (Int_t)fRunHeader->GetNumSamplesLoGain()-1;
365
366	if (lastdesired > lastavailable)
367	{
368	const Int_t diff = lastdesired - lastavailable;
369	*fLog << endl;
370	*fLog << warn << GetDescriptor()
371	<< Form("%s%2i%s%2i%s%2i%s",": Selected Lo Gain FADC Window [",
372	(int)fLoGainFirst,",",lastdesired,
373	"] ranges out of the available limits: [0,",lastavailable,"].") << endl;
374	*fLog << GetDescriptor() << ": Will reduce the upper edge to " << (int)(fLoGainLast - diff) << endl;
375	SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff);
376	}
377
378	lastdesired = (Int_t)fHiGainLast;
379	lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
380
381	if (lastdesired > lastavailable)
382	{
383	const Int_t diff = lastdesired - lastavailable;
384	*fLog << endl;
385	*fLog << warn << GetDescriptor()
386	<< Form("%s%2i%s%2i%s%2i%s",": Selected Hi Gain Range [",
387	(int)fHiGainFirst,",",lastdesired,
388	"] ranges out of the available limits: [0,",lastavailable,"].") << endl;
389	*fLog << warn << GetDescriptor()
390	<< Form("%s%2i%s",": Will possibly use ",diff," samples from the Low-Gain for the High-Gain range")
391	<< endl;
392	fHiGainLast -= diff;
393	fHiLoLast = diff;
394	}
395
396	lastdesired = (Int_t)fHiGainFirst+fWindowSizeHiGain-1;
397	lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
398
399	if (lastdesired > lastavailable)
400	{
401	const Int_t diff = lastdesired - lastavailable;
402	*fLog << endl;
403	*fLog << warn << GetDescriptor()
404	<< Form("%s%2i%s%2i%s",": Selected Hi Gain FADC Window size ",
405	(int)fWindowSizeHiGain,
406	" ranges out of the available limits: [0,",lastavailable,"].") << endl;
407	*fLog << warn << GetDescriptor()
408	<< Form("%s%2i%s",": Will use ",diff," samples from the Low-Gain for the High-Gain extraction")
409	<< endl;
410
411	if ((Int_t)fWindowSizeHiGain > diff)
412	{
413	fWindowSizeHiGain -= diff;
414	fWindowSizeLoGain += diff;
415	}
416	else
417	{
418	fWindowSizeLoGain += fWindowSizeHiGain;
419	fLoGainFirst = diff-fWindowSizeHiGain;
420	fWindowSizeHiGain = 0;
421	}
422	}
423
424
425	const Int_t npixels = fPedestals->GetSize();
426	const Int_t areas = fPedestals->GetAverageAreas();
427	const Int_t sectors = fPedestals->GetAverageSectors();
428
429	if (fSumx.GetSize()==0)
430	{
431	fSumx. Set(npixels);
432	fSumx2. Set(npixels);
433	fSumAB0.Set(npixels);
434	fSumAB1.Set(npixels);
435
436	fAreaSumx. Set(areas);
437	fAreaSumx2. Set(areas);
438	fAreaSumAB0.Set(areas);
439	fAreaSumAB1.Set(areas);
440	fAreaValid. Set(areas);
441
442	fSectorSumx. Set(sectors);
443	fSectorSumx2. Set(sectors);
444	fSectorSumAB0.Set(sectors);
445	fSectorSumAB1.Set(sectors);
446	fSectorValid. Set(sectors);
447
448	fSumx. Reset();
449	fSumx2.Reset();
450	}
451
452	if (fWindowSizeHiGain == 0 && fWindowSizeLoGain == 0)
453	{
454	*fLog << err << GetDescriptor()
455	<< ": Number of extracted Slices is 0, abort ... " << endl;
456	return kFALSE;
457	}
458
459
460	*fLog << endl;
461	*fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeHiGain
462	<< " HiGain FADC samples starting with slice: " << (int)fHiGainFirst << endl;
463	*fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeLoGain
464	<< " LoGain FADC samples starting with slice: " << (int)fLoGainFirst << endl;
465
466	return kTRUE;
467	}
468
469	// --------------------------------------------------------------------------
470	//
471	// Fill the MPedestalCam container with the signal mean and rms for the event.
472	// Store the measured signal in arrays fSumx and fSumx2 so that we can
473	// calculate the overall mean and rms in the PostProcess()
474	//
475	Int_t MPedCalcPedRun::Process()
476	{
477
478	MRawEvtPixelIter pixel(fRawEvt);
479
480	while (pixel.Next())
481	{
482	const UInt_t idx = pixel.GetPixelId();
483	const UInt_t aidx = (*fGeom)[idx].GetAidx();
484	const UInt_t sector = (*fGeom)[idx].GetSector();
485
486	Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst;
487	Byte_t *end = ptr + fWindowSizeHiGain;
488
489	UInt_t sum = 0;
490	UInt_t sqr = 0;
491	UInt_t ab0 = 0;
492	UInt_t ab1 = 0;
493	Int_t cnt = 0;
494
495	if (fWindowSizeHiGain != 0)
496	{
497	do
498	{
499	sum += *ptr;
500	sqr += ptr *ptr;
501
502	if (pixel.IsABFlagValid())
503	{
504	const Int_t abFlag = (fHiGainFirst + pixel.HasABFlag() + cnt) & 0x1;
505	if (abFlag)
506	ab1 += *ptr;
507	else
508	ab0 += *ptr;
509
510	cnt++;
511	}
512	}
513	while (++ptr != end);
514	}
515
516	cnt = 0;
517
518	if (fWindowSizeLoGain != 0)
519	{
520
521	ptr = pixel.GetLoGainSamples() + fLoGainFirst;
522	end = ptr + fWindowSizeLoGain;
523
524	do
525	{
526	sum += *ptr;
527	sqr += ptr *ptr;
528
529	if (pixel.IsABFlagValid())
530	{
531	const Int_t abFlag = (fLoGainFirst + pixel.GetNumHiGainSamples() + pixel.HasABFlag() + cnt)
532	& 0x1;
533	if (abFlag)
534	ab1 += *ptr;
535	else
536	ab0 += *ptr;
537
538	cnt++;
539	}
540
541	}
542	while (++ptr != end);
543
544	}
545
546	const Float_t msum = (Float_t)sum;
547
548	//
549	// These three lines have been uncommented by Markus Gaug
550	// If anybody needs them, please contact me!!
551	//
552	// const Float_t higainped = msum/fNumHiGainSlices;
553	// const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSlices)/(fNumHiGainSlices-1.));
554	// (*fPedestals)[idx].Set(higainped, higainrms);
555
556	fSumx[idx] += msum;
557	fAreaSumx[aidx] += msum;
558	fSectorSumx[sector] += msum;
559
560	//
561	// The old version:
562	//
563	// const Float_t msqr = (Float_t)sqr;
564	// fSumx2[idx] += msqr;
565	//
566	// The new version:
567	//
568	const Float_t sqrsum = msum*msum;
569	fSumx2[idx] += sqrsum;
570	fAreaSumx2[aidx] += sqrsum;
571	fSectorSumx2[sector] += sqrsum;
572
573	//
574	// Now, the sums separated for AB0 and AB1
575	//
576	fSumAB0[idx] += ab0;
577	fSumAB1[idx] += ab1;
578
579	fAreaSumAB0[aidx] += ab0;
580	fAreaSumAB1[aidx] += ab1;
581
582	fSectorSumAB0[aidx] += ab0;
583	fSectorSumAB1[aidx] += ab1;
584	}
585
586	fPedestals->SetReadyToSave();
587	fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain;
588
589	return kTRUE;
590	}
591
592	// --------------------------------------------------------------------------
593	//
594	// Compute signal mean and rms in the whole run and store it in MPedestalCam
595	//
596	Int_t MPedCalcPedRun::PostProcess()
597	{
598
599	// Compute pedestals and rms from the whole run
600	const ULong_t n = fNumSamplesTot;
601	const ULong_t nevts = GetNumExecutions();
602
603	MRawEvtPixelIter pixel(fRawEvt);
604
605	while (pixel.Next())
606	{
607
608	const Int_t pixid = pixel.GetPixelId();
609	const UInt_t aidx = (*fGeom)[pixid].GetAidx();
610	const UInt_t sector = (*fGeom)[pixid].GetSector();
611
612	fAreaValid [aidx]++;
613	fSectorValid[sector]++;
614
615	const Float_t sum = fSumx.At(pixid);
616	const Float_t sum2 = fSumx2.At(pixid);
617	const Float_t higainped = sum/n;
618	//
619	// The old version:
620	//
621	// const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.));
622	//
623	// The new version:
624	//
625	// 1. Calculate the Variance of the sums:
626	Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.);
627	// 2. Scale the variance to the number of slices:
628	higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain);
629	// 3. Calculate the RMS from the Variance:
630	const Float_t rms = higainVar<0 ? 0. : TMath::Sqrt(higainVar);
631	// 4. Calculate the amplitude of the 150MHz "AB" noise
632	const Float_t abOffs = (fSumAB0.At(pixid) - fSumAB1.At(pixid)) / n;
633
634	(*fPedestals)[pixid].Set(higainped,rms,abOffs);
635	}
636
637	//
638	// Loop over the (two) area indices to get the averaged pedestal per aidx
639	//
640	for (Int_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
641	{
642
643	const Int_t napix = fAreaValid.At(aidx);
644
645	if (napix == 0)
646	continue;
647
648	const Float_t sum = fAreaSumx.At(aidx);
649	const Float_t sum2 = fAreaSumx2.At(aidx);
650	const ULong_t an = napix * n;
651	const ULong_t aevts = napix * nevts;
652
653	const Float_t higainped = sum/an;
654
655	// 1. Calculate the Variance of the sums:
656	Float_t higainVar = (sum2-sum*sum/aevts)/(aevts-1.);
657	// 2. Scale the variance to the number of slices:
658	higainVar /= fWindowSizeHiGain+fWindowSizeLoGain;
659	// 3. Calculate the RMS from the Variance:
660	Float_t higainrms = TMath::Sqrt(higainVar);
661	// 4. Re-scale it with the square root of the number of involved pixels
662	// in order to be comparable to the mean of pedRMS of that area
663	higainrms *= TMath::Sqrt((Float_t)napix);
664	// 5. Calculate the amplitude of the 150MHz "AB" noise
665	const Float_t abOffs = (fAreaSumAB0.At(aidx) - fAreaSumAB1.At(aidx)) / an;
666
667	fPedestals->GetAverageArea(aidx).Set(higainped, higainrms,abOffs);
668	}
669
670	//
671	// Loop over the (six) sector indices to get the averaged pedestal per sector
672	//
673	for (Int_t sector=0; sector<fSectorValid.GetSize(); sector++)
674	{
675
676	const Int_t nspix = fSectorValid.At(sector);
677
678	if (nspix == 0)
679	continue;
680
681	const Float_t sum = fSectorSumx.At(sector);
682	const Float_t sum2 = fSectorSumx2.At(sector);
683	const ULong_t sn = nspix * n;
684	const ULong_t sevts = nspix * nevts;
685
686	const Float_t higainped = sum/sn;
687
688	// 1. Calculate the Variance of the sums:
689	Float_t higainVar = (sum2-sum*sum/sevts)/(sevts-1.);
690	// 2. Scale the variance to the number of slices:
691	higainVar /= fWindowSizeHiGain+fWindowSizeLoGain;
692	// 3. Calculate the RMS from the Variance:
693	Float_t higainrms = TMath::Sqrt(higainVar);
694	// 4. Re-scale it with the square root of the number of involved pixels
695	// in order to be comparable to the mean of pedRMS of that sector
696	higainrms *= TMath::Sqrt((Float_t)nspix);
697	// 5. Calculate the amplitude of the 150MHz "AB" noise
698	const Float_t abOffs = (fSectorSumAB0.At(sector) - fSectorSumAB1.At(sector)) / sn;
699
700	fPedestals->GetAverageSector(sector).Set(higainped, higainrms, abOffs);
701	}
702
703	fPedestals->SetTotalEntries(fNumSamplesTot);
704	fPedestals->SetReadyToSave();
705
706	return kTRUE;
707	}
708
709	Int_t MPedCalcPedRun::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
710	{
711	if (MExtractor::ReadEnv(env, prefix, print)==kERROR)
712	return kERROR;
713
714	Byte_t hw = fWindowSizeHiGain;
715	Byte_t lw = fWindowSizeLoGain;
716
717	Bool_t rc = kFALSE;
718
719	if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print))
720	{
721	hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw);
722	rc=kTRUE;
723	}
724
725	if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print))
726	{
727	lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw);
728	rc=kTRUE;
729	}
730
731	if (rc)
732	SetWindowSize(hw, lw);
733
734	return rc;
735	}

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

Download in other formats: