source: trunk/Mars/msignal/MExtractTimeAndChargeSlidingWindow.cc@ 11876

Last change on this file since 11876 was 7069, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 14.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 analyzing 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! Author(s): Markus Gaug, 09/2004 <mailto:markus@ifae.es>
18! Author(s): Hendrik Bartko, 01/2004 <mailto:hbartko@mppmu.mpg.de>
19!
20! Copyright: MAGIC Software Development, 2002-2005
21!
22!
23\* ======================================================================== */
24//////////////////////////////////////////////////////////////////////////////
25//
26// MExtractSlidingWindow
27//
28// Extracts the signal from a sliding window of size fHiGainWindowSize and
29// fLoGainWindowSize, respectively. The signal is the one which maximizes
30// the clock-noise and pedestal-corrected integral contents.
31//
32// The amplitude-weighted arrival time is calculated from the window with
33// the highest integral using the following formula:
34//
35// t = sum(s(i) * i) / sum(i)
36//
37// where i denotes the FADC slice index and s(i) the clock-noise and
38/// pedestal-corrected FADC value at slice index i. The sum runs over the
39// extraction window size.
40//
41// Call: SetRange(higainfirst, higainlast, logainfirst, logainlast)
42// to modify the ranges in which the window is allowed to move.
43//
44// Defaults are:
45//
46// fHiGainFirst = fgHiGainFirst = 0
47// fHiGainLast = fgHiGainLast = 14
48// fLoGainFirst = fgLoGainFirst = 2
49// fLoGainLast = fgLoGainLast = 14
50//
51// Call: SetWindowSize(windowhigain, windowlogain)
52// to modify the sliding window widths. Windows have to be an even number.
53// Odd numbers are possible since the clock-noise is corrected for.
54//
55// Defaults are:
56//
57// fHiGainWindowSize = 6
58// fLoGainWindowSize = 6
59//
60//////////////////////////////////////////////////////////////////////////////
61#include "MExtractTimeAndChargeSlidingWindow.h"
62
63#include "MPedestalPix.h"
64
65#include "MLog.h"
66#include "MLogManip.h"
67
68ClassImp(MExtractTimeAndChargeSlidingWindow);
69
70using namespace std;
71
72const Byte_t MExtractTimeAndChargeSlidingWindow::fgHiGainFirst = 0;
73const Byte_t MExtractTimeAndChargeSlidingWindow::fgHiGainLast = 16;
74const Byte_t MExtractTimeAndChargeSlidingWindow::fgLoGainFirst = 2;
75const Byte_t MExtractTimeAndChargeSlidingWindow::fgLoGainLast = 14;
76const Byte_t MExtractTimeAndChargeSlidingWindow::fgHiGainWindowSize = 6;
77const Byte_t MExtractTimeAndChargeSlidingWindow::fgLoGainWindowSize = 8;
78
79// --------------------------------------------------------------------------
80//
81// Default constructor.
82//
83// Calls:
84// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
85//
86// Initializes:
87// - fWindowSizeHiGain to fgHiGainWindowSize
88// - fWindowSizeLoGain to fgLoGainWindowSize
89//
90MExtractTimeAndChargeSlidingWindow::MExtractTimeAndChargeSlidingWindow(const char *name, const char *title)
91{
92
93 fName = name ? name : "MExtractTimeAndChargeSlidingWindow";
94 fTitle = title ? title : "Calculate arrival times and charges using a sliding window";
95
96 fWindowSizeHiGain = fgHiGainWindowSize;
97 fWindowSizeLoGain = fgLoGainWindowSize;
98
99 SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
100}
101
102//-------------------------------------------------------------------
103//
104// Set the ranges
105//
106// Calls:
107// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
108// - SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
109//
110void MExtractTimeAndChargeSlidingWindow::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
111{
112
113 MExtractor::SetRange(hifirst, hilast, lofirst, lolast);
114
115 //
116 // Redo the checks if the window is still inside the ranges
117 //
118 SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
119
120}
121
122// -----------------------------------------------------------------------------------------
123//
124// Checks:
125// - if a window is bigger than the one defined by the ranges, set it to the available range
126// - if a window is smaller than 2, set it to 2
127//
128// Sets:
129// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain
130// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain
131// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
132// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)
133//
134void MExtractTimeAndChargeSlidingWindow::SetWindowSize(Int_t windowh, Int_t windowl)
135{
136
137 fWindowSizeHiGain = windowh;
138 fWindowSizeLoGain = windowl;
139
140 const Int_t availhirange = (Int_t)(fHiGainLast-fHiGainFirst+1);
141
142 if (fWindowSizeHiGain > availhirange)
143 {
144 *fLog << warn << GetDescriptor()
145 << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",(int)fWindowSizeHiGain,
146 " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl;
147 *fLog << warn << GetDescriptor()
148 << ": Will set window size to: " << (int)availhirange << endl;
149 fWindowSizeHiGain = availhirange;
150 }
151
152 if (fWindowSizeHiGain<1)
153 {
154 fWindowSizeHiGain = 1;
155 *fLog << warn << GetDescriptor() << ": High Gain window size too small, set to one sample" << endl;
156 }
157
158 if (fLoGainLast != 0 && fWindowSizeLoGain != 0)
159 {
160 const Int_t availlorange = (Int_t)(fLoGainLast-fLoGainFirst+1);
161
162 if (fWindowSizeLoGain > availlorange)
163 {
164 *fLog << warn << GetDescriptor()
165 << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",(int)fWindowSizeLoGain,
166 " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl;
167 *fLog << warn << GetDescriptor()
168 << ": Will set window size to: " << (int)availlorange << endl;
169 fWindowSizeLoGain = availlorange;
170 }
171 }
172
173 fNumHiGainSamples = (Float_t)fWindowSizeHiGain;
174 fNumLoGainSamples = fLoGainLast ? (Float_t)fWindowSizeLoGain : 0.;
175
176 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
177 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
178
179 switch (fWindowSizeHiGain)
180 {
181 case 2:
182 SetResolutionPerPheHiGain(0.050);
183 break;
184 case 4:
185 SetResolutionPerPheHiGain(0.039);
186 break;
187 case 6:
188 case 8:
189 SetResolutionPerPheHiGain(0.011);
190 break;
191 case 14:
192 SetResolutionPerPheHiGain(0.009);
193 break;
194 default:
195 *fLog << warn << GetDescriptor() << ": Could not set the high-gain extractor resolution per phe for window size "
196 << fWindowSizeHiGain << endl;
197 }
198
199 switch (fWindowSizeLoGain)
200 {
201 case 2:
202 SetResolutionPerPheLoGain(0.028);
203 break;
204 case 4:
205 SetResolutionPerPheLoGain(0.013);
206 break;
207 case 6:
208 SetResolutionPerPheLoGain(0.008);
209 break;
210 case 8:
211 case 10:
212 SetResolutionPerPheLoGain(0.005);
213 break;
214 default:
215 *fLog << warn << GetDescriptor() << ": Could not set the low-gain extractor resolution per phe for window size "
216 << fWindowSizeLoGain << endl;
217 }
218}
219
220// --------------------------------------------------------------------------
221//
222// InitArrays
223//
224// Gets called in the ReInit() and initialized the arrays
225//
226Bool_t MExtractTimeAndChargeSlidingWindow::InitArrays()
227{
228 Int_t range = (Int_t)(fHiGainLast - fHiGainFirst + 1 + fHiLoLast);
229 fHiGainSignal.Set(range);
230 range = (Int_t)(fLoGainLast - fLoGainFirst + 1);
231 fLoGainSignal.Set(range);
232
233 return kTRUE;
234
235}
236
237// --------------------------------------------------------------------------
238//
239// Calculates the arrival time for each pixel
240//
241void MExtractTimeAndChargeSlidingWindow::FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum,
242 Float_t &time, Float_t &dtime,
243 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
244{
245
246 Int_t range = fHiGainLast - fHiGainFirst + 1;
247 const Byte_t *end = first + range;
248 Byte_t *p = first;
249
250 Float_t max = 0; // highest integral content of all windows
251 sat = 0;
252
253 const Float_t pedes = ped.GetPedestal();
254 const Float_t ABoffs = ped.GetPedestalABoffset();
255
256 const Float_t PedMean[2] = { pedes + ABoffs, pedes - ABoffs };
257
258 fMaxBinContent = 0;
259 //
260 // Check for saturation in all other slices
261 //
262 Int_t ids = fHiGainFirst;
263 Float_t *sample = fHiGainSignal.GetArray();
264
265 while (p<first+fWindowSizeHiGain)
266 {
267
268 const Float_t signal = (Float_t)*p - PedMean[(ids++ + abflag) & 0x1];
269 sum += signal;
270 *sample++ = signal;
271
272 if (*p > fMaxBinContent)
273 fMaxBinContent = *p;
274
275 if (*p++ >= fSaturationLimit)
276 if (!sat)
277 sat = ids-2;
278 }
279
280 if (IsNoiseCalculation())
281 return;
282
283 //
284 // Check for saturation in all other slices
285 //
286 while (p<end)
287 {
288 if (*p > fMaxBinContent)
289 fMaxBinContent = *p;
290
291 if (*p++ >= fSaturationLimit)
292 if (!sat)
293 sat = ids-2;
294 }
295
296 //
297 // Calculate the i-th sum as
298 // sum_i+1 = sum_i + slice[i+8] - slice[i]
299 // This is fast and accurate (because we are using int's)
300 //
301 Int_t count = 0;
302 max = sum;
303 Int_t idx = 0; // idx of the first slice of the maximum window
304
305 for (p=first; p+fWindowSizeHiGain<end; p++)
306 {
307
308 const Float_t signal = (Float_t)*(p+fWindowSizeHiGain) - PedMean[(ids++ + abflag) & 0x1];
309 sum += signal - *(sample-fWindowSizeHiGain);
310 *sample++ = signal;
311
312 if (sum>max)
313 {
314 max = sum;
315 idx = count+1;
316 }
317 count++;
318 }
319
320 //
321 // overlap bins
322 //
323 Byte_t *l = logain;
324 while (l < logain+fHiLoLast)
325 {
326
327 const Float_t signal = (Float_t)*l - PedMean[(ids++ + abflag) & 0x1];
328 sum += signal - *(sample-fWindowSizeHiGain);
329 *sample++ = signal;
330
331 if (*l > fMaxBinContent)
332 fMaxBinContent = *logain;
333
334 if (*l++ >= fSaturationLimit)
335 if (!sat)
336 sat = ids-2;
337
338 if (sum>max)
339 {
340 max = sum;
341 idx = count+1;
342 }
343 count++;
344 } /* while (l<logain+fHiLoLast) */
345
346 //
347 // now calculate the time for the maximum window
348 //
349 Float_t timesignalsum = 0.;
350 Int_t timesquaredsum = 0;
351
352 for (Int_t i=idx; i<idx+fWindowSizeHiGain; i++)
353 {
354 timesignalsum += fHiGainSignal[i]*i;
355 timesquaredsum += i*i;
356 }
357
358 sum = max;
359
360 time = max > 0.1 ? timesignalsum / max + Float_t(fHiGainFirst) : 1.;
361 dtime = max > 0.1 ? ped.GetPedestalRms() / max * sqrt(timesquaredsum - fWindowSizeHiGain*time) : 1.;
362
363}
364
365
366// --------------------------------------------------------------------------
367//
368// Calculates the arrival time for each pixel
369//
370void MExtractTimeAndChargeSlidingWindow::FindTimeAndChargeLoGain(Byte_t *first, Float_t &sum, Float_t &dsum,
371 Float_t &time, Float_t &dtime,
372 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
373{
374
375 Int_t range = fLoGainLast - fLoGainFirst + 1;
376 const Byte_t *end = first + range;
377 Byte_t *p = first;
378
379 Float_t max = 0; // highest integral content of all windows
380 Int_t count = 0; // counter to recognize the AB-flag
381
382 Float_t pedes = ped.GetPedestal();
383 const Float_t ABoffs = ped.GetPedestalABoffset();
384
385 Float_t PedMean[2] = { pedes + ABoffs, pedes - ABoffs };
386 //
387 // Check for saturation in all other slices
388 //
389 Int_t ids = fLoGainFirst;
390
391 while (p<first+fWindowSizeLoGain)
392 {
393 const Float_t signal = (Float_t)*p - PedMean[(ids++ + abflag) & 0x1];
394 sum += signal;
395 fLoGainSignal[count] = signal;
396
397 if (*p++ >= fSaturationLimit)
398 sat++;
399
400 count++;
401 }
402
403 //
404 // Check for saturation in all other slices
405 //
406 while (p<end)
407 if (*p++ >= fSaturationLimit)
408 sat++;
409
410 if (IsNoiseCalculation())
411 return;
412
413 //
414 // Calculate the i-th sum as
415 // sum_i+1 = sum_i + slice[i+8] - slice[i]
416 // This is fast and accurate (because we are using int's)
417 //
418 count = 0;
419 max = sum;
420 Int_t idx = 0; // idx of the first slice of the maximum window
421
422 for (p=first; p+fWindowSizeLoGain<end; p++)
423 {
424
425 const Float_t signal = (Float_t)*(p+fWindowSizeLoGain) - PedMean[(ids++ + abflag) & 0x1];
426 sum += signal - fLoGainSignal[count];
427 fLoGainSignal[count + fWindowSizeLoGain] = signal;
428
429 if (sum>max)
430 {
431 max = sum;
432 idx = count+1;
433 }
434 count++;
435 }
436
437 //
438 // now calculate the time for the maximum window
439 //
440 Float_t timesignalsum = 0;
441 Int_t timesquaredsum = 0;
442
443 for (Int_t i=idx; i<idx+fWindowSizeLoGain; i++)
444 {
445 timesignalsum += fLoGainSignal[i]*i;
446 timesquaredsum += i*i;
447 }
448
449 sum = max;
450
451 time = max > 0.1 ? timesignalsum / max + Float_t(fLoGainFirst) : 1.;
452 dtime = max > 0.1 ? ped.GetPedestalRms() / max * sqrt(timesquaredsum - fWindowSizeLoGain*time) : 1.;
453}
454
455// --------------------------------------------------------------------------
456//
457// In addition to the resources of the base-class MExtractor:
458// MJPedestal.MExtractor.WindowSizeHiGain: 6
459// MJPedestal.MExtractor.WindowSizeLoGain: 6
460//
461Int_t MExtractTimeAndChargeSlidingWindow::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
462{
463
464 Byte_t hw = fWindowSizeHiGain;
465 Byte_t lw = fWindowSizeLoGain;
466
467 Bool_t rc = kFALSE;
468
469 if (IsEnvDefined(env, prefix, "HiGainWindowSize", print))
470 {
471 hw = GetEnvValue(env, prefix, "HiGainWindowSize", hw);
472 rc = kTRUE;
473 }
474 if (IsEnvDefined(env, prefix, "LoGainWindowSize", print))
475 {
476 lw = GetEnvValue(env, prefix, "LoGainWindowSize", lw);
477 rc = kTRUE;
478 }
479
480 if (rc)
481 SetWindowSize(hw, lw);
482
483 return MExtractTimeAndCharge::ReadEnv(env, prefix, print) ? kTRUE : rc;
484
485}
486
Note: See TracBrowser for help on using the repository browser.