| 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-2006
|
|---|
| 21 | !
|
|---|
| 22 | !
|
|---|
| 23 | \* ======================================================================== */
|
|---|
| 24 |
|
|---|
| 25 | //////////////////////////////////////////////////////////////////////////////
|
|---|
| 26 | //
|
|---|
| 27 | // MExtralgoSpline
|
|---|
| 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 "MExtralgoSpline.h"
|
|---|
| 144 |
|
|---|
| 145 | using namespace std;
|
|---|
| 146 |
|
|---|
| 147 | void MExtralgoSpline::InitDerivatives() const
|
|---|
| 148 | {
|
|---|
| 149 | fDer1[0] = 0.;
|
|---|
| 150 | fDer2[0] = 0.;
|
|---|
| 151 |
|
|---|
| 152 | for (Int_t i=1; i<fNum-1; i++)
|
|---|
| 153 | {
|
|---|
| 154 | const Float_t pp = fDer2[i-1] + 4.;
|
|---|
| 155 |
|
|---|
| 156 | fDer2[i] = -1.0/pp;
|
|---|
| 157 |
|
|---|
| 158 | const Float_t d1 = fVal[i+1] - 2*fVal[i] + fVal[i-1];
|
|---|
| 159 | fDer1[i] = (6.0*d1-fDer1[i-1])/pp;
|
|---|
| 160 | }
|
|---|
| 161 |
|
|---|
| 162 | fDer2[fNum-1] = 0.;
|
|---|
| 163 |
|
|---|
| 164 | for (Int_t k=fNum-2; k>=0; k--)
|
|---|
| 165 | fDer2[k] = fDer2[k]*fDer2[k+1] + fDer1[k];
|
|---|
| 166 |
|
|---|
| 167 | for (Int_t k=fNum-2; k>=0; k--)
|
|---|
| 168 | fDer2[k] /= 6.;
|
|---|
| 169 | }
|
|---|
| 170 |
|
|---|
| 171 | Float_t MExtralgoSpline::CalcIntegral(Float_t start, Float_t range) const
|
|---|
| 172 | {
|
|---|
| 173 | // The number of steps is calculated directly from the integration
|
|---|
| 174 | // window. This is the only way to ensure we are not dealing with
|
|---|
| 175 | // numerical rounding uncertanties, because we always get the same
|
|---|
| 176 | // value under the same conditions -- it might still be different on
|
|---|
| 177 | // other machines!
|
|---|
| 178 | const Float_t step = 0.2;
|
|---|
| 179 | const Float_t width = fRiseTime+fFallTime;
|
|---|
| 180 | const Float_t max = range-1 - (width+step);
|
|---|
| 181 | const Int_t num = TMath::Nint(width/step);
|
|---|
| 182 |
|
|---|
| 183 | // The order is important. In some cases (limlo-/limup-check) it can
|
|---|
| 184 | // happen that max<0. In this case we start at 0
|
|---|
| 185 | if (start > max)
|
|---|
| 186 | start = max;
|
|---|
| 187 | if (start < 0)
|
|---|
| 188 | start = 0;
|
|---|
| 189 |
|
|---|
| 190 | start += step/2;
|
|---|
| 191 |
|
|---|
| 192 | Double_t sum = 0.;
|
|---|
| 193 | for (Int_t i=0; i<num; i++)
|
|---|
| 194 | {
|
|---|
| 195 | const Float_t x = start+i*step;
|
|---|
| 196 | const Int_t klo = (Int_t)TMath::Floor(x);
|
|---|
| 197 | // Note: if x is close to one integer number (= a FADC sample)
|
|---|
| 198 | // we get the same result by using that sample as klo, and the
|
|---|
| 199 | // next one as khi, or using the sample as khi and the previous
|
|---|
| 200 | // one as klo (the spline is of course continuous). So we do not
|
|---|
| 201 | // expect problems from rounding issues in the argument of
|
|---|
| 202 | // Floor() above (we have noticed differences in roundings
|
|---|
| 203 | // depending on the compilation options).
|
|---|
| 204 |
|
|---|
| 205 | sum += Eval(x, klo);
|
|---|
| 206 |
|
|---|
| 207 | // FIXME? Perhaps the integral should be done analitically
|
|---|
| 208 | // between every two FADC slices, instead of numerically
|
|---|
| 209 | }
|
|---|
| 210 | sum *= step; // Transform sum in integral
|
|---|
| 211 | return sum;
|
|---|
| 212 | }
|
|---|
| 213 |
|
|---|
| 214 | Float_t MExtralgoSpline::ExtractNoise(Int_t iter)
|
|---|
| 215 | {
|
|---|
| 216 | const Float_t nsx = iter * fResolution;
|
|---|
| 217 |
|
|---|
| 218 | if (fExtractionType == kAmplitude)
|
|---|
| 219 | return Eval(nsx+1, 1);
|
|---|
| 220 | else
|
|---|
| 221 | return CalcIntegral(2. + nsx, fNum);
|
|---|
| 222 | }
|
|---|
| 223 |
|
|---|
| 224 | void MExtralgoSpline::Extract(Byte_t sat, Int_t maxpos)
|
|---|
| 225 | {
|
|---|
| 226 | fSignal = 0;
|
|---|
| 227 | fTime = 0;
|
|---|
| 228 | fSignalDev = -1;
|
|---|
| 229 | fTimeDev = -1;
|
|---|
| 230 |
|
|---|
| 231 | //
|
|---|
| 232 | // Allow no saturated slice and
|
|---|
| 233 | // Don't start if the maxpos is too close to the limits.
|
|---|
| 234 | //
|
|---|
| 235 | const Bool_t limlo = maxpos < TMath::Ceil(fRiseTime);
|
|---|
| 236 | const Bool_t limup = maxpos > fNum-TMath::Ceil(fFallTime)-1;
|
|---|
| 237 | if (sat || limlo || limup)
|
|---|
| 238 | {
|
|---|
| 239 | fTimeDev = 1.0;
|
|---|
| 240 | if (fExtractionType == kAmplitude)
|
|---|
| 241 | {
|
|---|
| 242 | fSignal = fVal[maxpos];
|
|---|
| 243 | fTime = maxpos;
|
|---|
| 244 | fSignalDev = 0; // means: is valid
|
|---|
| 245 | return;
|
|---|
| 246 | }
|
|---|
| 247 |
|
|---|
| 248 | fSignal = CalcIntegral(limlo ? 0 : fNum, fNum);
|
|---|
| 249 | fTime = maxpos - 1;
|
|---|
| 250 | fSignalDev = 0; // means: is valid
|
|---|
| 251 | return;
|
|---|
| 252 | }
|
|---|
| 253 |
|
|---|
| 254 | fTimeDev = fResolution;
|
|---|
| 255 |
|
|---|
| 256 | //
|
|---|
| 257 | // Now find the maximum
|
|---|
| 258 | //
|
|---|
| 259 | Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one
|
|---|
| 260 |
|
|---|
| 261 | Int_t klo = maxpos-1;
|
|---|
| 262 | Float_t fAbMaxPos = maxpos;//! Current position of the maximum of the spline
|
|---|
| 263 | Float_t fAbMax = fVal[maxpos];//! Current maximum of the spline
|
|---|
| 264 |
|
|---|
| 265 | //
|
|---|
| 266 | // Search for the maximum, starting in interval maxpos-1 in steps of 0.2 till maxpos-0.2.
|
|---|
| 267 | // If no maximum is found, go to interval maxpos+1.
|
|---|
| 268 | //
|
|---|
| 269 | for (Int_t i=0; i<TMath::Nint(TMath::Ceil((1-0.3)/step)); i++)
|
|---|
| 270 | {
|
|---|
| 271 | const Float_t x = klo + step*(i+1);
|
|---|
| 272 | const Float_t y = Eval(x, klo);
|
|---|
| 273 | if (y > fAbMax)
|
|---|
| 274 | {
|
|---|
| 275 | fAbMax = y;
|
|---|
| 276 | fAbMaxPos = x;
|
|---|
| 277 | }
|
|---|
| 278 | }
|
|---|
| 279 |
|
|---|
| 280 | //
|
|---|
| 281 | // Search for the absolute maximum from maxpos to maxpos+1 in steps of 0.2
|
|---|
| 282 | //
|
|---|
| 283 | if (fAbMaxPos > maxpos - 0.1)
|
|---|
| 284 | {
|
|---|
| 285 | klo = maxpos;
|
|---|
| 286 |
|
|---|
| 287 | for (Int_t i=0; i<TMath::Nint(TMath::Ceil((1-0.3)/step)); i++)
|
|---|
| 288 | {
|
|---|
| 289 | const Float_t x = klo + step*(i+1);
|
|---|
| 290 | const Float_t y = Eval(x, klo);
|
|---|
| 291 | if (y > fAbMax)
|
|---|
| 292 | {
|
|---|
| 293 | fAbMax = y;
|
|---|
| 294 | fAbMaxPos = x;
|
|---|
| 295 | }
|
|---|
| 296 | }
|
|---|
| 297 | }
|
|---|
| 298 |
|
|---|
| 299 | //
|
|---|
| 300 | // Now, the time, abmax and khicont and klocont are set correctly within the previous precision.
|
|---|
| 301 | // Try a better precision.
|
|---|
| 302 | //
|
|---|
| 303 | const Float_t up = fAbMaxPos+step - 3.0*fResolution;
|
|---|
| 304 | const Float_t lo = fAbMaxPos-step + 3.0*fResolution;
|
|---|
| 305 | const Float_t abmaxpos = fAbMaxPos;
|
|---|
| 306 |
|
|---|
| 307 | step = 2.*fResolution; // step size of 0.1 FADC slices
|
|---|
| 308 |
|
|---|
| 309 | for (int i=0; i<TMath::Nint(TMath::Ceil((up-abmaxpos)/step)); i++)
|
|---|
| 310 | {
|
|---|
| 311 | const Float_t x = abmaxpos + (i+1)*step;
|
|---|
| 312 | const Float_t y = Eval(x, klo);
|
|---|
| 313 | if (y > fAbMax)
|
|---|
| 314 | {
|
|---|
| 315 | fAbMax = y;
|
|---|
| 316 | fAbMaxPos = x;
|
|---|
| 317 | }
|
|---|
| 318 | }
|
|---|
| 319 |
|
|---|
| 320 | //
|
|---|
| 321 | // Second, try from time down to time-0.2 in steps of fResolution.
|
|---|
| 322 | //
|
|---|
| 323 |
|
|---|
| 324 | //
|
|---|
| 325 | // Test the possibility that the absolute maximum has not been found between
|
|---|
| 326 | // maxpos and maxpos+0.05, then we have to look between maxpos-0.05 and maxpos
|
|---|
| 327 | // which requires new setting of klocont and khicont
|
|---|
| 328 | //
|
|---|
| 329 | if (abmaxpos < klo + fResolution)
|
|---|
| 330 | klo--;
|
|---|
| 331 |
|
|---|
| 332 | for (int i=TMath::Nint(TMath::Ceil((abmaxpos-lo)/step))-1; i>=0; i--)
|
|---|
| 333 | {
|
|---|
| 334 | const Float_t x = abmaxpos - (i+1)*step;
|
|---|
| 335 | const Float_t y = Eval(x, klo);
|
|---|
| 336 | if (y > fAbMax)
|
|---|
| 337 | {
|
|---|
| 338 | fAbMax = y;
|
|---|
| 339 | fAbMaxPos = x;
|
|---|
| 340 | }
|
|---|
| 341 | }
|
|---|
| 342 |
|
|---|
| 343 | if (fExtractionType == kAmplitude)
|
|---|
| 344 | {
|
|---|
| 345 | fTime = fAbMaxPos;
|
|---|
| 346 | fSignal = fAbMax;
|
|---|
| 347 | fSignalDev = 0; // means: is valid
|
|---|
| 348 | return;
|
|---|
| 349 | }
|
|---|
| 350 |
|
|---|
| 351 | Float_t fHalfMax = fAbMax/2.;//! Current half maximum of the spline
|
|---|
| 352 |
|
|---|
| 353 | //
|
|---|
| 354 | // Now, loop from the maximum bin leftward down in order to find the
|
|---|
| 355 | // position of the half maximum. First, find the right FADC slice:
|
|---|
| 356 | //
|
|---|
| 357 | klo = maxpos;
|
|---|
| 358 | while (klo > 0)
|
|---|
| 359 | {
|
|---|
| 360 | if (fVal[--klo] < fHalfMax)
|
|---|
| 361 | break;
|
|---|
| 362 | }
|
|---|
| 363 |
|
|---|
| 364 | //
|
|---|
| 365 | // Loop from the beginning of the slice upwards to reach the fHalfMax:
|
|---|
| 366 | // With means of bisection:
|
|---|
| 367 | //
|
|---|
| 368 | step = 0.5;
|
|---|
| 369 |
|
|---|
| 370 | Int_t maxcnt = 20;
|
|---|
| 371 | Int_t cnt = 0;
|
|---|
| 372 |
|
|---|
| 373 | Float_t y = Eval(klo, klo); // FIXME: IS THIS CORRECT???????
|
|---|
| 374 | Float_t x = klo;
|
|---|
| 375 | Bool_t back = kFALSE;
|
|---|
| 376 |
|
|---|
| 377 | while (TMath::Abs(y-fHalfMax) > fResolution)
|
|---|
| 378 | {
|
|---|
| 379 | x += back ? -step : +step;
|
|---|
| 380 |
|
|---|
| 381 | const Float_t y = Eval(x, klo);
|
|---|
| 382 |
|
|---|
| 383 | back = y > fHalfMax;
|
|---|
| 384 |
|
|---|
| 385 | if (++cnt > maxcnt)
|
|---|
| 386 | break;
|
|---|
| 387 |
|
|---|
| 388 | step /= 2.;
|
|---|
| 389 | }
|
|---|
| 390 |
|
|---|
| 391 | //
|
|---|
| 392 | // Now integrate the whole thing!
|
|---|
| 393 | //
|
|---|
| 394 | fTime = x;
|
|---|
| 395 | fSignal = CalcIntegral(fAbMaxPos - fRiseTime, fNum);
|
|---|
| 396 | fSignalDev = 0; // means: is valid
|
|---|
| 397 | }
|
|---|
