Index: trunk/MagicSoft/Mars/mpedestal/MExtractPedestal.cc =================================================================== --- trunk/MagicSoft/Mars/mpedestal/MExtractPedestal.cc (revision 5476) +++ trunk/MagicSoft/Mars/mpedestal/MExtractPedestal.cc (revision 5476) @@ -0,0 +1,621 @@ +/* ======================================================================== *\ +! +! * +! * 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): Markus Gaug 01/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MExtractPedestal +// +// Pedestal Extractor base class +// +// Input Containers: +// MRawEvtData +// MRawRunHeader +// MRawEvtHeader +// MGeomCam +// MPedestalCam +// +// Output Containers: +// MPedestalCam +// +// This class should be used for pedestal extractors with the following facilities: +// a) Standardized calculation of AB-noise, mean pedestals and RMS +// b) Standardized treatment of area- and sector averaged pedestals values +// c) Possibility to use a signal extractor to be applied on the pedestals +// d) Possibility to handle two MPedestalCams: one for the signal extractor and +// a second to be filled during the pedestal calculating process. +// +// ad a): Every calculated value is referred to one FADC slice (e.g. Mean pedestal per slice), +// RMS per slice. +// MExtractPedestal applies the following formula (1): +// +// Pedestal per slice = sum(x_i) / n / slices +// PedRMS per slice = Sqrt( ( sum(x_i^2) - sum(x_i)^2/n ) / n-1 / slices ) +// AB-Offset per slice = (sumAB0 - sumAB1) / n / slices +// +// where x_i is the sum of "slices" FADC slices and sum means the sum over all +// events. "n" is the number of events, "slices" is the number of summed FADC samples. +// +// Note that the slice-to-slice fluctuations are not Gaussian, but Poissonian, thus +// asymmetric and they are correlated. +// +// It is important to know that the Pedestal per slice and PedRMS per slice depend +// on the number of used FADC slices, as seen in the following plots: +// +//Begin_Html +/* + +*/ +//End_Html +// +//Begin_Html +/* + +*/ +// +// The plots show the inner and outer pixels, respectivly and have the following meaning: +// +// 1) The calculated mean pedestal per slice (from MPedCalcPedRun) +// 2) The fitted mean pedestal per slice (from MHPedestalCam) +// 3) The calculated pedestal RMS per slice (from MPedCalcPedRun) +// 4) The fitted sigma of the pedestal distribution per slice +// (from MHPedestalCam) +// 5) The relative difference between calculation and histogram fit +// for the mean +// 6) The relative difference between calculation and histogram fit +// for the sigma or RMS, respectively. +// +// The calculated means do not change significantly except for the case of 2 slices, +// however the RMS changes from 5.7 per slice in the case of 2 extracted slices +// to 8.3 per slice in the case of 26 extracted slices. This change is very significant. +// +// ad b) Every calculated value is referred to one FADC slice and one (averaged) pixel, +// (e.g. Mean Pedestal per area index per slice per pixel, etc. ) +// +// MExtractPedestal applies the following formula (2): +// +// Averaged Pedestal per slice = sum(x_i) / n / slices / n_pix +// PedRMS per slice = Sqrt( ( sum(x_i^2) - sum(x_i)^2/n ) / n-1 / slices / n_pix ) +// AB-Offset per slice = (sumAB0 - sumAB1) / n / slices / n_pix +// +// where x_i is the sum of "slices" FADC slices and sum means the sum over all +// events and all concerned pixels. +// "n" is the number of events, "slices" is the number of summed FADC samples and +// "n_pix" is the number of pixels belonging to the specific area index or camera sector. +// +// Calculating these averaged event-by-event values is very important to trace coherent +// fluctuations. An example is given in the following plots: +// +//Begin_Html +/* + +*/ +//End_Html +// +// The plots show the extracted pedestals of the inner pixels (obtained +// with MHPedestalCam), averaged on an event-by-event basis from +// run 13428 with switched off camera LV. +// The meaning of the four plots is: +// +// 1) The distribution of the averaged pedestals +// 2) The averaged pedestals vs. time. +// One can see clearly the oscillation pattern +// 3) The fourier transform of the averaged pedestals vs. time. +// One can see clearly a peak at a certain frequency +// 4) The projection of the fourier components with the non-exponential +// (and therefore significant) outlier. +// +// ad c) Many signal extractors, especially those using a sliding window +// have biases and their resolutions for zero-signals do not agree +// with the pedestal RMS. For the F-Factor method in the calibration +// and the image cleaning, however, both have to be known and measured. +// +// For this reason, a signal extractor can be handed over to the +// pedestal extractor and applied on the pedestal events with the +// function SetExtractor(). +// The results will get stored in an MPedestalCam. +// +// Note that only extractors deriving from MExtractTimeAndCharge +// can be used. +// +// ad d) The signal extractors themselves need a pedestal to be subtracted +// from the FADC slices. +// If the user wishes that the pededestals do not get overwritten by +// the results from the signal extractor, a different named MPedestalCam +// can be created with the function: SetNamePedestalOut(). +// +// See also: MPedestalCam, MPedestalPix, MPedCalcPedRun, MPedCalcFromLoGain +// +///////////////////////////////////////////////////////////////////////////// +#include "MExtractPedestal.h" + +#include "MParList.h" + +#include "MLog.h" +#include "MLogManip.h" + +#include "MRawRunHeader.h" +#include "MRawEvtHeader.h" +#include "MRawEvtPixelIter.h" +#include "MRawEvtData.h" + +#include "MPedestalPix.h" +#include "MPedestalCam.h" + +#include "MGeomPix.h" +#include "MGeomCam.h" + +ClassImp(MExtractPedestal); + +using namespace std; + +const TString MExtractPedestal::fgNamePedestalCam = "MPedestalCam"; +// -------------------------------------------------------------------------- +// +// Default constructor: +// +// Sets: +// - all pointers to NULL +// +// Calls: +// - AddToBranchList("fHiGainPixId"); +// - AddToBranchList("fHiGainFadcSamples"); +// - Clear() +// +MExtractPedestal::MExtractPedestal(const char *name, const char *title) + : fGeom(NULL), fPedestalsIn(NULL), fPedestalsOut(NULL), fExtractor(NULL), + fExtractWinFirst(0), fExtractWinSize(0) +{ + + fName = name ? name : "MExtractPedestal"; + fTitle = title ? title : "Base class to calculate pedestals"; + + AddToBranchList("fHiGainPixId"); + AddToBranchList("fLoGainPixId"); + AddToBranchList("fHiGainFadcSamples"); + AddToBranchList("fLoGainFadcSamples"); + + SetPedestalUpdate(kTRUE); + + SetNamePedestalCamIn(); + SetNamePedestalCamOut(); + SetNumEventsDump(); + + Clear(); +} + +void MExtractPedestal::ResetArrays() +{ + // Reset contents of arrays. + fSumx.Reset(); + fSumx2.Reset(); + fSumAB0.Reset(); + fSumAB1.Reset(); + fAreaSumx.Reset(); + fAreaSumx2.Reset(); + fAreaSumAB0.Reset(); + fAreaSumAB1.Reset(); + fAreaFilled.Reset(); + fAreaValid.Reset(); + fSectorSumx.Reset(); + fSectorSumx2.Reset(); + fSectorSumAB0.Reset(); + fSectorSumAB1.Reset(); + fSectorFilled.Reset(); + fSectorValid.Reset(); + +} + +// -------------------------------------------------------------------------- +// +// Resets Arrays: +// +// Sets: +// - fRawEvt to NULL +// - fRunHeader to NULL +// - fEvtHeader to NULL +// - fPedestalsIn to NULL +// - fPedestalsOut to NULL +// +void MExtractPedestal::Clear(const Option_t *o) +{ + + fRawEvt = NULL; + fRunHeader = NULL; + fEvtHeader = NULL; + fPedestalsIn = NULL; + fPedestalsOut = NULL; + + // If the size is yet set, set the size + if (fSumx.GetSize()>0) + ResetArrays(); + +} + +// -------------------------------------------------------------------------- +// +// Checks: +// - if a window is odd +// +Bool_t MExtractPedestal::SetExtractWindow(UShort_t windowf, UShort_t windows) +{ + + Bool_t rc = kTRUE; + + const Int_t odd = windows & 0x1; + + if (odd) + { + *fLog << warn << GetDescriptor(); + *fLog << " - WARNING: Window size in SetExtraxtWindow has to be even... " + " raising from " << windows << " to " << flush; + windows += 1; + *fLog << windows+1 << "!" << endl; + rc = kFALSE; + } + + if (windows==0) + { + *fLog << warn << GetDescriptor(); + *fLog << " - WARNING: Window size in SetExtraxtWindow has to be > 0... adjusting to 2!" << endl; + windows = 2; + rc = kFALSE; + } + + fExtractWinSize = windows; + fExtractWinFirst = windowf; + fExtractWinLast = fExtractWinFirst+fExtractWinSize-1; + + return rc; +} + +// -------------------------------------------------------------------------- +// +// Look for the following input containers: +// +// - MRawEvtData +// - MRawRunHeader +// - MRawEvtHeader +// - MGeomCam +// +// The following output containers are also searched and created if +// they were not found: +// +// - MPedestalCam with the name fPedContainerName +// +Int_t MExtractPedestal::PreProcess(MParList *pList) +{ + + Clear(); + + fRawEvt = (MRawEvtData*)pList->FindObject(AddSerialNumber("MRawEvtData")); + if (!fRawEvt) + { + *fLog << err << AddSerialNumber("MRawEvtData") << " not found... aborting." << endl; + return kFALSE; + } + + fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader")); + if (!fRunHeader) + { + *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl; + return kFALSE; + } + + fEvtHeader = (MRawEvtHeader*)pList->FindObject(AddSerialNumber("MRawEvtHeader")); + if (!fEvtHeader) + { + *fLog << err << AddSerialNumber("MRawEvtHeader") << " not found... aborting." << endl; + return kFALSE; + } + + fGeom = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam")); + if (!fGeom) + { + *fLog << err << AddSerialNumber("MGeomCam") << " not found... aborting." << endl; + return kFALSE; + } + + fPedestalsIn = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", AddSerialNumber(fNamePedestalCamIn.Data())); + if (!fPedestalsIn) + { + *fLog << err << fNamePedestalCamIn.Data() << " could not be found nor created... aborting" << endl; + return kFALSE; + } + + fPedestalsOut = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", AddSerialNumber(fNamePedestalCamOut.Data())); + if (!fPedestalsOut) + { + *fLog << err << fNamePedestalCamOut.Data() << " could not be found nor created... aborting" << endl; + return kFALSE; + } + + *fLog << inf; + Print(); + + return kTRUE; +} + +// --------------------------------------------------------------------------------- +// +// Sets the size (from MPedestalCam::GetSize() ) and resets the following arrays: +// - fSumx +// - fSumx2 +// - fSumAB0 +// - fSumAB1 +// - fAreaSumx +// - fAreaSumx2 +// - fAreaSumAB0 +// - fAreaSumAB1 +// - fAreaFilled +// - fAreaValid +// - fSectorSumx +// - fSectorSumx2 +// - fSectorSumAB0 +// - fSectorSumAB1 +// - fSectorFilled +// - fSectorValid +// +Bool_t MExtractPedestal::ReInit(MParList *pList) +{ + + // If the size is not yet set, set the size + if (fSumx.GetSize()==0) + { + const Int_t npixels = fPedestalsOut->GetSize(); + const Int_t areas = fPedestalsOut->GetAverageAreas(); + const Int_t sectors = fPedestalsOut->GetAverageSectors(); + + fSumx. Set(npixels); + fSumx2. Set(npixels); + fSumAB0.Set(npixels); + fSumAB1.Set(npixels); + + fAreaSumx. Set(areas); + fAreaSumx2. Set(areas); + fAreaSumAB0.Set(areas); + fAreaSumAB1.Set(areas); + fAreaFilled.Set(areas); + fAreaValid .Set(areas); + + fSectorSumx. Set(sectors); + fSectorSumx2. Set(sectors); + fSectorSumAB0.Set(sectors); + fSectorSumAB1.Set(sectors); + fSectorFilled.Set(sectors); + fSectorValid .Set(sectors); + + for (Int_t i=0; iGetAverageArea(aidx).Set(ped, rms,abOffs,nevts); + +} + +// --------------------------------------------------------------------------------- +// +// Calculates for sector idx "sector" with "nspix" valid pixels: +// +// Ped per slice = sum / nevts / fExtractWinSize / nspix; +// RMS per slice = sqrt { (sum2 - sum*sum/nevts) / (nevts-1) / fExtractWinSize / nspix } +// ABOffset per slice = (fSumAB0[idx] - fSumAB1[idx]) / nevts / fExtractWinSize / nspix; +// +// Stores the results in MPedestalCam::GetAverageSector(sector) +// +void MExtractPedestal::CalcSectorResults(const UInt_t nevts, const UInt_t nspix, const UInt_t sector) +{ + + const Float_t sum = fSectorSumx.At(sector); + const Float_t sum2 = fSectorSumx2.At(sector); + + // 1. Calculate the mean of the sums: + Float_t ped = sum/nevts; + // + // 2. Calculate the Variance of the sums: + // + Float_t var = (sum2-sum*sum/nevts)/(nevts-1.); + // + // 3. Calculate the amplitude of the 150MHz "AB" noise + // + Float_t abOffs = (fSectorSumAB0[sector] - fSectorSumAB1[sector]) / nevts; + // + // 4. Scale the mean, variance and AB-noise to the number of slices: + // + ped /= fExtractWinSize; + var /= fExtractWinSize; + abOffs /= fExtractWinSize; + // + // 5. Scale the mean, variance and AB-noise to the number of pixels: + // + ped /= nspix; + var /= nspix; + abOffs /= nspix; + // + // 6. Calculate the RMS from the Variance: + // + const Float_t rms = var<0 ? 0 : TMath::Sqrt(var); + + fPedestalsOut->GetAverageSector(sector).Set(ped, rms,abOffs,nevts); +} + + + +void MExtractPedestal::Print(Option_t *o) const +{ + + *fLog << GetDescriptor() << ":" << endl; + *fLog << "Name of incoming MPedestalCam Container: " << fNamePedestalCamIn.Data() << endl; + *fLog << "Name of outgoing MPedestalCam Container: " << fNamePedestalCamOut.Data() << endl; + *fLog << "Number events for pedestal calculation: " << fNumEventsDump << endl; + *fLog << "Number events for av. areas calculation: " << fNumAreasDump << endl; + *fLog << "Number events for av. sectors calculation: " << fNumSectorsDump << endl; + *fLog << "Pedestal Update is " << (fPedestalUpdate?"on":"off") << endl; + *fLog << "ExtractWindow from slice " << fExtractWinFirst << " to " << fExtractWinLast << " incl." << endl; +}