source: trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc@ 4742

Last change on this file since 4742 was 4732, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 15.8 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// MExtractTimeFastSpline
27//
28// Fast arrival Time extractor using a cubic spline algorithm of Numerical Recipes.
29// It returns the position of the half maximum between absolute maximum
30// and pedestal of the spline that interpolates the FADC slices.
31//
32// The precision of the half-maximum searches can be chosen by:
33// SetPrecision().
34//
35// The precision of the maximum-finder is fixed to 0.025 FADC units.
36//
37//////////////////////////////////////////////////////////////////////////////
38#include "MExtractTimeFastSpline.h"
39
40#include "MPedestalPix.h"
41
42#include "MLog.h"
43#include "MLogManip.h"
44
45
46ClassImp(MExtractTimeFastSpline);
47
48using namespace std;
49
50const Byte_t MExtractTimeFastSpline::fgHiGainFirst = 2;
51const Byte_t MExtractTimeFastSpline::fgHiGainLast = 14;
52const Byte_t MExtractTimeFastSpline::fgLoGainFirst = 3;
53const Byte_t MExtractTimeFastSpline::fgLoGainLast = 14;
54const Float_t MExtractTimeFastSpline::fgResolution = 0.003;
55const Float_t MExtractTimeFastSpline::fgRiseTime = 1.5;
56
57// --------------------------------------------------------------------------
58//
59// Default constructor.
60//
61// Calls:
62// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
63//
64// Initializes:
65// - fResolution to fgResolution
66//
67MExtractTimeFastSpline::MExtractTimeFastSpline(const char *name, const char *title)
68 : fHiGainFirstDeriv(NULL), fLoGainFirstDeriv(NULL),
69 fHiGainSecondDeriv(NULL), fLoGainSecondDeriv(NULL)
70{
71
72 fName = name ? name : "MExtractTimeFastSpline";
73 fTitle = title ? title : "Calculate photons arrival time using a fast spline";
74
75 SetResolution();
76 SetRiseTime ();
77 SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
78
79}
80
81MExtractTimeFastSpline::~MExtractTimeFastSpline()
82{
83
84 if (fHiGainFirstDeriv)
85 delete fHiGainFirstDeriv;
86 if (fLoGainFirstDeriv)
87 delete fLoGainFirstDeriv;
88 if (fHiGainSecondDeriv)
89 delete fHiGainSecondDeriv;
90 if (fLoGainSecondDeriv)
91 delete fLoGainSecondDeriv;
92
93}
94
95
96// --------------------------------------------------------------------------
97//
98// SetRange:
99//
100// Calls:
101// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
102// - Deletes x, if not NULL
103// - Creates x according to the range
104//
105void MExtractTimeFastSpline::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
106{
107
108 MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
109
110 if (fHiGainFirstDeriv)
111 delete fHiGainFirstDeriv;
112 if (fLoGainFirstDeriv)
113 delete fLoGainFirstDeriv;
114 if (fHiGainSecondDeriv)
115 delete fHiGainSecondDeriv;
116 if (fLoGainSecondDeriv)
117 delete fLoGainSecondDeriv;
118
119 Int_t range = fHiGainLast - fHiGainFirst + 1;
120
121 if (range < 2)
122 {
123 *fLog << warn << GetDescriptor()
124 << Form("%s%2i%s%2i%s",": Hi-Gain Extraction range [",(int)fHiGainFirst,","
125 ,fHiGainLast,"] too small, ") << endl;
126 *fLog << warn << GetDescriptor()
127 << " will move higher limit to obtain 4 slices " << endl;
128 SetRange(fHiGainFirst, fHiGainLast+4-range,fLoGainFirst,fLoGainLast);
129 range = fHiGainLast - fHiGainFirst + 1;
130 }
131
132
133 fHiGainFirstDeriv = new Float_t[range];
134 memset(fHiGainFirstDeriv,0,range*sizeof(Float_t));
135 fHiGainSecondDeriv = new Float_t[range];
136 memset(fHiGainSecondDeriv,0,range*sizeof(Float_t));
137
138 range = fLoGainLast - fLoGainFirst + 1;
139
140 if (range >= 2)
141 {
142
143 fLoGainFirstDeriv = new Float_t[range];
144 memset(fLoGainFirstDeriv,0,range*sizeof(Float_t));
145 fLoGainSecondDeriv = new Float_t[range];
146 memset(fLoGainSecondDeriv,0,range*sizeof(Float_t));
147
148 }
149
150}
151
152
153// --------------------------------------------------------------------------
154//
155// Calculates the arrival time for each pixel
156//
157void MExtractTimeFastSpline::FindTimeHiGain(Byte_t *first, Float_t &time, Float_t &dtime,
158 Byte_t &sat, const MPedestalPix &ped) const
159{
160
161 const Int_t range = fHiGainLast - fHiGainFirst + 1;
162 const Byte_t *end = first + range;
163 Byte_t *p = first;
164 Byte_t max = 0;
165 Byte_t maxpos = 0;
166
167 //
168 // Check for saturation in all other slices
169 //
170 while (++p<end)
171 {
172 if (*p > max)
173 {
174 max = *p;
175 maxpos = p-first;
176 }
177
178 if (*p >= fSaturationLimit)
179 {
180 sat++;
181 break;
182 }
183 }
184
185 if (sat)
186 return;
187
188 if (maxpos < 2)
189 return;
190
191 Float_t pp;
192
193 p = first;
194 fHiGainSecondDeriv[0] = 0.;
195 fHiGainFirstDeriv[0] = 0.;
196
197 for (Int_t i=1;i<range-1;i++)
198 {
199 p++;
200 pp = fHiGainSecondDeriv[i-1] + 4.;
201 fHiGainSecondDeriv[i] = -1.0/pp;
202 fHiGainFirstDeriv [i] = *(p+1) - 2.* *(p) + *(p-1);
203 fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
204 }
205
206 fHiGainSecondDeriv[range-1] = 0.;
207
208 for (Int_t k=range-2;k>0;k--)
209 fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k];
210 for (Int_t k=range-2;k>0;k--)
211 fHiGainSecondDeriv[k] /= 6.;
212
213 //
214 // Now find the maximum
215 //
216 Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
217 Float_t lower = (Float_t)maxpos-1.;
218 Float_t upper = (Float_t)maxpos;
219 Float_t x = lower;
220 Float_t y = 0.;
221 Float_t a = 1.;
222 Float_t b = 0.;
223 Int_t klo = maxpos-1;
224 Int_t khi = maxpos;
225 Float_t klocont = (Float_t)*(first+klo);
226 Float_t khicont = (Float_t)*(first+khi);
227 time = upper;
228 Float_t abmax = khicont;
229
230 //
231 // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to
232 // interval maxpos+1.
233 //
234 while (x<upper-0.3)
235 {
236
237 x += step;
238 a -= step;
239 b += step;
240
241 y = a*klocont
242 + b*khicont
243 + (a*a*a-a)*fHiGainSecondDeriv[klo]
244 + (b*b*b-b)*fHiGainSecondDeriv[khi];
245
246 if (y > abmax)
247 {
248 abmax = y;
249 time = x;
250 }
251 }
252
253
254 if (time > upper-0.1)
255 {
256
257 upper = (Float_t)maxpos+1.;
258 lower = (Float_t)maxpos;
259 x = lower;
260 a = 1.;
261 b = 0.;
262 khi = maxpos+1;
263 klo = maxpos;
264 klocont = (Float_t)*(first+klo);
265 khicont = (Float_t)*(first+khi);
266
267 while (x<upper-0.3)
268 {
269
270 x += step;
271 a -= step;
272 b += step;
273
274 y = a* klocont
275 + b* khicont
276 + (a*a*a-a)*fHiGainSecondDeriv[klo]
277 + (b*b*b-b)*fHiGainSecondDeriv[khi];
278
279 if (y > abmax)
280 {
281 abmax = y;
282 time = x;
283 }
284 }
285 }
286
287 const Float_t up = time+step-0.055;
288 const Float_t lo = time-step+0.055;
289 const Float_t maxpossave = time;
290
291 x = time;
292 a = upper - x;
293 b = x - lower;
294
295 step = 0.04; // step size of 83 ps
296
297 while (x<up)
298 {
299
300 x += step;
301 a -= step;
302 b += step;
303
304 y = a* klocont
305 + b* khicont
306 + (a*a*a-a)*fHiGainSecondDeriv[klo]
307 + (b*b*b-b)*fHiGainSecondDeriv[khi];
308
309 if (y > abmax)
310 {
311 abmax = y;
312 time = x;
313 }
314
315 }
316
317 if (time < klo + 0.02)
318 {
319 klo--;
320 khi--;
321 klocont = (Float_t)*(first+klo);
322 khicont = (Float_t)*(first+khi);
323 upper--;
324 lower--;
325 }
326
327 x = maxpossave;
328 a = upper - x;
329 b = x - lower;
330
331 while (x>lo)
332 {
333
334 x -= step;
335 a += step;
336 b -= step;
337
338 y = a* klocont
339 + b* khicont
340 + (a*a*a-a)*fHiGainSecondDeriv[klo]
341 + (b*b*b-b)*fHiGainSecondDeriv[khi];
342
343 if (y > abmax)
344 {
345 abmax = y;
346 time = x;
347 }
348 }
349
350 const Float_t pedes = ped.GetPedestal();
351 const Float_t halfmax = pedes + (abmax - pedes)/2.;
352
353 //
354 // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
355 // First, find the right FADC slice:
356 //
357 klo = maxpos;
358 while (klo > maxpos-fStartBeforeMax)
359 {
360 if (*(first+klo) < (Byte_t)halfmax)
361 break;
362 klo--;
363 }
364
365 //
366 // Loop from the beginning of the slice upwards to reach the halfmax:
367 // With means of bisection:
368 //
369 x = (Float_t)klo;
370 a = 1.;
371 b = 0.;
372 klocont = (Float_t)*(first+klo);
373 khicont = (Float_t)*(first+klo+1);
374 time = x;
375
376 step = 0.5;
377 Bool_t back = kFALSE;
378
379 while (step > fResolution)
380 {
381
382 if (back)
383 {
384 x -= step;
385 a += step;
386 b -= step;
387 }
388 else
389 {
390 x += step;
391 a -= step;
392 b += step;
393 }
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 >= halfmax)
401 back = kTRUE;
402 else
403 back = kFALSE;
404
405 step /= 2.;
406
407 }
408
409 time = (Float_t)fHiGainFirst + x;
410 dtime = fResolution;
411}
412
413
414// --------------------------------------------------------------------------
415//
416// Calculates the arrival time for each pixel
417//
418void MExtractTimeFastSpline::FindTimeLoGain(Byte_t *first, Float_t &time, Float_t &dtime,
419 Byte_t &sat, const MPedestalPix &ped) const
420{
421
422 const Int_t range = fLoGainLast - fLoGainFirst + 1;
423 const Byte_t *end = first + range;
424 Byte_t *p = first;
425 Byte_t max = 0;
426 Byte_t maxpos = 0;
427
428 //
429 // Check for saturation in all other slices
430 //
431 while (++p<end)
432 {
433 if (*p > max)
434 {
435 max = *p;
436 maxpos = p-first;
437 }
438
439 if (*p >= fSaturationLimit)
440 {
441 sat++;
442 break;
443 }
444 }
445
446 if (sat)
447 return;
448
449 if (maxpos < 2)
450 return;
451
452 Float_t pp;
453
454 p = first;
455 fLoGainSecondDeriv[0] = 0.;
456 fLoGainFirstDeriv[0] = 0.;
457
458 for (Int_t i=1;i<range-1;i++)
459 {
460 p++;
461 pp = fLoGainSecondDeriv[i-1] + 4.;
462 fLoGainSecondDeriv[i] = -1.0/pp;
463 fLoGainFirstDeriv [i] = *(p+1) - 2.* *(p) + *(p-1);
464 fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
465 }
466
467 fLoGainSecondDeriv[range-1] = 0.;
468
469 for (Int_t k=range-2;k>0;k--)
470 fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k];
471 for (Int_t k=range-2;k>0;k--)
472 fLoGainSecondDeriv[k] /= 6.;
473
474 //
475 // Now find the maximum
476 //
477 Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
478 Float_t lower = (Float_t)maxpos-1.;
479 Float_t upper = (Float_t)maxpos;
480 Float_t x = lower;
481 Float_t y = 0.;
482 Float_t a = 1.;
483 Float_t b = 0.;
484 Int_t klo = maxpos-1;
485 Int_t khi = maxpos;
486 Float_t klocont = (Float_t)*(first+klo);
487 Float_t khicont = (Float_t)*(first+khi);
488 time = upper;
489 Float_t abmax = khicont;
490
491 //
492 // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to
493 // interval maxpos+1.
494 //
495 while (x<upper-0.3)
496 {
497
498 x += step;
499 a -= step;
500 b += step;
501
502 y = a*klocont
503 + b*khicont
504 + (a*a*a-a)*fLoGainSecondDeriv[klo]
505 + (b*b*b-b)*fLoGainSecondDeriv[khi];
506
507 if (y > abmax)
508 {
509 abmax = y;
510 time = x;
511 }
512
513 }
514
515 if (time > upper-0.1)
516 {
517
518 upper = (Float_t)maxpos+1.;
519 lower = (Float_t)maxpos;
520 x = lower;
521 a = 1.;
522 b = 0.;
523 khi = maxpos+1;
524 klo = maxpos;
525 klocont = (Float_t)*(first+klo);
526 khicont = (Float_t)*(first+khi);
527
528 while (x<upper-0.3)
529 {
530
531 x += step;
532 a -= step;
533 b += step;
534
535 y = a* klocont
536 + b* khicont
537 + (a*a*a-a)*fLoGainSecondDeriv[klo]
538 + (b*b*b-b)*fLoGainSecondDeriv[khi];
539
540 if (y > abmax)
541 {
542 abmax = y;
543 time = x;
544 }
545 }
546 }
547
548 const Float_t up = time+step-0.055;
549 const Float_t lo = time-step+0.055;
550 const Float_t maxpossave = time;
551
552 x = time;
553 a = upper - x;
554 b = x - lower;
555
556 step = 0.025; // step size of 165 ps
557
558 while (x<up)
559 {
560
561 x += step;
562 a -= step;
563 b += step;
564
565 y = a* klocont
566 + b* khicont
567 + (a*a*a-a)*fLoGainSecondDeriv[klo]
568 + (b*b*b-b)*fLoGainSecondDeriv[khi];
569
570 if (y > abmax)
571 {
572 abmax = y;
573 time = x;
574 }
575
576 }
577
578 if (time < klo + 0.01)
579 {
580 klo--;
581 khi--;
582 klocont = (Float_t)*(first+klo);
583 khicont = (Float_t)*(first+khi);
584 upper--;
585 lower--;
586 }
587
588 x = maxpossave;
589 a = upper - x;
590 b = x - lower;
591
592 while (x>lo)
593 {
594
595 x -= step;
596 a += step;
597 b -= step;
598
599 y = a* klocont
600 + b* khicont
601 + (a*a*a-a)*fLoGainSecondDeriv[klo]
602 + (b*b*b-b)*fLoGainSecondDeriv[khi];
603
604 if (y > abmax)
605 {
606 abmax = y;
607 time = x;
608 }
609
610 }
611
612 const Float_t pedes = ped.GetPedestal();
613 const Float_t halfmax = pedes + (abmax - pedes)/2.;
614
615 //
616 // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
617 // First, find the right FADC slice:
618 //
619 klo = maxpos;
620 while (klo > maxpos-4)
621 {
622 if (*(first+klo) < (Byte_t)halfmax)
623 break;
624 klo--;
625 }
626
627 //
628 // Loop from the beginning of the slice upwards to reach the halfmax:
629 // With means of bisection:
630 //
631 x = (Float_t)klo;
632 a = 1.;
633 b = 0.;
634 klocont = (Float_t)*(first+klo);
635 khicont = (Float_t)*(first+klo+1);
636 time = x;
637
638 step = 0.5;
639 Bool_t back = kFALSE;
640
641 while (step > fResolution)
642 {
643
644 if (back)
645 {
646 x -= step;
647 a += step;
648 b -= step;
649 }
650 else
651 {
652 x += step;
653 a -= step;
654 b += step;
655 }
656
657 y = a*klocont
658 + b*khicont
659 + (a*a*a-a)*fLoGainSecondDeriv[klo]
660 + (b*b*b-b)*fLoGainSecondDeriv[khi];
661
662 if (y >= halfmax)
663 back = kTRUE;
664 else
665 back = kFALSE;
666
667 step /= 2.;
668
669 }
670
671 time = (Float_t)fLoGainFirst + x;
672 dtime = fResolution;
673}
674
675// --------------------------------------------------------------------------
676//
677// In addition to the resources of the base-class MExtractor:
678// MJPedestal.MExtractor.Resolution: 0.003
679// MJPedestal.MExtractor.RiseTime: 1.5
680//
681Int_t MExtractTimeFastSpline::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
682{
683 Bool_t rc = kFALSE;
684
685 if (IsEnvDefined(env, prefix, "HiGainWindowSize", print))
686 {
687 SetResolution(GetEnvValue(env, prefix, "Resolution", fResolution));
688 rc = kTRUE;
689 }
690 if (IsEnvDefined(env, prefix, "LoGainWindowSize", print))
691 {
692 SetRiseTime(GetEnvValue(env, prefix, "RiseTime", fRiseTime));
693 rc = kTRUE;
694 }
695
696 return MExtractTime::ReadEnv(env, prefix, print) ? kTRUE : rc;
697}
698
699void MExtractTimeFastSpline::Print(Option_t *o) const
700{
701 *fLog << all;
702 *fLog << GetDescriptor() << ":" << endl;
703 *fLog << " Resolution: " << fResolution << endl;
704 *fLog << " RiseTime: " << fRiseTime << endl;
705 MExtractTime::Print(o);
706}
Note: See TracBrowser for help on using the repository browser.