source: trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSlidingWindow.cc@ 5712

Last change on this file since 5712 was 5511, checked in by gaug, 20 years ago
*** empty log message ***
File size: 14.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 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): Thomas Bretz, 02/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
18! Author(s): Hendrik Bartko, 01/2004 <mailto:hbartko@mppmu.mpg.de>
19! Author(s): Markus Gaug, 09/2004 <mailto:markus@ifae.es>
20!
21! Copyright: MAGIC Software Development, 2002-2004
22!
23!
24\* ======================================================================== */
25//////////////////////////////////////////////////////////////////////////////
26//
27// MExtractSlidingWindow
28//
29// Extracts the signal from a sliding window of size fHiGainWindowSize and
30// fLoGainWindowSize, respectively. The signal is the one which maximizes
31// the clock-noise and pedestal-corrected integral contents.
32//
33// The amplitude-weighted arrival time is calculated from the window with
34// the highest integral using the following formula:
35//
36// t = sum(s(i) * i) / sum(i)
37//
38// where i denotes the FADC slice index and s(i) the clock-noise and
39/// pedestal-corrected FADC value at slice index i. The sum runs over the
40// extraction window size.
41//
42// Call: SetRange(higainfirst, higainlast, logainfirst, logainlast)
43// to modify the ranges in which the window is allowed to move.
44//
45// Defaults are:
46//
47// fHiGainFirst = fgHiGainFirst = 0
48// fHiGainLast = fgHiGainLast = 14
49// fLoGainFirst = fgLoGainFirst = 2
50// fLoGainLast = fgLoGainLast = 14
51//
52// Call: SetWindowSize(windowhigain, windowlogain)
53// to modify the sliding window widths. Windows have to be an even number.
54// Odd numbers are possible since the clock-noise is corrected for.
55//
56// Defaults are:
57//
58// fHiGainWindowSize = 6
59// fLoGainWindowSize = 6
60//
61//////////////////////////////////////////////////////////////////////////////
62#include "MExtractTimeAndChargeSlidingWindow.h"
63
64#include "MPedestalPix.h"
65
66#include "MLog.h"
67#include "MLogManip.h"
68
69ClassImp(MExtractTimeAndChargeSlidingWindow);
70
71using namespace std;
72
73const Byte_t MExtractTimeAndChargeSlidingWindow::fgHiGainFirst = 2;
74const Byte_t MExtractTimeAndChargeSlidingWindow::fgHiGainLast = 16;
75const Byte_t MExtractTimeAndChargeSlidingWindow::fgLoGainFirst = 2;
76const Byte_t MExtractTimeAndChargeSlidingWindow::fgLoGainLast = 14;
77const Byte_t MExtractTimeAndChargeSlidingWindow::fgHiGainWindowSize = 6;
78const Byte_t MExtractTimeAndChargeSlidingWindow::fgLoGainWindowSize = 6;
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}
180
181// --------------------------------------------------------------------------
182//
183// InitArrays
184//
185// Gets called in the ReInit() and initialized the arrays
186//
187Bool_t MExtractTimeAndChargeSlidingWindow::InitArrays()
188{
189 Int_t range = (Int_t)(fHiGainLast - fHiGainFirst + 1 + fHiLoLast);
190 fHiGainSignal.Set(range);
191 range = (Int_t)(fLoGainLast - fLoGainFirst + 1);
192 fLoGainSignal.Set(range);
193
194 return kTRUE;
195
196}
197
198// --------------------------------------------------------------------------
199//
200// Calculates the arrival time for each pixel
201//
202void MExtractTimeAndChargeSlidingWindow::FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum,
203 Float_t &time, Float_t &dtime,
204 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
205{
206
207 Int_t range = fHiGainLast - fHiGainFirst + 1;
208 const Byte_t *end = first + range;
209 Byte_t *p = first;
210
211 Float_t max = 0; // highest integral content of all windows
212 Int_t count = 0; // counter to recognize the AB-flag
213
214 const Float_t pedes = ped.GetPedestal();
215 const Float_t ABoffs = ped.GetPedestalABoffset();
216
217 Float_t PedMean[2];
218 PedMean[0] = pedes + ABoffs;
219 PedMean[1] = pedes - ABoffs;
220
221 //
222 // Check for saturation in all other slices
223 //
224 while (p<first+fWindowSizeHiGain)
225 {
226 const Int_t ids = fHiGainFirst + count;
227 const Float_t signal = (Float_t)*p - PedMean[(ids+abflag) & 0x1];
228 sum += signal;
229 fHiGainSignal[count] = signal;
230
231 if (*p++ >= fSaturationLimit)
232 sat++;
233
234 count++;
235 }
236
237 //
238 // Check for saturation in all other slices
239 //
240 while (p<end)
241 if (*p++ >= fSaturationLimit)
242 sat++;
243
244 if (IsNoiseCalculation())
245 return;
246
247 //
248 // Calculate the i-th sum as
249 // sum_i+1 = sum_i + slice[i+8] - slice[i]
250 // This is fast and accurate (because we are using int's)
251 //
252 count = 0;
253 max = sum;
254 Int_t idx = 0; // idx of the first slice of the maximum window
255
256 for (p=first; p+fWindowSizeHiGain<end; p++)
257 {
258
259 const Int_t ids = fHiGainFirst + count + fWindowSizeHiGain;
260 const Float_t signal = (Float_t)*(p+fWindowSizeHiGain) - PedMean[(ids+abflag) & 0x1];
261 sum += signal - fHiGainSignal[count];
262 fHiGainSignal[count + fWindowSizeHiGain] = signal;
263
264 if (sum>max)
265 {
266 max = sum;
267 idx = count+1;
268 }
269 count++;
270 }
271
272 if (fHiLoLast != 0)
273 {
274
275 //
276 // overlap bins
277 //
278 Byte_t *l = logain;
279
280 while (p < end && l < logain+fHiLoLast)
281 {
282
283 const Int_t ids = fHiGainFirst + count + fWindowSizeHiGain;
284 const Float_t signal = (Float_t)*l - PedMean[(ids+abflag) & 0x1];
285 sum += signal - fHiGainSignal[count];
286 fHiGainSignal[count + fWindowSizeHiGain] = signal;
287 count++;
288 if (*l++ >= fSaturationLimit)
289 sat++;
290
291 if (sum>max)
292 {
293 max = sum;
294 idx = count+1;
295 }
296 p++;
297 }
298
299 if (fHiLoLast > fWindowSizeHiGain)
300 {
301 while (l < logain + fHiLoLast)
302 {
303 const Int_t ids = fHiGainFirst + count + fWindowSizeHiGain;
304 const Float_t signal = (Float_t)*l - PedMean[(ids+abflag) & 0x1];
305 sum += signal - fHiGainSignal[count];
306 fHiGainSignal[count+fWindowSizeHiGain] = signal;
307 count++;
308 if (*(l+fWindowSizeHiGain) >= fSaturationLimit)
309 sat++;
310
311 if (sum>max)
312 {
313 max = sum;
314 idx = count+1;
315 }
316 l++;
317 } /* while (l < logain + fHiLoLast) */
318 } /* if (fHiLoLast > fWindowSizeHiGain) */
319 } /* if (fHiLoLast != 0) */
320
321 //
322 // now calculate the time for the maximum window
323 //
324 Float_t timesignalsum = 0.;
325 Int_t timesquaredsum = 0;
326
327 for (Int_t i=idx; i<idx+fWindowSizeHiGain; i++)
328 {
329 timesignalsum += fHiGainSignal[i]*i;
330 timesquaredsum += i*i;
331 }
332
333 sum = max;
334
335 time = sum != 0 ? timesignalsum / max + Float_t(fHiGainFirst) : 1.;
336 dtime = sum != 0 ? ped.GetPedestalRms() / max * sqrt(timesquaredsum - fWindowSizeHiGain*time) : 1.;
337
338}
339
340
341// --------------------------------------------------------------------------
342//
343// Calculates the arrival time for each pixel
344//
345void MExtractTimeAndChargeSlidingWindow::FindTimeAndChargeLoGain(Byte_t *first, Float_t &sum, Float_t &dsum,
346 Float_t &time, Float_t &dtime,
347 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
348{
349
350 Int_t range = fLoGainLast - fLoGainFirst + 1;
351 const Byte_t *end = first + range;
352 Byte_t *p = first;
353
354 Float_t max = 0; // highest integral content of all windows
355 Int_t count = 0; // counter to recognize the AB-flag
356
357 Float_t pedes = ped.GetPedestal();
358 const Float_t ABoffs = ped.GetPedestalABoffset();
359
360 Float_t PedMean[2];
361 PedMean[0] = pedes + ABoffs;
362 PedMean[1] = pedes - ABoffs;
363
364 //
365 // Check for saturation in all other slices
366 //
367 while (p<first+fWindowSizeLoGain)
368 {
369 const Int_t ids = fLoGainFirst + count;
370 const Float_t signal = (Float_t)*p - PedMean[(ids+abflag) & 0x1];
371 sum += signal;
372 fLoGainSignal[count] = signal;
373
374 if (*p++ >= fSaturationLimit)
375 sat++;
376
377 count++;
378 }
379
380 //
381 // Check for saturation in all other slices
382 //
383 while (p<end)
384 if (*p++ >= fSaturationLimit)
385 sat++;
386
387 if (IsNoiseCalculation())
388 return;
389
390 //
391 // Calculate the i-th sum as
392 // sum_i+1 = sum_i + slice[i+8] - slice[i]
393 // This is fast and accurate (because we are using int's)
394 //
395 count = 0;
396 max = sum;
397 Int_t idx = 0; // idx of the first slice of the maximum window
398
399 for (p=first; p+fWindowSizeLoGain<end; p++)
400 {
401
402 const Int_t ids = fLoGainFirst + count + fWindowSizeLoGain;
403 const Float_t signal = (Float_t)*(p+fWindowSizeLoGain) - PedMean[(ids+abflag) & 0x1];
404 sum += signal - fLoGainSignal[count];
405 fLoGainSignal[count + fWindowSizeLoGain] = signal;
406
407 if (sum>max)
408 {
409 max = sum;
410 idx = count+1;
411 }
412 count++;
413 }
414
415 //
416 // now calculate the time for the maximum window
417 //
418 Float_t timesignalsum = 0;
419 Int_t timesquaredsum = 0;
420
421 for (Int_t i=idx; i<idx+fWindowSizeLoGain; i++)
422 {
423 timesignalsum += fLoGainSignal[i]*i;
424 timesquaredsum += i*i;
425 }
426
427 sum = max;
428
429 time = sum != 0 ? timesignalsum / max + Float_t(fLoGainFirst) : 1.;
430 dtime = sum != 0 ? ped.GetPedestalRms() / max * sqrt(timesquaredsum - fWindowSizeLoGain*time) : 1.;
431}
432
433// --------------------------------------------------------------------------
434//
435// In addition to the resources of the base-class MExtractor:
436// MJPedestal.MExtractor.WindowSizeHiGain: 6
437// MJPedestal.MExtractor.WindowSizeLoGain: 6
438//
439Int_t MExtractTimeAndChargeSlidingWindow::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
440{
441
442 Byte_t hw = fWindowSizeHiGain;
443 Byte_t lw = fWindowSizeLoGain;
444
445 Bool_t rc = kFALSE;
446
447 if (IsEnvDefined(env, prefix, "HiGainWindowSize", print))
448 {
449 hw = GetEnvValue(env, prefix, "HiGainWindowSize", hw);
450 rc = kTRUE;
451 }
452 if (IsEnvDefined(env, prefix, "LoGainWindowSize", print))
453 {
454 lw = GetEnvValue(env, prefix, "LoGainWindowSize", lw);
455 rc = kTRUE;
456 }
457
458 if (rc)
459 SetWindowSize(hw, lw);
460
461 return MExtractTime::ReadEnv(env, prefix, print) ? kTRUE : rc;
462
463}
464
Note: See TracBrowser for help on using the repository browser.