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

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