Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 5145) +++ /trunk/MagicSoft/Mars/Changelog (revision 5146) @@ -19,4 +19,18 @@ -*-*- END OF LINE -*-*- + + 2004/09/29: Hendrik Bartko, Diego Tescaro, Markus Gaug + + * msignal/MExtractTimeAndChargeDigitalFilter.[h,cc] + - Hendriks digital filter finally implemented + + * msignal/Makefile + * msignal/SignalLinkDef.h + - added MExtractTimeAndChargeDigitalFilter and MExtractTimeAndCharge + + * msignal/MExtractTime.h + - made PreProcess, ReInit and Process protected + + 2004/09/29: Thomas Bretz Index: /trunk/MagicSoft/Mars/msignal/MExtractTime.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTime.h (revision 5145) +++ /trunk/MagicSoft/Mars/msignal/MExtractTime.h (revision 5146) @@ -12,11 +12,9 @@ { private: + static const Float_t fgOffsetLoGain; // Default for fOffsetLoGain (now set to 1.51 (= 5ns) - Int_t PreProcess(MParList *pList); - Bool_t ReInit(MParList *pList); - Int_t Process(); - protected: + Float_t fOffsetLoGain; // Offset of the low-gain signal w.r.t. the High-Gain slices @@ -27,4 +25,8 @@ virtual void FindTimeLoGain(Byte_t *firstused, Float_t &time, Float_t &dtime, Byte_t &sat, const MPedestalPix &ped) const; + + Int_t PreProcess( MParList *pList ); + Bool_t ReInit ( MParList *pList ); + Int_t Process (); Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.cc (revision 5145) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.cc (revision 5146) @@ -126,5 +126,5 @@ void MExtractTimeAndCharge::FindTimeAndChargeHiGain(Byte_t *firstused, Byte_t *logain, Float_t &sum, Float_t &dsum, Float_t &time, Float_t &dtime, - Byte_t &sat, const MPedestalPix &ped) + Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag) { return; @@ -133,5 +133,5 @@ void MExtractTimeAndCharge::FindTimeAndChargeLoGain(Byte_t *firstused, Float_t &sum, Float_t &dsum, Float_t &time, Float_t &dtime, - Byte_t &sat, const MPedestalPix &ped) + Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag) { return; @@ -170,5 +170,5 @@ sumhi, deltasumhi, timehi, deltatimehi, - sathi, ped); + sathi, ped, pixel.HasABFlag()); Float_t sumlo =0., deltasumlo =0.; @@ -179,5 +179,5 @@ sumlo, deltasumlo, timelo, deltatimelo, - satlo, ped); + satlo, ped, pixel.HasABFlag()); MExtractedSignalPix &pix = (*fSignals)[pixid]; @@ -186,7 +186,4 @@ pix.SetExtractedSignal(sumhi, deltasumhi,sumlo, deltasumlo); pix.SetGainSaturation(sathi, sathi, satlo); - - pix.SetNumHiGainSlices(fNumHiGainSamples); - pix.SetNumLoGainSlices(fNumLoGainSamples); tix.SetArrivalTime(timehi, deltatimehi, timelo-fOffsetLoGain, deltatimelo); Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.h (revision 5145) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.h (revision 5146) @@ -10,4 +10,5 @@ { private: + Int_t PreProcess(MParList *pList); Bool_t ReInit(MParList *pList); @@ -15,10 +16,11 @@ protected: + virtual void FindTimeAndChargeHiGain(Byte_t *firstused, Byte_t *logain, Float_t &sum, Float_t &dsum, Float_t &time, Float_t &dtime, - Byte_t &sat, const MPedestalPix &ped); + Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag); virtual void FindTimeAndChargeLoGain(Byte_t *firstused, Float_t &sum, Float_t &dsum, Float_t &time, Float_t &dtime, - Byte_t &sat, const MPedestalPix &ped); + Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag); public: Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc (revision 5146) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.cc (revision 5146) @@ -0,0 +1,471 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MExtractTimeAndChargeDigitalFilter +// +// 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 "MExtractTimeAndChargeDigitalFilter.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MPedestalPix.h" + +#include "TFile.h" // remove!! +#include "TH1F.h" // remove!! +#include "TString.h" + +ClassImp(MExtractTimeAndChargeDigitalFilter); + +using namespace std; + +const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainFirst = 0; +const Byte_t MExtractTimeAndChargeDigitalFilter::fgHiGainLast = 14; +const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainFirst = 3; +const Byte_t MExtractTimeAndChargeDigitalFilter::fgLoGainLast = 14; +const Int_t MExtractTimeAndChargeDigitalFilter::fgHiGainWindowSize = 6; +const Int_t MExtractTimeAndChargeDigitalFilter::fgLoGainWindowSize = 6; +const Int_t MExtractTimeAndChargeDigitalFilter::fgBinningResolution = 10; +// -------------------------------------------------------------------------- +// +// Default constructor. +// +// Calls: +// - SetWindowSize(); +// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast) +// - SetBinningResolution(); +// +MExtractTimeAndChargeDigitalFilter::MExtractTimeAndChargeDigitalFilter(const char *name, const char *title) +{ + + fName = name ? name : "MExtractTimeAndChargeDigitalFilter"; + fTitle = title ? title : "Digital Filter"; + + SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast); + SetWindowSize(); + SetBinningResolution(); + + for (int i=0; i<60; i++){ + fw_amp[i] = 1; + fw_time[i] = 1; + } +} + + +// --------------------------------------------------------------------------------------- +// +// Checks: +// - if a window is bigger than the one defined by the ranges, set it to the available range +// +// Sets: +// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain +// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain +// +void MExtractTimeAndChargeDigitalFilter::SetWindowSize(Int_t windowh, Int_t windowl) +{ + + fWindowSizeHiGain = windowh; + fWindowSizeLoGain = windowl; + + const Int_t availhirange = (Int_t)(fHiGainLast-fHiGainFirst+1); + + if (fWindowSizeHiGain > availhirange) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",fWindowSizeHiGain, + " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl; + fHiGainLast = fHiGainFirst + fWindowSizeHiGain; + *fLog << warn << GetDescriptor() + << ": Will set the upper range to: " << (int)fHiGainLast << endl; + } + + if (fWindowSizeHiGain < 2) + { + fWindowSizeHiGain = 2; + *fLog << warn << GetDescriptor() << ": High Gain window size set to two samples" << endl; + } + + if (fLoGainLast != 0 && fWindowSizeLoGain != 0) + { + const Int_t availlorange = (Int_t)(fLoGainLast-fLoGainFirst+1); + + if (fWindowSizeLoGain > availlorange) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",fWindowSizeLoGain, + " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl; + fLoGainLast = fLoGainFirst + fWindowSizeLoGain; + *fLog << warn << GetDescriptor() + << ": Will set the upper range to: " << (int)fLoGainLast << endl; + } + + if (fWindowSizeLoGain<2) + { + fWindowSizeLoGain = 2; + *fLog << warn << GetDescriptor() << ": Low Gain window size set to two samples" << endl; + } + } + + fNumHiGainSamples = (Float_t)fWindowSizeHiGain; + fNumLoGainSamples = (Float_t)fWindowSizeLoGain; + +} + +//---------------------------------------------------------------------------- +// +// Read a pre-defined weights file into the class. +// This is mandatory for the extraction +// +void MExtractTimeAndChargeDigitalFilter::SetWeightsFile(TString filename) +{ + + TFile *f = new TFile(filename); + TH1F *hw_amp = (TH1F*)f->Get("hw_amp"); + TH1F *hw_time = (TH1F*)f->Get("hw_time"); + + if (!hw_amp) + { + *fLog << warn << GetDescriptor() + << ": Could not find hw_amp in " << filename << endl; + return; + } + + if (!hw_time) + { + *fLog << warn << GetDescriptor() + << ": Could not find hw_time in " << filename << endl; + return; + } + + for (int i=0; iGetBinContent(i+1); + fw_time[i] = hw_time->GetBinContent(i+1); + // y[i]=hshape->GetBinContent(i+1); + // derivative[i]=hderivative->GetBinContent(i+1); + + // cout << "w_amp[" << i << "] = " << fw_amp[i] << endl; + } + delete f; +} + +void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeHiGain(Byte_t *ptr, 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) +{ + + const Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1; + + Float_t time_sum = 0.; + Float_t fmax = 0.; + Float_t ftime_max = 0.; + + Float_t pedes = ped.GetPedestal(); + + const Float_t ABoffs = ped.GetPedestalABoffset(); + + Float_t PedMean[2]; + PedMean[0] = pedes + ABoffs; + PedMean[1] = pedes - ABoffs; + + // + // Calculate the sum of the first fWindowSize slices + // + Byte_t *p = ptr; + Byte_t *max_p = ptr; + + while (p < end-fWindowSizeHiGain) + { + sum = 0.; + time_sum = 0.; + + // + // Slide with a window of size fWindowSizeHiGain over the sample + // and multiply the entries with the corresponding weights + // + Byte_t *p1 = p; + for (Int_t sample=0; sample < fWindowSizeHiGain; sample++) + { + const Int_t idx = fBinningResolution*sample+fBinningResolutionHalf; + const Float_t pex = (Float_t)*p1-PedMean[(Int_t(p1-ptr)+abflag) & 0x1]; + sum += fw_amp [idx]*pex; + time_sum += fw_time[idx]*pex; + p1++; + } + + if (sum>fmax) + { + fmax = sum; + ftime_max = time_sum; + max_p = p; + } + p++; + } /* while (pfmax) + { + fmax = sum; + ftime_max = time_sum; + max_p = p; + } + + } +#endif + + if (fmax!=0) + { + ftime_max /= fmax; + Int_t t_iter = Int_t(ftime_max*fBinningResolution); + Int_t sample_iter = 0; + p = max_p; + + while((t_iter)>fBinningResolutionHalf-1 || t_iter<-fBinningResolutionHalf) + { + if (t_iter>fBinningResolutionHalf-1) + { + t_iter -= fBinningResolution; + p--; + sample_iter--; + } + if (t_iter < -fBinningResolutionHalf) + { + t_iter += fBinningResolution; + p++; + sample_iter++; + } + } + + sum = 0.; + Byte_t *p2 = p; + + // + // Slide with a window of size fWindowSizeHiGain over the sample + // and multiply the entries with the corresponding weights + // + for (Int_t sample=0; sample < fWindowSizeHiGain; sample++) + { + const Int_t idx = fBinningResolution*sample + fBinningResolutionHalf + t_iter; + const Float_t pex = (Float_t)*p2-PedMean[(Int_t(p2-ptr)+abflag) & 0x1]; + sum += fw_amp [idx]*pex; + // cout << idx << " " << fw_amp[idx] << " " << pex << " " << t_iter << " " << sum << " " << ((Int_t(p2-ptr)+abflag) & 0x1) << " " << PedMean[(Int_t(p2-ptr)+abflag) & 0x1] << endl; + time_sum += fw_time[idx]*pex; + p2++; + } + + if (sum != 0) + time = (Float_t)(max_p - ptr) + 1. + sample_iter + - (Float_t)t_iter/fBinningResolution - 0.4 - time_sum/sum ; + else + time = 0.; + } /* if (max!=0) */ + else + time=0.; + + return; +} + +void MExtractTimeAndChargeDigitalFilter::FindTimeAndChargeLoGain(Byte_t *ptr, Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag) +{ + + const Byte_t *end = ptr + fLoGainLast - fLoGainFirst + 1; + + Float_t time_sum = 0; + Float_t fmax = 0; + Float_t ftime_max = 0; + + Float_t pedes = ped.GetPedestal(); + + const Float_t ABoffs = ped.GetPedestalABoffset(); + + Float_t PedMean[2]; + PedMean[0] = pedes + ABoffs; + PedMean[1] = pedes - ABoffs; + + // + // Calculate the sum of the first fWindowSize slices + // + sat = 0; + Byte_t *p = ptr; + Byte_t *max_p = ptr; + + while (p < end-fWindowSizeLoGain) + { + sum = 0; + time_sum = 0; + + // + // Slide with a window of size fWindowSizeLoGain over the sample + // and multiply the entries with the corresponding weights + // + Byte_t *p1 = p++; + for (Int_t sample=0; sample < fWindowSizeLoGain; sample++) + { + const Int_t idx = fBinningResolution*sample+fBinningResolutionHalf; + const Int_t adx = (Int_t(p1-ptr)+abflag) & 0x1; + sum += fw_amp [idx]*(*p1-PedMean[adx]); + time_sum += fw_time[idx]*(*p1-PedMean[adx]); + p1++; + } + + if (sum>fmax) + { + fmax = sum; + ftime_max = time_sum; + max_p = p; + } + } /* while (pfBinningResolutionHalf-1 || t_iter<-fBinningResolutionHalf) + { + if (t_iter>fBinningResolutionHalf-1) + { + t_iter -= fBinningResolution; + p--; + sample_iter--; + } + if (t_iter < -fBinningResolutionHalf) + { + t_iter += fBinningResolution; + p++; + sample_iter++; + } + } + + sum = 0; + Byte_t *p2 = p; + + // + // Slide with a window of size fWindowSizeLoGain over the sample + // and multiply the entries with the corresponding weights + // + for (Int_t sample=0; sample < fWindowSizeLoGain; sample++) + { + const Int_t idx = fBinningResolution*sample+fBinningResolutionHalf+t_iter; + const Int_t adx = (Int_t(p2-ptr)+abflag) & 0x1; + sum += fw_amp [idx]*(*p2-PedMean[adx]); + time_sum += fw_time[idx]*(*p2-PedMean[adx]); + p2++; + } + + if (sum != 0) + time = (Float_t)(max_p - ptr) + 1. + sample_iter + - (Float_t)t_iter/fBinningResolution - 0.4 - time_sum/sum ; + else + time = 0.; + } /* if (max!=0) */ + else + time=0.; + + return; +} + +Int_t MExtractTimeAndChargeDigitalFilter::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + + Byte_t hw = fWindowSizeHiGain; + Byte_t lw = fWindowSizeLoGain; + Bool_t rc = kFALSE; + + if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print)) + { + hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw); + rc = kTRUE; + } + if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print)) + { + lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw); + rc = kTRUE; + } + + if (rc) + SetWindowSize(hw, lw); + + if (IsEnvDefined(env, prefix, "BinningResolution", print)) + { + SetBinningResolution(GetEnvValue(env, prefix, "BinningResolution", fBinningResolution)); + rc = kTRUE; + } + + rc = MExtractor::ReadEnv(env, prefix, print) ? kTRUE : rc; + + return rc; +} + Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h (revision 5146) +++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeDigitalFilter.h (revision 5146) @@ -0,0 +1,60 @@ +#ifndef MARS_MExtractTimeAndChargeDigitalFilter +#define MARS_MExtractTimeAndChargeDigitalFilter + +#ifndef MARS_MExtractTimeAndCharge +#include "MExtractTimeAndCharge.h" +#endif + +class MPedestalPix; +class MExtractTimeAndChargeDigitalFilter : public MExtractTimeAndCharge +{ +private: + + static const Byte_t fgHiGainFirst; + static const Byte_t fgHiGainLast; + static const Byte_t fgLoGainFirst; + static const Byte_t fgLoGainLast; + static const Int_t fgHiGainWindowSize; + static const Int_t fgLoGainWindowSize; + static const Int_t fgBinningResolution; + + Int_t fWindowSizeHiGain; + Int_t fWindowSizeLoGain; + + Int_t fBinningResolution; + Int_t fBinningResolutionHalf; + + Float_t fw_amp [60]; + Float_t fw_time[60]; + + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); + +protected: + + virtual void FindTimeAndChargeHiGain(Byte_t *firstused, 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); + virtual void FindTimeAndChargeLoGain(Byte_t *firstused, Float_t &sum, Float_t &dsum, + Float_t &time, Float_t &dtime, + Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag); + +public: + + MExtractTimeAndChargeDigitalFilter(const char *name=NULL, const char *title=NULL); + + void SetWeightsFile(TString filename="pulpo_weights.root"); + void SetWindowSize(Int_t windowh=fgHiGainWindowSize, + Int_t windowl=fgLoGainWindowSize); + + void SetBinningResolution(Int_t r=fgBinningResolution) { + fBinningResolution = r & ~1; + fBinningResolutionHalf = fBinningResolution/2; } + + + ClassDef(MExtractTimeAndChargeDigitalFilter, 0) // Hendrik's digital filter +}; + +#endif + + + Index: /trunk/MagicSoft/Mars/msignal/Makefile =================================================================== --- /trunk/MagicSoft/Mars/msignal/Makefile (revision 5145) +++ /trunk/MagicSoft/Mars/msignal/Makefile (revision 5146) @@ -48,4 +48,6 @@ MExtractTimeFastSpline.cc \ MExtractTimeHighestIntegral.cc \ + MExtractTimeAndCharge.cc \ + MExtractTimeAndChargeDigitalFilter.cc \ MArrivalTime.cc \ MArrivalTimeCalc.cc \ Index: /trunk/MagicSoft/Mars/msignal/SignalLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/msignal/SignalLinkDef.h (revision 5145) +++ /trunk/MagicSoft/Mars/msignal/SignalLinkDef.h (revision 5146) @@ -29,4 +29,8 @@ #pragma link C++ class MExtractTimeHighestIntegral+; +#pragma link C++ class MExtractTimeAndCharge+; +#pragma link C++ class MExtractTimeAndChargeDigitalFilter+; + + #pragma link C++ class MArrivalTimeCam+; #pragma link C++ class MArrivalTimePix+;