source: branches/removing_cpp11_features/msignal/MExtractTimeAndChargeSpline.cc@ 19240

Last change on this file since 19240 was 13003, checked in by tbretz, 14 years ago
Improved the possible combinations of how to extract the leading edge or position of maximum and the intgeral or amplitude.
File size: 15.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): Thomas Bretz <mailto:tbretz@astro.uni-wuerzbrug.de>
18! Author(s): Markus Gaug 09/2004 <mailto:markus@ifae.es>
19!
20! Copyright: MAGIC Software Development, 2002-2007
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MExtractTimeAndChargeSpline
28//
29// Fast Spline extractor using a cubic spline algorithm, adapted from
30// Numerical Recipes in C++, 2nd edition, pp. 116-119.
31//
32// The coefficients "ya" are here denoted as "fHiGainSignal" and "fLoGainSignal"
33// which means the FADC value subtracted by the clock-noise corrected pedestal.
34//
35// The coefficients "y2a" get immediately divided 6. and are called here
36// "fHiGainSecondDeriv" and "fLoGainSecondDeriv" although they are now not exactly
37// the second derivative coefficients any more.
38//
39// The calculation of the cubic-spline interpolated value "y" on a point
40// "x" along the FADC-slices axis becomes:
41//
42// y = a*fHiGainSignal[klo] + b*fHiGainSignal[khi]
43// + (a*a*a-a)*fHiGainSecondDeriv[klo] + (b*b*b-b)*fHiGainSecondDeriv[khi]
44//
45// with:
46// a = (khi - x)
47// b = (x - klo)
48//
49// and "klo" being the lower bin edge FADC index and "khi" the upper bin edge FADC index.
50// fHiGainSignal[klo] and fHiGainSignal[khi] are the FADC values at "klo" and "khi".
51//
52// An analogues formula is used for the low-gain values.
53//
54// The coefficients fHiGainSecondDeriv and fLoGainSecondDeriv are calculated with the
55// following simplified algorithm:
56//
57// for (Int_t i=1;i<range-1;i++) {
58// pp = fHiGainSecondDeriv[i-1] + 4.;
59// fHiGainFirstDeriv[i] = fHiGainSignal[i+1] - 2.*fHiGainSignal[i] + fHiGainSignal[i-1]
60// fHiGainFirstDeriv[i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
61// }
62//
63// for (Int_t k=range-2;k>=0;k--)
64// fHiGainSecondDeriv[k] = (fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k])/6.;
65//
66//
67// This algorithm takes advantage of the fact that the x-values are all separated by exactly 1
68// which simplifies the Numerical Recipes algorithm.
69// (Note that the variables "fHiGainFirstDeriv" are not real first derivative coefficients.)
70//
71// The algorithm to search the time proceeds as follows:
72//
73// 1) Calculate all fHiGainSignal from fHiGainFirst to fHiGainLast
74// (note that an "overlap" to the low-gain arrays is possible: i.e. fHiGainLast>14 in the case of
75// the MAGIC FADCs).
76// 2) Remember the position of the slice with the highest content "fAbMax" at "fAbMaxPos".
77// 3) If one or more slices are saturated or fAbMaxPos is less than 2 slices from fHiGainFirst,
78// return fAbMaxPos as time and fAbMax as charge (note that the pedestal is subtracted here).
79// 4) Calculate all fHiGainSecondDeriv from the fHiGainSignal array
80// 5) Search for the maximum, starting in interval fAbMaxPos-1 in steps of 0.2 till fAbMaxPos-0.2.
81// If no maximum is found, go to interval fAbMaxPos+1.
82// --> 4 function evaluations
83// 6) Search for the absolute maximum from fAbMaxPos to fAbMaxPos+1 in steps of 0.2
84// --> 4 function evaluations
85// 7) Try a better precision searching from new max. position fAbMaxPos-0.2 to fAbMaxPos+0.2
86// in steps of 0.025 (83 psec. in the case of the MAGIC FADCs).
87// --> 14 function evaluations
88// 8) If Time Extraction Type kMaximum has been chosen, the position of the found maximum is
89// returned, else:
90// 9) The Half Maximum is calculated.
91// 10) fHiGainSignal is called beginning from fAbMaxPos-1 backwards until a value smaller than fHalfMax
92// is found at "klo".
93// 11) Then, the spline value between "klo" and "klo"+1 is halfed by means of bisection as long as
94// the difference between fHalfMax and spline evaluation is less than fResolution (default: 0.01).
95// --> maximum 12 interations.
96//
97// The algorithm to search the charge proceeds as follows:
98//
99// 1) If Charge Type: kAmplitude was chosen, return the Maximum of the spline, found during the
100// time search.
101// 2) If Charge Type: kIntegral was chosen, sum the fHiGainSignal between:
102// (Int_t)(fAbMaxPos - fRiseTimeHiGain) and
103// (Int_t)(fAbMaxPos + fFallTimeHiGain)
104// (default: fRiseTime: 1.5, fFallTime: 4.5)
105// sum the fLoGainSignal between:
106// (Int_t)(fAbMaxPos - fRiseTimeHiGain*fLoGainStretch) and
107// (Int_t)(fAbMaxPos + fFallTimeHiGain*fLoGainStretch)
108// (default: fLoGainStretch: 1.5)
109//
110// The values: fNumHiGainSamples and fNumLoGainSamples are set to:
111// 1) If Charge Type: kAmplitude was chosen: 1.
112// 2) If Charge Type: kIntegral was chosen: fRiseTimeHiGain + fFallTimeHiGain
113// or: fNumHiGainSamples*fLoGainStretch in the case of the low-gain
114//
115// Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast)
116// to modify the ranges.
117//
118// Defaults:
119// fHiGainFirst = 2
120// fHiGainLast = 14
121// fLoGainFirst = 2
122// fLoGainLast = 14
123//
124// Call: SetResolution() to define the resolution of the half-maximum search.
125// Default: 0.01
126//
127// Call: SetRiseTime() and SetFallTime() to define the integration ranges
128// for the case, the extraction type kIntegral has been chosen.
129//
130// Call: - SetChargeType(MExtractTimeAndChargeSpline::kAmplitude) for the
131// computation of the amplitude at the maximum (default) and extraction
132// the position of the maximum (default)
133// --> no further function evaluation needed
134// - SetChargeType(MExtractTimeAndChargeSpline::kIntegral) for the
135// computation of the integral beneith the spline between fRiseTimeHiGain
136// from the position of the maximum to fFallTimeHiGain after the position of
137// the maximum. The Low Gain is computed with half a slice more at the rising
138// edge and half a slice more at the falling edge.
139// The time of the half maximum is returned.
140// --> needs one function evaluations but is more precise
141//
142//////////////////////////////////////////////////////////////////////////////
143#include "MExtractTimeAndChargeSpline.h"
144
145#include "MPedestalPix.h"
146
147#include "MLog.h"
148#include "MLogManip.h"
149
150ClassImp(MExtractTimeAndChargeSpline);
151
152using namespace std;
153
154const Byte_t MExtractTimeAndChargeSpline::fgHiGainFirst = 0;
155const Byte_t MExtractTimeAndChargeSpline::fgHiGainLast = 14;
156const Int_t MExtractTimeAndChargeSpline::fgLoGainFirst = 1;
157const Byte_t MExtractTimeAndChargeSpline::fgLoGainLast = 14;
158const Float_t MExtractTimeAndChargeSpline::fgResolution = 0.05;
159const Float_t MExtractTimeAndChargeSpline::fgRiseTimeHiGain = 0.64;
160const Float_t MExtractTimeAndChargeSpline::fgFallTimeHiGain = 0.76;
161const Float_t MExtractTimeAndChargeSpline::fgLoGainStretch = 1.5;
162const Float_t MExtractTimeAndChargeSpline::fgOffsetLoGain = 1.3;
163
164// --------------------------------------------------------------------------
165//
166// Default constructor.
167//
168// Calls:
169// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
170//
171// Initializes:
172// - fResolution to fgResolution
173// - fRiseTimeHiGain to fgRiseTimeHiGain
174// - fFallTimeHiGain to fgFallTimeHiGain
175// - Charge Extraction Type to kAmplitude
176// - fLoGainStretch to fgLoGainStretch
177//
178MExtractTimeAndChargeSpline::MExtractTimeAndChargeSpline(const char *name, const char *title)
179 : fRiseTimeHiGain(0), fFallTimeHiGain(0), fHeightTm(0.5), fExtractionType(MExtralgoSpline::kIntegralRel)
180{
181
182 fName = name ? name : "MExtractTimeAndChargeSpline";
183 fTitle = title ? title : "Calculate photons arrival time using a fast spline";
184
185 SetResolution();
186 SetLoGainStretch();
187 SetOffsetLoGain(fgOffsetLoGain);
188
189 SetRiseTimeHiGain();
190 SetFallTimeHiGain();
191
192 SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
193}
194
195
196//-------------------------------------------------------------------
197//
198// Set the ranges
199// In order to set the fNum...Samples variables correctly for the case,
200// the integral is computed, have to overwrite this function and make an
201// explicit call to SetChargeType().
202//
203void MExtractTimeAndChargeSpline::SetRange(UShort_t hifirst, UShort_t hilast, Int_t lofirst, Byte_t lolast)
204{
205 MExtractor::SetRange(hifirst, hilast, lofirst, lolast);
206
207 SetChargeType(fExtractionType);
208}
209
210//-------------------------------------------------------------------
211//
212// Set the Charge Extraction type. Possible are:
213// - kAmplitude: Search the value of the spline at the maximum
214// - kIntegral: Integral the spline from fHiGainFirst to fHiGainLast,
215// by counting the edge bins only half and setting the
216// second derivative to zero, there.
217//
218void MExtractTimeAndChargeSpline::SetChargeType(MExtralgoSpline::ExtractionType_t typ)
219{
220 fExtractionType = typ;
221
222 InitArrays(fHiGainFirstDeriv.GetSize());
223
224 switch (fExtractionType)
225 {
226 case MExtralgoSpline::kAmplitude:
227 case MExtralgoSpline::kAmplitudeRel:
228 case MExtralgoSpline::kAmplitudeAbs:
229 SetResolutionPerPheHiGain(0.053);
230 SetResolutionPerPheLoGain(0.016);
231 return;
232
233 case MExtralgoSpline::kIntegralRel:
234 case MExtralgoSpline::kIntegralAbs:
235 switch (fWindowSizeHiGain)
236 {
237 case 1:
238 SetResolutionPerPheHiGain(0.041);
239 break;
240 case 2:
241 SetResolutionPerPheHiGain(0.064);
242 break;
243 case 3:
244 case 4:
245 SetResolutionPerPheHiGain(0.050);
246 break;
247 case 5:
248 case 6:
249 SetResolutionPerPheHiGain(0.030);
250 break;
251 default:
252 *fLog << warn << GetDescriptor() << ": Could not set the high-gain extractor resolution per phe for window size "
253 << fWindowSizeHiGain << "... using default!" << endl;
254 SetResolutionPerPheHiGain(0.050);
255 break;
256 }
257
258 switch (fWindowSizeLoGain)
259 {
260 case 1:
261 case 2:
262 SetResolutionPerPheLoGain(0.005);
263 break;
264 case 3:
265 case 4:
266 SetResolutionPerPheLoGain(0.017);
267 break;
268 case 5:
269 case 6:
270 case 7:
271 SetResolutionPerPheLoGain(0.005);
272 break;
273 case 8:
274 case 9:
275 SetResolutionPerPheLoGain(0.005);
276 break;
277 default:
278 *fLog << warn << "Could not set the low-gain extractor resolution per phe for window size "
279 << fWindowSizeLoGain << "... using default!" << endl;
280 SetResolutionPerPheLoGain(0.005);
281 break;
282 }
283 }
284}
285
286// --------------------------------------------------------------------------
287//
288// InitArrays
289//
290// Gets called in the ReInit() and initialized the arrays
291//
292Bool_t MExtractTimeAndChargeSpline::InitArrays(Int_t n)
293{
294 // Initialize arrays to the maximum number of entries necessary
295 fHiGainFirstDeriv .Set(n);
296 fHiGainSecondDeriv.Set(n);
297
298 fLoGainFirstDeriv .Set(n);
299 fLoGainSecondDeriv.Set(n);
300
301 fRiseTimeLoGain = fRiseTimeHiGain * fLoGainStretch;
302 fFallTimeLoGain = fFallTimeHiGain * fLoGainStretch;
303
304 switch (fExtractionType)
305 {
306 case MExtralgoSpline::kAmplitude:
307 case MExtralgoSpline::kAmplitudeRel:
308 case MExtralgoSpline::kAmplitudeAbs:
309 fNumHiGainSamples = 1.;
310 fNumLoGainSamples = fLoGainLast ? 1. : 0.;
311 fSqrtHiGainSamples = 1.;
312 fSqrtLoGainSamples = 1.;
313 fWindowSizeHiGain = 1;
314 fWindowSizeLoGain = 1;
315 fRiseTimeHiGain = 0.5;
316 break;
317
318 case MExtralgoSpline::kIntegralAbs:
319 case MExtralgoSpline::kIntegralRel:
320 fNumHiGainSamples = fRiseTimeHiGain + fFallTimeHiGain;
321 fNumLoGainSamples = fLoGainLast ? fRiseTimeLoGain + fFallTimeLoGain : 0.;
322 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
323 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
324 fWindowSizeHiGain = TMath::CeilNint(fRiseTimeHiGain + fFallTimeHiGain);
325 fWindowSizeLoGain = TMath::CeilNint(fRiseTimeLoGain + fFallTimeLoGain);
326 break;
327 }
328
329 return kTRUE;
330}
331
332void MExtractTimeAndChargeSpline::FindTimeAndChargeHiGain2(const Float_t *ptr, Int_t num,
333 Float_t &sum, Float_t &dsum,
334 Float_t &time, Float_t &dtime,
335 Byte_t sat, Int_t maxpos) const
336{
337 // Do some handling if maxpos is last slice!
338 MExtralgoSpline s(ptr, num, fHiGainFirstDeriv.GetArray(), fHiGainSecondDeriv.GetArray());
339
340 s.SetExtractionType(fExtractionType);
341 s.SetHeightTm(fHeightTm);
342 s.SetRiseFallTime(fRiseTimeHiGain, fFallTimeHiGain);
343
344 if (IsNoiseCalculation())
345 {
346 sum = s.ExtractNoise();
347 dsum = 1;
348 return;
349 }
350
351 s.Extract(maxpos);
352 s.GetTime(time, dtime);
353 s.GetSignal(sum, dsum);
354}
355
356void MExtractTimeAndChargeSpline::FindTimeAndChargeLoGain2(const Float_t *ptr, Int_t num,
357 Float_t &sum, Float_t &dsum,
358 Float_t &time, Float_t &dtime,
359 Byte_t sat, Int_t maxpos) const
360{
361 MExtralgoSpline s(ptr, num, fLoGainFirstDeriv.GetArray(), fLoGainSecondDeriv.GetArray());
362
363 s.SetExtractionType(fExtractionType);
364 s.SetHeightTm(fHeightTm);
365 s.SetRiseFallTime(fRiseTimeLoGain, fFallTimeLoGain);
366
367 if (IsNoiseCalculation())
368 {
369 sum = s.ExtractNoise();
370 return;
371 }
372
373 s.Extract(maxpos);
374 s.GetTime(time, dtime);
375 s.GetSignal(sum, dsum);
376}
377
378// --------------------------------------------------------------------------
379//
380// In addition to the resources of the base-class MExtractor:
381// Resolution
382// RiseTimeHiGain
383// FallTimeHiGain
384// LoGainStretch
385// ExtractionType: amplitude, integral
386//
387Int_t MExtractTimeAndChargeSpline::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
388{
389
390 Bool_t rc = kFALSE;
391
392 if (IsEnvDefined(env, prefix, "Resolution", print))
393 {
394 SetResolution(GetEnvValue(env, prefix, "Resolution",fResolution));
395 rc = kTRUE;
396 }
397 if (IsEnvDefined(env, prefix, "RiseTimeHiGain", print))
398 {
399 SetRiseTimeHiGain(GetEnvValue(env, prefix, "RiseTimeHiGain", fRiseTimeHiGain));
400 rc = kTRUE;
401 }
402 if (IsEnvDefined(env, prefix, "FallTimeHiGain", print))
403 {
404 SetFallTimeHiGain(GetEnvValue(env, prefix, "FallTimeHiGain", fFallTimeHiGain));
405 rc = kTRUE;
406 }
407 if (IsEnvDefined(env, prefix, "LoGainStretch", print))
408 {
409 SetLoGainStretch(GetEnvValue(env, prefix, "LoGainStretch", fLoGainStretch));
410 rc = kTRUE;
411 }
412 if (IsEnvDefined(env, prefix, "HeightTm", print))
413 {
414 fHeightTm = GetEnvValue(env, prefix, "HeightTm", fHeightTm);
415 rc = kTRUE;
416 }
417
418 if (IsEnvDefined(env, prefix, "ExtractionType", print))
419 {
420 TString type = GetEnvValue(env, prefix, "ExtractionType", "");
421 type.ToLower();
422 type = type.Strip(TString::kBoth);
423 if (type==(TString)"amplitude")
424 SetChargeType(MExtralgoSpline::kAmplitude);
425 if (type==(TString)"integralabsolute")
426 SetChargeType(MExtralgoSpline::kIntegralAbs);
427 if (type==(TString)"integralrelative")
428 SetChargeType(MExtralgoSpline::kIntegralRel);
429 rc=kTRUE;
430 }
431
432 return MExtractTimeAndCharge::ReadEnv(env, prefix, print) ? kTRUE : rc;
433}
Note: See TracBrowser for help on using the repository browser.