Index: trunk/MagicSoft/Mars/mextralgo/ExtralgoIncl.h =================================================================== --- trunk/MagicSoft/Mars/mextralgo/ExtralgoIncl.h (revision 7942) +++ trunk/MagicSoft/Mars/mextralgo/ExtralgoIncl.h (revision 7942) @@ -0,0 +1,3 @@ +#ifndef __CINT__ + +#endif // __CINT__ Index: trunk/MagicSoft/Mars/mextralgo/ExtralgoLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mextralgo/ExtralgoLinkDef.h (revision 7942) +++ trunk/MagicSoft/Mars/mextralgo/ExtralgoLinkDef.h (revision 7942) @@ -0,0 +1,10 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MExtralgoSpline+; +#pragma link C++ class MExtralgoDigitalFilter+; + +#endif Index: trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.cc =================================================================== --- trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.cc (revision 7942) +++ trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.cc (revision 7942) @@ -0,0 +1,132 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Hendrik Bartko, 09/2004 +! Author(s): Thomas Bretz, 08/2006 +! +! Copyright: MAGIC Software Development, 2000-2006 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MExtralgoDigitalFilter +// +////////////////////////////////////////////////////////////////////////////// +#include "MExtralgoDigitalFilter.h" + +using namespace std; + +Float_t MExtralgoDigitalFilter::ExtractNoise(Int_t iter, Int_t windowsize) const +{ + Double_t sum, dummy; + Eval(sum, dummy, windowsize, 0, iter-fWeightsPerBin/2); + return sum; +} + +void MExtralgoDigitalFilter::Extract(Int_t windowsize, Float_t timeshift) +{ + fSignal = 0; // default is: no pulse found + fTime = -1; // default is: out if range (--> random) + fSignalDev = 0; // default is: valid + fTimeDev = 0; // default is: valid + + // FIXME: How to handle saturation? + + Float_t maxamp = -FLT_MAX; + Float_t maxtime = 0; + Int_t maxp = -1; + + // + // Calculate the sum of the first fWindowSize slices + // + for (Int_t i=0; imaxamp) + { + maxamp = sumamp; + maxtime = sumtime; + maxp = i; + } + } + + // The number of available slices were smaller than the + // extraction window size of the extractor + if (maxp<0) + { + fSignalDev = -1; // means: is invalid + fTimeDev = -1; // means: is invalid + return; + } + + // For some reason (by chance or because all slices contain 0) + // maxamp is 0. This means the signal is zero and no arrival + // time can be extracted (but both informations are valid) + if (maxamp==0) + return; + + // This is the time offset from the extraction position + const Float_t time = maxtime/maxamp; + + // This is used to align the weights to bins between + // -0.5/fWeightsPerBin and 0.5/fWeightsPerBin instead of + // 0 and 1./fWeightsPerBin + const Float_t halfres = 0.5/fWeightsPerBin; + + // move the extractor by an offset number of slices + // determined by the extracted time + const Int_t integ = TMath::FloorNint(time+halfres+0.5); + maxp -= integ; + + // Convert the residual fraction of one slice into an + // offset position in the extraction weights + Int_t frac = TMath::FloorNint((time+halfres-integ)*fWeightsPerBin); + + // Align maxp into available range + if (maxp < 0) + { + maxp = 0; + frac = -fWeightsPerBin/2; // Assume peak at the beginning of the first slice + } + if (maxp+windowsize > fNum) + { + maxp = fNum-windowsize; + frac = fWeightsPerBin/2-1; // Assume peak at the end of the last slice + } + + Double_t sumamp, sumtime; + Eval(sumamp, sumtime, windowsize, maxp, frac); + + fSignal = sumamp; + if (sumamp == 0) + return; + + // here, the first high-gain slice is already included + // in the fTimeShiftHiGain this shifts the time to the + // start of the rising edge + fTime = maxp - Float_t(frac)/fWeightsPerBin + timeshift; + + const Float_t timefineadjust = sumtime/sumamp; + + // FIXME: Is 3 a good value? + if (TMath::Abs(timefineadjust)*fWeightsPerBin < 3.) + fTime -= timefineadjust; +} Index: trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.h =================================================================== --- trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.h (revision 7942) +++ trunk/MagicSoft/Mars/mextralgo/MExtralgoDigitalFilter.h (revision 7942) @@ -0,0 +1,74 @@ +#ifndef MARS_MExtralgoDigitalFilter +#define MARS_MExtralgoDigitalFilter + +#ifndef ROOT_TROOT +#include +#endif + +class MExtralgoDigitalFilter +{ +private: + // Input + Float_t *fVal; + Int_t fNum; + + Float_t *fWeightsAmp; + Float_t *fWeightsTime; + + Int_t fWeightsPerBin; // Number of weights per data bin + + // Result + Float_t fTime; + Float_t fTimeDev; + Float_t fSignal; + Float_t fSignalDev; + + inline void Eval(Double_t &sumamp, Double_t &sumtime, const Int_t windowsize, const Int_t startv, const Int_t startw=0) const + { + // + // Slide with a window of size windowsize over the sample + // and multiply the entries with the corresponding weights + // + sumamp = 0; + sumtime = 0; + + // Shift the start of the weight to the center of sample 0 + Float_t const *wa = fWeightsAmp + fWeightsPerBin/2 + startw; + Float_t const *wt = fWeightsTime + fWeightsPerBin/2 + startw; + Float_t *const beg = fVal+startv; + + for (Float_t const *pex=beg; pex +! Author(s): Markus Gaug 09/2004 +! +! Copyright: MAGIC Software Development, 2002-2006 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MExtralgoSpline +// +// Fast Spline extractor using a cubic spline algorithm, adapted from +// Numerical Recipes in C++, 2nd edition, pp. 116-119. +// +// The coefficients "ya" are here denoted as "fHiGainSignal" and "fLoGainSignal" +// which means the FADC value subtracted by the clock-noise corrected pedestal. +// +// The coefficients "y2a" get immediately divided 6. and are called here +// "fHiGainSecondDeriv" and "fLoGainSecondDeriv" although they are now not exactly +// the second derivative coefficients any more. +// +// The calculation of the cubic-spline interpolated value "y" on a point +// "x" along the FADC-slices axis becomes: +// +// y = a*fHiGainSignal[klo] + b*fHiGainSignal[khi] +// + (a*a*a-a)*fHiGainSecondDeriv[klo] + (b*b*b-b)*fHiGainSecondDeriv[khi] +// +// with: +// a = (khi - x) +// b = (x - klo) +// +// and "klo" being the lower bin edge FADC index and "khi" the upper bin edge FADC index. +// fHiGainSignal[klo] and fHiGainSignal[khi] are the FADC values at "klo" and "khi". +// +// An analogues formula is used for the low-gain values. +// +// The coefficients fHiGainSecondDeriv and fLoGainSecondDeriv are calculated with the +// following simplified algorithm: +// +// for (Int_t i=1;i=0;k--) +// fHiGainSecondDeriv[k] = (fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k])/6.; +// +// +// This algorithm takes advantage of the fact that the x-values are all separated by exactly 1 +// which simplifies the Numerical Recipes algorithm. +// (Note that the variables "fHiGainFirstDeriv" are not real first derivative coefficients.) +// +// The algorithm to search the time proceeds as follows: +// +// 1) Calculate all fHiGainSignal from fHiGainFirst to fHiGainLast +// (note that an "overlap" to the low-gain arrays is possible: i.e. fHiGainLast>14 in the case of +// the MAGIC FADCs). +// 2) Remember the position of the slice with the highest content "fAbMax" at "fAbMaxPos". +// 3) If one or more slices are saturated or fAbMaxPos is less than 2 slices from fHiGainFirst, +// return fAbMaxPos as time and fAbMax as charge (note that the pedestal is subtracted here). +// 4) Calculate all fHiGainSecondDeriv from the fHiGainSignal array +// 5) Search for the maximum, starting in interval fAbMaxPos-1 in steps of 0.2 till fAbMaxPos-0.2. +// If no maximum is found, go to interval fAbMaxPos+1. +// --> 4 function evaluations +// 6) Search for the absolute maximum from fAbMaxPos to fAbMaxPos+1 in steps of 0.2 +// --> 4 function evaluations +// 7) Try a better precision searching from new max. position fAbMaxPos-0.2 to fAbMaxPos+0.2 +// in steps of 0.025 (83 psec. in the case of the MAGIC FADCs). +// --> 14 function evaluations +// 8) If Time Extraction Type kMaximum has been chosen, the position of the found maximum is +// returned, else: +// 9) The Half Maximum is calculated. +// 10) fHiGainSignal is called beginning from fAbMaxPos-1 backwards until a value smaller than fHalfMax +// is found at "klo". +// 11) Then, the spline value between "klo" and "klo"+1 is halfed by means of bisection as long as +// the difference between fHalfMax and spline evaluation is less than fResolution (default: 0.01). +// --> maximum 12 interations. +// +// The algorithm to search the charge proceeds as follows: +// +// 1) If Charge Type: kAmplitude was chosen, return the Maximum of the spline, found during the +// time search. +// 2) If Charge Type: kIntegral was chosen, sum the fHiGainSignal between: +// (Int_t)(fAbMaxPos - fRiseTimeHiGain) and +// (Int_t)(fAbMaxPos + fFallTimeHiGain) +// (default: fRiseTime: 1.5, fFallTime: 4.5) +// sum the fLoGainSignal between: +// (Int_t)(fAbMaxPos - fRiseTimeHiGain*fLoGainStretch) and +// (Int_t)(fAbMaxPos + fFallTimeHiGain*fLoGainStretch) +// (default: fLoGainStretch: 1.5) +// +// The values: fNumHiGainSamples and fNumLoGainSamples are set to: +// 1) If Charge Type: kAmplitude was chosen: 1. +// 2) If Charge Type: kIntegral was chosen: fRiseTimeHiGain + fFallTimeHiGain +// or: fNumHiGainSamples*fLoGainStretch in the case of the low-gain +// +// Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast) +// to modify the ranges. +// +// Defaults: +// fHiGainFirst = 2 +// fHiGainLast = 14 +// fLoGainFirst = 2 +// fLoGainLast = 14 +// +// Call: SetResolution() to define the resolution of the half-maximum search. +// Default: 0.01 +// +// Call: SetRiseTime() and SetFallTime() to define the integration ranges +// for the case, the extraction type kIntegral has been chosen. +// +// Call: - SetChargeType(MExtractTimeAndChargeSpline::kAmplitude) for the +// computation of the amplitude at the maximum (default) and extraction +// the position of the maximum (default) +// --> no further function evaluation needed +// - SetChargeType(MExtractTimeAndChargeSpline::kIntegral) for the +// computation of the integral beneith the spline between fRiseTimeHiGain +// from the position of the maximum to fFallTimeHiGain after the position of +// the maximum. The Low Gain is computed with half a slice more at the rising +// edge and half a slice more at the falling edge. +// The time of the half maximum is returned. +// --> needs one function evaluations but is more precise +// +////////////////////////////////////////////////////////////////////////////// +#include "MExtralgoSpline.h" + +using namespace std; + +void MExtralgoSpline::InitDerivatives() const +{ + fDer1[0] = 0.; + fDer2[0] = 0.; + + for (Int_t i=1; i=0; k--) + fDer2[k] = fDer2[k]*fDer2[k+1] + fDer1[k]; + + for (Int_t k=fNum-2; k>=0; k--) + fDer2[k] /= 6.; +} + +Float_t MExtralgoSpline::CalcIntegral(Float_t start, Float_t range) const +{ + // The number of steps is calculated directly from the integration + // window. This is the only way to ensure we are not dealing with + // numerical rounding uncertanties, because we always get the same + // value under the same conditions -- it might still be different on + // other machines! + const Float_t step = 0.2; + const Float_t width = fRiseTime+fFallTime; + const Float_t max = range-1 - (width+step); + const Int_t num = TMath::Nint(width/step); + + // The order is important. In some cases (limlo-/limup-check) it can + // happen that max<0. In this case we start at 0 + if (start > max) + start = max; + if (start < 0) + start = 0; + + start += step/2; + + Double_t sum = 0.; + for (Int_t i=0; i fNum-TMath::Ceil(fFallTime)-1; + if (sat || limlo || limup) + { + fTimeDev = 1.0; + if (fExtractionType == kAmplitude) + { + fSignal = fVal[maxpos]; + fTime = maxpos; + fSignalDev = 0; // means: is valid + return; + } + + fSignal = CalcIntegral(limlo ? 0 : fNum, fNum); + fTime = maxpos - 1; + fSignalDev = 0; // means: is valid + return; + } + + fTimeDev = fResolution; + + // + // Now find the maximum + // + Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one + + Int_t klo = maxpos-1; + Float_t fAbMaxPos = maxpos;//! Current position of the maximum of the spline + Float_t fAbMax = fVal[maxpos];//! Current maximum of the spline + + // + // Search for the maximum, starting in interval maxpos-1 in steps of 0.2 till maxpos-0.2. + // If no maximum is found, go to interval maxpos+1. + // + for (Int_t i=0; i fAbMax) + { + fAbMax = y; + fAbMaxPos = x; + } + } + + // + // Search for the absolute maximum from maxpos to maxpos+1 in steps of 0.2 + // + if (fAbMaxPos > maxpos - 0.1) + { + klo = maxpos; + + for (Int_t i=0; i fAbMax) + { + fAbMax = y; + fAbMaxPos = x; + } + } + } + + // + // Now, the time, abmax and khicont and klocont are set correctly within the previous precision. + // Try a better precision. + // + const Float_t up = fAbMaxPos+step - 3.0*fResolution; + const Float_t lo = fAbMaxPos-step + 3.0*fResolution; + const Float_t abmaxpos = fAbMaxPos; + + step = 2.*fResolution; // step size of 0.1 FADC slices + + for (int i=0; i fAbMax) + { + fAbMax = y; + fAbMaxPos = x; + } + } + + // + // Second, try from time down to time-0.2 in steps of fResolution. + // + + // + // Test the possibility that the absolute maximum has not been found between + // maxpos and maxpos+0.05, then we have to look between maxpos-0.05 and maxpos + // which requires new setting of klocont and khicont + // + if (abmaxpos < klo + fResolution) + klo--; + + for (int i=TMath::Nint(TMath::Ceil((abmaxpos-lo)/step))-1; i>=0; i--) + { + const Float_t x = abmaxpos - (i+1)*step; + const Float_t y = Eval(x, klo); + if (y > fAbMax) + { + fAbMax = y; + fAbMaxPos = x; + } + } + + if (fExtractionType == kAmplitude) + { + fTime = fAbMaxPos; + fSignal = fAbMax; + fSignalDev = 0; // means: is valid + return; + } + + Float_t fHalfMax = fAbMax/2.;//! Current half maximum of the spline + + // + // Now, loop from the maximum bin leftward down in order to find the + // position of the half maximum. First, find the right FADC slice: + // + klo = maxpos; + while (klo > 0) + { + if (fVal[--klo] < fHalfMax) + break; + } + + // + // Loop from the beginning of the slice upwards to reach the fHalfMax: + // With means of bisection: + // + step = 0.5; + + Int_t maxcnt = 20; + Int_t cnt = 0; + + Float_t y = Eval(klo, klo); // FIXME: IS THIS CORRECT??????? + Float_t x = klo; + Bool_t back = kFALSE; + + while (TMath::Abs(y-fHalfMax) > fResolution) + { + x += back ? -step : +step; + + const Float_t y = Eval(x, klo); + + back = y > fHalfMax; + + if (++cnt > maxcnt) + break; + + step /= 2.; + } + + // + // Now integrate the whole thing! + // + fTime = x; + fSignal = CalcIntegral(fAbMaxPos - fRiseTime, fNum); + fSignalDev = 0; // means: is valid +} Index: trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h =================================================================== --- trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h (revision 7942) +++ trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h (revision 7942) @@ -0,0 +1,75 @@ +#ifndef MARS_MExtralgoSpline +#define MARS_MExtralgoSpline + +#ifndef ROOT_TROOT +#include +#endif + +class MExtralgoSpline +{ +public: + enum ExtractionType_t { kAmplitude, kIntegral }; //! Possible time and charge extraction types + +private: + ExtractionType_t fExtractionType; + +private: + //Bool_t fIsOwner; // Owner of derivatives.... + + // Input + Float_t *fVal; + Int_t fNum; + + Float_t *fDer1; + Float_t *fDer2; + + Float_t fRiseTime; + Float_t fFallTime; + + Float_t fResolution; + + // Result + Float_t fTime; + Float_t fTimeDev; + Float_t fSignal; + Float_t fSignalDev; + + inline Float_t Eval(Float_t val, Float_t a, Float_t deriv) const + { + return a*val + (a*a*a-a)*deriv; + } + + inline Float_t Eval(const Float_t x, const Int_t i) const + { + const Float_t b = x-i; + return Eval(fVal[i], 1-b, fDer2[i]) + Eval(fVal[i+1], b, fDer2[i+1]); + } + + void InitDerivatives() const; + Float_t CalcIntegral(Float_t start, Float_t range) const; + +public: + MExtralgoSpline(Float_t *val, Int_t n, Float_t *der1, Float_t *der2) + : fExtractionType(kIntegral), fVal(val), fNum(n), fDer1(der1), fDer2(der2), fTime(0), fTimeDev(-1), fSignal(0), fSignalDev(-1) + { + InitDerivatives(); + } + + void SetRiseFallTime(Float_t rise, Float_t fall) { fRiseTime=rise; fFallTime=fall; } + void SetExtrationType(ExtractionType_t typ) { fExtractionType = typ; } + void SetResolution(Float_t res) { fResolution=res; } + + Float_t GetTime() const { return fTime; } + Float_t GetSignal() const { return fSignal; } + + Float_t GetTimeDev() const { return fTimeDev; } + Float_t GetSignalDev() const { return fSignalDev; } + + void GetSignal(Float_t &sig, Float_t &dsig) const { sig=fSignal; dsig=fSignalDev; } + void GetTime(Float_t &sig, Float_t &dsig) const { sig=fTime; dsig=fTimeDev; } + + Float_t ExtractNoise(Int_t iter); + void Extract(Byte_t sat, Int_t maxpos); +}; + +#endif Index: trunk/MagicSoft/Mars/mextralgo/Makefile =================================================================== --- trunk/MagicSoft/Mars/mextralgo/Makefile (revision 7942) +++ trunk/MagicSoft/Mars/mextralgo/Makefile (revision 7942) @@ -0,0 +1,32 @@ +################################################################## +# +# makefile +# +# for the MARS software +# +################################################################## +include ../Makefile.conf.$(OSTYPE) +include ../Makefile.conf.general + +#------------------------------------------------------------------------------ + +# +# Handling name of the Root Dictionary Files +# +CINT = Extralgo + +# +# connect the include files defined in the config.mk file +# +INCLUDES = -I. + +SRCFILES = MExtralgoSpline.cc \ + MExtralgoDigitalFilter.cc + +############################################################ + +all: $(OBJS) + +include ../Makefile.rules + +mrproper: clean rmbak