/* ======================================================================== *\ ! ! * ! * 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 // // MPedCalcPedRun 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 ) // // 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. // // The plots have been produced on run 20123. You can reproduce them using // the macro pedestalstudies.C // // Usage of this class: // ==================== // // 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 = 29 // 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 // // // ToDo: // - Take a setup file (ReadEnv-implementation) as input // // // Input Containers: // MRawEvtData // MRawRunHeader // MRawEvtHeader // MGeomCam // // Output Containers: // MPedestalCam // // See also: MPedestalCam, MPedestalPix, MHPedestalCam, MExtractor // ///////////////////////////////////////////////////////////////////////////// #include "MPedCalcPedRun.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" #include "MExtractPedestal.h" #include "MExtractTimeAndCharge.h" ClassImp(MPedCalcPedRun); using namespace std; const UShort_t MPedCalcPedRun::fgExtractWinFirst = 0; const UShort_t MPedCalcPedRun::fgExtractWinSize = 6; const UInt_t MPedCalcPedRun::gkFirstRunWithFinalBits = 45589; // -------------------------------------------------------------------------- // // Default constructor: // // Sets: // - all pointers to NULL // - fWindowSizeHiGain to fgHiGainWindowSize // - fWindowSizeLoGain to fgLoGainWindowSize // // Calls: // - AddToBranchList("fHiGainPixId"); // - AddToBranchList("fHiGainFadcSamples"); // - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast) // MPedCalcPedRun::MPedCalcPedRun(const char *name, const char *title) : fOverlap(0) { fName = name ? name : "MPedCalcPedRun"; fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data"; SetExtractWindow(fgExtractWinFirst, fgExtractWinSize); } // -------------------------------------------------------------------------- // // In case that the variables fExtractLast is greater than the number of // High-Gain FADC samples obtained from the run header, the flag // fOverlap is set to the difference and fExtractLast is set back by the // same number of slices. // void MPedCalcPedRun::CheckExtractionWindow() { const UShort_t lastavailable = fRunHeader->GetNumSamplesHiGain()-1; if (fExtractWinLast <= lastavailable) return; const UShort_t diff = fExtractWinLast - lastavailable; *fLog << warn << endl; *fLog << "Selected ExtractWindow [" << fExtractWinFirst << "," << fExtractWinLast; *fLog << "] ranges out of range [0," << lastavailable << "]." << endl; *fLog << "Using " << diff << " samples from the 'Lo-Gain' slices for the pedestal extraction" << endl; fExtractWinLast -= diff; fOverlap = diff; } // -------------------------------------------------------------------------- // // Set fIsFirstPedRun=kTRUE // Int_t MPedCalcPedRun::PreProcess(MParList *pList) { fUsedEvents = 0; fIsFirstPedRun = kTRUE; return MExtractPedestal::PreProcess(pList); } // -------------------------------------------------------------------------- // // The run type is checked for "kRTPedestal" // and the variable fSkip is set in that case // Bool_t MPedCalcPedRun::ReInit(MParList *pList) { if (!MExtractPedestal::ReInit(pList)) return kFALSE; CheckExtractionWindow(); // If this is the first ped run, the next run (call to ReInit) // is not the first anymore if (fRunHeader->GetRunType()==MRawRunHeader::kRTPedestal) { fIsFirstPedRun = kFALSE; return kTRUE; } // If this is the first call to ReInit (before reading the first file) // nothing should be done if (fIsFirstPedRun) return kTRUE; // In any other case some kind of finaliazation must be done *fLog << inf << "Finalizing pedestal calculations..." << flush; if (!Finalize()) return kFALSE; Reset(); return kTRUE; } // -------------------------------------------------------------------------- // // Return kTRUE (without doing anything) in case that the run type is not // equal to 1 (pedestal run) // // 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() { if (!IsPedBitSet()) { *fLog << err << GetDescriptor() << " - ERROR: IsPedBitSet() returned kFALSE... abort." << endl; return kERROR; } if (fRunHeader->GetRunType() != MRawRunHeader::kRTPedestal) return kTRUE; fUsedEvents++; 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(); Float_t sum = 0.; UInt_t ab0 = 0; UInt_t ab1 = 0; if (fExtractor) CalcExtractor(pixel, sum, (*fPedestalsIn)[idx]); else CalcSums(pixel, sum, ab0, ab1); fSumx[idx] += sum; fAreaSumx[aidx] += sum; fSectorSumx[sector] += sum; const Float_t sqrsum = sum*sum; fSumx2[idx] += sqrsum; fAreaSumx2[aidx] += sqrsum; fSectorSumx2[sector] += sqrsum; fSumAB0[idx] += ab0; fSumAB1[idx] += ab1; fAreaSumAB0[aidx] += ab0; fAreaSumAB1[aidx] += ab1; fSectorSumAB0[aidx] += ab0; fSectorSumAB1[aidx] += ab1; } fPedestalsOut->SetReadyToSave(); return kTRUE; } void MPedCalcPedRun::CalcExtractor(const MRawEvtPixelIter &pixel, Float_t &sum, MPedestalPix &ped) { const Bool_t abflag = pixel.HasABFlag(); Byte_t *first = pixel.GetHiGainSamples() + fExtractWinFirst; Byte_t *logain = pixel.GetLoGainSamples(); Byte_t sat = 0; Float_t dummy; fExtractor->SetNoiseCalculation(kTRUE); fExtractor->FindTimeAndChargeHiGain(first,logain,sum,dummy,dummy,dummy,sat,ped,abflag); fExtractor->SetNoiseCalculation(kFALSE); } void MPedCalcPedRun::CalcSums(const MRawEvtPixelIter &pixel, Float_t &sum, UInt_t &ab0, UInt_t &ab1) { Byte_t *higain = pixel.GetHiGainSamples() + fExtractWinFirst; Byte_t *ptr = higain; Byte_t *end = ptr + fExtractWinSize; const Bool_t abflag = pixel.HasABFlag(); Int_t sumi = 0; Int_t cnt = 0; do { sumi += *ptr; if (!pixel.IsABFlagValid()) continue; const Int_t abFlag = (fExtractWinFirst + abflag + cnt) & 0x1; if (abFlag) ab1 += *ptr; else ab0 += *ptr; cnt++; } while (++ptr != end); if (fOverlap != 0) { ptr = pixel.GetLoGainSamples(); end = ptr + fOverlap; do { sumi += *ptr; if (!pixel.IsABFlagValid()) continue; const Int_t abFlag = (fExtractWinFirst + abflag + cnt) & 0x1; if (abFlag) ab1 += *ptr; else ab0 += *ptr; cnt++; } while (++ptr != end); } sum = (Float_t)sumi; } // -------------------------------------------------------------------------- // // Compute signal mean and rms in the whole run and store it in MPedestalCam // Int_t MPedCalcPedRun::Finalize() { if (fUsedEvents == 0) return kTRUE; MRawEvtPixelIter pixel(fRawEvt); while (pixel.Next()) CalcPixResults(fUsedEvents, pixel.GetPixelId()); // // Loop over the (two) area indices to get the averaged pedestal per aidx // for (UInt_t aidx=0; aidx0) CalcAreaResults(fUsedEvents, fAreaValid[aidx], aidx); // // Loop over the (six) sector indices to get the averaged pedestal per sector // for (UInt_t sector=0; sector0) CalcSectorResults(fUsedEvents, fSectorValid[sector], sector); fPedestalsOut->SetTotalEntries(fUsedEvents*fExtractWinSize); fPedestalsOut->SetReadyToSave(); return kTRUE; } Int_t MPedCalcPedRun::PostProcess() { if (!Finalize()) return kFALSE; return MExtractPedestal::PostProcess(); } //------------------------------------------------------------- // // Return if the pedestal bit was set from the calibration trigger box. // The last but one bit is used for the "pedestal-bit". // // This bit is set since run gkFinalizationTriggerBit // Bool_t MPedCalcPedRun::IsPedBitSet() { if (fRunHeader->GetRunNumber() == 38996) return kTRUE; if (fRunHeader->GetRunNumber() < gkFirstRunWithFinalBits) return kTRUE; return (fEvtHeader->GetTriggerID() & BIT(3)) ? kTRUE : kFALSE; } void MPedCalcPedRun::Print(Option_t *o) const { *fLog << GetDescriptor() << ":" << endl; *fLog << "Name of input MPedestalCam: " << (fPedestalsIn?fPedestalsIn->GetName():fNamePedestalCamIn.Data()) << " (" << fPedestalsIn << ")" << endl; *fLog << "Name of output MPedestalCam: " << (fPedestalsOut?fPedestalsOut->GetName():fNamePedestalCamOut.Data()) << " (" << fPedestalsOut << ")" << endl; *fLog << "ExtractWindow from slice " << fExtractWinFirst << " to " << fExtractWinLast << " incl." << endl; *fLog << "Number overlap lo-gain slices: " << fOverlap << endl; *fLog << "First ped run out of sequence: " << (fIsFirstPedRun?"yes":"no") << endl; *fLog << "Number of used events so far: " << fUsedEvents << endl; if (fExtractor) *fLog << "Extractor used: " << fExtractor->ClassName() << endl; *fLog << endl; }