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

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