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

Last change on this file since 5245 was 4752, checked in by tbretz, 21 years ago
*** empty log message ***
File size: 16.0 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 const Double_t deriv = *(p+1) - 2.* *(p) + *(p-1);
203 fHiGainFirstDeriv [i] = (6.0*deriv-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
235 Float_t higainklo = fHiGainSecondDeriv[klo];
236 Float_t higainkhi = fHiGainSecondDeriv[khi];
237 while (x<upper-0.3)
238 {
239
240 x += step;
241 a -= step;
242 b += step;
243
244 y = a*klocont
245 + b*khicont
246 + (a*a*a-a)*higainklo
247 + (b*b*b-b)*higainkhi;
248
249 if (y > abmax)
250 {
251 abmax = y;
252 time = x;
253 }
254 }
255
256
257 if (time > upper-0.1)
258 {
259
260 upper = (Float_t)maxpos+1.;
261 lower = (Float_t)maxpos;
262 x = lower;
263 a = 1.;
264 b = 0.;
265 khi = maxpos+1;
266 klo = maxpos;
267 klocont = (Float_t)*(first+klo);
268 khicont = (Float_t)*(first+khi);
269
270 higainklo = fHiGainSecondDeriv[klo];
271 higainkhi = fHiGainSecondDeriv[khi];
272 while (x<upper-0.3)
273 {
274
275 x += step;
276 a -= step;
277 b += step;
278
279 y = a* klocont
280 + b* khicont
281 + (a*a*a-a)*higainklo
282 + (b*b*b-b)*higainkhi;
283
284 if (y > abmax)
285 {
286 abmax = y;
287 time = x;
288 }
289 }
290 }
291
292 const Float_t up = time+step-0.055;
293 const Float_t lo = time-step+0.055;
294 const Float_t maxpossave = time;
295
296 x = time;
297 a = upper - x;
298 b = x - lower;
299
300 step = 0.04; // step size of 83 ps
301
302 higainklo = fHiGainSecondDeriv[klo];
303 higainkhi = fHiGainSecondDeriv[khi];
304 while (x<up)
305 {
306
307 x += step;
308 a -= step;
309 b += step;
310
311 y = a* klocont
312 + b* khicont
313 + (a*a*a-a)*higainklo
314 + (b*b*b-b)*higainkhi;
315
316 if (y > abmax)
317 {
318 abmax = y;
319 time = x;
320 }
321
322 }
323
324 if (time < klo + 0.02)
325 {
326 klo--;
327 khi--;
328 klocont = (Float_t)*(first+klo);
329 khicont = (Float_t)*(first+khi);
330 upper--;
331 lower--;
332 }
333
334 x = maxpossave;
335 a = upper - x;
336 b = x - lower;
337
338 higainklo = fHiGainSecondDeriv[klo];
339 higainkhi = fHiGainSecondDeriv[khi];
340 while (x>lo)
341 {
342
343 x -= step;
344 a += step;
345 b -= step;
346
347 y = a* klocont
348 + b* khicont
349 + (a*a*a-a)*higainklo
350 + (b*b*b-b)*higainkhi;
351
352 if (y > abmax)
353 {
354 abmax = y;
355 time = x;
356 }
357 }
358
359 const Float_t pedes = ped.GetPedestal();
360 const Float_t halfmax = pedes + (abmax - pedes)/2.;
361
362 //
363 // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
364 // First, find the right FADC slice:
365 //
366 klo = maxpos;
367 while (klo > maxpos-fStartBeforeMax)
368 {
369 if (*(first+klo) < (Byte_t)halfmax)
370 break;
371 klo--;
372 }
373
374 //
375 // Loop from the beginning of the slice upwards to reach the halfmax:
376 // With means of bisection:
377 //
378 x = (Float_t)klo;
379 a = 1.;
380 b = 0.;
381 klocont = (Float_t)*(first+klo);
382 khicont = (Float_t)*(first+klo+1);
383 time = x;
384
385 step = 0.5;
386 Bool_t back = kFALSE;
387
388 while (step > fResolution)
389 {
390
391 if (back)
392 {
393 x -= step;
394 a += step;
395 b -= step;
396 }
397 else
398 {
399 x += step;
400 a -= step;
401 b += step;
402 }
403
404 y = a*klocont
405 + b*khicont
406 + (a*a*a-a)*fHiGainSecondDeriv[klo]
407 + (b*b*b-b)*fHiGainSecondDeriv[khi];
408
409 if (y >= halfmax)
410 back = kTRUE;
411 else
412 back = kFALSE;
413
414 step /= 2.;
415
416 }
417
418 time = (Float_t)fHiGainFirst + x;
419 dtime = fResolution;
420}
421
422
423// --------------------------------------------------------------------------
424//
425// Calculates the arrival time for each pixel
426//
427void MExtractTimeFastSpline::FindTimeLoGain(Byte_t *first, Float_t &time, Float_t &dtime,
428 Byte_t &sat, const MPedestalPix &ped) const
429{
430
431 const Int_t range = fLoGainLast - fLoGainFirst + 1;
432 const Byte_t *end = first + range;
433 Byte_t *p = first;
434 Byte_t max = 0;
435 Byte_t maxpos = 0;
436
437 //
438 // Check for saturation in all other slices
439 //
440 while (++p<end)
441 {
442 if (*p > max)
443 {
444 max = *p;
445 maxpos = p-first;
446 }
447
448 if (*p >= fSaturationLimit)
449 {
450 sat++;
451 break;
452 }
453 }
454
455 if (sat)
456 return;
457
458 if (maxpos < 2)
459 return;
460
461 Float_t pp;
462
463 p = first;
464 fLoGainSecondDeriv[0] = 0.;
465 fLoGainFirstDeriv[0] = 0.;
466
467 for (Int_t i=1;i<range-1;i++)
468 {
469 p++;
470 pp = fLoGainSecondDeriv[i-1] + 4.;
471 fLoGainSecondDeriv[i] = -1.0/pp;
472 fLoGainFirstDeriv [i] = *(p+1) - 2.* *(p) + *(p-1);
473 fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
474 }
475
476 fLoGainSecondDeriv[range-1] = 0.;
477
478 for (Int_t k=range-2;k>0;k--)
479 fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k];
480 for (Int_t k=range-2;k>0;k--)
481 fLoGainSecondDeriv[k] /= 6.;
482
483 //
484 // Now find the maximum
485 //
486 Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
487 Float_t lower = (Float_t)maxpos-1.;
488 Float_t upper = (Float_t)maxpos;
489 Float_t x = lower;
490 Float_t y = 0.;
491 Float_t a = 1.;
492 Float_t b = 0.;
493 Int_t klo = maxpos-1;
494 Int_t khi = maxpos;
495 Float_t klocont = (Float_t)*(first+klo);
496 Float_t khicont = (Float_t)*(first+khi);
497 time = upper;
498 Float_t abmax = khicont;
499
500 //
501 // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to
502 // interval maxpos+1.
503 //
504 while (x<upper-0.3)
505 {
506
507 x += step;
508 a -= step;
509 b += step;
510
511 y = a*klocont
512 + b*khicont
513 + (a*a*a-a)*fLoGainSecondDeriv[klo]
514 + (b*b*b-b)*fLoGainSecondDeriv[khi];
515
516 if (y > abmax)
517 {
518 abmax = y;
519 time = x;
520 }
521
522 }
523
524 if (time > upper-0.1)
525 {
526
527 upper = (Float_t)maxpos+1.;
528 lower = (Float_t)maxpos;
529 x = lower;
530 a = 1.;
531 b = 0.;
532 khi = maxpos+1;
533 klo = maxpos;
534 klocont = (Float_t)*(first+klo);
535 khicont = (Float_t)*(first+khi);
536
537 while (x<upper-0.3)
538 {
539
540 x += step;
541 a -= step;
542 b += step;
543
544 y = a* klocont
545 + b* khicont
546 + (a*a*a-a)*fLoGainSecondDeriv[klo]
547 + (b*b*b-b)*fLoGainSecondDeriv[khi];
548
549 if (y > abmax)
550 {
551 abmax = y;
552 time = x;
553 }
554 }
555 }
556
557 const Float_t up = time+step-0.055;
558 const Float_t lo = time-step+0.055;
559 const Float_t maxpossave = time;
560
561 x = time;
562 a = upper - x;
563 b = x - lower;
564
565 step = 0.025; // step size of 165 ps
566
567 while (x<up)
568 {
569
570 x += step;
571 a -= step;
572 b += step;
573
574 y = a* klocont
575 + b* khicont
576 + (a*a*a-a)*fLoGainSecondDeriv[klo]
577 + (b*b*b-b)*fLoGainSecondDeriv[khi];
578
579 if (y > abmax)
580 {
581 abmax = y;
582 time = x;
583 }
584
585 }
586
587 if (time < klo + 0.01)
588 {
589 klo--;
590 khi--;
591 klocont = (Float_t)*(first+klo);
592 khicont = (Float_t)*(first+khi);
593 upper--;
594 lower--;
595 }
596
597 x = maxpossave;
598 a = upper - x;
599 b = x - lower;
600
601 while (x>lo)
602 {
603
604 x -= step;
605 a += step;
606 b -= step;
607
608 y = a* klocont
609 + b* khicont
610 + (a*a*a-a)*fLoGainSecondDeriv[klo]
611 + (b*b*b-b)*fLoGainSecondDeriv[khi];
612
613 if (y > abmax)
614 {
615 abmax = y;
616 time = x;
617 }
618
619 }
620
621 const Float_t pedes = ped.GetPedestal();
622 const Float_t halfmax = pedes + (abmax - pedes)/2.;
623
624 //
625 // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
626 // First, find the right FADC slice:
627 //
628 klo = maxpos;
629 while (klo > maxpos-4)
630 {
631 if (*(first+klo) < (Byte_t)halfmax)
632 break;
633 klo--;
634 }
635
636 //
637 // Loop from the beginning of the slice upwards to reach the halfmax:
638 // With means of bisection:
639 //
640 x = (Float_t)klo;
641 a = 1.;
642 b = 0.;
643 klocont = (Float_t)*(first+klo);
644 khicont = (Float_t)*(first+klo+1);
645 time = x;
646
647 step = 0.5;
648 Bool_t back = kFALSE;
649
650 while (step > fResolution)
651 {
652
653 if (back)
654 {
655 x -= step;
656 a += step;
657 b -= step;
658 }
659 else
660 {
661 x += step;
662 a -= step;
663 b += step;
664 }
665
666 y = a*klocont
667 + b*khicont
668 + (a*a*a-a)*fLoGainSecondDeriv[klo]
669 + (b*b*b-b)*fLoGainSecondDeriv[khi];
670
671 if (y >= halfmax)
672 back = kTRUE;
673 else
674 back = kFALSE;
675
676 step /= 2.;
677
678 }
679
680 time = (Float_t)fLoGainFirst + x;
681 dtime = fResolution;
682}
683
684// --------------------------------------------------------------------------
685//
686// In addition to the resources of the base-class MExtractor:
687// MJPedestal.MExtractor.Resolution: 0.003
688// MJPedestal.MExtractor.RiseTime: 1.5
689//
690Int_t MExtractTimeFastSpline::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
691{
692 Bool_t rc = kFALSE;
693
694 if (IsEnvDefined(env, prefix, "HiGainWindowSize", print))
695 {
696 SetResolution(GetEnvValue(env, prefix, "Resolution", fResolution));
697 rc = kTRUE;
698 }
699 if (IsEnvDefined(env, prefix, "LoGainWindowSize", print))
700 {
701 SetRiseTime(GetEnvValue(env, prefix, "RiseTime", fRiseTime));
702 rc = kTRUE;
703 }
704
705 return MExtractTime::ReadEnv(env, prefix, print) ? kTRUE : rc;
706}
707
708void MExtractTimeFastSpline::Print(Option_t *o) const
709{
710 *fLog << all;
711 *fLog << GetDescriptor() << ":" << endl;
712 *fLog << " Resolution: " << fResolution << endl;
713 *fLog << " RiseTime: " << fRiseTime << endl;
714 MExtractTime::Print(o);
715}
Note: See TracBrowser for help on using the repository browser.