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

Last change on this file since 4013 was 4013, checked in by gaug, 21 years ago
*** empty log message ***
File size: 20.1 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 = 14
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// Input Containers:
104// MRawEvtData
105// MRawRunHeader
106// MGeomCam
107//
108// Output Containers:
109// MPedestalCam
110//
111// See also: MPedestalCam, MPedestalPix, MHPedestalCam, MExtractor
112//
113/////////////////////////////////////////////////////////////////////////////
114#include "MPedCalcPedRun.h"
115#include "MExtractor.h"
116
117#include "MParList.h"
118
119#include "MLog.h"
120#include "MLogManip.h"
121
122#include "MRawRunHeader.h"
123#include "MRawEvtPixelIter.h"
124#include "MRawEvtData.h"
125
126#include "MPedestalPix.h"
127#include "MPedestalCam.h"
128
129#include "MGeomPix.h"
130#include "MGeomCam.h"
131
132#include "MBadPixelsPix.h"
133#include "MBadPixelsCam.h"
134
135#include "MGeomCamMagic.h"
136
137ClassImp(MPedCalcPedRun);
138
139using namespace std;
140
141const Byte_t MPedCalcPedRun::fgHiGainFirst = 0;
142const Byte_t MPedCalcPedRun::fgHiGainLast = 14;
143const Byte_t MPedCalcPedRun::fgLoGainFirst = 0;
144const Byte_t MPedCalcPedRun::fgLoGainLast = 14;
145const Byte_t MPedCalcPedRun::fgHiGainWindowSize = 14;
146const Byte_t MPedCalcPedRun::fgLoGainWindowSize = 0;
147// --------------------------------------------------------------------------
148//
149// Default constructor:
150//
151// Sets:
152// - all pointers to NULL
153// - fWindowSizeHiGain to fgHiGainWindowSize
154// - fWindowSizeLoGain to fgLoGainWindowSize
155//
156// Calls:
157// - AddToBranchList("fHiGainPixId");
158// - AddToBranchList("fHiGainFadcSamples");
159// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
160// - Clear()
161//
162MPedCalcPedRun::MPedCalcPedRun(const char *name, const char *title)
163 : fWindowSizeHiGain(fgHiGainWindowSize),
164 fWindowSizeLoGain(fgLoGainWindowSize),
165 fGeom(NULL),fBad(NULL)
166{
167 fName = name ? name : "MPedCalcPedRun";
168 fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data";
169
170 AddToBranchList("fHiGainPixId");
171 AddToBranchList("fHiGainFadcSamples");
172
173 SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
174 Clear();
175}
176
177// --------------------------------------------------------------------------
178//
179// Sets:
180// - fNumSamplesTot to 0
181// - fRawEvt to NULL
182// - fRunHeader to NULL
183// - fPedestals to NULL
184//
185void MPedCalcPedRun::Clear(const Option_t *o)
186{
187
188 fNumSamplesTot = 0;
189
190 fRawEvt = NULL;
191 fRunHeader = NULL;
192 fPedestals = NULL;
193}
194
195// --------------------------------------------------------------------------
196//
197// SetRange:
198//
199// Calls:
200// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
201// - SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
202//
203void MPedCalcPedRun::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
204{
205
206 MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
207
208 //
209 // Redo the checks if the window is still inside the ranges
210 //
211 SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
212
213}
214
215
216// --------------------------------------------------------------------------
217//
218// Checks:
219// - if a window is odd, subtract one
220// - if a window is bigger than the one defined by the ranges, set it to the available range
221// - if a window is smaller than 2, set it to 2
222//
223// Sets:
224// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain
225// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain
226// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
227// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)
228//
229void MPedCalcPedRun::SetWindowSize(Byte_t windowh, Byte_t windowl)
230{
231
232 fWindowSizeHiGain = windowh & ~1;
233 fWindowSizeLoGain = windowl & ~1;
234
235 if (fWindowSizeHiGain != windowh)
236 *fLog << warn << GetDescriptor() << ": Hi Gain window size has to be even, set to: "
237 << int(fWindowSizeHiGain) << " samples " << endl;
238
239 if (fWindowSizeLoGain != windowl)
240 *fLog << warn << GetDescriptor() << ": Lo Gain window size has to be even, set to: "
241 << int(fWindowSizeLoGain) << " samples " << endl;
242
243 const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1;
244 const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1;
245
246 if (fWindowSizeHiGain > availhirange)
247 {
248 *fLog << warn << GetDescriptor()
249 << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",(int)fWindowSizeHiGain,
250 " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl;
251 *fLog << warn << GetDescriptor()
252 << ": Will set window size to: " << (int)availhirange << endl;
253 fWindowSizeHiGain = availhirange;
254 }
255
256 if (fWindowSizeLoGain > availlorange)
257 {
258 *fLog << warn << GetDescriptor()
259 << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",(int)fWindowSizeLoGain,
260 " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl;
261 *fLog << warn << GetDescriptor()
262 << ": Will set window size to: " << (int)availlorange << endl;
263 fWindowSizeLoGain = availlorange;
264 }
265
266
267 fNumHiGainSamples = (Float_t)fWindowSizeHiGain;
268 fNumLoGainSamples = (Float_t)fWindowSizeLoGain;
269
270 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
271 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
272
273}
274
275
276
277// --------------------------------------------------------------------------
278//
279// Look for the following input containers:
280//
281// - MRawEvtData
282// - MRawRunHeader
283// - MGeomCam
284// - MBadPixelsCam
285//
286// The following output containers are also searched and created if
287// they were not found:
288//
289// - MPedestalCam
290//
291Int_t MPedCalcPedRun::PreProcess( MParList *pList )
292{
293
294 Clear();
295
296 fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
297 if (!fRawEvt)
298 {
299 *fLog << err << "MRawEvtData not found... aborting." << endl;
300 return kFALSE;
301 }
302
303 fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
304 if (!fRunHeader)
305 {
306 *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl;
307 return kFALSE;
308 }
309
310 fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
311 if (!fGeom)
312 {
313 *fLog << err << "MGeomCam not found... aborting." << endl;
314 return kFALSE;
315 }
316
317 fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam");
318 if (!fPedestals)
319 return kFALSE;
320
321 fBad = (MBadPixelsCam*)pList->FindObject("MBadPixelsCam");
322
323 return kTRUE;
324}
325
326// --------------------------------------------------------------------------
327//
328// The ReInit searches for:
329// - MRawRunHeader::GetNumSamplesHiGain()
330// - MRawRunHeader::GetNumSamplesLoGain()
331//
332// In case that the variables fHiGainLast and fLoGainLast are smaller than
333// the even part of the number of samples obtained from the run header, a
334// warning is given an the range is set back accordingly. A call to:
335// - SetRange(fHiGainFirst, fHiGainLast-diff, fLoGainFirst, fLoGainLast) or
336// - SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff)
337// is performed in that case. The variable diff means here the difference
338// between the requested range (fHiGainLast) and the available one. Note that
339// the functions SetRange() are mostly overloaded and perform more checks,
340// modifying the ranges again, if necessary.
341//
342// A loop over the MBadPixelsCam is performed and bad pixels are set
343// to MPedestalPix::SetValid(kFALSE);
344//
345Bool_t MPedCalcPedRun::ReInit(MParList *pList)
346{
347
348 Int_t lastdesired = (Int_t)fLoGainLast;
349 Int_t lastavailable = (Int_t)fRunHeader->GetNumSamplesLoGain()-1;
350
351 if (lastdesired > lastavailable)
352 {
353 const Int_t diff = lastdesired - lastavailable;
354 *fLog << endl;
355 *fLog << warn << GetDescriptor()
356 << Form("%s%2i%s%2i%s%2i%s",": Selected Lo Gain FADC Window [",
357 (int)fLoGainFirst,",",lastdesired,
358 "] ranges out of the available limits: [0,",lastavailable,"].") << endl;
359 *fLog << GetDescriptor() << ": Will reduce the upper edge to " << (int)(fLoGainLast - diff) << endl;
360 SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff);
361 }
362
363 lastdesired = (Int_t)fHiGainLast;
364 lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
365
366 if (lastdesired > lastavailable)
367 {
368 const Int_t diff = lastdesired - lastavailable;
369 *fLog << endl;
370 *fLog << warn << GetDescriptor()
371 << Form("%s%2i%s%2i%s%2i%s",": Selected Hi Gain Range [",
372 (int)fHiGainFirst,",",lastdesired,
373 "] ranges out of the available limits: [0,",lastavailable,"].") << endl;
374 *fLog << warn << GetDescriptor()
375 << Form("%s%2i%s",": Will possibly use ",diff," samples from the Low-Gain for the High-Gain range")
376 << endl;
377 fHiGainLast -= diff;
378 fHiLoLast = diff;
379 }
380
381 lastdesired = (Int_t)fHiGainFirst+fWindowSizeHiGain-1;
382 lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
383
384 if (lastdesired > lastavailable)
385 {
386 const Int_t diff = lastdesired - lastavailable;
387 *fLog << endl;
388 *fLog << warn << GetDescriptor()
389 << Form("%s%2i%s%2i%s",": Selected Hi Gain FADC Window size ",
390 (int)fWindowSizeHiGain,
391 " ranges out of the available limits: [0,",lastavailable,"].") << endl;
392 *fLog << warn << GetDescriptor()
393 << Form("%s%2i%s",": Will use ",diff," samples from the Low-Gain for the High-Gain extraction")
394 << endl;
395
396 if ((Int_t)fWindowSizeHiGain > diff)
397 {
398 fWindowSizeHiGain -= diff;
399 fWindowSizeLoGain += diff;
400 }
401 else
402 {
403 fWindowSizeLoGain += fWindowSizeHiGain;
404 fLoGainFirst = diff-fWindowSizeHiGain;
405 fWindowSizeHiGain = 0;
406 }
407 }
408
409
410 Int_t npixels = fPedestals->GetSize();
411 Int_t areas = fPedestals->GetAverageAreas();
412 Int_t sectors = fPedestals->GetAverageSectors();
413
414 if (fSumx.GetSize()==0)
415 {
416 fSumx. Set(npixels);
417 fSumx2.Set(npixels);
418
419 fAreaSumx. Set(areas);
420 fAreaSumx2.Set(areas);
421 fAreaValid.Set(areas);
422
423 fSectorSumx. Set(sectors);
424 fSectorSumx2.Set(sectors);
425 fSectorValid.Set(sectors);
426
427 fSumx.Reset();
428 fSumx2.Reset();
429 }
430
431 if (fWindowSizeHiGain == 0 && fWindowSizeLoGain == 0)
432 {
433 *fLog << err << GetDescriptor()
434 << ": Number of extracted Slices is 0, abort ... " << endl;
435 return kFALSE;
436 }
437
438
439 *fLog << endl;
440 *fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeHiGain
441 << " HiGain FADC samples starting with slice: " << (int)fHiGainFirst << endl;
442 *fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeLoGain
443 << " LoGain FADC samples starting with slice: " << (int)fLoGainFirst << endl;
444
445
446 if (fBad)
447 {
448 const Int_t nbads = fBad->GetSize();
449 for (Int_t i=0; i<(nbads>npixels?npixels:nbads);i++)
450 if ((*fBad)[i].IsBad())
451 (*fPedestals)[i].SetValid(kFALSE);
452 }
453
454 return kTRUE;
455
456}
457// --------------------------------------------------------------------------
458//
459// Fill the MPedestalCam container with the signal mean and rms for the event.
460// Store the measured signal in arrays fSumx and fSumx2 so that we can
461// calculate the overall mean and rms in the PostProcess()
462//
463Int_t MPedCalcPedRun::Process()
464{
465
466 MRawEvtPixelIter pixel(fRawEvt);
467
468 while (pixel.Next())
469 {
470
471 const UInt_t idx = pixel.GetPixelId();
472 const UInt_t aidx = (*fGeom)[idx].GetAidx();
473 const UInt_t sector = (*fGeom)[idx].GetSector();
474
475 Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst;
476 Byte_t *end = ptr + fWindowSizeHiGain;
477
478 UInt_t sum = 0;
479 UInt_t sqr = 0;
480
481 if (fWindowSizeHiGain != 0)
482 {
483 do
484 {
485 sum += *ptr;
486 sqr += *ptr * *ptr;
487 }
488 while (++ptr != end);
489 }
490
491 if (fWindowSizeLoGain != 0)
492 {
493
494 ptr = pixel.GetLoGainSamples() + fLoGainFirst;
495 end = ptr + fWindowSizeLoGain;
496
497 do
498 {
499 sum += *ptr;
500 sqr += *ptr * *ptr;
501 }
502 while (++ptr != end);
503
504 }
505
506 const Float_t msum = (Float_t)sum;
507
508 //
509 // These three lines have been uncommented by Markus Gaug
510 // If anybody needs them, please contact me!!
511 //
512 // const Float_t higainped = msum/fNumHiGainSlices;
513 // const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSlices)/(fNumHiGainSlices-1.));
514 // (*fPedestals)[idx].Set(higainped, higainrms);
515
516 fSumx[idx] += msum;
517 fAreaSumx[aidx] += msum;
518 fSectorSumx[sector] += msum;
519 //
520 // The old version:
521 //
522 // const Float_t msqr = (Float_t)sqr;
523 // fSumx2[idx] += msqr;
524 //
525 // The new version:
526 //
527 const Float_t sqrsum = msum*msum;
528 fSumx2[idx] += sqrsum;
529 fAreaSumx2[aidx] += sqrsum;
530 fSectorSumx2[sector] += sqrsum;
531 }
532
533 fPedestals->SetReadyToSave();
534 fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain;
535
536 return kTRUE;
537}
538
539// --------------------------------------------------------------------------
540//
541// Compute signal mean and rms in the whole run and store it in MPedestalCam
542//
543Int_t MPedCalcPedRun::PostProcess()
544{
545
546 // Compute pedestals and rms from the whole run
547 const ULong_t n = fNumSamplesTot;
548 const ULong_t nevts = GetNumExecutions();
549
550 MRawEvtPixelIter pixel(fRawEvt);
551
552 while (pixel.Next())
553 {
554
555 const Int_t pixid = pixel.GetPixelId();
556 const UInt_t aidx = (*fGeom)[pixid].GetAidx();
557 const UInt_t sector = (*fGeom)[pixid].GetSector();
558
559 fAreaValid [aidx]++;
560 fSectorValid[sector]++;
561
562 const Float_t sum = fSumx.At(pixid);
563 const Float_t sum2 = fSumx2.At(pixid);
564 const Float_t higainped = sum/n;
565 //
566 // The old version:
567 //
568 // const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.));
569 //
570 // The new version:
571 //
572 // 1. Calculate the Variance of the sums:
573 Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.);
574 // 2. Scale the variance to the number of slices:
575 higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain);
576 // 3. Calculate the RMS from the Variance:
577 (*fPedestals)[pixid].Set(higainped, higainVar < 0 ? 0. : TMath::Sqrt(higainVar));
578
579 }
580
581 //
582 // Loop over the (two) area indices to get the averaged pedestal per aidx
583 //
584 for (Int_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
585 {
586
587 const Int_t napix = fAreaValid.At(aidx);
588 const Float_t sum = fAreaSumx.At(aidx);
589 const Float_t sum2 = fAreaSumx2.At(aidx);
590 const ULong_t an = napix * n;
591 const ULong_t aevts = napix * nevts;
592
593 const Float_t higainped = sum/an;
594
595 // 1. Calculate the Variance of the sums:
596 Float_t higainVar = (sum2-sum*sum/aevts)/(aevts-1.);
597 // 2. Scale the variance to the number of slices:
598 higainVar /= fWindowSizeHiGain+fWindowSizeLoGain;
599 // 3. Calculate the RMS from the Variance:
600 Float_t higainrms = TMath::Sqrt(higainVar);
601 // 4. Re-scale it with the square root of the number of involved pixels
602 // in order to be comparable to the mean of pedRMS of that area
603 higainrms *= TMath::Sqrt((Float_t)napix);
604
605 fPedestals->GetAverageArea(aidx).Set(higainped, higainrms);
606 }
607
608 //
609 // Loop over the (six) sector indices to get the averaged pedestal per sector
610 //
611 for (Int_t sector=0; sector<fSectorValid.GetSize(); sector++)
612 {
613
614 const Int_t nspix = fSectorValid.At(sector);
615 const Float_t sum = fSectorSumx.At(sector);
616 const Float_t sum2 = fSectorSumx2.At(sector);
617 const ULong_t sn = nspix * n;
618 const ULong_t sevts = nspix * nevts;
619
620 const Float_t higainped = sum/sn;
621
622 // 1. Calculate the Variance of the sums:
623 Float_t higainVar = (sum2-sum*sum/sevts)/(sevts-1.);
624 // 2. Scale the variance to the number of slices:
625 higainVar /= fWindowSizeHiGain+fWindowSizeLoGain;
626 // 3. Calculate the RMS from the Variance:
627 Float_t higainrms = TMath::Sqrt(higainVar);
628 // 4. Re-scale it with the square root of the number of involved pixels
629 // in order to be comparable to the mean of pedRMS of that sector
630 higainrms *= TMath::Sqrt((Float_t)nspix);
631
632 fPedestals->GetAverageSector(sector).Set(higainped, higainrms);
633 }
634
635 fPedestals->SetTotalEntries(fNumSamplesTot);
636 fPedestals->SetReadyToSave();
637
638 return kTRUE;
639}
Note: See TracBrowser for help on using the repository browser.