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 | !
|
---|
25 | ! Copyright: MAGIC Software Development, 2000-2004
|
---|
26 | !
|
---|
27 | !
|
---|
28 | \* ======================================================================== */
|
---|
29 |
|
---|
30 | /////////////////////////////////////////////////////////////////////////////
|
---|
31 | //
|
---|
32 | // MPedCalcLoGain
|
---|
33 | //
|
---|
34 | //
|
---|
35 | // This task is devide form MPedCalcPedRun, described below. However, It
|
---|
36 | // calculates the pedstals 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 | //
|
---|
76 | // The plots show the inner and outer pixels, respectivly and have the following meaning:
|
---|
77 | //
|
---|
78 | // 1) The calculated mean pedestal per slice (from MPedCalcFromLoGain)
|
---|
79 | // 2) The fitted mean pedestal per slice (from MHPedestalCam)
|
---|
80 | // 3) The calculated pedestal RMS per slice (from MPedCalcFromLoGain)
|
---|
81 | // 4) The fitted sigma of the pedestal distribution per slice
|
---|
82 | // (from MHPedestalCam)
|
---|
83 | // 5) The relative difference between calculation and histogram fit
|
---|
84 | // for the mean
|
---|
85 | // 6) The relative difference between calculation and histogram fit
|
---|
86 | // for the sigma or RMS, respectively.
|
---|
87 | //
|
---|
88 | // The calculated means do not change significantly except for the case of 2 slices,
|
---|
89 | // however the RMS changes from 5.7 per slice in the case of 2 extracted slices
|
---|
90 | // to 8.3 per slice in the case of 26 extracted slices. This change is very significant.
|
---|
91 | //
|
---|
92 | // The plots have been produced on run 20123. You can reproduce them using
|
---|
93 | // the macro pedestalstudies.C
|
---|
94 | //
|
---|
95 | // Usage of this class:
|
---|
96 | // ====================
|
---|
97 | //
|
---|
98 | // Call: SetRange(higainfirst, higainlast, logainfirst, logainlast)
|
---|
99 | // to modify the ranges in which the window is allowed to move.
|
---|
100 | // Defaults are:
|
---|
101 | //
|
---|
102 | // fHiGainFirst = fgHiGainFirst = 0
|
---|
103 | // fHiGainLast = fgHiGainLast = 29
|
---|
104 | // fLoGainFirst = fgLoGainFirst = 0
|
---|
105 | // fLoGainLast = fgLoGainLast = 14
|
---|
106 | //
|
---|
107 | // Call: SetWindowSize(windowhigain, windowlogain)
|
---|
108 | // to modify the sliding window widths. Windows have to be an even number.
|
---|
109 | // In case of odd numbers, the window will be modified.
|
---|
110 | //
|
---|
111 | // Defaults are:
|
---|
112 | //
|
---|
113 | // fHiGainWindowSize = fgHiGainWindowSize = 14
|
---|
114 | // fLoGainWindowSize = fgLoGainWindowSize = 0
|
---|
115 | //
|
---|
116 | // Input Containers:
|
---|
117 | // MRawEvtData
|
---|
118 | // MRawRunHeader
|
---|
119 | // MGeomCam
|
---|
120 | //
|
---|
121 | // Output Containers:
|
---|
122 | // MPedestalCam
|
---|
123 | //
|
---|
124 | // See also: MPedestalCam, MPedestalPix, MHPedestalCam, MExtractor
|
---|
125 | //
|
---|
126 | /////////////////////////////////////////////////////////////////////////////
|
---|
127 | #include "MPedCalcFromLoGain.h"
|
---|
128 | #include "MExtractor.h"
|
---|
129 |
|
---|
130 | #include "MParList.h"
|
---|
131 |
|
---|
132 | #include "MLog.h"
|
---|
133 | #include "MLogManip.h"
|
---|
134 |
|
---|
135 | #include "MRawRunHeader.h"
|
---|
136 | #include "MRawEvtPixelIter.h"
|
---|
137 | #include "MRawEvtData.h"
|
---|
138 |
|
---|
139 | #include "MPedestalPix.h"
|
---|
140 | #include "MPedestalCam.h"
|
---|
141 |
|
---|
142 | #include "MGeomPix.h"
|
---|
143 | #include "MGeomCam.h"
|
---|
144 |
|
---|
145 | ClassImp(MPedCalcFromLoGain);
|
---|
146 |
|
---|
147 | using namespace std;
|
---|
148 |
|
---|
149 | const Byte_t MPedCalcFromLoGain::fgHiGainFirst = 0;
|
---|
150 | const Byte_t MPedCalcFromLoGain::fgHiGainLast = 11;
|
---|
151 | const Byte_t MPedCalcFromLoGain::fgLoGainFirst = 1;
|
---|
152 | const Byte_t MPedCalcFromLoGain::fgLoGainLast = 14;
|
---|
153 | const Byte_t MPedCalcFromLoGain::fgHiGainWindowSize = 12;
|
---|
154 | const Byte_t MPedCalcFromLoGain::fgLoGainWindowSize = 14;
|
---|
155 | const Byte_t MPedCalcFromLoGain::fgMaxHiGainVar = 20;
|
---|
156 |
|
---|
157 | // --------------------------------------------------------------------------
|
---|
158 | //
|
---|
159 | // Default constructor:
|
---|
160 | //
|
---|
161 | // Sets:
|
---|
162 | // - all pointers to NULL
|
---|
163 | // - fWindowSizeHiGain to fgHiGainWindowSize
|
---|
164 | // - fWindowSizeLoGain to fgLoGainWindowSize
|
---|
165 | //
|
---|
166 | // Calls:
|
---|
167 | // - AddToBranchList("fHiGainPixId");
|
---|
168 | // - AddToBranchList("fHiGainFadcSamples");
|
---|
169 | // - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
|
---|
170 | // - Clear()
|
---|
171 | //
|
---|
172 | MPedCalcFromLoGain::MPedCalcFromLoGain(const char *name, const char *title)
|
---|
173 | : fWindowSizeHiGain(fgHiGainWindowSize),
|
---|
174 | fWindowSizeLoGain(fgLoGainWindowSize),
|
---|
175 | fGeom(NULL)
|
---|
176 | {
|
---|
177 | fName = name ? name : "MPedCalcFromLoGain";
|
---|
178 | fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data";
|
---|
179 |
|
---|
180 | AddToBranchList("fHiGainPixId");
|
---|
181 | AddToBranchList("fHiGainFadcSamples");
|
---|
182 |
|
---|
183 | SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
|
---|
184 |
|
---|
185 | SetMaxHiGainVar(fgMaxHiGainVar);
|
---|
186 | SetPedestalUpdate(kTRUE);
|
---|
187 | Clear();
|
---|
188 | }
|
---|
189 |
|
---|
190 | // --------------------------------------------------------------------------
|
---|
191 | //
|
---|
192 | // Sets:
|
---|
193 | // - fNumSamplesTot to 0
|
---|
194 | // - fRawEvt to NULL
|
---|
195 | // - fRunHeader to NULL
|
---|
196 | // - fPedestals to NULL
|
---|
197 | //
|
---|
198 | void MPedCalcFromLoGain::Clear(const Option_t *o)
|
---|
199 | {
|
---|
200 | fRawEvt = NULL;
|
---|
201 | fRunHeader = NULL;
|
---|
202 | fPedestals = NULL;
|
---|
203 | }
|
---|
204 |
|
---|
205 | // --------------------------------------------------------------------------
|
---|
206 | //
|
---|
207 | // SetRange:
|
---|
208 | //
|
---|
209 | // Calls:
|
---|
210 | // - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
|
---|
211 | // - SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
|
---|
212 | //
|
---|
213 | void MPedCalcFromLoGain::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
|
---|
214 | {
|
---|
215 | MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
|
---|
216 |
|
---|
217 | //
|
---|
218 | // Redo the checks if the window is still inside the ranges
|
---|
219 | //
|
---|
220 | SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
|
---|
221 | }
|
---|
222 |
|
---|
223 | // --------------------------------------------------------------------------
|
---|
224 | //
|
---|
225 | void MPedCalcFromLoGain::SetMaxHiGainVar(Byte_t maxvar)
|
---|
226 | {
|
---|
227 | fMaxHiGainVar = maxvar;
|
---|
228 | }
|
---|
229 |
|
---|
230 | // --------------------------------------------------------------------------
|
---|
231 | //
|
---|
232 | // Checks:
|
---|
233 | // - if a window is odd, subtract one
|
---|
234 | // - if a window is bigger than the one defined by the ranges, set it to the available range
|
---|
235 | // - if a window is smaller than 2, set it to 2
|
---|
236 | //
|
---|
237 | // Sets:
|
---|
238 | // - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain
|
---|
239 | // - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain
|
---|
240 | // - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
|
---|
241 | // - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)
|
---|
242 | //
|
---|
243 | void MPedCalcFromLoGain::SetWindowSize(Byte_t windowh, Byte_t windowl)
|
---|
244 | {
|
---|
245 | fWindowSizeHiGain = windowh & ~1;
|
---|
246 | fWindowSizeLoGain = windowl & ~1;
|
---|
247 |
|
---|
248 | if (fWindowSizeHiGain != windowh)
|
---|
249 | *fLog << warn << GetDescriptor() << ": Hi Gain window size has to be even, set to: "
|
---|
250 | << int(fWindowSizeHiGain) << " samples " << endl;
|
---|
251 |
|
---|
252 | if (fWindowSizeLoGain != windowl)
|
---|
253 | *fLog << warn << GetDescriptor() << ": Lo Gain window size has to be even, set to: "
|
---|
254 | << int(fWindowSizeLoGain) << " samples " << endl;
|
---|
255 |
|
---|
256 | const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1;
|
---|
257 | const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1;
|
---|
258 |
|
---|
259 | if (fWindowSizeHiGain > availhirange)
|
---|
260 | {
|
---|
261 | *fLog << warn << GetDescriptor()
|
---|
262 | << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",(int)fWindowSizeHiGain,
|
---|
263 | " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl;
|
---|
264 | *fLog << warn << GetDescriptor()
|
---|
265 | << ": Will set window size to: " << (int)availhirange << endl;
|
---|
266 | fWindowSizeHiGain = availhirange;
|
---|
267 | }
|
---|
268 |
|
---|
269 | if (fWindowSizeLoGain > availlorange)
|
---|
270 | {
|
---|
271 | *fLog << warn << GetDescriptor()
|
---|
272 | << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",(int)fWindowSizeLoGain,
|
---|
273 | " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl;
|
---|
274 | *fLog << warn << GetDescriptor()
|
---|
275 | << ": Will set window size to: " << (int)availlorange << endl;
|
---|
276 | fWindowSizeLoGain = availlorange;
|
---|
277 | }
|
---|
278 |
|
---|
279 | fNumHiGainSamples = (Float_t)fWindowSizeHiGain;
|
---|
280 | fNumLoGainSamples = (Float_t)fWindowSizeLoGain;
|
---|
281 |
|
---|
282 | fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
|
---|
283 | fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
|
---|
284 | }
|
---|
285 |
|
---|
286 | // --------------------------------------------------------------------------
|
---|
287 | //
|
---|
288 | // Look for the following input containers:
|
---|
289 | //
|
---|
290 | // - MRawEvtData
|
---|
291 | // - MRawRunHeader
|
---|
292 | // - MGeomCam
|
---|
293 | //
|
---|
294 | // The following output containers are also searched and created if
|
---|
295 | // they were not found:
|
---|
296 | //
|
---|
297 | // - MPedestalCam
|
---|
298 | //
|
---|
299 | Int_t MPedCalcFromLoGain::PreProcess( MParList *pList )
|
---|
300 | {
|
---|
301 | Clear();
|
---|
302 |
|
---|
303 | fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
|
---|
304 | if (!fRawEvt)
|
---|
305 | {
|
---|
306 | *fLog << err << "MRawEvtData not found... aborting." << endl;
|
---|
307 | return kFALSE;
|
---|
308 | }
|
---|
309 |
|
---|
310 | fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader"));
|
---|
311 | if (!fRunHeader)
|
---|
312 | {
|
---|
313 | *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl;
|
---|
314 | return kFALSE;
|
---|
315 | }
|
---|
316 |
|
---|
317 | fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
|
---|
318 | if (!fGeom)
|
---|
319 | {
|
---|
320 | *fLog << err << "MGeomCam not found... aborting." << endl;
|
---|
321 | return kFALSE;
|
---|
322 | }
|
---|
323 |
|
---|
324 | fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam");
|
---|
325 | if (!fPedestals)
|
---|
326 | return kFALSE;
|
---|
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 | Bool_t MPedCalcFromLoGain::ReInit(MParList *pList)
|
---|
348 | {
|
---|
349 |
|
---|
350 | Int_t lastdesired = (Int_t)fLoGainLast;
|
---|
351 | Int_t lastavailable = (Int_t)fRunHeader->GetNumSamplesLoGain()-1;
|
---|
352 |
|
---|
353 | if (lastdesired > lastavailable)
|
---|
354 | {
|
---|
355 | const Int_t diff = lastdesired - lastavailable;
|
---|
356 | *fLog << endl;
|
---|
357 | *fLog << warn << GetDescriptor()
|
---|
358 | << Form("%s%2i%s%2i%s%2i%s",": Selected Lo Gain FADC Window [",
|
---|
359 | (int)fLoGainFirst,",",lastdesired,
|
---|
360 | "] ranges out of the available limits: [0,",lastavailable,"].") << endl;
|
---|
361 | *fLog << GetDescriptor() << ": Will reduce the upper edge to " << (int)(fLoGainLast - diff) << endl;
|
---|
362 | SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff);
|
---|
363 | }
|
---|
364 |
|
---|
365 | lastdesired = (Int_t)fHiGainLast;
|
---|
366 | lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
|
---|
367 |
|
---|
368 | if (lastdesired > lastavailable)
|
---|
369 | {
|
---|
370 | const Int_t diff = lastdesired - lastavailable;
|
---|
371 | *fLog << endl;
|
---|
372 | *fLog << warn << GetDescriptor()
|
---|
373 | << Form("%s%2i%s%2i%s%2i%s",": Selected Hi Gain Range [",
|
---|
374 | (int)fHiGainFirst,",",lastdesired,
|
---|
375 | "] ranges out of the available limits: [0,",lastavailable,"].") << endl;
|
---|
376 | *fLog << warn << GetDescriptor()
|
---|
377 | << Form("%s%2i%s",": Will possibly use ",diff," samples from the Low-Gain for the High-Gain range")
|
---|
378 | << endl;
|
---|
379 | fHiGainLast -= diff;
|
---|
380 | fHiLoLast = diff;
|
---|
381 | }
|
---|
382 |
|
---|
383 | lastdesired = (Int_t)fHiGainFirst+fWindowSizeHiGain-1;
|
---|
384 | lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1;
|
---|
385 |
|
---|
386 | if (lastdesired > lastavailable)
|
---|
387 | {
|
---|
388 | const Int_t diff = lastdesired - lastavailable;
|
---|
389 | *fLog << endl;
|
---|
390 | *fLog << warn << GetDescriptor()
|
---|
391 | << Form("%s%2i%s%2i%s",": Selected Hi Gain FADC Window size ",
|
---|
392 | (int)fWindowSizeHiGain,
|
---|
393 | " ranges out of the available limits: [0,",lastavailable,"].") << endl;
|
---|
394 | *fLog << warn << GetDescriptor()
|
---|
395 | << Form("%s%2i%s",": Will use ",diff," samples from the Low-Gain for the High-Gain extraction")
|
---|
396 | << endl;
|
---|
397 |
|
---|
398 | if ((Int_t)fWindowSizeHiGain > diff)
|
---|
399 | {
|
---|
400 | fWindowSizeHiGain -= diff;
|
---|
401 | fWindowSizeLoGain += diff;
|
---|
402 | }
|
---|
403 | else
|
---|
404 | {
|
---|
405 | fWindowSizeLoGain += fWindowSizeHiGain;
|
---|
406 | fLoGainFirst = diff-fWindowSizeHiGain;
|
---|
407 | fWindowSizeHiGain = 0;
|
---|
408 | }
|
---|
409 | }
|
---|
410 |
|
---|
411 |
|
---|
412 | Int_t npixels = fPedestals->GetSize();
|
---|
413 |
|
---|
414 | if (fSumx.GetSize()==0)
|
---|
415 | {
|
---|
416 | fSumx. Set(npixels);
|
---|
417 | fSumx2.Set(npixels);
|
---|
418 | fSumAB0.Set(npixels);
|
---|
419 | fSumAB1.Set(npixels);
|
---|
420 | fNumEventsUsed.Set(npixels);
|
---|
421 | fTotalCounter.Set(npixels);
|
---|
422 |
|
---|
423 | fSumx.Reset();
|
---|
424 | fSumx2.Reset();
|
---|
425 | fSumAB0.Reset();
|
---|
426 | fSumAB1.Reset();
|
---|
427 | fNumEventsUsed.Reset();
|
---|
428 | fTotalCounter.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 | return kTRUE;
|
---|
446 | }
|
---|
447 |
|
---|
448 | // --------------------------------------------------------------------------
|
---|
449 | //
|
---|
450 | // Fill the MPedestalCam container with the signal mean and rms for the event.
|
---|
451 | // Store the measured signal in arrays fSumx and fSumx2 so that we can
|
---|
452 | // calculate the overall mean and rms in the PostProcess()
|
---|
453 | //
|
---|
454 | Int_t MPedCalcFromLoGain::Process()
|
---|
455 | {
|
---|
456 | MRawEvtPixelIter pixel(fRawEvt);
|
---|
457 |
|
---|
458 | while (pixel.Next()) {
|
---|
459 |
|
---|
460 | const UInt_t idx = pixel.GetPixelId();
|
---|
461 |
|
---|
462 | Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst;
|
---|
463 | Byte_t *end = ptr + fWindowSizeHiGain;
|
---|
464 |
|
---|
465 | UInt_t sum = 0;
|
---|
466 | UInt_t sqr = 0;
|
---|
467 |
|
---|
468 | UInt_t max = 0;
|
---|
469 | UInt_t min = 255;
|
---|
470 |
|
---|
471 | // Find the maximum and minimum signal per slice in the high gain window
|
---|
472 | do {
|
---|
473 | if (*ptr > max) {
|
---|
474 | max = *ptr;
|
---|
475 | }
|
---|
476 | if (*ptr < min) {
|
---|
477 | min = *ptr;
|
---|
478 | }
|
---|
479 | } while (++ptr != end);
|
---|
480 |
|
---|
481 | // If the maximum in the high gain window is smaller than
|
---|
482 | if ((max-min < fMaxHiGainVar) && (max < 255)) {
|
---|
483 |
|
---|
484 | Byte_t *firstSlice = pixel.GetLoGainSamples() + fLoGainFirst;
|
---|
485 | Byte_t *lastSlice = firstSlice + fWindowSizeLoGain;
|
---|
486 |
|
---|
487 | Byte_t *slice = firstSlice;
|
---|
488 | do {
|
---|
489 | sum += *slice;
|
---|
490 | sqr += *slice * *slice;
|
---|
491 | } while (++slice != lastSlice);
|
---|
492 |
|
---|
493 | const Float_t msum = (Float_t)sum;
|
---|
494 | const Float_t sqrsum = msum*msum;
|
---|
495 |
|
---|
496 | fSumx[idx] += msum;
|
---|
497 | fSumx2[idx] += sqrsum;
|
---|
498 | fNumEventsUsed[idx]++;
|
---|
499 |
|
---|
500 | // Calculate the amplitude of the 150MHz "AB" noise
|
---|
501 |
|
---|
502 | Int_t abFlag = (fRunHeader->GetNumSamplesHiGain()
|
---|
503 | + fLoGainFirst
|
---|
504 | + pixel.HasABFlag()) & 0x1;
|
---|
505 | for (Int_t islice=0; islice<fWindowSizeLoGain; islice+=2) {
|
---|
506 | Int_t sliceAB0 = islice + abFlag;
|
---|
507 | Int_t sliceAB1 = islice + 1 - abFlag;
|
---|
508 | fSumAB0[idx] += firstSlice[sliceAB0];
|
---|
509 | fSumAB1[idx] += firstSlice[sliceAB1];
|
---|
510 | }
|
---|
511 |
|
---|
512 | if (fPedestalUpdate && (fNumEventsUsed[idx] == fNumEventsDump)) {
|
---|
513 |
|
---|
514 | const ULong_t n = fNumEventsDump;
|
---|
515 |
|
---|
516 | const ULong_t nsamplestot = n*fWindowSizeLoGain;
|
---|
517 |
|
---|
518 | const Float_t sum = fSumx.At(idx);
|
---|
519 | const Float_t sum2 = fSumx2.At(idx);
|
---|
520 | const Float_t ped = sum/(nsamplestot);
|
---|
521 |
|
---|
522 | // 1. Calculate the Variance of the sums:
|
---|
523 | Float_t var = (sum2-sum*sum/n)/(n-1.);
|
---|
524 |
|
---|
525 | // 2. Scale the variance to the number of slices:
|
---|
526 | var /= (Float_t)(fWindowSizeLoGain);
|
---|
527 |
|
---|
528 | // 3. Calculate the RMS from the Variance:
|
---|
529 | Float_t rms = var < 0 ? 0. : TMath::Sqrt(var);
|
---|
530 |
|
---|
531 | // 4. Calculate the amplitude of the 150MHz "AB" noise
|
---|
532 | Float_t abOffs = (fSumAB0[idx] - fSumAB1[idx]) / nsamplestot;
|
---|
533 |
|
---|
534 | (*fPedestals)[idx].Set(ped, rms, abOffs, n);
|
---|
535 |
|
---|
536 | fTotalCounter[idx]++;
|
---|
537 | fNumEventsUsed[idx]=0;
|
---|
538 | fSumx[idx]=0;
|
---|
539 | fSumx2[idx]=0;
|
---|
540 | fSumAB0[idx]=0;
|
---|
541 | fSumAB1[idx]=0;
|
---|
542 | }
|
---|
543 | }
|
---|
544 | }
|
---|
545 |
|
---|
546 | fPedestals->SetReadyToSave();
|
---|
547 |
|
---|
548 | return kTRUE;
|
---|
549 | }
|
---|
550 |
|
---|
551 | // --------------------------------------------------------------------------
|
---|
552 | //
|
---|
553 | // Compute signal mean and rms in the whole run and store it in MPedestalCam
|
---|
554 | //
|
---|
555 | Int_t MPedCalcFromLoGain::PostProcess()
|
---|
556 | {
|
---|
557 |
|
---|
558 | // Compute pedestals and rms from the whole run
|
---|
559 |
|
---|
560 | if (!fPedestalUpdate) {
|
---|
561 |
|
---|
562 | MRawEvtPixelIter pixel(fRawEvt);
|
---|
563 |
|
---|
564 | while (pixel.Next()) {
|
---|
565 |
|
---|
566 | const Int_t idx = pixel.GetPixelId();
|
---|
567 |
|
---|
568 | const ULong_t n = fNumEventsUsed[idx];
|
---|
569 |
|
---|
570 | if (n < 2)
|
---|
571 | continue;
|
---|
572 |
|
---|
573 | const ULong_t nsamplestot = n*fWindowSizeLoGain;
|
---|
574 |
|
---|
575 | const Float_t sum = fSumx.At(idx);
|
---|
576 | const Float_t sum2 = fSumx2.At(idx);
|
---|
577 | const Float_t ped = sum/(nsamplestot);
|
---|
578 |
|
---|
579 | // 1. Calculate the Variance of the sums:
|
---|
580 | Float_t var = (sum2-sum*sum/n)/(n-1.);
|
---|
581 |
|
---|
582 | // 2. Scale the variance to the number of slices:
|
---|
583 | var /= (Float_t)(fWindowSizeLoGain);
|
---|
584 |
|
---|
585 | // 3. Calculate the RMS from the Variance:
|
---|
586 | Float_t rms = var < 0 ? 0. : TMath::Sqrt(var);
|
---|
587 |
|
---|
588 | // 4. Calculate the amplitude of the 150MHz "AB" noise
|
---|
589 | Float_t abOffs = (fSumAB0[idx] - fSumAB1[idx]) / nsamplestot;
|
---|
590 |
|
---|
591 | (*fPedestals)[idx].Set(ped, rms, abOffs, n);
|
---|
592 |
|
---|
593 | fTotalCounter[idx]++;
|
---|
594 | }
|
---|
595 |
|
---|
596 | fPedestals->SetReadyToSave();
|
---|
597 | }
|
---|
598 |
|
---|
599 | return kTRUE;
|
---|
600 | }
|
---|