source: trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc@ 6292

Last change on this file since 6292 was 6280, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 13.8 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!
24! Copyright: MAGIC Software Development, 2000-2004
25!
26!
27\* ======================================================================== */
28/////////////////////////////////////////////////////////////////////////////
29//
30// MPedCalcPedRun
31//
32// This task takes a pedestal run file and fills MPedestalCam during
33// the Process() with the pedestal and rms computed in an event basis.
34// In the PostProcess() MPedestalCam is finally filled with the pedestal
35// mean and rms computed in a run basis.
36// More than one run (file) can be merged
37//
38// MPedCalcPedRun applies the following formula (1):
39//
40// Pedestal per slice = sum(x_i) / n / slices
41// PedRMS per slice = Sqrt( ( sum(x_i^2) - sum(x_i)^2/n ) / n-1 / slices )
42//
43// where x_i is the sum of "slices" FADC slices and sum means the sum over all
44// events. "n" is the number of events, "slices" is the number of summed FADC samples.
45//
46// Note that the slice-to-slice fluctuations are not Gaussian, but Poissonian, thus
47// asymmetric and they are correlated.
48//
49// It is important to know that the Pedestal per slice and PedRMS per slice depend
50// on the number of used FADC slices, as seen in the following plots:
51//
52//Begin_Html
53/*
54<img src="images/PedestalStudyInner.gif">
55*/
56//End_Html
57//
58//Begin_Html
59/*
60<img src="images/PedestalStudyOuter.gif">
61*/
62//
63// The plots show the inner and outer pixels, respectivly and have the following meaning:
64//
65// 1) The calculated mean pedestal per slice (from MPedCalcPedRun)
66// 2) The fitted mean pedestal per slice (from MHPedestalCam)
67// 3) The calculated pedestal RMS per slice (from MPedCalcPedRun)
68// 4) The fitted sigma of the pedestal distribution per slice
69// (from MHPedestalCam)
70// 5) The relative difference between calculation and histogram fit
71// for the mean
72// 6) The relative difference between calculation and histogram fit
73// for the sigma or RMS, respectively.
74//
75// The calculated means do not change significantly except for the case of 2 slices,
76// however the RMS changes from 5.7 per slice in the case of 2 extracted slices
77// to 8.3 per slice in the case of 26 extracted slices. This change is very significant.
78//
79// The plots have been produced on run 20123. You can reproduce them using
80// the macro pedestalstudies.C
81//
82// Usage of this class:
83// ====================
84//
85// Call: SetRange(higainfirst, higainlast, logainfirst, logainlast)
86// to modify the ranges in which the window is allowed to move.
87// Defaults are:
88//
89// fHiGainFirst = fgHiGainFirst = 0
90// fHiGainLast = fgHiGainLast = 29
91// fLoGainFirst = fgLoGainFirst = 0
92// fLoGainLast = fgLoGainLast = 14
93//
94// Call: SetWindowSize(windowhigain, windowlogain)
95// to modify the sliding window widths. Windows have to be an even number.
96// In case of odd numbers, the window will be modified.
97//
98// Defaults are:
99//
100// fHiGainWindowSize = fgHiGainWindowSize = 14
101// fLoGainWindowSize = fgLoGainWindowSize = 0
102//
103//
104// ToDo:
105// - Take a setup file (ReadEnv-implementation) as input
106//
107//
108// Input Containers:
109// MRawEvtData
110// MRawRunHeader
111// MRawEvtHeader
112// MGeomCam
113//
114// Output Containers:
115// MPedestalCam
116//
117// See also: MPedestalCam, MPedestalPix, MHPedestalCam, MExtractor
118//
119/////////////////////////////////////////////////////////////////////////////
120#include "MPedCalcPedRun.h"
121
122#include "MParList.h"
123
124#include "MLog.h"
125#include "MLogManip.h"
126
127#include "MRawRunHeader.h"
128#include "MRawEvtHeader.h"
129#include "MRawEvtPixelIter.h"
130#include "MRawEvtData.h"
131
132#include "MPedestalPix.h"
133#include "MPedestalCam.h"
134
135#include "MGeomPix.h"
136#include "MGeomCam.h"
137
138#include "MExtractPedestal.h"
139#include "MExtractTimeAndCharge.h"
140
141#include "MTriggerPattern.h"
142
143ClassImp(MPedCalcPedRun);
144
145using namespace std;
146
147const UShort_t MPedCalcPedRun::fgExtractWinFirst = 0;
148const UShort_t MPedCalcPedRun::fgExtractWinSize = 6;
149const UInt_t MPedCalcPedRun::gkFirstRunWithFinalBits = 45605;
150// --------------------------------------------------------------------------
151//
152// Default constructor:
153//
154// Sets:
155// - all pointers to NULL
156// - fWindowSizeHiGain to fgHiGainWindowSize
157// - fWindowSizeLoGain to fgLoGainWindowSize
158//
159// Calls:
160// - AddToBranchList("fHiGainPixId");
161// - AddToBranchList("fHiGainFadcSamples");
162// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
163//
164MPedCalcPedRun::MPedCalcPedRun(const char *name, const char *title)
165 : fOverlap(0), fIsFirstPedRun(kFALSE), fUsedEvents(0), fTrigPattern(NULL)
166{
167 fName = name ? name : "MPedCalcPedRun";
168 fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data";
169
170 SetExtractWindow(fgExtractWinFirst, fgExtractWinSize);
171
172 // SetNumEventsDump(1001);
173 // SetNumAreasDump(1001);
174 // SetNumSectorsDump(1001);
175}
176
177// --------------------------------------------------------------------------
178//
179// In case that the variables fExtractLast is greater than the number of
180// High-Gain FADC samples obtained from the run header, the flag
181// fOverlap is set to the difference and fExtractLast is set back by the
182// same number of slices.
183//
184void MPedCalcPedRun::CheckExtractionWindow()
185{
186 const UShort_t lastavailable = fRunHeader->GetNumSamplesHiGain()-1;
187
188 if (fExtractWinLast <= lastavailable)
189 return;
190
191 const UShort_t diff = fExtractWinLast - lastavailable;
192
193 *fLog << warn << endl;
194 *fLog << "Selected ExtractWindow [" << fExtractWinFirst << "," << fExtractWinLast;
195 *fLog << "] ranges out of range [0," << lastavailable << "]." << endl;
196 *fLog << "Using " << diff << " samples from the 'Lo-Gain' slices for the pedestal extraction" << endl;
197
198 fExtractWinLast -= diff;
199 fOverlap = diff;
200}
201
202void MPedCalcPedRun::Reset()
203{
204
205 MExtractPedestal::ResetArrays();
206 fUsedEvents = 0;
207}
208
209// --------------------------------------------------------------------------
210//
211// Set fIsFirstPedRun=kTRUE
212//
213Int_t MPedCalcPedRun::PreProcess(MParList *pList)
214{
215
216 fUsedEvents = 0;
217 fIsFirstPedRun = kTRUE;
218 fIsNotPedRun = kFALSE;
219
220 fTrigPattern = (MTriggerPattern*)pList->FindObject("MTriggerPattern");
221 if (!fTrigPattern)
222 *fLog << warn << "MTriggerPattern not found... Cannot use interlaced pedestal events." << endl;
223
224 return MExtractPedestal::PreProcess(pList);
225}
226
227// --------------------------------------------------------------------------
228//
229// The run type is checked for "kRTPedestal"
230// and the variable fIsNotPedRun is set in that case
231//
232Bool_t MPedCalcPedRun::ReInit(MParList *pList)
233{
234
235 if (!MExtractPedestal::ReInit(pList))
236 return kFALSE;
237
238 CheckExtractionWindow();
239
240 //
241 // If this is the first ped run, the next run (call to ReInit)
242 // is not the first anymore
243 //
244 if (fRunHeader->GetRunType()==MRawRunHeader::kRTPedestal)
245 {
246 fIsFirstPedRun = kFALSE;
247 fIsNotPedRun = kFALSE;
248 return kTRUE;
249 }
250
251 if (fRunHeader->GetRunType()==MRawRunHeader::kRTCalibration)
252 {
253 TString proj(fRunHeader->GetProjectName());
254 proj.ToLower();
255
256 // test if a continuous light run has been declared as calibration...
257 if (proj.Contains("cl"))
258 {
259 fIsFirstPedRun = kFALSE;
260 fIsNotPedRun = kFALSE;
261 return kTRUE;
262 }
263 }
264
265
266
267 fIsNotPedRun = kTRUE;
268 //
269 // If this is the first call to ReInit (before reading the first file)
270 // nothing should be done. It occurs for the case that the first treated
271 // file is not of pedestal type.
272 //
273 if (fIsFirstPedRun)
274 return kTRUE;
275
276 //
277 // In the other case, produce the MPedestalCam to be use the subsequent (non-pedestal) events
278 //
279 *fLog << inf << "Finalizing pedestal calculations..." << flush;
280
281 if (!Finalize())
282 return kFALSE;
283
284 Reset();
285
286 return kTRUE;
287}
288
289// --------------------------------------------------------------------------
290//
291// Return kTRUE (without doing anything) in case that the run type is not
292// equal to 1 (pedestal run)
293//
294// Fill the MPedestalCam container with the signal mean and rms for the event.
295// Store the measured signal in arrays fSumx and fSumx2 so that we can
296// calculate the overall mean and rms in the PostProcess()
297//
298Int_t MPedCalcPedRun::Calc()
299{
300
301 if (fIsNotPedRun && !IsPedBitSet())
302 return kTRUE;
303
304 // if (fUsedEvents == fNumEventsDump)
305 // {
306 // *fLog << endl;
307 // *fLog << inf << GetDescriptor() << " : Finalize pedestal calculations..." << flush;
308 // Finalize();
309 // Reset();
310 // }
311
312 fUsedEvents++;
313
314 MRawEvtPixelIter pixel(fRawEvt);
315
316 while (pixel.Next())
317 {
318 const UInt_t idx = pixel.GetPixelId();
319 const UInt_t aidx = (*fGeom)[idx].GetAidx();
320 const UInt_t sector = (*fGeom)[idx].GetSector();
321
322 Float_t sum = 0.;
323 UInt_t ab0 = 0;
324 UInt_t ab1 = 0;
325
326 if (fExtractor)
327 CalcExtractor(pixel, sum, (*fPedestalsIn)[idx]);
328 else
329 CalcSums(pixel, sum, ab0, ab1);
330
331 fSumx[idx] += sum;
332 fAreaSumx[aidx] += sum;
333 fSectorSumx[sector] += sum;
334
335 const Float_t sqrsum = sum*sum;
336 fSumx2[idx] += sqrsum;
337 fAreaSumx2[aidx] += sqrsum;
338 fSectorSumx2[sector] += sqrsum;
339
340 fSumAB0[idx] += ab0;
341 fSumAB1[idx] += ab1;
342
343 fAreaSumAB0[aidx] += ab0;
344 fAreaSumAB1[aidx] += ab1;
345
346 fSectorSumAB0[aidx] += ab0;
347 fSectorSumAB1[aidx] += ab1;
348 }
349
350 fPedestalsOut->SetReadyToSave();
351
352 return kTRUE;
353}
354
355
356void MPedCalcPedRun::CalcExtractor(const MRawEvtPixelIter &pixel, Float_t &sum, MPedestalPix &ped)
357{
358 const Bool_t abflag = pixel.HasABFlag();
359
360 Byte_t *first = pixel.GetHiGainSamples() + fExtractWinFirst;
361 Byte_t *logain = pixel.GetLoGainSamples();
362
363 Byte_t sat = 0;
364
365 Float_t dummy;
366 fExtractor->FindTimeAndChargeHiGain(first,logain,sum,dummy,dummy,dummy,sat,ped,abflag);
367}
368
369void MPedCalcPedRun::CalcSums(const MRawEvtPixelIter &pixel, Float_t &sum, UInt_t &ab0, UInt_t &ab1)
370{
371
372 Byte_t *higain = pixel.GetHiGainSamples() + fExtractWinFirst;
373 Byte_t *ptr = higain;
374 Byte_t *end = ptr + fExtractWinSize;
375
376 const Bool_t abflag = pixel.HasABFlag();
377
378 Int_t sumi = 0;
379
380 Int_t cnt = 0;
381 do
382 {
383 sumi += *ptr;
384 if (!pixel.IsABFlagValid())
385 continue;
386
387 const Int_t abFlag = (fExtractWinFirst + abflag + cnt) & 0x1;
388 if (abFlag)
389 ab1 += *ptr;
390 else
391 ab0 += *ptr;
392
393 cnt++;
394 } while (++ptr != end);
395
396 if (fOverlap != 0)
397 {
398 ptr = pixel.GetLoGainSamples();
399 end = ptr + fOverlap;
400
401 do
402 {
403 sumi += *ptr;
404 if (!pixel.IsABFlagValid())
405 continue;
406
407 const Int_t abFlag = (fExtractWinFirst + abflag + cnt) & 0x1;
408 if (abFlag)
409 ab1 += *ptr;
410 else
411 ab0 += *ptr;
412
413 cnt++;
414 } while (++ptr != end);
415 }
416
417 sum = (Float_t)sumi;
418}
419
420
421// --------------------------------------------------------------------------
422//
423// Compute signal mean and rms in the whole run and store it in MPedestalCam
424//
425Int_t MPedCalcPedRun::Finalize()
426{
427
428 if (fUsedEvents == 0)
429 return kTRUE;
430
431 MRawEvtPixelIter pixel(fRawEvt);
432 while (pixel.Next())
433 CalcPixResults(fUsedEvents, pixel.GetPixelId());
434
435 //
436 // Loop over the (two) area indices to get the averaged pedestal per aidx
437 //
438 for (UInt_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
439 if (fAreaValid[aidx]>0)
440 CalcAreaResults(fUsedEvents, fAreaValid[aidx], aidx);
441
442 //
443 // Loop over the (six) sector indices to get the averaged pedestal per sector
444 //
445 for (UInt_t sector=0; sector<fSectorValid.GetSize(); sector++)
446 if (fSectorValid[sector]>0)
447 CalcSectorResults(fUsedEvents, fSectorValid[sector], sector);
448
449 fPedestalsOut->SetTotalEntries(fUsedEvents*fExtractWinSize);
450 fPedestalsOut->SetReadyToSave();
451
452 return kTRUE;
453}
454
455Int_t MPedCalcPedRun::PostProcess()
456{
457 if (!Finalize())
458 return kFALSE;
459
460 return MExtractPedestal::PostProcess();
461}
462
463//-------------------------------------------------------------
464//
465// Return if the pedestal bit was set from the calibration trigger box.
466// The last but one bit is used for the "pedestal-bit".
467//
468// This bit is set since run gkFirstRunWithFinalBits
469//
470Bool_t MPedCalcPedRun::IsPedBitSet()
471{
472
473 if (!fTrigPattern)
474 return kTRUE;
475
476 if (fRunHeader->GetRunNumber() < gkFirstRunWithFinalBits)
477 return kTRUE;
478
479 return (fTrigPattern->GetPrescaled() & MTriggerPattern::kPedestal) ? kTRUE : kFALSE;
480}
481
482void MPedCalcPedRun::Print(Option_t *o) const
483{
484
485 MExtractPedestal::Print(o);
486
487 const Int_t last = fExtractor
488 ? fExtractWinFirst + fExtractor->GetWindowSizeHiGain() -1
489 : fExtractWinLast;
490
491 *fLog << "ExtractWindow from slice " << fExtractWinFirst << " to " << last << " incl." << endl;
492 *fLog << "Num overlap lo-gain slices: " << fOverlap << endl;
493 *fLog << "First pedrun out of sequence: " << (fIsFirstPedRun?"yes":"no") << endl;
494 *fLog << "Number of used events so far: " << fUsedEvents << endl;
495}
Note: See TracBrowser for help on using the repository browser.