/* ======================================================================== *\ ! ! * ! * 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 /* */ //End_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" #include "MTaskEnv.h" #include "MExtractTimeAndCharge.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), fPedestalsInter(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"); SetIntermediateStorage( kFALSE ); SetPedestalUpdate ( kTRUE ); SetRandomCalculation ( kTRUE ); SetNamePedestalCamIn(); SetNamePedestalCamOut(); SetNumDump(); 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 // void MExtractPedestal::Clear(const Option_t *o) { fRawEvt = NULL; fRunHeader = NULL; fEvtHeader = 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 && !fExtractor) { *fLog << warn << GetDescriptor(); *fLog << " - WARNING: Window size in SetExtractWindow has to be even... "; *fLog << " raising from " << windows << " to "; windows += 1; *fLog << windows << "!" << endl; rc = kFALSE; } if (windows==0) { *fLog << warn << GetDescriptor(); *fLog << " - WARNING: Window size in SetExtractWindow 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; } if (fExtractor && !fPedestalsIn) { fPedestalsIn = (MPedestalCam*)pList->FindObject(AddSerialNumber(fNamePedestalCamIn), "MPedestalCam"); if (!fPedestalsIn) { *fLog << err << AddSerialNumber(fNamePedestalCamIn) << " not found... aborting." << endl; return kFALSE; } } if (!fPedestalsInter && fIntermediateStorage) { fPedestalsInter = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", AddSerialNumber(fNamePedestalCamInter)); if (!fPedestalsInter) return kFALSE; } if (!fPedestalsOut) { fPedestalsOut = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", AddSerialNumber(fNamePedestalCamOut)); if (!fPedestalsOut) return kFALSE; } *fLog << inf; Print(); return fExtractor ? fExtractor->CallPreProcess(pList) : kTRUE; } Int_t MExtractPedestal::Process() { if (fExtractor) fExtractor->SetNoiseCalculation(fRandomCalculation); const Int_t rc = Calc(); if (fExtractor) fExtractor->SetNoiseCalculation(kFALSE); return rc; } // --------------------------------------------------------------------------------- // // 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->GetNumAverageArea(); const Int_t sectors = fPedestalsOut->GetNumAverageSector(); 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; iReInit(pList)) return kFALSE; if (!fExtractor->InitArrays()) return kFALSE; SetExtractWindow(fExtractor->GetHiGainFirst(), (Int_t)TMath::Nint(fExtractor->GetNumHiGainSamples())); } return kTRUE; } Int_t MExtractPedestal::PostProcess() { fPedestalsIn = NULL; return fExtractor ? fExtractor->CallPostProcess() : kTRUE; } // -------------------------------------------------------------------------- // // The following resources are available: // ExtractWindowFirst: 15 // ExtractWindowSize: 6 // NumEventsDump: 500 // PedestalUpdate: yes // RandomCalculation: yes // Int_t MExtractPedestal::ReadEnv(const TEnv &env, TString prefix, Bool_t print) { Bool_t rc=kFALSE; // find resource for numeventsdump if (IsEnvDefined(env, prefix, "NumEventsDump", print)) { SetNumEventsDump(GetEnvValue(env, prefix, "NumEventsDump", (Int_t)fNumEventsDump)); rc = kTRUE; } // find resource for numeventsdump if (IsEnvDefined(env, prefix, "NumAreasDump", print)) { SetNumAreasDump(GetEnvValue(env, prefix, "NumAreasDump", (Int_t)fNumAreasDump)); rc = kTRUE; } // find resource for numeventsdump if (IsEnvDefined(env, prefix, "NumSectorsDump", print)) { SetNumSectorsDump(GetEnvValue(env, prefix, "NumSectorsDump", (Int_t)fNumSectorsDump)); rc = kTRUE; } // find resource for pedestal update if (IsEnvDefined(env, prefix, "PedestalUpdate", print)) { SetPedestalUpdate(GetEnvValue(env, prefix, "PedestalUpdate", fPedestalUpdate)); rc = kTRUE; } if (IsEnvDefined(env, prefix, "IntermediateStorage", print)) { SetIntermediateStorage(GetEnvValue(env, prefix, "IntermediateStorage", fIntermediateStorage)); rc = kTRUE; } // find resource for random calculation if (IsEnvDefined(env, prefix, "RandomCalculation", print)) { SetRandomCalculation(GetEnvValue(env, prefix, "RandomCalculation", fRandomCalculation)); rc = kTRUE; } // Find resources for ExtractWindow Int_t ef = fExtractWinFirst; Int_t es = fExtractWinSize; if (IsEnvDefined(env, prefix, "ExtractWinFirst", print)) { ef = GetEnvValue(env, prefix, "ExtractWinFirst", ef); rc = kTRUE; } if (IsEnvDefined(env, prefix, "ExtractWinSize", print)) { es = GetEnvValue(env, prefix, "ExtractWinSize", es); rc = kTRUE; } SetExtractWindow(ef,es); // find resource for MPedestalCam if (IsEnvDefined(env, prefix, "NamePedestalCamIn", print)) { SetNamePedestalCamIn(GetEnvValue(env, prefix, "NamePedestalCamIn", fNamePedestalCamIn)); rc = kTRUE; } if (IsEnvDefined(env, prefix, "NamePedestalCamInter", print)) { SetNamePedestalCamInter(GetEnvValue(env, prefix, "NamePedestalCamInter", fNamePedestalCamInter)); rc = kTRUE; } if (IsEnvDefined(env, prefix, "NamePedestalCamOut", print)) { SetNamePedestalCamOut(GetEnvValue(env, prefix, "NamePedestalCamOut", fNamePedestalCamOut)); rc = kTRUE; } return rc; } // --------------------------------------------------------------------------------- // // Calculates for pixel "idx": // // Ped per slice = sum / n / fExtractWinSize; // RMS per slice = sqrt { (sum2 - sum*sum/n) / (n-1) / fExtractWinSize } // ABOffset per slice = (fSumAB0[idx] - fSumAB1[idx]) / n / fExtractWinSize; // // Stores the results in MPedestalCam[pixid] // void MExtractPedestal::CalcPixResults(const UInt_t nevts, const UInt_t pixid) { const Float_t sum = fSumx[pixid]; const Float_t sum2 = fSumx2[pixid]; // 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 = (fSumAB0[pixid] - fSumAB1[pixid]) / nevts; // 4. Scale the mean, variance and AB-noise to the number of slices: ped /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize; var /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize; abOffs /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize; // 5. Calculate the RMS from the Variance: const Float_t rms = var<0 ? 0 : TMath::Sqrt(var); (*fPedestalsOut)[pixid].Set(ped, rms, abOffs, nevts); } // --------------------------------------------------------------------------------- // // Calculates for area idx "aidx" with "napix" valid pixels: // // Ped per slice = sum / nevts / fExtractWinSize / napix; // RMS per slice = sqrt { (sum2 - sum*sum/nevts) / (nevts-1) / fExtractWinSize / napix } // ABOffset per slice = (fSumAB0[idx] - fSumAB1[idx]) / nevts / fExtractWinSize / napix; // // Stores the results in MPedestalCam::GetAverageArea(aidx) // void MExtractPedestal::CalcAreaResults(const UInt_t nevts, const UInt_t napix, const UInt_t aidx) { const Float_t sum = fAreaSumx[aidx]; const Float_t sum2 = fAreaSumx2[aidx]; // 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 = (fAreaSumAB0[aidx] - fAreaSumAB1[aidx]) / nevts; // 4. Scale the mean, variance and AB-noise to the number of slices: ped /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize; var /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize; abOffs /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize; // 5. Scale the mean, variance and AB-noise to the number of pixels: ped /= napix; var /= napix; abOffs /= napix; // 6. Calculate the RMS from the Variance: const Float_t rms = var<0 ? 0 : TMath::Sqrt(var); fPedestalsOut->GetAverageArea(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[sector]; const Float_t sum2 = fSectorSumx2[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 /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize; var /= fExtractor ? fExtractor->GetNumHiGainSamples() : fExtractWinSize; abOffs /= fExtractor ? fExtractor->GetNumHiGainSamples() : 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 input MPedestalCam: " << (fPedestalsIn?fPedestalsIn->GetName():fNamePedestalCamIn.Data()) << " (" << fPedestalsIn << ")" << endl; *fLog << "Name of interm. MPedestalCam: " << (fPedestalsInter?fPedestalsInter->GetName():fNamePedestalCamInter.Data()) << " (" << fPedestalsInter << ")" << endl; *fLog << "Name of output MPedestalCam: " << (fPedestalsOut?fPedestalsOut->GetName():fNamePedestalCamOut.Data()) << " (" << fPedestalsOut << ")" << endl; *fLog << "Intermediate Storage is " << (fIntermediateStorage?"on":"off") << endl; *fLog << "Pedestal Update is " << (fPedestalUpdate?"on":"off") << endl; if (fPedestalUpdate) { *fLog << "Num evts for pedestal calc: " << fNumEventsDump << endl; *fLog << "Num evts for avg.areas calc: " << fNumAreasDump << endl; *fLog << "Num evts for avg.sector calc: " << fNumSectorsDump << endl; } if (fExtractor) { *fLog << "Extractor used: " << fExtractor->ClassName() << " ("; *fLog << (fRandomCalculation?"":"non-") << "random)" << endl; } }