source: trunk/Mars/msignal/MExtractTimeAndChargeSpline.cc@ 10394

Last change on this file since 10394 was 10166, checked in by tbretz, 14 years ago
Removed the old obsolete cvs header line.
File size: 15.7 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(Byte_t hifirst, Byte_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 SetResolutionPerPheHiGain(0.053);
228 SetResolutionPerPheLoGain(0.016);
229 return;
230
231 case MExtralgoSpline::kIntegralRel:
232 case MExtralgoSpline::kIntegralAbs:
233 switch (fWindowSizeHiGain)
234 {
235 case 1:
236 SetResolutionPerPheHiGain(0.041);
237 break;
238 case 2:
239 SetResolutionPerPheHiGain(0.064);
240 break;
241 case 3:
242 case 4:
243 SetResolutionPerPheHiGain(0.050);
244 break;
245 case 5:
246 case 6:
247 SetResolutionPerPheHiGain(0.030);
248 break;
249 default:
250 *fLog << warn << GetDescriptor() << ": Could not set the high-gain extractor resolution per phe for window size "
251 << fWindowSizeHiGain << "... using default!" << endl;
252 SetResolutionPerPheHiGain(0.050);
253 break;
254 }
255
256 switch (fWindowSizeLoGain)
257 {
258 case 1:
259 case 2:
260 SetResolutionPerPheLoGain(0.005);
261 break;
262 case 3:
263 case 4:
264 SetResolutionPerPheLoGain(0.017);
265 break;
266 case 5:
267 case 6:
268 case 7:
269 SetResolutionPerPheLoGain(0.005);
270 break;
271 case 8:
272 case 9:
273 SetResolutionPerPheLoGain(0.005);
274 break;
275 default:
276 *fLog << warn << "Could not set the low-gain extractor resolution per phe for window size "
277 << fWindowSizeLoGain << "... using default!" << endl;
278 SetResolutionPerPheLoGain(0.005);
279 break;
280 }
281 }
282}
283
284// --------------------------------------------------------------------------
285//
286// InitArrays
287//
288// Gets called in the ReInit() and initialized the arrays
289//
290Bool_t MExtractTimeAndChargeSpline::InitArrays(Int_t n)
291{
292 // Initialize arrays to the maximum number of entries necessary
293 fHiGainFirstDeriv .Set(n);
294 fHiGainSecondDeriv.Set(n);
295
296 fLoGainFirstDeriv .Set(n);
297 fLoGainSecondDeriv.Set(n);
298
299 fRiseTimeLoGain = fRiseTimeHiGain * fLoGainStretch;
300 fFallTimeLoGain = fFallTimeHiGain * fLoGainStretch;
301
302 switch (fExtractionType)
303 {
304 case MExtralgoSpline::kAmplitude:
305 fNumHiGainSamples = 1.;
306 fNumLoGainSamples = fLoGainLast ? 1. : 0.;
307 fSqrtHiGainSamples = 1.;
308 fSqrtLoGainSamples = 1.;
309 fWindowSizeHiGain = 1;
310 fWindowSizeLoGain = 1;
311 fRiseTimeHiGain = 0.5;
312 break;
313
314 case MExtralgoSpline::kIntegralAbs:
315 case MExtralgoSpline::kIntegralRel:
316 fNumHiGainSamples = fRiseTimeHiGain + fFallTimeHiGain;
317 fNumLoGainSamples = fLoGainLast ? fRiseTimeLoGain + fFallTimeLoGain : 0.;
318 fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
319 fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
320 fWindowSizeHiGain = TMath::CeilNint(fRiseTimeHiGain + fFallTimeHiGain);
321 fWindowSizeLoGain = TMath::CeilNint(fRiseTimeLoGain + fFallTimeLoGain);
322 break;
323 }
324
325 return kTRUE;
326}
327
328void MExtractTimeAndChargeSpline::FindTimeAndChargeHiGain2(const Float_t *ptr, Int_t num,
329 Float_t &sum, Float_t &dsum,
330 Float_t &time, Float_t &dtime,
331 Byte_t sat, Int_t maxpos) const
332{
333 // Do some handling if maxpos is last slice!
334 MExtralgoSpline s(ptr, num, fHiGainFirstDeriv.GetArray(), fHiGainSecondDeriv.GetArray());
335
336 s.SetExtractionType(fExtractionType);
337 s.SetHeightTm(fHeightTm);
338 s.SetRiseFallTime(fRiseTimeHiGain, fFallTimeHiGain);
339
340 if (IsNoiseCalculation())
341 {
342 sum = s.ExtractNoise();
343 return;
344 }
345
346 s.Extract(maxpos);
347 s.GetTime(time, dtime);
348 s.GetSignal(sum, dsum);
349
350}
351
352void MExtractTimeAndChargeSpline::FindTimeAndChargeLoGain2(const Float_t *ptr, Int_t num,
353 Float_t &sum, Float_t &dsum,
354 Float_t &time, Float_t &dtime,
355 Byte_t sat, Int_t maxpos) const
356{
357 MExtralgoSpline s(ptr, num, fLoGainFirstDeriv.GetArray(), fLoGainSecondDeriv.GetArray());
358
359 s.SetExtractionType(fExtractionType);
360 s.SetHeightTm(fHeightTm);
361 s.SetRiseFallTime(fRiseTimeLoGain, fFallTimeLoGain);
362
363 if (IsNoiseCalculation())
364 {
365 sum = s.ExtractNoise();
366 return;
367 }
368
369 s.Extract(maxpos);
370 s.GetTime(time, dtime);
371 s.GetSignal(sum, dsum);
372}
373
374// --------------------------------------------------------------------------
375//
376// In addition to the resources of the base-class MExtractor:
377// Resolution
378// RiseTimeHiGain
379// FallTimeHiGain
380// LoGainStretch
381// ExtractionType: amplitude, integral
382//
383Int_t MExtractTimeAndChargeSpline::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
384{
385
386 Bool_t rc = kFALSE;
387
388 if (IsEnvDefined(env, prefix, "Resolution", print))
389 {
390 SetResolution(GetEnvValue(env, prefix, "Resolution",fResolution));
391 rc = kTRUE;
392 }
393 if (IsEnvDefined(env, prefix, "RiseTimeHiGain", print))
394 {
395 SetRiseTimeHiGain(GetEnvValue(env, prefix, "RiseTimeHiGain", fRiseTimeHiGain));
396 rc = kTRUE;
397 }
398 if (IsEnvDefined(env, prefix, "FallTimeHiGain", print))
399 {
400 SetFallTimeHiGain(GetEnvValue(env, prefix, "FallTimeHiGain", fFallTimeHiGain));
401 rc = kTRUE;
402 }
403 if (IsEnvDefined(env, prefix, "LoGainStretch", print))
404 {
405 SetLoGainStretch(GetEnvValue(env, prefix, "LoGainStretch", fLoGainStretch));
406 rc = kTRUE;
407 }
408 if (IsEnvDefined(env, prefix, "HeightTm", print))
409 {
410 fHeightTm = GetEnvValue(env, prefix, "HeightTm", fHeightTm);
411 rc = kTRUE;
412 }
413
414 if (IsEnvDefined(env, prefix, "ExtractionType", print))
415 {
416 TString type = GetEnvValue(env, prefix, "ExtractionType", "");
417 type.ToLower();
418 type = type.Strip(TString::kBoth);
419 if (type==(TString)"amplitude")
420 SetChargeType(MExtralgoSpline::kAmplitude);
421 if (type==(TString)"integralabsolute")
422 SetChargeType(MExtralgoSpline::kIntegralAbs);
423 if (type==(TString)"integralrelative")
424 SetChargeType(MExtralgoSpline::kIntegralRel);
425 rc=kTRUE;
426 }
427
428 return MExtractTimeAndCharge::ReadEnv(env, prefix, print) ? kTRUE : rc;
429}
Note: See TracBrowser for help on using the repository browser.