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

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