Context Navigation

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

Visit:

Last change on this file since 5298 was 5220, checked in by gaug, 20 years ago
* empty log message *
File size: 23.1 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
262	if (fWindowSizeHiGain == 0)
263	{
264	*fLog << warn;
265	*fLog << GetDescriptor() << ": HiGain window currently set to 0, will set it to 2 samples " << endl;
266	fWindowSizeHiGain = 2;
267	}
268
269	if (fWindowSizeLoGain == 0)
270	{
271	*fLog << warn;
272	*fLog << GetDescriptor() << ": LoGain window currently set to 0, will set it to 2 samples " << endl;
273	fWindowSizeLoGain = 2;
274	}
275
276	const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1;
277	const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1;
278
279	if (fWindowSizeHiGain > availhirange)
280	{
281	*fLog << warn << GetDescriptor()
282	<< Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",(int)fWindowSizeHiGain,
283	" is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl;
284	*fLog << warn << GetDescriptor()
285	<< ": Will set window size to: " << (int)availhirange << endl;
286	fWindowSizeHiGain = availhirange;
287	}
288
289	if (fWindowSizeLoGain > availlorange)
290	{
291	*fLog << warn << GetDescriptor()
292	<< Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",(int)fWindowSizeLoGain,
293	" is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl;
294	*fLog << warn << GetDescriptor()
295	<< ": Will set window size to: " << (int)availlorange << endl;
296	fWindowSizeLoGain = availlorange;
297	}
298
299
300	fNumHiGainSamples = (Float_t)fWindowSizeHiGain;
301	fNumLoGainSamples = (Float_t)fWindowSizeLoGain;
302
303	fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
304	fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
305	}
306
307	// --------------------------------------------------------------------------
308	//
309	// Look for the following input containers:
310	//
311	// - MRawEvtData
312	// - MRawRunHeader
313	// - MGeomCam
314	//
315	// The following output containers are also searched and created if
316	// they were not found:
317	//
318	// - MPedestalCam
319	//
320	Int_t MPedCalcPedRun::PreProcess( MParList *pList )
321	{
322	Clear();
323
324	fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
325	if (!fRawEvt)
326	{
327	*fLog << err << "MRawEvtData not found... aborting." << endl;
328	return kFALSE;
329	}
330
331	fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
332	if (!fRunHeader)
333	{
334	*fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl;
335	return kFALSE;
336	}
337
338	fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
339	if (!fGeom)
340	{
341	*fLog << err << "MGeomCam not found... aborting." << endl;
342	return kFALSE;
343	}
344
345	fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam");
346	if (!fPedestals)
347	return kFALSE;
348
349	return kTRUE;
350	}
351
352	// --------------------------------------------------------------------------
353	//
354	// The ReInit searches for:
355	// - MRawRunHeader::GetNumSamplesHiGain()
356	// - MRawRunHeader::GetNumSamplesLoGain()
357	//
358	// In case that the variables fHiGainLast and fLoGainLast are smaller than
359	// the even part of the number of samples obtained from the run header, a
360	// warning is given an the range is set back accordingly. A call to:
361	// - SetRange(fHiGainFirst, fHiGainLast-diff, fLoGainFirst, fLoGainLast) or
362	// - SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff)
363	// is performed in that case. The variable diff means here the difference
364	// between the requested range (fHiGainLast) and the available one. Note that
365	// the functions SetRange() are mostly overloaded and perform more checks,
366	// modifying the ranges again, if necessary.
367	//
368	Bool_t MPedCalcPedRun::ReInit(MParList *pList)
369	{
370
371	Int_t lastdesired = (Int_t)fLoGainLast;
372	Int_t lastavailable = (Int_t)fRunHeader->GetNumSamplesLoGain()-1;
373
374	if (lastdesired > lastavailable)
375	{
376	const Int_t diff = lastdesired - lastavailable;
377	*fLog << endl;
378	*fLog << warn << GetDescriptor()
379	<< Form("%s%2i%s%2i%s%2i%s",": Selected Lo Gain FADC Window [",
380	(int)fLoGainFirst,",",lastdesired,
381	"] ranges out of the available limits: [0,",lastavailable,"].") << endl;
382	*fLog << GetDescriptor() << ": Will reduce the upper edge to " << (int)(fLoGainLast - diff) << endl;
383	SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff);
384	}
385
386	lastdesired = (Int_t)fHiGainLast;
387	lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
388
389	if (lastdesired > lastavailable)
390	{
391	const Int_t diff = lastdesired - lastavailable;
392	*fLog << endl;
393	*fLog << warn << GetDescriptor()
394	<< Form("%s%2i%s%2i%s%2i%s",": Selected Hi Gain Range [",
395	(int)fHiGainFirst,",",lastdesired,
396	"] ranges out of the available limits: [0,",lastavailable,"].") << endl;
397	*fLog << warn << GetDescriptor()
398	<< Form("%s%2i%s",": Will possibly use ",diff," samples from the Low-Gain for the High-Gain range")
399	<< endl;
400	fHiGainLast -= diff;
401	fHiLoLast = diff;
402	}
403
404	lastdesired = (Int_t)fHiGainFirst+fWindowSizeHiGain-1;
405	lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
406
407	if (lastdesired > lastavailable)
408	{
409	const Int_t diff = lastdesired - lastavailable;
410	*fLog << endl;
411	*fLog << warn << GetDescriptor()
412	<< Form("%s%2i%s%2i%s",": Selected Hi Gain FADC Window size ",
413	(int)fWindowSizeHiGain,
414	" ranges out of the available limits: [0,",lastavailable,"].") << endl;
415	*fLog << warn << GetDescriptor()
416	<< Form("%s%2i%s",": Will use ",diff," samples from the Low-Gain for the High-Gain extraction")
417	<< endl;
418
419	if ((Int_t)fWindowSizeHiGain > diff)
420	{
421	fWindowSizeHiGain -= diff;
422	fWindowSizeLoGain += diff;
423	}
424	else
425	{
426	fWindowSizeLoGain += fWindowSizeHiGain;
427	fLoGainFirst = diff-fWindowSizeHiGain;
428	fWindowSizeHiGain = 0;
429	}
430	}
431
432
433	const Int_t npixels = fPedestals->GetSize();
434	const Int_t areas = fPedestals->GetAverageAreas();
435	const Int_t sectors = fPedestals->GetAverageSectors();
436
437	if (fSumx.GetSize()==0)
438	{
439	fSumx. Set(npixels);
440	fSumx2. Set(npixels);
441	fSumAB0.Set(npixels);
442	fSumAB1.Set(npixels);
443
444	fAreaSumx. Set(areas);
445	fAreaSumx2. Set(areas);
446	fAreaSumAB0.Set(areas);
447	fAreaSumAB1.Set(areas);
448	fAreaValid. Set(areas);
449
450	fSectorSumx. Set(sectors);
451	fSectorSumx2. Set(sectors);
452	fSectorSumAB0.Set(sectors);
453	fSectorSumAB1.Set(sectors);
454	fSectorValid. Set(sectors);
455
456	fSumx. Reset();
457	fSumx2.Reset();
458	}
459
460	if (fWindowSizeHiGain == 0 && fWindowSizeLoGain == 0)
461	{
462	*fLog << err << GetDescriptor()
463	<< ": Number of extracted Slices is 0, abort ... " << endl;
464	return kFALSE;
465	}
466
467
468	*fLog << endl;
469	*fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeHiGain
470	<< " HiGain FADC samples starting with slice: " << (int)fHiGainFirst << endl;
471	*fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeLoGain
472	<< " LoGain FADC samples starting with slice: " << (int)fLoGainFirst << endl;
473
474	return kTRUE;
475	}
476
477	// --------------------------------------------------------------------------
478	//
479	// Fill the MPedestalCam container with the signal mean and rms for the event.
480	// Store the measured signal in arrays fSumx and fSumx2 so that we can
481	// calculate the overall mean and rms in the PostProcess()
482	//
483	Int_t MPedCalcPedRun::Process()
484	{
485
486	MRawEvtPixelIter pixel(fRawEvt);
487
488	while (pixel.Next())
489	{
490	const UInt_t idx = pixel.GetPixelId();
491	const UInt_t aidx = (*fGeom)[idx].GetAidx();
492	const UInt_t sector = (*fGeom)[idx].GetSector();
493
494	Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst;
495	Byte_t *end = ptr + fWindowSizeHiGain;
496
497	UInt_t sum = 0;
498	UInt_t sqr = 0;
499	UInt_t ab0 = 0;
500	UInt_t ab1 = 0;
501	Int_t cnt = 0;
502
503	if (fWindowSizeHiGain != 0)
504	{
505	do
506	{
507	sum += *ptr;
508	sqr += ptr *ptr;
509
510	if (pixel.IsABFlagValid())
511	{
512	const Int_t abFlag = (fHiGainFirst + pixel.HasABFlag() + cnt) & 0x1;
513	if (abFlag)
514	ab1 += *ptr;
515	else
516	ab0 += *ptr;
517
518	cnt++;
519	}
520	}
521	while (++ptr != end);
522	}
523
524	cnt = 0;
525
526	if (fWindowSizeLoGain != 0)
527	{
528
529	ptr = pixel.GetLoGainSamples() + fLoGainFirst;
530	end = ptr + fWindowSizeLoGain;
531
532	do
533	{
534	sum += *ptr;
535	sqr += ptr *ptr;
536
537	if (pixel.IsABFlagValid())
538	{
539	const Int_t abFlag = (fLoGainFirst + pixel.GetNumHiGainSamples() + pixel.HasABFlag() + cnt)
540	& 0x1;
541	if (abFlag)
542	ab1 += *ptr;
543	else
544	ab0 += *ptr;
545
546	cnt++;
547	}
548
549	}
550	while (++ptr != end);
551
552	}
553
554	const Float_t msum = (Float_t)sum;
555
556	//
557	// These three lines have been uncommented by Markus Gaug
558	// If anybody needs them, please contact me!!
559	//
560	// const Float_t higainped = msum/fNumHiGainSlices;
561	// const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSlices)/(fNumHiGainSlices-1.));
562	// (*fPedestals)[idx].Set(higainped, higainrms);
563
564	fSumx[idx] += msum;
565	fAreaSumx[aidx] += msum;
566	fSectorSumx[sector] += msum;
567
568	//
569	// The old version:
570	//
571	// const Float_t msqr = (Float_t)sqr;
572	// fSumx2[idx] += msqr;
573	//
574	// The new version:
575	//
576	const Float_t sqrsum = msum*msum;
577	fSumx2[idx] += sqrsum;
578	fAreaSumx2[aidx] += sqrsum;
579	fSectorSumx2[sector] += sqrsum;
580
581	//
582	// Now, the sums separated for AB0 and AB1
583	//
584	fSumAB0[idx] += ab0;
585	fSumAB1[idx] += ab1;
586
587	fAreaSumAB0[aidx] += ab0;
588	fAreaSumAB1[aidx] += ab1;
589
590	fSectorSumAB0[aidx] += ab0;
591	fSectorSumAB1[aidx] += ab1;
592	}
593
594	fPedestals->SetReadyToSave();
595	fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain;
596
597	return kTRUE;
598	}
599
600	// --------------------------------------------------------------------------
601	//
602	// Compute signal mean and rms in the whole run and store it in MPedestalCam
603	//
604	Int_t MPedCalcPedRun::PostProcess()
605	{
606
607	// Compute pedestals and rms from the whole run
608	const ULong_t n = fNumSamplesTot;
609	const ULong_t nevts = GetNumExecutions();
610
611	MRawEvtPixelIter pixel(fRawEvt);
612
613	while (pixel.Next())
614	{
615
616	const Int_t pixid = pixel.GetPixelId();
617	const UInt_t aidx = (*fGeom)[pixid].GetAidx();
618	const UInt_t sector = (*fGeom)[pixid].GetSector();
619
620	fAreaValid [aidx]++;
621	fSectorValid[sector]++;
622
623	const Float_t sum = fSumx.At(pixid);
624	const Float_t sum2 = fSumx2.At(pixid);
625	const Float_t higainped = sum/n;
626	//
627	// The old version:
628	//
629	// const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.));
630	//
631	// The new version:
632	//
633	// 1. Calculate the Variance of the sums:
634	Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.);
635	// 2. Scale the variance to the number of slices:
636	higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain);
637	// 3. Calculate the RMS from the Variance:
638	const Float_t rms = higainVar<0 ? 0. : TMath::Sqrt(higainVar);
639	// 4. Calculate the amplitude of the 150MHz "AB" noise
640	const Float_t abOffs = (fSumAB0.At(pixid) - fSumAB1.At(pixid)) / n;
641
642	(*fPedestals)[pixid].Set(higainped,rms,abOffs);
643	}
644
645	//
646	// Loop over the (two) area indices to get the averaged pedestal per aidx
647	//
648	for (Int_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
649	{
650
651	const Int_t napix = fAreaValid.At(aidx);
652
653	if (napix == 0)
654	continue;
655
656	const Float_t sum = fAreaSumx.At(aidx);
657	const Float_t sum2 = fAreaSumx2.At(aidx);
658	const ULong_t an = napix * n;
659	const ULong_t aevts = napix * nevts;
660
661	const Float_t higainped = sum/an;
662
663	// 1. Calculate the Variance of the sums:
664	Float_t higainVar = (sum2-sum*sum/aevts)/(aevts-1.);
665	// 2. Scale the variance to the number of slices:
666	higainVar /= fWindowSizeHiGain+fWindowSizeLoGain;
667	// 3. Calculate the RMS from the Variance:
668	Float_t higainrms = TMath::Sqrt(higainVar);
669	// 4. Re-scale it with the square root of the number of involved pixels
670	// in order to be comparable to the mean of pedRMS of that area
671	higainrms *= TMath::Sqrt((Float_t)napix);
672	// 5. Calculate the amplitude of the 150MHz "AB" noise
673	const Float_t abOffs = (fAreaSumAB0.At(aidx) - fAreaSumAB1.At(aidx)) / an;
674
675	fPedestals->GetAverageArea(aidx).Set(higainped, higainrms,abOffs);
676	}
677
678	//
679	// Loop over the (six) sector indices to get the averaged pedestal per sector
680	//
681	for (Int_t sector=0; sector<fSectorValid.GetSize(); sector++)
682	{
683
684	const Int_t nspix = fSectorValid.At(sector);
685
686	if (nspix == 0)
687	continue;
688
689	const Float_t sum = fSectorSumx.At(sector);
690	const Float_t sum2 = fSectorSumx2.At(sector);
691	const ULong_t sn = nspix * n;
692	const ULong_t sevts = nspix * nevts;
693
694	const Float_t higainped = sum/sn;
695
696	// 1. Calculate the Variance of the sums:
697	Float_t higainVar = (sum2-sum*sum/sevts)/(sevts-1.);
698	// 2. Scale the variance to the number of slices:
699	higainVar /= fWindowSizeHiGain+fWindowSizeLoGain;
700	// 3. Calculate the RMS from the Variance:
701	Float_t higainrms = TMath::Sqrt(higainVar);
702	// 4. Re-scale it with the square root of the number of involved pixels
703	// in order to be comparable to the mean of pedRMS of that sector
704	higainrms *= TMath::Sqrt((Float_t)nspix);
705	// 5. Calculate the amplitude of the 150MHz "AB" noise
706	const Float_t abOffs = (fSectorSumAB0.At(sector) - fSectorSumAB1.At(sector)) / sn;
707
708	fPedestals->GetAverageSector(sector).Set(higainped, higainrms, abOffs);
709	}
710
711	fPedestals->SetTotalEntries(fNumSamplesTot);
712	fPedestals->SetReadyToSave();
713
714	return kTRUE;
715	}
716
717	Int_t MPedCalcPedRun::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
718	{
719	if (MExtractor::ReadEnv(env, prefix, print)==kERROR)
720	return kERROR;
721
722	Byte_t hw = fWindowSizeHiGain;
723	Byte_t lw = fWindowSizeLoGain;
724
725	Bool_t rc = kFALSE;
726
727	if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print))
728	{
729	hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw);
730	rc=kTRUE;
731	}
732
733	if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print))
734	{
735	lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw);
736	rc=kTRUE;
737	}
738
739	if (rc)
740	SetWindowSize(hw, lw);
741
742	return rc;
743	}

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

Download in other formats: