Index: /trunk/Mars/mpedestal/MDrsCalibApply.cc =================================================================== --- /trunk/Mars/mpedestal/MDrsCalibApply.cc (revision 14920) +++ /trunk/Mars/mpedestal/MDrsCalibApply.cc (revision 14920) @@ -0,0 +1,260 @@ +/* ======================================================================== *\ +! +! * +! * 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, 10/2006 +! +! Copyright: MAGIC Software Development, 2000-2008 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MDrsCalibApply +// +// This class merges hi- and lo-gain samples into one array and +// subtracts the pedestal (including the AB-offset) from the +// data and stores the result in MDrsCalibApplyedEvt. +// +// Input Containers: +// MRawEvtData +// MRawRunHeader +// MPedestalCam +// +// Output Containers: +// MDrsCalibApplyedEvt +// +////////////////////////////////////////////////////////////////////////////// +#include "MDrsCalibApply.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MArrayB.h" + +#include "MRawRunHeader.h" +#include "MRawEvtData.h" +#include "MRawEvtPixelIter.h" + +#include "MPedestalCam.h" +#include "MPedestalPix.h" + +#include "../mhist/MHDrsCalib.h" +#include "MPedestalSubtractedEvt.h" + +#include "MExtractedSignalCam.h" +#include "MExtractedSignalPix.h" + +ClassImp(MDrsCalibResult); +ClassImp(MDrsCalibApply); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. +// +MDrsCalibApply::MDrsCalibApply(const char *name, const char *title) + : fRawEvt(NULL), fDrsCalib(NULL), fSignal(NULL), fResult(0), + fMaxNumPrevEvents(5), fRemoveSpikes(3), fSlidingAverage(0) +{ + fName = name ? name : "MDrsCalibApply"; + fTitle = title ? title : "Class to subtract pedestal"; +} + +// -------------------------------------------------------------------------- +// +// 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: +// +// - MDrsCalibApplyedEvt +// +Int_t MDrsCalibApply::PreProcess(MParList *pList) +{ + fRawEvt = (MRawEvtData*)pList->FindObject(AddSerialNumber("MRawEvtData")); + if (!fRawEvt) + { + *fLog << err << AddSerialNumber("MRawEvtData") << " not found... aborting." << endl; + return kFALSE; + } + + fDrsCalib = (MDrsCalibration*)pList->FindObject("MDrsCalibration"); + if (!fDrsCalib) + *fLog << warn << "No DRS calibration container (MDrsCalibration) found... skipping." << endl; + + else + *fLog << inf << "DRS calibration will be applied." << endl; + + fResult = 0; + if (fMaxNumPrevEvents>0) + { + fPrevStart.clear(); + fResult = (MDrsCalibResult*)pList->FindCreateObj("MDrsCalibResult"); + if (!fResult) + return kFALSE; + } + + fSignal = (MPedestalSubtractedEvt*)pList->FindCreateObj("MPedestalSubtractedEvt");//, AddSerialNumber(fNamePedestalSubtractedEvt)); + if (!fSignal) + return kFALSE; + + return kTRUE; +} + +Bool_t MDrsCalibApply::ReInit(MParList *pList) +{ + fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader")); + if (!fRunHeader) + { + *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl; + return kFALSE; + } + + if (!fDrsCalib) + return kTRUE; + + if (fDrsCalib->fRoi>0 && fDrsCalib->fRoi != fRunHeader->GetNumSamplesHiGain()) + { + *fLog << warn << "Slice number of DRS calibration " << fDrsCalib->fRoi << " and run-header " << fRunHeader->GetNumSamplesHiGain() << " don't match...\n"; + *fLog << "secondary baseline correction will not be applied." << endl; + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// +Int_t MDrsCalibApply::Process() +{ + // Please note: + // - for data with only hi-gain samples numl is 0 + // - for data with hi- and lo-gain samples + // numl is 0 if read from a raw-data file or a new MC root-file(?) + // numl is not 0 if read from an old MC root-file + + // Total number of samples + const Int_t roi = fRawEvt->GetNumHiGainSamples(); + + // Check if event is empty (presumably MC event -- sanity check) + if (roi==0) + return kCONTINUE; + + // Check for consistency (our simulation can do weird things!) + if (roi!=fRunHeader->GetNumSamplesHiGain()) + { + *fLog << err << "MDrsCalibApply::Process: ERROR - Number of samples in event "; + *fLog << "(hi=" << roi << ")" << endl << " doesn't match number in run-header "; + *fLog << "(" << fRunHeader->GetNumSamplesHiGain() << ")" << endl; + *fLog << " In case you are processing real data you have just found a bug." << endl; + *fLog << " If you process Monte Carlo data, it means probably that the length" << endl; + *fLog << " of the arrays in MRawEvtData are inconsistent with the run-header." << endl; + return kERROR; + } + + if (fRawEvt->GetNumLoGainSamples()+fRunHeader->GetNumSamplesLoGain()>0) + { + *fLog << err << "MDrsCalibApply::Process: ERROR - Data must not contain lo gain slices." << endl; + return kERROR; + } + + // initialize fSignal + fSignal->InitSamples(roi); + + if (!fDrsCalib) // FIXME: Do a fake (none) calibration + return kTRUE; + + const int16_t *val = reinterpret_cast(fRawEvt->GetSamples()); + + const int32_t *offset = fDrsCalib->fOffset.data(); + const int64_t *gain = fDrsCalib->fGain.data(); + const int64_t *trgoff = fDrsCalib->fTrgOff.data(); + + const uint64_t scaleabs = fDrsCalib->fNumOffset; + const uint64_t scalegain = fDrsCalib->fNumGain; + const uint64_t scalerel = fDrsCalib->fNumTrgOff; + + const int16_t *start = reinterpret_cast(fRawEvt->GetStartCells()); + + Float_t *vec = fSignal->GetSamples(0); + + const UShort_t *idx = fRawEvt->GetPixelIds(); + for (size_t ch=0; chGetNumPixels(); ch++) + { + const size_t drs = ch*1024; + const size_t hw = ch*roi; + const size_t sw = idx[ch]*roi; + + DrsCalibrate::ApplyCh(vec+sw, val+hw, start[ch], roi, + offset+drs, scaleabs, + gain +drs, scalegain, + trgoff+hw, scalerel); + } + + if (fResult) + { + fResult->fData.resize(fPrevStart.size()*2); + + deque>::const_iterator it = fPrevStart.begin(); + + for (size_t i=0; ifData[i] = DrsCalibrate::Step(); + + const int32_t off = i%2 ? 3 : roi+10; + + fResult->fData[i] = + DrsCalibrate::CorrectStep(vec, fRawEvt->GetNumPixels(), roi, + it->data(), start, off, idx); + + if (i%2==1) + it++; + } + + fPrevStart.push_front(vector(start, start+fRawEvt->GetNumPixels())); + if (fPrevStart.size()>fMaxNumPrevEvents) + fPrevStart.pop_back(); + } + + switch (fRemoveSpikes) + { + case 1: + DrsCalibrate::RemoveSpikes(vec, roi); + break; + case 2: + DrsCalibrate::RemoveSpikes2(vec, roi); + break; + case 3: + DrsCalibrate::RemoveSpikes3(vec, roi); + break; + } + + + if (fSlidingAverage) + DrsCalibrate::SlidingAverage(vec, roi, fSlidingAverage); + + + return kTRUE; +} Index: /trunk/Mars/mpedestal/MDrsCalibApply.h =================================================================== --- /trunk/Mars/mpedestal/MDrsCalibApply.h (revision 14920) +++ /trunk/Mars/mpedestal/MDrsCalibApply.h (revision 14920) @@ -0,0 +1,65 @@ +#ifndef MARS_MDrsCalibApply +#define MARS_MDrsCalibApply + +#ifndef MARS_MTask +#include "MTask.h" +#endif + +#ifndef MARS_DrsCalib +#include "DrsCalib.h" +#endif + +#include + +class MRawRunHeader; +class MRawEvtData; +class MPedestalSubtractedEvt; +class MDrsCalibration; + +class MDrsCalibResult: public MParContainer +{ + friend class MDrsCalibApply; + +public: + std::vector fData; + + MDrsCalibResult() + { + fName = "MDrsCalibResult"; + } + + const std::vector &GetData() const { return fData; } + + ClassDef(MDrsCalibResult, 0) +}; + +class MDrsCalibApply : public MTask +{ +private: + MRawRunHeader *fRunHeader; //! Run Header + MRawEvtData *fRawEvt; //! Input Raw data + MDrsCalibration *fDrsCalib; //! DRS4 calibration data + MPedestalSubtractedEvt *fSignal; //! Output container + MDrsCalibResult *fResult; //! Output container + + std::deque> fPrevStart; //! History for start cells of previous events + + UShort_t fMaxNumPrevEvents; + UShort_t fRemoveSpikes; + UShort_t fSlidingAverage; + + Int_t PreProcess(MParList *pList); + Bool_t ReInit(MParList *pList); + Int_t Process(); + +public: + MDrsCalibApply(const char *name=NULL, const char *title=NULL); + + void SetMaxNumPrevEvents(UShort_t n) { fMaxNumPrevEvents = n; } + void SetRemoveSpikes(int typ=1) { fRemoveSpikes = typ; } + void SetSlidingAverage(UShort_t w) { fSlidingAverage = w; } + + ClassDef(MDrsCalibApply, 0) // Class to subtract pedestal including AB-noise +}; + +#endif