Index: trunk/Mars/msignal/MExtractFACT.cc =================================================================== --- trunk/Mars/msignal/MExtractFACT.cc (revision 17830) +++ trunk/Mars/msignal/MExtractFACT.cc (revision 17830) @@ -0,0 +1,262 @@ +/* ======================================================================== *\ +! +! * +! * 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): Thomas Bretz, 2014 +! +! Copyright: MAGIC Software Development, 2014 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MExtractFACT +// +////////////////////////////////////////////////////////////////////////////// +#include "MExtractFACT.h" + +#include + +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MRawRunHeader.h" + +#include "MArrivalTimeCam.h" +#include "MArrivalTimePix.h" + +#include "MExtractedSignalCam.h" +#include "MExtractedSignalPix.h" + +#include "MPedestalSubtractedEvt.h" + +ClassImp(MExtractFACT); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +// Sets: +// - fWindowSizeHiGain and fWindowSizeLoGain to 0 +// - fLoGainStartShift to fgLoGainStartShift +// - fLoGainSwitch to fgLoGainSwitch +// +MExtractFACT::MExtractFACT(const char *name, const char *title) +{ + fName = name ? name : "MExtractFACT"; + fTitle = title ? title : "Base class for signal and time extractors"; +} + +// -------------------------------------------------------------------------- +// +// The PreProcess searches for the following input containers: +// - MRawEvtData +// - MRawRunHeader +// - MPedestalCam +// +// The following output containers are also searched and created if +// they were not found: +// +// - MExtractedSignalCam +// - MArrivalTimeCam +// +Int_t MExtractFACT::PreProcess(MParList *pList) +{ + if (!MExtractTime::PreProcess(pList)) + return kFALSE; + + fSignals = (MExtractedSignalCam*)pList->FindCreateObj("MExtractedSignalCam",AddSerialNumber(fNameSignalCam)); + if (!fSignals) + return kFALSE; + + *fLog << flush << inf; + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// The ReInit calls: +// - MExtractor::ReInit() +// +// Call: +// - MArrivalTimeCam::SetUsedFADCSlices(fHiGainFirst, fHiGainLast, fNumHiGainSamples, +// fLoGainFirst, fLoGainLast, fNumLoGainSamples); +// - InitArrays(); +// +Bool_t MExtractFACT::ReInit(MParList *pList) +{ + if (!MExtractTime::ReInit(pList)) + return kFALSE; + + //if (fSignals) + // fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast, fNumHiGainSamples, + // fLoGainFirst, fLoGainLast, fNumLoGainSamples); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Calculate the integral of the FADC time slices and store them as a new +// pixel in the MArrivalTimeCam container. +// Calculate the integral of the FADC time slices and store them as a new +// pixel in the MExtractedSignalCam container. +// The functions FindTimeAndChargeHiGain() and FindTimeAndChargeLoGain() are +// supposed to extract the signal themselves. +// +Int_t MExtractFACT::Process() +{ + // FIXME: Check ranges + //const UInt_t nroi = fSignal->GetNumSamples(); + const UInt_t npix = fSignal->GetNumPixels(); + + const UInt_t rangehi = fHiGainLast - fHiGainFirst - 1; + + for (UInt_t ch=0; chGetSamples(ch); + + const Float_t *pbeg = ptr+fHiGainFirst+1; + const Float_t *pend = ptr+fHiGainLast-1; + + // Find maximum + const Float_t *pmax = pbeg; + + if (!IsNoiseCalculation()) + { + for (const Float_t *p=pbeg; p*pmax) + pmax = p; + } + + // ------------------------------------------------------------------------- + // Calculate a parabola through this and the surrounding points + // on logarithmic values (that's a gaussian) + // + // Quadratic interpolation + // + // calculate the parameters of a parabula such that + // y(i) = a + b*x(i) + c*x(i)^2 + // y'(i) = b + 2*c*x(i) + // + // y = exp( - (x-k)^2 / s^2 / 2 ) + // + // -2*s^2 * log(y) = x^2 - 2*k*x + k^2 + // + // a = (k/s0)^2/2 + // b = k/s^2 + // c = -1/(2s^2) --> s = sqrt(-1/2c) + // + Float_t max = -1; + if (pmax>pbeg && pmax0 && v2>0 && v3>0) + { + const double y1 = log(v1); + const double y2 = log(v2); + const double y3 = log(v3); + + const double a = y2; + const double b = (y3-y1)/2; + const double c = (y1+y3)/2 - y2; + if (c<0) + { + const double w = -1./(2*c); + const double dx = b*w; + + if (dx>=-1 && dx<=1) + { + max = a + b*dx + c*dx*dx; + //time = dx; + //time += Int_t(pmax-ptr); + } + } + } + } + + Float_t time = -1; + if (max>=0) + { + // -10: maximum search window + pend = std::max(pbeg, pmax-10); + for (;pmax>=pend; pmax--) + if (*pmaxpend && pmax[0]!=pmax[1]) + { + time = (max/2-pmax[0])/(pmax[1]-pmax[0]); + time += Int_t(pmax-ptr); + } + } + + if (time<0) + { + time = gRandom->Uniform(rangehi)+fHiGainFirst+1; + max = pbeg[uint16_t(time)]; + } + + // Now store the result in the corresponding containers + MExtractedSignalPix &pix = (*fSignals)[ch]; + MArrivalTimePix &tix = (*fArrTime)[ch]; + pix.SetExtractedSignal(max); + pix.SetGainSaturation(0); + + tix.SetArrivalTime(time); + tix.SetGainSaturation(0); + } + + fArrTime->SetReadyToSave(); + fSignals->SetReadyToSave(); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// In addition to the resources of the base-class MExtractor: +// MJPedestal.MExtractor.LoGainStartShift: -2.8 +// +Int_t MExtractFACT::ReadEnv(const TEnv &env, TString prefix, Bool_t print) +{ + Bool_t rc = MExtractTime::ReadEnv(env, prefix, print); + /* + if (IsEnvDefined(env, prefix, "LoGainStartShift", print)) + { + fLoGainStartShift = GetEnvValue(env, prefix, "LoGainStartShift", fgLoGainStartShift); + SetLoGainStartShift(fLoGainStartShift); + rc = kTRUE; + } + + if (IsEnvDefined(env, prefix, "LoGainSwitch", print)) + { + fLoGainSwitch = GetEnvValue(env, prefix, "LoGainSwitch", (Int_t)fLoGainSwitch); + rc = kTRUE; + } + */ + return rc; +} + Index: trunk/Mars/msignal/MExtractFACT.h =================================================================== --- trunk/Mars/msignal/MExtractFACT.h (revision 17830) +++ trunk/Mars/msignal/MExtractFACT.h (revision 17830) @@ -0,0 +1,26 @@ +#ifndef MARS_MExtractFACT +#define MARS_MExtractFACT + +#ifndef MARS_MExtractTime +#include "MExtractTime.h" +#endif + +class MExtractFACT : public MExtractTime +{ +private: + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); + +public: + MExtractFACT(const char *name=NULL, const char *title=NULL); + + virtual Bool_t InitArrays(Int_t) { return kTRUE; } + + // For MExtractPedestal + Bool_t ReInit(MParList *pList); + + ClassDef(MExtractFACT, 0) // Time And Charge Extractor Base Class +}; + +#endif Index: trunk/Mars/msignal/MFilterData.cc =================================================================== --- trunk/Mars/msignal/MFilterData.cc (revision 17830) +++ trunk/Mars/msignal/MFilterData.cc (revision 17830) @@ -0,0 +1,146 @@ +/* ======================================================================== *\ +! +! * +! * 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): Thomas Bretz, 2013 +! +! Copyright: MAGIC Software Development, 2013 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MFilterData +// +////////////////////////////////////////////////////////////////////////////// +#include "MFilterData.h" + +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MPedestalSubtractedEvt.h" + +ClassImp(MFilterData); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MFilterData::MFilterData(const char *name, const char *title) + : fSignalIn(NULL), fSignalOut(NULL), + fNameSignalIn("MPedestalSubtractedEvt"), + fNameSignalOut("FilteredEvt") +{ + fName = name ? name : "MFilterData"; + fTitle = title ? title : "Class to filter the data"; +} + +// -------------------------------------------------------------------------- +// +Int_t MFilterData::PreProcess(MParList *pList) +{ + fSignalIn = (MPedestalSubtractedEvt*)pList->FindObject(fNameSignalIn); + if (!fSignalIn) + { + *fLog << err << fNameSignalIn << " [MPedestalSubtractedEvt] not found... aborting." << endl; + return kFALSE; + } + + fSignalOut = (MPedestalSubtractedEvt*)pList->FindCreateObj("MPedestalSibtractedEvt", fNameSignalOut); + if (!fSignalOut) + return kFALSE; + + if (fWeights.size()==0) + { + fWeights.resize(14); + fWeights[ 0] = -0.217305; + fWeights[ 1] = -0.213277; + fWeights[ 2] = -0.193537; + fWeights[ 3] = -0.181686; + fWeights[ 4] = -0.15356; + fWeights[ 5] = -0.129926; + fWeights[ 6] = -0.0792033; + fWeights[ 7] = -0.0219311; + fWeights[ 8] = 0.0550301; + fWeights[ 9] = 0.127364; + fWeights[10] = 0.206711; + fWeights[11] = 0.246864; + fWeights[12] = 0.271012; + fWeights[13] = 0.283444; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// +Int_t MFilterData::Process() +{ + const uint16_t nroi = fSignalIn->GetNumSamples(); + const uint16_t npix = fSignalIn->GetNumPixels(); + const uint16_t nw = fWeights.size(); + //const uint16_t last = nroi-1; + + if (fSignalIn!=fSignalOut) + fSignalOut->InitSamples(nroi); // contains setting to 0 + + const float *begw = fWeights.data(); + const float *endw = fWeights.data()+nw; + + const Float_t *beg = fSignalIn->GetSamples(); + + Float_t *out = fSignalOut->GetSamples(); + for (const Float_t *in=beg; in + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +class MPedestalSubtractedEvt; + +class MFilterData : public MTask +{ +private: + MPedestalSubtractedEvt *fSignalIn; //! Input container + MPedestalSubtractedEvt *fSignalOut; //! Output container + + TString fNameSignalIn; + TString fNameSignalOut; + + std::vector fWeights; + + Int_t PreProcess(MParList *pList); + //Bool_t ReInit(MParList *pList); + Int_t Process(); + +public: + MFilterData(const char *name=NULL, const char *title=NULL); + + void SetNameSignalIn(const char *name="MPedestalSubtractedEvt") { fNameSignalIn=name; } + void SetNameSignalOut(const char *name="MPedestalSubtractedEvt") { fNameSignalOut=name; } + + void SetWeights(const std::vector &w) { fWeights = w; } + void SetAverage(const UInt_t n) { fWeights.assign(n, 1./n); } + + ClassDef(MFilterData, 0) // Class to filter the calibrated data +}; + +#endif