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

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