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

Last change on this file since 6126 was 6117, checked in by gaug, 20 years ago
*** empty log message ***
File size: 11.9 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 = 1;
80const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgHiGainLast = 15;
81const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgLoGainFirst = 3;
82const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgLoGainLast = 14;
83const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgOffsetLeftFromPeak = 1;
84const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgOffsetRightFromPeak = 2;
85const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgPeakSearchWindowSize = 2;
86// --------------------------------------------------------------------------
87//
88// Default constructor.
89//
90// Sets:
91// - fOffsetLeftFromPeak to fgOffsetLeftFromPeak
92// - fOffsetRightFromPeak to fgOffsetRightFromPeak
93//
94MExtractTimeAndChargeDigitalFilterPeakSearch::MExtractTimeAndChargeDigitalFilterPeakSearch(const char *name, const char *title)
95{
96 fName = name ? name : "MExtractTimeAndChargeDigitalFilterPeakSearch";
97 fTitle = title ? title : "Digital Filter";
98
99 SetOffsetLeftFromPeak();
100 SetOffsetRightFromPeak();
101 SetPeakSearchWindowSize();
102}
103
104// --------------------------------------------------------------------------
105//
106// FindPeak
107//
108// Finds highest sum of "window" consecutive FADC slices in a pixel, and store
109// in "startslice" the first slice of the group which yields the maximum sum.
110// In "max" the value of the maximum sum is stored, in "sat" the number of
111// saturated slices.
112//
113void MExtractTimeAndChargeDigitalFilterPeakSearch::FindPeak(Byte_t *ptr, Byte_t &startslice, Int_t &max,
114 Int_t &sat, Byte_t &satpos) const
115{
116
117 const Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1;
118
119 sat = 0;
120 satpos = 0;
121
122 startslice = 0;
123 Int_t sum=0;
124
125 //
126 // Calculate the sum of the first "fPeakSearchWindowSize" slices
127 //
128 sat = 0;
129 Byte_t *p = ptr;
130
131 while (p<ptr+fPeakSearchWindowSize)
132 {
133 sum += *p;
134 if (*p++ >= fSaturationLimit)
135 {
136 if (sat == 0)
137 satpos = p-ptr;
138 sat++;
139 }
140 }
141
142 //
143 // Check for saturation in all other slices
144 //
145 while (p<end)
146 if (*p++ >= fSaturationLimit)
147 {
148 if (sat == 0)
149 satpos = p-ptr;
150 sat++;
151 }
152
153 //
154 // Calculate the i-th sum of n as
155 // sum_i+1 = sum_i + slice[i+fPeakSearchWindowSize] - slice[i]
156 // This is fast and accurate (because we are using int's)
157 //
158 max=sum;
159 for (p=ptr; p+fPeakSearchWindowSize<end; p++)
160 {
161 sum += *(p+fPeakSearchWindowSize) - *p;
162 if (sum > max)
163 {
164 max = sum;
165 startslice = p-ptr+1;
166 }
167 }
168
169 return;
170}
171
172
173// --------------------------------------------------------------------------
174//
175// Process
176//
177// First find the pixel with highest sum of fPeakSearchWindowSize slices (default:2)
178// "startslice" will mark the slice at which the highest sum begins for that pixel.
179//
180// Then define the beginning of the digital filter search window for ALL pixels as the slice
181// before that: startslice-fOffsetLeftFromPeak until: startslice+fOffsetRightFromPeak
182//
183Int_t MExtractTimeAndChargeDigitalFilterPeakSearch::Process()
184{
185
186 MRawEvtPixelIter pixel(fRawEvt);
187
188 Int_t sat;
189 Byte_t satpos;
190 ULong_t gsatpos = 0;
191
192 Int_t maxsumhi = -1000000;
193 Int_t numsat = 0;
194 Byte_t startslice;
195
196 Byte_t hiGainFirstsave = fHiGainFirst;
197 Byte_t hiGainLastsave = fHiGainLast;
198 Byte_t loGainFirstsave = fLoGainFirst;
199 Byte_t loGainLastsave = fLoGainLast;
200
201 Byte_t higainfirst = fHiGainFirst;
202
203 while (pixel.Next())
204 {
205
206 Int_t sumhi;
207 sat = 0;
208
209 FindPeak(pixel.GetHiGainSamples()+fHiGainFirst+fOffsetLeftFromPeak, startslice, sumhi, sat, satpos);
210
211 if (sumhi > maxsumhi && sat == 0)
212 {
213 maxsumhi = sumhi;
214 higainfirst = fHiGainFirst + startslice;
215 }
216 else if (sat)
217 {
218 numsat++;
219 gsatpos += satpos;
220 }
221 }
222
223 //
224 // Check necessary for calibration events
225 //
226 if (numsat > fSignals->GetSize()*0.9)
227 higainfirst = gsatpos/numsat - 1;
228
229 //
230 // Shift the start slice to the left:
231 //
232 if (higainfirst > fOffsetLeftFromPeak)
233 fHiGainFirst = higainfirst - fOffsetLeftFromPeak;
234 else
235 *fLog << warn << " High Gain ranges out of limits to the left!!! " << endl;
236
237 //
238 // Shift the last slice to the right:
239 //
240 if (higainfirst + fOffsetRightFromPeak + fWindowSizeHiGain > hiGainLastsave)
241 fHiGainLast = higainfirst + fOffsetRightFromPeak + fWindowSizeHiGain;
242 else
243 *fLog << warn << " High Gain ranges out of limits to the right!!! " << endl;
244
245 if ( fHiGainFirst+(Int_t)fOffsetLoGain > fLoGainFirst )
246 fLoGainFirst = fHiGainFirst + (Int_t)fOffsetLoGain;
247 else
248 *fLog << warn << " High Gain ranges out of limits to the left!!! " << endl;
249
250 //
251 // Make sure we will not integrate beyond the lo gain limit:
252 //
253 if (fLoGainFirst+fWindowSizeLoGain+fOffsetRightFromPeak <= pixel.GetNumLoGainSamples())
254 fLoGainLast = fLoGainFirst+fWindowSizeLoGain+fOffsetRightFromPeak;
255 else
256 *fLog << warn << " Low Gain ranges out of limits to the right!!! " << endl;
257
258 pixel.Reset();
259
260 while (pixel.Next())
261 {
262 //
263 // Find signal in hi- and lo-gain
264 //
265 Float_t sumhi =0., deltasumhi =0; // Set hi-gain of MExtractedSignalPix valid
266 Float_t timehi=0., deltatimehi=0; // Set hi-gain of MArrivalTimePix valid
267 Byte_t sathi=0;
268
269 // Initialize fMaxBinContent for the case, it gets not set by the derived class
270 fMaxBinContent = fLoGainSwitch + 1;
271
272 const Int_t pixidx = pixel.GetPixelId();
273 const MPedestalPix &ped = (*fPedestals)[pixidx];
274 const Bool_t higainabflag = pixel.HasABFlag();
275
276 FindTimeAndChargeHiGain(pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(),
277 sumhi, deltasumhi,
278 timehi, deltatimehi,
279 sathi, ped, higainabflag);
280
281 //
282 // Make sure that in cases the time couldn't be correctly determined
283 // more meaningfull default values are assigned
284 //
285 if (timehi<0)
286 timehi = -1;
287 if (timehi>pixel.GetNumHiGainSamples())
288 timehi = pixel.GetNumHiGainSamples();
289
290 Float_t sumlo =0., deltasumlo =-1.; // invalidate logain of MExtractedSignalPix
291 Float_t timelo=0., deltatimelo=-1; // invalidate logain of MArrivalTimePix
292 Byte_t satlo=0;
293
294 //
295 // Adapt the low-gain extraction range from the obtained high-gain time
296 //
297 if (pixel.HasLoGain() && (fMaxBinContent > fLoGainSwitch) )
298 {
299 deltasumlo = 0; // make logain of MExtractedSignalPix valid
300 deltatimelo = 0; // make logain of MArrivalTimePix valid
301
302 fLoGainFirstSave = fLoGainFirst;
303 const Byte_t logainstart = sathi
304 ? sathi + (Int_t)fLoGainStartShift
305 : (Byte_t)(timehi + fLoGainStartShift);
306
307 fLoGainFirst = logainstart > fLoGainFirstSave ? logainstart : fLoGainFirstSave;
308
309 if ( fLoGainFirst < fLoGainLast )
310 {
311 const Bool_t logainabflag = (higainabflag + pixel.GetNumHiGainSamples()) & 0x1;
312 FindTimeAndChargeLoGain(pixel.GetLoGainSamples()+fLoGainFirst,
313 sumlo, deltasumlo,
314 timelo, deltatimelo,
315 satlo, ped, logainabflag);
316 }
317 fLoGainFirst = fLoGainFirstSave;
318
319 // Make sure that in cases the time couldn't be correctly determined
320 // more meaningfull default values are assigned
321 if (timelo<0)
322 timelo = -1;
323 if (timelo>pixel.GetNumLoGainSamples())
324 timelo = pixel.GetNumLoGainSamples();
325 }
326
327 MExtractedSignalPix &pix = (*fSignals)[pixidx];
328 MArrivalTimePix &tix = (*fArrTime)[pixidx];
329 pix.SetExtractedSignal(sumhi, deltasumhi,sumlo, deltasumlo);
330 pix.SetGainSaturation(sathi, sathi, satlo);
331
332 tix.SetArrivalTime(timehi, deltatimehi, timelo-fOffsetLoGain, deltatimelo);
333 tix.SetGainSaturation(sathi, sathi, satlo);
334
335 } /* while (pixel.Next()) */
336
337 fArrTime->SetReadyToSave();
338 fSignals->SetReadyToSave();
339
340 fHiGainFirst = hiGainFirstsave;
341 fHiGainLast = hiGainLastsave ;
342 fLoGainFirst = loGainFirstsave;
343 fLoGainLast = loGainLastsave ;
344
345 return kTRUE;
346}
347
348
349
350// --------------------------------------------------------------------------
351//
352// Read the setup from a TEnv, eg:
353// MJPedestal.MExtractor.WindowSizeHiGain: 6
354// MJPedestal.MExtractor.WindowSizeLoGain: 6
355// MJPedestal.MExtractor.BinningResolutionHiGain: 10
356// MJPedestal.MExtractor.BinningResolutionLoGain: 10
357// MJPedestal.MExtractor.WeightsFile: filename
358//
359Int_t MExtractTimeAndChargeDigitalFilterPeakSearch::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
360{
361
362 Bool_t rc = kFALSE;
363
364 if (IsEnvDefined(env, prefix, "OffsetLeftFromPeak", print))
365 {
366 fOffsetLeftFromPeak = GetEnvValue(env, prefix, fOffsetLeftFromPeak);
367 rc = kTRUE;
368 }
369
370 if (IsEnvDefined(env, prefix, "OffsetRightFromPeak", print))
371 {
372 fOffsetRightFromPeak = GetEnvValue(env, prefix, fOffsetRightFromPeak);
373 rc = kTRUE;
374 }
375
376 if (IsEnvDefined(env, prefix, "PeakSearchWindowSize", print))
377 {
378 fPeakSearchWindowSize = GetEnvValue(env, prefix, fPeakSearchWindowSize);
379 rc = kTRUE;
380 }
381
382 return MExtractTimeAndChargeDigitalFilter::ReadEnv(env, prefix, print) ? kTRUE : rc;
383}
384
385
386void MExtractTimeAndChargeDigitalFilterPeakSearch::Print(Option_t *o) const
387{
388 if (IsA()==Class())
389 *fLog << GetDescriptor() << ":" << endl;
390
391 MExtractTimeAndChargeDigitalFilter::Print(o);
392 *fLog << " Offset from Peak left: " << fOffsetLeftFromPeak << endl;
393 *fLog << " Offset from Peak right: " << fOffsetRightFromPeak << endl;
394 *fLog << " Peak search window size: " << fPeakSearchWindowSize << endl;
395}
Note: See TracBrowser for help on using the repository browser.