Index: trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc =================================================================== --- trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc (revision 3888) +++ trunk/MagicSoft/Mars/mpedestal/MPedCalcPedRun.cc (revision 3889) @@ -26,5 +26,4 @@ ! \* ======================================================================== */ - ///////////////////////////////////////////////////////////////////////////// // @@ -70,7 +69,27 @@ // (and moreover asymmetric) and that they are also correlated.) // +// Call: SetRange(higainfirst, higainlast, logainfirst, logainlast) +// to modify the ranges in which the window is allowed to move. +// Defaults are: +// +// fHiGainFirst = fgHiGainFirst = 0 +// fHiGainLast = fgHiGainLast = 14 +// fLoGainFirst = fgLoGainFirst = 0 +// fLoGainLast = fgLoGainLast = 14 +// +// Call: SetWindowSize(windowhigain, windowlogain) +// to modify the sliding window widths. Windows have to be an even number. +// In case of odd numbers, the window will be modified. +// +// Defaults are: +// +// fHiGainWindowSize = fgHiGainWindowSize = 14 +// fLoGainWindowSize = fgLoGainWindowSize = 0 +// // // Input Containers: // MRawEvtData +// MRawRunHeader +// MGeomCam // // Output Containers: @@ -80,4 +99,5 @@ ///////////////////////////////////////////////////////////////////////////// #include "MPedCalcPedRun.h" +#include "MExtractor.h" #include "MParList.h" @@ -93,7 +113,4 @@ #include "MPedestalCam.h" -#include "MExtractedSignalPix.h" -#include "MExtractedSignalCam.h" - #include "MGeomPix.h" #include "MGeomCam.h" @@ -105,21 +122,46 @@ using namespace std; -// -------------------------------------------------------------------------- -// -// default constructor +const Byte_t MPedCalcPedRun::fgHiGainFirst = 0; +const Byte_t MPedCalcPedRun::fgHiGainLast = 14; +const Byte_t MPedCalcPedRun::fgLoGainFirst = 0; +const Byte_t MPedCalcPedRun::fgLoGainLast = 14; +const Byte_t MPedCalcPedRun::fgHiGainWindowSize = 14; +const Byte_t MPedCalcPedRun::fgLoGainWindowSize = 0; +// -------------------------------------------------------------------------- +// +// Default constructor: +// +// Sets: +// - fWindowSizeHiGain to fgHiGainWindowSize +// - fWindowSizeLoGain to fgLoGainWindowSize +// +// Calls: +// - AddToBranchList("fHiGainPixId"); +// - AddToBranchList("fHiGainFadcSamples"); +// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast) +// - Clear() // MPedCalcPedRun::MPedCalcPedRun(const char *name, const char *title) - : fRawEvt(NULL), fPedestals(NULL) -{ - fName = name ? name : "MPedCalcPedRun"; - fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data"; - - AddToBranchList("fHiGainPixId"); - AddToBranchList("fHiGainFadcSamples"); - - fNumHiGainSamples = 0; - Clear(); -} - + : fWindowSizeHiGain(fgHiGainWindowSize), + fWindowSizeLoGain(fgLoGainWindowSize) +{ + fName = name ? name : "MPedCalcPedRun"; + fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data"; + + AddToBranchList("fHiGainPixId"); + AddToBranchList("fHiGainFadcSamples"); + + SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast); + Clear(); +} + +// -------------------------------------------------------------------------- +// +// Sets: +// - fNumSamplesTot to 0 +// - fRawEvt to NULL +// - fRunHeader to NULL +// - fPedestals to NULL +// void MPedCalcPedRun::Clear(const Option_t *o) { @@ -128,8 +170,90 @@ fRawEvt = NULL; + fRunHeader = NULL; fPedestals = NULL; } - +// -------------------------------------------------------------------------- +// +// SetRange: +// +// Calls: +// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast); +// - SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain); +// +void MPedCalcPedRun::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast) +{ + + MExtractor::SetRange(hifirst,hilast,lofirst,lolast); + + // + // Redo the checks if the window is still inside the ranges + // + SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain); + +} + + +// -------------------------------------------------------------------------- +// +// Checks: +// - if a window is odd, subtract one +// - if a window is bigger than the one defined by the ranges, set it to the available range +// - if a window is smaller than 2, set it to 2 +// +// Sets: +// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain +// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain +// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples) +// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples) +// +void MPedCalcPedRun::SetWindowSize(Byte_t windowh, Byte_t windowl) +{ + + fWindowSizeHiGain = windowh & ~1; + fWindowSizeLoGain = windowl & ~1; + + if (fWindowSizeHiGain != windowh) + *fLog << warn << GetDescriptor() << ": Hi Gain window size has to be even, set to: " + << int(fWindowSizeHiGain) << " samples " << endl; + + if (fWindowSizeLoGain != windowl) + *fLog << warn << GetDescriptor() << ": Lo Gain window size has to be even, set to: " + << int(fWindowSizeLoGain) << " samples " << endl; + + const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1; + const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1; + + if (fWindowSizeHiGain > availhirange) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s%2i%s",": Hi Gain window size: ",(int)fWindowSizeHiGain, + " is bigger than available range: [",(int)fHiGainFirst,",",(int)fHiGainLast,"]") << endl; + *fLog << warn << GetDescriptor() + << ": Will set window size to: " << (int)availhirange << endl; + fWindowSizeHiGain = availhirange; + } + + if (fWindowSizeLoGain > availlorange) + { + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s%2i%s",": Lo Gain window size: ",(int)fWindowSizeLoGain, + " is bigger than available range: [",(int)fLoGainFirst,",",(int)fLoGainLast,"]") << endl; + *fLog << warn << GetDescriptor() + << ": Will set window size to: " << (int)availlorange << endl; + fWindowSizeLoGain = availlorange; + } + + + fNumHiGainSamples = (Float_t)fWindowSizeHiGain; + fNumLoGainSamples = (Float_t)fWindowSizeLoGain; + + fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples); + fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples); + +} + + + // -------------------------------------------------------------------------- // @@ -137,4 +261,6 @@ // // - MRawEvtData +// - MRawRunHeader +// - MGeomCam // // The following output containers are also searched and created if @@ -155,4 +281,11 @@ } + fRunHeader = (MRawRunHeader*)pList->FindObject(AddSerialNumber("MRawRunHeader")); + if (!fRunHeader) + { + *fLog << err << AddSerialNumber("MRawRunHeader") << " not found... aborting." << endl; + return kFALSE; + } + fGeom = (MGeomCam*)pList->FindObject("MGeomCam"); if (!fGeom) @@ -165,5 +298,5 @@ if (!fPedestals) return kFALSE; - + return kTRUE; } @@ -171,49 +304,88 @@ // -------------------------------------------------------------------------- // -// The ReInit searches for the following input containers: -// - MRawRunHeader -// -// It also initializes the data arrays fSumx and fSumx2 -// (only for the first read file) -// +// The ReInit searches for: +// - MRawRunHeader::GetNumSamplesHiGain() +// - MRawRunHeader::GetNumSamplesLoGain() +// +// In case that the variables fHiGainLast and fLoGainLast are smaller than +// the even part of the number of samples obtained from the run header, a +// warning is given an the range is set back accordingly. A call to: +// - SetRange(fHiGainFirst, fHiGainLast-diff, fLoGainFirst, fLoGainLast) or +// - SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff) +// is performed in that case. The variable diff means here the difference +// between the requested range (fHiGainLast) and the available one. Note that +// the functions SetRange() are mostly overloaded and perform more checks, +// modifying the ranges again, if necessary. +// Bool_t MPedCalcPedRun::ReInit(MParList *pList) { - const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader"); - if (!runheader) - { - *fLog << warn << dbginf; - *fLog << "Warning - cannot check file type, MRawRunHeader not found." << endl; - } - else - if (runheader->IsMonteCarloRun()) - return kTRUE; - - Int_t npixels = fPedestals->GetSize(); - Int_t areas = fPedestals->GetAverageAreas(); - Int_t sectors = fPedestals->GetAverageSectors(); - - if (fSumx.GetSize()==0) - { - fSumx. Set(npixels); - fSumx2.Set(npixels); - - fAreaSumx. Set(areas); - fAreaSumx2.Set(areas); - fAreaValid.Set(areas); - - fSectorSumx. Set(sectors); - fSectorSumx2.Set(sectors); - fSectorValid.Set(sectors); - - fSumx.Reset(); - fSumx2.Reset(); - } - - // Calculate an even number for the hi gain samples to avoid - // biases due to the fluctuation in pedestal from one slice to - // the other one - fNumHiGainSamples = runheader->GetNumSamplesHiGain() & ~1; - - return kTRUE; + + Int_t lastdesired = (Int_t)fHiGainLast; + Int_t lastavailable = (Int_t)fRunHeader->GetNumSamplesHiGain()-1; + + if (lastdesired > lastavailable) + { + const Int_t diff = lastdesired - lastavailable; + *fLog << endl; + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s%2i%s",": Selected Hi Gain FADC Window [", + (int)fHiGainFirst,",",lastdesired, + "] ranges out of the available limits: [0,",lastavailable,"].") << endl; + *fLog << GetDescriptor() << ": Will reduce the upper edge to " << (int)(fHiGainLast - diff) << endl; + SetRange(fHiGainFirst, fHiGainLast-diff, fLoGainFirst, fLoGainLast); + } + + lastdesired = (Int_t)(fLoGainLast); + lastavailable = (Int_t)fRunHeader->GetNumSamplesLoGain()-1; + + if (lastdesired > lastavailable) + { + const Int_t diff = lastdesired - lastavailable; + *fLog << endl; + *fLog << warn << GetDescriptor() + << Form("%s%2i%s%2i%s%2i%s",": Selected Lo Gain FADC Window [", + (int)fLoGainFirst,",",lastdesired, + "] ranges out of the available limits: [0,",lastavailable,"].") << endl; + *fLog << GetDescriptor() << ": Will reduce the upper edge to " << (int)(fLoGainLast - diff) << endl; + SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff); + } + + Int_t npixels = fPedestals->GetSize(); + Int_t areas = fPedestals->GetAverageAreas(); + Int_t sectors = fPedestals->GetAverageSectors(); + + if (fSumx.GetSize()==0) + { + fSumx. Set(npixels); + fSumx2.Set(npixels); + + fAreaSumx. Set(areas); + fAreaSumx2.Set(areas); + fAreaValid.Set(areas); + + fSectorSumx. Set(sectors); + fSectorSumx2.Set(sectors); + fSectorValid.Set(sectors); + + fSumx.Reset(); + fSumx2.Reset(); + } + + if (fWindowSizeHiGain == 0 && fWindowSizeLoGain == 0) + { + *fLog << err << GetDescriptor() + << ": Number of extracted Slices is 0, abort ... " << endl; + return kFALSE; + } + + + *fLog << endl; + *fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeHiGain + << " HiGain FADC samples starting with slice: " << (int)fHiGainFirst << endl; + *fLog << inf << GetDescriptor() << ": Taking " << (int)fWindowSizeLoGain + << " LoGain FADC samples starting with slice: " << (int)fLoGainFirst << endl; + + return kTRUE; + } // -------------------------------------------------------------------------- @@ -235,6 +407,6 @@ const UInt_t sector = (*fGeom)[idx].GetSector(); - Byte_t *ptr = pixel.GetHiGainSamples(); - const Byte_t *end = ptr + fNumHiGainSamples; + Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst; + Byte_t *end = ptr + fWindowSizeHiGain; UInt_t sum = 0; @@ -247,4 +419,19 @@ } while (++ptr != end); + + if (fWindowSizeLoGain != 0) + { + + ptr = pixel.GetLoGainSamples() + fLoGainFirst; + end = ptr + fWindowSizeLoGain; + + do + { + sum += *ptr; + sqr += *ptr * *ptr; + } + while (++ptr != end); + + } const Float_t msum = (Float_t)sum; @@ -254,6 +441,6 @@ // If anybody needs them, please contact me!! // - // const Float_t higainped = msum/fNumHiGainSamples; - // const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSamples)/(fNumHiGainSamples-1.)); + // const Float_t higainped = msum/fNumHiGainSlices; + // const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSlices)/(fNumHiGainSlices-1.)); // (*fPedestals)[idx].Set(higainped, higainrms); @@ -276,6 +463,5 @@ fPedestals->SetReadyToSave(); - fNumSamplesTot += fNumHiGainSamples; - + fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain; return kTRUE; @@ -319,5 +505,5 @@ Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.); // 2. Scale the variance to the number of slices: - higainVar /= (Float_t)fNumHiGainSamples; + higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain); // 3. Calculate the RMS from the Variance: const Float_t higainrms = TMath::Sqrt(higainVar); @@ -344,5 +530,5 @@ Float_t higainVar = (sum2-sum*sum/aevts)/(aevts-1.); // 2. Scale the variance to the number of slices: - higainVar /= (Float_t)fNumHiGainSamples; + higainVar /= fWindowSizeHiGain+fWindowSizeLoGain; // 3. Calculate the RMS from the Variance: Float_t higainrms = TMath::Sqrt(higainVar); @@ -371,5 +557,5 @@ Float_t higainVar = (sum2-sum*sum/sevts)/(sevts-1.); // 2. Scale the variance to the number of slices: - higainVar /= (Float_t)fNumHiGainSamples; + higainVar /= fWindowSizeHiGain+fWindowSizeLoGain; // 3. Calculate the RMS from the Variance: Float_t higainrms = TMath::Sqrt(higainVar);