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

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