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

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