source: trunk/MagicSoft/Mars/mpedestal/MExtractPedestal.cc@ 8203

Last change on this file since 8203 was 8151, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 30.8 KB
Line 
1/* ======================================================================== *\
2! $Name: not supported by cvs2svn $:$Id: MExtractPedestal.cc,v 1.23 2006-10-24 07:58:13 tbretz Exp $
3! --------------------------------------------------------------------------
4!
5! *
6! * This file is part of MARS, the MAGIC Analysis and Reconstruction
7! * Software. It is distributed to you in the hope that it can be a useful
8! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
9! * It is distributed WITHOUT ANY WARRANTY.
10! *
11! * Permission to use, copy, modify and distribute this software and its
12! * documentation for any purpose is hereby granted without fee,
13! * provided that the above copyright notice appear in all copies and
14! * that both that copyright notice and this permission notice appear
15! * in supporting documentation. It is provided "as is" without express
16! * or implied warranty.
17! *
18!
19!
20! Author(s): Markus Gaug 01/2004 <mailto:markus@ifae.es>
21! Author(s): Thomas Bretz 01/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
22!
23! Copyright: MAGIC Software Development, 2000-2006
24!
25!
26\* ======================================================================== */
27
28/////////////////////////////////////////////////////////////////////////////
29//
30// MExtractPedestal
31//
32// Pedestal Extractor base class
33//
34// Input Containers:
35// MRawEvtData
36// MRawRunHeader
37// MRawEvtHeader
38// MGeomCam
39// MPedestalCam
40//
41// Output Containers:
42// MPedestalCam
43//
44// This class should be used for pedestal extractors with the following facilities:
45// a) Standardized calculation of AB-noise, mean pedestals and RMS
46// b) Standardized treatment of area- and sector averaged pedestals values
47// c) Possibility to use a signal extractor to be applied on the pedestals
48// d) Possibility to handle two MPedestalCams: one for the signal extractor and
49// a second to be filled during the pedestal calculating process.
50//
51// ad a): Every calculated value is referred to one FADC slice (e.g. Mean pedestal per slice),
52// RMS per slice.
53// MExtractPedestal applies the following formula (1):
54//
55// Pedestal per slice = sum(x_i) / n / slices
56// PedRMS per slice = Sqrt( ( sum(x_i^2) - sum(x_i)^2/n ) / n-1 / slices )
57// AB-Offset per slice = (sumAB0 - sumAB1) / n / slices
58//
59// where x_i is the sum of "slices" FADC slices and sum means the sum over all
60// events. "n" is the number of events, "slices" is the number of summed FADC samples.
61//
62// Note that the slice-to-slice fluctuations are not Gaussian, but Poissonian, thus
63// asymmetric and they are correlated.
64//
65// It is important to know that the Pedestal per slice and PedRMS per slice depend
66// on the number of used FADC slices, as seen in the following plots:
67//
68//Begin_Html
69/*
70<img src="images/PedestalStudyInner.gif">
71*/
72//End_Html
73//
74//Begin_Html
75/*
76<img src="images/PedestalStudyOuter.gif">
77*/
78//End_Html
79//
80// The plots show the inner and outer pixels, respectivly and have the following meaning:
81//
82// 1) The calculated mean pedestal per slice (from MPedCalcPedRun)
83// 2) The fitted mean pedestal per slice (from MHPedestalCam)
84// 3) The calculated pedestal RMS per slice (from MPedCalcPedRun)
85// 4) The fitted sigma of the pedestal distribution per slice
86// (from MHPedestalCam)
87// 5) The relative difference between calculation and histogram fit
88// for the mean
89// 6) The relative difference between calculation and histogram fit
90// for the sigma or RMS, respectively.
91//
92// The calculated means do not change significantly except for the case of 2 slices,
93// however the RMS changes from 5.7 per slice in the case of 2 extracted slices
94// to 8.3 per slice in the case of 26 extracted slices. This change is very significant.
95//
96// ad b) Every calculated value is referred to one FADC slice and one (averaged) pixel,
97// (e.g. Mean Pedestal per area index per slice per pixel, etc. )
98//
99// MExtractPedestal applies the following formula (2):
100//
101// Averaged Pedestal per slice = sum(x_i) / n / slices / n_pix
102// PedRMS per slice = Sqrt( ( sum(x_i^2) - sum(x_i)^2/n ) / n-1 / slices / n_pix )
103// AB-Offset per slice = (sumAB0 - sumAB1) / n / slices / n_pix
104//
105// where x_i is the sum of "slices" FADC slices and sum means the sum over all
106// events and all concerned pixels.
107// "n" is the number of events, "slices" is the number of summed FADC samples and
108// "n_pix" is the number of pixels belonging to the specific area index or camera sector.
109//
110// Calculating these averaged event-by-event values is very important to trace coherent
111// fluctuations. An example is given in the following plots:
112//
113//Begin_Html
114/*
115<img src="images/PedestalOscillations.gif">
116*/
117//End_Html
118//
119// The plots show the extracted pedestals of the inner pixels (obtained
120// with MHPedestalCam), averaged on an event-by-event basis from
121// run 13428 with switched off camera LV.
122// The meaning of the four plots is:
123//
124// 1) The distribution of the averaged pedestals
125// 2) The averaged pedestals vs. time.
126// One can see clearly the oscillation pattern
127// 3) The fourier transform of the averaged pedestals vs. time.
128// One can see clearly a peak at a certain frequency
129// 4) The projection of the fourier components with the non-exponential
130// (and therefore significant) outlier.
131//
132// ad c) Many signal extractors, especially those using a sliding window
133// have biases and their resolutions for zero-signals do not agree
134// with the pedestal RMS. For the F-Factor method in the calibration
135// and the image cleaning, however, both have to be known and measured.
136//
137// For this reason, a signal extractor can be handed over to the
138// pedestal extractor and applied on the pedestal events with the
139// function SetExtractor().
140// The results will get stored in an MPedestalCam.
141//
142// Note that only extractors deriving from MExtractTimeAndCharge
143// can be used.
144//
145// ad d) The signal extractors themselves need a pedestal to be subtracted
146// from the FADC slices.
147// If the user wishes that the pededestals do not get overwritten by
148// the results from the signal extractor, a different named MPedestalCam
149// can be created with the function: SetNamePedestalOut().
150//
151// See also: MPedestalCam, MPedestalPix, MPedCalcPedRun, MPedCalcFromLoGain
152//
153/////////////////////////////////////////////////////////////////////////////
154#include "MExtractPedestal.h"
155
156#include "MParList.h"
157
158#include "MLog.h"
159#include "MLogManip.h"
160
161#include "MRawRunHeader.h"
162#include "MRawEvtHeader.h"
163#include "MRawEvtPixelIter.h"
164
165#include "MPedestalPix.h"
166#include "MPedestalCam.h"
167
168#include "MGeomPix.h"
169#include "MGeomCam.h"
170
171#include "MExtractTimeAndCharge.h"
172#include "MPedestalSubtractedEvt.h"
173
174ClassImp(MExtractPedestal);
175
176using namespace std;
177
178const TString MExtractPedestal::fgNamePedestalCam = "MPedestalCam";
179const TString MExtractPedestal::fgNameRawEvtData = "MRawEvtData";
180
181const UShort_t MExtractPedestal::fgCheckWinFirst = 0;
182const UShort_t MExtractPedestal::fgCheckWinLast = 29;
183const UShort_t MExtractPedestal::fgMaxSignalVar = 40;
184
185// --------------------------------------------------------------------------
186//
187// Default constructor:
188//
189// Sets:
190// - all pointers to NULL
191//
192// Calls:
193// - Clear()
194//
195MExtractPedestal::MExtractPedestal(const char *name, const char *title)
196 : fGeom(NULL), fPedestalsIn(NULL), fPedestalsInter(NULL),
197 fPedestalsOut(NULL), fExtractor(NULL), fSignal(0),
198 fExtractWinFirst(0), fExtractWinSize(0), fUseSpecialPixels(kFALSE)
199{
200 fName = name ? name : "MExtractPedestal";
201 fTitle = title ? title : "Base class to calculate pedestals";
202
203 SetIntermediateStorage( kFALSE );
204 SetRandomCalculation ( kTRUE );
205
206 SetNamePedestalCamIn();
207 SetNamePedestalCamOut();
208 SetNamePedestalCamInter();
209 SetNameRawEvtData();
210
211 SetCheckRange(fgCheckWinFirst, fgCheckWinLast);
212 SetMaxSignalVar(fgMaxSignalVar);
213
214 Clear();
215}
216
217// --------------------------------------------------------------------------
218//
219// Call reset() of all Arays
220//
221void MExtractPedestal::ResetArrays()
222{
223 // Reset contents of arrays.
224 fSumx.Reset();
225 fSumx2.Reset();
226 fSumAB0.Reset();
227 fSumAB1.Reset();
228 fAreaSumx.Reset();
229 fAreaSumx2.Reset();
230 fAreaSumAB0.Reset();
231 fAreaSumAB1.Reset();
232 fAreaFilled.Reset();
233 fAreaValid.Reset();
234 fSectorSumx.Reset();
235 fSectorSumx2.Reset();
236 fSectorSumAB0.Reset();
237 fSectorSumAB1.Reset();
238 fSectorFilled.Reset();
239 fSectorValid.Reset();
240
241}
242
243// --------------------------------------------------------------------------
244//
245// Resets Arrays:
246//
247// Sets:
248// - fRawEvt to NULL
249// - fRunHeader to NULL
250//
251void MExtractPedestal::Clear(const Option_t *o)
252{
253
254 fRawEvt = NULL;
255 fRunHeader = NULL;
256
257 // If the size is yet set, set the size
258 if (fSumx.GetSize()>0)
259 ResetArrays();
260
261}
262
263// --------------------------------------------------------------------------
264//
265// Checks:
266// - if a window is odd
267//
268Bool_t MExtractPedestal::SetExtractWindow(UShort_t windowf, UShort_t windows)
269{
270
271 Bool_t rc = kTRUE;
272
273 const Int_t odd = TMath::Odd(windows);
274
275 if (odd && !fExtractor)
276 {
277 *fLog << warn << GetDescriptor();
278 *fLog << " - WARNING: Window size in SetExtractWindow has to be even... ";
279 *fLog << " lowering from " << windows << " to ";
280 windows -= 1;
281 *fLog << windows << "!" << endl;
282 rc = kFALSE;
283 }
284
285 if (windows==0)
286 {
287 *fLog << warn << GetDescriptor();
288 *fLog << " - WARNING: Window size in SetExtractWindow has to be > 0... adjusting to 2!" << endl;
289 windows = 2;
290 rc = kFALSE;
291 }
292
293 fExtractWinSize = windows;
294 fExtractWinFirst = windowf;
295 fExtractWinLast = fExtractWinFirst+fExtractWinSize-1;
296
297 return rc;
298}
299
300// --------------------------------------------------------------------------
301//
302// SetCheckRange:
303//
304// Exits, if the first argument is smaller than 0
305// Exits, if the the last argument is smaller than the first
306//
307Bool_t MExtractPedestal::SetCheckRange(UShort_t chfirst, UShort_t chlast)
308{
309
310 Bool_t rc = kTRUE;
311
312 if (chlast<=chfirst)
313 {
314 *fLog << warn << GetDescriptor();
315 *fLog << " - WARNING: Last slice in SetCheckRange smaller than first slice... set to first+2" << endl;
316 chlast = chfirst+1;
317 rc = kFALSE;
318 }
319
320 fCheckWinFirst = chfirst;
321 fCheckWinLast = chlast;
322
323 return rc;
324}
325
326
327// --------------------------------------------------------------------------
328//
329// Check (and if neccesary: correct) the extraction and check ranges.
330//
331void MExtractPedestal::CheckExtractionWindow()
332{
333 // fSignal->GetNumSamples() not yet initialized!!!
334 const UInt_t num = fRunHeader->GetNumSamplesHiGain()+fRunHeader->GetNumSamplesLoGain();
335
336 if (fCheckWinLast >= num)
337 {
338 *fLog << warn << endl;
339 *fLog << "Selected CheckWindow [" << fCheckWinFirst << "," << fCheckWinLast;
340 *fLog << "] ranges out of range [0," << num-1 << "]." << endl;
341 *fLog << "Resetting upper edge." << endl;
342
343 fCheckWinLast = num;
344 }
345
346 if (fExtractWinLast >= num)
347 {
348 *fLog << warn << endl;
349 *fLog << "Selected ExtractWindow [" << fExtractWinFirst << "," << fExtractWinLast;
350 *fLog << "] ranges out of range [0," << num-1 << "]." << endl;
351 *fLog << "Resetting upper edge." << endl;
352
353 fExtractWinLast = num;
354 }
355
356 fExtractWinSize = fExtractWinLast-fExtractWinFirst+1;
357
358 if (fExtractor || TMath::Even(fExtractWinSize))
359 return;
360
361 fExtractWinLast += fExtractWinLast==num-1 ? -1 : +1;
362
363 *fLog << warn << endl;
364 *fLog << "ExtractionWindow odd...";
365 *fLog << "New extraction window [" << fExtractWinFirst << "," << fExtractWinLast << "]" << endl;
366
367 fExtractWinSize = fExtractWinLast-fExtractWinFirst+1;
368}
369
370// --------------------------------------------------------------------------
371//
372// Look for the following input containers:
373//
374// - MRawEvtData
375// - MRawRunHeader
376// - MRawEvtHeader
377// - MGeomCam
378//
379// The following output containers are also searched and created if
380// they were not found:
381//
382// - MPedestalCam with the name fPedContainerName
383//
384Int_t MExtractPedestal::PreProcess(MParList *pList)
385{
386
387 Clear();
388
389 fRawEvt = (MRawEvtData*)pList->FindObject(AddSerialNumber(fNameRawEvtData));
390 if (!fRawEvt)
391 if (!fUseSpecialPixels)
392 {
393 *fLog << err << AddSerialNumber(fNameRawEvtData) << " not found... aborting." << endl;
394 return kFALSE;
395 }
396
397
398 fRawEvt = (MRawEvtData*)pList->FindObject(AddSerialNumber("MRawEvtData"));
399 if (!fRawEvt)
400 {
401 *fLog << err << AddSerialNumber("MRawEvtData") << " not found... aborting." << endl;
402 return kFALSE;
403 }
404
405 fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
406 if (!fRunHeader)
407 {
408 *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl;
409 return kFALSE;
410 }
411
412 fGeom = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
413 if (!fGeom)
414 {
415 *fLog << err << AddSerialNumber("MGeomCam") << " not found... aborting." << endl;
416 return kFALSE;
417 }
418
419 fSignal = (MPedestalSubtractedEvt*)pList->FindObject(AddSerialNumber("MPedestalSubtractedEvt"));
420 if (!fSignal)
421 {
422 *fLog << err << AddSerialNumber("MPedestalSubtractedEvt") << " not found... aborting." << endl;
423 return kFALSE;
424 }
425
426 if (fExtractor && !fPedestalsIn)
427 {
428 fPedestalsIn = (MPedestalCam*)pList->FindObject(AddSerialNumber(fNamePedestalCamIn), "MPedestalCam");
429 if (!fPedestalsIn)
430 {
431 *fLog << err << AddSerialNumber(fNamePedestalCamIn) << " not found... aborting." << endl;
432 return kFALSE;
433 }
434 }
435
436 if (!fPedestalsInter && fIntermediateStorage)
437 {
438 fPedestalsInter = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", AddSerialNumber(fNamePedestalCamInter));
439 if (!fPedestalsInter)
440 return kFALSE;
441 }
442
443 if (!fPedestalsOut)
444 {
445 fPedestalsOut = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", AddSerialNumber(fNamePedestalCamOut));
446 if (!fPedestalsOut)
447 return kFALSE;
448 }
449
450 return fExtractor ? fExtractor->CallPreProcess(pList) : kTRUE;
451}
452
453//-----------------------------------------------------------------------
454//
455// Call Calc(). If fExtractor!=NULL enclose the call in setting the
456// NoiseCalculation to fRandomCalculation
457//
458Int_t MExtractPedestal::Process()
459{
460 //
461 // Necessary check for extraction of special pixels
462 // together with data which does not yet have them
463 //
464 if (fSumx.GetSize()==0)
465 return kTRUE;
466
467 if (fExtractor)
468 fExtractor->SetNoiseCalculation(fRandomCalculation);
469
470 const Int_t rc = Calc();
471
472 if (fExtractor)
473 fExtractor->SetNoiseCalculation(kFALSE);
474
475 return rc;
476}
477
478// ---------------------------------------------------------------------------------
479//
480// Sets the size (from MPedestalCam::GetSize() ) and resets the following arrays:
481// - fSumx
482// - fSumx2
483// - fSumAB0
484// - fSumAB1
485// - fAreaSumx
486// - fAreaSumx2
487// - fAreaSumAB0
488// - fAreaSumAB1
489// - fAreaFilled
490// - fAreaValid
491// - fSectorSumx
492// - fSectorSumx2
493// - fSectorSumAB0
494// - fSectorSumAB1
495// - fSectorFilled
496// - fSectorValid
497//
498Bool_t MExtractPedestal::ReInit(MParList *pList)
499{
500 // Necessary check for special pixels which might not yet have existed
501 if (!fRawEvt)
502 {
503 if (fRunHeader->GetFormatVersion() > 3)
504 return kTRUE;
505
506 *fLog << err << "ERROR - " << fNameRawEvtData << " [MRawEvtData] has not ";
507 *fLog << "been found and format version > 3... abort." << endl;
508 return kFALSE;
509 }
510
511 // If the size is not yet set, set the size
512 if (fSumx.GetSize()==0)
513 {
514 // Initialize the normal pixels (size of MPedestalCam already set by MGeomApply)
515 const Int_t npixels = fPedestalsOut->GetSize();
516
517 fSumx. Set(npixels);
518 fSumx2. Set(npixels);
519 fSumAB0.Set(npixels);
520 fSumAB1.Set(npixels);
521
522 if (fUseSpecialPixels)
523 {
524 // Initialize size of MPedestalCam in case of special pixels (not done by MGeomApply)
525 const UShort_t nspecial = fRunHeader->GetNumSpecialPixels();
526 if (nspecial == 0)
527 {
528 *fLog << warn << "WARNING - Number of special pixels is 0." << endl;
529 return kTRUE;
530 }
531
532 fPedestalsOut->InitSize((UInt_t)nspecial);
533 }
534 else
535 {
536 // Initialize the averaged areas and sectors (do not exist for special pixels)
537 const Int_t areas = fPedestalsOut->GetNumAverageArea();
538 const Int_t sectors = fPedestalsOut->GetNumAverageSector();
539
540 fAreaSumx. Set(areas);
541 fAreaSumx2. Set(areas);
542 fAreaSumAB0.Set(areas);
543 fAreaSumAB1.Set(areas);
544 fAreaFilled.Set(areas);
545 fAreaValid .Set(areas);
546
547 fSectorSumx. Set(sectors);
548 fSectorSumx2. Set(sectors);
549 fSectorSumAB0.Set(sectors);
550 fSectorSumAB1.Set(sectors);
551 fSectorFilled.Set(sectors);
552 fSectorValid .Set(sectors);
553
554 for (Int_t i=0; i<npixels; i++)
555 {
556 const UInt_t aidx = (*fGeom)[i].GetAidx();
557 const UInt_t sector = (*fGeom)[i].GetSector();
558
559 fAreaValid [aidx] ++;
560 fSectorValid[sector]++;
561 }
562 }
563 }
564
565 if (fExtractor)
566 {
567 if (!fExtractor->ReInit(pList))
568 return kFALSE;
569
570 // If an extractor is set determin the window size automatically!
571 fExtractWinFirst = fExtractor->GetHiGainFirst();
572 fExtractWinLast = fExtractor->GetHiGainLast();
573
574 // fSignal->GetNumSamples() not yet initialized!!!
575 const UInt_t num = fRunHeader->GetNumSamplesHiGain()+fRunHeader->GetNumSamplesLoGain();
576 if (fExtractWinLast >= num)
577 {
578 *fLog << err << GetDescriptor();
579 *fLog << " - ERROR: Selected fExtractWinLast out of range." << endl;
580 return kFALSE;
581 }
582 }
583
584 CheckExtractionWindow();
585
586 return kTRUE;
587}
588
589// ---------------------------------------------------------------------------------
590//
591// PostProcess the extractor if available
592//
593Int_t MExtractPedestal::PostProcess()
594{
595 fPedestalsIn = NULL;
596 return fExtractor ? fExtractor->CallPostProcess() : kTRUE;
597}
598
599// ---------------------------------------------------------------------------------
600//
601// Check whether the signal variation between fCheckWinFirst and fCheckWinLast
602// exceeds fMaxSignalVar or the signal is greater than 250
603//
604Bool_t MExtractPedestal::CheckVariation(UInt_t idx) const
605{
606 // This is the fast workaround to put hi- and lo-gains together
607 Byte_t *slices = fSignal->GetSamplesRaw(idx);//pixel.GetSamples();
608
609 // Start 'real' work
610 UShort_t max = 0;
611 UShort_t min = (UShort_t)-1;
612
613 // Find the maximum and minimum signal per slice in the high gain window
614 for (Byte_t *slice=slices+fCheckWinFirst; slice<=slices+fCheckWinLast; slice++)
615 {
616 if (*slice > max)
617 max = *slice;
618 if (*slice < min)
619 min = *slice;
620 }
621
622 // If the maximum in the high gain window is smaller than
623 return max-min<fMaxSignalVar && max<250;
624}
625
626// ---------------------------------------------------------------------------------
627//
628// Invoke the hi-gain extraction starting at fExtractWinFirst+offset
629// for fExtractWinLast-fExtractWinFirst+1 slices. If Noise calculation
630// is set it is up to the signal extractor to do the right thing.
631//
632// Returns the extracted signal.
633//
634Float_t MExtractPedestal::CalcExtractor(const MRawEvtPixelIter &pixel, Int_t offset) const
635{
636 // Use the same extraction window as for signal extraction
637 const Int_t first = fExtractWinFirst;
638 const Int_t last = fExtractWinLast;
639
640 const Int_t start = first+offset;
641
642 const Int_t range = last-first+1;
643
644 // This check is already done in CheckExtractionWindow
645 // if (range>pixel.GetNumSamples()-start)
646 // range = pixel.GetNumSamples()-start;
647
648 const Int_t idx = pixel.GetPixelId();
649
650 // Do some handling if maxpos is last slice?
651 const Int_t maxposhi = fSignal->GetMax(idx, start, start+range-1);
652
653 const Float_t *sig = fSignal->GetSamples(idx);
654
655 // The pedestal is extracted with the hi-gain extractor (eg. digital
656 // filter weights) but from the lo-gains
657 Float_t dummy[3];
658 Float_t sum = 0;
659 fExtractor->FindTimeAndChargeHiGain2(sig+start, range, sum,
660 dummy[0], dummy[1], dummy[2],
661 0, maxposhi);
662 return sum;
663}
664
665// ---------------------------------------------------------------------------------
666//
667// Sum slices from fExtractWinFirst to fExtractWinLast. The total sum is
668// returned. ab0 and ab1 will contain the total sum splitted by the
669// AB-flag. If the AB-flag is invalid ab0=ab1=0 is returned.
670//
671UInt_t MExtractPedestal::CalcSums(const MRawEvtPixelIter &pixel, Int_t offset, UInt_t &ab0, UInt_t &ab1) const
672{
673 const Int_t first = fExtractWinFirst+offset;
674
675 Byte_t *ptr = fSignal->GetSamplesRaw(pixel.GetPixelId())+first;
676 Byte_t *end = ptr + fExtractWinSize;
677
678 Int_t abflag = pixel.HasABFlag() + first;
679
680 UInt_t ab[2] = { 0, 0 };
681 while (ptr<end)
682 ab[abflag++ & 0x1] += *ptr++;
683
684 // This check if for old data without AB-Flag in the data
685 const Bool_t valid = pixel.IsABFlagValid();
686
687 ab0 = valid ? ab[0] : 0;
688 ab1 = valid ? ab[1] : 0;
689
690 return ab0+ab1;
691}
692
693// --------------------------------------------------------------------------
694//
695// The following resources are available:
696// ExtractWindowFirst: 15
697// ExtractWindowSize: 6
698// PedestalUpdate: yes
699// RandomCalculation: yes
700//
701Int_t MExtractPedestal::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
702{
703 Bool_t rc=kFALSE;
704
705 // find resource for fUseSpecialPixels
706 if (IsEnvDefined(env, prefix, "UseSpecialPixels", print))
707 {
708 SetUseSpecialPixels(GetEnvValue(env, prefix, "UseSpecialPixels", fUseSpecialPixels));
709 rc = kTRUE;
710 }
711
712 if (IsEnvDefined(env, prefix, "IntermediateStorage", print))
713 {
714 SetIntermediateStorage(GetEnvValue(env, prefix, "IntermediateStorage", fIntermediateStorage));
715 rc = kTRUE;
716 }
717
718 // find resource for random calculation
719 if (IsEnvDefined(env, prefix, "RandomCalculation", print))
720 {
721 SetRandomCalculation(GetEnvValue(env, prefix, "RandomCalculation", fRandomCalculation));
722 rc = kTRUE;
723 }
724
725 // Find resources for ExtractWindow
726 Int_t ef = fExtractWinFirst;
727 Int_t es = fExtractWinSize;
728 if (IsEnvDefined(env, prefix, "ExtractWinFirst", print))
729 {
730 ef = GetEnvValue(env, prefix, "ExtractWinFirst", ef);
731 rc = kTRUE;
732 }
733 if (IsEnvDefined(env, prefix, "ExtractWinSize", print))
734 {
735 es = GetEnvValue(env, prefix, "ExtractWinSize", es);
736 rc = kTRUE;
737 }
738
739 SetExtractWindow(ef,es);
740
741 // Find resources for CheckWindow
742 Int_t cfs = fCheckWinFirst;
743 Int_t cls = fCheckWinLast;
744 if (IsEnvDefined(env, prefix, "CheckWinFirst", print))
745 {
746 cfs = GetEnvValue(env, prefix, "CheckWinFirst", cfs);
747 rc = kTRUE;
748 }
749 if (IsEnvDefined(env, prefix, "CheckWinLast", print))
750 {
751 cls = GetEnvValue(env, prefix, "CheckWinLast", cls);
752 rc = kTRUE;
753 }
754
755 SetCheckRange(cfs,cls);
756
757 // find resource for maximum signal variation
758 if (IsEnvDefined(env, prefix, "MaxSignalVar", print))
759 {
760 SetMaxSignalVar(GetEnvValue(env, prefix, "MaxSignalVar", fMaxSignalVar));
761 rc = kTRUE;
762 }
763
764 // find resource for MPedestalCam
765 if (IsEnvDefined(env, prefix, "NamePedestalCamIn", print))
766 {
767 SetNamePedestalCamIn(GetEnvValue(env, prefix, "NamePedestalCamIn", fNamePedestalCamIn));
768 rc = kTRUE;
769 }
770
771 if (IsEnvDefined(env, prefix, "NamePedestalCamInter", print))
772 {
773 SetNamePedestalCamInter(GetEnvValue(env, prefix, "NamePedestalCamInter", fNamePedestalCamInter));
774 rc = kTRUE;
775 }
776
777 if (IsEnvDefined(env, prefix, "NamePedestalCamOut", print))
778 {
779 SetNamePedestalCamOut(GetEnvValue(env, prefix, "NamePedestalCamOut", fNamePedestalCamOut));
780 rc = kTRUE;
781 }
782
783 return rc;
784}
785
786// ---------------------------------------------------------------------------------
787//
788// Calculates for pixel "idx":
789//
790// Ped per slice = sum / n / fExtractWinSize;
791// RMS per slice = sqrt { (sum2 - sum*sum/n) / (n-1) / fExtractWinSize }
792// ABOffset per slice = (fSumAB0[idx] - fSumAB1[idx]) / n / fExtractWinSize;
793//
794// Stores the results in MPedestalCam[pixid]
795//
796void MExtractPedestal::CalcPixResults(const UInt_t nevts, const UInt_t pixid)
797{
798 const Double_t sum = fSumx[pixid];
799 const Double_t sum2 = fSumx2[pixid];
800
801 // 1. Calculate the mean of the sums:
802 Double_t ped = sum/nevts;
803
804 // 2. Calculate the Variance of the sums:
805 Double_t var = (sum2-sum*sum/nevts)/(nevts-1.);
806
807 // 3. Calculate the amplitude of the 150MHz "AB" noise
808 Double_t abOffs = (fSumAB0[pixid] - fSumAB1[pixid]) / nevts;
809
810 // 4. Scale the mean, variance and AB-noise to the number of slices:
811 ped /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
812 var /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
813 abOffs /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
814 // The pedestal extracted with the extractor is divided by
815 // the number of hi-gain samples because the calibration
816 // multiplies by this number
817
818 // 5. Calculate the RMS from the Variance:
819 const Double_t rms = var<0 ? 0 : TMath::Sqrt(var);
820
821 // abOffs contains only half of the signal as ped.
822 // Therefor abOffs is not the full, but the half amplitude
823 (*fPedestalsOut)[pixid].Set(ped, rms, abOffs, nevts);
824}
825
826// ---------------------------------------------------------------------------------
827//
828// Calculates for area idx "aidx" with "napix" valid pixels:
829//
830// Ped per slice = sum / nevts / fExtractWinSize / napix;
831// RMS per slice = sqrt { (sum2 - sum*sum/nevts) / (nevts-1) / fExtractWinSize / napix }
832// ABOffset per slice = (fSumAB0[idx] - fSumAB1[idx]) / nevts / fExtractWinSize / napix;
833//
834// Stores the results in MPedestalCam::GetAverageArea(aidx)
835//
836void MExtractPedestal::CalcAreaResults(const UInt_t nevts, const UInt_t napix, const UInt_t aidx)
837{
838 const Double_t sum = fAreaSumx[aidx];
839 const Double_t sum2 = fAreaSumx2[aidx];
840
841 // 1. Calculate the mean of the sums:
842 Double_t ped = sum/nevts;
843
844 // 2. Calculate the Variance of the sums:
845 Double_t var = (sum2/napix-sum*sum/nevts)/(nevts-1.);
846
847 // 3. Calculate the amplitude of the 150MHz "AB" noise
848 Double_t abOffs = (fAreaSumAB0[aidx] - fAreaSumAB1[aidx]) / nevts;
849
850 // 4. Scale the mean, variance and AB-noise to the number of slices:
851 ped /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
852 var /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
853 abOffs /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
854 // The pedestal extracted with the extractor is divided by
855 // the number of hi-gain samples because the calibration
856 // multiplies by this number
857
858 // 5. Scale the mean, variance and AB-noise to the number of pixels:
859 ped /= napix;
860 var /= napix;
861 abOffs /= napix;
862
863 // 6. Calculate the RMS from the Variance:
864 const Double_t rms = var<0 ? 0 : TMath::Sqrt(var);
865
866 // abOffs contains only half of the signal as ped.
867 // Therefor abOffs is not the full, but the half amplitude
868 fPedestalsOut->GetAverageArea(aidx).Set(ped, rms, abOffs, nevts);
869}
870
871// ---------------------------------------------------------------------------------
872//
873// Calculates for sector idx "sector" with "nspix" valid pixels:
874//
875// Ped per slice = sum / nevts / fExtractWinSize / nspix;
876// RMS per slice = sqrt { (sum2 - sum*sum/nevts) / (nevts-1) / fExtractWinSize / nspix }
877// ABOffset per slice = (fSumAB0[idx] - fSumAB1[idx]) / nevts / fExtractWinSize / nspix;
878//
879// Stores the results in MPedestalCam::GetAverageSector(sector)
880//
881void MExtractPedestal::CalcSectorResults(const UInt_t nevts, const UInt_t nspix, const UInt_t sector)
882{
883 const Double_t sum = fSectorSumx[sector];
884 const Double_t sum2 = fSectorSumx2[sector];
885
886 // 1. Calculate the mean of the sums:
887 Double_t ped = sum/nevts;
888
889 // 2. Calculate the Variance of the sums:
890 Double_t var = (sum2/nspix-sum*sum/nevts)/(nevts-1.);
891
892 // 3. Calculate the amplitude of the 150MHz "AB" noise
893 Double_t abOffs = (fSectorSumAB0[sector] - fSectorSumAB1[sector]) / nevts;
894
895 // 4. Scale the mean, variance and AB-noise to the number of slices:
896 ped /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
897 var /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
898 abOffs /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
899 // The pedestal extracted with the extractor is divided by
900 // the number of hi-gain samples because the calibration
901 // multiplies by this number
902
903 // 5. Scale the mean, variance and AB-noise to the number of pixels:
904 ped /= nspix;
905 var /= nspix;
906 abOffs /= nspix;
907
908 // 6. Calculate the RMS from the Variance:
909 const Double_t rms = var<0 ? 0 : TMath::Sqrt(var);
910
911 // abOffs contains only half of the signal as ped.
912 // Therefor abOffs is not the full, but the half amplitude
913 fPedestalsOut->GetAverageSector(sector).Set(ped, rms, abOffs, nevts);
914}
915
916//-----------------------------------------------------------------------
917//
918void MExtractPedestal::Print(Option_t *o) const
919{
920 *fLog << GetDescriptor() << ":" << endl;
921 *fLog << "Name of input MPedestalCam: " << (fPedestalsIn?fPedestalsIn->GetName():fNamePedestalCamIn.Data()) << " (" << fPedestalsIn << ")" << endl;
922 *fLog << "Name of interm. MPedestalCam: " << (fPedestalsInter?fPedestalsInter->GetName():fNamePedestalCamInter.Data()) << " (" << fPedestalsInter << ")" << endl;
923 *fLog << "Name of output MPedestalCam: " << (fPedestalsOut?fPedestalsOut->GetName():fNamePedestalCamOut.Data()) << " (" << fPedestalsOut << ")" << endl;
924 *fLog << "Intermediate Storage is " << (fIntermediateStorage?"on":"off") << endl;
925 *fLog << "Special pixel mode " << (fUseSpecialPixels?"on":"off") << endl;
926 if (fExtractor)
927 {
928 *fLog << "Extractor used: " << fExtractor->ClassName() << " (";
929 *fLog << (fRandomCalculation?"":"non-") << "random)" << endl;
930 }
931 *fLog << "ExtractWindow from slice " << fExtractWinFirst << " to " << fExtractWinLast << " incl." << endl;
932 *fLog << "CheckWindow from slice " << fCheckWinFirst << " to " << fCheckWinLast << " incl." << endl;
933 *fLog << "Max.allowed signal variation: " << fMaxSignalVar << endl;
934}
Note: See TracBrowser for help on using the repository browser.