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

Last change on this file since 5633 was 5604, checked in by gaug, 20 years ago
*** empty log message ***
File size: 11.5 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 else
271 {
272
273 //
274 // There are two overlap slices which we have to treat sepatately
275 //
276 // First, the last high-gain slice as center
277 //
278 summ += *ptr;
279 i++;
280
281 if (*logain >= fSaturationLimit)
282 sat++;
283
284 if (fHiLoLast == 1)
285 sum += (Float_t)*logain/2;
286 else
287 {
288 pp = fHiGainSecondDeriv[i-1] + 4.;
289 fHiGainSecondDeriv[i] = -1.0/pp;
290 fHiGainFirstDeriv [i] = *(logain) - 2.* *(ptr) + *(ptr-1);
291 fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
292
293 //
294 // Second, the first low-gain slice as center
295 //
296 summ += *logain;
297 i++;
298
299 pp = fHiGainSecondDeriv[i-1] + 4.;
300 fHiGainSecondDeriv[i] = -1.0/pp;
301 fHiGainFirstDeriv [i] = *(logain+1) - 2.* *(logain) + *(ptr);
302 fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
303
304 if (*logain++ >= fSaturationLimit)
305 sat++;
306
307 end = logain + fHiLoLast - 2;
308
309 while (logain<end)
310 {
311
312 summ += *logain;
313 i++;
314
315 pp = fHiGainSecondDeriv[i-1] + 4.;
316 fHiGainSecondDeriv[i] = -1.0/pp;
317 fHiGainFirstDeriv [i] = *(logain+1) - 2.* *(logain) + *(logain-1);
318 fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
319
320 if (*logain++ >= fSaturationLimit)
321 sat++;
322
323 }
324 }
325 }
326
327 fHiGainSecondDeriv[range-1] = 0.;
328
329 for (Int_t k=range-2;k>0;k--)
330 {
331 fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k];
332 sum += 0.25*fHiGainSecondDeriv[k];
333 }
334
335 sum += (Float_t)summ;
336}
337
338// --------------------------------------------------------------------------
339//
340// FindSignalLoGain:
341//
342// - Loop from ptr to (ptr+fLoGainLast-fLoGainFirst)
343// - Sum up contents of *ptr
344// - If *ptr is greater than fSaturationLimit, raise sat by 1
345//
346void MExtractFixedWindowSpline::FindSignalLoGain(Byte_t *ptr, Float_t &sum, Byte_t &sat) const
347{
348
349 const Byte_t *end = ptr + fLoGainLast - fLoGainFirst;
350 Int_t range = fLoGainLast - fLoGainFirst + 1;
351
352 Float_t pp;
353 Int_t i = 0;
354
355 fLoGainSecondDeriv[0] = 0.;
356 fLoGainFirstDeriv[0] = 0.;
357
358 Int_t summ = 0;
359 sum = (Float_t)*ptr++/2.;
360 //
361 // Check for saturation in all other slices
362 //
363 while (ptr<end)
364 {
365
366 summ += *ptr;
367 i++;
368
369 pp = fLoGainSecondDeriv[i-1] + 4.;
370 fLoGainSecondDeriv[i] = -1.0/pp;
371 fLoGainFirstDeriv [i] = *(ptr+1) - 2.* *(ptr) + *(ptr-1);
372 fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
373
374 if (*ptr++ >= fSaturationLimit)
375 sat++;
376 }
377
378 if (*ptr++ >= fSaturationLimit)
379 sat++;
380
381 sum += (Float_t)*ptr/2.;
382 fLoGainSecondDeriv[++i] = 0.;
383
384 for (Int_t k=range-2;k>0;k--)
385 {
386 fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k];
387 sum += 0.25*fLoGainSecondDeriv[k];
388 }
389
390 sum += (Float_t)summ;
391}
392
Note: See TracBrowser for help on using the repository browser.