Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.cc (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.cc (revision 3636) @@ -53,7 +53,16 @@ #include "MLogManip.h" +#include "MCalibrationPix.h" +#include "MCalibrationCam.h" + #include "MHGausEvents.h" +#include "MBadPixelsPix.h" +#include "MBadPixelsCam.h" + #include "MGeomCam.h" +#include "MGeomPix.h" + +#include "MParList.h" ClassImp(MHCalibrationCam); @@ -61,8 +70,12 @@ using namespace std; +const Int_t MHCalibrationCam::fgAverageNbins = 2000; const Int_t MHCalibrationCam::fgPulserFrequency = 500; // -------------------------------------------------------------------------- // // Default Constructor. +// +// Sets: +// - all pointers to NULL // // Initializes and sets owner of: @@ -75,4 +88,5 @@ // MHCalibrationCam::MHCalibrationCam(const char *name, const char *title) + : fGeom(NULL), fBadPixels(NULL), fCam(NULL) { fName = name ? name : "MHCalibrationCam"; @@ -97,4 +111,5 @@ fAverageLoGainSectors->SetOwner(); + SetAverageNbins(); SetPulserFrequency(); } @@ -306,4 +321,472 @@ // -------------------------------------------------------------------------- // +// Gets the pointers to: +// - MGeomCam +// +// Calls SetupHists(const MParList *pList) +// +// Calls Delete-Function of: +// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray +// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas +// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors +// +Bool_t MHCalibrationCam::SetupFill(const MParList *pList) +{ + + fGeom = (MGeomCam*)pList->FindObject("MGeomCam"); + if (!fGeom) + { + *fLog << err << "MGeomCam not found... aborting." << endl; + return kFALSE; + } + + fHiGainArray->Delete(); + fLoGainArray->Delete(); + + fAverageHiGainAreas->Delete(); + fAverageLoGainAreas->Delete(); + + fAverageHiGainSectors->Delete(); + fAverageLoGainSectors->Delete(); + + return SetupHists(pList); +} + + +Bool_t MHCalibrationCam::SetupHists(const MParList *pList) +{ + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Gets or creates the pointers to: +// - MBadPixelsCam +// +// Searches pointer to: +// - MArrivalTimeCam +// +// Initializes, if empty to MArrivalTimeCam::GetSize() for: +// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray +// +// Initializes, if empty to MGeomCam::GetNumAreas() for: +// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas +// +// Initializes, if empty to MGeomCam::GetNumSectors() for: +// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors +// +// Initializes TArray's to MGeomCam::GetNumAreas and MGeomCam::GetNumSectors, respectively +// Fills with number of valid pixels (if !MBadPixelsPix::IsBad()): +// - MHCalibrationCam::fAverageAreaNum[area index] +// - MHCalibrationCam::fAverageSectorNum[area index] +// +// Calls InitializeHists() for every entry in: +// - MHCalibrationCam::fHiGainArray +// - MHCalibrationCam::fAverageHiGainAreas +// - MHCalibrationCam::fAverageHiGainSectors +// +// Sets Titles and Names for the Histograms +// - MHCalibrationCam::fAverageHiGainAreas +// - MHCalibrationCam::fAverageHiGainSectors +// +Bool_t MHCalibrationCam::ReInit(MParList *pList) +{ + + fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam"); + if (!fBadPixels) + return kFALSE; + + + const Int_t npixels = fGeom->GetNumPixels(); + const Int_t nsectors = fGeom->GetNumSectors(); + const Int_t nareas = fGeom->GetNumAreas(); + + // + // The function TArrayF::Set() already sets all entries to 0. + // + fAverageAreaNum. Set(nareas); + fAverageAreaSat. Set(nareas); + fAverageAreaSigma. Set(nareas); + fAverageAreaSigmaErr. Set(nareas); + fAverageAreaRelSigma. Set(nareas); + fAverageAreaRelSigmaErr.Set(nareas); + fAverageSectorNum. Set(nsectors); + + for (Int_t i=0; iGetNumPixels(); + const Int_t nareas = fGeom->GetNumAreas(); + const Int_t nsectors = fGeom->GetNumSectors(); + + if (fHiGainArray->GetEntries() != npixels) + { + gLog << err << "ERROR - Size mismatch... abort." << endl; + return kFALSE; + } + + if (fLoGainArray->GetEntries() != npixels) + { + gLog << err << "ERROR - Size mismatch... abort." << endl; + return kFALSE; + } + + if (fAverageHiGainAreas->GetEntries() != nareas) + { + *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl; + return kFALSE; + } + + if (fAverageLoGainAreas->GetEntries() != nareas) + { + *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl; + return kFALSE; + } + + if (fAverageHiGainSectors->GetEntries() != nsectors) + { + *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl; + return kFALSE; + } + + if (fAverageLoGainSectors->GetEntries() != nsectors) + { + *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl; + return kFALSE; + } + + return FillHists(par, w); +} + +Bool_t MHCalibrationCam::FillHists(const MParContainer *par, const Stat_t w) +{ + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// 1) FinalizeHists() +// 2) FinalizeBadPixels() +// 3) CalcAverageSigma() +// +Bool_t MHCalibrationCam::Finalize() +{ + if (!FinalizeHists()) + return kFALSE; + + FinalizeBadPixels(); + CalcAverageSigma(); + + return kTRUE; +} + +Bool_t MHCalibrationCam::FinalizeHists() +{ + return kTRUE; +} + +void MHCalibrationCam::FinalizeBadPixels() +{ +} + + +// ------------------------------------------------------------- +// +// If MBadPixelsPix::IsBad(): +// - calls MHGausEvents::SetExcluded() +// +// Calls: +// - MHGausEvents::InitBins() +// - MHGausEvents::ChangeHistId(i) +// - MHGausEvents::SetEventFrequency(fPulserFrequency) +// +void MHCalibrationCam::InitHists(MHGausEvents &hist, MBadPixelsPix &bad, const Int_t i) +{ + + if (bad.IsBad()) + hist.SetExcluded(); + + hist.InitBins(); + hist.ChangeHistId(i); + hist.SetEventFrequency(fPulserFrequency); + +} + +void MHCalibrationCam::FitHiGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp) +{ + + for (Int_t i=0; iGetSize(); i++) + { + + MHGausEvents &hist = (*this)[i]; + + if (hist.IsExcluded()) + continue; + + MCalibrationPix &pix = calcam[i]; + MBadPixelsPix &bad = badcam[i]; + + FitHiGainHists(hist,pix,bad,fittyp,osctyp); + + } + + for (Int_t j=0; jGetSize(); j++) + { + + MHGausEvents &hist = GetAverageHiGainArea(j); + MCalibrationPix &pix = calcam.GetAverageArea(j); + MBadPixelsPix &bad = calcam.GetAverageBadArea(j); + + FitHiGainHists(hist,pix,bad,fittyp,osctyp); + } + + + for (Int_t j=0; jGetSize(); j++) + { + + MHGausEvents &hist = GetAverageHiGainSector(j); + MCalibrationPix &pix = calcam.GetAverageSector(j); + MBadPixelsPix &bad = calcam.GetAverageBadSector(j); + + FitHiGainHists(hist,pix,bad,fittyp,osctyp); + } + +} + +void MHCalibrationCam::FitLoGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp) +{ + + for (Int_t i=0; iGetSize(); i++) + { + + MHGausEvents &hist = (*this)(i); + + if (hist.IsExcluded()) + continue; + + MCalibrationPix &pix = calcam[i]; + MBadPixelsPix &bad = badcam[i]; + + FitLoGainHists(hist,pix,bad,fittyp,osctyp); + + } + + for (Int_t j=0; jGetSize(); j++) + { + + MHGausEvents &hist = GetAverageLoGainArea(j); + MCalibrationPix &pix = calcam.GetAverageArea(j); + MBadPixelsPix &bad = calcam.GetAverageBadArea(j); + + FitLoGainHists(hist,pix,bad,fittyp,osctyp); + } + + + for (Int_t j=0; jGetSize(); j++) + { + + MHGausEvents &hist = GetAverageLoGainSector(j); + MCalibrationPix &pix = calcam.GetAverageSector(j); + MBadPixelsPix &bad = calcam.GetAverageBadSector(j); + + FitLoGainHists(hist,pix,bad,fittyp,osctyp); + } +} + +//------------------------------------------------------------ +// +// For all averaged areas, the fitted sigma is multiplied with the square root of +// the number involved pixels +// +void MHCalibrationCam::CalcAverageSigma() +{ + + for (Int_t j=0; jGetNumAreas(); j++) + { + + MCalibrationPix &pix = (*fCam).GetAverageArea(j); + + fAverageAreaSigma[j] = pix.GetSigma () * TMath::Sqrt((Float_t)fAverageAreaNum[j]); + fAverageAreaSigmaErr[j] = pix.GetSigmaErr () * TMath::Sqrt((Float_t)fAverageAreaNum[j]); + + pix.SetSigma (fAverageAreaSigma[j]); + pix.SetSigmaErr(fAverageAreaSigmaErr[j]); + + fAverageAreaRelSigma[j] = fAverageAreaSigma[j] / pix.GetMean(); + + Float_t relsigmaerr = fAverageAreaSigmaErr[j]*fAverageAreaSigmaErr[j] + / (fAverageAreaSigma[j] *fAverageAreaSigma[j] ); + relsigmaerr += pix.GetMeanErr()*pix.GetMeanErr() + / (pix.GetMean() *pix.GetMean() ); + relsigmaerr *= fAverageAreaRelSigma[j]; + fAverageAreaRelSigmaErr[j] = TMath::Sqrt(relsigmaerr); + } +} + +// --------------------------------------------------------------------------- +// +// Returns if the histogram is empty and sets the following flag: +// - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun) +// +// Fits the histograms with a Gaussian, in case of failure +// calls MHGausEvents::RepeatFit(), in case of repeated failure +// calls MHGausEvents::BypassFit() and sets the following flags: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t fittyp ) +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK(). +// In case no, sets the following flags: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t osctyp ) +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// Retrieves the results and store them in MCalibrationPix +// +void MHCalibrationCam::FitHiGainHists(MHGausEvents &hist, + MCalibrationPix &pix, + MBadPixelsPix &bad, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp) +{ + + if (hist.IsEmpty()) + return; + + // + // 2) Fit the Hi Gain histograms with a Gaussian + // + if (!hist.FitGaus()) + // + // 3) In case of failure set the bit Fitted to false and take histogram means and RMS + // + if (!hist.RepeatFit()) + { + hist.BypassFit(); + bad.SetUncalibrated( fittyp ); + } + + // + // 4) Check for oscillations + // + hist.CreateFourierSpectrum(); + + if (!hist.IsFourierSpectrumOK()) + bad.SetUncalibrated( osctyp ); + + // + // 5) Retrieve the results and store them in this class + // + pix.SetHiGainMean ( hist.GetMean() ); + pix.SetHiGainMeanErr ( hist.GetMeanErr() ); + pix.SetHiGainSigma ( hist.GetSigma() ); + pix.SetHiGainSigmaErr ( hist.GetSigmaErr() ); + pix.SetHiGainProb ( hist.GetProb() ); + pix.SetHiGainNumPickup ( hist.GetPickup() ); + +} + + +// --------------------------------------------------------------------------- +// +// Returns if the histogram is empty and sets the following flag: +// - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun) +// +// Fits the histograms with a Gaussian, in case of failure +// calls MHGausEvents::RepeatFit(), in case of repeated failure +// calls MHGausEvents::BypassFit() and sets the following flags: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t fittyp ) +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK(). +// In case no, sets the following flags: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::UncalibratedType_t osctyp ) +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// Retrieves the results and store them in MCalibrationPix +// +void MHCalibrationCam::FitLoGainHists(MHGausEvents &hist, + MCalibrationPix &pix, + MBadPixelsPix &bad, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp) +{ + + if (hist.IsEmpty()) + return; + + + // + // 2) Fit the Hi Gain histograms with a Gaussian + // + if (!hist.FitGaus()) + // + // 3) In case of failure set the bit Fitted to false and take histogram means and RMS + // + if (!hist.RepeatFit()) + { + hist.BypassFit(); + bad.SetUncalibrated( fittyp ); + } + + // + // 4) Check for oscillations + // + hist.CreateFourierSpectrum(); + + if (!hist.IsFourierSpectrumOK()) + bad.SetUncalibrated( osctyp ); + + // + // 5) Retrieve the results and store them in this class + // + pix.SetLoGainMean ( hist.GetMean() ); + pix.SetLoGainMeanErr ( hist.GetMeanErr() ); + pix.SetLoGainSigma ( hist.GetSigma() ); + pix.SetLoGainSigmaErr ( hist.GetSigmaErr() ); + pix.SetLoGainProb ( hist.GetProb() ); + pix.SetLoGainNumPickup ( hist.GetPickup() ); + +} + + + +// -------------------------------------------------------------------------- +// // Dummy, needed by MCamEvent // @@ -343,4 +826,5 @@ for (Int_t i=0; icd(2*(i+1)-1); GetAverageHiGainArea(i).Draw(opt); Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.h (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationCam.h (revision 3636) @@ -21,14 +21,28 @@ #endif +#ifndef MARS_MBadPixelsPix +#include "MBadPixelsPix.h" +#endif + class TText; class TArrayI; class TArrayF; class MHGausEvents; +class MGeomCam; +class MCalibrationCam; +class MCalibrationPix; +class MBadPixelsCam; +class MBadPixelsPix; class MHCalibrationCam : public MH, public MCamEvent { +private: + + static const Int_t fgAverageNbins; //! The default for fAverageNbins (now set to: 2000) + static const Int_t fgPulserFrequency; //! The default for fPulserFrequency (now set to: 500) + protected: - static const Int_t fgPulserFrequency; // The default for fPulserFrequency (now set to: 500) + Int_t fAverageNbins; // Number of bins for the average histograms Int_t fPulserFrequency; // Light pulser frequency @@ -40,4 +54,8 @@ TObjArray *fAverageLoGainSectors; //-> Array of calibration pixels, one per camera sector + MGeomCam *fGeom; //! Camera geometry + MBadPixelsCam *fBadPixels; //! Bad Pixels storage container + MCalibrationCam *fCam; //! Calibration Cam with the results + TArrayI fAverageAreaNum; // Number of pixels in average pixels per area TArrayI fAverageAreaSat; // Number of saturated slices in average pixels per area @@ -48,4 +66,34 @@ TArrayI fAverageSectorNum; // Number of pixels in average pixels per sector + virtual Bool_t SetupHists(const MParList *pList); + virtual Bool_t ReInitHists(MParList *pList); + virtual Bool_t FillHists(const MParContainer *par, const Stat_t w=1); + virtual Bool_t FinalizeHists(); + virtual void FinalizeBadPixels(); + + void InitHists(MHGausEvents &hist, MBadPixelsPix &bad, const Int_t i); + + void FitHiGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp); + + void FitLoGainArrays(MCalibrationCam &calcam, MBadPixelsCam &badcam, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp); + + void FitHiGainHists(MHGausEvents &hist, + MCalibrationPix &pix, + MBadPixelsPix &bad, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp); + + void FitLoGainHists(MHGausEvents &hist, + MCalibrationPix &pix, + MBadPixelsPix &bad, + MBadPixelsPix::UncalibratedType_t fittyp, + MBadPixelsPix::UncalibratedType_t osctyp); + + void CalcAverageSigma(); + void DrawAverageSigma(Bool_t sat, Bool_t inner, Float_t sigma, Float_t sigmaerr, @@ -57,4 +105,5 @@ ~MHCalibrationCam(); + void SetAverageNbins( const Int_t bins=fgAverageNbins ) { fAverageNbins = bins; } void SetPulserFrequency(const Int_t f=fgPulserFrequency) { fPulserFrequency = f; } @@ -76,4 +125,9 @@ MHGausEvents &GetAverageLoGainSector(UInt_t i); const MHGausEvents &GetAverageLoGainSector(UInt_t i) const; + + virtual Bool_t SetupFill(const MParList *pList); + virtual Bool_t ReInit ( MParList *pList); + virtual Bool_t Fill (const MParContainer *par, const Stat_t w=1); + virtual Bool_t Finalize ( ); // Clone Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.cc (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.cc (revision 3636) @@ -127,4 +127,6 @@ #include "MGeomPix.h" +#include "MHGausEvents.h" + #include "MBadPixelsCam.h" #include "MBadPixelsPix.h" @@ -140,5 +142,4 @@ using namespace std; -const Int_t MHCalibrationChargeCam::fgAverageNbins = 4000; const Float_t MHCalibrationChargeCam::fgNumHiGainSaturationLimit = 0.01; const Float_t MHCalibrationChargeCam::fgNumLoGainSaturationLimit = 0.005; @@ -150,16 +151,10 @@ // - all pointers to NULL // -// Initializes and sets owner of (done in MHCalibrationCam): -// - fHiGainArray, fLoGainArray -// - fAverageHiGainAreas, fAverageLoGainAreas -// - fAverageHiGainSectors, fAverageLoGainSectors -// // Initializes: // - fNumHiGainSaturationLimit to fgNumHiGainSaturationLimit // - fNumLoGainSaturationLimit to fgNumLoGainSaturationLimit -// - fPulserFrequency to fgPulserFrequency (in MHCalibrationCam) // MHCalibrationChargeCam::MHCalibrationChargeCam(const char *name, const char *title) - : fCam(NULL), fRawEvt(NULL), fGeom(NULL), fBadPixels(NULL) + : fRawEvt(NULL) { fName = name ? name : "MHCalibrationChargeCam"; @@ -174,12 +169,6 @@ // Gets the pointers to: // - MRawEvtData -// - MGeomCam -// -// Calls Delete-Function of: -// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray -// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas -// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors -// -Bool_t MHCalibrationChargeCam::SetupFill(const MParList *pList) +// +Bool_t MHCalibrationChargeCam::SetupHists(const MParList *pList) { @@ -191,20 +180,4 @@ } - fGeom = (MGeomCam*)pList->FindObject("MGeomCam"); - if (!fGeom) - { - *fLog << err << "MGeomCam not found... aborting." << endl; - return kFALSE; - } - - fHiGainArray->Delete(); - fLoGainArray->Delete(); - - fAverageHiGainAreas->Delete(); - fAverageLoGainAreas->Delete(); - - fAverageHiGainSectors->Delete(); - fAverageLoGainSectors->Delete(); - return kTRUE; } @@ -213,5 +186,4 @@ // // Gets or creates the pointers to: -// - MBadPixelsCam // - MCalibrationChargeCam // @@ -219,19 +191,5 @@ // - MExtractedSignalCam // -// Initializes, if empty to MExtractedSignalCam::GetSize() for: -// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray -// -// Initializes, if empty to MGeomCam::GetNumAreas() for: -// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas -// -// Initializes, if empty to MGeomCam::GetNumSectors() for: -// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors -// -// Initializes TArrays to MGeomCam::GetNumAreas and MGeomCam::GetNumSectors, respectively -// Fills with number of valid pixels (if !MBadPixelsPix::IsBad()): -// - MHCalibrationCam::fAverageAreaNum[area index] -// - MHCalibrationCam::fAverageSectorNum[area index] -// -// Calls InitializeHiGainHists() and InitializeLoGainHists() for every entry in: +// Calls InitializeHists() for every entry in: // - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray // - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas @@ -239,17 +197,12 @@ // // Sets Titles and Names for the Charge Histograms and -// Sets number of bins to fAverageNbins for: +// Sets number of bins to MHCalibrationCam::fAverageNbins for: // - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas // - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors // -Bool_t MHCalibrationChargeCam::ReInit(MParList *pList) +Bool_t MHCalibrationChargeCam::ReInitHists(MParList *pList) { - fBadPixels = (MBadPixelsCam*)pList->FindCreateObj("MBadPixelsCam"); - if (!fBadPixels) - return kFALSE; - - - fCam = (MCalibrationChargeCam*)pList->FindCreateObj("MCalibrationChargeCam"); + fCam = (MCalibrationCam*)pList->FindCreateObj("MCalibrationChargeCam"); if (!fCam) return kFALSE; @@ -262,46 +215,26 @@ } - const Int_t n = signal->GetSize(); + const Int_t npixels = fGeom->GetNumPixels(); const Int_t nsectors = fGeom->GetNumSectors(); const Int_t nareas = fGeom->GetNumAreas(); - // - // The function TArrayF::Set() already sets all entries to 0. - // - fAverageAreaNum. Set(nareas); - fAverageAreaSat. Set(nareas); - fAverageAreaSigma. Set(nareas); - fAverageAreaSigmaErr. Set(nareas); - fAverageAreaRelSigma. Set(nareas); - fAverageAreaRelSigmaErr.Set(nareas); - fAverageSectorNum. Set(nsectors); - - for (Int_t i=0; iGetEntries()==0) { - fHiGainArray->Expand(n); - for (Int_t i=0; iExpand(npixels); + for (Int_t i=0; iGetEntries()==0) { - fLoGainArray->Expand(n); - - for (Int_t i=0; iExpand(npixels); + + for (Int_t i=0; iGetAverageBadArea(j),j); + MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageHiGainArea(j); - InitializeHiGainHists(hist,fCam->GetAverageBadArea(j),j); - hist.GetHGausHist()->SetTitle("Summed FADC slices average HiGain Area Idx "); - hist.SetChargeNbins(fAverageNbins); + hist.SetNbins(fAverageNbins); hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average HiGain Area Idx "); } @@ -340,8 +273,8 @@ MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageLoGainArea(j); - InitializeLoGainHists(hist,fCam->GetAverageBadArea(j),j); + InitHists(hist,fCam->GetAverageBadArea(j),j); hist.GetHGausHist()->SetTitle("Summed FADC slices average LoGain Area Idx "); - hist.SetChargeNbins(fAverageNbins); + hist.SetNbins(fAverageNbins); hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average LoGain Area Idx "); @@ -361,8 +294,8 @@ MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageHiGainSector(j); - InitializeHiGainHists(hist,fCam->GetAverageBadSector(j),j); + InitHists(hist,fCam->GetAverageBadSector(j),j); hist.GetHGausHist()->SetTitle("Summed FADC slices average HiGain Sector "); - hist.SetChargeNbins(fAverageNbins); + hist.SetNbins(fAverageNbins); hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average HiGain Sector "); @@ -382,8 +315,8 @@ MHCalibrationChargePix &hist = (MHCalibrationChargePix&)GetAverageLoGainSector(j); - InitializeLoGainHists(hist,fCam->GetAverageBadSector(j),j); + InitHists(hist,fCam->GetAverageBadSector(j),j); hist.GetHGausHist()->SetTitle("Summed FADC slices average LoGain Sector "); - hist.SetChargeNbins(fAverageNbins); + hist.SetNbins(fAverageNbins); hist.GetHAbsTime()->SetTitle("Absolute Arrival Time average LoGain Sector "); @@ -394,58 +327,12 @@ } -// ------------------------------------------------------------- -// -// If MBadPixelsPix::IsBad(): -// - calls MHCalibrationChargePix::SetExcluded() -// -// Calls: -// - MHCalibrationChargePix::Init() -// - MHCalibrationChargePix::ChangeHistId(i) -// - MHGausEvents::SetEventFrequency(fPulserFrequency) -// -void MHCalibrationChargeCam::InitializeHiGainHists(MHCalibrationChargePix &hist, MBadPixelsPix &bad, const Int_t i) -{ - - if (bad.IsBad()) - { - *fLog << warn << "Excluded pixel: " << i << " from calibration " << endl; - hist.SetExcluded(); - } - - hist.Init(); - hist.ChangeHistId(i); - hist.SetEventFrequency(fPulserFrequency); - -} - -// ------------------------------------------------------------- -// -// If MBadPixelsPix::IsBad(): -// - calls MHCalibrationChargePix::SetExcluded() -// -// Calls: -// - MHCalibrationChargePix::Init() -// - MHCalibrationChargePix::ChangeHistId(i) -// - MHGausEvents::SetEventFrequency(fPulserFrequency) -// -void MHCalibrationChargeCam::InitializeLoGainHists(MHCalibrationChargePix &hist,MBadPixelsPix &bad, const Int_t i) -{ - - if (bad.IsBad()) - hist.SetExcluded(); - - hist.Init(); - hist.ChangeHistId(i); - hist.SetEventFrequency(fPulserFrequency); -} - // -------------------------------------------------------------------------- // // Retrieves from MExtractedSignalCam: +// - first used LoGain FADC slice +// +// Retrieves from MGeomCam: // - number of pixels -// - first used LoGain FADC slice -// -// Retrieves from MGeomCam: // - number of pixel areas // - number of sectors @@ -453,19 +340,17 @@ // For all TObjArray's (including the averaged ones), the following steps are performed: // -// 1) Test size and return kFALSE if not matching -// -// 2) Fill Charges histograms (MHGausEvents::FillHistAndArray()) with: +// 1) Fill Charges histograms (MHGausEvents::FillHistAndArray()) with: // - MExtractedSignalPix::GetExtractedSignalHiGain(); // - MExtractedSignalPix::GetExtractedSignalLoGain(); // -// 3) Set number of saturated slices (MHCalibrationChargePix::SetSaturated()) with: +// 2) Set number of saturated slices (MHCalibrationChargePix::SetSaturated()) with: // - MExtractedSignalPix::GetNumHiGainSaturated(); // - MExtractedSignalPix::GetNumLoGainSaturated(); // -// 4) Fill AbsTime histograms (MHCalibrationChargePix::FillAbsTime()) with: +// 3) Fill AbsTime histograms (MHCalibrationChargePix::FillAbsTime()) with: // - MRawEvtPixelIter::GetIdxMaxHiGainSample(); // - MRawEvtPixelIter::GetIdxMaxLoGainSample(first slice); // -Bool_t MHCalibrationChargeCam::Fill(const MParContainer *par, const Stat_t w) +Bool_t MHCalibrationChargeCam::FillHists(const MParContainer *par, const Stat_t w) { @@ -477,44 +362,8 @@ } - const Int_t npixels = signal->GetSize(); - const Int_t lofirst = signal->GetFirstUsedSliceLoGain(); + const Int_t npixels = fGeom->GetNumPixels(); const Int_t nareas = fGeom->GetNumAreas(); const Int_t nsectors = fGeom->GetNumSectors(); - - if (fHiGainArray->GetEntries() != npixels) - { - *fLog << err << "ERROR - Size mismatch in number of pixels ... abort." << endl; - return kFALSE; - } - - if (fLoGainArray->GetEntries() != npixels) - { - *fLog << err << "ERROR - Size mismatch in number of pixels ... abort." << endl; - return kFALSE; - } - - if (fAverageHiGainAreas->GetEntries() != nareas) - { - *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl; - return kFALSE; - } - - if (fAverageLoGainAreas->GetEntries() != nareas) - { - *fLog << err << "ERROR - Size mismatch in number of areas ... abort." << endl; - return kFALSE; - } - - if (fAverageHiGainSectors->GetEntries() != nsectors) - { - *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl; - return kFALSE; - } - - if (fAverageLoGainSectors->GetEntries() != nsectors) - { - *fLog << err << "ERROR - Size mismatch in number of sectors ... abort." << endl; - return kFALSE; - } + const Int_t lofirst = signal->GetFirstUsedSliceLoGain(); Float_t sumhiarea [nareas], sumloarea [nareas], timehiarea [nareas], timeloarea [nareas]; @@ -526,9 +375,14 @@ { sumhiarea [j] = sumloarea [j] = timehiarea [j] = timeloarea [j] = 0.; + sathiarea [j] = satloarea [j] = 0; + } + + for (UInt_t j=0; j fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries()) + { + pix.SetHiGainSaturation(); + histhi.CreateFourierSpectrum(); + continue; + } + + pix.SetAbsTimeMean ( histhi.GetAbsTimeMean()); + pix.SetAbsTimeRms ( histhi.GetAbsTimeRms() ); + } + + for (Int_t i=0; iGetSize(); i++) + { + + MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(i); MBadPixelsPix &bad = (*fBadPixels)[i]; - - FinalizeHiGainHists(histhi,pix,bad); - - } - - for (Int_t i=0; iGetSize(); i++) - { - - MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)(*this)(i); if (histlo.IsExcluded()) continue; - MCalibrationChargePix &pix = (*fCam)[i]; - MBadPixelsPix &bad = (*fBadPixels)[i]; - - FinalizeLoGainHists(histlo,pix,bad); - - } - + if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries()) + { + *fLog << warn << "Saturated Lo Gain histogram in pixel: " << i << endl; + bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation ); + histlo.CreateFourierSpectrum(); + continue; + } + + MCalibrationChargePix &pix = (MCalibrationChargePix&)(*fCam)[i]; + + if (pix.IsHiGainSaturation()) + { + pix.SetAbsTimeMean ( histlo.GetAbsTimeMean()); + pix.SetAbsTimeRms ( histlo.GetAbsTimeRms() ); + } + } + for (Int_t j=0; jGetSize(); j++) { MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)GetAverageHiGainArea(j); - MCalibrationChargePix &pix = fCam->GetAverageArea(j); - MBadPixelsPix &bad = fCam->GetAverageBadArea(j); - - FinalizeHiGainHists(histhi,pix,bad); + MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageArea(j); + + if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries()) + { + pix.SetHiGainSaturation(); + histhi.CreateFourierSpectrum(); + continue; + } + + pix.SetAbsTimeMean ( histhi.GetAbsTimeMean()); + pix.SetAbsTimeRms ( histhi.GetAbsTimeRms() ); } @@ -696,8 +579,18 @@ MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainArea(j); - MCalibrationChargePix &pix = fCam->GetAverageArea(j); - MBadPixelsPix &bad = fCam->GetAverageBadArea(j); - - FinalizeLoGainHists(histlo,pix,bad); + MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageArea(j); + + if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries()) + { + *fLog << warn << "Saturated Lo Gain histogram in area idx: " << j << endl; + histlo.CreateFourierSpectrum(); + continue; + } + + if (pix.IsHiGainSaturation()) + { + pix.SetAbsTimeMean ( histlo.GetAbsTimeMean()); + pix.SetAbsTimeRms ( histlo.GetAbsTimeRms() ); + } } @@ -706,239 +599,98 @@ MHCalibrationChargePix &histhi = (MHCalibrationChargePix&)GetAverageHiGainSector(j); - MCalibrationChargePix &pix = fCam->GetAverageSector(j); + MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageSector(j); + + if (histhi.GetSaturated() > fNumHiGainSaturationLimit*histhi.GetHGausHist()->GetEntries()) + { + pix.SetHiGainSaturation(); + histhi.CreateFourierSpectrum(); + continue; + } + + pix.SetAbsTimeMean ( histhi.GetAbsTimeMean()); + pix.SetAbsTimeRms ( histhi.GetAbsTimeRms() ); + } + + for (Int_t j=0; jGetSize(); j++) + { + + MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainSector(j); + MCalibrationChargePix &pix = (MCalibrationChargePix&)fCam->GetAverageSector(j); MBadPixelsPix &bad = fCam->GetAverageBadSector(j); - FinalizeHiGainHists(histhi,pix,bad); - } - - for (Int_t j=0; jGetSize(); j++) - { - - MHCalibrationChargePix &histlo = (MHCalibrationChargePix&)GetAverageLoGainSector(j); - MCalibrationChargePix &pix = fCam->GetAverageSector(j); - MBadPixelsPix &bad = fCam->GetAverageBadSector(j); - - FinalizeLoGainHists(histlo,pix,bad); - } - - for (Int_t j=0; jGetNumAreas();j++) - { - - MCalibrationChargePix &pix = fCam->GetAverageArea(j); + if (histlo.GetSaturated() > fNumLoGainSaturationLimit*histlo.GetHGausHist()->GetEntries()) + { + *fLog << warn << "Saturated Lo Gain histogram in sector: " << j << endl; + bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation ); + histlo.CreateFourierSpectrum(); + continue; + } if (pix.IsHiGainSaturation()) - fAverageAreaSat[j]++; - - fAverageAreaSigma[j] = pix.GetSigmaCharge () * TMath::Sqrt((Float_t)fAverageAreaNum[j]); - fAverageAreaSigmaErr[j] = pix.GetSigmaChargeErr () * TMath::Sqrt((Float_t)fAverageAreaNum[j]); - - pix.SetSigmaCharge (fAverageAreaSigma[j]); - pix.SetSigmaChargeErr(fAverageAreaSigmaErr[j]); - - fAverageAreaRelSigma[j] = fAverageAreaSigma[j] / pix.GetMeanCharge(); - - Float_t relsigmaerr = fAverageAreaSigmaErr[j]*fAverageAreaSigmaErr[j] - / (fAverageAreaSigma[j] *fAverageAreaSigma[j] ); - relsigmaerr += pix.GetMeanChargeErr()*pix.GetMeanChargeErr() - / (pix.GetMeanCharge() *pix.GetMeanCharge() ); - relsigmaerr *= fAverageAreaRelSigma[j]; - fAverageAreaRelSigmaErr[j] = TMath::Sqrt(relsigmaerr); - - } - + { + pix.SetAbsTimeMean ( histlo.GetAbsTimeMean()); + pix.SetAbsTimeRms ( histlo.GetAbsTimeRms() ); + } + } + + // + // Perform the fitting for the High Gain (done in MHCalibrationCam) + // + FitHiGainArrays((MCalibrationCam&)(*fCam),(*fBadPixels), + MBadPixelsPix::kHiGainNotFitted, + MBadPixelsPix::kHiGainOscillating); + // + // Perform the fitting for the Low Gain (done in MHCalibrationCam) + // + FitLoGainArrays((MCalibrationCam&)(*fCam),(*fBadPixels), + MBadPixelsPix::kLoGainNotFitted, + MBadPixelsPix::kLoGainOscillating); + return kTRUE; } -// --------------------------------------------------------------------------- -// -// Returns if the histogram is empty and sets the following flag: -// - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun) -// -// Returns if the number of saturated slices (MHCalibrationChargePix::GetSaturated()) is -// higher than the allowed limit (fNumHiGainSaturationLimit*histogram-entries), creates -// the fourier spectrum and sets the following flags: -// - MCalibrationChargePix::SetHiGainSaturation() -// -// Fits the histograms with a Gaussian, in case of failure -// calls MHCalibrationChargePix::RepeatFit(), in case of repeated failure -// calls MHGausEvents::BypassFit() and sets the following flags: -// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) -// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) -// -// Counts the number of pickup events -// -// Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK(). -// In case no, sets the following flags: -// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating ) -// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) -// -// Retrieves the results and store them in MCalibrationChargePix -// -void MHCalibrationChargeCam::FinalizeHiGainHists(MHCalibrationChargePix &hist, - MCalibrationChargePix &pix, - MBadPixelsPix &bad) +// -------------------------------------------------------------------------- +// +// Takes the decisions under which conditions a pixel is declared: +// MBadPixelsPix::kUnsuitableRun or MBadPixelsPix::kUnreliableRun, namely: +// * if MBadPixelsPix::kHiGainNotFitted and !MCalibrationPix::IsHiGainSaturation() +// sets MBadPixelsPix::kUnreliableRun +// * if MBadPixelsPix::kHiGainOscillating and !MCalibrationPix::IsHiGainSaturation() +// sets MBadPixelsPix::kUnreliableRun +// * if MBadPixelsPix::kLoGainNotFitted and MCalibrationPix::IsLoGainSaturation() +// sets MBadPixelsPix::kUnreliableRun +// * if MBadPixelsPix::kLoGainOscillating and MCalibrationPix::IsLoGainSaturation() +// sets MBadPixelsPix::kUnreliableRun +// * if MBadPixelsPix::kLoGainSaturation +// sets MBadPixelsPix::kUnsuitableRun +// +void MHCalibrationChargeCam::FinalizeBadPixels() { - - if (hist.IsEmpty()) - { - *fLog << warn << "Empty Hi Gain histogram in pixel: " << pix.GetPixId() << endl; - bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun); - return; - } - - if (hist.GetSaturated() > fNumHiGainSaturationLimit*hist.GetHGausHist()->GetEntries()) - { - *fLog << warn << "Saturated Hi Gain histogram in pixel: " << pix.GetPixId() << endl; - pix.SetHiGainSaturation(); - hist.CreateFourierSpectrum(); - return; - } - - // - // 2) Fit the Hi Gain histograms with a Gaussian - // - if (!hist.FitGaus()) - // - // 3) In case of failure set the bit Fitted to false and take histogram means and RMS - // - if (!hist.RepeatFit()) - { - hist.BypassFit(); - *fLog << warn << "Hi Gain could not be fitted in pixel: " << pix.GetPixId() << endl; - bad.SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ); - bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); - } - - // - // 4) Check for pickup - // - hist.CountPickup(); - - // - // 5) Check for oscillations - // - hist.CreateFourierSpectrum(); - - if (!hist.IsFourierSpectrumOK()) - { - *fLog << warn << "Oscillating Hi Gain in pixel: " << pix.GetPixId() << endl; - bad.SetUncalibrated( MBadPixelsPix::kHiGainOscillating ); - bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); - } - - // - // 6) Retrieve the results and store them in this class - // - pix.SetHiGainMeanCharge( hist.GetMean() ); - pix.SetHiGainMeanChargeErr( hist.GetMeanErr() ); - pix.SetHiGainSigmaCharge( hist.GetSigma() ); - pix.SetHiGainSigmaChargeErr( hist.GetSigmaErr() ); - pix.SetHiGainChargeProb ( hist.GetProb() ); - - pix.SetAbsTimeMean ( hist.GetAbsTimeMean()); - pix.SetAbsTimeRms ( hist.GetAbsTimeRms() ); - - pix.SetHiGainNumPickup ( hist.GetPickup() ); - + + for (Int_t i=0; iGetSize(); i++) + { + + MBadPixelsPix &bad = (*fBadPixels)[i]; + MCalibrationPix &pix = (*fCam)[i]; + + if (bad.IsUncalibrated( MBadPixelsPix::kHiGainNotFitted )) + if (!pix.IsHiGainSaturation()) + bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); + + if (bad.IsUncalibrated( MBadPixelsPix::kHiGainOscillating )) + bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); + + if (bad.IsUncalibrated( MBadPixelsPix::kLoGainNotFitted )) + if (pix.IsHiGainSaturation()) + bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); + + if (bad.IsUncalibrated( MBadPixelsPix::kLoGainOscillating )) + if (pix.IsHiGainSaturation()) + bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); + + if (bad.IsUncalibrated( MBadPixelsPix::kLoGainSaturation )) + bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun ); + } } - -// --------------------------------------------------------------------------- -// -// Returns if the histogram is empty. If yes and the flag MCalibrationChargePix::IsHiGainSaturation() -// is set, sets the following flag: -// - MBadPixelsPix::SetUnsuitable(MBadPixelsPix::kUnsuitableRun) -// -// Returns if the number of saturated slices (MHCalibrationChargePix::GetSaturated()) is -// higher than the allowed limit (fNumLoGainSaturationLimit*histogram-entries), creates -// the fourier spectrum and sets the following flags: -// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturation ); -// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnsuitableRun ); -// -// Fits the histograms with a Gaussian, in case of failure calls MHCalibrationChargePix::RepeatFit(), -// in case of repeated failure calls MHGausEvents::BypassFit() and sets the following flags: -// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted ) -// If the flag MCalibrationChargePix::IsHiGainSaturation() is set, sets additionally: -// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) -// -// Counts the number of pickup events -// -// Creates the fourier spectrum and tests MHGausEvents::IsFourierSpectrumOK(). -// In case no, sets the following flags: -// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating ) -// If the flag MCalibrationChargePix::IsHiGainSaturation() is set, sets additionally: -// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) -// -// Retrieves the results and store them in MCalibrationChargePix -// -void MHCalibrationChargeCam::FinalizeLoGainHists(MHCalibrationChargePix &hist, - MCalibrationChargePix &pix, - MBadPixelsPix &bad) -{ - - if (hist.IsEmpty()) - { - if (pix.IsHiGainSaturation()) - bad.SetUnsuitable(MBadPixelsPix::kUnsuitableRun); - return; - } - - if (hist.GetSaturated() > fNumLoGainSaturationLimit*hist.GetHGausHist()->GetEntries()) - { - *fLog << warn << "Saturated Lo Gain histogram in pixel: " << pix.GetPixId() << endl; - bad.SetUncalibrated( MBadPixelsPix::kLoGainSaturation ); - bad.SetUnsuitable( MBadPixelsPix::kUnsuitableRun ); - hist.CreateFourierSpectrum(); - return; - } - - // - // 2) Fit the Lo Gain histograms with a Gaussian - // - if (!hist.FitGaus()) - // - // 3) In case of failure set the bit kFitted to false and take histogram means and RMS - // - if (!hist.RepeatFit()) - { - hist.BypassFit(); - bad.SetUncalibrated( MBadPixelsPix::kLoGainNotFitted ); - if (pix.IsHiGainSaturation()) - bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun); - } - - // - // 4) Check for pickup - // - hist.CountPickup(); - - // - // 5) Check for oscillations - // - hist.CreateFourierSpectrum(); - - if (!hist.IsFourierSpectrumOK()) - { - bad.SetUncalibrated( MBadPixelsPix::kLoGainOscillating ); - if (pix.IsHiGainSaturation()) - bad.SetUnsuitable(MBadPixelsPix::kUnreliableRun); - } - // - // 6) Retrieve the results and store them in this class - // - pix.SetLoGainMeanCharge( hist.GetMean() ); - pix.SetLoGainMeanChargeErr( hist.GetMeanErr() ); - pix.SetLoGainSigmaCharge( hist.GetSigma() ); - pix.SetLoGainSigmaChargeErr( hist.GetSigmaErr() ); - pix.SetLoGainChargeProb ( hist.GetProb() ); - - if (pix.IsHiGainSaturation()) - { - pix.SetAbsTimeMean ( hist.GetAbsTimeMean()); - pix.SetAbsTimeRms ( hist.GetAbsTimeRms() ); - } - - pix.SetLoGainNumPickup ( hist.GetPickup() ); - -} // -------------------------------------------------------------------------- Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.h (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeCam.h (revision 3636) @@ -7,6 +7,4 @@ class MRawEvtData; -class MGeomCam; -class MBadPixelsCam; class MBadPixelsPix; class MCalibrationChargeCam; @@ -17,20 +15,18 @@ private: - static const Int_t fgAverageNbins; // The default for fAverageNbins static const Float_t fgNumHiGainSaturationLimit; // The default for fNumHiGainSaturationLimit static const Float_t fgNumLoGainSaturationLimit; // The default for fNumLoGainSaturationLimit - Int_t fAverageNbins; // Number of bins for the average histograms Float_t fNumHiGainSaturationLimit; // Rel. amount sat. higain FADC slices until pixel is called saturated Float_t fNumLoGainSaturationLimit; // Rel. amount sat. logain FADC slices until pixel is called saturated - MCalibrationChargeCam *fCam; //! Calibration Cam with the results - MRawEvtData *fRawEvt; //! Raw event data - MGeomCam *fGeom; //! Camera geometry - MBadPixelsCam *fBadPixels; //! Bad Pixels storage container + MRawEvtData *fRawEvt; //! Raw event data - void InitializeHiGainHists(MHCalibrationChargePix &hist, MBadPixelsPix &bad, const Int_t i); - void InitializeLoGainHists(MHCalibrationChargePix &hist, MBadPixelsPix &bad, const Int_t i); - + Bool_t SetupHists(const MParList *pList); + Bool_t ReInitHists(MParList *pList); + Bool_t FillHists(const MParContainer *par, const Stat_t w=1); + Bool_t FinalizeHists(); + void FinalizeBadPixels(); + void FinalizeHiGainHists(MHCalibrationChargePix &hist, MCalibrationChargePix &pix, MBadPixelsPix &bad); void FinalizeLoGainHists(MHCalibrationChargePix &hist, MCalibrationChargePix &pix, MBadPixelsPix &bad); @@ -41,6 +37,4 @@ ~MHCalibrationChargeCam() {} - - void SetAverageNbins( const Int_t bins=fgAverageNbins ) { fAverageNbins = bins; } void SetNumLoGainSaturationLimit( const Float_t lim=fgNumLoGainSaturationLimit) { fNumLoGainSaturationLimit = lim; } void SetNumHiGainSaturationLimit( const Float_t lim=fgNumHiGainSaturationLimit) { fNumHiGainSaturationLimit = lim; } @@ -49,13 +43,8 @@ Float_t GetNumLoGainSaturationLimit() const { return fNumLoGainSaturationLimit; } - Bool_t SetupFill(const MParList *pList); - Bool_t ReInit ( MParList *pList); - Bool_t Fill (const MParContainer *par, const Stat_t w=1); - Bool_t Finalize ( ); - Bool_t GetPixelContent ( Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; void DrawPixelContent( Int_t num ) const; - ClassDef(MHCalibrationChargeCam, 1) // Histogram class for camera calibration + ClassDef(MHCalibrationChargeCam, 1) // Histogram class for camera charge calibration }; Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeHiGainPix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeHiGainPix.cc (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeHiGainPix.cc (revision 3636) @@ -60,7 +60,7 @@ // // Sets: -// - the default number for fChargeNbins (fgChargeNbins) -// - the default number for fChargeFirst (fgChargeFirst) -// - the default number for fChargeLast (fgChargeLast) +// - the default number for fNbins (fgChargeNbins) +// - the default number for fFirst (fgChargeFirst) +// - the default number for fLast (fgChargeLast) // - the default number for fAbsTimeNbins (fgAbstTimeNbins) // - the default number for fAbsTimeFirst (fgAbsTimeFirst) @@ -79,7 +79,7 @@ fTitle = title ? title : "Fill the FADC sums of the HiGainPix events and perform the fits Pixel "; - SetChargeNbins(); - SetChargeFirst(); - SetChargeLast(); + SetNbins ( fgChargeNbins ); + SetFirst ( fgChargeFirst ); + SetLast ( fgChargeLast ); SetAbsTimeNbins(); @@ -87,8 +87,8 @@ SetAbsTimeLast(); - fHGausHist.SetName("HCalibrationChargeHiGainPix"); + fHGausHist.SetName ("HCalibrationChargeHiGainPix"); fHGausHist.SetTitle("Distribution of Summed Hi Gain FADC slices Pixel "); - fHAbsTime.SetName("HAbsTimeHiGainPix"); + fHAbsTime.SetName ("HAbsTimeHiGainPix"); fHAbsTime.SetTitle("Distribution of Absolute Arrival Times Hi Gain Pixel "); } Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeHiGainPix.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeHiGainPix.h (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeHiGainPix.h (revision 3636) @@ -12,10 +12,10 @@ private: - static const Int_t fgChargeNbins; // Default for fChargeNBins (now set to: 2000 ) - static const Axis_t fgChargeFirst; // Default for fChargeFirst (now set to: -0.5 ) - static const Axis_t fgChargeLast; // Default for fChargeLast (now set to: 1999.5) - static const Int_t fgAbsTimeNbins; // Default for fAbsTimeNbins (now set to: 20 ) - static const Axis_t fgAbsTimeFirst; // Default for fAbsTimeFirst (now set to: -0.5 ) - static const Axis_t fgAbsTimeLast; // Default for fAbsTimeLast (now set to: 19.5 ) + static const Int_t fgChargeNbins; // Default for MHGausEvents::fNbins (now set to: 2000 ) + static const Axis_t fgChargeFirst; // Default for MHGausEvents::fFirst (now set to: -0.5 ) + static const Axis_t fgChargeLast; // Default for MHGausEvents::fLast (now set to: 1999.5) + static const Int_t fgAbsTimeNbins; // Default for fAbsTimeNbins (now set to: 20 ) + static const Axis_t fgAbsTimeFirst; // Default for fAbsTimeFirst (now set to: -0.5 ) + static const Axis_t fgAbsTimeLast; // Default for fAbsTimeLast (now set to: 19.5 ) public: @@ -25,8 +25,4 @@ // Setters - void SetChargeNbins(const Int_t bins =fgChargeNbins) { fChargeNbins = bins; } - void SetChargeFirst(const Axis_t first=fgChargeFirst) { fChargeFirst = first; } - void SetChargeLast (const Axis_t last =fgChargeLast) { fChargeLast = last; } - void SetAbsTimeNbins(const Int_t bins =fgAbsTimeNbins) { fAbsTimeNbins = bins; } void SetAbsTimeFirst(const Axis_t first=fgAbsTimeFirst) { fAbsTimeFirst = first; } Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeLoGainPix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeLoGainPix.cc (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeLoGainPix.cc (revision 3636) @@ -60,7 +60,7 @@ // // Sets: -// - the default number for fChargeNbins (fgChargeNbins) -// - the default number for fChargeFirst (fgChargeFirst) -// - the default number for fChargeLast (fgChargeLast) +// - the default number for fNbins (fgChargeNbins) +// - the default number for fFirst (fgChargeFirst) +// - the default number for fLast (fgChargeLast) // - the default number for fAbsTimeNbins (fgAbstTimeNbins) // - the default number for fAbsTimeFirst (fgAbsTimeFirst) @@ -79,7 +79,7 @@ fTitle = title ? title : "Fill the FADC sums of the Low Gain events and perform the fits Pixel "; - SetChargeNbins(); - SetChargeFirst(); - SetChargeLast(); + SetNbins ( fgChargeNbins ); + SetFirst ( fgChargeFirst ); + SetLast ( fgChargeLast ); SetAbsTimeNbins(); @@ -87,8 +87,8 @@ SetAbsTimeLast(); - fHGausHist.SetName("HCalibrationChargeLoGainPix"); + fHGausHist.SetName ("HCalibrationChargeLoGainPix"); fHGausHist.SetTitle("Distribution of Summed Lo Gain FADC slices Pixel "); - fHAbsTime.SetName("HAbsTimeLoGainPix"); + fHAbsTime.SetName ("HAbsTimeLoGainPix"); fHAbsTime.SetTitle("Distribution of Absolute Arrival Times Lo Gain Pixel "); } Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeLoGainPix.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeLoGainPix.h (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeLoGainPix.h (revision 3636) @@ -11,10 +11,10 @@ private: - static const Int_t fgChargeNbins; // Default for fChargeNBins (now set to: 200 ) - static const Axis_t fgChargeFirst; // Default for fChargeFirst (now set to: -0.5 ) - static const Axis_t fgChargeLast; // Default for fChargeLast (now set to: 199.5 ) - static const Int_t fgAbsTimeNbins; // Default for fAbsTimeNbins (now set to: 15 ) - static const Axis_t fgAbsTimeFirst; // Default for fAbsTimeFirst (now set to: -0.5 ) - static const Axis_t fgAbsTimeLast; // Default for fAbsTimeLast (now set to: 14.5 ) + static const Int_t fgChargeNbins; // Default for MHGausEvents::fNbins (now set to: 200 ) + static const Axis_t fgChargeFirst; // Default for MHGausEvents::fFirst (now set to: -0.5 ) + static const Axis_t fgChargeLast; // Default for MHGausEvents::fLast (now set to: 199.5 ) + static const Int_t fgAbsTimeNbins; // Default for fAbsTimeNbins (now set to: 15 ) + static const Axis_t fgAbsTimeFirst; // Default for fAbsTimeFirst (now set to: -0.5 ) + static const Axis_t fgAbsTimeLast; // Default for fAbsTimeLast (now set to: 14.5 ) public: @@ -24,8 +24,4 @@ // Setters - void SetChargeNbins(const Int_t bins =fgChargeNbins) { fChargeNbins = bins; } - void SetChargeFirst(const Axis_t first=fgChargeFirst) { fChargeFirst = first; } - void SetChargeLast (const Axis_t last =fgChargeLast) { fChargeLast = last; } - void SetAbsTimeNbins(const Int_t bins =fgAbsTimeNbins) { fAbsTimeNbins = bins; } void SetAbsTimeFirst(const Axis_t first=fgAbsTimeFirst) { fAbsTimeFirst = first; } Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargePix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationChargePix.cc (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationChargePix.cc (revision 3636) @@ -56,5 +56,4 @@ const Axis_t MHCalibrationChargePix::fgAbsTimeFirst = -0.5; const Axis_t MHCalibrationChargePix::fgAbsTimeLast = 14.5; -const Float_t MHCalibrationChargePix::fgPickupLimit = 5.; // -------------------------------------------------------------------------- // @@ -62,11 +61,10 @@ // // Sets: -// - the default number for fChargeNbins (fgChargeNbins) -// - the default number for fChargeFirst (fgChargeFirst) -// - the default number for fChargeLast (fgChargeLast) +// - the default number for fNbins (fgChargeNbins) +// - the default number for fFirst (fgChargeFirst) +// - the default number for fLast (fgChargeLast) // - the default number for fAbsTimeNbins (fgAbstTimeNbins) // - the default number for fAbsTimeFirst (fgAbsTimeFirst) // - the default number for fAbsTimeLast (fgAbsTimeLast) -// - the default number for fPickupLimit (fgPickupLimit) // // - the default name of the fHGausHist ("HCalibrationCharge") @@ -88,5 +86,5 @@ // MHCalibrationChargePix::MHCalibrationChargePix(const char *name, const char *title) - : fHAbsTime(), fFlags(0) + : fHAbsTime() { @@ -94,13 +92,11 @@ fTitle = title ? title : "Statistics of the FADC sums of calibration events"; - SetChargeNbins(); - SetChargeFirst(); - SetChargeLast(); + SetNbins ( fgChargeNbins ); + SetFirst ( fgChargeFirst ); + SetLast ( fgChargeLast ); SetAbsTimeNbins(); SetAbsTimeFirst(); SetAbsTimeLast(); - - SetPickupLimit(); fHGausHist.SetName("HCalibrationCharge"); @@ -124,11 +120,10 @@ // // Sets: -// - fHGausHist.SetBins(fChargeNbins,fChargeFirst,fChargeLast); +// - fHGausHist.SetBins(fNbins,fFirst,fLast); // - fHAbsTime.SetBins(fAbsTimeNbins,fAbsTimeFirst,fAbsTimeLast); // -void MHCalibrationChargePix::Init() -{ - - fHGausHist.SetBins(fChargeNbins,fChargeFirst,fChargeLast); +void MHCalibrationChargePix::InitBins() +{ + MHGausEvents::InitBins(); fHAbsTime.SetBins(fAbsTimeNbins,fAbsTimeFirst,fAbsTimeLast); } @@ -137,6 +132,5 @@ // // Sets: -// - all variables to 0., except fPixId to -1 -// - all flags to kFALSE +// - all variables to 0. // // Executes: @@ -146,10 +140,6 @@ { - fPixId = -1; fSaturated = 0.; - fPickup = 0.; - - SetExcluded ( kFALSE ); - + MHGausEvents::Clear(); return; @@ -166,9 +156,7 @@ // -------------------------------------------------------------------------- // -// Set fPixId to id +// Calls MHGausEvents::ChangeHistId() // // Add id to names and titles of: -// - this -// - fHGausHist // - fHAbsTime // @@ -176,28 +164,9 @@ { - fPixId = id; - - fHGausHist.SetName(Form("%s%d", fHGausHist.GetName(), id)); - fHGausHist.SetTitle(Form("%s%d", fHGausHist.GetTitle(), id)); - - fHAbsTime.SetName(Form("%s%d", fHAbsTime.GetName(), id)); + MHGausEvents::ChangeHistId(id); + + fHAbsTime.SetName (Form("%s%d", fHAbsTime.GetName(), id)); fHAbsTime.SetTitle(Form("%s%d", fHAbsTime.GetTitle(), id)); - fName = Form("%s%d", fName.Data(), id); - fTitle = Form("%s%d", fTitle.Data(), id); -} - -// -------------------------------------------------------------------------- -// -// Set Excluded bit from outside -// -void MHCalibrationChargePix::SetExcluded(const Bool_t b) -{ - b ? SETBIT(fFlags,kExcluded) : CLRBIT(fFlags,kExcluded); -} - -const Bool_t MHCalibrationChargePix::IsExcluded() const -{ - return TESTBIT(fFlags,kExcluded); } @@ -317,63 +286,2 @@ } - -// ----------------------------------------------------------------------------- -// -// Repeats the Gauss fit in a smaller range, defined by: -// -// min = GetMean() - fPickupLimit * GetSigma(); -// max = GetMean() + fPickupLimit * GetSigma(); -// -Bool_t MHCalibrationChargePix::RepeatFit(const Option_t *option) -{ - - // - // Get new fitting ranges - // - Axis_t rmin = GetMean() - fPickupLimit * GetSigma(); - Axis_t rmax = GetMean() + fPickupLimit * GetSigma(); - - GetFGausFit()->SetRange(rmin,rmax); - - GetHGausHist()->Fit(GetFGausFit(),option); - - SetMean ( GetFGausFit()->GetParameter(1) ); - SetMeanErr ( GetFGausFit()->GetParameter(2) ); - SetSigma ( GetFGausFit()->GetParError(1) ); - SetSigmaErr ( GetFGausFit()->GetParError(2) ); - SetProb ( GetFGausFit()->GetProb() ); - - // - // The fit result is accepted under condition: - // 1) The results are not nan's - // 2) The NDF is not smaller than fNDFLimit (5) - // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%) - // - if ( TMath::IsNaN ( GetMean() ) - || TMath::IsNaN ( GetMeanErr() ) - || TMath::IsNaN ( GetProb() ) - || TMath::IsNaN ( GetSigma() ) - || TMath::IsNaN ( GetSigmaErr() ) - || GetFGausFit()->GetNDF() < fNDFLimit - || GetProb() < fProbLimit ) - return kFALSE; - - SetGausFitOK(kTRUE); - return kTRUE; - -} - - - -void MHCalibrationChargePix::CountPickup() -{ - fPickup = (Int_t)GetHGausHist()->Integral(GetHGausHist()->GetXaxis()->FindBin(GetMean()+fPickupLimit*GetSigma()), - GetHGausHist()->GetXaxis()->GetLast(), - "width"); -} - - - - - - Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargePix.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationChargePix.h (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationChargePix.h (revision 3636) @@ -12,33 +12,20 @@ private: - static const Int_t fgChargeNbins; // Default for fChargeNBins (now set to: 2000 ) - static const Axis_t fgChargeFirst; // Default for fChargeFirst (now set to: -0.5 ) - static const Axis_t fgChargeLast; // Default for fChargeLast (now set to: 1999.5) - static const Int_t fgAbsTimeNbins; // Default for fAbsTimeNbins (now set to: 15 ) - static const Axis_t fgAbsTimeFirst; // Default for fAbsTimeFirst (now set to: -0.5 ) - static const Axis_t fgAbsTimeLast; // Default for fAbsTimeLast (now set to: 14.5 ) - static const Float_t fgPickupLimit; // Default for fPickupLimit (now set to: 5. ) + static const Int_t fgChargeNbins; // Default for MHGausEvents::fNBins (now set to: 2000 ) + static const Axis_t fgChargeFirst; // Default for MHGausEvents::fFirst (now set to: -0.5 ) + static const Axis_t fgChargeLast; // Default for MHGausEvents::fLast (now set to: 1999.5) + static const Int_t fgAbsTimeNbins; // Default for fAbsTimeNbins (now set to: 15 ) + static const Axis_t fgAbsTimeFirst; // Default for fAbsTimeFirst (now set to: -0.5 ) + static const Axis_t fgAbsTimeLast; // Default for fAbsTimeLast (now set to: 14.5 ) protected: - Int_t fPixId; // The pixel ID - TH1F fHAbsTime; // Histogram containing the absolute arrival times - Int_t fChargeNbins; // Number of bins used for the fHGausHist - Axis_t fChargeFirst; // Lower bound bin used for the fHGausHist - Axis_t fChargeLast; // Upper bound bin used for the fHGausHist - Int_t fAbsTimeNbins; // Number of bins used for the fHAbsTime Axis_t fAbsTimeFirst; // Lower bound bin used for the fHAbsTime Axis_t fAbsTimeLast; // Upper bound bin used for the fHAbsTime - Float_t fPickupLimit; // Upper number of sigmas from the fitted mean above which events are considered as pickup - Float_t fSaturated; // Number of events classified as saturated - Float_t fPickup; // Number of events classified as pick-up - - Byte_t fFlags; // Bit-field for the flags - enum { kExcluded }; // Possible bits to be set public: @@ -49,20 +36,12 @@ virtual void Clear(Option_t *o=""); virtual void Reset(); - virtual void Init(); - virtual void ChangeHistId(Int_t i); + virtual void InitBins(); // Setters - virtual void SetChargeNbins(const Int_t bins =fgChargeNbins) { fChargeNbins = bins; } - virtual void SetChargeFirst(const Axis_t first=fgChargeFirst) { fChargeFirst = first; } - virtual void SetChargeLast( const Axis_t last =fgChargeLast) { fChargeLast = last; } - virtual void SetAbsTimeNbins(const Int_t bins =fgAbsTimeNbins) { fAbsTimeNbins = bins; } virtual void SetAbsTimeFirst(const Axis_t first=fgAbsTimeFirst) { fAbsTimeFirst = first; } virtual void SetAbsTimeLast( const Axis_t last =fgAbsTimeLast) { fAbsTimeLast = last; } - virtual void SetPickupLimit( const Float_t lim =fgPickupLimit) { fPickupLimit = lim; } - void SetSaturated ( const Float_t f ) { fSaturated += f; } - void SetExcluded ( const Bool_t b=kTRUE ); // Getters @@ -70,26 +49,18 @@ const TH1F *GetHAbsTime() const { return &fHAbsTime; } - const Float_t GetIntegral() const; - const Float_t GetAbsTimeMean( ) const; const Float_t GetAbsTimeRms() const; - + const Float_t GetIntegral() const; const Float_t GetSaturated() const { return fSaturated; } - const Float_t GetPickup() const { return fPickup; } - - const Bool_t IsExcluded() const; // Fill histos Bool_t FillAbsTime(const Float_t t); - // Fits - Bool_t RepeatFit(const Option_t *option="RQ0"); - // Draws virtual void Draw(Option_t *opt=""); // Miscelleaneous - void CountPickup(); - + void ChangeHistId(Int_t id); + ClassDef(MHCalibrationChargePix, 1) // Base Histogram class for a Charge Calibration Pixel }; Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.cc (revision 3636) @@ -22,13 +22,60 @@ ! \* ======================================================================== */ - ///////////////////////////////////////////////////////////////////////////// // // -// MHCalibrationRelTimeCam // +// MHCalibrationRelTimeCam // // // -// Hold the CalibrationRelTime information for all pixels in the camera // -// // +// Fills the extracted relative arrival times of MArrivalTimeCam into +// the MHGausEvents-classes MHCalibrationRelTimePix for every: +// +// - pixel, stored in the TObjArray's MHCalibrationCam::fHiGainArray +// +// - average pixel per area index (e.g. inner and outer for the MAGIC camera), +// stored in the TObjArray's MHCalibrationCam::fAverageHiGainAreas +// +// - average pixel per camera sector (e.g. sectors 1-6 for the MAGIC camera), +// stored in the TObjArray's MHCalibrationCam::fAverageHiGainSectors +// +// Every time is filled into a histogram and an array, in order to perform +// a Fourier analysis (see MHGausEvents). The signals are moreover averaged on an +// event-by-event basis and written into the corresponding average pixels. +// +// The histograms are fitted to a Gaussian, mean and sigma with its errors +// and the fit probability are extracted. If none of these values are NaN's and +// if the probability is bigger than fProbLimit (default: 0.5%), the fit is valid. +// Otherwise, the fit is repeated within ranges of the previous mean +- 5 sigma. +// In case this does not make the fit valid, the histogram means and RMS's are +// taken directly and the following flags are set: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kRelTimeNotFitted ) and +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// Outliers of more than MHCalibrationChargePix::fPickupLimit (default: 5) sigmas +// from the mean are counted as PickUp events (stored in MHCalibrationRelTimePix::fPickup) +// +// The class also fills arrays with the signal vs. event number, creates a fourier +// spectrum and investigates if the projected fourier components follow an exponential +// distribution. In case that the probability of the exponential fit is less than +// fProbLimit (default: 0.5%), the following flags are set: +// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kRelTimeOscillating ) and +// - MBadPixelsPix::SetUnsuitable( MBadPixelsPix::kUnreliableRun ) +// +// This same procedure is performed for the average pixels. +// +// The following results are written into MCalibrationRelTimeCam: +// +// - MCalibrationPix::SetMean() +// - MCalibrationPix::SetMeanErr() +// - MCalibrationPix::SetSigma() +// - MCalibrationPix::SetSigmaErr() +// - MCalibrationPix::SetProb() +// - MCalibrationPix::SetNumPickup() +// +// For all averaged areas, the fitted sigma is multiplied with the square root of +// the number involved pixels in order to be able to compare it to the average of +// sigmas in the camera. +// ///////////////////////////////////////////////////////////////////////////// #include "MHCalibrationRelTimeCam.h" +#include "MHCalibrationRelTimePix.h" #include "MLog.h" @@ -37,17 +84,22 @@ #include "MParList.h" -#include "MHCalibrationRelTimePix.h" +#include "MCalibrationRelTimeCam.h" +#include "MCalibrationPix.h" #include "MArrivalTimeCam.h" #include "MArrivalTimePix.h" +#include "MGeomCam.h" +#include "MGeomPix.h" + +#include "MBadPixelsCam.h" +#include "MBadPixelsPix.h" + ClassImp(MHCalibrationRelTimeCam); using namespace std; -const Float_t MHCalibrationRelTimeCam::fgTimeSliceWidth = 3.3; -const Int_t MHCalibrationRelTimeCam::fgPulserFrequency = 500; // -------------------------------------------------------------------------- // -// Default constructor. Creates a MHCalibrationRelTimePix object for each pixel +// Default Constructor. // MHCalibrationRelTimeCam::MHCalibrationRelTimeCam(const char *name, const char *title) @@ -55,97 +107,160 @@ fName = name ? name : "MHCalibrationRelTimeCam"; - fTitle = title ? title : "Histogram container for the relative time calibration of the camera"; - - // - // loop over all Pixels and create two histograms - // one for the Low and one for the High gain - // connect all the histogram with the container fHist - // - fArray = new TObjArray; - fArray->SetOwner(); - - - SetTimeSliceWidth(); - SetPulserFrequency(); + fTitle = title ? title : "Histogram class for the relative time calibration of the camera"; + } // -------------------------------------------------------------------------- // -// Delete the array conatining the pixel pedest information -// -MHCalibrationRelTimeCam::~MHCalibrationRelTimeCam() -{ - delete fArray; -} - -// -------------------------------------- -// -void MHCalibrationRelTimeCam::Clear(Option_t *o) -{ - - fArray->ForEach(TObject, Clear)(); - return; -} - -// -------------------------------------------------------------------------- -// -// Get i-th pixel (pixel number) -// -MHCalibrationRelTimePix &MHCalibrationRelTimeCam::operator[](UInt_t i) -{ - return *(MHCalibrationRelTimePix*)(fArray->At(i)); -} - -// -------------------------------------------------------------------------- -// -// Get i-th pixel (pixel number) -// -const MHCalibrationRelTimePix &MHCalibrationRelTimeCam::operator[](UInt_t i) const -{ - return *(MHCalibrationRelTimePix*)(fArray->At(i)); -} - - -// -------------------------------------------------------------------------- -// -// Our own clone function is necessary since root 3.01/06 or Mars 0.4 -// I don't know the reason -// -TObject *MHCalibrationRelTimeCam::Clone(const char *) const -{ - - const Int_t n = fArray->GetSize(); - - // - // FIXME, this might be done faster and more elegant, by direct copy. - // - MHCalibrationRelTimeCam *cam = new MHCalibrationRelTimeCam; - - cam->fArray->Expand(n); - - for (int i=0; ifArray)[i]; - (*cam->fArray)[i] = (*fArray)[i]->Clone(); - } - - return cam; -} - - - -// -------------------------------------------------------------------------- -// -// -Bool_t MHCalibrationRelTimeCam::SetupFill(const MParList *pList) -{ +// Gets or creates the pointers to: +// - MCalibrationRelTimeCam +// +// Searches pointer to: +// - MArrivalTimeCam +// +// Initializes, if empty to MArrivalTimeCam::GetSize() for: +// - MHCalibrationCam::fHiGainArray, MHCalibrationCam::fLoGainArray +// +// Initializes, if empty to MGeomCam::GetNumAreas() for: +// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas +// +// Initializes, if empty to MGeomCam::GetNumSectors() for: +// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors +// +// Calls InitializeHists() for every entry in: +// - MHCalibrationCam::fHiGainArray +// - MHCalibrationCam::fAverageHiGainAreas +// - MHCalibrationCam::fAverageHiGainSectors +// +// Sets Titles and Names for the Histograms +// - MHCalibrationCam::fAverageHiGainAreas +// - MHCalibrationCam::fAverageHiGainSectors +// +// Sets number of bins to MHCalibrationCam::fAverageNbins for: +// - MHCalibrationCam::fAverageHiGainAreas, MHCalibrationCam::fAverageLoGainAreas +// - MHCalibrationCam::fAverageHiGainSectors, MHCalibrationCam::fAverageLoGainSectors +// +Bool_t MHCalibrationRelTimeCam::ReInitHists(MParList *pList) +{ + + fCam = (MCalibrationCam*)pList->FindCreateObj("MCalibrationRelTimeCam"); + if (!fCam) + return kFALSE; + + MArrivalTimeCam *signal = (MArrivalTimeCam*)pList->FindObject("MArrivalTimeCam"); + if (!signal) + { + *fLog << err << "MArrivalTimeCam not found... abort." << endl; + return kFALSE; + } + + const Int_t npixels = fGeom->GetNumPixels(); + const Int_t nsectors = fGeom->GetNumSectors(); + const Int_t nareas = fGeom->GetNumAreas(); + + if (fHiGainArray->GetEntries()==0) + { + fHiGainArray->Expand(npixels); + for (Int_t i=0; iGetEntries()==0) + { + fLoGainArray->Expand(npixels); + for (Int_t i=0; iGetEntries()==0) + { + fAverageHiGainAreas->Expand(nareas); + + for (Int_t j=0; jGetAverageBadArea(j),j); + + GetAverageHiGainArea(j).GetHGausHist()->SetTitle("Relative Arrival Times average Area Idx "); + GetAverageHiGainArea(j).SetNbins(fAverageNbins); + } + } + + if (fAverageLoGainAreas->GetEntries()==0) + { + fAverageLoGainAreas->Expand(nareas); + + for (Int_t j=0; jGetAverageBadArea(j),j); + + GetAverageLoGainArea(j).GetHGausHist()->SetTitle("Relative Arrival Times average Area Idx "); + GetAverageLoGainArea(j).SetNbins(fAverageNbins); + } + } + + if (fAverageHiGainSectors->GetEntries()==0) + { + fAverageHiGainSectors->Expand(nsectors); + + for (Int_t j=0; jSetTitle("Relative Arrival Times average Sector "); + GetAverageHiGainSector(j).SetNbins(fAverageNbins); + } + } + + if (fAverageLoGainSectors->GetEntries()==0) + { + fAverageLoGainSectors->Expand(nsectors); + + for (Int_t j=0; jGetAverageBadSector(j),j); + + GetAverageLoGainSector(j).GetHGausHist()->SetTitle("Relative Arrival Times average Sector "); + GetAverageLoGainSector(j).SetNbins(fAverageNbins); + } + } return kTRUE; - -} - -// -------------------------------------------------------------------------- -// -Bool_t MHCalibrationRelTimeCam::Fill(const MParContainer *par, const Stat_t w) +} + + +//-------------------------------------------------------------------------------- +// +// Retrieves pointer to MArrivalTimeCam: +// +// Retrieves from MGeomCam: +// - number of pixels +// - number of pixel areas +// - number of sectors +// +// Fill RelTimes histograms (MHGausEvents::FillHistAndArray()) with: +// - MArrivalTimePix::GetArrivalTime(pixid) - MArrivalTimePix::GetArrivalTime(1); +// (i.e. the time difference between pixel i and pixel 1 (hardware number: 2) +// +Bool_t MHCalibrationRelTimeCam::FillHists(const MParContainer *par, const Stat_t w) { @@ -157,79 +272,128 @@ } - const Int_t n = arrtime->GetSize(); - - if (fArray->GetEntries()==0) - { - fArray->Expand(n); + const Int_t npixels = fGeom->GetNumPixels(); + const Int_t nareas = fGeom->GetNumAreas(); + const Int_t nsectors = fGeom->GetNumSectors(); + + Float_t sumareahi [nareas], sumarealo [nareas]; + Float_t sumsectorhi[nsectors], sumsectorlo[nsectors]; + Int_t numareahi [nareas], numarealo [nareas]; + Int_t numsectorhi[nsectors], numsectorlo[nsectors]; + + for (UInt_t j=0; jGetEntries() != n) - { - gLog << err << "ERROR - Size mismatch... abort." << endl; - return kFALSE; - } - - const MArrivalTimePix &refpix = (*arrtime)[1]; - const Float_t reftime = refpix.IsLoGainUsed() ? refpix.GetArrivalTimeLoGain() : refpix.GetArrivalTimeHiGain(); - - for (Int_t i=0; iGetSize(); i++) +// -------------------------------------------------------------------------- +// +// Calls: +// - MHCalibrationCam::FitHiGainArrays() with flags: +// MBadPixelsPix::kRelTimeNotFitted and MBadPixelsPix::kRelTimeOscillating +// - MHCalibrationCam::FitLoGainArrays() with flags: +// MBadPixelsPix::kRelTimeNotFitted and MBadPixelsPix::kRelTimeOscillating +// +Bool_t MHCalibrationRelTimeCam::FinalizeHists() +{ + + FitHiGainArrays((MCalibrationCam&)(*fCam),*fBadPixels, + MBadPixelsPix::kRelTimeNotFitted, + MBadPixelsPix::kRelTimeOscillating); + + FitLoGainArrays((MCalibrationCam&)(*fCam),*fBadPixels, + MBadPixelsPix::kRelTimeNotFitted, + MBadPixelsPix::kRelTimeOscillating); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Sets all pixels to MBadPixelsPix::kUnreliableRun, if following flags are set: +// - MBadPixelsPix::kRelTimeNotFitted +// - MBadPixelsPix::kRelTimeOscillating +// +void MHCalibrationRelTimeCam::FinalizeBadPixels() +{ + + for (Int_t i=0; iGetSize(); i++) { - MHCalibrationRelTimePix &hist = (*this)[i]; + MBadPixelsPix &bad = (*fBadPixels)[i]; + + if (bad.IsUncalibrated( MBadPixelsPix::kRelTimeNotFitted )) + bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); + + if (bad.IsUncalibrated( MBadPixelsPix::kRelTimeOscillating)) + bad.SetUnsuitable( MBadPixelsPix::kUnreliableRun ); - // - // 1) Return if the charge distribution is already succesfully fitted - // or if the histogram is empty - // - if (hist.IsGausFitOK() || hist.IsEmpty()) - continue; - - // - // 2) Fit the Hi Gain histograms with a Gaussian - // - hist.FitGaus(); - - // - // 3) If fit does not succeed , bypass the fit and take the histogram means and sigmas - // - if (!hist.IsGausFitOK()) - hist.BypassFit(); - - // - // 4) Create the fourier transform of the arrays - // - hist.CreateFourierSpectrum(); - - // - // 5) Renormalize to the real time in ns. - // - hist.Renorm(fTimeSliceWidth); - - } - return kTRUE; + } } @@ -250,5 +414,5 @@ // ============================================= // -// 4: Returned probability of Gauss fit to Charge distribution +// 4: Returned probability of Gauss fit to RelTime distribution // // Localized defects: @@ -261,30 +425,32 @@ { - if (fArray->GetSize() <= idx) + if (fHiGainArray->GetSize() <= idx) return kFALSE; + MHCalibrationRelTimePix &hist = (MHCalibrationRelTimePix&)(*this)[idx]; + switch (type) { case 0: - val = (*this)[idx].GetMean(); + val = hist.GetMean(); break; case 1: - val = (*this)[idx].GetMeanErr(); + val = hist.GetMeanErr(); break; case 2: - val = (*this)[idx].GetSigma(); + val = hist.GetSigma(); break; case 3: - val = (*this)[idx].GetSigmaErr(); + val = hist.GetSigmaErr(); break; case 4: - val = (*this)[idx].GetProb(); + val = hist.GetProb(); break; case 5: - if (!(*this)[idx].IsGausFitOK()) + if (!hist.IsGausFitOK()) val = 1.; break; case 6: - if (!(*this)[idx].IsFourierSpectrumOK()) + if (!hist.IsFourierSpectrumOK()) val = 1.; break; @@ -297,4 +463,5 @@ void MHCalibrationRelTimeCam::DrawPixelContent(Int_t idx) const { - (*this)[idx].DrawClone(); -} + const MHCalibrationRelTimePix &hist = (MHCalibrationRelTimePix&)(*this)[idx]; + hist.DrawClone(); +} Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.h (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimeCam.h (revision 3636) @@ -2,49 +2,24 @@ #define MARS_MHCalibrationRelTimeCam -#ifndef ROOT_TObjArray -#include +#ifndef MARS_MHCalibrationCam +#include "MHCalibrationCam.h" #endif -#ifndef MARS_MH -#include "MH.h" -#endif -#ifndef MARS_MCamEvent -#include "MCamEvent.h" -#endif - -class MHCalibrationRelTimePix; -class MHCalibrationRelTimeCam : public MH, public MCamEvent +class MGeomCam; +class MHCalibrationRelTimeCam : public MHCalibrationCam { private: - static const Float_t fgTimeSliceWidth; // Default for fTimeSliceWidth - static const Int_t fgPulserFrequency; // Default for fPulserFrequency - - Float_t fTimeSliceWidth; // FADC slice time width - Int_t fPulserFrequency; // The pulser frequency + Bool_t ReInitHists(MParList *pList); + Bool_t FillHists(const MParContainer *par, const Stat_t w=1); + Bool_t FinalizeHists(); + void FinalizeBadPixels(); - TObjArray *fArray; //-> List of MHCalibrationRelTimePix's - public: MHCalibrationRelTimeCam(const char *name=NULL, const char *title=NULL); - ~MHCalibrationRelTimeCam(); + ~MHCalibrationRelTimeCam() {} - void Clear(Option_t *o=""); - - MHCalibrationRelTimePix &operator[](UInt_t i); - const MHCalibrationRelTimePix &operator[](UInt_t i) const; - - Bool_t SetupFill(const MParList *pList); - Bool_t Fill(const MParContainer *par, const Stat_t w=1); - Bool_t Finalize(); - - // Setters - void SetTimeSliceWidth( const Float_t width=fgTimeSliceWidth) { fTimeSliceWidth = width; } - void SetPulserFrequency( const Int_t f=fgPulserFrequency) { fPulserFrequency = f; } - - TObject *Clone(const char *) const; - Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; void DrawPixelContent(Int_t idx) const; Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimePix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimePix.cc (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimePix.cc (revision 3636) @@ -40,7 +40,7 @@ using namespace std; // -const Int_t MHCalibrationRelTimePix::fgChargeNbins = 900; -const Axis_t MHCalibrationRelTimePix::fgChargeFirst = -13.; -const Axis_t MHCalibrationRelTimePix::fgChargeLast = 13.; +const Int_t MHCalibrationRelTimePix::fgRelTimeNbins = 900; +const Axis_t MHCalibrationRelTimePix::fgRelTimeFirst = -13.; +const Axis_t MHCalibrationRelTimePix::fgRelTimeLast = 13.; // -------------------------------------------------------------------------- // @@ -48,7 +48,7 @@ // // Sets: -// - the default number for fChargeNbins (fgChargeNbins) -// - the default number for fChargeFirst (fgChargeFirst) -// - the default number for fChargeLast (fgChargeLast) +// - the default number for fNbins (fgRelTimeNbins) +// - the default number for fFirst (fgRelTimeFirst) +// - the default number for fLast (fgRelTimeLast) // // - the default name of the fHGausHist ("HCalibrationRelTime") @@ -64,5 +64,4 @@ // MHCalibrationRelTimePix::MHCalibrationRelTimePix(const char *name, const char *title) - : fPixId(-1) { @@ -70,7 +69,7 @@ fTitle = title ? title : "Histogrammed Calibration Relative Arrival Time events"; - SetChargeNbins(); - SetChargeFirst(); - SetChargeLast(); + SetNbins ( fgRelTimeNbins ); + SetFirst ( fgRelTimeFirst ); + SetLast ( fgRelTimeLast ); // Create a large number of bins, later we will rebin @@ -96,4 +95,5 @@ fPixId = -1; + MHGausEvents::Clear(); return; @@ -109,47 +109,3 @@ } -// -------------------------------------------------------------------------- -// -// Sets: -// - fHGausHist.SetBins(fChargeNbins,fChargeFirst,fChargeLast); -// -void MHCalibrationRelTimePix::Init() -{ - fHGausHist.SetBins(fChargeNbins,fChargeFirst,fChargeLast); - -} - -// -------------------------------------------------------------------------- -// -// - Set fPixId to id -// -// Add id to names and titles of: -// - fHGausHist -// -void MHCalibrationRelTimePix::ChangeHistId(Int_t id) -{ - - fPixId = id; - - fHGausHist.SetName( Form("%s%d", fHGausHist.GetName(), id)); - fHGausHist.SetTitle( Form("%s%d", fHGausHist.GetTitle(), id)); - -} - - -// ---------------------------------------------------------------------- -// -// Renorm the results from FADC slices to times. -// The parameters slicewidth is the inverse of the FADC frequency and has the unit ns. -// -void MHCalibrationRelTimePix::Renorm(const Float_t slicewidth) -{ - - SetMean( GetMean() * slicewidth ); - SetMeanErr( GetMeanErr() * slicewidth ); - SetSigma( GetSigma() * slicewidth ); - SetSigmaErr( GetSigmaErr()* slicewidth ); - -} - Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimePix.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimePix.h (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationRelTimePix.h (revision 3636) @@ -11,13 +11,7 @@ private: - static const Int_t fgChargeNbins; // Default for fChargeNBins (now set to: 900 ) - static const Axis_t fgChargeFirst; // Default for fChargeFirst (now set to: -13.5 ) - static const Axis_t fgChargeLast; // Default for fChargeLast (now set to: 13.5 ) - - Int_t fChargeNbins; // Number of bins used for the fHGausHist - Axis_t fChargeFirst; // Lower bound bin used for the fHGausHist - Axis_t fChargeLast; // Upper bound bin used for the fHGausHist - - Int_t fPixId; // The pixel ID + static const Int_t fgRelTimeNbins; //! Default for MHGausEvents::fNBins (now set to: 900 ) + static const Axis_t fgRelTimeFirst; //! Default for MHGausEvents::fFirst (now set to: -13.5 ) + static const Axis_t fgRelTimeLast; //! Default for MHGausEvents::fLast (now set to: 13.5 ) public: @@ -26,17 +20,6 @@ ~MHCalibrationRelTimePix() {} - void Clear(Option_t *o=""); void Reset(); - void Init(); - - // Setters - void SetChargeNbins(const Int_t bins =fgChargeNbins) { fChargeNbins = bins; } - void SetChargeFirst(const Axis_t first=fgChargeFirst) { fChargeFirst = first; } - void SetChargeLast( const Axis_t last =fgChargeLast) { fChargeLast = last; } - - // Others - void ChangeHistId(Int_t i); - void Renorm(const Float_t slicewidth); ClassDef(MHCalibrationRelTimePix, 1) // Histogram class for Relative Time Calibration pixel Index: trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc (revision 3636) @@ -40,5 +40,6 @@ // It is created with a default name and title and resides in directory NULL. // - Any class deriving from MHGausEvents needs to apply a binning to fHGausHist -// (e.g. with the function TH1F::SetBins(..) ) +// (e.g. by setting the variables fNbins, fFirst, fLast and calling the function +// InitBins() or by directly calling fHGausHist.SetBins(..) ) // - The histogram is filled with the functions FillHist() or FillHistAndArray(). // - The histogram can be fitted with the function FitGaus(). This involves stripping @@ -48,4 +49,8 @@ // The NDF is compared to fNDFLimit and a check is made whether results are NaNs. // Accordingly, a flag IsGausFitOK() is set. +// - One can repeat the fit within a given amount of sigmas from the previous mean +// (specified by the variable fPickupLimit) with the function RepeatFit() +// - One can completely skip the fitting to set mean, sigma and its errors directly +// from the histograms with the function BypassFit() // // 2) The TArrayF fEvents: @@ -96,4 +101,5 @@ const Float_t MHGausEvents::fgProbLimit = 0.001; const Int_t MHGausEvents::fgNDFLimit = 2; +const Float_t MHGausEvents::fgPickupLimit = 5.; const Int_t MHGausEvents::fgPowerProbabilityBins = 20; const Int_t MHGausEvents::fgBinsAfterStripping = 40; @@ -105,8 +111,12 @@ // - the default Probability Limit for fProbLimit (fgProbLimit) // - the default NDF Limit for fNDFLimit (fgNDFLimit) +// - the default number for fPickupLimit (fgPickupLimit) // - the default number of bins after stripping for fBinsAfterStipping (fgBinsAfterStipping) // - the default name of the fHGausHist ("HGausHist") // - the default title of the fHGausHist ("Histogram of Events with Gaussian Distribution") // - the default directory of the fHGausHist (NULL) +// - the default for fNbins (100) +// - the default for fFirst (0.) +// - the default for fLast (100.) // // Initializes: @@ -121,4 +131,5 @@ fPowerSpectrum(NULL), fGraphEvents(NULL), fGraphPowerSpectrum(NULL), + fNbins(100), fFirst(0.), fLast(100.), fPixId(-1), fHGausHist(), fEvents(0), fFGausFit(NULL), fFExpFit(NULL) @@ -133,4 +144,5 @@ SetProbLimit(); SetNDFLimit(); + SetPickupLimit(); SetBinsAfterStripping(); @@ -183,7 +195,9 @@ // -------------------------------------------------------------------------- // +// Default Clear(), can be overloaded. +// // Sets: // - all other pointers to NULL -// - all variables to 0., except fEventFrequency +// - all variables to 0., except fPixId to -1 and keep fEventFrequency // - all flags to kFALSE // @@ -197,4 +211,5 @@ SetExpFitOK ( kFALSE ); SetFourierSpectrumOK( kFALSE ); + SetExcluded ( kFALSE ); fMean = 0.; @@ -205,4 +220,5 @@ fCurrentSize = 0; + fPixId = -1; if (fHPowerProbability) @@ -247,4 +263,6 @@ // -------------------------------------------------------------------------- +// +// Default Reset(), can be overloaded. // // Executes: @@ -264,4 +282,16 @@ // -------------------------------------------------------------------------- // +// Default InitBins, can be overloaded. +// +// Executes: +// - fHGausHist.SetBins(fNbins,fFirst,fLast) +// +void MHGausEvents::InitBins() +{ + fHGausHist.SetBins(fNbins,fFirst,fLast); +} + +// -------------------------------------------------------------------------- +// // Executes: // - FillArray() @@ -303,4 +333,13 @@ } +// -------------------------------------------------------------------------- +// +// Set Excluded bit from outside +// +void MHGausEvents::SetExcluded(const Bool_t b) +{ + b ? SETBIT(fFlags,kExcluded) : CLRBIT(fFlags,kExcluded); +} + const Double_t MHGausEvents::GetChiSquare() const @@ -366,4 +405,9 @@ { return TESTBIT(fFlags,kExpFitOK); +} + +const Bool_t MHGausEvents::IsExcluded() const +{ + return TESTBIT(fFlags,kExcluded); } @@ -409,4 +453,12 @@ if (fFExpFit) return; + + if (fEvents.GetSize() < 8) + { + *fLog << warn << "Cannot create Fourier spectrum in pixel: " << fPixId + << ". Number of events smaller than 8 " << endl; + return; + } + // @@ -424,4 +476,12 @@ // This cuts only the non-used zero's, but MFFT will later cut the rest MArray::StripZeros(fEvents); + + if (fEvents.GetSize() < 8) + { + *fLog << warn << "Cannot create Fourier spectrum. " << endl; + *fLog << warn << "Number of events (after stripping of zeros) is smaller than 8 " + << "in pixel: " << fPixId << endl; + return; + } MFFT fourier; @@ -506,5 +566,6 @@ if (!fFGausFit) { - *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit" << endl; + *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit " + << "in pixel: " << fPixId << endl; return kFALSE; } @@ -528,6 +589,6 @@ // The fit result is accepted under condition: // 1) The results are not nan's - // 2) The NDF is not smaller than fNDFLimit (5) - // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%) + // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit) + // 3) The Probability is greater than fProbLimit (default: fgProbLimit) // if ( TMath::IsNaN(fMean) @@ -545,4 +606,70 @@ // ----------------------------------------------------------------------------- +// +// If flag IsGausFitOK() is set (histogram already successfully fitted), +// returns kTRUE +// +// If both fMean and fSigma are still zero, call FitGaus() +// +// Repeats the Gauss fit in a smaller range, defined by: +// +// min = GetMean() - fPickupLimit * GetSigma(); +// max = GetMean() + fPickupLimit * GetSigma(); +// +// The fit results are retrieved and stored in class-own variables. +// +// A flag IsGausFitOK() is set according to whether the fit probability +// is smaller or bigger than fProbLimit, whether the NDF is bigger than +// fNDFLimit and whether results are NaNs. +// +Bool_t MHGausEvents::RepeatFit(const Option_t *option) +{ + + if (IsGausFitOK()) + return kTRUE; + + if ((fMean == 0.) && (fSigma == 0.)) + return FitGaus(); + + // + // Get new fitting ranges + // + Axis_t rmin = fMean - fPickupLimit * fSigma; + Axis_t rmax = fMean + fPickupLimit * fSigma; + + Axis_t hmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst()); + Axis_t hmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast()) ; + + fFGausFit->SetRange(hmin < rmin ? rmin : hmin , hmax > rmax ? rmax : hmax); + + fHGausHist.Fit(fFGausFit,option); + + fMean = fFGausFit->GetParameter(1); + fMeanErr = fFGausFit->GetParameter(2); + fSigma = fFGausFit->GetParError(1) ; + fSigmaErr = fFGausFit->GetParError(2) ; + fProb = fFGausFit->GetProb() ; + + // + // The fit result is accepted under condition: + // 1) The results are not nan's + // 2) The NDF is not smaller than fNDFLimit (default: fgNDFLimit) + // 3) The Probability is greater than fProbLimit (default: fgProbLimit) + // + if ( TMath::IsNaN ( fMean ) + || TMath::IsNaN ( fMeanErr ) + || TMath::IsNaN ( fProb ) + || TMath::IsNaN ( fSigma ) + || TMath::IsNaN ( fSigmaErr ) + || fFGausFit->GetNDF() < fNDFLimit + || fProb < fProbLimit ) + return kFALSE; + + SetGausFitOK(kTRUE); + return kTRUE; + +} + +// ----------------------------------------------------------------------------- // // Bypasses the Gauss fit by taking mean and RMS from the histogram @@ -559,12 +686,13 @@ if (entries <= 0.) { - *fLog << warn << GetDescriptor() << ": Cannot bypass fit. Number of entries smaller or equal 0" << endl; + *fLog << warn << GetDescriptor() + << ": Cannot bypass fit. Number of entries smaller or equal 0 in pixel: " << fPixId << endl; return; } - SetMean ( fHGausHist.GetMean() ); - SetMeanErr ( fHGausHist.GetRMS() / TMath::Sqrt(entries) ); - SetSigma ( fHGausHist.GetRMS() ); - SetSigmaErr ( fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2. ); + fMean = fHGausHist.GetMean(); + fMeanErr = fHGausHist.GetRMS() / TMath::Sqrt(entries); + fSigma = fHGausHist.GetRMS() ; + fSigmaErr = fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2.; } @@ -577,5 +705,5 @@ *fLog << all << endl; - *fLog << all << "Results of the Gauss Fit: " << endl; + *fLog << all << "Results of the Gauss Fit in pixel: " << fPixId << endl; *fLog << all << "Mean: " << GetMean() << endl; *fLog << all << "Sigma: " << GetSigma() << endl; @@ -585,4 +713,24 @@ *fLog << all << endl; +} + +// -------------------------------------------------------------------------- +// +// - Set fPixId to id +// +// Add id to names and titles of: +// - fHGausHist +// +void MHGausEvents::ChangeHistId(Int_t id) +{ + + fPixId = id; + + fHGausHist.SetName( Form("%s%d", fHGausHist.GetName(), id)); + fHGausHist.SetTitle( Form("%s%d", fHGausHist.GetTitle(), id)); + + fName = Form("%s%d", fName.Data(), id); + fTitle = Form("%s%d", fTitle.Data(), id); + } @@ -766,3 +914,18 @@ } - +const Double_t MHGausEvents::GetPickup() const +{ + + if ((fMean == 0.) && (fSigma == 0.)) + return -1.; + + const Int_t first = fHGausHist.GetXaxis()->FindBin(fMean+fPickupLimit*fSigma); + const Int_t last = fHGausHist.GetXaxis()->GetLast(); + + if (first >= last) + return 0.; + + return fHGausHist.Integral(first, last, "width"); +} + + Index: trunk/MagicSoft/Mars/mcalib/MHGausEvents.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHGausEvents.h (revision 3635) +++ trunk/MagicSoft/Mars/mcalib/MHGausEvents.h (revision 3636) @@ -20,9 +20,10 @@ private: - const static Float_t fgProbLimit; // Default for fProbLimit (now 0.001) - const static Int_t fgNDFLimit; // Default for fNDFLimit (now 2) - const static Int_t fgPowerProbabilityBins; // Default for fPowerProbabilityBins (now 20) - const static Int_t fgBinsAfterStripping; // Default for fBinsAfterStripping (now 40) - + const static Float_t fgProbLimit; //! Default for fProbLimit (now set to: 0.001) + const static Int_t fgNDFLimit; //! Default for fNDFLimit (now set to: 2) + const static Float_t fgPickupLimit; //! Default for fPickupLimit (now set to: 5. ) + const static Int_t fgPowerProbabilityBins; //! Default for fPowerProbabilityBins (now set to: 20) + const static Int_t fgBinsAfterStripping; //! Default for fBinsAfterStripping (now set to: 40) + Int_t fPowerProbabilityBins; // Bins for the projected power spectrum Int_t fBinsAfterStripping; // Bins for the Gauss Histogram after stripping off the zeros at both ends @@ -41,12 +42,18 @@ Double_t fProb; // Probability of the Gauss fit (derived from Chi-Square and NDF - enum { kGausFitOK, kExpFitOK, kFourierSpectrumOK }; // Bits to hold information about fit results + enum { kGausFitOK, kExpFitOK, kFourierSpectrumOK, + kExcluded }; // Bits to hold information about fit results - Byte_t fFlags; // Byte to hold the bits fit result bits + Byte_t fFlags; // Byte to hold the bits fit result bits - Int_t fCurrentSize; // Current size of the array fEvents - + Int_t fCurrentSize; // Current size of the array fEvents + protected: + Int_t fNbins; // Number histogram bins for fHGausHist (used by Init()) + Axis_t fFirst; // Lower histogram edge for fHGausHist (used by Init()) + Axis_t fLast; // Upper histogram edge for fHGausHist (used by Init()) + Int_t fPixId; // Pixel ID + TH1F fHGausHist; // Histogram which should hold the Gaussian distribution TArrayF fEvents; // Array which holds the entries of GausHist @@ -57,4 +64,5 @@ Float_t fProbLimit; // Probability limit for judgement if fit is OK Int_t fNDFLimit; // NDF limit for judgement if fit is OK + Float_t fPickupLimit; // Upper number of sigmas from the mean until events are considered as pickup Float_t *CreateEventXaxis(Int_t n); // Create an x-axis for the Event TGraphs @@ -75,5 +83,6 @@ virtual void Clear(Option_t *o=""); virtual void Reset(); - + virtual void InitBins(); + // Setters void SetEventFrequency(const Float_t f) { fEventFrequency = f; } @@ -84,12 +93,19 @@ void SetSigmaErr( const Double_t d ) { fSigmaErr = d; } void SetProb ( const Double_t d ) { fProb = d; } + void SetPixId ( const Int_t i ) { fPixId = i; } + + void SetNbins ( const Int_t i ) { fNbins = i; } + void SetFirst ( const Double_t d ) { fFirst = d; } + void SetLast ( const Double_t d ) { fLast = d; } + void SetNDFLimit( const Int_t lim=fgNDFLimit ) { fNDFLimit = lim; } + void SetPickupLimit( const Float_t lim =fgPickupLimit) { fPickupLimit = lim; } void SetProbLimit( const Float_t lim=fgProbLimit ) { fProbLimit = lim; } - void SetNDFLimit( const Int_t lim=fgNDFLimit ) { fNDFLimit = lim; } // Setters ONLY for MC: - void SetGausFitOK( const Bool_t b ); - void SetExpFitOK( const Bool_t b ); - void SetFourierSpectrumOK( const Bool_t b ); + void SetGausFitOK( const Bool_t b=kTRUE ); + void SetExpFitOK( const Bool_t b=kTRUE ); + void SetFourierSpectrumOK( const Bool_t b=kTRUE ); + void SetExcluded ( const Bool_t b=kTRUE ); // Getters @@ -101,4 +117,6 @@ const Double_t GetProb() const { return fProb; } const Int_t GetNdf() const; + const Double_t GetPickup() const; + const Int_t GetPixId() const { return fPixId; } const Double_t GetSlope() const; @@ -136,4 +154,5 @@ const Bool_t IsEmpty() const; const Bool_t IsFourierSpectrumOK() const; + const Bool_t IsExcluded () const; // Fill @@ -145,6 +164,6 @@ Bool_t FitGaus(Option_t *option="RQ0", const Double_t xmin=0., const Double_t xmax=0.); // Fit the histogram HGausHist with a Gaussian - void BypassFit(); // Take mean and RMS from the histogram - + Bool_t RepeatFit(const Option_t *option="RQ0"); // Repeat fit within limits defined by fPickupLimit + void BypassFit(); // Take mean and RMS from the histogram // Draws @@ -155,4 +174,5 @@ // Miscelleaneous + virtual void ChangeHistId(Int_t id); // Changes names and titles of the histogram void CreateFourierSpectrum(); // Create the fourier spectrum out of fEvents void CreateGraphEvents(); // Create the TGraph fGraphEvents of fEvents