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

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