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

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