/* ======================================================================== *\ ! ! * ! * 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): Josep Flix 04/2001 ! Author(s): Thomas Bretz 05/2001 ! Author(s): Sebastian Commichau 12/2003 ! Author(s): Javier Rico 01/2004 ! Author(s): Markus Gaug 01/2004 ! ! Copyright: MAGIC Software Development, 2000-2004 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MPedCalcPedRun // // This task takes a pedestal run file and fills MPedestalCam during // the Process() with the pedestal and rms computed in an event basis. // In the PostProcess() MPedestalCam is finally filled with the pedestal // mean and rms computed in a run basis. // More than one run (file) can be merged // // // Actually, MPedCalcPedRun applies the following formula (1): // // PedRMS = Sqrt( (sum(x_i^2) - sum(x_i)^2/n) / n-1 / 14 ) // // where x_i is the sum of 14 FADC slices and sum means the sum over all // events, n is the number of events. // // For a high number of events, this formula is equivalent to formula (2): // // PedRMS = Sqrt( ( - *) / 14 ) // // where <> is the mean over all events and x_i again the sum over the 14 // slices. // // If you assume statistical equivalence of all slices (say, all have equal // offset and are not correlated and fluctuate Gaussian), it should also be // equivalent to (old formula) (3): // // PedRMS = Sqrt( ( - *) ) // // which is the RMS per slice of a single slice (p_i) and // <> the mean over the total number of measurements, i.e. n*14. // // If we assume that at least our pairs fluctuate independently and Gaussian, // then we can use the actual formula (1) in order to get // fluctuations of pairs by the transformation: // // PedRMS/pair = PedRMS (form. (3)) * Sqrt(2) // // (However, we know that our slice-to-slice fluctuations are not Gaussian // (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: // MPedestalCam // // ///////////////////////////////////////////////////////////////////////////// #include "MPedCalcPedRun.h" #include "MExtractor.h" #include "MParList.h" #include "MLog.h" #include "MLogManip.h" #include "MRawRunHeader.h" #include "MRawEvtPixelIter.h" #include "MRawEvtData.h" #include "MPedestalPix.h" #include "MPedestalCam.h" #include "MGeomPix.h" #include "MGeomCam.h" #include "MBadPixelsPix.h" #include "MBadPixelsCam.h" #include "MGeomCamMagic.h" ClassImp(MPedCalcPedRun); using namespace std; 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: // - all pointers to NULL // - 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) : fWindowSizeHiGain(fgHiGainWindowSize), fWindowSizeLoGain(fgLoGainWindowSize), fGeom(NULL),fBad(NULL) { 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) { fNumSamplesTot = 0; 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); } // -------------------------------------------------------------------------- // // Look for the following input containers: // // - MRawEvtData // - MRawRunHeader // - MGeomCam // - MBadPixelsCam // // The following output containers are also searched and created if // they were not found: // // - MPedestalCam // Int_t MPedCalcPedRun::PreProcess( MParList *pList ) { Clear(); fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); if (!fRawEvt) { *fLog << err << "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; } fGeom = (MGeomCam*)pList->FindObject("MGeomCam"); if (!fGeom) { *fLog << err << "MGeomCam not found... aborting." << endl; return kFALSE; } fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam"); if (!fPedestals) return kFALSE; fBad = (MBadPixelsCam*)pList->FindObject("MBadPixelsCam"); return kTRUE; } // -------------------------------------------------------------------------- // // 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. // // A loop over the MBadPixelsCam is performed and bad pixels are set // to MPedestalPix::SetValid(kFALSE); // Bool_t MPedCalcPedRun::ReInit(MParList *pList) { 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; if (fBad) { const Int_t nbads = fBad->GetSize(); for (Int_t i=0; i<(nbads>npixels?npixels:nbads);i++) if ((*fBad)[i].IsBad()) (*fPedestals)[i].SetValid(kFALSE); } return kTRUE; } // -------------------------------------------------------------------------- // // Fill the MPedestalCam container with the signal mean and rms for the event. // Store the measured signal in arrays fSumx and fSumx2 so that we can // calculate the overall mean and rms in the PostProcess() // Int_t MPedCalcPedRun::Process() { MRawEvtPixelIter pixel(fRawEvt); while (pixel.Next()) { const UInt_t idx = pixel.GetPixelId(); const UInt_t aidx = (*fGeom)[idx].GetAidx(); const UInt_t sector = (*fGeom)[idx].GetSector(); Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst; Byte_t *end = ptr + fWindowSizeHiGain; UInt_t sum = 0; UInt_t sqr = 0; do { sum += *ptr; sqr += *ptr * *ptr; } 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; // // These three lines have been uncommented by Markus Gaug // If anybody needs them, please contact me!! // // const Float_t higainped = msum/fNumHiGainSlices; // const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSlices)/(fNumHiGainSlices-1.)); // (*fPedestals)[idx].Set(higainped, higainrms); fSumx[idx] += msum; fAreaSumx[aidx] += msum; fSectorSumx[sector] += msum; // // The old version: // // const Float_t msqr = (Float_t)sqr; // fSumx2[idx] += msqr; // // The new version: // const Float_t sqrsum = msum*msum; fSumx2[idx] += sqrsum; fAreaSumx2[aidx] += sqrsum; fSectorSumx2[sector] += sqrsum; } fPedestals->SetReadyToSave(); fNumSamplesTot += fWindowSizeHiGain + fWindowSizeLoGain; return kTRUE; } // -------------------------------------------------------------------------- // // Compute signal mean and rms in the whole run and store it in MPedestalCam // Int_t MPedCalcPedRun::PostProcess() { // Compute pedestals and rms from the whole run const ULong_t n = fNumSamplesTot; const ULong_t nevts = GetNumExecutions(); MRawEvtPixelIter pixel(fRawEvt); while (pixel.Next()) { const Int_t pixid = pixel.GetPixelId(); const UInt_t aidx = (*fGeom)[pixid].GetAidx(); const UInt_t sector = (*fGeom)[pixid].GetSector(); fAreaValid [aidx]++; fSectorValid[sector]++; const Float_t sum = fSumx.At(pixid); const Float_t sum2 = fSumx2.At(pixid); const Float_t higainped = sum/n; // // The old version: // // const Float_t higainrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.)); // // The new version: // // 1. Calculate the Variance of the sums: Float_t higainVar = (sum2-sum*sum/nevts)/(nevts-1.); // 2. Scale the variance to the number of slices: higainVar /= (Float_t)(fWindowSizeHiGain+fWindowSizeLoGain); // 3. Calculate the RMS from the Variance: const Float_t higainrms = TMath::Sqrt(higainVar); (*fPedestals)[pixid].Set(higainped, higainrms); } // // Loop over the (two) area indices to get the averaged pedestal per aidx // for (Int_t aidx=0; aidxGetAverageArea(aidx).Set(higainped, higainrms); } // // Loop over the (six) sector indices to get the averaged pedestal per sector // for (Int_t sector=0; sectorGetAverageSector(sector).Set(higainped, higainrms); } fPedestals->SetTotalEntries(fNumSamplesTot); fPedestals->SetReadyToSave(); return kTRUE; }