Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 5231) +++ /trunk/MagicSoft/Mars/Changelog (revision 5232) @@ -28,4 +28,11 @@ * msignal/MExtractTimeAndChargeSpline.[h,cc] - fixed class documentation and some last bugs. + + * msignal/MExtractTimeAndChargeSlidingWindow.[h,cc] + * msignal/Makefile + * msignal/SignalLinkDef.h + - the combined extractor class sliding window and highest integral + with pedestal-AB-flag corrected, to be used for the TDAS-extractor + 2004/10/08: Markus Meyer and Keiichi Mase Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSlidingWindow.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSlidingWindow.cc (revision 5232) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSlidingWindow.cc (revision 5232) @@ -0,0 +1,482 @@ +/* ======================================================================== *\ +! +! * +! * 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): Thomas Bretz, 02/2004 +! Author(s): Hendrik Bartko, 01/2004 +! Author(s): Markus Gaug, 09/2004 +! +! Copyright: MAGIC Software Development, 2002-2004 +! +! +\* ======================================================================== */ +////////////////////////////////////////////////////////////////////////////// +// +// MExtractSlidingWindow +// +// Extracts the signal from a sliding window of size fHiGainWindowSize and +// fLoGainWindowSize. The signal is the one which maximizes the integral +// contents. +// +// Call: SetRange(higainfirst, higainlast, logainfirst, logainlast) +// to modify the ranges in which the window is allowed to move. +// Defaults are: +// +// fHiGainFirst = fgHiGainFirst = 0 +// fHiGainLast = fgHiGainLast = 14 +// fLoGainFirst = fgLoGainFirst = 3 +// fLoGainLast = fgLoGainLast = 14 +// +// Call: SetWindowSize(windowhigain, windowlogain) +// to modify the sliding window widths. Windows have to be an even number. +// In case of odd numbers, the window will be modified. +// +// Defaults are: +// +// fHiGainWindowSize = fgHiGainWindowSize = 6 +// fLoGainWindowSize = fgLoGainWindowSize = 6 +// +////////////////////////////////////////////////////////////////////////////// +#include "MExtractTimeAndChargeSlidingWindow.h" + +#include "MPedestalPix.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MExtractTimeAndChargeSlidingWindow); + +using namespace std; + +const Byte_t MExtractTimeAndChargeSlidingWindow::fgHiGainFirst = 2; +const Byte_t MExtractTimeAndChargeSlidingWindow::fgHiGainLast = 14; +const Byte_t MExtractTimeAndChargeSlidingWindow::fgLoGainFirst = 3; +const Byte_t MExtractTimeAndChargeSlidingWindow::fgLoGainLast = 14; +const Byte_t MExtractTimeAndChargeSlidingWindow::fgHiGainWindowSize = 6; +const Byte_t MExtractTimeAndChargeSlidingWindow::fgLoGainWindowSize = 6; +// -------------------------------------------------------------------------- +// +// Default constructor. +// +// Calls: +// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast) +// +// Initializes: +// - fWindowSizeHiGain to fgHiGainWindowSize +// - fWindowSizeLoGain to fgLoGainWindowSize +// +MExtractTimeAndChargeSlidingWindow::MExtractTimeAndChargeSlidingWindow(const char *name, const char *title) + : fWindowSizeHiGain(fgHiGainWindowSize), + fWindowSizeLoGain(fgLoGainWindowSize), + fHiGainSignal(NULL), fLoGainSignal(NULL) +{ + + fName = name ? name : "MExtractTimeAndChargeSlidingWindow"; + fTitle = title ? title : "Calculate arrival times and charges using a sliding window"; + + SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast); +} + + +// -------------------------------------------------------------------------- +// +// Destructor: Deletes the arrays +// +MExtractTimeAndChargeSlidingWindow::~MExtractTimeAndChargeSlidingWindow() +{ + + if (fHiGainSignal) + delete [] fHiGainSignal; + if (fLoGainSignal) + delete [] fLoGainSignal; + +} + +//------------------------------------------------------------------- +// +// Set the ranges +// +// Calls: +// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast); +// - SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain); +// +void MExtractTimeAndChargeSlidingWindow::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast) +{ + + MExtractor::SetRange(hifirst, hilast, lofirst, lolast); + + // + // Redo the checks if the window is still inside the ranges + // + SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain); + +} + +// ----------------------------------------------------------------------------------------- +// +// Checks: +// - if a window is bigger than the one defined by the ranges, set it to the available range +// - if a window is smaller than 2, set it to 2 +// +// Sets: +// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain +// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain +// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples) +// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples) +// +void MExtractTimeAndChargeSlidingWindow::SetWindowSize(Byte_t windowh, Byte_t windowl) +{ + + const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1; + + if (fWindowSizeHiGain > availhirange) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",(int)fWindowSizeHiGain, + " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl; + *fLog << warn << GetDescriptor() + << ": Will set window size to: " << (int)availhirange << endl; + fWindowSizeHiGain = availhirange; + } + + if (fWindowSizeHiGain<2) + { + fWindowSizeHiGain = 2; + *fLog << warn << GetDescriptor() << ": High Gain window size too small, set to two samples" << endl; + } + + if (fLoGainLast != 0 && fWindowSizeLoGain != 0) + { + const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1; + + if (fWindowSizeLoGain > availlorange) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",(int)fWindowSizeLoGain, + " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl; + *fLog << warn << GetDescriptor() + << ": Will set window size to: " << (int)availlorange << endl; + fWindowSizeLoGain = availlorange; + } + + if (fWindowSizeLoGain<2) + { + fWindowSizeLoGain = 2; + *fLog << warn << GetDescriptor() << ": Low Gain window size too small set to two samples" << endl; + } + } + + fNumHiGainSamples = (Float_t)fWindowSizeHiGain; + fNumLoGainSamples = fLoGainLast ? (Float_t)fWindowSizeLoGain : 0.; + + fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples); + fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples); + +} + +// -------------------------------------------------------------------------- +// +// ReInit +// +// Calls: +// - MExtractTimeAndCharge::ReInit(pList); +// - Deletes all arrays, if not NULL +// - Creates new arrays according to the extraction range +// +Bool_t MExtractTimeAndChargeSlidingWindow::ReInit(MParList *pList) +{ + + if (!MExtractTimeAndCharge::ReInit(pList)) + return kFALSE; + + if (fHiGainSignal) + delete [] fHiGainSignal; + if (fLoGainSignal) + delete [] fLoGainSignal; + + Int_t range = fHiGainLast - fHiGainFirst + 1 + fHiLoLast; + + fHiGainSignal = new Float_t[range]; + memset(fHiGainSignal,0,range*sizeof(Float_t)); + + range = fLoGainLast - fLoGainFirst + 1; + + fLoGainSignal = new Float_t[range]; + memset(fLoGainSignal,0,range*sizeof(Float_t)); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calculates the arrival time for each pixel +// +void MExtractTimeAndChargeSlidingWindow::FindTimeAndChargeHiGain(Byte_t *first, Byte_t *logain, Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag) +{ + + Int_t range = fHiGainLast - fHiGainFirst + 1; + const Byte_t *end = first + range; + Byte_t *p = first; + + Float_t max = 0; // highest integral content of all windows + Int_t count = 0; // counter to recognize the AB-flag + + Float_t pedes = ped.GetPedestal(); + const Float_t ABoffs = ped.GetPedestalABoffset(); + + Float_t PedMean[2]; + PedMean[0] = pedes + ABoffs; + PedMean[1] = pedes - ABoffs; + + // + // Check for saturation in all other slices + // + while (p= fSaturationLimit) + sat++; + + count++; + } + + // + // Check for saturation in all other slices + // + while (p= fSaturationLimit) + sat++; + + // + // Calculate the i-th sum as + // sum_i+1 = sum_i + slice[i+8] - slice[i] + // This is fast and accurate (because we are using int's) + // + count = 0; + max = sum; + Int_t idx = 0; // idx of the first slice of the maximum window + + for (p=first; p+fWindowSizeHiGainmax) + { + max = sum; + idx = count+1; + } + count++; + } + + if (fHiLoLast != 0) + { + + // + // overlap bins + // + Byte_t *l = logain; + + while (p < end && l < logain+fHiLoLast) + { + + const Int_t ids = fHiGainFirst + count + fWindowSizeHiGain; + const Float_t signal = (Float_t)*l - PedMean[(ids+abflag) & 0x1]; + sum += signal - fHiGainSignal[count]; + fHiGainSignal[count + fWindowSizeHiGain] = signal; + count++; + if (*l++ >= fSaturationLimit) + sat++; + + if (sum>max) + { + max = sum; + idx = count+1; + } + p++; + } + + if (fHiLoLast > fWindowSizeHiGain) + { + while (l < logain + fHiLoLast) + { + const Int_t ids = fHiGainFirst + count + fWindowSizeHiGain; + const Float_t signal = (Float_t)*l - PedMean[(ids+abflag) & 0x1]; + sum += signal - fHiGainSignal[count]; + fHiGainSignal[count+fWindowSizeHiGain] = signal; + count++; + if (*(l+fWindowSizeHiGain) >= fSaturationLimit) + sat++; + + if (sum>max) + { + max = sum; + idx = count+1; + } + l++; + } + } + } + + // + // now calculate the time for the maximum window + // + Float_t timesignalsum = 0.; + Int_t timesquaredsum = 0; + + for (Int_t i=idx; i= fSaturationLimit) + sat++; + + count++; + } + + // + // Check for saturation in all other slices + // + while (p= fSaturationLimit) + sat++; + + // + // Calculate the i-th sum as + // sum_i+1 = sum_i + slice[i+8] - slice[i] + // This is fast and accurate (because we are using int's) + // + count = 0; + max = sum; + Int_t idx = 0; // idx of the first slice of the maximum window + + for (p=first; p+fWindowSizeLoGainmax) + { + max = sum; + idx = count+1; + } + count++; + } + + // + // now calculate the time for the maximum window + // + Float_t timesignalsum = 0; + Int_t timesquaredsum = 0; + + for (Int_t i=idx; i