source: branches/Corsika7500Compatibility/msignal/MExtractFixedWindowPeakSearch.cc@ 18558

Last change on this file since 18558 was 7810, checked in by tbretz, 18 years ago
*** empty log message ***
File size: 17.7 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 analysing 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!
18! Author(s): Abelardo Moralejo, 04/2004 <mailto:moralejo@pd.infn.it>
19! Author(s): Markus Gaug, 04/2004 <mailto:markus@ifae.es>
20!
21! Copyright: MAGIC Software Development, 2000-2004
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27//
28// MExtractFixedWindowPeakSearch
29//
30// Calculate the signal integrating fWindowSize time slices, the same for
31// all pixels. The integrated slices change from event to event, since the
32// pulse positions in the FADC jump between events, but apparently in a
33// "coherent" fashion. We first loop over all pixels and find the one
34// which has the highest sum of fPeakSearchWindowSize (default: 4) consecutive
35// non-saturated high gain slices. The position of the peak in this pixel
36// determines the integration window position for all pixels of this event.
37// For the moment we neglect time delays between pixels (which are in most
38// cases small). The class is based on MExtractSlidingWindow.
39//
40// Call: SetRange(higainfirst, higainlast, logainfirst, logainlast)
41// to modify the ranges in which the window is allowed to move.
42// Defaults are:
43//
44// fHiGainFirst = fgHiGainFirst = 0
45// fHiGainLast = fgHiGainLast = 14
46// fLoGainFirst = fgLoGainFirst = 3
47// fLoGainLast = fgLoGainLast = 14
48//
49// Call: SetWindows(windowhigain, windowlogain,peaksearchwindow)
50// to modify the sliding window widths. Windows have to be an even number.
51// In case of odd numbers, the window will be modified.
52//
53// Call: SetOffsetFromWindow() to adjust the positionning of the extraction
54// window w.r.t. the peak search window. An fOffsetFromWindow of 0 means that
55// the starting slice of the extraction window is equal to the starting slice
56// of the maximizing peak search window. fOffsetFromWindow equal to 1 (default)
57// means that the extraction window starts one slice earlier than the peak
58// search window result. It is recommanded to always use a smaller peak search
59// window than the extraction window.
60//
61// Defaults are:
62//
63// fHiGainWindowSize = fgHiGainWindowSize = 6
64// fLoGainWindowSize = fgLoGainWindowSize = 6
65// fPeakSearchWindowSize = fgPeakSearchWindowSize = 4
66// fLoGainPeakShift = fgLoGainPeakShift = 0
67// fOffsetFromWindow = fgOffsetFromWindow = 1
68//
69//////////////////////////////////////////////////////////////////////////////
70#include "MExtractFixedWindowPeakSearch.h"
71#include "MExtractor.h"
72
73#include "MLog.h"
74#include "MLogManip.h"
75
76#include "MRawEvtPixelIter.h"
77
78#include "MExtractedSignalCam.h"
79#include "MExtractedSignalPix.h"
80
81#include "MPedestalCam.h"
82#include "MPedestalPix.h"
83
84ClassImp(MExtractFixedWindowPeakSearch);
85
86using namespace std;
87
88const Byte_t MExtractFixedWindowPeakSearch::fgHiGainFirst = 0;
89const Byte_t MExtractFixedWindowPeakSearch::fgHiGainLast = 16;
90const Byte_t MExtractFixedWindowPeakSearch::fgLoGainFirst = 3;
91const Byte_t MExtractFixedWindowPeakSearch::fgLoGainLast = 14;
92const Byte_t MExtractFixedWindowPeakSearch::fgHiGainWindowSize = 6;
93const Byte_t MExtractFixedWindowPeakSearch::fgLoGainWindowSize = 6;
94const Byte_t MExtractFixedWindowPeakSearch::fgPeakSearchWindowSize = 4;
95const Byte_t MExtractFixedWindowPeakSearch::fgOffsetFromWindow = 1;
96const Byte_t MExtractFixedWindowPeakSearch::fgLoGainPeakShift = 1;
97// --------------------------------------------------------------------------
98//
99// Default constructor.
100//
101// Sets:
102// - fHiGainWindowSize to fgHiGainWindowSize
103// - fLoGainWindowSize to fgLoGainWindowSize
104// - fPeakSearchWindowSize to fgPeakSearchWindowSize
105// - fLoGainPeakShift to fgLoGainPeakShift
106//
107// Calls:
108// - SetOffsetFromWindow()
109// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
110//
111MExtractFixedWindowPeakSearch::MExtractFixedWindowPeakSearch(const char *name, const char *title)
112 : fHiGainWindowSize(fgHiGainWindowSize),
113 fLoGainWindowSize(fgLoGainWindowSize),
114 fPeakSearchWindowSize(fgPeakSearchWindowSize),
115 fLoGainPeakShift(fgLoGainPeakShift)
116{
117
118 fName = name ? name : "MExtractFixedWindowPeakSearch";
119 fTitle = title ? title : "Task to extract the signal from the FADC slices";
120
121 SetOffsetFromWindow();
122 SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
123}
124
125// --------------------------------------------------------------------------
126//
127// SetRange:
128//
129// Calls:
130// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
131// - SetWindows(fHiGainWindowSize,fLoGainWindowSize,fPeakSearchWindowSize);
132//
133void MExtractFixedWindowPeakSearch::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
134{
135
136 MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
137
138 //
139 // Redo the checks if the window is still inside the ranges
140 //
141 SetWindows(fHiGainWindowSize,fLoGainWindowSize,fPeakSearchWindowSize);
142
143}
144
145// --------------------------------------------------------------------------
146//
147// Checks:
148// - if a window is odd, subtract one
149// - if a window is bigger than the one defined by the ranges, set it to the available range
150// - if a window is smaller than 2, set it to 2
151//
152// Sets:
153// - fNumHiGainSamples to: (Float_t)fHiGainWindowSize
154// - fNumLoGainSamples to: (Float_t)fLoGainWindowSize
155// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
156// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)
157//
158void MExtractFixedWindowPeakSearch::SetWindows(Byte_t windowh, Byte_t windowl, Byte_t peaksearchwindow)
159{
160
161 fHiGainWindowSize = windowh & ~1;
162 fLoGainWindowSize = windowl & ~1;
163 fPeakSearchWindowSize = peaksearchwindow & ~1;
164
165 if (fHiGainWindowSize != windowh)
166 *fLog << warn << GetDescriptor() << ": Hi Gain window size has to be even, set to: "
167 << int(fHiGainWindowSize) << " samples " << endl;
168
169 if (fLoGainWindowSize != windowl)
170 *fLog << warn << GetDescriptor() << ": Lo Gain window size has to be even, set to: "
171 << int(fLoGainWindowSize) << " samples " << endl;
172
173 if (fPeakSearchWindowSize != peaksearchwindow)
174 *fLog << warn << GetDescriptor() << ": Peak Search window size has to be even, set to: "
175 << int(fPeakSearchWindowSize) << " samples " << endl;
176
177 const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1;
178 const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1;
179
180 if (fHiGainWindowSize > availhirange)
181 {
182 *fLog << warn << GetDescriptor()
183 << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",(int)fHiGainWindowSize,
184 " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl;
185 *fLog << warn << GetDescriptor()
186 << ": Will set window size to: " << (int)availhirange << endl;
187 fHiGainWindowSize = availhirange;
188 }
189
190
191 if (fLoGainWindowSize > availlorange)
192 {
193 *fLog << warn << GetDescriptor()
194 << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",(int)fLoGainWindowSize,
195 " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl;
196 *fLog << warn << GetDescriptor()
197 << ": Will set window size to: " << (int)availlorange << endl;
198 fLoGainWindowSize = availlorange;
199 }
200
201 if (fHiGainWindowSize<2)
202 {
203 fHiGainWindowSize = 2;
204 *fLog << warn << GetDescriptor() << ": Hi Gain window size set to two samples" << endl;
205 }
206
207 if (fLoGainWindowSize<2)
208 {
209 fLoGainWindowSize = 2;
210 *fLog << warn << GetDescriptor() << ": Lo Gain window size set to two samples" << endl;
211 }
212
213 if (fPeakSearchWindowSize<2)
214 {
215 fPeakSearchWindowSize = 2;
216 *fLog << warn << GetDescriptor()
217 << ": Peak Search window size set to two samples" << endl;
218 }
219
220 fNumHiGainSamples = (Float_t)fHiGainWindowSize;
221 fNumLoGainSamples = (Float_t)fLoGainWindowSize;
222
223 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
224 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
225}
226
227
228// --------------------------------------------------------------------------
229//
230// The ReInit calls:
231// - MExtractor::ReInit()
232// - fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
233// fLoGainFirst, fLoGainLast, fNumLoGainSamples);
234//
235Bool_t MExtractFixedWindowPeakSearch::ReInit(MParList *pList)
236{
237
238 MExtractor::ReInit(pList);
239
240 fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
241 fLoGainFirst, fLoGainLast, fNumLoGainSamples);
242
243 *fLog << dec << endl;
244 *fLog << inf << "Taking " << fNumHiGainSamples
245 << " HiGain samples starting with slice " << (Int_t)fHiGainFirst
246 << " to " << (Int_t)(fHiGainLast+fHiLoLast) << " incl" << endl;
247 *fLog << inf << "Taking " << fNumLoGainSamples
248 << " LoGain samples starting with slice " << (Int_t)fLoGainFirst
249 << " to " << (Int_t)fLoGainLast << " incl" << endl;
250
251 return kTRUE;
252}
253
254// --------------------------------------------------------------------------
255//
256// FindPeak
257// Finds highest sum of "window" consecutive FADC slices in a pixel, and store
258// in "startslice" the first slice of the group which yields the maximum sum.
259// In "max" the value of the maximum sum is stored, in "sat" the number of
260// saturated slices.
261//
262void MExtractFixedWindowPeakSearch::FindPeak(Byte_t *ptr, Byte_t window, Byte_t &startslice, Int_t &max,
263 Int_t &sat, Byte_t &satpos) const
264{
265
266 const Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1;
267
268 sat = 0;
269 satpos = 0;
270
271 startslice = 0;
272 Int_t sum=0;
273
274 //
275 // Calculate the sum of the first "window" slices
276 //
277 sat = 0;
278 Byte_t *p = ptr;
279
280 while (p<ptr+window)
281 {
282 sum += *p;
283 if (*p++ >= fSaturationLimit)
284 {
285 if (sat == 0)
286 satpos = p-ptr;
287 sat++;
288 }
289 }
290
291 //
292 // Check for saturation in all other slices
293 //
294 while (p<end)
295 if (*p++ >= fSaturationLimit)
296 {
297 if (sat == 0)
298 satpos = p-ptr;
299 sat++;
300 }
301
302 //
303 // Calculate the i-th sum of n as
304 // sum_i+1 = sum_i + slice[i+window] - slice[i]
305 // This is fast and accurate (because we are using int's)
306 //
307 max=sum;
308 for (p=ptr; p+window<end; p++)
309 {
310 sum += *(p+window) - *p;
311 if (sum > max)
312 {
313 max = sum;
314 startslice = p-ptr+1;
315 }
316 }
317
318 return;
319}
320
321// --------------------------------------------------------------------------
322//
323// FindSignalHiGain:
324//
325// - Loop from ptr to (ptr+fHiGainWindowSize)
326// - Sum up contents of *ptr
327// - If *ptr is greater than fSaturationLimit, raise sat by 1
328//
329void MExtractFixedWindowPeakSearch::FindSignalHiGain(Byte_t *ptr, Byte_t *logain, Float_t &sum, Byte_t &sat) const
330{
331
332 Byte_t *end = ptr + fHiGainWindowSize-fHiLoLast;
333
334 Int_t summ = 0;
335 //
336 // Calculate the sum of the "window" slices starting in ptr
337 //
338 while (ptr<end)
339 {
340 summ += *ptr;
341 if (*ptr++ >= fSaturationLimit)
342 sat++;
343 }
344
345 //
346 // If part of the "lo-Gain" slices are used,
347 // repeat steps one and two for the logain part until fHiLoLast
348 //
349 Byte_t *p = logain;
350 end = logain + fHiLoLast;
351 while (p<end)
352 {
353
354 summ += *p;
355 if (*p++ >= fSaturationLimit)
356 sat++;
357 }
358
359 sum = (Float_t)summ;
360 return;
361}
362
363// --------------------------------------------------------------------------
364//
365// FindSignalLoGain:
366//
367// - Loop from ptr to (ptr+fLoGainWindowSize)
368// - Sum up contents of *ptr
369// - If *ptr is greater than fSaturationLimit, raise sat by 1
370//
371void MExtractFixedWindowPeakSearch::FindSignalLoGain(Byte_t *ptr, Float_t &sum, Byte_t &sat) const
372{
373 //
374 // Calculate the sum of the "window" slices starting in ptr
375 //
376 Byte_t *p = ptr;
377 Int_t summ = 0;
378
379 while (p<ptr+fLoGainWindowSize)
380 {
381 summ += *p;
382 if (*p++ >= fSaturationLimit)
383 sat++;
384 }
385
386 sum = (Float_t)summ;
387 return;
388}
389
390// --------------------------------------------------------------------------
391//
392// Process
393// First find the pixel with highest sum of fPeakSearchWindowSize slices (default:4)
394// "startslice" will mark the slice at which the highest sum begins for that pixel.
395// Then define the beginning of the integration window for ALL pixels as the slice
396// before that: startslice-fOffsetFromWindow, unless of course startslice-fOffsetFromWindow<=0,
397// in which case we start at 0. We will also check that the integration window does not
398// go beyond the FADC limits.
399//
400Int_t MExtractFixedWindowPeakSearch::Process()
401{
402
403 MRawEvtPixelIter pixel(fRawEvt);
404
405 Int_t sat;
406 Byte_t satpos;
407 ULong_t gsatpos = 0;
408
409 Int_t maxsumhi = -1000000;
410 Int_t numsat = 0;
411 Byte_t startslice;
412 Byte_t hiGainFirst = 0;
413 Byte_t loGainFirst = 0;
414 fHiLoLast = 0;
415
416 while (pixel.Next())
417 {
418 Int_t sumhi;
419 sat = 0;
420
421 FindPeak(pixel.GetHiGainSamples()+fHiGainFirst, fPeakSearchWindowSize, startslice, sumhi, sat, satpos);
422
423 if (sumhi > maxsumhi && sat == 0)
424 {
425 maxsumhi = sumhi;
426 if (startslice > fOffsetFromWindow)
427 hiGainFirst = fHiGainFirst + startslice - fOffsetFromWindow;
428 else
429 hiGainFirst = fHiGainFirst;
430 }
431 else if (sat)
432 {
433 numsat++;
434 gsatpos += satpos;
435 }
436 }
437
438 // Check necessary for calibration events
439 if (numsat > fSignals->GetSize()*0.9)
440 hiGainFirst = gsatpos/numsat - 1;
441
442 loGainFirst = ( hiGainFirst+fLoGainPeakShift > fLoGainFirst )
443 ? hiGainFirst+fLoGainPeakShift
444 : fLoGainFirst;
445
446 // Make sure we will not integrate beyond the hi gain limit:
447 if (hiGainFirst+fHiGainWindowSize > pixel.GetNumHiGainSamples())
448 fHiLoLast = hiGainFirst+fHiGainWindowSize - pixel.GetNumHiGainSamples();
449 // hiGainFirst = pixel.GetNumHiGainSamples()-fHiGainWindowSize;
450
451 // Make sure we will not integrate beyond the lo gain limit:
452 if (loGainFirst+fLoGainWindowSize > pixel.GetNumLoGainSamples())
453 loGainFirst = pixel.GetNumLoGainSamples()-fLoGainWindowSize;
454
455 pixel.Reset();
456
457 while (pixel.Next())
458 {
459 //
460 // Find signal in hi- and lo-gain
461 //
462 Float_t sumhi=0.;
463 Byte_t sathi=0;
464
465 FindSignalHiGain(pixel.GetHiGainSamples()+hiGainFirst, pixel.GetLoGainSamples(), sumhi, sathi);
466
467 Float_t sumlo=0.;
468 Byte_t satlo=0;
469 if (pixel.HasLoGain())
470 FindSignalLoGain(pixel.GetLoGainSamples()+loGainFirst, sumlo, satlo);
471
472 //
473 // Take corresponding pedestal
474 //
475 const Int_t pixid = pixel.GetPixelId();
476
477 const MPedestalPix &ped = (*fPedestals)[pixid];
478 MExtractedSignalPix &pix = (*fSignals)[pixid];
479
480 const Float_t pedes = ped.GetPedestal();
481 const Float_t pedrms = ped.GetPedestalRms();
482 //
483 // Set extracted signal with pedestal substracted
484 //
485 pix.SetExtractedSignal(sumhi - pedes*fNumHiGainSamples, pedrms*fSqrtHiGainSamples,
486 sumlo - pedes*fNumLoGainSamples, pedrms*fSqrtLoGainSamples);
487
488 pix.SetGainSaturation(sathi, satlo);
489
490 // pix.SetNumHiGainSlices(fNumHiGainSamples);
491 // pix.SetNumLoGainSlices(fNumLoGainSamples);
492
493 } /* while (pixel.Next()) */
494
495
496 fSignals->SetReadyToSave();
497
498 return kTRUE;
499}
500
501// --------------------------------------------------------------------------
502//
503// In addition to the resources of the base-class MExtractor:
504// MJPedestal.MExtractor.WindowSizeHiGain: 6
505// MJPedestal.MExtractor.WindowSizeLoGain: 6
506// MJPedestal.MExtractor.PeakSearchWindow: 4
507// MJPedestal.MExtractor.OffsetFromWindow: 1
508// MJPedestal.MExtractor.LoGainPeakShift: 1
509//
510Int_t MExtractFixedWindowPeakSearch::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
511{
512 Byte_t hw = fHiGainWindowSize;
513 Byte_t lw = fLoGainWindowSize;
514 Byte_t pw = fPeakSearchWindowSize;
515
516 Bool_t rc = kFALSE;
517
518 if (IsEnvDefined(env, prefix, "PeakSearchWindow", print))
519 {
520 pw = GetEnvValue(env, prefix, "PeakSearchWindow", pw);
521 rc = kTRUE;
522 }
523 if (IsEnvDefined(env, prefix, "HiGainWindowSize", print))
524 {
525 hw = GetEnvValue(env, prefix, "HiGainWindowSize", hw);
526 rc = kTRUE;
527 }
528 if (IsEnvDefined(env, prefix, "LoGainWindowSize", print))
529 {
530 lw = GetEnvValue(env, prefix, "LoGainWindowSize", lw);
531 rc = kTRUE;
532 }
533
534 if (rc)
535 SetWindows(hw, lw, pw);
536
537
538 if (IsEnvDefined(env, prefix, "OffsetFromWindow", print))
539 {
540 SetOffsetFromWindow(GetEnvValue(env, prefix, "OffsetFromWindow", fOffsetFromWindow));
541 rc = kTRUE;
542 }
543
544 if (IsEnvDefined(env, prefix, "LoGainPeakShift", print))
545 {
546 SetLoGainPeakShift(GetEnvValue(env, prefix, "LoGainPeakShift", fLoGainPeakShift));
547 rc = kTRUE;
548 }
549
550 rc = MExtractor::ReadEnv(env, prefix, print) ? kTRUE : rc;
551
552 return rc;
553}
554
555void MExtractFixedWindowPeakSearch::Print(Option_t *o) const
556{
557 *fLog << all;
558 *fLog << GetDescriptor() << ":" << endl;
559 *fLog << " Windows: Hi-Gain=" << (int)fHiGainWindowSize << " Lo-Gain=" << (int)fLoGainWindowSize;
560 *fLog << " Peak-Search=" << (int)fPeakSearchWindowSize << endl;
561 *fLog << " Offset From Window: " << (int)fOffsetFromWindow << endl;
562 *fLog << " Lo-Gain peak Shift: " << (int)fLoGainPeakShift << endl;
563 MExtractor::Print(o);
564}
Note: See TracBrowser for help on using the repository browser.