source: tags/Mars-V1.0/mpedestal/MExtractPedestal.cc

Last change on this file was 8305, 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.26 2007-02-04 15:49:08 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 Bool_t rc = kTRUE;
271
272 if (windows==0)
273 {
274 *fLog << warn << GetDescriptor();
275 *fLog << " - WARNING: Window size in SetExtractWindow has to be > 0... adjusting to 2!" << endl;
276 windows = 2;
277 rc = kFALSE;
278 }
279
280 fExtractWinSize = windows;
281 fExtractWinFirst = windowf;
282 fExtractWinLast = fExtractWinFirst+fExtractWinSize-1;
283
284 return rc;
285}
286
287// --------------------------------------------------------------------------
288//
289// SetCheckRange:
290//
291// Exits, if the first argument is smaller than 0
292// Exits, if the the last argument is smaller than the first
293//
294Bool_t MExtractPedestal::SetCheckRange(UShort_t chfirst, UShort_t chlast)
295{
296
297 Bool_t rc = kTRUE;
298
299 if (chlast<=chfirst)
300 {
301 *fLog << warn << GetDescriptor();
302 *fLog << " - WARNING: Last slice in SetCheckRange smaller than first slice... set to first+2" << endl;
303 chlast = chfirst+1;
304 rc = kFALSE;
305 }
306
307 fCheckWinFirst = chfirst;
308 fCheckWinLast = chlast;
309
310 return rc;
311}
312
313
314// --------------------------------------------------------------------------
315//
316// Check (and if neccesary: correct) the extraction and check ranges.
317//
318void MExtractPedestal::CheckExtractionWindow(UInt_t offset)
319{
320 *fLog << inf;
321 *fLog << "Requested CheckWindow is [" << fCheckWinFirst << "," << fCheckWinLast << "]." <<endl;
322 *fLog << "Requested ExtractWindow is [" << fExtractWinFirst+offset << "," << fExtractWinLast+offset << "]." <<endl;
323
324 // fSignal->GetNumSamples() not yet initialized!!!
325 const UInt_t num = fRunHeader->GetNumSamplesHiGain()+fRunHeader->GetNumSamplesLoGain();
326
327 if (fCheckWinLast >= num)
328 {
329 *fLog << "CheckWindow [" << fCheckWinFirst << "," << fCheckWinLast;
330 *fLog << "] out of range [0," << num-1 << "]... ";
331 *fLog << "reset upper edge." << endl;
332
333 fCheckWinLast = num-1;
334 }
335
336 if (offset+fExtractWinLast >= num)
337 {
338 *fLog << "ExtractWindow [" << fExtractWinFirst+offset << "," << fExtractWinLast+offset;
339 *fLog << "] out of range [0," << num-1 << "]... ";
340 *fLog << "reset upper edge." << endl;
341
342 fExtractWinLast = num-offset-1;
343 }
344
345 fExtractWinSize = fExtractWinLast-fExtractWinFirst+1;
346
347 if (fExtractor || TMath::Even(fExtractWinSize))
348 return;
349
350 fExtractWinLast += offset+fExtractWinLast==num-1 ? -1 : +1;
351
352 *fLog << "ExtractionWindow odd... set to [";
353 *fLog << fExtractWinFirst+offset << "," << fExtractWinLast+offset << "]" << endl;
354
355 fExtractWinSize = fExtractWinLast-fExtractWinFirst+1;
356}
357
358// --------------------------------------------------------------------------
359//
360// Look for the following input containers:
361//
362// - MRawEvtData
363// - MRawRunHeader
364// - MRawEvtHeader
365// - MGeomCam
366//
367// The following output containers are also searched and created if
368// they were not found:
369//
370// - MPedestalCam with the name fPedContainerName
371//
372Int_t MExtractPedestal::PreProcess(MParList *pList)
373{
374
375 Clear();
376
377 fRawEvt = (MRawEvtData*)pList->FindObject(AddSerialNumber(fNameRawEvtData));
378 if (!fRawEvt)
379 if (!fUseSpecialPixels)
380 {
381 *fLog << err << AddSerialNumber(fNameRawEvtData) << " not found... aborting." << endl;
382 return kFALSE;
383 }
384
385
386 fRawEvt = (MRawEvtData*)pList->FindObject(AddSerialNumber("MRawEvtData"));
387 if (!fRawEvt)
388 {
389 *fLog << err << AddSerialNumber("MRawEvtData") << " not found... aborting." << endl;
390 return kFALSE;
391 }
392
393 fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
394 if (!fRunHeader)
395 {
396 *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl;
397 return kFALSE;
398 }
399
400 fGeom = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
401 if (!fGeom)
402 {
403 *fLog << err << AddSerialNumber("MGeomCam") << " not found... aborting." << endl;
404 return kFALSE;
405 }
406
407 fSignal = (MPedestalSubtractedEvt*)pList->FindObject(AddSerialNumber("MPedestalSubtractedEvt"));
408 if (!fSignal)
409 {
410 *fLog << err << AddSerialNumber("MPedestalSubtractedEvt") << " not found... aborting." << endl;
411 return kFALSE;
412 }
413
414 if (fExtractor && !fPedestalsIn)
415 {
416 fPedestalsIn = (MPedestalCam*)pList->FindObject(AddSerialNumber(fNamePedestalCamIn), "MPedestalCam");
417 if (!fPedestalsIn)
418 {
419 *fLog << err << AddSerialNumber(fNamePedestalCamIn) << " not found... aborting." << endl;
420 return kFALSE;
421 }
422 }
423
424 if (!fPedestalsInter && fIntermediateStorage)
425 {
426 fPedestalsInter = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", AddSerialNumber(fNamePedestalCamInter));
427 if (!fPedestalsInter)
428 return kFALSE;
429 }
430
431 if (!fPedestalsOut)
432 {
433 fPedestalsOut = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", AddSerialNumber(fNamePedestalCamOut));
434 if (!fPedestalsOut)
435 return kFALSE;
436 }
437
438 return fExtractor ? fExtractor->CallPreProcess(pList) : kTRUE;
439}
440
441//-----------------------------------------------------------------------
442//
443// Call Calc(). If fExtractor!=NULL enclose the call in setting the
444// NoiseCalculation to fRandomCalculation
445//
446Int_t MExtractPedestal::Process()
447{
448 //
449 // Necessary check for extraction of special pixels
450 // together with data which does not yet have them
451 //
452 if (fSumx.GetSize()==0)
453 return kTRUE;
454
455 if (fExtractor)
456 fExtractor->SetNoiseCalculation(fRandomCalculation);
457
458 const Int_t rc = Calc();
459
460 if (fExtractor)
461 fExtractor->SetNoiseCalculation(kFALSE);
462
463 return rc;
464}
465
466// ---------------------------------------------------------------------------------
467//
468// Sets the size (from MPedestalCam::GetSize() ) and resets the following arrays:
469// - fSumx
470// - fSumx2
471// - fSumAB0
472// - fSumAB1
473// - fAreaSumx
474// - fAreaSumx2
475// - fAreaSumAB0
476// - fAreaSumAB1
477// - fAreaFilled
478// - fAreaValid
479// - fSectorSumx
480// - fSectorSumx2
481// - fSectorSumAB0
482// - fSectorSumAB1
483// - fSectorFilled
484// - fSectorValid
485//
486Bool_t MExtractPedestal::ReInit(MParList *pList)
487{
488 // Necessary check for special pixels which might not yet have existed
489 if (!fRawEvt)
490 {
491 if (fRunHeader->GetFormatVersion() > 3)
492 return kTRUE;
493
494 *fLog << err << "ERROR - " << fNameRawEvtData << " [MRawEvtData] has not ";
495 *fLog << "been found and format version > 3... abort." << endl;
496 return kFALSE;
497 }
498
499 // If the size is not yet set, set the size
500 if (fSumx.GetSize()==0)
501 {
502 // Initialize the normal pixels (size of MPedestalCam already set by MGeomApply)
503 const Int_t npixels = fPedestalsOut->GetSize();
504
505 fSumx. Set(npixels);
506 fSumx2. Set(npixels);
507 fSumAB0.Set(npixels);
508 fSumAB1.Set(npixels);
509
510 if (fUseSpecialPixels)
511 {
512 // Initialize size of MPedestalCam in case of special pixels (not done by MGeomApply)
513 const UShort_t nspecial = fRunHeader->GetNumSpecialPixels();
514 if (nspecial == 0)
515 {
516 *fLog << warn << "WARNING - Number of special pixels is 0." << endl;
517 return kTRUE;
518 }
519
520 fPedestalsOut->InitSize((UInt_t)nspecial);
521 }
522 else
523 {
524 // Initialize the averaged areas and sectors (do not exist for special pixels)
525 const Int_t areas = fPedestalsOut->GetNumAverageArea();
526 const Int_t sectors = fPedestalsOut->GetNumAverageSector();
527
528 fAreaSumx. Set(areas);
529 fAreaSumx2. Set(areas);
530 fAreaSumAB0.Set(areas);
531 fAreaSumAB1.Set(areas);
532 fAreaFilled.Set(areas);
533 fAreaValid .Set(areas);
534
535 fSectorSumx. Set(sectors);
536 fSectorSumx2. Set(sectors);
537 fSectorSumAB0.Set(sectors);
538 fSectorSumAB1.Set(sectors);
539 fSectorFilled.Set(sectors);
540 fSectorValid .Set(sectors);
541
542 for (Int_t i=0; i<npixels; i++)
543 {
544 const UInt_t aidx = (*fGeom)[i].GetAidx();
545 const UInt_t sector = (*fGeom)[i].GetSector();
546
547 fAreaValid [aidx] ++;
548 fSectorValid[sector]++;
549 }
550 }
551 }
552
553 if (fExtractor)
554 {
555 if (!fExtractor->ReInit(pList))
556 return kFALSE;
557
558 // If an extractor is set determin the window size automatically!
559 fExtractWinFirst = fExtractor->GetHiGainFirst();
560 fExtractWinLast = fExtractor->GetHiGainLast();
561
562 // fSignal->GetNumSamples() not yet initialized!!!
563 const UInt_t num = fRunHeader->GetNumSamplesHiGain()+fRunHeader->GetNumSamplesLoGain();
564 if (fExtractWinLast >= num)
565 {
566 *fLog << err << GetDescriptor();
567 *fLog << " - ERROR: Selected fExtractWinLast " << fExtractWinLast;
568 *fLog << " out of range (>=" << num<< ")." << endl;
569 return kFALSE;
570 }
571 }
572
573 //CheckExtractionWindow();
574
575 return kTRUE;
576}
577
578// ---------------------------------------------------------------------------------
579//
580// PostProcess the extractor if available
581//
582Int_t MExtractPedestal::PostProcess()
583{
584 fPedestalsIn = NULL;
585 return fExtractor ? fExtractor->CallPostProcess() : kTRUE;
586}
587
588// ---------------------------------------------------------------------------------
589//
590// Check whether the signal variation between fCheckWinFirst and fCheckWinLast
591// exceeds fMaxSignalVar or the signal is greater than 250
592//
593Bool_t MExtractPedestal::CheckVariation(UInt_t idx) const
594{
595 // This is the fast workaround to put hi- and lo-gains together
596 Byte_t *slices = fSignal->GetSamplesRaw(idx);//pixel.GetSamples();
597
598 // Start 'real' work
599 UShort_t max = 0;
600 UShort_t min = (UShort_t)-1;
601
602 // Find the maximum and minimum signal per slice in the high gain window
603 for (Byte_t *slice=slices+fCheckWinFirst; slice<=slices+fCheckWinLast; slice++)
604 {
605 if (*slice > max)
606 max = *slice;
607 if (*slice < min)
608 min = *slice;
609 }
610
611 // If the maximum in the high gain window is smaller than
612 // FIXME: Precompiled value!
613 return max-min<fMaxSignalVar && max<250;
614}
615
616// ---------------------------------------------------------------------------------
617//
618// Invoke the hi-gain extraction starting at fExtractWinFirst+offset
619// for fExtractWinLast-fExtractWinFirst+1 slices. If Noise calculation
620// is set it is up to the signal extractor to do the right thing.
621//
622// Returns the extracted signal.
623//
624Float_t MExtractPedestal::CalcExtractor(const MRawEvtPixelIter &pixel, Int_t offset) const
625{
626 // Use the same extraction window as for signal extraction
627 const Int_t first = fExtractWinFirst;
628 const Int_t last = fExtractWinLast;
629
630 const Int_t start = first+offset;
631
632 const Int_t range = last-first+1;
633
634 // This check is already done in CheckExtractionWindow
635 // if (range>pixel.GetNumSamples()-start)
636 // range = pixel.GetNumSamples()-start;
637
638 const Int_t idx = pixel.GetPixelId();
639
640 // Do some handling if maxpos is last slice?
641 const Int_t maxposhi = fSignal->GetMax(idx, start, start+range-1);
642
643 const Float_t *sig = fSignal->GetSamples(idx);
644
645 // The pedestal is extracted with the hi-gain extractor (eg. digital
646 // filter weights) but from the lo-gains
647 Float_t dummy[3];
648 Float_t sum = 0;
649 fExtractor->FindTimeAndChargeHiGain2(sig+start, range, sum,
650 dummy[0], dummy[1], dummy[2],
651 0, maxposhi);
652 return sum;
653}
654
655// ---------------------------------------------------------------------------------
656//
657// Sum slices from fExtractWinFirst to fExtractWinLast. The total sum is
658// returned. ab0 and ab1 will contain the total sum splitted by the
659// AB-flag. If the AB-flag is invalid ab0=ab1=0 is returned.
660//
661UInt_t MExtractPedestal::CalcSums(const MRawEvtPixelIter &pixel, Int_t offset, UInt_t &ab0, UInt_t &ab1) const
662{
663 const Int_t first = fExtractWinFirst+offset;
664
665 Byte_t *ptr = fSignal->GetSamplesRaw(pixel.GetPixelId())+first;
666 Byte_t *end = ptr + fExtractWinSize;
667
668 Int_t abflag = pixel.HasABFlag() + first;
669
670 UInt_t ab[2] = { 0, 0 };
671 while (ptr<end)
672 ab[abflag++ & 0x1] += *ptr++;
673
674 // This check if for old data without AB-Flag in the data
675 const Bool_t valid = pixel.IsABFlagValid();
676
677 ab0 = valid ? ab[0] : 0;
678 ab1 = valid ? ab[1] : 0;
679
680 return ab[0]+ab[1];
681}
682
683// --------------------------------------------------------------------------
684//
685// The following resources are available:
686// ExtractWindowFirst: 15
687// ExtractWindowSize: 6
688// PedestalUpdate: yes
689// RandomCalculation: yes
690//
691Int_t MExtractPedestal::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
692{
693 Bool_t rc=kFALSE;
694
695 // find resource for fUseSpecialPixels
696 if (IsEnvDefined(env, prefix, "UseSpecialPixels", print))
697 {
698 SetUseSpecialPixels(GetEnvValue(env, prefix, "UseSpecialPixels", fUseSpecialPixels));
699 rc = kTRUE;
700 }
701
702 if (IsEnvDefined(env, prefix, "IntermediateStorage", print))
703 {
704 SetIntermediateStorage(GetEnvValue(env, prefix, "IntermediateStorage", fIntermediateStorage));
705 rc = kTRUE;
706 }
707
708 // find resource for random calculation
709 if (IsEnvDefined(env, prefix, "RandomCalculation", print))
710 {
711 SetRandomCalculation(GetEnvValue(env, prefix, "RandomCalculation", fRandomCalculation));
712 rc = kTRUE;
713 }
714
715 // Find resources for ExtractWindow
716 Int_t ef = fExtractWinFirst;
717 Int_t es = fExtractWinSize;
718 if (IsEnvDefined(env, prefix, "ExtractWinFirst", print))
719 {
720 ef = GetEnvValue(env, prefix, "ExtractWinFirst", ef);
721 rc = kTRUE;
722 }
723 if (IsEnvDefined(env, prefix, "ExtractWinSize", print))
724 {
725 es = GetEnvValue(env, prefix, "ExtractWinSize", es);
726 rc = kTRUE;
727 }
728
729 SetExtractWindow(ef,es);
730
731 // Find resources for CheckWindow
732 Int_t cfs = fCheckWinFirst;
733 Int_t cls = fCheckWinLast;
734 if (IsEnvDefined(env, prefix, "CheckWinFirst", print))
735 {
736 cfs = GetEnvValue(env, prefix, "CheckWinFirst", cfs);
737 rc = kTRUE;
738 }
739 if (IsEnvDefined(env, prefix, "CheckWinLast", print))
740 {
741 cls = GetEnvValue(env, prefix, "CheckWinLast", cls);
742 rc = kTRUE;
743 }
744
745 SetCheckRange(cfs,cls);
746
747 // find resource for maximum signal variation
748 if (IsEnvDefined(env, prefix, "MaxSignalVar", print))
749 {
750 SetMaxSignalVar(GetEnvValue(env, prefix, "MaxSignalVar", fMaxSignalVar));
751 rc = kTRUE;
752 }
753
754 // find resource for MPedestalCam
755 if (IsEnvDefined(env, prefix, "NamePedestalCamIn", print))
756 {
757 SetNamePedestalCamIn(GetEnvValue(env, prefix, "NamePedestalCamIn", fNamePedestalCamIn));
758 rc = kTRUE;
759 }
760
761 if (IsEnvDefined(env, prefix, "NamePedestalCamInter", print))
762 {
763 SetNamePedestalCamInter(GetEnvValue(env, prefix, "NamePedestalCamInter", fNamePedestalCamInter));
764 rc = kTRUE;
765 }
766
767 if (IsEnvDefined(env, prefix, "NamePedestalCamOut", print))
768 {
769 SetNamePedestalCamOut(GetEnvValue(env, prefix, "NamePedestalCamOut", fNamePedestalCamOut));
770 rc = kTRUE;
771 }
772
773 return rc;
774}
775
776// ---------------------------------------------------------------------------------
777//
778// Calculates for pixel "idx":
779//
780// Ped per slice = sum / n / fExtractWinSize;
781// RMS per slice = sqrt { (sum2 - sum*sum/n) / (n-1) / fExtractWinSize }
782// ABOffset per slice = (fSumAB0[idx] - fSumAB1[idx]) / n / fExtractWinSize;
783//
784// Stores the results in MPedestalCam[pixid]
785//
786void MExtractPedestal::CalcPixResults(const UInt_t nevts, const UInt_t pixid)
787{
788 const Double_t sum = fSumx[pixid];
789 const Double_t sum2 = fSumx2[pixid];
790
791 // 1. Calculate the mean of the sums:
792 Double_t ped = sum/nevts;
793
794 // 2. Calculate the Variance of the sums:
795 Double_t var = (sum2-sum*sum/nevts)/(nevts-1.);
796
797 // 3. Calculate the amplitude of the 150MHz "AB" noise
798 Double_t abOffs = (fSumAB0[pixid] - fSumAB1[pixid]) / nevts;
799
800 // 4. Scale the mean, variance and AB-noise to the number of slices:
801 ped /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
802 var /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
803 abOffs /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
804 // The pedestal extracted with the extractor is divided by
805 // the number of hi-gain samples because the calibration
806 // multiplies by this number
807
808 // 5. Calculate the RMS from the Variance:
809 const Double_t rms = var<0 ? 0 : TMath::Sqrt(var);
810
811 // abOffs contains only half of the signal as ped.
812 // Therefor abOffs is not the full, but the half amplitude
813 (*fPedestalsOut)[pixid].Set(ped, rms, abOffs, nevts);
814}
815
816// ---------------------------------------------------------------------------------
817//
818// Calculates for area idx "aidx" with "napix" valid pixels:
819//
820// Ped per slice = sum / nevts / fExtractWinSize / napix;
821// RMS per slice = sqrt { (sum2 - sum*sum/nevts) / (nevts-1) / fExtractWinSize / napix }
822// ABOffset per slice = (fSumAB0[idx] - fSumAB1[idx]) / nevts / fExtractWinSize / napix;
823//
824// Stores the results in MPedestalCam::GetAverageArea(aidx)
825//
826void MExtractPedestal::CalcAreaResults(const UInt_t nevts, const UInt_t napix, const UInt_t aidx)
827{
828 const Double_t sum = fAreaSumx[aidx];
829 const Double_t sum2 = fAreaSumx2[aidx];
830
831 // 1. Calculate the mean of the sums:
832 Double_t ped = sum/nevts;
833
834 // 2. Calculate the Variance of the sums:
835 Double_t var = (sum2/napix-sum*sum/nevts)/(nevts-1.);
836
837 // 3. Calculate the amplitude of the 150MHz "AB" noise
838 Double_t abOffs = (fAreaSumAB0[aidx] - fAreaSumAB1[aidx]) / nevts;
839
840 // 4. Scale the mean, variance and AB-noise to the number of slices:
841 ped /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
842 var /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
843 abOffs /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
844 // The pedestal extracted with the extractor is divided by
845 // the number of hi-gain samples because the calibration
846 // multiplies by this number
847
848 // 5. Scale the mean, variance and AB-noise to the number of pixels:
849 ped /= napix;
850 var /= napix;
851 abOffs /= napix;
852
853 // 6. Calculate the RMS from the Variance:
854 const Double_t rms = var<0 ? 0 : TMath::Sqrt(var);
855
856 // abOffs contains only half of the signal as ped.
857 // Therefor abOffs is not the full, but the half amplitude
858 fPedestalsOut->GetAverageArea(aidx).Set(ped, rms, abOffs, nevts);
859}
860
861// ---------------------------------------------------------------------------------
862//
863// Calculates for sector idx "sector" with "nspix" valid pixels:
864//
865// Ped per slice = sum / nevts / fExtractWinSize / nspix;
866// RMS per slice = sqrt { (sum2 - sum*sum/nevts) / (nevts-1) / fExtractWinSize / nspix }
867// ABOffset per slice = (fSumAB0[idx] - fSumAB1[idx]) / nevts / fExtractWinSize / nspix;
868//
869// Stores the results in MPedestalCam::GetAverageSector(sector)
870//
871void MExtractPedestal::CalcSectorResults(const UInt_t nevts, const UInt_t nspix, const UInt_t sector)
872{
873 const Double_t sum = fSectorSumx[sector];
874 const Double_t sum2 = fSectorSumx2[sector];
875
876 // 1. Calculate the mean of the sums:
877 Double_t ped = sum/nevts;
878
879 // 2. Calculate the Variance of the sums:
880 Double_t var = (sum2/nspix-sum*sum/nevts)/(nevts-1.);
881
882 // 3. Calculate the amplitude of the 150MHz "AB" noise
883 Double_t abOffs = (fSectorSumAB0[sector] - fSectorSumAB1[sector]) / nevts;
884
885 // 4. Scale the mean, variance and AB-noise to the number of slices:
886 ped /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
887 var /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
888 abOffs /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize;
889 // The pedestal extracted with the extractor is divided by
890 // the number of hi-gain samples because the calibration
891 // multiplies by this number
892
893 // 5. Scale the mean, variance and AB-noise to the number of pixels:
894 ped /= nspix;
895 var /= nspix;
896 abOffs /= nspix;
897
898 // 6. Calculate the RMS from the Variance:
899 const Double_t rms = var<0 ? 0 : TMath::Sqrt(var);
900
901 // abOffs contains only half of the signal as ped.
902 // Therefor abOffs is not the full, but the half amplitude
903 fPedestalsOut->GetAverageSector(sector).Set(ped, rms, abOffs, nevts);
904}
905
906//-----------------------------------------------------------------------
907//
908void MExtractPedestal::Print(Option_t *o) const
909{
910 *fLog << GetDescriptor() << ":" << endl;
911 *fLog << "Name of input MPedestalCam: " << (fPedestalsIn?fPedestalsIn->GetName():fNamePedestalCamIn.Data()) << " (" << fPedestalsIn << ")" << endl;
912 *fLog << "Name of interm. MPedestalCam: " << (fPedestalsInter?fPedestalsInter->GetName():fNamePedestalCamInter.Data()) << " (" << fPedestalsInter << ")" << endl;
913 *fLog << "Name of output MPedestalCam: " << (fPedestalsOut?fPedestalsOut->GetName():fNamePedestalCamOut.Data()) << " (" << fPedestalsOut << ")" << endl;
914 *fLog << "Intermediate Storage is " << (fIntermediateStorage?"on":"off") << endl;
915 *fLog << "Special pixel mode " << (fUseSpecialPixels?"on":"off") << endl;
916 if (fExtractor)
917 {
918 *fLog << "Extractor used: " << fExtractor->ClassName() << " (";
919 *fLog << (fRandomCalculation?"":"non-") << "random)" << endl;
920 }
921 *fLog << "ExtractWindow from slice " << fExtractWinFirst << " to " << fExtractWinLast << " incl." << endl;
922 *fLog << "CheckWindow from slice " << fCheckWinFirst << " to " << fCheckWinLast << " incl." << endl;
923 *fLog << "Max.allowed signal variation: " << fMaxSignalVar << endl;
924}
Note: See TracBrowser for help on using the repository browser.