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

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