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

Last change on this file since 6290 was 6290, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 16.3 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 = 0;
274 Byte_t satpos = 0;
275 ULong_t gsatpos = 0;
276
277 Int_t maxsumhi = -1000000;
278 Int_t numsat = 0;
279 Byte_t startslice = 0;
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 {
322 // *fLog << err << (Int_t)higainfirst << " " << (int)fHiGainFirst << " " << (int)fOffsetLeftFromPeak << endl;
323 fHiGainOutOfRangeLeft++;
324 }
325
326 //
327 // Shift the last slice to the right:
328 //
329 const Byte_t rlim = higainfirst + fOffsetRightFromPeak + fWindowSizeHiGain;
330 if (rlim <= fHiGainLast+fHiLoLast)
331 if (rlim > fHiGainLast)
332 fHiLoLast = rlim - fHiGainLast;
333 else
334 {
335 fHiLoLast = 0;
336 fHiGainLast = rlim;
337 }
338 else
339 {
340 fHiGainOutOfRangeRight++;
341 // *fLog << err << (Int_t)higainfirst << endl;
342 }
343
344 const Byte_t llim = fHiGainFirst + (Int_t)fOffsetLoGain;
345 if ( llim >= fLoGainFirst )
346 fLoGainFirst = llim;
347 else
348 fLoGainOutOfRangeLeft++;
349
350 //
351 // Make sure we will not integrate beyond the lo gain limit:
352 //
353 if (fLoGainFirst+fWindowSizeLoGain+fOffsetRightFromPeak <= pixel.GetNumLoGainSamples())
354 fLoGainLast = fLoGainFirst+fWindowSizeLoGain+fOffsetRightFromPeak;
355 else
356 fLoGainOutOfRangeRight++;
357
358 pixel.Reset();
359
360 while (pixel.Next())
361 {
362
363 //
364 // Find signal in hi- and lo-gain
365 //
366 Float_t sumhi =0., deltasumhi =0; // Set hi-gain of MExtractedSignalPix valid
367 Float_t timehi=0., deltatimehi=0; // Set hi-gain of MArrivalTimePix valid
368 Byte_t sathi=0;
369
370 // Initialize fMaxBinContent for the case, it gets not set by the derived class
371 fMaxBinContent = fLoGainSwitch + 1;
372
373 const Int_t pixidx = pixel.GetPixelId();
374 const MPedestalPix &ped = (*fPedestals)[pixidx];
375 const Bool_t higainabflag = pixel.HasABFlag();
376
377 FindTimeAndChargeHiGain(pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(),
378 sumhi, deltasumhi,
379 timehi, deltatimehi,
380 sathi, ped, higainabflag);
381
382 //
383 // Make sure that in cases the time couldn't be correctly determined
384 // more meaningfull default values are assigned
385 //
386 if (timehi<0)
387 timehi = -1;
388 if (timehi>pixel.GetNumHiGainSamples())
389 timehi = pixel.GetNumHiGainSamples();
390
391 Float_t sumlo =0., deltasumlo =-1.; // invalidate logain of MExtractedSignalPix
392 Float_t timelo=0., deltatimelo=-1; // invalidate logain of MArrivalTimePix
393 Byte_t satlo=0;
394
395 //
396 // Adapt the low-gain extraction range from the obtained high-gain time
397 //
398 if (pixel.HasLoGain() && (fMaxBinContent > fLoGainSwitch) )
399 {
400 deltasumlo = 0; // make logain of MExtractedSignalPix valid
401 deltatimelo = 0; // make logain of MArrivalTimePix valid
402
403 fLoGainFirstSave = fLoGainFirst;
404 const Byte_t logainstart = sathi
405 ? sathi + (Int_t)fLoGainStartShift
406 : (Byte_t)(timehi + fLoGainStartShift);
407
408 fLoGainFirst = logainstart > fLoGainFirstSave ? logainstart : fLoGainFirstSave;
409
410 if ( fLoGainFirst < fLoGainLast )
411 {
412 const Bool_t logainabflag = (higainabflag + pixel.GetNumHiGainSamples()) & 0x1;
413 FindTimeAndChargeLoGain(pixel.GetLoGainSamples()+fLoGainFirst,
414 sumlo, deltasumlo,
415 timelo, deltatimelo,
416 satlo, ped, logainabflag);
417 }
418 fLoGainFirst = fLoGainFirstSave;
419
420 // Make sure that in cases the time couldn't be correctly determined
421 // more meaningfull default values are assigned
422 if (timelo<0)
423 timelo = -1;
424 if (timelo>pixel.GetNumLoGainSamples())
425 timelo = pixel.GetNumLoGainSamples();
426 }
427
428 MExtractedSignalPix &pix = (*fSignals)[pixidx];
429 MArrivalTimePix &tix = (*fArrTime)[pixidx];
430 pix.SetExtractedSignal(sumhi, deltasumhi,sumlo, deltasumlo);
431 pix.SetGainSaturation(sathi, sathi, satlo);
432
433 tix.SetArrivalTime(timehi, deltatimehi, timelo-fOffsetLoGain, deltatimelo);
434 tix.SetGainSaturation(sathi, sathi, satlo);
435
436 } /* while (pixel.Next()) */
437
438 fArrTime->SetReadyToSave();
439 fSignals->SetReadyToSave();
440
441 fHiGainFirst = hiGainFirstsave;
442 fHiGainLast = hiGainLastsave ;
443 fLoGainFirst = loGainFirstsave;
444 fLoGainLast = loGainLastsave ;
445 fHiLoLast = hiLoLastsave;
446
447 return kTRUE;
448}
449
450Int_t MExtractTimeAndChargeDigitalFilterPeakSearch::PostProcess()
451{
452 if (GetNumExecutions() == 0)
453 return kTRUE;
454
455 *fLog << warn << endl;
456
457 const Int_t higainfailure = ((fHiGainOutOfRangeLeft+fHiGainOutOfRangeRight)*100)/GetNumExecutions();
458 const Int_t logainfailure = ((fLoGainOutOfRangeLeft+fLoGainOutOfRangeRight)*100)/GetNumExecutions();
459
460 if (fHiGainOutOfRangeLeft > 0)
461 *fLog << Form("%5.1f",((Float_t)fHiGainOutOfRangeLeft*100)/GetNumExecutions())
462 << "% ranging out of hi-gain window to the left!" << endl;
463 if (fHiGainOutOfRangeRight > 0)
464 *fLog << Form("%5.1f",((Float_t)fHiGainOutOfRangeRight*100)/GetNumExecutions())
465 << "% ranging out of hi-gain window to the right!" << endl;
466 if (fLoGainOutOfRangeLeft > 0)
467 *fLog << Form("%5.1f",((Float_t)fLoGainOutOfRangeLeft*100)/GetNumExecutions())
468 << "% ranging out of lo-gain window to the left!" << endl;
469 if (fHiGainOutOfRangeRight > 0)
470 *fLog << Form("%5.1f",((Float_t)fLoGainOutOfRangeRight*100)/GetNumExecutions())
471 << "% ranging out of lo-gain window to the right!" << endl;
472
473
474 if (higainfailure > fHiGainFailureLimit)
475 {
476 *fLog << err << higainfailure << "% range failures in high gain above limit of: " << fHiGainFailureLimit << "%." << endl;
477 return kFALSE;
478 }
479
480 if (logainfailure > fLoGainFailureLimit)
481 {
482 *fLog << err << logainfailure << "% range failures in low gain above limit of: " << fLoGainFailureLimit << "%." << endl;
483 return kFALSE;
484 }
485
486 return kTRUE;
487
488}
489
490// --------------------------------------------------------------------------
491//
492// Read the setup from a TEnv, eg:
493// MJPedestal.MExtractor.WindowSizeHiGain: 6
494// MJPedestal.MExtractor.WindowSizeLoGain: 6
495// MJPedestal.MExtractor.BinningResolutionHiGain: 10
496// MJPedestal.MExtractor.BinningResolutionLoGain: 10
497// MJPedestal.MExtractor.WeightsFile: filename
498//
499Int_t MExtractTimeAndChargeDigitalFilterPeakSearch::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
500{
501
502 Bool_t rc = kFALSE;
503
504 if (IsEnvDefined(env, prefix, "OffsetLeftFromPeak", print))
505 {
506 fOffsetLeftFromPeak = GetEnvValue(env, prefix, "OffsetLeftFromPeak", fOffsetLeftFromPeak);
507 rc = kTRUE;
508 }
509
510 if (IsEnvDefined(env, prefix, "OffsetRightFromPeak", print))
511 {
512 fOffsetRightFromPeak = GetEnvValue(env, prefix, "OffsetRightFromPeak", fOffsetRightFromPeak);
513 rc = kTRUE;
514 }
515
516 if (IsEnvDefined(env, prefix, "PeakSearchWindowSize", print))
517 {
518 fPeakSearchWindowSize = GetEnvValue(env, prefix, "PeakSearchWindowSize", fPeakSearchWindowSize);
519 rc = kTRUE;
520 }
521
522 if (IsEnvDefined(env, prefix, "HiGainFailureLimit", print))
523 {
524 fHiGainFailureLimit = GetEnvValue(env, prefix, "HiGainFailureLimit", fHiGainFailureLimit);
525 rc = kTRUE;
526 }
527
528 if (IsEnvDefined(env, prefix, "LoGainFailureLimit", print))
529 {
530 fLoGainFailureLimit = GetEnvValue(env, prefix, "LoGainFailureLimit", fLoGainFailureLimit);
531 rc = kTRUE;
532 }
533
534 return MExtractTimeAndChargeDigitalFilter::ReadEnv(env, prefix, print) ? kTRUE : rc;
535}
536
537
538void MExtractTimeAndChargeDigitalFilterPeakSearch::Print(Option_t *o) const
539{
540 if (IsA()==Class())
541 *fLog << GetDescriptor() << ":" << endl;
542
543 MExtractTimeAndChargeDigitalFilter::Print(o);
544 *fLog << " Offset from Peak left: " << (int)fOffsetLeftFromPeak << endl;
545 *fLog << " Offset from Peak right: " << (int)fOffsetRightFromPeak << endl;
546 *fLog << " Peak search window size: " << (int)fPeakSearchWindowSize << endl;
547 *fLog << " High Gain Failure limit: " << fHiGainFailureLimit << endl;
548 *fLog << " Low Gain Failure limit: " << fLoGainFailureLimit << endl;
549}
Note: See TracBrowser for help on using the repository browser.