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

Last change on this file since 6893 was 6840, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 13.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): 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 = 6;
78// --------------------------------------------------------------------------
79//
80// Default constructor.
81//
82// Calls:
83// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
84//
85// Initializes:
86// - fWindowSizeHiGain to fgHiGainWindowSize
87// - fWindowSizeLoGain to fgLoGainWindowSize
88//
89MExtractTimeAndChargeSlidingWindow::MExtractTimeAndChargeSlidingWindow(const char *name, const char *title)
90{
91
92 fName = name ? name : "MExtractTimeAndChargeSlidingWindow";
93 fTitle = title ? title : "Calculate arrival times and charges using a sliding window";
94
95 fWindowSizeHiGain = fgHiGainWindowSize;
96 fWindowSizeLoGain = fgLoGainWindowSize;
97
98 SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
99}
100
101//-------------------------------------------------------------------
102//
103// Set the ranges
104//
105// Calls:
106// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
107// - SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
108//
109void MExtractTimeAndChargeSlidingWindow::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
110{
111
112 MExtractor::SetRange(hifirst, hilast, lofirst, lolast);
113
114 //
115 // Redo the checks if the window is still inside the ranges
116 //
117 SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
118
119}
120
121// -----------------------------------------------------------------------------------------
122//
123// Checks:
124// - if a window is bigger than the one defined by the ranges, set it to the available range
125// - if a window is smaller than 2, set it to 2
126//
127// Sets:
128// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain
129// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain
130// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
131// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)
132//
133void MExtractTimeAndChargeSlidingWindow::SetWindowSize(Int_t windowh, Int_t windowl)
134{
135
136 fWindowSizeHiGain = windowh;
137 fWindowSizeLoGain = windowl;
138
139 const Int_t availhirange = (Int_t)(fHiGainLast-fHiGainFirst+1);
140
141 if (fWindowSizeHiGain > availhirange)
142 {
143 *fLog << warn << GetDescriptor()
144 << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",(int)fWindowSizeHiGain,
145 " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl;
146 *fLog << warn << GetDescriptor()
147 << ": Will set window size to: " << (int)availhirange << endl;
148 fWindowSizeHiGain = availhirange;
149 }
150
151 if (fWindowSizeHiGain<1)
152 {
153 fWindowSizeHiGain = 1;
154 *fLog << warn << GetDescriptor() << ": High Gain window size too small, set to one sample" << endl;
155 }
156
157 if (fLoGainLast != 0 && fWindowSizeLoGain != 0)
158 {
159 const Int_t availlorange = (Int_t)(fLoGainLast-fLoGainFirst+1);
160
161 if (fWindowSizeLoGain > availlorange)
162 {
163 *fLog << warn << GetDescriptor()
164 << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",(int)fWindowSizeLoGain,
165 " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl;
166 *fLog << warn << GetDescriptor()
167 << ": Will set window size to: " << (int)availlorange << endl;
168 fWindowSizeLoGain = availlorange;
169 }
170 }
171
172 fNumHiGainSamples = (Float_t)fWindowSizeHiGain;
173 fNumLoGainSamples = fLoGainLast ? (Float_t)fWindowSizeLoGain : 0.;
174
175 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
176 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
177
178}
179
180// --------------------------------------------------------------------------
181//
182// InitArrays
183//
184// Gets called in the ReInit() and initialized the arrays
185//
186Bool_t MExtractTimeAndChargeSlidingWindow::InitArrays()
187{
188 Int_t range = (Int_t)(fHiGainLast - fHiGainFirst + 1 + fHiLoLast);
189 fHiGainSignal.Set(range);
190 range = (Int_t)(fLoGainLast - fLoGainFirst + 1);
191 fLoGainSignal.Set(range);
192
193 return kTRUE;
194
195}
196
197// --------------------------------------------------------------------------
198//
199// Calculates the arrival time for each pixel
200//
201void MExtractTimeAndChargeSlidingWindow::FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum,
202 Float_t &time, Float_t &dtime,
203 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
204{
205
206 Int_t range = fHiGainLast - fHiGainFirst + 1;
207 const Byte_t *end = first + range;
208 Byte_t *p = first;
209
210 Float_t max = 0; // highest integral content of all windows
211 sat = 0;
212
213 const Float_t pedes = ped.GetPedestal();
214 const Float_t ABoffs = ped.GetPedestalABoffset();
215
216 const Float_t PedMean[2] = { pedes + ABoffs, pedes - ABoffs };
217
218 fMaxBinContent = 0;
219 //
220 // Check for saturation in all other slices
221 //
222 Int_t ids = fHiGainFirst;
223 Float_t *sample = fHiGainSignal.GetArray();
224
225 while (p<first+fWindowSizeHiGain)
226 {
227
228 const Float_t signal = (Float_t)*p - PedMean[(ids++ + abflag) & 0x1];
229 sum += signal;
230 *sample++ = signal;
231
232 if (*p > fMaxBinContent)
233 fMaxBinContent = *p;
234
235 if (*p++ >= fSaturationLimit)
236 if (!sat)
237 sat = ids-2;
238 }
239
240 if (IsNoiseCalculation())
241 return;
242
243 //
244 // Check for saturation in all other slices
245 //
246 while (p<end)
247 {
248 if (*p > fMaxBinContent)
249 fMaxBinContent = *p;
250
251 if (*p++ >= fSaturationLimit)
252 if (!sat)
253 sat = ids-2;
254 }
255
256 //
257 // Calculate the i-th sum as
258 // sum_i+1 = sum_i + slice[i+8] - slice[i]
259 // This is fast and accurate (because we are using int's)
260 //
261 Int_t count = 0;
262 max = sum;
263 Int_t idx = 0; // idx of the first slice of the maximum window
264
265 for (p=first; p+fWindowSizeHiGain<end; p++)
266 {
267
268 const Float_t signal = (Float_t)*(p+fWindowSizeHiGain) - PedMean[(ids++ + abflag) & 0x1];
269 sum += signal - *(sample-fWindowSizeHiGain);
270 *sample++ = signal;
271
272 if (sum>max)
273 {
274 max = sum;
275 idx = count+1;
276 }
277 count++;
278 }
279
280 //
281 // overlap bins
282 //
283 Byte_t *l = logain;
284 while (l < logain+fHiLoLast)
285 {
286
287 const Float_t signal = (Float_t)*l - PedMean[(ids++ + abflag) & 0x1];
288 sum += signal - *(sample-fWindowSizeHiGain);
289 *sample++ = signal;
290
291 if (*l > fMaxBinContent)
292 fMaxBinContent = *logain;
293
294 if (*l++ >= fSaturationLimit)
295 if (!sat)
296 sat = ids-2;
297
298 if (sum>max)
299 {
300 max = sum;
301 idx = count+1;
302 }
303 count++;
304 } /* while (l<logain+fHiLoLast) */
305
306 //
307 // now calculate the time for the maximum window
308 //
309 Float_t timesignalsum = 0.;
310 Int_t timesquaredsum = 0;
311
312 for (Int_t i=idx; i<idx+fWindowSizeHiGain; i++)
313 {
314 timesignalsum += fHiGainSignal[i]*i;
315 timesquaredsum += i*i;
316 }
317
318 sum = max;
319
320 time = max > 0.1 ? timesignalsum / max + Float_t(fHiGainFirst) : 1.;
321 dtime = max > 0.1 ? ped.GetPedestalRms() / max * sqrt(timesquaredsum - fWindowSizeHiGain*time) : 1.;
322
323}
324
325
326// --------------------------------------------------------------------------
327//
328// Calculates the arrival time for each pixel
329//
330void MExtractTimeAndChargeSlidingWindow::FindTimeAndChargeLoGain(Byte_t *first, Float_t &sum, Float_t &dsum,
331 Float_t &time, Float_t &dtime,
332 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
333{
334
335 Int_t range = fLoGainLast - fLoGainFirst + 1;
336 const Byte_t *end = first + range;
337 Byte_t *p = first;
338
339 Float_t max = 0; // highest integral content of all windows
340 Int_t count = 0; // counter to recognize the AB-flag
341
342 Float_t pedes = ped.GetPedestal();
343 const Float_t ABoffs = ped.GetPedestalABoffset();
344
345 Float_t PedMean[2] = { pedes + ABoffs, pedes - ABoffs };
346 //
347 // Check for saturation in all other slices
348 //
349 Int_t ids = fLoGainFirst;
350
351 while (p<first+fWindowSizeLoGain)
352 {
353 const Float_t signal = (Float_t)*p - PedMean[(ids++ + abflag) & 0x1];
354 sum += signal;
355 fLoGainSignal[count] = signal;
356
357 if (*p++ >= fSaturationLimit)
358 sat++;
359
360 count++;
361 }
362
363 //
364 // Check for saturation in all other slices
365 //
366 while (p<end)
367 if (*p++ >= fSaturationLimit)
368 sat++;
369
370 if (IsNoiseCalculation())
371 return;
372
373 //
374 // Calculate the i-th sum as
375 // sum_i+1 = sum_i + slice[i+8] - slice[i]
376 // This is fast and accurate (because we are using int's)
377 //
378 count = 0;
379 max = sum;
380 Int_t idx = 0; // idx of the first slice of the maximum window
381
382 for (p=first; p+fWindowSizeLoGain<end; p++)
383 {
384
385 const Float_t signal = (Float_t)*(p+fWindowSizeLoGain) - PedMean[(ids++ + abflag) & 0x1];
386 sum += signal - fLoGainSignal[count];
387 fLoGainSignal[count + fWindowSizeLoGain] = signal;
388
389 if (sum>max)
390 {
391 max = sum;
392 idx = count+1;
393 }
394 count++;
395 }
396
397 //
398 // now calculate the time for the maximum window
399 //
400 Float_t timesignalsum = 0;
401 Int_t timesquaredsum = 0;
402
403 for (Int_t i=idx; i<idx+fWindowSizeLoGain; i++)
404 {
405 timesignalsum += fLoGainSignal[i]*i;
406 timesquaredsum += i*i;
407 }
408
409 sum = max;
410
411 time = max > 0.1 ? timesignalsum / max + Float_t(fLoGainFirst) : 1.;
412 dtime = max > 0.1 ? ped.GetPedestalRms() / max * sqrt(timesquaredsum - fWindowSizeLoGain*time) : 1.;
413}
414
415// --------------------------------------------------------------------------
416//
417// In addition to the resources of the base-class MExtractor:
418// MJPedestal.MExtractor.WindowSizeHiGain: 6
419// MJPedestal.MExtractor.WindowSizeLoGain: 6
420//
421Int_t MExtractTimeAndChargeSlidingWindow::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
422{
423
424 Byte_t hw = fWindowSizeHiGain;
425 Byte_t lw = fWindowSizeLoGain;
426
427 Bool_t rc = kFALSE;
428
429 if (IsEnvDefined(env, prefix, "HiGainWindowSize", print))
430 {
431 hw = GetEnvValue(env, prefix, "HiGainWindowSize", hw);
432 rc = kTRUE;
433 }
434 if (IsEnvDefined(env, prefix, "LoGainWindowSize", print))
435 {
436 lw = GetEnvValue(env, prefix, "LoGainWindowSize", lw);
437 rc = kTRUE;
438 }
439
440 if (rc)
441 SetWindowSize(hw, lw);
442
443 return MExtractTimeAndCharge::ReadEnv(env, prefix, print) ? kTRUE : rc;
444
445}
446
Note: See TracBrowser for help on using the repository browser.