source: branches/MarsMoreSimulationTruth/msignal/MExtractFixedWindowSpline.cc@ 18656

Last change on this file since 18656 was 6188, checked in by gaug, 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 firstderiv++;
233 secondderiv++;
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 firstderiv--;
356 secondderiv--;
357 // fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k];
358 // sum += 0.25*fHiGainSecondDeriv[k];
359 }
360
361 sum += (Float_t)summ;
362}
363
364// --------------------------------------------------------------------------
365//
366// FindSignalLoGain:
367//
368// - Loop from ptr to (ptr+fLoGainLast-fLoGainFirst)
369// - Sum up contents of *ptr
370// - If *ptr is greater than fSaturationLimit, raise sat by 1
371//
372void MExtractFixedWindowSpline::FindSignalLoGain(Byte_t *ptr, Float_t &sum, Byte_t &sat) const
373{
374
375 const Byte_t *end = ptr + fLoGainLast - fLoGainFirst;
376 Int_t range = fLoGainLast - fLoGainFirst + 1;
377
378 Float_t pp;
379 // Int_t i = 0;
380
381 Int_t summ = 0;
382 //
383 // Take half of the first slice content
384 //
385 Float_t *firstderiv = fLoGainFirstDeriv.GetArray();
386 Float_t *secondderiv = fLoGainSecondDeriv.GetArray();
387 sum = (Float_t)*ptr/2.;
388 //
389 // The first slice has already been treated now!
390 //
391 ptr++; // i++;
392 secondderiv++;
393 firstderiv++;
394 //
395 // Check for saturation in all other slices
396 //
397 while (ptr<end)
398 {
399
400 summ += *ptr;
401 // i++;
402
403 // pp = fLoGainSecondDeriv[i-1] + 4.;
404 // fLoGainSecondDeriv[i] = -1.0/pp;
405 // fLoGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
406 // fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
407
408 pp = *(secondderiv-1) + 4.;
409 *secondderiv = -1.0/pp;
410 *firstderiv = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
411 *firstderiv = (6.0* *(firstderiv) - *(firstderiv-1))/pp;
412
413 if (*ptr++ >= fSaturationLimit)
414 sat++;
415
416 secondderiv++;
417 firstderiv++;
418 }
419
420 sum += (Float_t)*ptr/2.;
421
422 //
423 // Go back to last but one element:
424 //
425 secondderiv--;
426 firstderiv--;
427
428 for (Int_t k=range-2;k>0;k--)
429 {
430 *secondderiv = *secondderiv * *(secondderiv+1) + *firstderiv;
431 sum += 0.25* *secondderiv;
432 firstderiv--;
433 secondderiv--;
434 // fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k];
435 // sum += 0.25*fLoGainSecondDeriv[k];
436 }
437
438 sum += (Float_t)summ;
439}
440
Note: See TracBrowser for help on using the repository browser.