source: trunk/MagicSoft/Mars/msignal/MExtractAmplitudeSpline.cc@ 4693

Last change on this file since 4693 was 4378, checked in by gaug, 20 years ago
*** empty log message ***
File size: 15.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 "MExtractAmplitudeSpline.h"
44
45#include "MExtractedSignalCam.h"
46
47#include "MLog.h"
48#include "MLogManip.h"
49
50ClassImp(MExtractAmplitudeSpline);
51
52using namespace std;
53
54const Byte_t MExtractAmplitudeSpline::fgHiGainFirst = 2;
55const Byte_t MExtractAmplitudeSpline::fgHiGainLast = 14;
56const Byte_t MExtractAmplitudeSpline::fgLoGainFirst = 3;
57const Byte_t MExtractAmplitudeSpline::fgLoGainLast = 14;
58const Float_t MExtractAmplitudeSpline::fgResolution = 0.003;
59// --------------------------------------------------------------------------
60//
61// Default constructor.
62//
63// Calls:
64// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
65//
66MExtractAmplitudeSpline::MExtractAmplitudeSpline(const char *name, const char *title)
67 : fHiGainSignal(NULL), fLoGainSignal(NULL),
68 fHiGainFirstDeriv(NULL), fLoGainFirstDeriv(NULL),
69 fHiGainSecondDeriv(NULL), fLoGainSecondDeriv(NULL)
70{
71
72 fName = name ? name : "MExtractAmplitudeSpline";
73 fTitle = title ? title : "Signal Extractor for a fixed FADC window using a fast spline";
74
75 SetResolution();
76 SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
77}
78
79MExtractAmplitudeSpline::~MExtractAmplitudeSpline()
80{
81
82 if (fHiGainSignal)
83 delete [] fHiGainSignal;
84 if (fLoGainSignal)
85 delete [] fLoGainSignal;
86 if (fHiGainFirstDeriv)
87 delete [] fHiGainFirstDeriv;
88 if (fLoGainFirstDeriv)
89 delete [] fLoGainFirstDeriv;
90 if (fHiGainSecondDeriv)
91 delete [] fHiGainSecondDeriv;
92 if (fLoGainSecondDeriv)
93 delete [] fLoGainSecondDeriv;
94
95}
96
97// --------------------------------------------------------------------------
98//
99// SetRange:
100//
101// Checks:
102// - if the window defined by (fHiGainLast-fHiGainFirst-1) are odd, subtract one
103// - if the window defined by (fLoGainLast-fLoGainFirst-1) are odd, subtract one
104// - if the Hi Gain window is smaller than 2, set fHiGainLast to fHiGainFirst+1
105// - if the Lo Gain window is smaller than 2, set fLoGainLast to fLoGainFirst+1
106//
107// Calls:
108// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
109//
110// Sets:
111// - fNumHiGainSamples to: (Float_t)(fHiGainLast-fHiGainFirst+1)
112// - fNumLoGainSamples to: (Float_t)(fLoGainLast-fLoGainFirst+1)
113// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
114// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)
115//
116void MExtractAmplitudeSpline::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
117{
118
119
120 MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
121
122 fNumHiGainSamples = 2.;
123 fNumLoGainSamples = 2.;
124
125 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
126 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
127
128 fHiLoLast = 0;
129
130}
131
132// --------------------------------------------------------------------------
133//
134// ReInit
135//
136// Calls:
137// - MExtractor::ReInit(pList);
138// - fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
139// fLoGainFirst, fLoGainLast, fNumLoGainSamples);
140//
141// Deletes all arrays, if not NULL
142// Creates new arrays according to the extraction range
143//
144Bool_t MExtractAmplitudeSpline::ReInit(MParList *pList)
145{
146
147 if (!MExtractor::ReInit(pList))
148 return kFALSE;
149
150 fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples,
151 fLoGainFirst, fLoGainLast, fNumLoGainSamples);
152
153 if (fHiGainSignal)
154 delete [] fHiGainSignal;
155 if (fHiGainFirstDeriv)
156 delete [] fHiGainFirstDeriv;
157 if (fHiGainSecondDeriv)
158 delete [] fHiGainSecondDeriv;
159 if (fLoGainSignal)
160 delete [] fLoGainSignal;
161 if (fLoGainFirstDeriv)
162 delete [] fLoGainFirstDeriv;
163 if (fLoGainSecondDeriv)
164 delete [] fLoGainSecondDeriv;
165
166 Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast;
167
168 fHiGainSignal = new Float_t[range];
169 memset(fHiGainSignal,0,range*sizeof(Float_t));
170 fHiGainFirstDeriv = new Float_t[range];
171 memset(fHiGainFirstDeriv,0,range*sizeof(Float_t));
172 fHiGainSecondDeriv = new Float_t[range];
173 memset(fHiGainSecondDeriv,0,range*sizeof(Float_t));
174
175 *fLog << endl;
176 *fLog << inf << GetDescriptor() << ": Using for High-Gain Extraction " << range
177 << " FADC samples from "
178 << Form("%s%2i%s"," High Gain slice ",(Int_t)fHiGainFirst," to (including) ")
179 << Form("%s%2i",fHiLoLast ? "Low Gain slice " : " High Gain slice ",
180 fHiLoLast ? (Int_t)fHiLoLast : (Int_t)fHiGainLast )
181 << endl;
182
183 range = fLoGainLast - fLoGainFirst + 1;
184
185 fLoGainSignal = new Float_t[range];
186 memset(fLoGainSignal,0,range*sizeof(Float_t));
187 fLoGainFirstDeriv = new Float_t[range];
188 memset(fLoGainFirstDeriv,0,range*sizeof(Float_t));
189 fLoGainSecondDeriv = new Float_t[range];
190 memset(fLoGainSecondDeriv,0,range*sizeof(Float_t));
191
192 *fLog << endl;
193 *fLog << inf << GetDescriptor() << ": Using for Low-Gain Extraction " << range
194 << " FADC samples from "
195 << Form("%s%2i%s%2i"," Low Gain slice ",(Int_t)fLoGainFirst,
196 " to (including) ",(Int_t)fLoGainLast) << endl;
197
198 return kTRUE;
199}
200
201
202// --------------------------------------------------------------------------
203//
204// FindSignalHiGain:
205//
206// - Loop from ptr to (ptr+fHiGainLast-fHiGainFirst)
207// - Sum up contents of *ptr
208// - If *ptr is greater than fSaturationLimit, raise sat by 1
209//
210// - If fHiLoLast is not 0, loop also from logain to (logain+fHiLoLast)
211// - Sum up contents of logain
212// - If *logain is greater than fSaturationLimit, raise sat by 1
213//
214void MExtractAmplitudeSpline::FindSignalHiGain(Byte_t *ptr, Byte_t *logain, Float_t &sum, Byte_t &sat) const
215{
216
217 Int_t count = 0;
218 Float_t abmaxpos = 0.;
219 Byte_t max = 0;
220 Byte_t maxpos = 0;
221
222 Int_t range = fHiGainLast - fHiGainFirst + 1;
223 Byte_t *end = ptr + range;
224 Byte_t *p = ptr;
225 //
226 // Check for saturation in all other slices
227 //
228 while (++p<end)
229 {
230
231 fHiGainSignal[count] = (Float_t)*p;
232
233 if (*p > max)
234 {
235 max = *p;
236 maxpos = count;
237 }
238
239 count++;
240
241 if (*p >= fSaturationLimit)
242 {
243 sat++;
244 break;
245 }
246 }
247
248 if (fHiLoLast != 0)
249 {
250
251 p = logain;
252 end = logain + fHiLoLast + 1;
253
254 while (p<end)
255 {
256
257 fHiGainSignal[count] = (Float_t)*p;
258
259 if (*p > max)
260 {
261 max = *p;
262 maxpos = count;
263 }
264
265 range++;
266 count++;
267
268 if (*p++ >= fSaturationLimit)
269 {
270 sat++;
271 break;
272 }
273 }
274 }
275
276 //
277 // allow one saturated slice
278 //
279 if (sat > 1)
280 return;
281
282 //
283 // Don't start if the maxpos is too close to the left limit.
284 //
285 if (maxpos < 2)
286 return;
287
288 Float_t pp;
289 fHiGainSecondDeriv[0] = 0.;
290 fHiGainFirstDeriv[0] = 0.;
291
292 for (Int_t i=1;i<range-1;i++)
293 {
294 pp = fHiGainSecondDeriv[i-1] + 4.;
295 fHiGainSecondDeriv[i] = -1.0/pp;
296 fHiGainFirstDeriv [i] = fHiGainSignal[i+1] - fHiGainSignal[i] - fHiGainSignal[i] + fHiGainSignal[i-1];
297 fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
298 p++;
299 }
300
301 fHiGainSecondDeriv[range-1] = 0.;
302 for (Int_t k=range-2;k>=0;k--)
303 fHiGainSecondDeriv[k] = (fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k])/6.;
304
305 //
306 // Now find the maximum
307 //
308 Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
309 Float_t lower = (Float_t)maxpos-1.;
310 Float_t upper = (Float_t)maxpos;
311 Float_t x = lower;
312 Float_t y = 0.;
313 Float_t a = 1.;
314 Float_t b = 0.;
315 Int_t klo = maxpos-1;
316 Int_t khi = maxpos;
317 Float_t klocont = fHiGainSignal[klo];
318 Float_t khicont = fHiGainSignal[khi];
319 sum = khicont;
320 abmaxpos = lower;
321
322 //
323 // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to
324 // interval maxpos+1.
325 //
326 while (x<upper-0.3)
327 {
328
329 x += step;
330 a -= step;
331 b += step;
332
333 y = a*klocont
334 + b*khicont
335 + (a*a*a-a)*fHiGainSecondDeriv[klo]
336 + (b*b*b-b)*fHiGainSecondDeriv[khi];
337
338 if (y > sum)
339 {
340 sum = y;
341 abmaxpos = x;
342 }
343 }
344
345 if (abmaxpos > upper-0.1)
346 {
347
348 upper = (Float_t)maxpos+1;
349 lower = (Float_t)maxpos;
350 x = lower;
351 a = 1.;
352 b = 0.;
353 khi = maxpos+1;
354 klo = maxpos;
355 klocont = fHiGainSignal[klo];
356 khicont = fHiGainSignal[khi];
357
358 while (x<upper-0.3)
359 {
360
361 x += step;
362 a -= step;
363 b += step;
364
365 y = a* klocont
366 + b* khicont
367 + (a*a*a-a)*fHiGainSecondDeriv[klo]
368 + (b*b*b-b)*fHiGainSecondDeriv[khi];
369
370 if (y > sum)
371 {
372 sum = y;
373 abmaxpos = x;
374 }
375 }
376 }
377
378 const Float_t up = abmaxpos+step-0.055;
379 const Float_t lo = abmaxpos-step+0.055;
380 const Float_t maxpossave = abmaxpos;
381
382 x = abmaxpos;
383 a = upper - x;
384 b = x - lower;
385
386 step = 0.04; // step size of 83 ps
387
388 while (x<up)
389 {
390
391 x += step;
392 a -= step;
393 b += step;
394
395 y = a* klocont
396 + b* khicont
397 + (a*a*a-a)*fHiGainSecondDeriv[klo]
398 + (b*b*b-b)*fHiGainSecondDeriv[khi];
399
400 if (y > sum)
401 {
402 sum = y;
403 abmaxpos = x;
404 }
405 }
406
407 if (abmaxpos < klo + 0.02)
408 {
409 klo--;
410 khi--;
411 klocont = fHiGainSignal[klo];
412 khicont = fHiGainSignal[khi];
413 upper--;
414 lower--;
415 }
416
417 x = maxpossave;
418 a = upper - x;
419 b = x - lower;
420
421 while (x>lo)
422 {
423
424 x -= step;
425 a += step;
426 b -= step;
427
428 y = a* klocont
429 + b* khicont
430 + (a*a*a-a)*fHiGainSecondDeriv[klo]
431 + (b*b*b-b)*fHiGainSecondDeriv[khi];
432
433 if (y > sum)
434 sum = y;
435 }
436}
437
438
439// --------------------------------------------------------------------------
440//
441// FindSignalLoGain:
442//
443// - Loop from ptr to (ptr+fLoGainLast-fLoGainFirst)
444// - Sum up contents of *ptr
445// - If *ptr is greater than fSaturationLimit, raise sat by 1
446//
447void MExtractAmplitudeSpline::FindSignalLoGain(Byte_t *ptr, Float_t &sum, Byte_t &sat) const
448{
449
450 Int_t count = 0;
451 Float_t abmaxpos = 0.;
452 Byte_t max = 0;
453 Byte_t maxpos = 0;
454
455 Int_t range = fLoGainLast - fLoGainFirst + 1;
456 Byte_t *end = ptr + range;
457 Byte_t *p = ptr;
458 //
459 // Check for saturation in all other slices
460 //
461 while (++p<end)
462 {
463
464 fLoGainSignal[count] = (Float_t)*p;
465
466 if (*p > max)
467 {
468 max = *p;
469 maxpos = count;
470 }
471
472 range++;
473 count++;
474
475 if (*p >= fSaturationLimit)
476 {
477 sat++;
478 break;
479 }
480 }
481
482
483 //
484 // allow one saturated slice
485 //
486 if (sat > 1)
487 return;
488
489 //
490 // Don't start if the maxpos is too close to the left limit.
491 //
492 if (maxpos < 2)
493 return;
494
495 Float_t pp;
496 fLoGainSecondDeriv[0] = 0.;
497 fLoGainFirstDeriv[0] = 0.;
498
499 for (Int_t i=1;i<range-1;i++)
500 {
501 pp = fLoGainSecondDeriv[i-1] + 4.;
502 fLoGainSecondDeriv[i] = -1.0/pp;
503 fLoGainFirstDeriv [i] = fLoGainSignal[i+1] - fLoGainSignal[i] - fLoGainSignal[i] + fLoGainSignal[i-1];
504 fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
505 p++;
506 }
507
508 fLoGainSecondDeriv[range-1] = 0.;
509 for (Int_t k=range-2;k>=0;k--)
510 fLoGainSecondDeriv[k] = (fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k])/6.;
511
512 //
513 // Now find the maximum
514 //
515 Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
516 Float_t lower = (Float_t)maxpos-1.;
517 Float_t upper = (Float_t)maxpos;
518 Float_t x = lower;
519 Float_t y = 0.;
520 Float_t a = 1.;
521 Float_t b = 0.;
522 Int_t klo = maxpos-1;
523 Int_t khi = maxpos;
524 Float_t klocont = fLoGainSignal[klo];
525 Float_t khicont = fLoGainSignal[khi];
526 sum = khicont;
527 abmaxpos = lower;
528
529 //
530 // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to
531 // interval maxpos+1.
532 //
533 while (x<upper-0.3)
534 {
535
536 x += step;
537 a -= step;
538 b += step;
539
540 y = a*klocont
541 + b*khicont
542 + (a*a*a-a)*fLoGainSecondDeriv[klo]
543 + (b*b*b-b)*fLoGainSecondDeriv[khi];
544
545 if (y > sum)
546 {
547 sum = y;
548 abmaxpos = x;
549 }
550 }
551
552 if (abmaxpos > upper-0.1)
553 {
554
555 upper = (Float_t)maxpos+1;
556 lower = (Float_t)maxpos;
557 x = lower;
558 a = 1.;
559 b = 0.;
560 khi = maxpos+1;
561 klo = maxpos;
562 klocont = fLoGainSignal[klo];
563 khicont = fLoGainSignal[khi];
564
565 while (x<upper-0.3)
566 {
567
568 x += step;
569 a -= step;
570 b += step;
571
572 y = a* klocont
573 + b* khicont
574 + (a*a*a-a)*fLoGainSecondDeriv[klo]
575 + (b*b*b-b)*fLoGainSecondDeriv[khi];
576
577 if (y > sum)
578 {
579 sum = y;
580 abmaxpos = x;
581 }
582 }
583 }
584
585 const Float_t up = abmaxpos+step-0.055;
586 const Float_t lo = abmaxpos-step+0.055;
587 const Float_t maxpossave = abmaxpos;
588
589 x = abmaxpos;
590 a = upper - x;
591 b = x - lower;
592
593 step = 0.04; // step size of 83 ps
594
595 while (x<up)
596 {
597
598 x += step;
599 a -= step;
600 b += step;
601
602 y = a* klocont
603 + b* khicont
604 + (a*a*a-a)*fLoGainSecondDeriv[klo]
605 + (b*b*b-b)*fLoGainSecondDeriv[khi];
606
607 if (y > sum)
608 {
609 sum = y;
610 abmaxpos = x;
611 }
612 }
613
614 if (abmaxpos < klo + 0.02)
615 {
616 klo--;
617 khi--;
618 klocont = fLoGainSignal[klo];
619 khicont = fLoGainSignal[khi];
620 upper--;
621 lower--;
622 }
623
624 x = maxpossave;
625 a = upper - x;
626 b = x - lower;
627
628 while (x>lo)
629 {
630
631 x -= step;
632 a += step;
633 b -= step;
634
635 y = a* klocont
636 + b* khicont
637 + (a*a*a-a)*fLoGainSecondDeriv[klo]
638 + (b*b*b-b)*fLoGainSecondDeriv[khi];
639
640 if (y > sum)
641 sum = y;
642 }
643}
644
Note: See TracBrowser for help on using the repository browser.