Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 4275) +++ /trunk/MagicSoft/Mars/Changelog (revision 4276) @@ -56,5 +56,16 @@ fixed. (The difference was usually a couple of pixels) - + * msignal/MExtractTime.h + - make members protected instead of private in order to allow + classes to derive from it + + * msignal/MExtractFixedWindow.cc + - added some documentation + + * msignal/MExtractFixedWindowSpline.[h,cc] + * msignal/Makefile + * msignal/SignalLinkDef.h + - new fast spline signal extractor + 2004/06/02: Antonio Stamerra Index: /trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.cc (revision 4276) +++ /trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.cc (revision 4276) @@ -0,0 +1,359 @@ +/* ======================================================================== *\ +! +! * +! * 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 analyzing 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): Markus Gaug 05/2004 +! +! Copyright: MAGIC Software Development, 2002-2004 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MExtractTimeAndChargeSpline +// +// Fast Spline extractor using a cubic spline algorithm of Numerical Recipes. +// It returns the integral below the interpolating spline. +// +// Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast) +// to modify the ranges. Ranges have to be an even number. In case of odd +// ranges, the last slice will be reduced by one. +// +// Defaults are: +// +// fHiGainFirst = fgHiGainFirst = 3 +// fHiGainLast = fgHiGainLast = 14 +// fLoGainFirst = fgLoGainFirst = 3 +// fLoGainLast = fgLoGainLast = 14 +// +////////////////////////////////////////////////////////////////////////////// +#include "MExtractFixedWindowSpline.h" + +#include "MExtractedSignalCam.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MExtractFixedWindowSpline); + +using namespace std; + +const Byte_t MExtractFixedWindowSpline::fgHiGainFirst = 2; +const Byte_t MExtractFixedWindowSpline::fgHiGainLast = 14; +const Byte_t MExtractFixedWindowSpline::fgLoGainFirst = 3; +const Byte_t MExtractFixedWindowSpline::fgLoGainLast = 14; +// -------------------------------------------------------------------------- +// +// Default constructor. +// +// Calls: +// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast) +// +MExtractFixedWindowSpline::MExtractFixedWindowSpline(const char *name, const char *title) + : fHiGainFirstDeriv(NULL), fLoGainFirstDeriv(NULL), + fHiGainSecondDeriv(NULL), fLoGainSecondDeriv(NULL) +{ + + fName = name ? name : "MExtractFixedWindowSpline"; + fTitle = title ? title : "Signal Extractor for a fixed FADC window using a fast spline"; + + SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast); +} + +MExtractFixedWindowSpline::~MExtractFixedWindowSpline() +{ + + if (fHiGainFirstDeriv) + delete fHiGainFirstDeriv; + if (fLoGainFirstDeriv) + delete fLoGainFirstDeriv; + if (fHiGainSecondDeriv) + delete fHiGainSecondDeriv; + if (fLoGainSecondDeriv) + delete fLoGainSecondDeriv; + +} + +// -------------------------------------------------------------------------- +// +// SetRange: +// +// Checks: +// - if the window defined by (fHiGainLast-fHiGainFirst-1) are odd, subtract one +// - if the window defined by (fLoGainLast-fLoGainFirst-1) are odd, subtract one +// - if the Hi Gain window is smaller than 2, set fHiGainLast to fHiGainFirst+1 +// - if the Lo Gain window is smaller than 2, set fLoGainLast to fLoGainFirst+1 +// +// Calls: +// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast); +// +// Sets: +// - fNumHiGainSamples to: (Float_t)(fHiGainLast-fHiGainFirst+1) +// - fNumLoGainSamples to: (Float_t)(fLoGainLast-fLoGainFirst+1) +// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples) +// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples) +// +void MExtractFixedWindowSpline::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast) +{ + + const Byte_t windowhi = hilast-hifirst+1; + const Byte_t windowlo = lolast-lofirst+1; + + const Byte_t whieven = windowhi & ~1; + const Byte_t wloeven = windowlo & ~1; + + if (whieven != windowhi) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i",": Hi Gain window size has to be even, set last slice from " + ,(int)hilast," to ",(int)(hilast-1)) << endl; + hilast -= 1; + } + + if (wloeven != windowlo) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i",": Lo Gain window size has to be even, set last slice from " + ,(int)lolast," to ",(int)(lolast-1)) << endl; + lolast -= 1; + } + + if (whieven<2) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i",": Hi Gain window is smaller than 2 FADC sampes, set last slice from" + ,(int)hilast," to ",(int)(hifirst+1)) << endl; + hilast = hifirst+1; + } + + if (wloeven<2) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i",": Lo Gain window is smaller than 2 FADC sampes, set last slice from" + ,(int)lolast," to ",(int)(lofirst+1)) << endl; + lolast = lofirst+1; + } + + + MExtractor::SetRange(hifirst,hilast,lofirst,lolast); + + // + // Very important: Because the spline interpolates between the slice, + // the number of samples for the pedestal subtraction + // is now 1 less than with e.g. MExtractFixedWindow + // + fNumHiGainSamples = (Float_t)(fHiGainLast-fHiGainFirst); + fNumLoGainSamples = (Float_t)(fLoGainLast-fLoGainFirst); + + fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples); + fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples); + +} + +// -------------------------------------------------------------------------- +// +// ReInit +// +// Calls: +// - MExtractor::ReInit(pList); +// - fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples, +// fLoGainFirst, fLoGainLast, fNumLoGainSamples); +// +// Deletes all arrays, if not NULL +// Creates new arrays according to the extraction range +// +Bool_t MExtractFixedWindowSpline::ReInit(MParList *pList) +{ + + if (!MExtractor::ReInit(pList)) + return kFALSE; + + fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast+fHiLoLast, fNumHiGainSamples, + fLoGainFirst, fLoGainLast, fNumLoGainSamples); + + if (fHiGainFirstDeriv) + delete fHiGainFirstDeriv; + if (fLoGainFirstDeriv) + delete fLoGainFirstDeriv; + if (fHiGainSecondDeriv) + delete fHiGainSecondDeriv; + if (fLoGainSecondDeriv) + delete fLoGainSecondDeriv; + + Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast; + + fHiGainFirstDeriv = new Float_t[range]; + memset(fHiGainFirstDeriv,0,range*sizeof(Float_t)); + fHiGainSecondDeriv = new Float_t[range]; + memset(fHiGainSecondDeriv,0,range*sizeof(Float_t)); + + range = fLoGainLast - fLoGainFirst + 1; + + fLoGainFirstDeriv = new Float_t[range]; + memset(fLoGainFirstDeriv,0,range*sizeof(Float_t)); + fLoGainSecondDeriv = new Float_t[range]; + memset(fLoGainSecondDeriv,0,range*sizeof(Float_t)); + + return kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// FindSignalHiGain: +// +// - Loop from ptr to (ptr+fHiGainLast-fHiGainFirst) +// - Sum up contents of *ptr +// - If *ptr is greater than fSaturationLimit, raise sat by 1 +// +// - If fHiLoLast is not 0, loop also from logain to (logain+fHiLoLast) +// - Sum up contents of logain +// - If *logain is greater than fSaturationLimit, raise sat by 1 +// +void MExtractFixedWindowSpline::FindSignalHiGain(Byte_t *ptr, Byte_t *logain, Int_t &sum, Byte_t &sat) const +{ + + const Byte_t *end = ptr + fHiGainLast - fHiGainFirst; + Int_t range = fHiGainLast - fHiGainFirst + fHiLoLast + 1; + + Float_t pp; + Int_t i = 0; + + fHiGainSecondDeriv[0] = 0.; + fHiGainFirstDeriv[0] = 0.; + + Float_t sumf = *ptr++/2.; + // + // Check for saturation in all other slices + // + while (ptr= fSaturationLimit) + sat++; + } + + if (*ptr++ >= fSaturationLimit) + sat++; + + if (fHiLoLast == 0) + { + sumf += *ptr/2.; + fHiGainSecondDeriv[++i] = 0.; + + } + else + { + sumf += *ptr; + i++; + + pp = fHiGainSecondDeriv[i-1] + 4.; + fHiGainSecondDeriv[i] = -1.0/pp; + fHiGainFirstDeriv [i] = *(logain) - 2.* *(ptr) + *(ptr-1); + fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp; + + end = logain + fHiLoLast - 1; + + while (logain= fSaturationLimit) + sat++; + + } + sumf += *logain/2; + fHiGainSecondDeriv[++i] = 0.; + } + + for (Int_t k=range-2;k>0;k--) + { + fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k]; + sumf += 0.25*fHiGainSecondDeriv[k]; + } + + sum = (Int_t)sumf; +} + +// -------------------------------------------------------------------------- +// +// FindSignalLoGain: +// +// - Loop from ptr to (ptr+fLoGainLast-fLoGainFirst) +// - Sum up contents of *ptr +// - If *ptr is greater than fSaturationLimit, raise sat by 1 +// +void MExtractFixedWindowSpline::FindSignalLoGain(Byte_t *ptr, Int_t &sum, Byte_t &sat) const +{ + + const Byte_t *end = ptr + fLoGainLast - fLoGainFirst; + Int_t range = fLoGainLast - fLoGainFirst + 1; + + Float_t pp; + Int_t i = 0; + + fLoGainSecondDeriv[0] = 0.; + fLoGainFirstDeriv[0] = 0.; + + Float_t sumf = *ptr++/2.; + // + // Check for saturation in all other slices + // + while (ptr= fSaturationLimit) + sat++; + } + + if (*ptr++ >= fSaturationLimit) + sat++; + + sumf += *ptr/2.; + fLoGainSecondDeriv[++i] = 0.; + + for (Int_t k=range-2;k>0;k--) + { + fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k]; + sumf += 0.25*fLoGainSecondDeriv[k]; + } + + sum = (Int_t)sumf; +} + Index: /trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.h (revision 4276) +++ /trunk/MagicSoft/Mars/msignal/MExtractFixedWindowSpline.h (revision 4276) @@ -0,0 +1,41 @@ +#ifndef MARS_MExtractFixedWindowSpline +#define MARS_MExtractFixedWindowSpline + +#ifndef MARS_MExtractor +#include "MExtractor.h" +#endif + +class MExtractFixedWindowSpline : public MExtractor +{ + +private: + + static const Byte_t fgHiGainFirst; // Default for fHiGainFirst (now set to: 2) + static const Byte_t fgHiGainLast; // Default for fHiGainLast (now set to: 14) + static const Byte_t fgLoGainFirst; // Default for fLOGainFirst (now set to: 3) + static const Byte_t fgLoGainLast; // Default for fLoGainLast (now set to: 14) + + Float_t *fHiGainFirstDeriv; + Float_t *fLoGainFirstDeriv; + Float_t *fHiGainSecondDeriv; + Float_t *fLoGainSecondDeriv; + + Bool_t ReInit (MParList *pList); + + void FindSignalHiGain(Byte_t *ptr, Byte_t *logain, Int_t &sum, Byte_t &sat) const; + void FindSignalLoGain(Byte_t *ptr, Int_t &sum, Byte_t &sat) const; + +public: + + MExtractFixedWindowSpline(const char *name=NULL, const char *title=NULL); + ~MExtractFixedWindowSpline(); + + void SetRange(Byte_t hifirst=0, Byte_t hilast=0, Byte_t lofirst=0, Byte_t lolast=0); + + ClassDef(MExtractFixedWindowSpline, 0) // Task to Extract the Arrival Times using a Fast Spline +}; + +#endif + + + Index: /trunk/MagicSoft/Mars/msignal/Makefile =================================================================== --- /trunk/MagicSoft/Mars/msignal/Makefile (revision 4275) +++ /trunk/MagicSoft/Mars/msignal/Makefile (revision 4276) @@ -36,4 +36,5 @@ MExtractSlidingWindow.cc \ MExtractFixedWindowPeakSearch.cc \ + MExtractFixedWindowSpline.cc \ MExtractSignal.cc \ MExtractSignal2.cc \ Index: /trunk/MagicSoft/Mars/msignal/SignalLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/SignalLinkDef.h (revision 4275) +++ /trunk/MagicSoft/Mars/msignal/SignalLinkDef.h (revision 4276) @@ -18,4 +18,5 @@ #pragma link C++ class MExtractSlidingWindow+; #pragma link C++ class MExtractFixedWindowPeakSearch+; +#pragma link C++ class MExtractFixedWindowSpline+; #pragma link C++ class MExtractPINDiode++; #pragma link C++ class MExtractBlindPixel++; @@ -25,4 +26,5 @@ #pragma link C++ class MExtractTimeFastSpline+; #pragma link C++ class MExtractTimeHighestIntegral+; + #pragma link C++ class MArrivalTimeCam+; #pragma link C++ class MArrivalTimePix+; @@ -31,5 +33,3 @@ #pragma link C++ class MArrivalTimeCalc2+; - - #endif