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

Last change on this file since 4147 was 4147, checked in by gaug, 20 years ago
*** empty log message ***
File size: 20.2 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
349 Int_t lastdesired = (Int_t)fLoGainLast;
350 Int_t lastavailable = (Int_t)fRunHeader->GetNumSamplesLoGain()-1;
351
352 if (lastdesired > lastavailable)
353 {
354 const Int_t diff = lastdesired - lastavailable;
355 *fLog << endl;
356 *fLog << warn << GetDescriptor()
357 << Form("%s%2i%s%2i%s%2i%s",": Selected Lo Gain FADC Window [",
358 (int)fLoGainFirst,",",lastdesired,
359 "] ranges out of the available limits: [0,",lastavailable,"].") << endl;
360 *fLog << GetDescriptor() << ": Will reduce the upper edge to " << (int)(fLoGainLast - diff) << endl;
361 SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff);
362 }
363
364 lastdesired = (Int_t)fHiGainLast;
365 lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
366
367 if (lastdesired > lastavailable)
368 {
369 const Int_t diff = lastdesired - lastavailable;
370 *fLog << endl;
371 *fLog << warn << GetDescriptor()
372 << Form("%s%2i%s%2i%s%2i%s",": Selected Hi Gain Range [",
373 (int)fHiGainFirst,",",lastdesired,
374 "] ranges out of the available limits: [0,",lastavailable,"].") << endl;
375 *fLog << warn << GetDescriptor()
376 << Form("%s%2i%s",": Will possibly use ",diff," samples from the Low-Gain for the High-Gain range")
377 << endl;
378 fHiGainLast -= diff;
379 fHiLoLast = diff;
380 }
381
382 lastdesired = (Int_t)fHiGainFirst+fWindowSizeHiGain-1;
383 lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
384
385 if (lastdesired > lastavailable)
386 {
387 const Int_t diff = lastdesired - lastavailable;
388 *fLog << endl;
389 *fLog << warn << GetDescriptor()
390 << Form("%s%2i%s%2i%s",": Selected Hi Gain FADC Window size ",
391 (int)fWindowSizeHiGain,
392 " ranges out of the available limits: [0,",lastavailable,"].") << endl;
393 *fLog << warn << GetDescriptor()
394 << Form("%s%2i%s",": Will use ",diff," samples from the Low-Gain for the High-Gain extraction")
395 << endl;
396
397 if ((Int_t)fWindowSizeHiGain > diff)
398 {
399 fWindowSizeHiGain -= diff;
400 fWindowSizeLoGain += diff;
401 }
402 else
403 {
404 fWindowSizeLoGain += fWindowSizeHiGain;
405 fLoGainFirst = diff-fWindowSizeHiGain;
406 fWindowSizeHiGain = 0;
407 }
408 }
409
410
411 Int_t npixels = fPedestals->GetSize();
412 Int_t areas = fPedestals->GetAverageAreas();
413 Int_t sectors = fPedestals->GetAverageSectors();
414
415 if (fSumx.GetSize()==0)
416 {
417 fSumx. Set(npixels);
418 fSumx2.Set(npixels);
419
420 fAreaSumx. Set(areas);
421 fAreaSumx2.Set(areas);
422 fAreaValid.Set(areas);
423
424 fSectorSumx. Set(sectors);
425 fSectorSumx2.Set(sectors);
426 fSectorValid.Set(sectors);
427
428 fSumx.Reset();
429 fSumx2.Reset();
430 }
431
432 if (fWindowSizeHiGain == 0 && fWindowSizeLoGain == 0)
433 {
434 *fLog << err << GetDescriptor()
435 << ": Number of extracted Slices is 0, abort ... " << endl;
436 return kFALSE;
437 }
438
439
440 *fLog << endl;
441 *fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeHiGain
442 << " HiGain FADC samples starting with slice: " << (int)fHiGainFirst << endl;
443 *fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeLoGain
444 << " LoGain FADC samples starting with slice: " << (int)fLoGainFirst << endl;
445
446
447 if (fBad)
448 {
449 const Int_t nbads = fBad->GetSize();
450 for (Int_t i=0; i<(nbads>npixels?npixels:nbads);i++)
451 if ((*fBad)[i].IsBad())
452 (*fPedestals)[i].SetValid(kFALSE);
453 }
454
455 return kTRUE;
456
457}
458// --------------------------------------------------------------------------
459//
460// Fill the MPedestalCam container with the signal mean and rms for the event.
461// Store the measured signal in arrays fSumx and fSumx2 so that we can
462// calculate the overall mean and rms in the PostProcess()
463//
464Int_t MPedCalcPedRun::Process()
465{
466
467 MRawEvtPixelIter pixel(fRawEvt);
468
469 while (pixel.Next())
470 {
471
472 const UInt_t idx = pixel.GetPixelId();
473 const UInt_t aidx = (*fGeom)[idx].GetAidx();
474 const UInt_t sector = (*fGeom)[idx].GetSector();
475
476 Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst;
477 Byte_t *end = ptr + fWindowSizeHiGain;
478
479 UInt_t sum = 0;
480 UInt_t sqr = 0;
481
482 if (fWindowSizeHiGain != 0)
483 {
484 do
485 {
486 sum += *ptr;
487 sqr += *ptr * *ptr;
488 }
489 while (++ptr != end);
490 }
491
492 if (fWindowSizeLoGain != 0)
493 {
494
495 ptr = pixel.GetLoGainSamples() + fLoGainFirst;
496 end = ptr + fWindowSizeLoGain;
497
498 do
499 {
500 sum += *ptr;
501 sqr += *ptr * *ptr;
502 }
503 while (++ptr != end);
504
505 }
506
507 const Float_t msum = (Float_t)sum;
508
509 //
510 // These three lines have been uncommented by Markus Gaug
511 // If anybody needs them, please contact me!!
512 //
513 // const Float_t higainped = msum/fNumHiGainSlices;
514 // const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSlices)/(fNumHiGainSlices-1.));
515 // (*fPedestals)[idx].Set(higainped, higainrms);
516
517 fSumx[idx] += msum;
518 fAreaSumx[aidx] += msum;
519 fSectorSumx[sector] += msum;
520 //
521 // The old version:
522 //
523 // const Float_t msqr = (Float_t)sqr;
524 // fSumx2[idx] += msqr;
525 //
526 // The new version:
527 //
528 const Float_t sqrsum = msum*msum;
529 fSumx2[idx] += sqrsum;
530 fAreaSumx2[aidx] += sqrsum;
531 fSectorSumx2[sector] += sqrsum;
532 }
533
534 fPedestals->SetReadyToSave();
535 fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain;
536
537 return kTRUE;
538}
539
540// --------------------------------------------------------------------------
541//
542// Compute signal mean and rms in the whole run and store it in MPedestalCam
543//
544Int_t MPedCalcPedRun::PostProcess()
545{
546
547 // Compute pedestals and rms from the whole run
548 const ULong_t n = fNumSamplesTot;
549 const ULong_t nevts = GetNumExecutions();
550
551 MRawEvtPixelIter pixel(fRawEvt);
552
553 while (pixel.Next())
554 {
555
556 const Int_t pixid = pixel.GetPixelId();
557 const UInt_t aidx = (*fGeom)[pixid].GetAidx();
558 const UInt_t sector = (*fGeom)[pixid].GetSector();
559
560 fAreaValid [aidx]++;
561 fSectorValid[sector]++;
562
563 const Float_t sum = fSumx.At(pixid);
564 const Float_t sum2 = fSumx2.At(pixid);
565 const Float_t higainped = sum/n;
566 //
567 // The old version:
568 //
569 // const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.));
570 //
571 // The new version:
572 //
573 // 1. Calculate the Variance of the sums:
574 Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.);
575 // 2. Scale the variance to the number of slices:
576 higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain);
577 // 3. Calculate the RMS from the Variance:
578 (*fPedestals)[pixid].Set(higainped, higainVar < 0 ? 0. : TMath::Sqrt(higainVar));
579
580 }
581
582 //
583 // Loop over the (two) area indices to get the averaged pedestal per aidx
584 //
585 for (Int_t aidx=0; aidx<fAreaValid.GetSize(); aidx++)
586 {
587
588 const Int_t napix = fAreaValid.At(aidx);
589
590 if (napix == 0)
591 continue;
592
593 const Float_t sum = fAreaSumx.At(aidx);
594 const Float_t sum2 = fAreaSumx2.At(aidx);
595 const ULong_t an = napix * n;
596 const ULong_t aevts = napix * nevts;
597
598 const Float_t higainped = sum/an;
599
600 // 1. Calculate the Variance of the sums:
601 Float_t higainVar = (sum2-sum*sum/aevts)/(aevts-1.);
602 // 2. Scale the variance to the number of slices:
603 higainVar /= fWindowSizeHiGain+fWindowSizeLoGain;
604 // 3. Calculate the RMS from the Variance:
605 Float_t higainrms = TMath::Sqrt(higainVar);
606 // 4. Re-scale it with the square root of the number of involved pixels
607 // in order to be comparable to the mean of pedRMS of that area
608 higainrms *= TMath::Sqrt((Float_t)napix);
609
610 fPedestals->GetAverageArea(aidx).Set(higainped, higainrms);
611 }
612
613 //
614 // Loop over the (six) sector indices to get the averaged pedestal per sector
615 //
616 for (Int_t sector=0; sector<fSectorValid.GetSize(); sector++)
617 {
618
619 const Int_t nspix = fSectorValid.At(sector);
620
621 if (nspix == 0)
622 continue;
623
624 const Float_t sum = fSectorSumx.At(sector);
625 const Float_t sum2 = fSectorSumx2.At(sector);
626 const ULong_t sn = nspix * n;
627 const ULong_t sevts = nspix * nevts;
628
629 const Float_t higainped = sum/sn;
630
631 // 1. Calculate the Variance of the sums:
632 Float_t higainVar = (sum2-sum*sum/sevts)/(sevts-1.);
633 // 2. Scale the variance to the number of slices:
634 higainVar /= fWindowSizeHiGain+fWindowSizeLoGain;
635 // 3. Calculate the RMS from the Variance:
636 Float_t higainrms = TMath::Sqrt(higainVar);
637 // 4. Re-scale it with the square root of the number of involved pixels
638 // in order to be comparable to the mean of pedRMS of that sector
639 higainrms *= TMath::Sqrt((Float_t)nspix);
640
641 fPedestals->GetAverageSector(sector).Set(higainped, higainrms);
642 }
643
644 fPedestals->SetTotalEntries(fNumSamplesTot);
645 fPedestals->SetReadyToSave();
646
647 return kTRUE;
648}
Note: See TracBrowser for help on using the repository browser.