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

Last change on this file since 6107 was 6107, checked in by mazin, 20 years ago
*** empty log message ***
File size: 13.2 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{
75
76 fName = name ? name : "MExtractFixedWindowSpline";
77 fTitle = title ? title : "Signal Extractor for a fixed FADC window using a fast spline";
78
79 SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
80}
81
82// --------------------------------------------------------------------------
83//
84// SetRange:
85//
86// Checks:
87// - if the window defined by (fHiGainLast-fHiGainFirst-1) are odd, subtract one
88// - if the window defined by (fLoGainLast-fLoGainFirst-1) are odd, subtract one
89// - if the Hi Gain window is smaller than 2, set fHiGainLast to fHiGainFirst+1
90// - if the Lo Gain window is smaller than 2, set fLoGainLast to fLoGainFirst+1
91//
92// Calls:
93// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
94//
95// Sets:
96// - fNumHiGainSamples to: (Float_t)(fHiGainLast-fHiGainFirst+1)
97// - fNumLoGainSamples to: (Float_t)(fLoGainLast-fLoGainFirst+1)
98// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
99// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)
100//
101void MExtractFixedWindowSpline::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
102{
103
104 const Byte_t windowhi = hilast-hifirst;
105 const Byte_t whieven = windowhi & ~1;
106
107 if (whieven != windowhi)
108 {
109 *fLog << warn << GetDescriptor()
110 << Form("%s%2i%s%2i",": Hi Gain window size has to be uneven, set last slice from "
111 ,(int)hilast," to ",(int)(hilast-1)) << endl;
112 hilast -= 1;
113 }
114
115 if (whieven<2)
116 {
117 *fLog << warn << GetDescriptor()
118 << Form("%s%2i%s%2i",": Hi Gain window is smaller than 2 FADC sampes, set last slice from"
119 ,(int)hilast," to ",(int)(hifirst+2)) << endl;
120 hilast = hifirst+2;
121 }
122
123 const Byte_t windowlo = lolast-lofirst;
124 const Byte_t wloeven = windowlo & ~1;
125
126 if (lolast != 0)
127 {
128 if (wloeven != windowlo)
129 {
130 *fLog << warn << GetDescriptor()
131 << Form("%s%2i%s%2i",": Lo Gain window size has to be uneven, set last slice from "
132 ,(int)lolast," to ",(int)(lolast-1)) << endl;
133 lolast -= 1;
134 }
135
136
137 if (wloeven<2)
138 {
139 *fLog << warn << GetDescriptor()
140 << Form("%s%2i%s%2i",": Lo Gain window is smaller than 2 FADC sampes, set last slice from"
141 ,(int)lolast," to ",(int)(lofirst+2)) << endl;
142 lolast = lofirst+2;
143 }
144 }
145
146
147 MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
148
149 //
150 // Very important: Because the spline interpolates between the slices,
151 // the number of samples for the pedestal subtraction
152 // is now 1 less than with e.g. MExtractFixedWindow
153 //
154 fNumHiGainSamples = (Float_t)(fHiGainLast-fHiGainFirst);
155 if (fLoGainLast != 0)
156 fNumLoGainSamples = (Float_t)(fLoGainLast-fLoGainFirst);
157 else
158 fNumLoGainSamples = 0.;
159
160 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
161 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
162
163}
164
165// --------------------------------------------------------------------------
166//
167// ReInit
168//
169// Calls:
170// - MExtractor::ReInit(pList);
171// - fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
172// fLoGainFirst, fLoGainLast, fNumLoGainSamples);
173//
174// Deletes all arrays, if not NULL
175// Creates new arrays according to the extraction range
176//
177Bool_t MExtractFixedWindowSpline::ReInit(MParList *pList)
178{
179
180 if (!MExtractor::ReInit(pList))
181 return kFALSE;
182
183 fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
184 fLoGainFirst, fLoGainLast, fNumLoGainSamples);
185
186 Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast;
187
188 fHiGainFirstDeriv.Set(range);
189 fHiGainSecondDeriv.Set(range);
190
191 range = fLoGainLast - fLoGainFirst + 1;
192
193 fLoGainFirstDeriv.Set(range);
194 fLoGainSecondDeriv.Set(range);
195
196 return kTRUE;
197}
198
199
200// --------------------------------------------------------------------------
201//
202// FindSignalHiGain:
203//
204// - Loop from ptr to (ptr+fHiGainLast-fHiGainFirst)
205// - Sum up contents of *ptr
206// - If *ptr is greater than fSaturationLimit, raise sat by 1
207//
208// - If fHiLoLast is not 0, loop also from logain to (logain+fHiLoLast)
209// - Sum up contents of logain
210// - If *logain is greater than fSaturationLimit, raise sat by 1
211//
212void MExtractFixedWindowSpline::FindSignalHiGain(Byte_t *ptr, Byte_t *logain, Float_t &sum, Byte_t &sat) const
213{
214
215 const Byte_t *end = ptr + fHiGainLast - fHiGainFirst;
216 Int_t range = fHiGainLast - fHiGainFirst + fHiLoLast + 1;
217
218 Float_t pp;
219 // Int_t i = 0;
220
221 Int_t summ = 0;
222 //
223 // Take half of the first slice content
224 //
225 Float_t *firstderiv = fHiGainFirstDeriv.GetArray();
226 Float_t *secondderiv = fHiGainSecondDeriv.GetArray();
227 sum = (Float_t)*ptr/2.;
228 //
229 // The first slice has already been treated now!
230 //
231 ptr++; // i++;
232 secondderiv++;
233 firstderiv++;
234 //
235 // Check for saturation in all other slices
236 //
237 while (ptr<end)
238 {
239
240 summ += *ptr;
241
242 // pp = fHiGainSecondDeriv[i-1] + 4.;
243 // fHiGainSecondDeriv[i] = -1.0/pp;
244 // fHiGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
245 // fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
246
247 pp = *(secondderiv-1) + 4.;
248 *secondderiv = -1.0/pp;
249 *firstderiv = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
250 *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
251
252 if (*ptr++ >= fSaturationLimit)
253 sat++;
254
255 secondderiv++;
256 firstderiv++;
257
258 // i++;
259 }
260
261 switch (fHiLoLast)
262 {
263 case 0:
264 // Treat the last slice of the high-gain as half slice:
265 sum += (Float_t)*ptr/2.;
266 break;
267 case 1:
268 // Treat the last slice of the high-gain as full slice:
269 summ += *ptr;
270 pp = *(secondderiv-1) + 4.;
271 *secondderiv = -1.0/pp;
272 *firstderiv = *(logain) - 2.* *(ptr) + *(ptr-1);
273 *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
274 secondderiv++;
275 firstderiv++;
276 if (*logain >= fSaturationLimit)
277 sat++;
278 // Treat the first slice of the low-gain as half slice:
279 sum += (Float_t)*logain/2;
280 break;
281 case 2:
282 // Treat the last slice of the high-gain as full slice:
283 summ += *ptr;
284 pp = *(secondderiv-1) + 4.;
285 *secondderiv = -1.0/pp;
286 *firstderiv = *(logain) - 2.* *(ptr) + *(ptr-1);
287 *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
288 secondderiv++;
289 firstderiv++;
290 // Treat the last first slice of the low-gain as full slice:
291 summ += *logain;
292 pp = *(secondderiv-1) + 4.;
293 *secondderiv = -1.0/pp;
294 *firstderiv = *(logain+1) - 2.* *(logain) + *(ptr);
295 *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
296 secondderiv++;
297 firstderiv++;
298 if (*logain++ >= fSaturationLimit)
299 sat++;
300 // Treat the second slice of the low-gain as half slice:
301 sum += (Float_t)*logain/2;
302 if (*logain >= fSaturationLimit)
303 sat++;
304 break;
305 default:
306 // Treat the last slice of the high-gain as full slice:
307 summ += *ptr;
308 pp = *(secondderiv-1) + 4.;
309 *secondderiv = -1.0/pp;
310 *firstderiv = *(logain) - 2.* *(ptr) + *(ptr-1);
311 *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
312 secondderiv++;
313 firstderiv++;
314 // Treat the last first slice of the low-gain as full slice:
315 summ += *logain;
316 pp = *(secondderiv-1) + 4.;
317 *secondderiv = -1.0/pp;
318 *firstderiv = *(logain+1) - 2.* *(logain) + *(ptr);
319 *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
320 secondderiv++;
321 firstderiv++;
322 if (*logain++ >= fSaturationLimit)
323 sat++;
324 // Treat the rest of the slices:
325 const Byte_t *end = logain+fHiLoLast;
326 while (logain<end)
327 {
328 summ += *logain;
329 pp = *(secondderiv-1) + 4.;
330 *secondderiv = -1.0/pp;
331 *firstderiv = *(logain+1) - 2.* *(logain) + *(logain-1);
332 *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
333 // pp = fHiGainSecondDeriv[i-1] + 4.;
334 // fHiGainSecondDeriv[i] = -1.0/pp;
335 // fHiGainFirstDeriv [i] = *(logain+1) - 2.* *(logain) + *(logain-1);
336 // fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
337 secondderiv++;
338 firstderiv++;
339 if (*logain++ >= fSaturationLimit)
340 sat++;
341 }
342 break;
343 }
344
345 //
346 // Go back to last but one element:
347 //
348 secondderiv--;
349 firstderiv--;
350
351 for (Int_t k=range-2;k>0;k--)
352 {
353 *secondderiv = *secondderiv * *(secondderiv+1) + *firstderiv;
354 sum += 0.25* *secondderiv;
355 secondderiv--;
356 // fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k];
357 // sum += 0.25*fHiGainSecondDeriv[k];
358 }
359
360 sum += (Float_t)summ;
361}
362
363// --------------------------------------------------------------------------
364//
365// FindSignalLoGain:
366//
367// - Loop from ptr to (ptr+fLoGainLast-fLoGainFirst)
368// - Sum up contents of *ptr
369// - If *ptr is greater than fSaturationLimit, raise sat by 1
370//
371void MExtractFixedWindowSpline::FindSignalLoGain(Byte_t *ptr, Float_t &sum, Byte_t &sat) const
372{
373
374 const Byte_t *end = ptr + fLoGainLast - fLoGainFirst;
375 Int_t range = fLoGainLast - fLoGainFirst + 1;
376
377 Float_t pp;
378 // Int_t i = 0;
379
380 Int_t summ = 0;
381 //
382 // Take half of the first slice content
383 //
384 Float_t *firstderiv = fLoGainFirstDeriv.GetArray();
385 Float_t *secondderiv = fLoGainSecondDeriv.GetArray();
386 sum = (Float_t)*ptr/2.;
387 //
388 // The first slice has already been treated now!
389 //
390 ptr++; // i++;
391 secondderiv++;
392 firstderiv++;
393 //
394 // Check for saturation in all other slices
395 //
396 while (ptr<end)
397 {
398
399 summ += *ptr;
400 // i++;
401
402 // pp = fLoGainSecondDeriv[i-1] + 4.;
403 // fLoGainSecondDeriv[i] = -1.0/pp;
404 // fLoGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
405 // fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
406
407 pp = *(secondderiv-1) + 4.;
408 *secondderiv = -1.0/pp;
409 *firstderiv = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
410 *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
411
412 if (*ptr++ >= fSaturationLimit)
413 sat++;
414
415 secondderiv++;
416 firstderiv++;
417 }
418
419 sum += (Float_t)*ptr/2.;
420
421 //
422 // Go back to last but one element:
423 //
424 secondderiv--;
425 firstderiv--;
426
427 for (Int_t k=range-2;k>0;k--)
428 {
429 *secondderiv = *secondderiv * *(secondderiv+1) + *firstderiv;
430 sum += 0.25* *secondderiv;
431 secondderiv--;
432 // fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k];
433 // sum += 0.25*fLoGainSecondDeriv[k];
434 }
435
436 sum += (Float_t)summ;
437}
438
Note: See TracBrowser for help on using the repository browser.