source: branches/Mars_McMismatchStudy/msignal/MExtractTimeFastSpline.cc@ 18459

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