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

Last change on this file since 4624 was 4613, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 23.0 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
137ClassImp(MPedCalcPedRun);
138
139using namespace std;
140
141const Byte_t MPedCalcPedRun::fgHiGainFirst = 0;
142const Byte_t MPedCalcPedRun::fgHiGainLast = 29;
143const Byte_t MPedCalcPedRun::fgLoGainFirst = 0;
144const Byte_t MPedCalcPedRun::fgLoGainLast = 14;
145const Byte_t MPedCalcPedRun::fgHiGainWindowSize = 14;
146const 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//
162MPedCalcPedRun::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//
185void MPedCalcPedRun::Clear(const Option_t *o)
186{
187 fNumSamplesTot = 0;
188
189 fRawEvt = NULL;
190 fRunHeader = NULL;
191 fPedestals = NULL;
192}
193
194// --------------------------------------------------------------------------
195//
196// SetRange:
197//
198// Calls:
199// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
200// - SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
201//
202void MPedCalcPedRun::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
203{
204 MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
205
206 //
207 // Redo the checks if the window is still inside the ranges
208 //
209 SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
210}
211
212// --------------------------------------------------------------------------
213//
214// Checks:
215// - if a window is odd, subtract one
216// - if a window is bigger than the one defined by the ranges, set it to the available range
217// - if a window is smaller than 2, set it to 2
218//
219// Sets:
220// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain
221// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain
222// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
223// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)
224//
225void MPedCalcPedRun::SetWindowSize(Byte_t windowh, Byte_t windowl)
226{
227
228 fWindowSizeHiGain = windowh & ~1;
229 fWindowSizeLoGain = windowl & ~1;
230
231 if (fWindowSizeHiGain != windowh)
232 *fLog << warn << GetDescriptor() << ": Hi Gain window size has to be even, set to: "
233 << int(fWindowSizeHiGain) << " samples " << endl;
234
235 if (fWindowSizeLoGain != windowl)
236 *fLog << warn << GetDescriptor() << ": Lo Gain window size has to be even, set to: "
237 << int(fWindowSizeLoGain) << " samples " << endl;
238
239 const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1;
240 const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1;
241
242 if (fWindowSizeHiGain > availhirange)
243 {
244 *fLog << warn << GetDescriptor()
245 << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",(int)fWindowSizeHiGain,
246 " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl;
247 *fLog << warn << GetDescriptor()
248 << ": Will set window size to: " << (int)availhirange << endl;
249 fWindowSizeHiGain = availhirange;
250 }
251
252 if (fWindowSizeLoGain > availlorange)
253 {
254 *fLog << warn << GetDescriptor()
255 << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",(int)fWindowSizeLoGain,
256 " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl;
257 *fLog << warn << GetDescriptor()
258 << ": Will set window size to: " << (int)availlorange << endl;
259 fWindowSizeLoGain = availlorange;
260 }
261
262
263 fNumHiGainSamples = (Float_t)fWindowSizeHiGain;
264 fNumLoGainSamples = (Float_t)fWindowSizeLoGain;
265
266 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
267 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
268}
269
270// --------------------------------------------------------------------------
271//
272// Look for the following input containers:
273//
274// - MRawEvtData
275// - MRawRunHeader
276// - MGeomCam
277//
278// The following output containers are also searched and created if
279// they were not found:
280//
281// - MPedestalCam
282//
283Int_t MPedCalcPedRun::PreProcess( MParList *pList )
284{
285 Clear();
286
287 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
288 if (!fRawEvt)
289 {
290 *fLog << err << "MRawEvtData not found... aborting." << endl;
291 return kFALSE;
292 }
293
294 fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
295 if (!fRunHeader)
296 {
297 *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl;
298 return kFALSE;
299 }
300
301 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
302 if (!fGeom)
303 {
304 *fLog << err << "MGeomCam not found... aborting." << endl;
305 return kFALSE;
306 }
307
308 fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam");
309 if (!fPedestals)
310 return kFALSE;
311
312 return kTRUE;
313}
314
315// --------------------------------------------------------------------------
316//
317// The ReInit searches for:
318// - MRawRunHeader::GetNumSamplesHiGain()
319// - MRawRunHeader::GetNumSamplesLoGain()
320//
321// In case that the variables fHiGainLast and fLoGainLast are smaller than
322// the even part of the number of samples obtained from the run header, a
323// warning is given an the range is set back accordingly. A call to:
324// - SetRange(fHiGainFirst, fHiGainLast-diff, fLoGainFirst, fLoGainLast) or
325// - SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff)
326// is performed in that case. The variable diff means here the difference
327// between the requested range (fHiGainLast) and the available one. Note that
328// the functions SetRange() are mostly overloaded and perform more checks,
329// modifying the ranges again, if necessary.
330//
331Bool_t MPedCalcPedRun::ReInit(MParList *pList)
332{
333 Int_t lastdesired = (Int_t)fLoGainLast;
334 Int_t lastavailable = (Int_t)fRunHeader->GetNumSamplesLoGain()-1;
335
336 if (lastdesired > lastavailable)
337 {
338 const Int_t diff = lastdesired - lastavailable;
339 *fLog << endl;
340 *fLog << warn << GetDescriptor()
341 << Form("%s%2i%s%2i%s%2i%s",": Selected Lo Gain FADC Window [",
342 (int)fLoGainFirst,",",lastdesired,
343 "] ranges out of the available limits: [0,",lastavailable,"].") << endl;
344 *fLog << GetDescriptor() << ": Will reduce the upper edge to " << (int)(fLoGainLast - diff) << endl;
345 SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff);
346 }
347
348 lastdesired = (Int_t)fHiGainLast;
349 lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
350
351 if (lastdesired > lastavailable)
352 {
353 const Int_t diff = lastdesired - lastavailable;
354 *fLog << endl;
355 *fLog << warn << GetDescriptor()
356 << Form("%s%2i%s%2i%s%2i%s",": Selected Hi Gain Range [",
357 (int)fHiGainFirst,",",lastdesired,
358 "] ranges out of the available limits: [0,",lastavailable,"].") << endl;
359 *fLog << warn << GetDescriptor()
360 << Form("%s%2i%s",": Will possibly use ",diff," samples from the Low-Gain for the High-Gain range")
361 << endl;
362 fHiGainLast -= diff;
363 fHiLoLast = diff;
364 }
365
366 lastdesired = (Int_t)fHiGainFirst+fWindowSizeHiGain-1;
367 lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
368
369 if (lastdesired > lastavailable)
370 {
371 const Int_t diff = lastdesired - lastavailable;
372 *fLog << endl;
373 *fLog << warn << GetDescriptor()
374 << Form("%s%2i%s%2i%s",": Selected Hi Gain FADC Window size ",
375 (int)fWindowSizeHiGain,
376 " ranges out of the available limits: [0,",lastavailable,"].") << endl;
377 *fLog << warn << GetDescriptor()
378 << Form("%s%2i%s",": Will use ",diff," samples from the Low-Gain for the High-Gain extraction")
379 << endl;
380
381 if ((Int_t)fWindowSizeHiGain > diff)
382 {
383 fWindowSizeHiGain -= diff;
384 fWindowSizeLoGain += diff;
385 }
386 else
387 {
388 fWindowSizeLoGain += fWindowSizeHiGain;
389 fLoGainFirst = diff-fWindowSizeHiGain;
390 fWindowSizeHiGain = 0;
391 }
392 }
393
394
395 const Int_t npixels = fPedestals->GetSize();
396
397 if (fSumx.GetSize()==0)
398 {
399 fSumx. Set(npixels);
400 fSumx2.Set(npixels);
401
402 fSumx.Reset();
403 fSumx2.Reset();
404 }
405
406 if (fWindowSizeHiGain == 0 && fWindowSizeLoGain == 0)
407 {
408 *fLog << err << GetDescriptor()
409 << ": Number of extracted Slices is 0, abort ... " << endl;
410 return kFALSE;
411 }
412
413
414 *fLog << endl;
415 *fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeHiGain
416 << " HiGain FADC samples starting with slice: " << (int)fHiGainFirst << endl;
417 *fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeLoGain
418 << " LoGain FADC samples starting with slice: " << (int)fLoGainFirst << endl;
419
420 return kTRUE;
421}
422
423// --------------------------------------------------------------------------
424//
425// Fill the MPedestalCam container with the signal mean and rms for the event.
426// Store the measured signal in arrays fSumx and fSumx2 so that we can
427// calculate the overall mean and rms in the PostProcess()
428//
429Int_t MPedCalcPedRun::Process()
430{
431 MRawEvtPixelIter pixel(fRawEvt);
432
433 while (pixel.Next())
434 {
435 const UInt_t idx = pixel.GetPixelId();
436
437 Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst;
438 Byte_t *end = ptr + fWindowSizeHiGain;
439
440 UInt_t sum = 0;
441 UInt_t sqr = 0;
442
443 if (fWindowSizeHiGain != 0)
444 {
445 do
446 {
447 sum += *ptr;
448 sqr += *ptr * *ptr;
449 }
450 while (++ptr != end);
451 }
452
453 if (fWindowSizeLoGain != 0)
454 {
455 ptr = pixel.GetLoGainSamples() + fLoGainFirst;
456 end = ptr + fWindowSizeLoGain;
457
458 do
459 {
460 sum += *ptr;
461 sqr += *ptr * *ptr;
462 }
463 while (++ptr != end);
464 }
465
466 const Float_t msum = (Float_t)sum;
467
468 //
469 // These three lines have been uncommented by Markus Gaug
470 // If anybody needs them, please contact me!!
471 //
472 // const Float_t higainped = msum/fNumHiGainSlices;
473 // const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSlices)/(fNumHiGainSlices-1.));
474 // (*fPedestals)[idx].Set(higainped, higainrms);
475
476 fSumx[idx] += msum;
477 //
478 // The old version:
479 //
480 // const Float_t msqr = (Float_t)sqr;
481 // fSumx2[idx] += msqr;
482 //
483 // The new version:
484 //
485 const Float_t sqrsum = msum*msum;
486 fSumx2[idx] += sqrsum;
487 }
488
489 fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain;
490
491 return kTRUE;
492/*
493 MRawEvtPixelIter pixel(fRawEvt);
494
495 while (pixel.Next())
496 {
497 const UInt_t idx = pixel.GetPixelId();
498 const UInt_t aidx = (*fGeom)[idx].GetAidx();
499 const UInt_t sector = (*fGeom)[idx].GetSector();
500
501 Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst;
502 Byte_t *end = ptr + fWindowSizeHiGain;
503
504 UInt_t sum = 0;
505 UInt_t sqr = 0;
506
507 if (fWindowSizeHiGain != 0)
508 {
509 do
510 {
511 sum += *ptr;
512 sqr += *ptr * *ptr;
513 }
514 while (++ptr != end);
515 }
516
517 if (fWindowSizeLoGain != 0)
518 {
519
520 ptr = pixel.GetLoGainSamples() + fLoGainFirst;
521 end = ptr + fWindowSizeLoGain;
522
523 do
524 {
525 sum += *ptr;
526 sqr += *ptr * *ptr;
527 }
528 while (++ptr != end);
529
530 }
531
532 const Float_t msum = (Float_t)sum;
533
534 //
535 // These three lines have been uncommented by Markus Gaug
536 // If anybody needs them, please contact me!!
537 //
538 // const Float_t higainped = msum/fNumHiGainSlices;
539 // const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSlices)/(fNumHiGainSlices-1.));
540 // (*fPedestals)[idx].Set(higainped, higainrms);
541
542 fSumx[idx] += msum;
543 fAreaSumx[aidx] += msum;
544 fSectorSumx[sector] += msum;
545 //
546 // The old version:
547 //
548 // const Float_t msqr = (Float_t)sqr;
549 // fSumx2[idx] += msqr;
550 //
551 // The new version:
552 //
553 const Float_t sqrsum = msum*msum;
554 fSumx2[idx] += sqrsum;
555 fAreaSumx2[aidx] += sqrsum;
556 fSectorSumx2[sector] += sqrsum;
557 }
558
559 fPedestals->SetReadyToSave();
560 fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain;
561*/
562 return kTRUE;
563}
564
565// --------------------------------------------------------------------------
566//
567// Compute signal mean and rms in the whole run and store it in MPedestalCam
568//
569Int_t MPedCalcPedRun::PostProcess()
570{
571 const ULong_t nevts = GetNumExecutions();
572 if (nevts<1)
573 {
574 *fLog << err << "ERROR - Not enough events recorded...abort." << endl;
575 return kFALSE;
576 }
577
578 *fLog << flush << inf << "Calculating pedestals..." << flush;
579
580 // Compute pedestals and rms from the whole run
581 const ULong_t n = fNumSamplesTot;
582 const ULong_t npix = fGeom->GetNumPixels();
583
584 for (UInt_t pixidx=0; pixidx<npix; pixidx++)
585 {
586 const Float_t sum = fSumx.At(pixidx);
587 const Float_t sum2 = fSumx2.At(pixidx);
588 const Float_t higainped = sum/n;
589 //
590 // The old version:
591 //
592 // const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.));
593 //
594 // The new version:
595 //
596 // 1. Calculate the Variance of the sums:
597 Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.);
598 // 2. Scale the variance to the number of slices:
599 higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain);
600 // 3. Calculate the RMS from the Variance:
601 (*fPedestals)[pixidx].Set(higainped, higainVar < 0 ? 0. : TMath::Sqrt(higainVar), 0, nevts);
602 }
603
604 *fLog << flush << inf << "Calculating means..." << flush;
605
606 fPedestals->SetTotalEntries(fNumSamplesTot);
607 fPedestals->ReCalc(*fGeom);
608 fPedestals->SetReadyToSave();
609
610 return kTRUE;
611
612 //
613 // Loop over the (two) area indices to get the averaged pedestal per aidx
614 //
615 /*
616
617 // Compute pedestals and rms from the whole run
618 const ULong_t n = fNumSamplesTot;
619 const ULong_t nevts = GetNumExecutions();
620
621 MRawEvtPixelIter pixel(fRawEvt);
622
623 while (pixel.Next())
624 {
625 const Int_t pixid = pixel.GetPixelId();
626 const UInt_t aidx = (*fGeom)[pixid].GetAidx();
627 const UInt_t sector = (*fGeom)[pixid].GetSector();
628
629 fAreaValid [aidx]++;
630 fSectorValid[sector]++;
631
632 const Float_t sum = fSumx.At(pixid);
633 const Float_t sum2 = fSumx2.At(pixid);
634 const Float_t higainped = sum/n;
635 //
636 // The old version:
637 //
638 // const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.));
639 //
640 // The new version:
641 //
642 // 1. Calculate the Variance of the sums:
643 Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.);
644 // 2. Scale the variance to the number of slices:
645 higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain);
646 // 3. Calculate the RMS from the Variance:
647 (*fPedestals)[pixid].Set(higainped, higainVar < 0 ? 0. : TMath::Sqrt(higainVar), 0, nevts);
648 }
649 for (Int_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
650 {
651
652 const Int_t napix = fAreaValid.At(aidx);
653
654 if (napix == 0)
655 continue;
656
657 const Float_t sum = fAreaSumx.At(aidx);
658 const Float_t sum2 = fAreaSumx2.At(aidx);
659 const ULong_t an = napix * n;
660 const ULong_t aevts = napix * nevts;
661
662 const Float_t higainped = sum/an;
663
664 // 1. Calculate the Variance of the sums:
665 Float_t higainVar = (sum2-sum*sum/aevts)/(aevts-1.);
666 // 2. Scale the variance to the number of slices:
667 higainVar /= fWindowSizeHiGain+fWindowSizeLoGain;
668 // 3. Calculate the RMS from the Variance:
669 Float_t higainrms = TMath::Sqrt(higainVar);
670 // 4. Re-scale it with the square root of the number of involved pixels
671 // in order to be comparable to the mean of pedRMS of that area
672 higainrms *= TMath::Sqrt((Float_t)napix);
673
674 fPedestals->GetAverageArea(aidx).Set(higainped, higainrms);
675 }
676 */
677 //
678 // Loop over the (six) sector indices to get the averaged pedestal per sector
679 //
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
706 fPedestals->GetAverageSector(sector).Set(higainped, higainrms);
707 }
708 fPedestals->SetTotalEntries(fNumSamplesTot);
709 fPedestals->SetReadyToSave();
710
711 return kTRUE;
712 */
713}
714
715Int_t MPedCalcPedRun::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
716{
717 if (MExtractor::ReadEnv(env, prefix, print)==kERROR)
718 return kERROR;
719
720 Byte_t hw = fWindowSizeHiGain;
721 Byte_t lw = fWindowSizeLoGain;
722
723 Bool_t rc = kFALSE;
724
725 if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print))
726 {
727 hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw);
728 rc=kTRUE;
729 }
730
731 if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print))
732 {
733 lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw);
734 rc=kTRUE;
735 }
736
737 if (rc)
738 SetWindowSize(hw, lw);
739
740 return rc;
741}
Note: See TracBrowser for help on using the repository browser.