source: trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.cc@ 6034

Last change on this file since 6034 was 5697, checked in by gaug, 20 years ago
*** empty log message ***
File size: 11.8 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 analyzing 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! Author(s): Markus Gaug 05/2004 <mailto:markus@ifae.es>
18!
19! Copyright: MAGIC Software Development, 2002-2004
20!
21!
22\* ======================================================================== */
23
24//////////////////////////////////////////////////////////////////////////////
25//
26// MExtractTimeAndChargeSpline
27//
28// Fast Spline extractor using a cubic spline algorithm of Numerical Recipes.
29// It returns the integral below the interpolating spline.
30//
31// Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast)
32// to modify the ranges.
33//
34// The spline will then be integrated from fHiGainFirst to fHiGainLast,
35// including half of the outer edges. The effective number of intergrated
36// slices ("range") is thus:
37//
38// range = fHiGainLast - fHiGainFirst
39//
40// Ranges have to be an even number. In case of odd ranges, the last slice
41// will be reduced by one.
42//
43// Defaults are:
44//
45// fHiGainFirst = fgHiGainFirst = 2
46// fHiGainLast = fgHiGainLast = 14
47// fLoGainFirst = fgLoGainFirst = 3
48// fLoGainLast = fgLoGainLast = 13
49//
50//////////////////////////////////////////////////////////////////////////////
51#include "MExtractFixedWindowSpline.h"
52
53#include "MExtractedSignalCam.h"
54
55#include "MLog.h"
56#include "MLogManip.h"
57
58ClassImp(MExtractFixedWindowSpline);
59
60using namespace std;
61
62const Byte_t MExtractFixedWindowSpline::fgHiGainFirst = 2;
63const Byte_t MExtractFixedWindowSpline::fgHiGainLast = 14;
64const Byte_t MExtractFixedWindowSpline::fgLoGainFirst = 3;
65const Byte_t MExtractFixedWindowSpline::fgLoGainLast = 13;
66// --------------------------------------------------------------------------
67//
68// Default constructor.
69//
70// Calls:
71// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
72//
73MExtractFixedWindowSpline::MExtractFixedWindowSpline(const char *name, const char *title)
74 : fHiGainFirstDeriv(NULL), fLoGainFirstDeriv(NULL),
75 fHiGainSecondDeriv(NULL), fLoGainSecondDeriv(NULL)
76{
77
78 fName = name ? name : "MExtractFixedWindowSpline";
79 fTitle = title ? title : "Signal Extractor for a fixed FADC window using a fast spline";
80
81 SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
82}
83
84MExtractFixedWindowSpline::~MExtractFixedWindowSpline()
85{
86
87 if (fHiGainFirstDeriv)
88 delete [] fHiGainFirstDeriv;
89 if (fLoGainFirstDeriv)
90 delete [] fLoGainFirstDeriv;
91 if (fHiGainSecondDeriv)
92 delete [] fHiGainSecondDeriv;
93 if (fLoGainSecondDeriv)
94 delete [] fLoGainSecondDeriv;
95
96}
97
98// --------------------------------------------------------------------------
99//
100// SetRange:
101//
102// Checks:
103// - if the window defined by (fHiGainLast-fHiGainFirst-1) are odd, subtract one
104// - if the window defined by (fLoGainLast-fLoGainFirst-1) are odd, subtract one
105// - if the Hi Gain window is smaller than 2, set fHiGainLast to fHiGainFirst+1
106// - if the Lo Gain window is smaller than 2, set fLoGainLast to fLoGainFirst+1
107//
108// Calls:
109// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
110//
111// Sets:
112// - fNumHiGainSamples to: (Float_t)(fHiGainLast-fHiGainFirst+1)
113// - fNumLoGainSamples to: (Float_t)(fLoGainLast-fLoGainFirst+1)
114// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
115// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)
116//
117void MExtractFixedWindowSpline::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
118{
119
120 const Byte_t windowhi = hilast-hifirst;
121 const Byte_t whieven = windowhi & ~1;
122
123 if (whieven != windowhi)
124 {
125 *fLog << warn << GetDescriptor()
126 << Form("%s%2i%s%2i",": Hi Gain window size has to be even, set last slice from "
127 ,(int)hilast," to ",(int)(hilast-1)) << endl;
128 hilast -= 1;
129 }
130
131 if (whieven<2)
132 {
133 *fLog << warn << GetDescriptor()
134 << Form("%s%2i%s%2i",": Hi Gain window is smaller than 2 FADC sampes, set last slice from"
135 ,(int)hilast," to ",(int)(hifirst+1)) << endl;
136 hilast = hifirst+1;
137 }
138
139 const Byte_t windowlo = lolast-lofirst;
140 const Byte_t wloeven = windowlo & ~1;
141
142 if (lolast != 0)
143 {
144 if (wloeven != windowlo)
145 {
146 *fLog << warn << GetDescriptor()
147 << Form("%s%2i%s%2i",": Lo Gain window size has to be even, set last slice from "
148 ,(int)lolast," to ",(int)(lolast-1)) << endl;
149 lolast -= 1;
150 }
151
152
153 if (wloeven<2)
154 {
155 *fLog << warn << GetDescriptor()
156 << Form("%s%2i%s%2i",": Lo Gain window is smaller than 2 FADC sampes, set last slice from"
157 ,(int)lolast," to ",(int)(lofirst+1)) << endl;
158 lolast = lofirst+1;
159 }
160 }
161
162
163 MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
164
165 //
166 // Very important: Because the spline interpolates between the slice,
167 // the number of samples for the pedestal subtraction
168 // is now 1 less than with e.g. MExtractFixedWindow
169 //
170 fNumHiGainSamples = (Float_t)(fHiGainLast-fHiGainFirst);
171 if (fLoGainLast != 0)
172 fNumLoGainSamples = (Float_t)(fLoGainLast-fLoGainFirst);
173 else
174 fNumLoGainSamples = 0.;
175
176 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
177 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
178
179}
180
181// --------------------------------------------------------------------------
182//
183// ReInit
184//
185// Calls:
186// - MExtractor::ReInit(pList);
187// - fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
188// fLoGainFirst, fLoGainLast, fNumLoGainSamples);
189//
190// Deletes all arrays, if not NULL
191// Creates new arrays according to the extraction range
192//
193Bool_t MExtractFixedWindowSpline::ReInit(MParList *pList)
194{
195
196 if (!MExtractor::ReInit(pList))
197 return kFALSE;
198
199 fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
200 fLoGainFirst, fLoGainLast, fNumLoGainSamples);
201
202 if (fHiGainFirstDeriv)
203 delete [] fHiGainFirstDeriv;
204 if (fLoGainFirstDeriv)
205 delete [] fLoGainFirstDeriv;
206 if (fHiGainSecondDeriv)
207 delete [] fHiGainSecondDeriv;
208 if (fLoGainSecondDeriv)
209 delete [] fLoGainSecondDeriv;
210
211 Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast;
212
213 fHiGainFirstDeriv = new Float_t[range];
214 memset(fHiGainFirstDeriv,0,range*sizeof(Float_t));
215 fHiGainSecondDeriv = new Float_t[range];
216 memset(fHiGainSecondDeriv,0,range*sizeof(Float_t));
217
218 range = fLoGainLast - fLoGainFirst + 1;
219
220 fLoGainFirstDeriv = new Float_t[range];
221 memset(fLoGainFirstDeriv,0,range*sizeof(Float_t));
222 fLoGainSecondDeriv = new Float_t[range];
223 memset(fLoGainSecondDeriv,0,range*sizeof(Float_t));
224
225 return kTRUE;
226}
227
228
229// --------------------------------------------------------------------------
230//
231// FindSignalHiGain:
232//
233// - Loop from ptr to (ptr+fHiGainLast-fHiGainFirst)
234// - Sum up contents of *ptr
235// - If *ptr is greater than fSaturationLimit, raise sat by 1
236//
237// - If fHiLoLast is not 0, loop also from logain to (logain+fHiLoLast)
238// - Sum up contents of logain
239// - If *logain is greater than fSaturationLimit, raise sat by 1
240//
241void MExtractFixedWindowSpline::FindSignalHiGain(Byte_t *ptr, Byte_t *logain, Float_t &sum, Byte_t &sat) const
242{
243
244 const Byte_t *end = ptr + fHiGainLast - fHiGainFirst;
245
246 Int_t range = fHiGainLast - fHiGainFirst + fHiLoLast + 1;
247
248 Float_t pp;
249 Int_t i = 0;
250
251 Int_t summ = 0;
252 sum = (Float_t)*ptr++/2.;
253 //
254 // Check for saturation in all other slices
255 //
256 while (ptr<end)
257 {
258
259 summ += *ptr;
260 i++;
261
262 pp = fHiGainSecondDeriv[i-1] + 4.;
263 fHiGainSecondDeriv[i] = -1.0/pp;
264 fHiGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
265 fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
266
267 if (*ptr++ >= fSaturationLimit)
268 sat++;
269 }
270
271 if (*ptr++ >= fSaturationLimit)
272 sat++;
273
274 if (fHiLoLast == 0)
275 {
276 sum += (Float_t)*ptr/2.;
277 fHiGainSecondDeriv[++i] = 0.;
278 }
279 else
280 {
281
282 //
283 // There are two overlap slices which we have to treat sepatately
284 //
285 // First, the last high-gain slice as center
286 //
287 summ += *ptr;
288 i++;
289
290 if (*logain >= fSaturationLimit)
291 sat++;
292
293 if (fHiLoLast == 1)
294 sum += (Float_t)*logain/2;
295 else
296 {
297 pp = fHiGainSecondDeriv[i-1] + 4.;
298 fHiGainSecondDeriv[i] = -1.0/pp;
299 fHiGainFirstDeriv [i] = *(logain) - 2.* *(ptr) + *(ptr-1);
300 fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
301
302 //
303 // Second, the first low-gain slice as center
304 //
305 summ += *logain;
306 i++;
307
308 pp = fHiGainSecondDeriv[i-1] + 4.;
309 fHiGainSecondDeriv[i] = -1.0/pp;
310 fHiGainFirstDeriv [i] = *(logain+1) - 2.* *(logain) + *(ptr);
311 fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
312
313 if (*logain++ >= fSaturationLimit)
314 sat++;
315
316 end = logain + fHiLoLast - 2;
317
318 while (logain<end)
319 {
320
321 summ += *logain;
322 i++;
323
324 pp = fHiGainSecondDeriv[i-1] + 4.;
325 fHiGainSecondDeriv[i] = -1.0/pp;
326 fHiGainFirstDeriv [i] = *(logain+1) - 2.* *(logain) + *(logain-1);
327 fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
328
329 if (*logain++ >= fSaturationLimit)
330 sat++;
331
332 }
333 }
334 }
335
336 fHiGainSecondDeriv[range-1] = 0.;
337
338 for (Int_t k=range-2;k>0;k--)
339 {
340 fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k];
341 sum += 0.25*fHiGainSecondDeriv[k];
342 }
343
344 sum += (Float_t)summ;
345}
346
347// --------------------------------------------------------------------------
348//
349// FindSignalLoGain:
350//
351// - Loop from ptr to (ptr+fLoGainLast-fLoGainFirst)
352// - Sum up contents of *ptr
353// - If *ptr is greater than fSaturationLimit, raise sat by 1
354//
355void MExtractFixedWindowSpline::FindSignalLoGain(Byte_t *ptr, Float_t &sum, Byte_t &sat) const
356{
357
358 const Byte_t *end = ptr + fLoGainLast - fLoGainFirst;
359 Int_t range = fLoGainLast - fLoGainFirst + 1;
360
361 Float_t pp;
362 Int_t i = 0;
363
364 fLoGainSecondDeriv[0] = 0.;
365 fLoGainFirstDeriv[0] = 0.;
366
367 Int_t summ = 0;
368 sum = (Float_t)*ptr++/2.;
369 //
370 // Check for saturation in all other slices
371 //
372 while (ptr<end)
373 {
374
375 summ += *ptr;
376 i++;
377
378 pp = fLoGainSecondDeriv[i-1] + 4.;
379 fLoGainSecondDeriv[i] = -1.0/pp;
380 fLoGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
381 fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
382
383 if (*ptr++ >= fSaturationLimit)
384 sat++;
385 }
386
387 if (*ptr++ >= fSaturationLimit)
388 sat++;
389
390 sum += (Float_t)*ptr/2.;
391 fLoGainSecondDeriv[++i] = 0.;
392
393 for (Int_t k=range-2;k>0;k--)
394 {
395 fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k];
396 sum += 0.25*fLoGainSecondDeriv[k];
397 }
398
399 sum += (Float_t)summ;
400}
401
Note: See TracBrowser for help on using the repository browser.