Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 3970) +++ /trunk/MagicSoft/Mars/Changelog (revision 3971) @@ -18,4 +18,14 @@ -*-*- END OF LINE -*-*- + + 2004/05/05: Markus Gaug + + * msignal/Makefile + * msignal/SignalLinkDef.h + * msignal/MExtractTimeFastSpline.[h,cc] + - new fast spline extractor for the equally spaced time slices. + Searches for the position of the half maximum between maximum and + pedestal. About 6 times faster than MArrivalTimeCalc + 2004/05/04: Raquel de los Reyes Index: /trunk/MagicSoft/Mars/NEWS =================================================================== --- /trunk/MagicSoft/Mars/NEWS (revision 3970) +++ /trunk/MagicSoft/Mars/NEWS (revision 3971) @@ -13,4 +13,7 @@ *** Version 0.8.4 (2004/04/19) + + - new fast arrival time extractor using cubic splines: + MExtractTimeFastSpline - implementes multi-argument support in MDataChain Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc (revision 3971) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc (revision 3971) @@ -0,0 +1,654 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MExtractTimeFastSpline +// +// Fast arrival Time extractor using a cubic spline algorithm of Numerical Recipes. +// It returns the position of the half maximum between absolute maximum +// and pedestal of the spline that interpolates the FADC slices. +// +// The precision of the half-maximum searches can be chosen by: +// SetPrecision(). +// +// The precision of the maximum-finder is fixed to 0.025 FADC units. +// +////////////////////////////////////////////////////////////////////////////// +#include "MExtractTimeFastSpline.h" + +#include "MPedestalPix.h" + +#include "MLog.h" +#include "MLogManip.h" + + +ClassImp(MExtractTimeFastSpline); + +using namespace std; + +const Byte_t MExtractTimeFastSpline::fgHiGainFirst = 2; +const Byte_t MExtractTimeFastSpline::fgHiGainLast = 14; +const Byte_t MExtractTimeFastSpline::fgLoGainFirst = 3; +const Byte_t MExtractTimeFastSpline::fgLoGainLast = 14; +const Float_t MExtractTimeFastSpline::fgResolution = 0.003; +// -------------------------------------------------------------------------- +// +// Default constructor. +// +// Calls: +// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast) +// +// Initializes: +// - fResolution to fgResolution +// +MExtractTimeFastSpline::MExtractTimeFastSpline(const char *name, const char *title) + : fHiGainFirstDeriv(NULL), fLoGainFirstDeriv(NULL), + fHiGainSecondDeriv(NULL), fLoGainSecondDeriv(NULL) +{ + + fName = name ? name : "MExtractTimeFastSpline"; + fTitle = title ? title : "Calculate photons arrival time using a fast spline"; + + SetResolution(); + SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast); + +} + +MExtractTimeFastSpline::~MExtractTimeFastSpline() +{ + + if (fHiGainFirstDeriv) + delete fHiGainFirstDeriv; + if (fLoGainFirstDeriv) + delete fLoGainFirstDeriv; + if (fHiGainSecondDeriv) + delete fHiGainSecondDeriv; + if (fLoGainSecondDeriv) + delete fLoGainSecondDeriv; + +} + + +// -------------------------------------------------------------------------- +// +// SetRange: +// +// Calls: +// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast); +// - Deletes x, if not NULL +// - Creates x according to the range +// +void MExtractTimeFastSpline::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast) +{ + + MExtractor::SetRange(hifirst,hilast,lofirst,lolast); + + if (fHiGainFirstDeriv) + delete fHiGainFirstDeriv; + if (fLoGainFirstDeriv) + delete fLoGainFirstDeriv; + if (fHiGainSecondDeriv) + delete fHiGainSecondDeriv; + if (fLoGainSecondDeriv) + delete fLoGainSecondDeriv; + + Int_t range = fHiGainLast - fHiGainFirst + 1; + + if (range < 2) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s",": Hi-Gain Extraction range [",(int)fHiGainFirst,"," + ,fHiGainLast,"] too small, ") << endl; + *fLog << warn << GetDescriptor() + << " will move higher limit to obtain 4 slices " << endl; + SetRange(fHiGainFirst, fHiGainLast+4-range,fLoGainFirst,fLoGainLast); + range = fHiGainLast - fHiGainFirst + 1; + } + + + 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; + + if (range >= 2) + { + + fLoGainFirstDeriv = new Float_t[range]; + memset(fLoGainFirstDeriv,0,range*sizeof(Float_t)); + fLoGainSecondDeriv = new Float_t[range]; + memset(fLoGainSecondDeriv,0,range*sizeof(Float_t)); + + } +} + + +// -------------------------------------------------------------------------- +// +// Calculates the arrival time for each pixel +// +void MExtractTimeFastSpline::FindTimeHiGain(Byte_t *first, Float_t &time, Float_t &dtime, + Byte_t &sat, const MPedestalPix &ped) const +{ + + const Int_t range = fHiGainLast - fHiGainFirst + 1; + const Byte_t *end = first + range; + Byte_t *p = first; + Byte_t max = 0; + Byte_t maxpos = 0; + + // + // Check for saturation in all other slices + // + while (++p max) + { + max = *p; + maxpos = p-first; + } + + if (*p >= fSaturationLimit) + { + sat++; + break; + } + } + + if (sat) + return; + + if (maxpos < 2) + return; + + Float_t pp; + + p = first; + fHiGainSecondDeriv[0] = 0.; + fHiGainFirstDeriv[0] = 0.; + + for (Int_t i=1;i=0;k--) + fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k]; + + // + // Now find the maximum + // + Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one + Float_t lower = (Float_t)maxpos-1.; + Float_t upper = (Float_t)maxpos; + Float_t abmax = 0.; + Float_t x = lower; + Float_t y = 0.; + Float_t a = 1.; + Float_t b = 0.; + Int_t klo = maxpos-1; + Int_t khi = maxpos; + Float_t klocont = (Float_t)*(first+klo); + Float_t khicont = (Float_t)*(first+khi); + time = lower; + + // + // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to + // interval maxpos+1. + // + while (x abmax) + { + abmax = y; + time = x; + } + + } + + if (time > upper-0.1) + { + + upper = (Float_t)maxpos+1.; + lower = (Float_t)maxpos; + x = lower; + a = 1.; + b = 0.; + khi = maxpos+1; + klo = maxpos; + klocont = (Float_t)*(first+klo); + khicont = (Float_t)*(first+khi); + + while (x abmax) + { + abmax = y; + time = x; + } + } + } + + const Float_t up = time+step+0.015; + + x = time; + a = upper - x; + b = x - lower; + + step = 0.025; // step size of 83 ps + + while (x abmax) + { + abmax = y; + time = x; + } + + } + + Float_t lo = time-0.225; + + x = time; + a = upper - x; + b = x - lower; + khi--; + klo--; + klocont = (Float_t)*(first+klo); + khicont = (Float_t)*(first+khi); + + while (x>lo) + { + + x -= step; + a += step; + b -= step; + + y = a* klocont + + b* khicont + + ( (a*a*a-a)*fHiGainSecondDeriv[klo] + +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0; + + if (y > abmax) + { + abmax = y; + time = x; + } + + } + + const Float_t pedes = ped.GetPedestal(); + const Float_t halfmax = pedes + (abmax - pedes)/2.; + + // + // 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 > maxpos-4) + { + if (*(first+klo) < (Byte_t)halfmax) + break; + klo--; + } + + // + // Loop from the beginning of the slice upwards to reach the halfmax: + // With means of bisection: + // + x = (Float_t)klo; + a = 1.; + b = 0.; + klocont = (Float_t)*(first+klo); + khicont = (Float_t)*(first+klo+1); + time = x; + + step = 0.5; + Bool_t back = kFALSE; + + while (step > fResolution) + { + + if (back) + { + x -= step; + a += step; + b -= step; + } + else + { + x += step; + a -= step; + b += step; + } + + y = a*klocont + + b*khicont + + ( (a*a*a-a)*fHiGainSecondDeriv[klo] + +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0; + + if (y >= halfmax) + back = kTRUE; + else + back = kFALSE; + + step /= 2.; + } + + time = x; + dtime = fResolution; +} + + +// -------------------------------------------------------------------------- +// +// Calculates the arrival time for each pixel +// +void MExtractTimeFastSpline::FindTimeLoGain(Byte_t *first, Float_t &time, Float_t &dtime, + Byte_t &sat, const MPedestalPix &ped) const +{ + + const Int_t range = fLoGainLast - fLoGainFirst + 1; + const Byte_t *end = first + range; + Byte_t *p = first; + Byte_t max = 0; + Byte_t maxpos = 0; + + // + // Check for saturation in all other slices + // + while (++p max) + { + max = *p; + maxpos = p-first; + } + + if (*p >= fSaturationLimit) + { + sat++; + break; + } + } + + if (sat) + return; + + if (maxpos < 2) + return; + + Float_t pp; + + p = first; + fLoGainSecondDeriv[0] = 0.; + fLoGainFirstDeriv[0] = 0.; + + for (Int_t i=1;i=0;k--) + fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k]; + + // + // Now find the maximum + // + Float_t step = 0.2; // start with step size of 1ns and loop again with the smaller one + Float_t lower = (Float_t)maxpos-1.; + Float_t upper = (Float_t)maxpos; + Float_t abmax = 0.; + Float_t x = lower; + Float_t y = 0.; + Float_t a = 1.; + Float_t b = 0.; + Int_t klo = maxpos-1; + Int_t khi = maxpos; + Float_t klocont = (Float_t)*(first+klo); + Float_t khicont = (Float_t)*(first+khi); + time = lower; + + // + // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to + // interval maxpos+1. + // + while (x abmax) + { + abmax = y; + time = x; + } + + } + + if (time > upper-0.1) + { + + upper = (Float_t)maxpos+1.; + lower = (Float_t)maxpos; + x = lower; + a = 1.; + b = 0.; + khi = maxpos+1; + klo = maxpos; + klocont = (Float_t)*(first+klo); + khicont = (Float_t)*(first+khi); + + while (x abmax) + { + abmax = y; + time = x; + } + } + } + + const Float_t up = time+step+0.015; + + x = time; + a = upper - x; + b = x - lower; + + step = 0.025; // step size of 165 ps + + while (x abmax) + { + abmax = y; + time = x; + } + + } + + Float_t lo = time-0.225; + + x = time; + a = upper - x; + b = x - lower; + khi--; + klo--; + klocont = (Float_t)*(first+klo); + khicont = (Float_t)*(first+khi); + + while (x>lo) + { + + x -= step; + a += step; + b -= step; + + y = a* klocont + + b* khicont + + ( (a*a*a-a)*fLoGainSecondDeriv[klo] + +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0; + + if (y > abmax) + { + abmax = y; + time = x; + } + + } + + const Float_t pedes = ped.GetPedestal(); + const Float_t halfmax = pedes + (abmax - pedes)/2.; + + // + // 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 > maxpos-4) + { + if (*(first+klo) < (Byte_t)halfmax) + break; + klo--; + } + + // + // Loop from the beginning of the slice upwards to reach the halfmax: + // With means of bisection: + // + x = (Float_t)klo; + a = 1.; + b = 0.; + klocont = (Float_t)*(first+klo); + khicont = (Float_t)*(first+klo+1); + time = x; + + step = 0.5; + Bool_t back = kFALSE; + + while (step > fResolution) + { + + if (back) + { + x -= step; + a += step; + b -= step; + } + else + { + x += step; + a -= step; + b += step; + } + + y = a*klocont + + b*khicont + + ( (a*a*a-a)*fLoGainSecondDeriv[klo] + +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0; + + if (y >= halfmax) + back = kTRUE; + else + back = kFALSE; + + step /= 2.; + } + + time = x; + dtime = fResolution; +} + + Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.h (revision 3971) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.h (revision 3971) @@ -0,0 +1,42 @@ +#ifndef MARS_MExtractTimeFastSpline +#define MARS_MExtractTimeFastSpline + +#ifndef MARS_MExtractTime +#include "MExtractTime.h" +#endif + +class MPedestalPix; +class MExtractTimeFastSpline : public MExtractTime +{ + + 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) + static const Float_t fgResolution; // Default for fResolution (now set to: 0.003) + + Float_t *fHiGainFirstDeriv; + Float_t *fLoGainFirstDeriv; + Float_t *fHiGainSecondDeriv; + Float_t *fLoGainSecondDeriv; + + Float_t fResolution; // The time resolution in FADC units + + void FindTimeHiGain(Byte_t *first, Float_t &time, Float_t &dtime, Byte_t &sat, const MPedestalPix &ped) const; + void FindTimeLoGain(Byte_t *first, Float_t &time, Float_t &dtime, Byte_t &sat, const MPedestalPix &ped) const; + +public: + + MExtractTimeFastSpline(const char *name=NULL, const char *title=NULL); + ~MExtractTimeFastSpline(); + + void SetRange(Byte_t hifirst=0, Byte_t hilast=0, Byte_t lofirst=0, Byte_t lolast=0); + void SetResolution(Float_t f=fgResolution) { fResolution = f; } + + ClassDef(MExtractTimeFastSpline, 0) // Task to Extract the Arrival Times using a Fast Spline +}; + +#endif + + + Index: /trunk/MagicSoft/Mars/msignal/Makefile =================================================================== --- /trunk/MagicSoft/Mars/msignal/Makefile (revision 3970) +++ /trunk/MagicSoft/Mars/msignal/Makefile (revision 3971) @@ -43,4 +43,5 @@ MExtractTime.cc \ MExtractTimeSpline.cc \ + MExtractTimeFastSpline.cc \ MExtractTimeHighestIntegral.cc \ MArrivalTime.cc \ Index: /trunk/MagicSoft/Mars/msignal/SignalLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/SignalLinkDef.h (revision 3970) +++ /trunk/MagicSoft/Mars/msignal/SignalLinkDef.h (revision 3971) @@ -23,4 +23,5 @@ #pragma link C++ class MExtractTime+; #pragma link C++ class MExtractTimeSpline+; +#pragma link C++ class MExtractTimeFastSpline+; #pragma link C++ class MExtractTimeHighestIntegral+; #pragma link C++ class MArrivalTimeCam;