Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 3632) +++ /trunk/MagicSoft/Mars/Changelog (revision 3633) @@ -19,4 +19,14 @@ -*-*- END OF LINE -*-*- + + 2004/04/02: Abelardo Moralejo + + * msignal/MExtractSignal3.[h,cc] + - added. Alternative methid for signal extraction. Uses sliding + window from event to event, but a common window for all pixels + in an event. + + * msignal/Makefile, SignalLinkDef.h + - added MExtractSignal3. 2004/04/02: Markus Gaug Index: /trunk/MagicSoft/Mars/msignal/MExtractSignal3.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractSignal3.cc (revision 3633) +++ /trunk/MagicSoft/Mars/msignal/MExtractSignal3.cc (revision 3633) @@ -0,0 +1,365 @@ +/* ======================================================================== *\ + ! + ! * + ! * 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): Abelardo Moralejo, 4/2004 + ! + ! Copyright: MAGIC Software Development, 2000-2004 + ! + ! + \* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MExtractSignal3 +// +// Calculate the signal integrating fWindowSize time slices, the same for +// all pixels. The integrated slices change from event to event, since the +// pulse positions in the FADC jump between events, but apparently in a +// "coherent" fashion. We first loop over all pixels and find the one +// which has the highest sum of fPeakSearchWindowSize (default: 4) consecutive +// non-saturated high gain slices. The position of the peak in this pixel +// determines the integration window position for all pixels of this event. +// For the moment we neglect time delays between pixels (which are in most +// cases small). The class is based on MExtractSignal2. +// +////////////////////////////////////////////////////////////////////////////// + +#include "MExtractSignal3.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MRawEvtData.h" +#include "MRawEvtPixelIter.h" + +#include "MPedestalCam.h" +#include "MPedestalPix.h" + +#include "MExtractedSignalCam.h" +#include "MExtractedSignalPix.h" + +ClassImp(MExtractSignal3); + +using namespace std; + +const Byte_t MExtractSignal3::fgSaturationLimit = 254; +const Byte_t MExtractSignal3::fgFirst = 0; +const Byte_t MExtractSignal3::fgLast = 14; +const Byte_t MExtractSignal3::fgWindowSize = 6; +const Byte_t MExtractSignal3::fgPeakSearchWindowSize = 4; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MExtractSignal3::MExtractSignal3(const char *name, const char *title) + : fNumHiGainSamples(15), fNumLoGainSamples(15), fSaturationLimit(fgSaturationLimit) +{ + + fName = name ? name : "MExtractSignal3"; + fTitle = title ? title : "Task to extract the signal from the FADC slices"; + + AddToBranchList("MRawEvtData.*"); + + SetRange(); +} + +void MExtractSignal3::SetRange(Byte_t windowh, Byte_t windowl, Byte_t peaksearchwindow) +{ + // + // Set windows to even number of slices due to clock noise (odd/even slice effect). + // + fWindowSizeHiGain = windowh & ~1; + fWindowSizeLoGain = windowl & ~1; + fPeakSearchWindowSize = peaksearchwindow & ~1; + + + if (fWindowSizeHiGain != windowh) + *fLog << endl << warn << + "MExtractSignal3::SetRange - Hi Gain window size has to be even, set to: " + << int(fWindowSizeHiGain) << " samples " << endl; + + if (fWindowSizeLoGain != windowl) + *fLog << endl << warn << + "MExtractSignal3::SetRange - Lo Gain window size has to be even, set to: " + << int(fWindowSizeLoGain) << " samples " << endl; + + if (fPeakSearchWindowSize != peaksearchwindow) + *fLog << endl << warn << + "MExtractSignal3::SetRange - Peak Search window size has to be even, set to: " + << int(fPeakSearchWindowSize) << " samples " << endl; + + + if (fWindowSizeHiGain<2) + { + fWindowSizeHiGain = 2; + *fLog << warn << "MExtractSignal3::SetRange - Hi Gain window size set to two samples" << endl; + } + + if (fWindowSizeLoGain<2) + { + fWindowSizeLoGain = 2; + *fLog << warn << "MExtractSignal3::SetRange - Lo Gain window size set to two samples" << endl; + } + + if (fPeakSearchWindowSize<2) + { + fPeakSearchWindowSize = 2; + *fLog << warn << "MExtractSignal3::SetRange - Peak Search window size set to two samples" << endl; + } + + + if (fWindowSizeHiGain > fNumHiGainSamples) + { + fWindowSizeHiGain = fNumHiGainSamples & ~1; + *fLog << warn << "MExtractSignal3::SetRange - Hi Gain window size set to " + << int(fWindowSizeHiGain) << " samples " << endl; + } + + if (fWindowSizeLoGain > fNumLoGainSamples) + { + fWindowSizeLoGain = fNumLoGainSamples & ~1; + *fLog << warn << "MExtractSignal3::SetRange - Lo Gain window size set to " + << int(fWindowSizeLoGain) << " samples " << endl; + } + + fWindowSqrtHiGain = TMath::Sqrt((Float_t)fWindowSizeHiGain); + fWindowSqrtLoGain = TMath::Sqrt((Float_t)fWindowSizeLoGain); + +} + + +// -------------------------------------------------------------------------- +// +// The PreProcess searches for the following input containers: +// - MRawEvtData +// - MPedestalCam +// +// The following output containers are also searched and created if +// they were not found: +// +// - MExtractedSignalCam +// +Int_t MExtractSignal3::PreProcess(MParList *pList) +{ + fRawEvt = (MRawEvtData*)pList->FindObject(AddSerialNumber("MRawEvtData")); + if (!fRawEvt) + { + *fLog << err << AddSerialNumber("MRawEvtData") << " not found... aborting." << endl; + return kFALSE; + } + + fSignals = (MExtractedSignalCam*)pList->FindCreateObj(AddSerialNumber("MExtractedSignalCam")); + if (!fSignals) + return kFALSE; + + // + // FIXME? We should keep track in MExtractedSignalCam of the signal extraction method used. + // + fSignals->SetUsedFADCSlices(0, fNumHiGainSamples-1, (Float_t)fWindowSizeHiGain, + 0, fNumLoGainSamples-1, (Float_t)fWindowSizeLoGain); + + fPedestals = (MPedestalCam*)pList->FindObject(AddSerialNumber("MPedestalCam")); + if (!fPedestals) + { + *fLog << err << AddSerialNumber("MPedestalCam") << " not found... aborting" << endl; + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// FindPeak +// Finds highest sum of "window" consecutive FADC slices in a pixel, and store +// in "startslice" the first slice of the group which yields the maximum sum. +// In "max" the value of the maximum sum is stored, in "sat" the number of +// saturated slices. +// +void MExtractSignal3::FindPeak(Byte_t *ptr, Byte_t window, Byte_t &startslice, + Int_t &max, Int_t &sat) const +{ + const Byte_t *end = ptr + fNumHiGainSamples; + + Int_t sum=0; + + // + // Calculate the sum of the first "window" slices + // + sat = 0; + Byte_t *p = ptr; + + while (p= fSaturationLimit) + sat++; + } + + // + // Check for saturation in all other slices + // + while (p= fSaturationLimit) + sat++; + + // + // Calculate the i-th sum of n as + // sum_i+1 = sum_i + slice[i+window] - slice[i] + // This is fast and accurate (because we are using int's) + // + max=sum; + for (p=ptr; p+window max) + { + max = sum; + startslice = p-ptr; + } + } + + return; +} + +// -------------------------------------------------------------------------- +// +// FindSignal +// Adds up "window" slices starting in slice to which "ptr" points. The result +// is stored in "sum", and the number of saturated slices in "sat". +// +void MExtractSignal3::FindSignal(Byte_t *ptr, Byte_t window, Int_t &sum, Int_t &sat) const +{ + // + // Calculate the sum of the "window" slices starting in ptr + // + sum=0; + sat = 0; + Byte_t *p = ptr; + + while (p= fSaturationLimit) + sat++; + } + + return; +} + +// -------------------------------------------------------------------------- +// +// Process +// First find the pixel with highest sum of fPeakSearchWindowSize slices (default:4) +// "startslice" will mark the slice at which the highest sum begins for that pixel. +// Then define the beginning of the integration window for ALL pixels as the slice +// before that: startslice-1 (this is somehow arbitrary), unless of course startslice=0, +// in which case we start at 0. We will also check that the integration window does not +// go beyond the FADC limits. +// +Int_t MExtractSignal3::Process() +{ + MRawEvtPixelIter pixel(fRawEvt); + + Int_t sat; + Int_t maxsumhi = 0; + Byte_t startslice; + Byte_t hiGainFirst = 0; + Byte_t loGainFirst = 0; + + while (pixel.Next()) + { + Int_t sumhi; + sat = 0; + + FindPeak(pixel.GetHiGainSamples(), fPeakSearchWindowSize, startslice, sumhi, sat); + if (sumhi > maxsumhi && sat == 0 ) + { + maxsumhi = sumhi; + if (startslice > 0) + hiGainFirst = startslice-1; + else + hiGainFirst = 0; + } + } + + + loGainFirst = hiGainFirst; + + // Make sure we will not integrate beyond the hi gain limit: + if (hiGainFirst+fWindowSizeHiGain > pixel.GetNumHiGainSamples()) + hiGainFirst = pixel.GetNumHiGainSamples()-fWindowSizeHiGain; + + // Make sure we will not integrate beyond the lo gain limit: + if (loGainFirst+fWindowSizeLoGain > pixel.GetNumLoGainSamples()) + loGainFirst = pixel.GetNumLoGainSamples()-fWindowSizeLoGain; + + pixel.Reset(); + fSignals->Clear(); + + sat = 0; + while (pixel.Next()) + { + // + // Find signal in hi- and lo-gain + // + Int_t sumhi, sathi; + + FindSignal(pixel.GetHiGainSamples()+hiGainFirst, fWindowSizeHiGain, sumhi, sathi); + + Int_t sumlo=0; + Int_t satlo=0; + if (pixel.HasLoGain()) + { + FindSignal(pixel.GetLoGainSamples()+loGainFirst, fWindowSizeLoGain, sumlo, satlo); + if (satlo) + sat++; + } + + // + // Take corresponding pedestal + // + const Int_t pixid = pixel.GetPixelId(); + + const MPedestalPix &ped = (*fPedestals)[pixid]; + MExtractedSignalPix &pix = (*fSignals)[pixid]; + + const Float_t pedes = ped.GetPedestal(); + const Float_t pedrms = ped.GetPedestalRms(); + + // + // Set extracted signal with pedestal substracted + // + pix.SetExtractedSignal(sumhi - pedes*fWindowSizeHiGain, pedrms*fWindowSqrtHiGain, + sumlo - pedes*fWindowSizeLoGain, pedrms*fWindowSqrtLoGain); + + pix.SetGainSaturation(sathi, sathi, satlo); + } /* while (pixel.Next()) */ + + // + // Print a warning if event has saturationg lo-gains + // + if (sat) + *fLog << warn << "WARNING - Lo Gain saturated in " << sat << " pixels." << endl; + + fSignals->SetReadyToSave(); + + return kTRUE; +} Index: /trunk/MagicSoft/Mars/msignal/MExtractSignal3.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractSignal3.h (revision 3633) +++ /trunk/MagicSoft/Mars/msignal/MExtractSignal3.h (revision 3633) @@ -0,0 +1,63 @@ +#ifndef MARS_MExtractSignal3 +#define MARS_MExtractSignal3 + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MRawEvtData; +class MRawRunHeader; + +class MPedestalCam; +class MExtractedSignalCam; + +class MExtractSignal3 : public MTask +{ +private: + static const Byte_t fgSaturationLimit; + static const Byte_t fgFirst; + static const Byte_t fgLast; + static const Byte_t fgWindowSize; + static const Byte_t fgPeakSearchWindowSize; + + MPedestalCam *fPedestals; // Pedestals of all pixels in the camera + MExtractedSignalCam *fSignals; // Extracted signal of all pixels in the camera + + MRawEvtData *fRawEvt; // raw event data (time slices) + MRawRunHeader *fRunHeader; // RunHeader information + + + Byte_t fPeakSearchWindowSize; // Size of FADC window in the search for the highest peak + // of all pixels. + + Byte_t fNumHiGainSamples; + Byte_t fNumLoGainSamples; + + Byte_t fWindowSizeHiGain; // Number of Hi Gain slices in window + Float_t fWindowSqrtHiGain; // Square root of number of Hi Gain slices in window + + Byte_t fWindowSizeLoGain; // Number of Lo Gain slices in window + Float_t fWindowSqrtLoGain; // Square root of number of Lo Gain slices in window + + Byte_t fSaturationLimit; + + void FindSignal(Byte_t *ptr, Byte_t window, Int_t &sum, Int_t &sat) const; + + void FindPeak(Byte_t *ptr, Byte_t window, Byte_t &startslice, Int_t &signal, Int_t &sat) const; + + + Int_t PreProcess(MParList *pList); + Int_t Process(); + +public: + MExtractSignal3(const char *name=NULL, const char *title=NULL); + + void SetRange(Byte_t windowh=fgWindowSize, Byte_t windowl=fgWindowSize, + Byte_t peaksearchwindow=fgPeakSearchWindowSize); + + void SetSaturationLimit(Byte_t lim) { fSaturationLimit = lim; } + + ClassDef(MExtractSignal3, 0) // Extracted Signal algorithm #3 +}; + +#endif Index: /trunk/MagicSoft/Mars/msignal/Makefile =================================================================== --- /trunk/MagicSoft/Mars/msignal/Makefile (revision 3632) +++ /trunk/MagicSoft/Mars/msignal/Makefile (revision 3633) @@ -39,4 +39,5 @@ MExtractSignal.cc \ MExtractSignal2.cc \ + MExtractSignal3.cc \ MExtractPINDiode.cc \ MExtractedSignalPINDiode.cc \ Index: /trunk/MagicSoft/Mars/msignal/SignalLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/SignalLinkDef.h (revision 3632) +++ /trunk/MagicSoft/Mars/msignal/SignalLinkDef.h (revision 3633) @@ -9,4 +9,5 @@ #pragma link C++ class MExtractSignal+; #pragma link C++ class MExtractSignal2+; +#pragma link C++ class MExtractSignal3+; #pragma link C++ class MExtractPINDiode++;