source: trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc@ 5218

Last change on this file since 5218 was 5215, checked in by gaug, 21 years ago
*** empty log message ***
File size: 21.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// MExtractTimeAndChargeSpline
27//
28// Fast Spline extractor using a cubic spline algorithm of Numerical Recipes.
29// It returns the integral below the interpolating spline.
30//
31// Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast)
32// to modify the ranges. Ranges have to be an even number. In case of odd
33// ranges, the last slice will be reduced by one.
34//
35// Defaults are:
36//
37// fHiGainFirst = fgHiGainFirst = 3
38// fHiGainLast = fgHiGainLast = 14
39// fLoGainFirst = fgLoGainFirst = 3
40// fLoGainLast = fgLoGainLast = 14
41//
42//////////////////////////////////////////////////////////////////////////////
43#include "MExtractTimeAndChargeSpline.h"
44
45#include "MPedestalPix.h"
46#include "MPedestalCam.h"
47
48#include "MLog.h"
49#include "MLogManip.h"
50
51#include "MParList.h"
52#include "MRawEvtPixelIter.h"
53
54ClassImp(MExtractTimeAndChargeSpline);
55
56using namespace std;
57
58const Byte_t MExtractTimeAndChargeSpline::fgHiGainFirst = 2;
59const Byte_t MExtractTimeAndChargeSpline::fgHiGainLast = 14;
60const Byte_t MExtractTimeAndChargeSpline::fgLoGainFirst = 3;
61const Byte_t MExtractTimeAndChargeSpline::fgLoGainLast = 14;
62const Float_t MExtractTimeAndChargeSpline::fgResolution = 0.001;
63const Float_t MExtractTimeAndChargeSpline::fgRiseTime = 2.0;
64const Float_t MExtractTimeAndChargeSpline::fgFallTime = 4.0;
65// --------------------------------------------------------------------------
66//
67// Default constructor.
68//
69// Calls:
70// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
71//
72// Initializes:
73// - fResolution to fgResolution
74// - fRatioMax2Fall to fgRatioMax2Fall
75// - fExtractCharges to kFALSE
76//
77MExtractTimeAndChargeSpline::MExtractTimeAndChargeSpline(const char *name, const char *title)
78 : fHiGainSignal(NULL), fLoGainSignal(NULL),
79 fHiGainFirstDeriv(NULL), fLoGainFirstDeriv(NULL),
80 fHiGainSecondDeriv(NULL), fLoGainSecondDeriv(NULL),
81 fAbMax(0.), fAbMaxPos(0.), fHalfMax(0.)
82{
83
84 fName = name ? name : "MExtractTimeAndChargeSpline";
85 fTitle = title ? title : "Calculate photons arrival time using a fast spline";
86
87 SetResolution();
88 SetRiseTime();
89 SetFallTime();
90
91 SetTimeType();
92 SetChargeType();
93
94 SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
95
96 fNumHiGainSamples = 1.;
97 fNumLoGainSamples = 1.;
98 fSqrtHiGainSamples = 1.;
99 fSqrtLoGainSamples = 1.;
100
101}
102
103MExtractTimeAndChargeSpline::~MExtractTimeAndChargeSpline()
104{
105
106 if (fHiGainSignal)
107 delete fHiGainSignal;
108 if (fLoGainSignal)
109 delete fLoGainSignal;
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}
120
121// --------------------------------------------------------------------------
122//
123// The PreProcess searches for the following input containers:
124// - MRawEvtData
125// - MRawRunHeader
126// - MPedestalCam
127//
128// The following output containers are also searched and created if
129// they were not found:
130//
131// - MArrivalTimeCam
132//
133// If the flag fExtractCharges is set, also following containers are searched:
134//
135// - MExtractedSignalCam
136//
137Int_t MExtractTimeAndChargeSpline::PreProcess(MParList *pList)
138{
139
140 if (!MExtractTimeAndCharge::PreProcess(pList))
141 return kFALSE;
142
143 return kTRUE;
144}
145
146// --------------------------------------------------------------------------
147//
148// ReInit
149//
150// Calls:
151// - MExtractTimeAndCharge::ReInit(pList);
152// - Deletes all arrays, if not NULL
153// - Creates new arrays according to the extraction range
154//
155Bool_t MExtractTimeAndChargeSpline::ReInit(MParList *pList)
156{
157
158 if (!MExtractTimeAndCharge::ReInit(pList))
159 return kFALSE;
160
161 if (fHiGainSignal)
162 delete fHiGainSignal;
163 if (fLoGainSignal)
164 delete fLoGainSignal;
165 if (fHiGainFirstDeriv)
166 delete fHiGainFirstDeriv;
167 if (fLoGainFirstDeriv)
168 delete fLoGainFirstDeriv;
169 if (fHiGainSecondDeriv)
170 delete fHiGainSecondDeriv;
171 if (fLoGainSecondDeriv)
172 delete fLoGainSecondDeriv;
173
174 Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast;
175
176 fHiGainSignal = new Float_t[range];
177 memset(fHiGainSignal,0,range*sizeof(Float_t));
178 fHiGainFirstDeriv = new Float_t[range];
179 memset(fHiGainFirstDeriv,0,range*sizeof(Float_t));
180 fHiGainSecondDeriv = new Float_t[range];
181 memset(fHiGainSecondDeriv,0,range*sizeof(Float_t));
182
183 range = fLoGainLast - fLoGainFirst + 1;
184
185 fLoGainSignal = new Float_t[range];
186 memset(fLoGainSignal,0,range*sizeof(Float_t));
187 fLoGainFirstDeriv = new Float_t[range];
188 memset(fLoGainFirstDeriv,0,range*sizeof(Float_t));
189 fLoGainSecondDeriv = new Float_t[range];
190 memset(fLoGainSecondDeriv,0,range*sizeof(Float_t));
191
192 return kTRUE;
193}
194
195
196// --------------------------------------------------------------------------
197//
198// Calculates the arrival time for each pixel
199//
200void MExtractTimeAndChargeSpline::FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum,
201 Float_t &time, Float_t &dtime,
202 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
203{
204
205 Int_t range = fHiGainLast - fHiGainFirst + 1;
206 const Byte_t *end = first + range;
207 Byte_t *p = first;
208 Int_t count = 0;
209
210 Float_t pedes = ped.GetPedestal();
211 const Float_t ABoffs = ped.GetPedestalABoffset();
212
213 Float_t PedMean[2];
214 PedMean[0] = pedes + ABoffs;
215 PedMean[1] = pedes - ABoffs;
216
217 fAbMax = 0.;
218 fAbMaxPos = 0.;
219 Byte_t maxpos = 0;
220
221 //
222 // Check for saturation in all other slices
223 //
224 while (p<end)
225 {
226
227 const Int_t ids = fHiGainFirst + count ;
228 const Float_t signal = (Float_t)*p - PedMean[(ids+abflag) & 0x1];
229 fHiGainSignal[count] = signal;
230
231 if (signal > fAbMax)
232 {
233 fAbMax = signal;
234 maxpos = p-first;
235 }
236
237 if (*p >= fSaturationLimit)
238 sat++;
239
240 p++;
241 count++;
242 }
243
244 if (fHiLoLast != 0)
245 {
246
247 end = logain + fHiLoLast;
248
249 while (logain<end)
250 {
251
252 const Int_t ids = fHiGainFirst + range ;
253 const Float_t signal = (Float_t)*logain - PedMean[(ids+abflag) & 0x1];
254 fHiGainSignal[range] = signal;
255 range++;
256
257 if (signal > fAbMax)
258 {
259 fAbMax = signal;
260 maxpos = logain-first;
261 }
262
263 if (*logain >= fSaturationLimit)
264 sat++;
265
266 logain++;
267 }
268 }
269
270 //
271 // Allow no saturated slice
272 // and
273 // Don't start if the maxpos is too close to the left limit.
274 //
275 if (sat || maxpos < 2)
276 {
277 time = IsTimeType(kMaximum)
278 ? (Float_t)(fHiGainFirst + maxpos)
279 : (Float_t)(fHiGainFirst + maxpos - 1);
280 return;
281 }
282
283 Float_t pp;
284
285 fHiGainSecondDeriv[0] = 0.;
286 fHiGainFirstDeriv[0] = 0.;
287
288 for (Int_t i=1;i<range-1;i++)
289 {
290 pp = fHiGainSecondDeriv[i-1] + 4.;
291 fHiGainSecondDeriv[i] = -1.0/pp;
292 fHiGainFirstDeriv [i] = fHiGainSignal[i+1] - fHiGainSignal[i] - fHiGainSignal[i] + fHiGainSignal[i-1];
293 fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
294 }
295
296 fHiGainSecondDeriv[range-1] = 0.;
297 for (Int_t k=range-2;k>=0;k--)
298 fHiGainSecondDeriv[k] = (fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k])/6.;
299
300 //
301 // Now find the maximum
302 //
303 Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
304 Float_t lower = (Float_t)maxpos-1.;
305 Float_t upper = (Float_t)maxpos;
306 fAbMaxPos = upper;
307 Float_t x = lower;
308 Float_t y = 0.;
309 Float_t a = 1.;
310 Float_t b = 0.;
311 Int_t klo = maxpos-1;
312 Int_t khi = maxpos;
313
314 //
315 // Search for the maximum, starting in interval maxpos-1 in steps of 0.2 till maxpos-0.2.
316 // If no maximum is found, go to interval maxpos+1.
317 //
318 while ( x < upper - 0.3 )
319 {
320
321 x += step;
322 a -= step;
323 b += step;
324
325 y = a*fHiGainSignal[klo]
326 + b*fHiGainSignal[khi]
327 + (a*a*a-a)*fHiGainSecondDeriv[klo]
328 + (b*b*b-b)*fHiGainSecondDeriv[khi];
329
330 if (y > fAbMax)
331 {
332 fAbMax = y;
333 fAbMaxPos = x;
334 }
335
336 // *fLog << err << x << " " << y << " " << fAbMaxPos<< endl;
337 }
338
339 //
340 // Test the possibility that the absolute maximum has not been found before the
341 // maxpos and search from maxpos to maxpos+1 in steps of 0.2
342 //
343 if (fAbMaxPos > upper-0.1)
344 {
345
346 upper = (Float_t)maxpos+1.;
347 lower = (Float_t)maxpos;
348 x = lower;
349 a = 1.;
350 b = 0.;
351 khi = maxpos+1;
352 klo = maxpos;
353
354 while (x<upper-0.3)
355 {
356
357 x += step;
358 a -= step;
359 b += step;
360
361 y = a*fHiGainSignal[klo]
362 + b*fHiGainSignal[khi]
363 + (a*a*a-a)*fHiGainSecondDeriv[klo]
364 + (b*b*b-b)*fHiGainSecondDeriv[khi];
365
366 if (y > fAbMax)
367 {
368 fAbMax = y;
369 fAbMaxPos = x;
370 }
371 // *fLog << inf << x << " " << y << " " << fAbMaxPos << endl;
372
373 }
374 }
375
376
377 //
378 // Now, the time, abmax and khicont and klocont are set correctly within the previous precision.
379 // Try a better precision.
380 //
381 const Float_t up = fAbMaxPos+step-0.055;
382 const Float_t lo = fAbMaxPos-step+0.055;
383 const Float_t maxpossave = fAbMaxPos;
384
385 x = fAbMaxPos;
386 a = upper - x;
387 b = x - lower;
388
389 step = 0.02; // step size of 42 ps
390
391 while (x<up)
392 {
393
394 x += step;
395 a -= step;
396 b += step;
397
398 y = a*fHiGainSignal[klo]
399 + b*fHiGainSignal[khi]
400 + (a*a*a-a)*fHiGainSecondDeriv[klo]
401 + (b*b*b-b)*fHiGainSecondDeriv[khi];
402
403 if (y > fAbMax)
404 {
405 fAbMax = y;
406 fAbMaxPos = x;
407 }
408 // *fLog << inf << x << " " << y << " " << fAbMaxPos << endl;
409 }
410
411 //
412 // Second, try from time down to time-0.2 in steps of 0.04.
413 //
414 x = maxpossave;
415
416 //
417 // Test the possibility that the absolute maximum has not been found between
418 // maxpos and maxpos+0.02, then we have to look between maxpos-0.02 and maxpos
419 // which requires new setting of klocont and khicont
420 //
421 if (x < klo + 0.02)
422 {
423 klo--;
424 khi--;
425 upper--;
426 lower--;
427 }
428
429 a = upper - x;
430 b = x - lower;
431
432 while (x>lo)
433 {
434
435 x -= step;
436 a += step;
437 b -= step;
438
439 y = a*fHiGainSignal[klo]
440 + b*fHiGainSignal[khi]
441 + (a*a*a-a)*fHiGainSecondDeriv[klo]
442 + (b*b*b-b)*fHiGainSecondDeriv[khi];
443
444 if (y > fAbMax)
445 {
446 fAbMax = y;
447 fAbMaxPos = x;
448 }
449 // *fLog << warn << x << " " << y << " " << fAbMaxPos << endl;
450 }
451
452 if (IsTimeType(kMaximum))
453 {
454 time = (Float_t)fHiGainFirst + fAbMaxPos;
455 dtime = 0.02;
456 }
457 else
458 {
459 fHalfMax = fAbMax/2.;
460
461 //
462 // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
463 // First, find the right FADC slice:
464 //
465 klo = maxpos - 1;
466 while (klo >= 0)
467 {
468 if (fHiGainSignal[klo] < fHalfMax)
469 break;
470 klo--;
471 }
472
473 //
474 // Loop from the beginning of the slice upwards to reach the fHalfMax:
475 // With means of bisection:
476 //
477 x = (Float_t)klo;
478 a = 1.;
479 b = 0.;
480
481 step = 0.5;
482 Bool_t back = kFALSE;
483
484 Int_t maxcnt = 1000;
485 Int_t cnt = 0;
486
487 while (TMath::Abs(y-fHalfMax) > fResolution)
488 {
489
490 if (back)
491 {
492 x -= step;
493 a += step;
494 b -= step;
495 }
496 else
497 {
498 x += step;
499 a -= step;
500 b += step;
501 }
502
503 y = a*fHiGainSignal[klo]
504 + b*fHiGainSignal[khi]
505 + (a*a*a-a)*fHiGainSecondDeriv[klo]
506 + (b*b*b-b)*fHiGainSecondDeriv[khi];
507
508 if (y > fHalfMax)
509 back = kTRUE;
510 else
511 back = kFALSE;
512
513 if (++cnt > maxcnt)
514 {
515 // *fLog << inf << x << " " << y << " " << fHalfMax << endl;
516 break;
517 }
518
519 step /= 2.;
520 }
521
522 time = (Float_t)fHiGainFirst + x;
523 dtime = fResolution;
524 }
525
526 if (IsChargeType(kAmplitude))
527 {
528 sum = fAbMax;
529 return;
530 }
531
532 if (IsChargeType(kIntegral))
533 {
534 //
535 // Now integrate the whole thing!
536 //
537 Int_t startslice = (Int_t)(fAbMaxPos - fRiseTime);
538 Int_t lastslice = (Int_t)(fAbMaxPos + fFallTime);
539
540 if (startslice < 0)
541 {
542 lastslice -= startslice;
543 startslice = 0;
544 }
545
546 Int_t i = startslice;
547 sum = 0.5*fHiGainSignal[i];
548
549 for (i=startslice+1; i<lastslice; i++)
550 sum += fHiGainSignal[i] + 1.5*fHiGainSecondDeriv[i];
551
552 sum += 0.5*fHiGainSignal[lastslice];
553 }
554
555}
556
557
558// --------------------------------------------------------------------------
559//
560// Calculates the arrival time for each pixel
561//
562void MExtractTimeAndChargeSpline::FindTimeAndChargeLoGain(Byte_t *first, Float_t &sum, Float_t &dsum,
563 Float_t &time, Float_t &dtime,
564 Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag)
565{
566
567 Int_t range = fLoGainLast - fLoGainFirst + 1;
568 const Byte_t *end = first + range;
569 Byte_t *p = first;
570 Int_t count = 0;
571
572 Float_t pedes = ped.GetPedestal();
573 const Float_t ABoffs = ped.GetPedestalABoffset();
574
575 Float_t PedMean[2];
576 PedMean[0] = pedes + ABoffs;
577 PedMean[1] = pedes - ABoffs;
578
579 fAbMax = 0.;
580 fAbMaxPos = 0.;
581 Byte_t maxpos = 0;
582
583 //
584 // Check for saturation in all other slices
585 //
586 while (p<end)
587 {
588
589 const Int_t ids = fLoGainFirst + count ;
590 const Float_t signal = (Float_t)*p - PedMean[(ids+abflag) & 0x1];
591 fLoGainSignal[count] = signal;
592
593 if (signal > fAbMax)
594 {
595 fAbMax = signal;
596 maxpos = p-first;
597 }
598
599 if (*p >= fSaturationLimit)
600 sat++;
601
602 p++;
603 count++;
604 }
605
606 //
607 // Allow no saturated slice
608 // and
609 // Don't start if the maxpos is too close to the left limit.
610 //
611 if (sat || maxpos < 2)
612 {
613 time = IsTimeType(kMaximum)
614 ? (Float_t)(fLoGainFirst + maxpos)
615 : (Float_t)(fLoGainFirst + maxpos - 1);
616 return;
617 }
618
619 Float_t pp;
620
621 fLoGainSecondDeriv[0] = 0.;
622 fLoGainFirstDeriv[0] = 0.;
623
624 for (Int_t i=1;i<range-1;i++)
625 {
626 pp = fLoGainSecondDeriv[i-1] + 4.;
627 fLoGainSecondDeriv[i] = -1.0/pp;
628 fLoGainFirstDeriv [i] = fLoGainSignal[i+1] - fLoGainSignal[i] - fLoGainSignal[i] + fLoGainSignal[i-1];
629 fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
630 }
631
632 fLoGainSecondDeriv[range-1] = 0.;
633 for (Int_t k=range-2;k>=0;k--)
634 fLoGainSecondDeriv[k] = (fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k])/6.;
635
636 //
637 // Now find the maximum
638 //
639 Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
640 Float_t lower = (Float_t)maxpos-1.;
641 Float_t upper = (Float_t)maxpos;
642 fAbMaxPos = upper;
643 Float_t x = lower;
644 Float_t y = 0.;
645 Float_t a = 1.;
646 Float_t b = 0.;
647 Int_t klo = maxpos-1;
648 Int_t khi = maxpos;
649
650 //
651 // Search for the maximum, starting in interval maxpos-1 in steps of 0.2 till maxpos-0.2.
652 // If no maximum is found, go to interval maxpos+1.
653 //
654 while ( x < upper - 0.3 )
655 {
656
657 x += step;
658 a -= step;
659 b += step;
660
661 y = a*fLoGainSignal[klo]
662 + b*fLoGainSignal[khi]
663 + (a*a*a-a)*fLoGainSecondDeriv[klo]
664 + (b*b*b-b)*fLoGainSecondDeriv[khi];
665
666 if (y > fAbMax)
667 {
668 fAbMax = y;
669 fAbMaxPos = x;
670 }
671
672 // *fLog << err << x << " " << y << " " << fAbMaxPos<< endl;
673 }
674
675 //
676 // Test the possibility that the absolute maximum has not been found before the
677 // maxpos and search from maxpos to maxpos+1 in steps of 0.2
678 //
679 if (fAbMaxPos > upper-0.1)
680 {
681
682 upper = (Float_t)maxpos+1.;
683 lower = (Float_t)maxpos;
684 x = lower;
685 a = 1.;
686 b = 0.;
687 khi = maxpos+1;
688 klo = maxpos;
689
690 while (x<upper-0.3)
691 {
692
693 x += step;
694 a -= step;
695 b += step;
696
697 y = a*fLoGainSignal[klo]
698 + b*fLoGainSignal[khi]
699 + (a*a*a-a)*fLoGainSecondDeriv[klo]
700 + (b*b*b-b)*fLoGainSecondDeriv[khi];
701
702 if (y > fAbMax)
703 {
704 fAbMax = y;
705 fAbMaxPos = x;
706 }
707 // *fLog << inf << x << " " << y << " " << fAbMaxPos << endl;
708
709 }
710 }
711
712
713 //
714 // Now, the time, abmax and khicont and klocont are set correctly within the previous precision.
715 // Try a better precision.
716 //
717 const Float_t up = fAbMaxPos+step-0.055;
718 const Float_t lo = fAbMaxPos-step+0.055;
719 const Float_t maxpossave = fAbMaxPos;
720
721 x = fAbMaxPos;
722 a = upper - x;
723 b = x - lower;
724
725 step = 0.02; // step size of 42 ps
726
727 while (x<up)
728 {
729
730 x += step;
731 a -= step;
732 b += step;
733
734 y = a*fLoGainSignal[klo]
735 + b*fLoGainSignal[khi]
736 + (a*a*a-a)*fLoGainSecondDeriv[klo]
737 + (b*b*b-b)*fLoGainSecondDeriv[khi];
738
739 if (y > fAbMax)
740 {
741 fAbMax = y;
742 fAbMaxPos = x;
743 }
744 // *fLog << inf << x << " " << y << " " << fAbMaxPos << endl;
745 }
746
747 //
748 // Second, try from time down to time-0.2 in steps of 0.04.
749 //
750 x = maxpossave;
751
752 //
753 // Test the possibility that the absolute maximum has not been found between
754 // maxpos and maxpos+0.02, then we have to look between maxpos-0.02 and maxpos
755 // which requires new setting of klocont and khicont
756 //
757 if (x < klo + 0.02)
758 {
759 klo--;
760 khi--;
761 upper--;
762 lower--;
763 }
764
765 a = upper - x;
766 b = x - lower;
767
768 while (x>lo)
769 {
770
771 x -= step;
772 a += step;
773 b -= step;
774
775 y = a*fLoGainSignal[klo]
776 + b*fLoGainSignal[khi]
777 + (a*a*a-a)*fLoGainSecondDeriv[klo]
778 + (b*b*b-b)*fLoGainSecondDeriv[khi];
779
780 if (y > fAbMax)
781 {
782 fAbMax = y;
783 fAbMaxPos = x;
784 }
785 // *fLog << warn << x << " " << y << " " << fAbMaxPos << endl;
786 }
787
788 if (IsTimeType(kMaximum))
789 {
790 time = (Float_t)fLoGainFirst + fAbMaxPos;
791 dtime = 0.02;
792 }
793 else
794 {
795 fHalfMax = fAbMax/2.;
796
797 //
798 // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
799 // First, find the right FADC slice:
800 //
801 klo = maxpos - 1;
802 while (klo >= 0)
803 {
804 if (fLoGainSignal[klo] < fHalfMax)
805 break;
806 klo--;
807 }
808
809 //
810 // Loop from the beginning of the slice upwards to reach the fHalfMax:
811 // With means of bisection:
812 //
813 x = (Float_t)klo;
814 a = 1.;
815 b = 0.;
816
817 step = 0.5;
818 Bool_t back = kFALSE;
819
820 Int_t maxcnt = 1000;
821 Int_t cnt = 0;
822
823 while (TMath::Abs(y-fHalfMax) > fResolution)
824 {
825
826 if (back)
827 {
828 x -= step;
829 a += step;
830 b -= step;
831 }
832 else
833 {
834 x += step;
835 a -= step;
836 b += step;
837 }
838
839 y = a*fLoGainSignal[klo]
840 + b*fLoGainSignal[khi]
841 + (a*a*a-a)*fLoGainSecondDeriv[klo]
842 + (b*b*b-b)*fLoGainSecondDeriv[khi];
843
844 if (y > fHalfMax)
845 back = kTRUE;
846 else
847 back = kFALSE;
848
849 if (++cnt > maxcnt)
850 {
851 // *fLog << inf << x << " " << y << " " << fHalfMax << endl;
852 break;
853 }
854
855 step /= 2.;
856 }
857
858 time = (Float_t)fLoGainFirst + x;
859 dtime = fResolution;
860 }
861
862 if (IsChargeType(kAmplitude))
863 {
864 sum = fAbMax;
865 return;
866 }
867
868 if (IsChargeType(kIntegral))
869 {
870 //
871 // Now integrate the whole thing!
872 //
873 Int_t startslice = (Int_t)(fAbMaxPos - fRiseTime);
874 Int_t lastslice = (Int_t)(fAbMaxPos + fFallTime);
875
876 if (startslice < 0)
877 {
878 lastslice -= startslice;
879 startslice = 0;
880 }
881
882 Int_t i = startslice;
883 sum = 0.5*fLoGainSignal[i];
884
885 for (i=startslice+1; i<lastslice; i++)
886 sum += fLoGainSignal[i] + 1.5*fLoGainSecondDeriv[i];
887
888 sum += 0.5*fLoGainSignal[lastslice];
889 }
890
891
892}
893
Note: See TracBrowser for help on using the repository browser.