source: trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilterPeakSearch.cc@ 6742

Last change on this file since 6742 was 6332, checked in by gaug, 20 years ago
*** empty log message ***
File size: 17.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 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): Markus Gaug, 02/2005 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24//////////////////////////////////////////////////////////////////////////////
25//
26// MExtractTimeAndChargeDigitalFilterPeakSearch
27//
28// An extractor using the digital filter in combination with a global
29// peak search, analogue to the class MExtractFixedWindowPeakSearch.
30//
31// The extractor returns kFALSE if too many events lie outside the
32// range: (peak - fOffsetLeftFromPeak, peak - fOffsetRightFromPeak)
33//
34// Input Containers:
35// MRawEvtData
36// MRawRunHeader
37// MPedestalCam
38//
39// Output Containers:
40// MArrivalTimeCam
41// MExtractedSignalCam
42//
43//////////////////////////////////////////////////////////////////////////////
44#include "MExtractTimeAndChargeDigitalFilterPeakSearch.h"
45
46#include <fstream>
47
48#include <TFile.h>
49#include <TH1F.h>
50#include <TH2F.h>
51#include <TString.h>
52#include <TMatrix.h>
53
54#include "MLog.h"
55#include "MLogManip.h"
56
57#include "MParList.h"
58
59#include "MPedestalPix.h"
60#include "MPedestalCam.h"
61
62#include "MRawEvtPixelIter.h"
63
64#include "MExtractedSignalCam.h"
65#include "MExtractedSignalPix.h"
66
67#include "MArrivalTimeCam.h"
68#include "MArrivalTimePix.h"
69
70#include "MBadPixelsCam.h"
71#include "MBadPixelsPix.h"
72
73ClassImp(MExtractTimeAndChargeDigitalFilterPeakSearch);
74
75using namespace std;
76
77const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgHiGainFirst = 0;
78const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgHiGainLast = 20;
79const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgLoGainFirst = 0;
80const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgLoGainLast = 14;
81const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgOffsetLeftFromPeak = 3;
82const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgOffsetRightFromPeak = 3;
83const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgPeakSearchWindowSize = 2;
84const Int_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgHiGainFailureLimit = 10;
85const Int_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgLoGainFailureLimit = 25;
86// --------------------------------------------------------------------------
87//
88// Default constructor.
89//
90// Sets:
91// - fOffsetLeftFromPeak to fgOffsetLeftFromPeak
92// - fOffsetRightFromPeak to fgOffsetRightFromPeak
93// - fPeakSearchWindowSize to fgPeakSearchWindowSize
94// - fHiGainFailureLimit to fgHiGainFailureLimit
95// - fLoGainFailureLimit to fgLoGainFailureLimit
96//
97MExtractTimeAndChargeDigitalFilterPeakSearch::MExtractTimeAndChargeDigitalFilterPeakSearch(const char *name, const char *title)
98 : fHiGainOutOfRangeLeft(0), fHiGainOutOfRangeRight(0),
99 fLoGainOutOfRangeLeft(0), fLoGainOutOfRangeRight(0),
100 fBadPixels(NULL)
101{
102 fName = name ? name : "MExtractTimeAndChargeDigitalFilterPeakSearch";
103 fTitle = title ? title : "Digital Filter";
104
105 SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
106
107 SetOffsetLeftFromPeak();
108 SetOffsetRightFromPeak();
109 SetPeakSearchWindowSize();
110 SetHiGainFailureLimit();
111 SetLoGainFailureLimit();
112}
113
114// --------------------------------------------------------------------------
115//
116// The PreProcess searches for the following input containers:
117// - MRawEvtData
118// - MRawRunHeader
119// - MPedestalCam
120//
121// The following output containers are also searched and created if
122// they were not found:
123//
124// - MExtractedSignalCam
125// - MArrivalTimeCam
126//
127// The following variables are set to 0:
128//
129// - fHiGainOutOfRangeLeft
130// - fHiGainOutOfRangeRight
131// - fLoGainOutOfRangeLeft
132// - fLoGainOutOfRangeRight
133//
134Int_t MExtractTimeAndChargeDigitalFilterPeakSearch::PreProcess(MParList *pList)
135{
136
137 if (!MExtractTimeAndCharge::PreProcess(pList))
138 return kFALSE;
139
140 fHiGainOutOfRangeLeft = 0;
141 fHiGainOutOfRangeRight = 0;
142 fLoGainOutOfRangeLeft = 0;
143 fLoGainOutOfRangeRight = 0;
144
145 fBadPixels = (MBadPixelsCam*)pList->FindObject(AddSerialNumber("MBadPixelsCam"));
146 if (!fBadPixels)
147 {
148 *fLog << err << "Cannot find MBadPixelsCam ... abort." << endl;
149 return kFALSE;
150 }
151
152 return kTRUE;
153}
154
155// --------------------------------------------------------------------------
156//
157// FindPeak
158//
159// Finds highest sum of "window" consecutive FADC slices in a pixel, and store
160// in "startslice" the first slice of the group which yields the maximum sum.
161// In "max" the value of the maximum sum is stored, in "sat" the number of
162// saturated slices.
163//
164void MExtractTimeAndChargeDigitalFilterPeakSearch::FindPeak(Byte_t *ptr, Byte_t *logain, Byte_t &startslice, Int_t &max,
165 Int_t &sat, Byte_t &satpos) const
166{
167
168 Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1;
169
170 sat = 0;
171 satpos = 0;
172
173 startslice = 0;
174 Int_t sum=0;
175
176 //
177 // Calculate the sum of the first "fPeakSearchWindowSize" slices
178 //
179 sat = 0;
180 Byte_t *p = ptr;
181
182 while (p<ptr+fPeakSearchWindowSize)
183 {
184 sum += *p;
185 if (*p++ >= fSaturationLimit)
186 {
187 if (sat == 0)
188 satpos = p-ptr;
189 sat++;
190 }
191 }
192
193 //
194 // Check for saturation in all other slices
195 //
196 while (p<end)
197 if (*p++ >= fSaturationLimit)
198 {
199 if (sat == 0)
200 satpos = p-ptr;
201 sat++;
202 }
203
204 //
205 // Calculate the i-th sum of n as
206 // sum_i+1 = sum_i + slice[i+fPeakSearchWindowSize] - slice[i]
207 // This is fast and accurate (because we are using int's)
208 //
209 max=sum;
210 for (p=ptr; p+fPeakSearchWindowSize<end; p++)
211 {
212 sum += *(p+fPeakSearchWindowSize) - *p;
213 if (sum > max)
214 {
215 max = sum;
216 startslice = p-ptr+1;
217 }
218 }
219
220 if (!fHiLoLast)
221 return;
222
223 Byte_t *l = logain;
224
225 while (++p < end && l < logain+fHiLoLast)
226 {
227 sum += *l - *p;
228 if (sum > max)
229 {
230 max = sum;
231 startslice = p-ptr+1;
232 }
233
234 if (*l++ >= fSaturationLimit)
235 {
236 if (sat == 0)
237 satpos = l-logain + fHiGainLast - fHiGainFirst;
238 sat++;
239 }
240 }
241
242 if (fHiLoLast <= fPeakSearchWindowSize)
243 return;
244
245 while (l < logain+fHiLoLast)
246 {
247 if (*l++ >= fSaturationLimit)
248 {
249 if (sat == 0)
250 satpos = l-logain+fHiGainLast-fHiGainFirst;
251 sat++;
252 }
253
254 sum += *l - *(l-fPeakSearchWindowSize);
255 if (sum > max)
256 {
257 max = sum;
258 startslice = l-logain + fHiGainLast - fHiGainFirst - fPeakSearchWindowSize + 1;
259 }
260 }
261
262 return;
263}
264
265
266// --------------------------------------------------------------------------
267//
268// Process
269//
270// First find the pixel with highest sum of fPeakSearchWindowSize slices (default:2)
271// "startslice" will mark the slice at which the highest sum begins for that pixel.
272//
273// Then define the beginning of the digital filter search window for ALL pixels as the slice
274// before that: startslice-fOffsetLeftFromPeak until: startslice+fOffsetRightFromPeak
275//
276Int_t MExtractTimeAndChargeDigitalFilterPeakSearch::Process()
277{
278
279 MRawEvtPixelIter pixel(fRawEvt);
280
281 Int_t sat = 0;
282 Byte_t satpos = 0;
283 ULong_t gsatpos = 0;
284
285 Int_t maxsumhi = -1000000;
286 Int_t numsat = 0;
287 Byte_t startslice = 0;
288
289 Byte_t hiGainFirstsave = fHiGainFirst;
290 Byte_t hiGainLastsave = fHiGainLast;
291 Byte_t loGainFirstsave = fLoGainFirst;
292 Byte_t loGainLastsave = fLoGainLast;
293 Byte_t hiLoLastsave = fHiLoLast;
294
295 Byte_t higainfirst = fHiGainFirst;
296 UInt_t peakpixel = 0;
297
298 while (pixel.Next())
299 {
300
301 Int_t sumhi;
302 sat = 0;
303
304 MBadPixelsPix &bad = (*fBadPixels)[pixel.GetPixelId()];
305
306 if (!bad.IsUnsuitable())
307 FindPeak(pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(), startslice, sumhi, sat, satpos);
308
309 if (sumhi > maxsumhi && sat == 0)
310 {
311 maxsumhi = sumhi;
312 higainfirst = fHiGainFirst + startslice;
313 peakpixel = pixel.GetPixelId();
314 }
315 else if (sat)
316 {
317 numsat++;
318 gsatpos += satpos;
319 }
320 }
321
322 //
323 // Check necessary for calibration events
324 //
325 if (numsat > fSignals->GetSize()*0.9)
326 higainfirst = gsatpos/numsat - 1;
327
328 //
329 // Shift the start slice to the left:
330 //
331 if (higainfirst >= fHiGainFirst + fOffsetLeftFromPeak)
332 fHiGainFirst = higainfirst - fOffsetLeftFromPeak;
333 else
334 {
335 // *fLog << err << (Int_t)higainfirst << " " << peakpixel << endl;
336 fHiGainOutOfRangeLeft++;
337 }
338
339 //
340 // Shift the last slice to the right:
341 //
342 const Byte_t rlim = higainfirst + fOffsetRightFromPeak + fWindowSizeHiGain - 1;
343 if (rlim <= fHiGainLast+fHiLoLast)
344 if (rlim > fHiGainLast)
345 {
346 fHiLoLast = rlim - fHiGainLast;
347 fHiGainLast = pixel.GetNumHiGainSamples() - 1;
348 }
349 else
350 {
351 fHiLoLast = 0;
352 fHiGainLast = rlim;
353 }
354 else
355 {
356 fHiGainOutOfRangeRight++;
357 // *fLog << err << (Int_t)higainfirst << " " << peakpixel << endl;
358 }
359
360 const Byte_t llim = higainfirst + (Int_t)fOffsetLoGain;
361 if ( llim >= fLoGainFirst + fOffsetLeftFromPeak)
362 fLoGainFirst = llim - fOffsetLeftFromPeak;
363 else
364 fLoGainOutOfRangeLeft++;
365
366 //
367 // Make sure we will not integrate beyond the lo gain limit:
368 //
369 const Byte_t lolast = fLoGainFirst+fWindowSizeLoGain+fOffsetRightFromPeak-1;
370 if (lolast < pixel.GetNumLoGainSamples())
371 fLoGainLast = lolast;
372 else
373 {
374 fLoGainOutOfRangeRight++;
375 // *fLog << err << (Int_t)higainfirst << " " << peakpixel << " " << (int)fLoGainFirst << endl;
376 }
377
378 // *fLog << inf << (int)fHiGainFirst << " " << (int)higainfirst << " " << (int)fHiGainLast
379 // << " " << (int)fLoGainFirst << " " << (int)fLoGainLast << endl;
380
381 pixel.Reset();
382
383 while (pixel.Next())
384 {
385
386 //
387 // Find signal in hi- and lo-gain
388 //
389 Float_t sumhi =0., deltasumhi =0; // Set hi-gain of MExtractedSignalPix valid
390 Float_t timehi=0., deltatimehi=0; // Set hi-gain of MArrivalTimePix valid
391 Byte_t sathi=0;
392
393 // Initialize fMaxBinContent for the case, it gets not set by the derived class
394 fMaxBinContent = fLoGainSwitch + 1;
395
396 const Int_t pixidx = pixel.GetPixelId();
397 const MPedestalPix &ped = (*fPedestals)[pixidx];
398 const Bool_t higainabflag = pixel.HasABFlag();
399
400 FindTimeAndChargeHiGain(pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(),
401 sumhi, deltasumhi,
402 timehi, deltatimehi,
403 sathi, ped, higainabflag);
404
405 //
406 // Make sure that in cases the time couldn't be correctly determined
407 // more meaningfull default values are assigned
408 //
409 if (timehi<0)
410 timehi = -1;
411 if (timehi>pixel.GetNumHiGainSamples())
412 timehi = pixel.GetNumHiGainSamples();
413
414 Float_t sumlo =0., deltasumlo =-1.; // invalidate logain of MExtractedSignalPix
415 Float_t timelo=0., deltatimelo=-1; // invalidate logain of MArrivalTimePix
416 Byte_t satlo=0;
417
418 //
419 // Adapt the low-gain extraction range from the obtained high-gain time
420 //
421 if (pixel.HasLoGain() && (fMaxBinContent > fLoGainSwitch) )
422 {
423 deltasumlo = 0; // make logain of MExtractedSignalPix valid
424 deltatimelo = 0; // make logain of MArrivalTimePix valid
425
426 fLoGainFirstSave = fLoGainFirst;
427 const Byte_t logainstart = sathi
428 ? sathi + (Int_t)fLoGainStartShift
429 : (Byte_t)(timehi + fLoGainStartShift);
430
431 fLoGainFirst = logainstart > fLoGainFirstSave ? logainstart : fLoGainFirstSave;
432
433 if ( fLoGainFirst < fLoGainLast )
434 {
435 const Bool_t logainabflag = (higainabflag + pixel.GetNumHiGainSamples()) & 0x1;
436 FindTimeAndChargeLoGain(pixel.GetLoGainSamples()+fLoGainFirst,
437 sumlo, deltasumlo,
438 timelo, deltatimelo,
439 satlo, ped, logainabflag);
440 }
441 fLoGainFirst = fLoGainFirstSave;
442
443 // Make sure that in cases the time couldn't be correctly determined
444 // more meaningfull default values are assigned
445 if (timelo<0)
446 timelo = -1;
447 if (timelo>pixel.GetNumLoGainSamples())
448 timelo = pixel.GetNumLoGainSamples();
449 }
450
451 MExtractedSignalPix &pix = (*fSignals)[pixidx];
452 MArrivalTimePix &tix = (*fArrTime)[pixidx];
453 pix.SetExtractedSignal(sumhi, deltasumhi,sumlo, deltasumlo);
454 pix.SetGainSaturation(sathi, sathi, satlo);
455
456 tix.SetArrivalTime(timehi, deltatimehi, timelo-fOffsetLoGain, deltatimelo);
457 tix.SetGainSaturation(sathi, sathi, satlo);
458
459 } /* while (pixel.Next()) */
460
461 fArrTime->SetReadyToSave();
462 fSignals->SetReadyToSave();
463
464 fHiGainFirst = hiGainFirstsave;
465 fHiGainLast = hiGainLastsave ;
466 fLoGainFirst = loGainFirstsave;
467 fLoGainLast = loGainLastsave ;
468 fHiLoLast = hiLoLastsave;
469
470 return kTRUE;
471}
472
473Int_t MExtractTimeAndChargeDigitalFilterPeakSearch::PostProcess()
474{
475
476 if (GetNumExecutions() < 1)
477 return kTRUE;
478
479 *fLog << warn << endl;
480
481 const Int_t higainfailure = ((fHiGainOutOfRangeLeft+fHiGainOutOfRangeRight)*100)/GetNumExecutions();
482 const Int_t logainfailure = ((fLoGainOutOfRangeLeft+fLoGainOutOfRangeRight)*100)/GetNumExecutions();
483
484 if (fHiGainOutOfRangeLeft > 0)
485 *fLog << Form("%5.1f",((Float_t)fHiGainOutOfRangeLeft*100)/GetNumExecutions())
486 << "% ranging out of hi-gain window to the left!" << endl;
487 if (fHiGainOutOfRangeRight > 0)
488 *fLog << Form("%5.1f",((Float_t)fHiGainOutOfRangeRight*100)/GetNumExecutions())
489 << "% ranging out of hi-gain window to the right!" << endl;
490 if (fLoGainOutOfRangeLeft > 0)
491 *fLog << Form("%5.1f",((Float_t)fLoGainOutOfRangeLeft*100)/GetNumExecutions())
492 << "% ranging out of lo-gain window to the left!" << endl;
493 if (fLoGainOutOfRangeRight > 0)
494 *fLog << Form("%5.1f",((Float_t)fLoGainOutOfRangeRight*100)/GetNumExecutions())
495 << "% ranging out of lo-gain window to the right!" << endl;
496
497
498 if (higainfailure > fHiGainFailureLimit)
499 {
500 *fLog << err << higainfailure << "% range failures in high gain above limit of: " << fHiGainFailureLimit << "%." << endl;
501 return kFALSE;
502 }
503
504 if (logainfailure > fLoGainFailureLimit)
505 {
506 *fLog << err << logainfailure << "% range failures in low gain above limit of: " << fLoGainFailureLimit << "%." << endl;
507 return kFALSE;
508 }
509
510 return kTRUE;
511
512}
513
514// --------------------------------------------------------------------------
515//
516// Read the setup from a TEnv, eg:
517// MJPedestal.MExtractor.WindowSizeHiGain: 6
518// MJPedestal.MExtractor.WindowSizeLoGain: 6
519// MJPedestal.MExtractor.BinningResolutionHiGain: 10
520// MJPedestal.MExtractor.BinningResolutionLoGain: 10
521// MJPedestal.MExtractor.WeightsFile: filename
522//
523Int_t MExtractTimeAndChargeDigitalFilterPeakSearch::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
524{
525
526 Bool_t rc = kFALSE;
527
528 if (IsEnvDefined(env, prefix, "OffsetLeftFromPeak", print))
529 {
530 fOffsetLeftFromPeak = GetEnvValue(env, prefix, "OffsetLeftFromPeak", fOffsetLeftFromPeak);
531 rc = kTRUE;
532 }
533
534 if (IsEnvDefined(env, prefix, "OffsetRightFromPeak", print))
535 {
536 fOffsetRightFromPeak = GetEnvValue(env, prefix, "OffsetRightFromPeak", fOffsetRightFromPeak);
537 rc = kTRUE;
538 }
539
540 if (IsEnvDefined(env, prefix, "PeakSearchWindowSize", print))
541 {
542 fPeakSearchWindowSize = GetEnvValue(env, prefix, "PeakSearchWindowSize", fPeakSearchWindowSize);
543 rc = kTRUE;
544 }
545
546 if (IsEnvDefined(env, prefix, "HiGainFailureLimit", print))
547 {
548 fHiGainFailureLimit = GetEnvValue(env, prefix, "HiGainFailureLimit", fHiGainFailureLimit);
549 rc = kTRUE;
550 }
551
552 if (IsEnvDefined(env, prefix, "LoGainFailureLimit", print))
553 {
554 fLoGainFailureLimit = GetEnvValue(env, prefix, "LoGainFailureLimit", fLoGainFailureLimit);
555 rc = kTRUE;
556 }
557
558 return MExtractTimeAndChargeDigitalFilter::ReadEnv(env, prefix, print) ? kTRUE : rc;
559}
560
561
562void MExtractTimeAndChargeDigitalFilterPeakSearch::Print(Option_t *o) const
563{
564 if (IsA()==Class())
565 *fLog << GetDescriptor() << ":" << endl;
566
567 MExtractTimeAndChargeDigitalFilter::Print(o);
568 *fLog << " Offset from Peak left: " << (int)fOffsetLeftFromPeak << endl;
569 *fLog << " Offset from Peak right: " << (int)fOffsetRightFromPeak << endl;
570 *fLog << " Peak search window size: " << (int)fPeakSearchWindowSize << endl;
571 *fLog << " High Gain Failure limit: " << fHiGainFailureLimit << endl;
572 *fLog << " Low Gain Failure limit: " << fLoGainFailureLimit << endl;
573}
Note: See TracBrowser for help on using the repository browser.