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

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