source: trunk/MagicSoft/Mars/mpedestal/MPedCalcFromLoGain.cc@ 5520

Last change on this file since 5520 was 5507, checked in by gaug, 20 years ago
*** empty log message ***
File size: 16.3 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 = 15
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#include "MExtractPedestal.h"
137
138#include "MExtractTimeAndCharge.h"
139
140#include "MParList.h"
141
142#include "MLog.h"
143#include "MLogManip.h"
144
145#include "MRawRunHeader.h"
146#include "MRawEvtPixelIter.h"
147#include "MRawEvtData.h"
148
149#include "MPedestalPix.h"
150#include "MPedestalCam.h"
151
152#include "MGeomPix.h"
153#include "MGeomCam.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 = 15;
162const UShort_t MPedCalcFromLoGain::fgExtractWinSize = 6;
163const UShort_t MPedCalcFromLoGain::fgMaxSignalVar = 40;
164// --------------------------------------------------------------------------
165//
166// Default constructor:
167//
168// Sets:
169// - all pointers to NULL
170//
171// Calls:
172// - AddToBranchList("fHiGainPixId");
173// - AddToBranchList("fHiGainFadcSamples");
174// - SetCheckRange(fgCheckWinFirst, fgCheckWinLast, fgExtractWinFirst, fgExtractWinSize)
175// - Clear()
176//
177MPedCalcFromLoGain::MPedCalcFromLoGain(const char *name, const char *title)
178{
179
180 fName = name ? name : "MPedCalcFromLoGain";
181 fTitle = title ? title : "Task to calculate pedestals from lo-gains";
182
183 SetCheckRange(fgCheckWinFirst, fgCheckWinLast);
184 SetExtractWindow(fgExtractWinFirst, fgExtractWinSize);
185
186 SetMaxSignalVar(fgMaxSignalVar);
187
188 Clear();
189}
190
191void MPedCalcFromLoGain::ResetArrays()
192{
193
194 MExtractPedestal::ResetArrays();
195
196 fNumEventsUsed.Reset();
197 fTotalCounter.Reset();
198}
199
200// --------------------------------------------------------------------------
201//
202// SetCheckRange:
203//
204// Exits, if the first argument is smaller than 0
205// Exits, if the the last argument is smaller than the first
206//
207Bool_t MPedCalcFromLoGain::SetCheckRange(UShort_t chfirst, UShort_t chlast)
208{
209
210 Bool_t rc = kTRUE;
211
212 if (chlast<=chfirst)
213 {
214 *fLog << warn << GetDescriptor();
215 *fLog << " - WARNING: Last slice in SetCheckRange smaller than first slice... set to first+2" << endl;
216 chlast = chfirst+1;
217 rc = kFALSE;
218 }
219
220 fCheckWinFirst = chfirst;
221 fCheckWinLast = chlast;
222
223 return rc;
224}
225
226// ---------------------------------------------------------------------------------
227//
228// Checks:
229// - if the number of available slices
230// (fRunHeader->GetNumSamplesHiGain()+(Int_t)fRunHeader->GetNumSamplesLoGain()-1)
231// is smaller than the number of used slices
232// (fExtractWinSize+ fExtractWinFirst-1) or
233// fCheckWinLast
234//
235Bool_t MPedCalcFromLoGain::ReInit(MParList *pList)
236{
237
238 const UShort_t hisamples = fRunHeader->GetNumSamplesHiGain();
239 const UShort_t losamples = fRunHeader->GetNumSamplesLoGain();
240
241 UShort_t lastavailable = hisamples+losamples-1;
242
243 if (fExtractor)
244 lastavailable = losamples-1;
245
246 // If the size is not yet set, set the size
247 if (fSumx.GetSize()==0)
248 {
249 const Int_t npixels = fPedestalsOut->GetSize();
250 fNumEventsUsed.Set(npixels);
251 fTotalCounter.Set(npixels);
252 }
253
254 if (fExtractWinLast > lastavailable) //changed to override check
255 {
256 const UShort_t diff = fExtractWinLast - lastavailable;
257 *fLog << warn << GetDescriptor();
258 *fLog << " - WARNING: Selected Extract Window ranges out of range...adjusting to last available slice ";
259 *fLog << lastavailable << endl;
260
261 fExtractWinLast -= diff;
262 fExtractWinSize -= diff;
263 }
264
265 lastavailable = hisamples + losamples -1;
266
267 if (fCheckWinLast > lastavailable) //changed to override check
268 {
269 *fLog << warn << GetDescriptor();
270 *fLog << " - WARNING: Last Check Window slice out of range...adjusting to last available slice ";
271 *fLog << lastavailable << endl;
272
273 fCheckWinLast = lastavailable;
274 }
275
276 MExtractPedestal::ReInit(pList);
277
278 return kTRUE;
279}
280
281// ---------------------------------------------------------------------------------
282//
283// Returns the pointer to slice "slice".
284//
285UShort_t MPedCalcFromLoGain::GetSlice(MRawEvtPixelIter *pixel, UInt_t slice)
286{
287 const UShort_t nh = (Int_t)fRunHeader->GetNumSamplesHiGain();
288
289 Byte_t *ptr;
290
291 if(slice<nh)
292 ptr = pixel->GetHiGainSamples() + slice;
293 else
294 ptr = pixel->GetLoGainSamples() + slice - nh;
295
296 return *ptr;
297}
298
299// --------------------------------------------------------------------------
300//
301// Fill the MPedestalCam container with the signal mean and rms for the event.
302// Store the measured signal in arrays fSumx and fSumx2 so that we can
303// calculate the overall mean and rms in the PostProcess()
304//
305Int_t MPedCalcFromLoGain::Process()
306{
307
308 MRawEvtPixelIter pixel(fRawEvt);
309
310 while (pixel.Next())
311 {
312 const UInt_t idx = pixel.GetPixelId();
313 const UInt_t aidx = (*fGeom)[idx].GetAidx();
314 const UInt_t sector = (*fGeom)[idx].GetSector();
315
316 UShort_t max = 0;
317 UShort_t min = (UShort_t)-1;
318
319 // Find the maximum and minimum signal per slice in the high gain window
320 for (Int_t slice=fCheckWinFirst; slice<=fCheckWinLast; slice++)
321 {
322 const UShort_t svalue = GetSlice(&pixel,slice);
323 if (svalue > max)
324 max = svalue;
325 if (svalue < min)
326 min = svalue;
327 }
328
329 // If the maximum in the high gain window is smaller than
330 if (max-min>=fMaxSignalVar || max>=250)
331 continue;
332
333 Float_t sum = 0.;
334 UInt_t sumi = 0;
335
336 //extract pedestal
337 if (fExtractor)
338 CalcExtractor( &pixel, sum, (*fPedestalsIn)[idx] );
339 else
340 {
341 for(Int_t slice=fExtractWinFirst; slice<=fExtractWinLast; slice++)
342 sumi += GetSlice(&pixel,slice);
343 sum = (Float_t)sumi;
344 }
345
346 const Float_t sqrsum = sum*sum;
347
348 fSumx[idx] += sum;
349 fSumx2[idx] += sqrsum;
350 fAreaSumx[aidx] += sum;
351 fAreaSumx2[aidx] += sqrsum;
352 fSectorSumx[sector] += sum;
353 fSectorSumx2[sector] += sqrsum;
354
355 fNumEventsUsed[idx] ++;
356 fAreaFilled [aidx] ++;
357 fSectorFilled [sector]++;
358
359 if (!fExtractor)
360 {
361 //
362 // Calculate the amplitude of the 150MHz "AB" noise
363 //
364 const UShort_t abFlag = (fExtractWinFirst + pixel.HasABFlag()) & 0x1;
365
366 for (UShort_t islice=fExtractWinFirst; islice<=fExtractWinLast; islice+=2)
367 {
368 const UShort_t sliceAB0 = islice + abFlag;
369 const UShort_t sliceAB1 = islice - abFlag + 1;
370 const UShort_t ab0 = GetSlice(&pixel, sliceAB0);
371 const UShort_t ab1 = GetSlice(&pixel, sliceAB1);
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 {
396 //
397 // Loop over the (two) area indices to get the averaged pedestal per aidx
398 //
399 for (UInt_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
400 {
401 const Int_t napix = fAreaValid.At(aidx);
402 if (napix == 0)
403 continue;
404
405 CalcAreaResults(fSectorFilled[aidx],napix,aidx);
406 }
407 }
408
409 if (!(GetNumExecutions() % fNumSectorsDump))
410 {
411 //
412 // Loop over the (two) sector indices to get the averaged pedestal per sector
413 //
414 for (UInt_t sector=0; sector<fSectorValid.GetSize(); sector++)
415 {
416 const Int_t nspix = fSectorValid.At(sector);
417 if (nspix == 0)
418 continue;
419
420 CalcSectorResults(fSectorFilled[sector],nspix,sector);
421 }
422 }
423
424 if (fPedestalUpdate)
425 fPedestalsOut->SetReadyToSave();
426
427 return kTRUE;
428}
429
430void MPedCalcFromLoGain::CalcExtractor( MRawEvtPixelIter *pixel, Float_t &sum, MPedestalPix &ped)
431{
432
433 Byte_t sat = 0;
434 Byte_t *logain = pixel->GetLoGainSamples() + fExtractWinFirst;
435 const Bool_t logainabflag = (pixel->HasABFlag() + pixel->GetNumHiGainSamples()) & 0x1;
436 Float_t dummy;
437 fExtractor->FindTimeAndChargeHiGain(logain,logain,sum,dummy,dummy,dummy,sat,ped,logainabflag);
438}
439
440
441// --------------------------------------------------------------------------
442//
443// Compute signal mean and rms in the whole run and store it in MPedestalCam
444//
445Int_t MPedCalcFromLoGain::PostProcess()
446{
447
448 // Compute pedestals and rms from the whole run
449 if (fPedestalUpdate)
450 return kTRUE;
451
452 *fLog << flush << inf << "Calculating Pedestals..." << flush;
453
454 const Int_t npix = fGeom->GetNumPixels();
455 for (Int_t idx=0; idx<npix; idx++)
456 {
457 const ULong_t n = fNumEventsUsed[idx];
458 if (n>1)
459 {
460 CalcPixResults(n, idx);
461 fTotalCounter[idx]++;
462 }
463 }
464
465 //
466 // Loop over the (two) area indices to get the averaged pedestal per aidx
467 //
468 for (UInt_t aidx=0; aidx<fAreaFilled.GetSize(); aidx++)
469 {
470 const Int_t napix = fAreaValid.At(aidx);
471 if (napix == 0)
472 continue;
473
474 CalcAreaResults(fAreaFilled[aidx],napix,aidx);
475 }
476
477 //
478 // Loop over the (six) sector indices to get the averaged pedestal per sector
479 //
480 for (UInt_t sector=0; sector<fSectorFilled.GetSize(); sector++)
481 {
482 const Int_t nspix = fSectorValid.At(sector);
483 if (nspix == 0)
484 continue;
485
486 CalcSectorResults(fSectorFilled[sector],nspix,sector);
487 }
488
489 fPedestalsOut->SetReadyToSave();
490
491 return kTRUE;
492}
493
494
495// --------------------------------------------------------------------------
496//
497// The following resources are available:
498// FirstCheckWindowSlice: 0
499// LastCheckWindowSlice: 29
500// MaxSignalVar: 40
501//
502Int_t MPedCalcFromLoGain::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
503{
504
505 Bool_t rc=kFALSE;
506
507 // Find resources for CheckWindow
508 Int_t fs = fCheckWinFirst;
509 Int_t ls = fCheckWinLast;
510 if (IsEnvDefined(env, prefix, "CheckWinFirst", print))
511 {
512 fs = GetEnvValue(env, prefix, "CheckWinFirst", fs);
513 rc = kTRUE;
514 }
515 if (IsEnvDefined(env, prefix, "CheckWinLast", print))
516 {
517 ls = GetEnvValue(env, prefix, "CheckWinLast", ls);
518 rc = kTRUE;
519 }
520
521 SetCheckRange(fs,ls);
522
523 // find resource for maximum signal variation
524 if (IsEnvDefined(env, prefix, "MaxSignalVar", print))
525 {
526 SetMaxSignalVar(GetEnvValue(env, prefix, "MaxSignalVar", fMaxSignalVar));
527 rc = kTRUE;
528 }
529
530 return MExtractPedestal::ReadEnv(env,prefix,print) ? kTRUE : rc;
531}
532
533void MPedCalcFromLoGain::Print(Option_t *o) const
534{
535
536 MExtractPedestal::Print(o);
537
538 *fLog << "CheckWindow from slice " << fCheckWinFirst << " to " << fCheckWinLast << " incl." << endl;
539 *fLog << "Max. allowed signal variation: " << fMaxSignalVar << endl;
540 *fLog << endl;
541}
Note: See TracBrowser for help on using the repository browser.