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

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