source: tags/Mars-V0.9.2/mpedestal/MPedCalcFromLoGain.cc

Last change on this file was 6325, checked in by gaug, 20 years ago
*** empty log message ***
File size: 16.4 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Josep Flix 04/2001 <mailto:jflix@ifae.es>
19! Author(s): Thomas Bretz 05/2001 <mailto:tbretz@astro.uni-wuerzburg.de>
20! Author(s): Sebastian Commichau 12/2003
21! Author(s): Javier Rico 01/2004 <mailto:jrico@ifae.es>
22! Author(s): Markus Gaug 01/2004 <mailto:markus@ifae.es>
23! Author(s): Florian Goebel 06/2004 <mailto:fgoebel@mppmu.mpg.de>
24! Author(s): Nepomuk Otte 10/2004 <mailto:otte@mppmu.mpg.de>
25!
26! Copyright: MAGIC Software Development, 2000-2004
27!
28!
29\* ======================================================================== */
30
31/////////////////////////////////////////////////////////////////////////////
32//
33// MPedCalcLoGain
34//
35// This task derives from MExtractPedestal.
36// It calculates the pedestals using the low gain slices, whenever the difference
37// between the highest and the lowest slice in the high gain
38// slices is below a given threshold (SetMaxHiGainVar). In this case the receiver
39// boards do not switch to lo gain and the so called lo gain slices are actually
40// high gain slices.
41//
42// MPedCalcLoGain also fills the ABoffset in MPedestalPix which allows to correct
43// the 150 MHz clock noise.
44//
45// This task takes a pedestal run file and fills MPedestalCam during
46// the Process() with the pedestal and rms computed in an event basis.
47// In the PostProcess() MPedestalCam is finally filled with the pedestal
48// mean and rms computed in a run basis.
49// More than one run (file) can be merged
50//
51// MPedCalcPedRun applies the following formula (1):
52//
53// Pedestal per slice = sum(x_i) / n / slices
54// PedRMS per slice = Sqrt( ( sum(x_i^2) - sum(x_i)^2/n ) / n-1 / slices )
55//
56// where x_i is the sum of "slices" FADC slices and sum means the sum over all
57// events. "n" is the number of events, "slices" is the number of summed FADC samples.
58//
59// Note that the slice-to-slice fluctuations are not Gaussian, but Poissonian, thus
60// asymmetric and they are correlated.
61//
62// It is important to know that the Pedestal per slice and PedRMS per slice depend
63// on the number of used FADC slices, as seen in the following plots:
64//
65//Begin_Html
66/*
67<img src="images/PedestalStudyInner.gif">
68*/
69//End_Html
70//
71//Begin_Html
72/*
73<img src="images/PedestalStudyOuter.gif">
74*/
75//End_Html
76//
77// The plots show the inner and outer pixels, respectivly and have the following meaning:
78//
79// 1) The calculated mean pedestal per slice (from MPedCalcFromLoGain)
80// 2) The fitted mean pedestal per slice (from MHPedestalCam)
81// 3) The calculated pedestal RMS per slice (from MPedCalcFromLoGain)
82// 4) The fitted sigma of the pedestal distribution per slice
83// (from MHPedestalCam)
84// 5) The relative difference between calculation and histogram fit
85// for the mean
86// 6) The relative difference between calculation and histogram fit
87// for the sigma or RMS, respectively.
88//
89// The calculated means do not change significantly except for the case of 2 slices,
90// however the RMS changes from 5.7 per slice in the case of 2 extracted slices
91// to 8.3 per slice in the case of 26 extracted slices. This change is very significant.
92//
93// The plots have been produced on run 20123. You can reproduce them using
94// the macro pedestalstudies.C
95//
96// Usage of this class:
97// ====================
98//
99//
100// fCheckWinFirst = fgCheckWinFirst = 0
101// fCheckWinLast = fgCheckWinLast = 29
102// fExtractWinFirst = fgExtractWinFirst = 17
103// fExtractWinSize = fgExtractWinSize = 6
104// fMaxSignalVar = fgMaxSignalVar = 40;
105//
106// Call:
107// SetCheckRange(fCheckWinFirst,fCheckWinLast);
108// to set the Window in which a signal is searched
109//
110// SetExtractWindow(fExtractWinFirst,fExtractWinSize);
111// to set the Window from which a signal is extracted
112//
113// SetMaxSignalVar(fMaxSignalVar);
114// set the maximum allowed difference between maximum and minimal signal in CheckWindow
115//
116// Variables:
117// fgCheckWinFirst; First FADC slice to check for signal (currently set to: 0)
118// fgCheckWinLast: Last FADC slice to check for signal (currently set to: 29)
119// fgExtractWinFirst: First FADC slice to be used for pedestal extraction (currently set to: 15)
120// fgExtractWinSize: Window size in slices used for pedestal extraction (currently set to: 6)
121// fgMaxSignalVar: The maximum difference between the highest and lowest slice
122// in the check window allowed in order to use event
123//
124// Input Containers:
125// MRawEvtData
126// MRawRunHeader
127// MGeomCam
128//
129// Output Containers:
130// MPedestalCam
131//
132// See also: MPedestalCam, MPedestalPix, MHPedestalCam, MExtractor
133//
134/////////////////////////////////////////////////////////////////////////////
135#include "MPedCalcFromLoGain.h"
136
137#include "MParList.h"
138
139#include "MLog.h"
140#include "MLogManip.h"
141
142#include "MRawRunHeader.h"
143#include "MRawEvtPixelIter.h"
144#include "MRawEvtData.h"
145
146#include "MPedestalPix.h"
147#include "MPedestalCam.h"
148
149#include "MGeomPix.h"
150#include "MGeomCam.h"
151
152#include "MExtractPedestal.h"
153#include "MExtractTimeAndCharge.h"
154
155ClassImp(MPedCalcFromLoGain);
156
157using namespace std;
158
159const UShort_t MPedCalcFromLoGain::fgCheckWinFirst = 0;
160const UShort_t MPedCalcFromLoGain::fgCheckWinLast = 29;
161const UShort_t MPedCalcFromLoGain::fgExtractWinFirst = 17;
162const UShort_t MPedCalcFromLoGain::fgExtractWinSize = 6;
163const UShort_t MPedCalcFromLoGain::fgMaxSignalVar = 40;
164
165// --------------------------------------------------------------------------
166//
167// Default constructor:
168//
169// Sets:
170// - all pointers to NULL
171//
172// Calls:
173// - AddToBranchList("fHiGainPixId");
174// - AddToBranchList("fHiGainFadcSamples");
175// - SetCheckRange(fgCheckWinFirst, fgCheckWinLast, fgExtractWinFirst, fgExtractWinSize)
176//
177MPedCalcFromLoGain::MPedCalcFromLoGain(const char *name, const char *title)
178{
179 fName = name ? name : "MPedCalcFromLoGain";
180 fTitle = title ? title : "Task to calculate pedestals from lo-gains";
181
182 SetCheckRange(fgCheckWinFirst, fgCheckWinLast);
183 SetExtractWindow(fgExtractWinFirst, fgExtractWinSize);
184
185 SetMaxSignalVar(fgMaxSignalVar);
186}
187
188void MPedCalcFromLoGain::ResetArrays()
189{
190 MExtractPedestal::ResetArrays();
191
192 fNumEventsUsed.Reset();
193 fTotalCounter.Reset();
194}
195
196// --------------------------------------------------------------------------
197//
198// SetCheckRange:
199//
200// Exits, if the first argument is smaller than 0
201// Exits, if the the last argument is smaller than the first
202//
203Bool_t MPedCalcFromLoGain::SetCheckRange(UShort_t chfirst, UShort_t chlast)
204{
205
206 Bool_t rc = kTRUE;
207
208 if (chlast<=chfirst)
209 {
210 *fLog << warn << GetDescriptor();
211 *fLog << " - WARNING: Last slice in SetCheckRange smaller than first slice... set to first+2" << endl;
212 chlast = chfirst+1;
213 rc = kFALSE;
214 }
215
216 fCheckWinFirst = chfirst;
217 fCheckWinLast = chlast;
218
219 return rc;
220}
221
222// ---------------------------------------------------------------------------------
223//
224// Checks:
225// - if the number of available slices
226// (fRunHeader->GetNumSamplesHiGain()+(Int_t)fRunHeader->GetNumSamplesLoGain()-1)
227// is smaller than the number of used slices
228// (fExtractWinSize+ fExtractWinFirst-1) or
229// fCheckWinLast
230//
231Bool_t MPedCalcFromLoGain::ReInit(MParList *pList)
232{
233
234 const UShort_t hisamples = fRunHeader->GetNumSamplesHiGain();
235 const UShort_t losamples = fRunHeader->GetNumSamplesLoGain();
236
237 fSlices.Set(hisamples+losamples);
238
239 UShort_t lastavailable = hisamples+losamples-1;
240
241 if (fExtractor)
242 fExtractWinLast = fExtractWinFirst + fExtractor->GetWindowSizeHiGain() - 1;
243
244 // If the size is not yet set, set the size
245 if (fSumx.GetSize()==0)
246 {
247 const Int_t npixels = fPedestalsOut->GetSize();
248 fNumEventsUsed.Set(npixels);
249 fTotalCounter.Set(npixels);
250 }
251
252 if (fCheckWinLast > lastavailable) //changed to override check
253 {
254 *fLog << warn << GetDescriptor();
255 *fLog << " - WARNING: Last Check Window slice out of range...adjusting to last available slice ";
256 *fLog << lastavailable << endl;
257
258 fCheckWinLast = lastavailable;
259 }
260
261 if (fExtractWinLast > lastavailable)
262 {
263 if (fExtractor)
264 {
265 *fLog << err << GetDescriptor();
266 *fLog << " - ERROR: Selected Last Extraction Window: "
267 << fExtractWinFirst + fExtractor->GetWindowSizeHiGain()-1
268 << "ranges out of range: " << lastavailable-1 << endl;
269 return kFALSE;
270 }
271 else
272 {
273 const UShort_t diff = fExtractWinLast - lastavailable;
274 *fLog << warn << GetDescriptor();
275 *fLog << " - WARNING: Selected Extract Window ranges out of range...adjusting to last available slice ";
276 *fLog << lastavailable << endl;
277
278 fExtractWinLast -= diff;
279 fExtractWinSize -= diff;
280 }
281 }
282
283 return MExtractPedestal::ReInit(pList);
284}
285
286// --------------------------------------------------------------------------
287//
288// Fill the MPedestalCam container with the signal mean and rms for the event.
289// Store the measured signal in arrays fSumx and fSumx2 so that we can
290// calculate the overall mean and rms in the PostProcess()
291//
292Int_t MPedCalcFromLoGain::Calc()
293{
294 // This is the workaround to put hi- and lo-gains together
295 const Int_t nhigain = fRunHeader->GetNumSamplesHiGain();
296 const Int_t nlogain = fRunHeader->GetNumSamplesLoGain();
297
298 Byte_t *slices = fSlices.GetArray();
299
300 // Real Process
301 MRawEvtPixelIter pixel(fRawEvt);
302
303 while (pixel.Next())
304 {
305 // This is the fast workaround to put hi- and lo-gains together
306 memcpy(slices, pixel.GetHiGainSamples(), nhigain);
307 memcpy(slices+nhigain, pixel.GetLoGainSamples(), nlogain);
308
309 // Start 'real' work
310 const UInt_t idx = pixel.GetPixelId();
311
312 const UInt_t aidx = (*fGeom)[idx].GetAidx();
313 const UInt_t sector = (*fGeom)[idx].GetSector();
314
315 UShort_t max = 0;
316 UShort_t min = (UShort_t)-1;
317
318 // Find the maximum and minimum signal per slice in the high gain window
319 for (Byte_t *slice=slices+fCheckWinFirst; slice<=slices+fCheckWinLast; slice++)
320 {
321 if (*slice > max)
322 max = *slice;
323 if (*slice < min)
324 min = *slice;
325 }
326
327 // If the maximum in the high gain window is smaller than
328 if (max-min>=fMaxSignalVar || max>=250)
329 continue;
330
331 Float_t sum = 0.;
332
333 //extract pedestal
334 if (fExtractor)
335 CalcExtractor(pixel, sum, (*fPedestalsIn)[idx]);
336 else
337 {
338 UInt_t sumi = 0;
339 for(Byte_t *slice=slices+fExtractWinFirst; slice<=slices+fExtractWinLast; slice++)
340 sumi += *slice;
341 sum = (Float_t)sumi;
342 }
343
344 const Float_t sqrsum = sum*sum;
345
346 fSumx[idx] += sum;
347 fSumx2[idx] += sqrsum;
348 fAreaSumx[aidx] += sum;
349 fAreaSumx2[aidx] += sqrsum;
350 fSectorSumx[sector] += sum;
351 fSectorSumx2[sector] += sqrsum;
352
353 if (fIntermediateStorage)
354 (*fPedestalsInter)[idx].Set(sum,0.,0.,fNumEventsUsed[idx]);
355
356 fNumEventsUsed[idx] ++;
357 fAreaFilled [aidx] ++;
358 fSectorFilled [sector]++;
359
360 if (!fExtractor)
361 {
362 //
363 // Calculate the amplitude of the 150MHz "AB" noise
364 //
365 const UShort_t abFlag = (fExtractWinFirst + pixel.HasABFlag()) & 0x1;
366
367 for (Byte_t *slice=slices+fExtractWinFirst; slice<=slices+fExtractWinLast; slice+=2)
368 {
369 const UShort_t ab0 = *(slice + abFlag);
370 const UShort_t ab1 = *(slice - abFlag + 1);
371
372 fSumAB0[idx] += ab0;
373 fSumAB1[idx] += ab1;
374 fAreaSumAB0[aidx] += ab0;
375 fAreaSumAB1[aidx] += ab1;
376 fSectorSumAB0[aidx] += ab0;
377 fSectorSumAB1[aidx] += ab1;
378 }
379 }
380
381 if (!fPedestalUpdate || (UInt_t)fNumEventsUsed[idx]<fNumEventsDump)
382 continue;
383
384 CalcPixResults(fNumEventsDump, idx);
385 fTotalCounter[idx]++;
386
387 fNumEventsUsed[idx]=0;
388 fSumx[idx]=0;
389 fSumx2[idx]=0;
390 fSumAB0[idx]=0;
391 fSumAB1[idx]=0;
392 }
393
394 if (!(GetNumExecutions() % fNumAreasDump))
395 CalcAreaResult();
396
397 if (!(GetNumExecutions() % fNumSectorsDump))
398 CalcSectorResult();
399
400 if (fPedestalUpdate)
401 fPedestalsOut->SetReadyToSave();
402
403 return kTRUE;
404}
405
406// --------------------------------------------------------------------------
407//
408// Loop over the sector indices to get the averaged pedestal per sector
409//
410void MPedCalcFromLoGain::CalcSectorResult()
411{
412 for (UInt_t sector=0; sector<fSectorFilled.GetSize(); sector++)
413 if (fSectorValid[sector]>0)
414 CalcSectorResults(fSectorFilled[sector], fSectorValid[sector], sector);
415}
416
417// --------------------------------------------------------------------------
418//
419// Loop over the (two) area indices to get the averaged pedestal per aidx
420//
421void MPedCalcFromLoGain::CalcAreaResult()
422{
423 for (UInt_t aidx=0; aidx<fAreaFilled.GetSize(); aidx++)
424 if (fAreaValid[aidx]>0)
425 CalcAreaResults(fAreaFilled[aidx], fAreaValid[aidx], aidx);
426
427}
428
429void MPedCalcFromLoGain::CalcExtractor(const MRawEvtPixelIter &pixel, Float_t &sum, MPedestalPix &ped)
430{
431 Byte_t sat = 0;
432 Byte_t *logain = pixel.GetLoGainSamples() + fExtractWinFirst;
433
434 const Bool_t logainabflag = (pixel.HasABFlag() + pixel.GetNumHiGainSamples()) & 0x1;
435
436 Float_t dummy;
437 fExtractor->FindTimeAndChargeHiGain(logain,logain,sum,dummy,dummy,dummy,sat,ped,logainabflag);
438}
439
440// --------------------------------------------------------------------------
441//
442// Compute signal mean and rms in the whole run and store it in MPedestalCam
443//
444Int_t MPedCalcFromLoGain::PostProcess()
445{
446 // Compute pedestals and rms from the whole run
447 if (fPedestalUpdate)
448 return kTRUE;
449
450 *fLog << flush << inf << "Calculating Pedestals..." << flush;
451
452 const Int_t npix = fGeom->GetNumPixels();
453 for (Int_t idx=0; idx<npix; idx++)
454 {
455 const ULong_t n = fNumEventsUsed[idx];
456 if (n>1)
457 {
458 CalcPixResults(n, idx);
459 fTotalCounter[idx]++;
460 }
461 }
462
463 CalcAreaResult();
464 CalcSectorResult();
465
466 fPedestalsOut->SetReadyToSave();
467
468 return MExtractPedestal::PostProcess();
469}
470
471
472// --------------------------------------------------------------------------
473//
474// The following resources are available:
475// FirstCheckWindowSlice: 0
476// LastCheckWindowSlice: 29
477// MaxSignalVar: 40
478//
479Int_t MPedCalcFromLoGain::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
480{
481 Bool_t rc=kFALSE;
482
483 // Find resources for CheckWindow
484 Int_t fs = fCheckWinFirst;
485 Int_t ls = fCheckWinLast;
486 if (IsEnvDefined(env, prefix, "CheckWinFirst", print))
487 {
488 fs = GetEnvValue(env, prefix, "CheckWinFirst", fs);
489 rc = kTRUE;
490 }
491 if (IsEnvDefined(env, prefix, "CheckWinLast", print))
492 {
493 ls = GetEnvValue(env, prefix, "CheckWinLast", ls);
494 rc = kTRUE;
495 }
496
497 SetCheckRange(fs,ls);
498
499 // find resource for maximum signal variation
500 if (IsEnvDefined(env, prefix, "MaxSignalVar", print))
501 {
502 SetMaxSignalVar(GetEnvValue(env, prefix, "MaxSignalVar", fMaxSignalVar));
503 rc = kTRUE;
504 }
505
506 return MExtractPedestal::ReadEnv(env,prefix,print) ? kTRUE : rc;
507}
508
509void MPedCalcFromLoGain::Print(Option_t *o) const
510{
511
512 MExtractPedestal::Print(o);
513
514 const Int_t last = fExtractor
515 ? fExtractWinFirst + fExtractor->GetWindowSizeHiGain() -1
516 : fExtractWinLast;
517
518 *fLog << "ExtractWindow from slice " << fExtractWinFirst << " to " << last << " incl." << endl;
519 *fLog << "Max.allowed signal variation: " << fMaxSignalVar << endl;
520 *fLog << "CheckWindow from slice " << fCheckWinFirst << " to " << fCheckWinLast << " incl." << endl;
521}
Note: See TracBrowser for help on using the repository browser.