Index: /trunk/MagicSoft/GRB-Proposal/GRB_proposal_2005.tex =================================================================== --- /trunk/MagicSoft/GRB-Proposal/GRB_proposal_2005.tex (revision 6116) +++ /trunk/MagicSoft/GRB-Proposal/GRB_proposal_2005.tex (revision 6117) @@ -159,4 +159,5 @@ \bibitem{MESZAROS94} Meszaros P., Rees M., MNRAS 289 L41 (1994) \bibitem{NICOLAGRB} http://www.pd.infn.it/magic/GRB/grb.html +\bibitem{GCN} http://gcn.gsfc.nasa.gov/ \bibitem{GCNARCHIVE} http://lheawww.gsfc.nasa.gov/docs/gamcosray/legr/bacodine/gcn3\_archive.html \bibitem{GOTZ} D. Gotz, S. Mereghetti 2002 Observation of Gamma-ray Bursts with INTEGRAL Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 6116) +++ /trunk/MagicSoft/Mars/Changelog (revision 6117) @@ -29,4 +29,8 @@ - take out call to fSignals->Clear() which is already done in MTaskList. + + * msignal/MExtractTimeAndChargeDigitalFilterPeakSearch.[h,cc] + - new extractor reducing the bias and solving the problem with the + jumping pulse readout. Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.h (revision 6116) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.h (revision 6117) @@ -14,8 +14,8 @@ static const Byte_t fgLoGainSwitch; //! Default for fLoGainSwitch (now set to: 100) +protected: + Byte_t fLoGainFirstSave; //! Temporary variable to store the original position of low-gain start slice Float_t fLoGainStartShift; // Shift to start searching the low-gain signal obtained from the high-gain times. - -protected: Byte_t fLoGainSwitch; // Limit for max. bin content before the low-gain gets extracted Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilterPeakSearch.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilterPeakSearch.cc (revision 6117) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilterPeakSearch.cc (revision 6117) @@ -0,0 +1,395 @@ +/* ======================================================================== *\ +! +! * +! * 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): Markus Gaug, 05/2004 +! Author(s): Diego Tescaro, 05/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ +////////////////////////////////////////////////////////////////////////////// +// +// MExtractTimeAndChargeDigitalFilterPeakSearch +// +// Hendrik has promised to write more documentation +// +// +// The following variables have to be set by the derived class and +// do not have defaults: +// - fNumHiGainSamples +// - fNumLoGainSamples +// - fSqrtHiGainSamples +// - fSqrtLoGainSamples +// +// Input Containers: +// MRawEvtData +// MRawRunHeader +// MPedestalCam +// +// Output Containers: +// MArrivalTimeCam +// MExtractedSignalCam +// +////////////////////////////////////////////////////////////////////////////// +#include "MExtractTimeAndChargeDigitalFilterPeakSearch.h" + +#include +#include + +#include +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MPedestalPix.h" +#include "MPedestalCam.h" + +#include "MRawEvtPixelIter.h" + +#include "MExtractedSignalCam.h" +#include "MExtractedSignalPix.h" + +#include "MArrivalTimeCam.h" +#include "MArrivalTimePix.h" + +ClassImp(MExtractTimeAndChargeDigitalFilterPeakSearch); + +using namespace std; + +const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgHiGainFirst = 1; +const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgHiGainLast = 15; +const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgLoGainFirst = 3; +const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgLoGainLast = 14; +const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgOffsetLeftFromPeak = 1; +const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgOffsetRightFromPeak = 2; +const Byte_t MExtractTimeAndChargeDigitalFilterPeakSearch::fgPeakSearchWindowSize = 2; +// -------------------------------------------------------------------------- +// +// Default constructor. +// +// Sets: +// - fOffsetLeftFromPeak to fgOffsetLeftFromPeak +// - fOffsetRightFromPeak to fgOffsetRightFromPeak +// +MExtractTimeAndChargeDigitalFilterPeakSearch::MExtractTimeAndChargeDigitalFilterPeakSearch(const char *name, const char *title) +{ + fName = name ? name : "MExtractTimeAndChargeDigitalFilterPeakSearch"; + fTitle = title ? title : "Digital Filter"; + + SetOffsetLeftFromPeak(); + SetOffsetRightFromPeak(); + SetPeakSearchWindowSize(); +} + +// -------------------------------------------------------------------------- +// +// 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 MExtractTimeAndChargeDigitalFilterPeakSearch::FindPeak(Byte_t *ptr, Byte_t &startslice, Int_t &max, + Int_t &sat, Byte_t &satpos) const +{ + + const Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1; + + sat = 0; + satpos = 0; + + startslice = 0; + Int_t sum=0; + + // + // Calculate the sum of the first "fPeakSearchWindowSize" slices + // + sat = 0; + Byte_t *p = ptr; + + while (p= fSaturationLimit) + { + if (sat == 0) + satpos = p-ptr; + sat++; + } + } + + // + // Check for saturation in all other slices + // + while (p= fSaturationLimit) + { + if (sat == 0) + satpos = p-ptr; + sat++; + } + + // + // Calculate the i-th sum of n as + // sum_i+1 = sum_i + slice[i+fPeakSearchWindowSize] - slice[i] + // This is fast and accurate (because we are using int's) + // + max=sum; + for (p=ptr; p+fPeakSearchWindowSize max) + { + max = sum; + startslice = p-ptr+1; + } + } + + return; +} + + +// -------------------------------------------------------------------------- +// +// Process +// +// First find the pixel with highest sum of fPeakSearchWindowSize slices (default:2) +// "startslice" will mark the slice at which the highest sum begins for that pixel. +// +// Then define the beginning of the digital filter search window for ALL pixels as the slice +// before that: startslice-fOffsetLeftFromPeak until: startslice+fOffsetRightFromPeak +// +Int_t MExtractTimeAndChargeDigitalFilterPeakSearch::Process() +{ + + MRawEvtPixelIter pixel(fRawEvt); + + Int_t sat; + Byte_t satpos; + ULong_t gsatpos = 0; + + Int_t maxsumhi = -1000000; + Int_t numsat = 0; + Byte_t startslice; + + Byte_t hiGainFirstsave = fHiGainFirst; + Byte_t hiGainLastsave = fHiGainLast; + Byte_t loGainFirstsave = fLoGainFirst; + Byte_t loGainLastsave = fLoGainLast; + + Byte_t higainfirst = fHiGainFirst; + + while (pixel.Next()) + { + + Int_t sumhi; + sat = 0; + + FindPeak(pixel.GetHiGainSamples()+fHiGainFirst+fOffsetLeftFromPeak, startslice, sumhi, sat, satpos); + + if (sumhi > maxsumhi && sat == 0) + { + maxsumhi = sumhi; + higainfirst = fHiGainFirst + startslice; + } + else if (sat) + { + numsat++; + gsatpos += satpos; + } + } + + // + // Check necessary for calibration events + // + if (numsat > fSignals->GetSize()*0.9) + higainfirst = gsatpos/numsat - 1; + + // + // Shift the start slice to the left: + // + if (higainfirst > fOffsetLeftFromPeak) + fHiGainFirst = higainfirst - fOffsetLeftFromPeak; + else + *fLog << warn << " High Gain ranges out of limits to the left!!! " << endl; + + // + // Shift the last slice to the right: + // + if (higainfirst + fOffsetRightFromPeak + fWindowSizeHiGain > hiGainLastsave) + fHiGainLast = higainfirst + fOffsetRightFromPeak + fWindowSizeHiGain; + else + *fLog << warn << " High Gain ranges out of limits to the right!!! " << endl; + + if ( fHiGainFirst+(Int_t)fOffsetLoGain > fLoGainFirst ) + fLoGainFirst = fHiGainFirst + (Int_t)fOffsetLoGain; + else + *fLog << warn << " High Gain ranges out of limits to the left!!! " << endl; + + // + // Make sure we will not integrate beyond the lo gain limit: + // + if (fLoGainFirst+fWindowSizeLoGain+fOffsetRightFromPeak <= pixel.GetNumLoGainSamples()) + fLoGainLast = fLoGainFirst+fWindowSizeLoGain+fOffsetRightFromPeak; + else + *fLog << warn << " Low Gain ranges out of limits to the right!!! " << endl; + + pixel.Reset(); + + while (pixel.Next()) + { + // + // Find signal in hi- and lo-gain + // + Float_t sumhi =0., deltasumhi =0; // Set hi-gain of MExtractedSignalPix valid + Float_t timehi=0., deltatimehi=0; // Set hi-gain of MArrivalTimePix valid + Byte_t sathi=0; + + // Initialize fMaxBinContent for the case, it gets not set by the derived class + fMaxBinContent = fLoGainSwitch + 1; + + const Int_t pixidx = pixel.GetPixelId(); + const MPedestalPix &ped = (*fPedestals)[pixidx]; + const Bool_t higainabflag = pixel.HasABFlag(); + + FindTimeAndChargeHiGain(pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(), + sumhi, deltasumhi, + timehi, deltatimehi, + sathi, ped, higainabflag); + + // + // Make sure that in cases the time couldn't be correctly determined + // more meaningfull default values are assigned + // + if (timehi<0) + timehi = -1; + if (timehi>pixel.GetNumHiGainSamples()) + timehi = pixel.GetNumHiGainSamples(); + + Float_t sumlo =0., deltasumlo =-1.; // invalidate logain of MExtractedSignalPix + Float_t timelo=0., deltatimelo=-1; // invalidate logain of MArrivalTimePix + Byte_t satlo=0; + + // + // Adapt the low-gain extraction range from the obtained high-gain time + // + if (pixel.HasLoGain() && (fMaxBinContent > fLoGainSwitch) ) + { + deltasumlo = 0; // make logain of MExtractedSignalPix valid + deltatimelo = 0; // make logain of MArrivalTimePix valid + + fLoGainFirstSave = fLoGainFirst; + const Byte_t logainstart = sathi + ? sathi + (Int_t)fLoGainStartShift + : (Byte_t)(timehi + fLoGainStartShift); + + fLoGainFirst = logainstart > fLoGainFirstSave ? logainstart : fLoGainFirstSave; + + if ( fLoGainFirst < fLoGainLast ) + { + const Bool_t logainabflag = (higainabflag + pixel.GetNumHiGainSamples()) & 0x1; + FindTimeAndChargeLoGain(pixel.GetLoGainSamples()+fLoGainFirst, + sumlo, deltasumlo, + timelo, deltatimelo, + satlo, ped, logainabflag); + } + fLoGainFirst = fLoGainFirstSave; + + // Make sure that in cases the time couldn't be correctly determined + // more meaningfull default values are assigned + if (timelo<0) + timelo = -1; + if (timelo>pixel.GetNumLoGainSamples()) + timelo = pixel.GetNumLoGainSamples(); + } + + MExtractedSignalPix &pix = (*fSignals)[pixidx]; + MArrivalTimePix &tix = (*fArrTime)[pixidx]; + pix.SetExtractedSignal(sumhi, deltasumhi,sumlo, deltasumlo); + pix.SetGainSaturation(sathi, sathi, satlo); + + tix.SetArrivalTime(timehi, deltatimehi, timelo-fOffsetLoGain, deltatimelo); + tix.SetGainSaturation(sathi, sathi, satlo); + + } /* while (pixel.Next()) */ + + fArrTime->SetReadyToSave(); + fSignals->SetReadyToSave(); + + fHiGainFirst = hiGainFirstsave; + fHiGainLast = hiGainLastsave ; + fLoGainFirst = loGainFirstsave; + fLoGainLast = loGainLastsave ; + + return kTRUE; +} + + + +// -------------------------------------------------------------------------- +// +// Read the setup from a TEnv, eg: +// MJPedestal.MExtractor.WindowSizeHiGain: 6 +// MJPedestal.MExtractor.WindowSizeLoGain: 6 +// MJPedestal.MExtractor.BinningResolutionHiGain: 10 +// MJPedestal.MExtractor.BinningResolutionLoGain: 10 +// MJPedestal.MExtractor.WeightsFile: filename +// +Int_t MExtractTimeAndChargeDigitalFilterPeakSearch::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + + Bool_t rc = kFALSE; + + if (IsEnvDefined(env, prefix, "OffsetLeftFromPeak", print)) + { + fOffsetLeftFromPeak = GetEnvValue(env, prefix, fOffsetLeftFromPeak); + rc = kTRUE; + } + + if (IsEnvDefined(env, prefix, "OffsetRightFromPeak", print)) + { + fOffsetRightFromPeak = GetEnvValue(env, prefix, fOffsetRightFromPeak); + rc = kTRUE; + } + + if (IsEnvDefined(env, prefix, "PeakSearchWindowSize", print)) + { + fPeakSearchWindowSize = GetEnvValue(env, prefix, fPeakSearchWindowSize); + rc = kTRUE; + } + + return MExtractTimeAndChargeDigitalFilter::ReadEnv(env, prefix, print) ? kTRUE : rc; +} + + +void MExtractTimeAndChargeDigitalFilterPeakSearch::Print(Option_t *o) const +{ + if (IsA()==Class()) + *fLog << GetDescriptor() << ":" << endl; + + MExtractTimeAndChargeDigitalFilter::Print(o); + *fLog << " Offset from Peak left: " << fOffsetLeftFromPeak << endl; + *fLog << " Offset from Peak right: " << fOffsetRightFromPeak << endl; + *fLog << " Peak search window size: " << fPeakSearchWindowSize << endl; +} Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilterPeakSearch.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilterPeakSearch.h (revision 6117) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilterPeakSearch.h (revision 6117) @@ -0,0 +1,54 @@ +#ifndef MARS_MExtractTimeAndChargeDigitalFilterPeakSearch +#define MARS_MExtractTimeAndChargeDigitalFilterPeakSearch + +#ifndef MARS_MExtractTimeAndChargeDigitalFilter +#include "MExtractTimeAndChargeDigitalFilter.h" +#endif + +#ifndef MARS_MArrayF +#include "MArrayF.h" +#endif + +#ifndef MARS_MArrayI +#include "MArrayI.h" +#endif + +class TH1F; +class TH2F; +class MPedestalPix; +class MExtractTimeAndChargeDigitalFilterPeakSearch : public MExtractTimeAndChargeDigitalFilter +{ +private: + + static const Byte_t fgHiGainFirst; //! Default for fHiGainFirst (now set to: 0) + static const Byte_t fgHiGainLast; //! Default for fHiGainLast (now set to:18) + static const Byte_t fgLoGainFirst; //! Default for fLoGainFirst (now set to: 2) + static const Byte_t fgLoGainLast; //! Default for fLoGainLast (now set to:14) + static const Byte_t fgOffsetLeftFromPeak; //! Default for fOffsetLeftFromPeak (now set to: 1) + static const Byte_t fgOffsetRightFromPeak; //! Default for fOffsetRightFromPeak (now set to: 2) + static const Byte_t fgPeakSearchWindowSize; //! Default for fPeakSearchWindowSize (now set to: 2) + + Byte_t fOffsetLeftFromPeak; // Number of slices to start extraction before peak slice + Byte_t fOffsetRightFromPeak; // Number of slices to stop extraction after peak slice + Byte_t fPeakSearchWindowSize; // Size of FADC window in the search for the highest peak of all pixels. + + void FindPeak(Byte_t *ptr, Byte_t &startslice, Int_t &signal, Int_t &sat, Byte_t &satpos) const; + + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); + Int_t Process(); + +public: + + MExtractTimeAndChargeDigitalFilterPeakSearch(const char *name=NULL, const char *title=NULL); + ~MExtractTimeAndChargeDigitalFilterPeakSearch() { } + + void SetOffsetLeftFromPeak ( Byte_t offset=fgOffsetLeftFromPeak ) { fOffsetLeftFromPeak = offset; } + void SetOffsetRightFromPeak( Byte_t offset=fgOffsetRightFromPeak ) { fOffsetRightFromPeak = offset; } + void SetPeakSearchWindowSize(Byte_t size =fgPeakSearchWindowSize ) { fPeakSearchWindowSize = size; } + + void Print(Option_t *o="") const; + + ClassDef(MExtractTimeAndChargeDigitalFilterPeakSearch, 1) // Digital filter with global Peak Search +}; + +#endif Index: /trunk/MagicSoft/Mars/msignal/Makefile =================================================================== --- /trunk/MagicSoft/Mars/msignal/Makefile (revision 6116) +++ /trunk/MagicSoft/Mars/msignal/Makefile (revision 6117) @@ -52,4 +52,5 @@ MExtractTimeAndChargeSpline.cc \ MExtractTimeAndChargeDigitalFilter.cc \ + MExtractTimeAndChargeDigitalFilterPeakSearch.cc \ MArrivalTime.cc \ MArrivalTimeCalc.cc \ Index: /trunk/MagicSoft/Mars/msignal/SignalLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/SignalLinkDef.h (revision 6116) +++ /trunk/MagicSoft/Mars/msignal/SignalLinkDef.h (revision 6117) @@ -33,4 +33,5 @@ #pragma link C++ class MExtractTimeAndChargeSpline+; #pragma link C++ class MExtractTimeAndChargeDigitalFilter+; +#pragma link C++ class MExtractTimeAndChargeDigitalFilterPeakSearch+;