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

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