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

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