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

Last change on this file since 5504 was 5504, 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 = 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 {
245 SetExtractWindow(fExtractor->GetHiGainFirst(),(Int_t)fExtractor->GetNumHiGainSamples());
246 lastavailable = losamples-1;
247 }
248
249 // If the size is not yet set, set the size
250 if (fSumx.GetSize()==0)
251 {
252 const Int_t npixels = fPedestalsOut->GetSize();
253 fNumEventsUsed.Set(npixels);
254 fTotalCounter.Set(npixels);
255 }
256
257 if (fExtractWinLast > lastavailable) //changed to override check
258 {
259 const UShort_t diff = fExtractWinLast - lastavailable;
260 *fLog << warn << GetDescriptor();
261 *fLog << " - WARNING: Selected Extract Window ranges out of range...adjusting to last available slice ";
262 *fLog << lastavailable << endl;
263
264 fExtractWinLast -= diff;
265 fExtractWinSize -= diff;
266 }
267
268 lastavailable = hisamples + losamples -1;
269
270 if (fCheckWinLast > lastavailable) //changed to override check
271 {
272 *fLog << warn << GetDescriptor();
273 *fLog << " - WARNING: Last Check Window slice out of range...adjusting to last available slice ";
274 *fLog << lastavailable << endl;
275
276 fCheckWinLast = lastavailable;
277 }
278
279 MExtractPedestal::ReInit(pList);
280
281 return kTRUE;
282}
283
284// ---------------------------------------------------------------------------------
285//
286// Returns the pointer to slice "slice".
287//
288UShort_t MPedCalcFromLoGain::GetSlice(MRawEvtPixelIter *pixel, UInt_t slice)
289{
290 const UShort_t nh = (Int_t)fRunHeader->GetNumSamplesHiGain();
291
292 Byte_t *ptr;
293
294 if(slice<nh)
295 ptr = pixel->GetHiGainSamples() + slice;
296 else
297 ptr = pixel->GetLoGainSamples() + slice - nh;
298
299 return *ptr;
300}
301
302// --------------------------------------------------------------------------
303//
304// Fill the MPedestalCam container with the signal mean and rms for the event.
305// Store the measured signal in arrays fSumx and fSumx2 so that we can
306// calculate the overall mean and rms in the PostProcess()
307//
308Int_t MPedCalcFromLoGain::Process()
309{
310
311 MRawEvtPixelIter pixel(fRawEvt);
312
313 while (pixel.Next())
314 {
315 const UInt_t idx = pixel.GetPixelId();
316 const UInt_t aidx = (*fGeom)[idx].GetAidx();
317 const UInt_t sector = (*fGeom)[idx].GetSector();
318
319 UShort_t max = 0;
320 UShort_t min = (UShort_t)-1;
321
322 // Find the maximum and minimum signal per slice in the high gain window
323 for (Int_t slice=fCheckWinFirst; slice<=fCheckWinLast; slice++)
324 {
325 const UShort_t svalue = GetSlice(&pixel,slice);
326 if (svalue > max)
327 max = svalue;
328 if (svalue < min)
329 min = svalue;
330 }
331
332 // If the maximum in the high gain window is smaller than
333 if (max-min>=fMaxSignalVar || max>=250)
334 continue;
335
336 Float_t sum = 0.;
337 UInt_t sumi = 0;
338
339 //extract pedestal
340 if (fExtractor)
341 CalcExtractor( &pixel, sum, (*fPedestalsIn)[idx] );
342 else
343 {
344 for(Int_t slice=fExtractWinFirst; slice<=fExtractWinLast; slice++)
345 sumi += GetSlice(&pixel,slice);
346 sum = (Float_t)sumi;
347 }
348
349 const Float_t sqrsum = sum*sum;
350
351 fSumx[idx] += sum;
352 fSumx2[idx] += sqrsum;
353 fAreaSumx[aidx] += sum;
354 fAreaSumx2[aidx] += sqrsum;
355 fSectorSumx[sector] += sum;
356 fSectorSumx2[sector] += sqrsum;
357
358 fNumEventsUsed[idx] ++;
359 fAreaFilled [aidx] ++;
360 fSectorFilled [sector]++;
361
362 if (!fExtractor)
363 {
364 //
365 // Calculate the amplitude of the 150MHz "AB" noise
366 //
367 const UShort_t abFlag = (fExtractWinFirst + pixel.HasABFlag()) & 0x1;
368
369 for (UShort_t islice=fExtractWinFirst; islice<=fExtractWinLast; islice+=2)
370 {
371 const UShort_t sliceAB0 = islice + abFlag;
372 const UShort_t sliceAB1 = islice - abFlag + 1;
373 const UShort_t ab0 = GetSlice(&pixel, sliceAB0);
374 const UShort_t ab1 = GetSlice(&pixel, sliceAB1);
375 fSumAB0[idx] += ab0;
376 fSumAB1[idx] += ab1;
377 fAreaSumAB0[aidx] += ab0;
378 fAreaSumAB1[aidx] += ab1;
379 fSectorSumAB0[aidx] += ab0;
380 fSectorSumAB1[aidx] += ab1;
381 }
382 }
383
384 if (!fPedestalUpdate || (UInt_t)fNumEventsUsed[idx]<fNumEventsDump)
385 continue;
386
387 CalcPixResults(fNumEventsDump, idx);
388 fTotalCounter[idx]++;
389
390 fNumEventsUsed[idx]=0;
391 fSumx[idx]=0;
392 fSumx2[idx]=0;
393 fSumAB0[idx]=0;
394 fSumAB1[idx]=0;
395 }
396
397 if (!(GetNumExecutions() % fNumAreasDump))
398 {
399 //
400 // Loop over the (two) area indices to get the averaged pedestal per aidx
401 //
402 for (UInt_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
403 {
404 const Int_t napix = fAreaValid.At(aidx);
405 if (napix == 0)
406 continue;
407
408 CalcAreaResults(fSectorFilled[aidx],napix,aidx);
409 }
410 }
411
412 if (!(GetNumExecutions() % fNumSectorsDump))
413 {
414 //
415 // Loop over the (two) sector indices to get the averaged pedestal per sector
416 //
417 for (UInt_t sector=0; sector<fSectorValid.GetSize(); sector++)
418 {
419 const Int_t nspix = fSectorValid.At(sector);
420 if (nspix == 0)
421 continue;
422
423 CalcSectorResults(fSectorFilled[sector],nspix,sector);
424 }
425 }
426
427 if (fPedestalUpdate)
428 fPedestalsOut->SetReadyToSave();
429
430 return kTRUE;
431}
432
433void MPedCalcFromLoGain::CalcExtractor( MRawEvtPixelIter *pixel, Float_t &sum, MPedestalPix &ped)
434{
435
436 Byte_t sat = 0;
437 Byte_t *logain = pixel->GetLoGainSamples() + fExtractWinFirst;
438 const Bool_t logainabflag = (pixel->HasABFlag() + pixel->GetNumHiGainSamples()) & 0x1;
439 Float_t dummy;
440 fExtractor->FindTimeAndChargeLoGain(logain,sum,dummy,dummy,dummy,sat,ped,logainabflag);
441}
442
443
444// --------------------------------------------------------------------------
445//
446// Compute signal mean and rms in the whole run and store it in MPedestalCam
447//
448Int_t MPedCalcFromLoGain::PostProcess()
449{
450
451 // Compute pedestals and rms from the whole run
452 if (fPedestalUpdate)
453 return kTRUE;
454
455 *fLog << flush << inf << "Calculating Pedestals..." << flush;
456
457 const Int_t npix = fGeom->GetNumPixels();
458 for (Int_t idx=0; idx<npix; idx++)
459 {
460 const ULong_t n = fNumEventsUsed[idx];
461 if (n>1)
462 {
463 CalcPixResults(n, idx);
464 fTotalCounter[idx]++;
465 }
466 }
467
468 //
469 // Loop over the (two) area indices to get the averaged pedestal per aidx
470 //
471 for (UInt_t aidx=0; aidx<fAreaFilled.GetSize(); aidx++)
472 {
473 const Int_t napix = fAreaValid.At(aidx);
474 if (napix == 0)
475 continue;
476
477 CalcAreaResults(fAreaFilled[aidx],napix,aidx);
478 }
479
480 //
481 // Loop over the (six) sector indices to get the averaged pedestal per sector
482 //
483 for (UInt_t sector=0; sector<fSectorFilled.GetSize(); sector++)
484 {
485 const Int_t nspix = fSectorValid.At(sector);
486 if (nspix == 0)
487 continue;
488
489 CalcSectorResults(fSectorFilled[sector],nspix,sector);
490 }
491
492 fPedestalsOut->SetReadyToSave();
493
494 return kTRUE;
495}
496
497
498// --------------------------------------------------------------------------
499//
500// The following resources are available:
501// FirstCheckWindowSlice: 0
502// LastCheckWindowSlice: 29
503// MaxSignalVar: 40
504//
505Int_t MPedCalcFromLoGain::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
506{
507
508 Bool_t rc=kFALSE;
509
510 // Find resources for CheckWindow
511 Int_t fs = fCheckWinFirst;
512 Int_t ls = fCheckWinLast;
513 if (IsEnvDefined(env, prefix, "CheckWinFirst", print))
514 {
515 fs = GetEnvValue(env, prefix, "CheckWinFirst", fs);
516 rc = kTRUE;
517 }
518 if (IsEnvDefined(env, prefix, "CheckWinLast", print))
519 {
520 ls = GetEnvValue(env, prefix, "CheckWinLast", ls);
521 rc = kTRUE;
522 }
523
524 SetCheckRange(fs,ls);
525
526 // find resource for maximum signal variation
527 if (IsEnvDefined(env, prefix, "MaxSignalVar", print))
528 {
529 SetMaxSignalVar(GetEnvValue(env, prefix, "MaxSignalVar", fMaxSignalVar));
530 rc = kTRUE;
531 }
532
533 return MExtractPedestal::ReadEnv(env,prefix,print) ? kTRUE : rc;
534}
535
536void MPedCalcFromLoGain::Print(Option_t *o) const
537{
538
539 MExtractPedestal::Print(o);
540
541 *fLog << "CheckWindow from slice " << fCheckWinFirst << " to " << fCheckWinLast << " incl." << endl;
542 *fLog << "Max. allowed signal variation: " << fMaxSignalVar << endl;
543 *fLog << endl;
544}
Note: See TracBrowser for help on using the repository browser.