Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 3153) +++ /trunk/MagicSoft/Mars/Changelog (revision 3154) @@ -20,4 +20,9 @@ - "CutEdges" renamed to "StripZeros" + * mcalib/MHGausEvent.[h,cc] + * mcalib/Makefile + * mcalib/CalibLinkDef.h + - replaced by the improved version: MHGausEvents.[h,cc] + 2004/02/14: Thomas Bretz Index: /trunk/MagicSoft/Mars/mcalib/CalibLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/CalibLinkDef.h (revision 3153) +++ /trunk/MagicSoft/Mars/mcalib/CalibLinkDef.h (revision 3154) @@ -20,5 +20,5 @@ #pragma link C++ class MMcCalibrationCalc++; -#pragma link C++ class MHGausEvent++; +#pragma link C++ class MHGausEvents++; #endif Index: unk/MagicSoft/Mars/mcalib/MHGausEvent.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHGausEvent.cc (revision 3153) +++ (revision ) @@ -1,317 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * This file is part of MARS, the MAGIC Analysis and Reconstruction -! * Software. It is distributed to you in the hope that it can be a useful -! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. -! * It is distributed WITHOUT ANY WARRANTY. -! * -! * Permission to use, copy, modify and distribute this software and its -! * documentation for any purpose is hereby granted without fee, -! * provided that the above copyright notice appear in all copies and -! * that both that copyright notice and this permission notice appear -! * in supporting documentation. It is provided "as is" without express -! * or implied warranty. -! * -! -! -! Author(s): Markus Gaug 11/2003 -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// -// MHGausEvent -// -// A base class for all kind of event which follow a Gaussian distribution -// with time, i.e. observables containing white noise. -// -// The class provides the basic tools for fitting, -// spectrum analysis, etc. -// -////////////////////////////////////////////////////////////////////////////// -#include "MHGausEvent.h" - -#include -#include - -#include "MFFT.h" -#include "MArray.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHGausEvent); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default Constructor. -// -MHGausEvent::MHGausEvent(const char *name, const char *title) - : fHGausHist(NULL), fHPowerProbability(NULL), - fFGausFit(NULL), fFExpFit(NULL), - fEvents(NULL), fPowerSpectrum(NULL) -{ - - fName = name ? name : "MHGausEvent"; - fTitle = title ? title : "Events which follow a Gaussian distribution"; - - Clear(); -} - - -MHGausEvent::~MHGausEvent() -{ - Clear(); - - if (fHGausHist) - delete fHGausHist; -} - - - -void MHGausEvent::Clear(Option_t *o) -{ - - fGausHistBins = 50; - fGausHistAxisFirst = 0.; - fGausHistAxisLast = 100.; - - fPowerProbabilityBins = 30; - - fProbLimit = 0.01; - fGausFitOK = kFALSE; - fExpFitOK = kFALSE; - fOscillating = kFALSE; - - fMean = 0.; - fSigma = 0.; - fMeanErr = 0.; - fSigmaErr = 0.; - - fProb = 0.; - - if (fHPowerProbability) - delete fHPowerProbability; - if (fFGausFit) - delete fFGausFit; - if (fEvents) - delete fEvents; - if (fPowerSpectrum) - delete fPowerSpectrum; -} - - -void MHGausEvent::Reset() -{ - - Clear(); - fHGausHist->Reset(); - -} - -const Double_t MHGausEvent::GetChiSquare() const -{ - return ( fFGausFit ? fFGausFit->GetChisquare() : 0.); -} - -const Int_t MHGausEvent::GetNdf() const -{ - return ( fFGausFit ? fFGausFit->GetNDF() : 0); -} - - -const Double_t MHGausEvent::GetExpChiSquare() const -{ - return ( fFExpFit ? fFExpFit->GetChisquare() : 0.); -} - - -const Int_t MHGausEvent::GetExpNdf() const -{ - return ( fFExpFit ? fFExpFit->GetNDF() : 0); -} - -const Double_t MHGausEvent::GetExpProb() const -{ - return ( fFExpFit ? fFExpFit->GetProb() : 0.); -} - -const Double_t MHGausEvent::GetOffset() const -{ - return ( fFExpFit ? fFExpFit->GetParameter(0) : 0.); -} - -const Double_t MHGausEvent::GetSlope() const -{ - return ( fFExpFit ? fFExpFit->GetParameter(1) : 0.); -} - - - -Bool_t MHGausEvent::CheckOscillations() -{ - - if (fFExpFit) - return IsOscillating(); - - if (!fEvents) - return kFALSE; - - // - // The number of entries HAS to be a potence of 2, - // so we can only cut out from the last potence of 2 to the rest. - // Another possibility would be to fill everything with - // zeros, but that gives a low frequency peak, which we would - // have to cut out later again. - // - // So, we have to live with the possibility that at the end - // of the calibration run, something has happened without noticing - // it... - // - - // This cuts only the non-used zero's, but MFFT will later cut the rest - MArray::StripZeros(*fEvents); - - MFFT fourier; - - fPowerSpectrum = fourier.PowerSpectrumDensity(fEvents); - fHPowerProbability = ProjectArray(*fPowerSpectrum, fPowerProbabilityBins, - "PowerProbability", - "Probability of Power occurrance"); - // - // First guesses for the fit (should be as close to reality as possible, - // - const Double_t xmax = fHPowerProbability->GetXaxis()->GetXmax(); - - fFExpFit = new TF1("FExpFit","exp([0]-[1]*x)",0.,xmax); - - const Double_t slope_guess = (TMath::Log(fHPowerProbability->GetEntries())+1.)/xmax; - const Double_t offset_guess = slope_guess*xmax; - - fFExpFit->SetParameters(offset_guess, slope_guess); - fFExpFit->SetParNames("Offset","Slope"); - fFExpFit->SetParLimits(0,offset_guess/2.,2.*offset_guess); - fFExpFit->SetParLimits(1,slope_guess/1.5,1.5*slope_guess); - fFExpFit->SetRange(0.,xmax); - - fHPowerProbability->Fit(fFExpFit,"RQL0"); - - if (GetExpProb() < fProbLimit) - fExpFitOK = kFALSE; - - // For the moment, this is the only check, later we can add more... - fOscillating = fExpFitOK; - - return fOscillating; -} - - -Bool_t MHGausEvent::IsEmpty() const -{ - return !(fHGausHist->GetEntries()); -} - -Bool_t MHGausEvent::IsOscillating() -{ - - if (fFExpFit) - return fOscillating; - - return CheckOscillations(); - -} - - -Bool_t MHGausEvent::FitGaus(Option_t *option) -{ - - if (IsGausFitOK()) - return kTRUE; - - // - // First, cut the edges which contain only zeros and rebin - // to about 20 bins. - // - // (ATTENTION: The Chisquare method is more sensitive, - // the _less_ bins, you have!) - // - Int_t newbins = 20; - MH::StripZeros(fHGausHist,newbins); - - // - // Get the fitting ranges - // - Axis_t rmin = fHGausHist->GetXaxis()->GetFirst(); - Axis_t rmax = fHGausHist->GetXaxis()->GetLast(); - - // - // First guesses for the fit (should be as close to reality as possible, - // - const Stat_t entries = fHGausHist->Integral("width"); - const Double_t mu_guess = fHGausHist->GetBinCenter(fHGausHist->GetMaximumBin()); - const Double_t sigma_guess = (rmax-rmin)/2.; - const Double_t area_guess = entries/TMath::Sqrt(TMath::TwoPi())/sigma_guess; - - fFGausFit = new TF1("GausFit","gaus",rmin,rmax); - - if (!fFGausFit) - { - *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit" << endl; - return kFALSE; - } - - fFGausFit->SetParameters(area_guess,mu_guess,sigma_guess); - fFGausFit->SetParNames("Area","#mu","#sigma"); - fFGausFit->SetParLimits(0,0.,entries); - fFGausFit->SetParLimits(1,rmin,rmax); - fFGausFit->SetParLimits(2,0.,rmax-rmin); - fFGausFit->SetRange(rmin,rmax); - - fHGausHist->Fit(fFGausFit,option); - - fMean = fFGausFit->GetParameter(1); - fSigma = fFGausFit->GetParameter(2); - fMeanErr = 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 (5) - // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%) - // - if ( TMath::IsNaN(fMean) - || TMath::IsNaN(fMeanErr) - || TMath::IsNaN(fProb) - || TMath::IsNaN(fSigma) - || TMath::IsNaN(fSigmaErr) ) - { - fGausFitOK = kFALSE; - return kFALSE; - } - - fGausFitOK = kTRUE; - return kTRUE; -} - -void MHGausEvent::Print(const Option_t *o) const -{ - - *fLog << all << endl; - *fLog << all << "Results of the Gauss Fit: " << endl; - *fLog << all << "Mean: " << GetMean() << endl; - *fLog << all << "Sigma: " << GetSigma() << endl; - *fLog << all << "Chisquare: " << GetChiSquare() << endl; - *fLog << all << "DoF: " << GetNdf() << endl; - *fLog << all << "Probability: " << GetProb() << endl; - *fLog << all << endl; - -} - Index: unk/MagicSoft/Mars/mcalib/MHGausEvent.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHGausEvent.h (revision 3153) +++ (revision ) @@ -1,97 +1,0 @@ -#ifndef MARS_MHGausEvent -#define MARS_MHGausEvent - -#ifndef MARS_MH -#include "MH.h" -#endif - -class TArrayF; -class TH1F; -class TH1I; -class TF1; -class MHGausEvent : public MH -{ -private: - - Bool_t fGausFitOK; - Bool_t fExpFitOK; - Bool_t fOscillating; - - Double_t fMean; - Double_t fSigma; - Double_t fMeanErr; - Double_t fSigmaErr; - - Double_t fProb; - -protected: - - TH1F* fHGausHist; - TH1I* fHPowerProbability; - - TF1* fFGausFit; - TF1* fFExpFit; - - TArrayF *fEvents; - TArrayF *fPowerSpectrum; - - Int_t fGausHistBins; - Axis_t fGausHistAxisFirst; - Axis_t fGausHistAxisLast; - - Int_t fPowerProbabilityBins; - Float_t fProbLimit; - - void CutArrayBorder(TArrayF *array); - -public: - - MHGausEvent(const char *name=NULL, const char *title=NULL); - ~MHGausEvent(); - - void Clear(Option_t *o=""); - void Reset(); - - const Double_t GetMean() const { return fMean; } - const Double_t GetMeanErr() const { return fMeanErr; } - const Double_t GetSigma() const { return fSigma; } - const Double_t GetSigmaErr() const { return fSigmaErr; } - const Double_t GetChiSquare() const; - const Double_t GetProb() const { return fProb; } - const Int_t GetNdf() const; - - const Double_t GetSlope() const; - const Double_t GetOffset() const; - const Double_t GetExpChiSquare() const; - const Double_t GetExpProb() const; - const Int_t GetExpNdf() const; - - const TH1F *GetHGausHist() { return fHGausHist; } - const TH1F *GetHGausHist() const { return fHGausHist; } - - const TF1 *GetFGausFit() { return fFGausFit; } - const TF1 *GetFGausFit() const { return fFGausFit; } - - const TH1I *GetHPowerProbability() { return fHPowerProbability; } - const TH1I *GetHPowerProbability() const { return fHPowerProbability; } - - const TF1 *GetFExpFit() { return fFExpFit; } - const TF1 *GetFExpFit() const { return fFExpFit; } - - Bool_t CheckOscillations(); - - virtual Bool_t IsGausFitOK() const { return fGausFitOK; } - virtual Bool_t IsExpFitOK() const { return fExpFitOK; } - virtual Bool_t IsEmpty() const; - virtual Bool_t IsOscillating(); - - // Fits - virtual Bool_t FitGaus(Option_t *option="RQ0"); - - // Prints - virtual void Print(const Option_t *o="") const; - - ClassDef(MHGausEvent, 1) // Histograms for events with Gaussian distributed values -}; - -#endif Index: /trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc (revision 3154) +++ /trunk/MagicSoft/Mars/mcalib/MHGausEvents.cc (revision 3154) @@ -0,0 +1,639 @@ +/* ======================================================================== *\ +! +! * +! * This file is part of MARS, the MAGIC Analysis and Reconstruction +! * Software. It is distributed to you in the hope that it can be a useful +! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. +! * It is distributed WITHOUT ANY WARRANTY. +! * +! * Permission to use, copy, modify and distribute this software and its +! * documentation for any purpose is hereby granted without fee, +! * provided that the above copyright notice appear in all copies and +! * that both that copyright notice and this permission notice appear +! * in supporting documentation. It is provided "as is" without express +! * or implied warranty. +! * +! +! +! Author(s): Markus Gaug 02/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHGausEvents +// +// A base class for events which are believed follow a Gaussian distribution +// with time, (e.g. calibration events, observables containing white noise, etc.) +// +// MHGausEvents derives from MH, thus all features of MH can be used by a class +// deriving from MHGausEvents, especially the filling functions. +// +// The central objects are: +// +// 1) The TH1F fHGausHist: +// ============================== +// +// 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(..) ) +// - The histogram is filled with the functions FillHist or FillHistAndArray. +// - The histogram can be fitted with the function FitGaus(). This involves stripping +// of all zeros at the lower and upper end of the histogram and re-binning to +// a new number of bins, specified in the variable fBinsAfterStripping. +// - The fit result's probability is compared to a reference probability fProbLimit +// The NDF is compared to fNDFLimit and a check is made whether results are NaNs. +// Accordingly the flag GausFitOK is set. +// +// 2) The TArrayF fEvents: +// ========================== +// +// It is created with 0 entries and not expanded unless FillArray or FillHistAndArray is called. +// - A first call to FillArray or FillHistAndArray initializes fEvents by default to 512 entries. +// - Any later call to FillArray or FillHistAndArray fills up the array. Reaching the limit, +// the array is expanded by a factor 2. +// - The array can be fourier-transformed into the array fPowerSpectrum. Note that any FFT accepts +// only number of events which are a power of 2. Thus, fEvents is cut to the next power of 2 smaller +// than its actual number of entries. You might lose information at the end of your analysis. +// - Calling the function CreateFourierTransform creates the array fPowerSpectrum and its projection +// fHPowerProbability which in turn is fit to an exponential. +// - The fit result's probability is compared to a referenc probability fProbLimit +// and accordingly the flag ExpFitOK is set. +// - The flag FourierSpectrumOK is set accordingly to ExpFitOK. Later, a closer check will be +// installed. +// - You can display all arrays by calls to: CreateGraphEvents() and CreateGraphPowerSpectrum() and +// successive calls to the corresponding Getters. +// +////////////////////////////////////////////////////////////////////////////// +#include "MHGausEvents.h" + +#include +#include +#include +#include +#include +#include +#include + +#include "MFFT.h" +#include "MArray.h" + +#include "MH.h" + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHGausEvents); + +using namespace std; + +const Float_t MHGausEvents::fgProbLimit = 0.01; +const Int_t MHGausEvents::fgNDFLimit = 2; +const Int_t MHGausEvents::fgPowerProbabilityBins = 25; +const Int_t MHGausEvents::fgBinsAfterStripping = 25; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MHGausEvents::MHGausEvents(const char *name, const char *title) + : fFGausFit(NULL), fFExpFit(NULL), + fHPowerProbability(NULL), + fPowerSpectrum(NULL), + fGraphEvents(NULL), fGraphPowerSpectrum(NULL), + fHGausHist(), + fEvents(0) +{ + + fName = name ? name : "MHGausEvents"; + fTitle = title ? title : "Events with expected Gaussian distributions"; + + Clear(); + + SetPowerProbabilityBins(); + SetEventFrequency(); + SetProbLimit(); + SetNDFLimit(); + SetBinsAfterStripping(); + + fHGausHist.SetName("HGausHist"); + fHGausHist.SetTitle("Histogram of Events with Gaussian Distribution"); + // important, other ROOT will not draw the axes: + fHGausHist.UseCurrentStyle(); + fHGausHist.SetDirectory(NULL); +} + + + + +MHGausEvents::~MHGausEvents() +{ + + // delete histograms + if (fHPowerProbability) + delete fHPowerProbability; + + // delete fits + if (fFGausFit) + delete fFGausFit; + if (fFExpFit) + delete fFExpFit; + + // delete arrays + if (fPowerSpectrum) + delete fPowerSpectrum; + + // delete graphs + if (fGraphEvents) + delete fGraphEvents; + if (fGraphPowerSpectrum) + delete fGraphPowerSpectrum; + +} + + +void MHGausEvents::Clear(Option_t *o) +{ + + SetGausFitOK ( kFALSE ); + SetExpFitOK ( kFALSE ); + SetFourierSpectrumOK( kFALSE ); + + fMean = 0.; + fSigma = 0.; + fMeanErr = 0.; + fSigmaErr = 0.; + fProb = 0.; + + fCurrentSize = 0; + + if (fHPowerProbability) + delete fHPowerProbability; + + // delete fits + if (fFGausFit) + delete fFGausFit; + if (fFExpFit) + delete fFExpFit; + + // delete arrays + if (fPowerSpectrum) + delete fPowerSpectrum; + + // delete graphs + if (fGraphEvents) + delete fGraphEvents; + if (fGraphPowerSpectrum) + delete fGraphPowerSpectrum; +} + + +void MHGausEvents::Reset() +{ + + Clear(); + fHGausHist.Reset(); + fEvents.Set(0); + +} + + +Bool_t MHGausEvents::FillHistAndArray(Float_t f) +{ + + FillArray(f); + return FillHist(f); +} + +Bool_t MHGausEvents::FillHist(Float_t f) +{ + + if (fHGausHist.Fill(f) == -1) + return kFALSE; + + return kTRUE; +} + +void MHGausEvents::FillArray(Float_t f) +{ + if (fEvents.GetSize() == 0) + fEvents.Set(512); + + if (fCurrentSize >= fEvents.GetSize()) + fEvents.Set(fEvents.GetSize()*2); + + fEvents.AddAt(f,fCurrentSize++); +} + + +const Double_t MHGausEvents::GetChiSquare() const +{ + return ( fFGausFit ? fFGausFit->GetChisquare() : 0.); +} + +const Int_t MHGausEvents::GetNdf() const +{ + return ( fFGausFit ? fFGausFit->GetNDF() : 0); +} + + +const Double_t MHGausEvents::GetExpChiSquare() const +{ + return ( fFExpFit ? fFExpFit->GetChisquare() : 0.); +} + + +const Int_t MHGausEvents::GetExpNdf() const +{ + return ( fFExpFit ? fFExpFit->GetNDF() : 0); +} + +const Double_t MHGausEvents::GetExpProb() const +{ + return ( fFExpFit ? fFExpFit->GetProb() : 0.); +} + +const Double_t MHGausEvents::GetOffset() const +{ + return ( fFExpFit ? fFExpFit->GetParameter(0) : 0.); +} + +const Double_t MHGausEvents::GetSlope() const +{ + return ( fFExpFit ? fFExpFit->GetParameter(1) : 0.); +} + +const Bool_t MHGausEvents::IsEmpty() const +{ + return !(fHGausHist.GetEntries()); +} + +const Bool_t MHGausEvents::IsFourierSpectrumOK() const +{ + return TESTBIT(fFlags,kFourierSpectrumOK); +} + +const Bool_t MHGausEvents::IsGausFitOK() const +{ + return TESTBIT(fFlags,kGausFitOK); +} + +const Bool_t MHGausEvents::IsExpFitOK() const +{ + return TESTBIT(fFlags,kExpFitOK); +} + +// ------------------------------------------------------------------- +// +// The flag setters are to be used ONLY for Monte-Carlo!! +// +void MHGausEvents::SetGausFitOK(const Bool_t b) +{ + b ? SETBIT(fFlags,kGausFitOK) : CLRBIT(fFlags,kGausFitOK); + +} +// ------------------------------------------------------------------- +// +// The flag setters are to be used ONLY for Monte-Carlo!! +// +void MHGausEvents::SetExpFitOK(const Bool_t b) +{ + + b ? SETBIT(fFlags,kExpFitOK) : CLRBIT(fFlags,kExpFitOK); +} + +// ------------------------------------------------------------------- +// +// The flag setters are to be used ONLY for Monte-Carlo!! +// +void MHGausEvents::SetFourierSpectrumOK(const Bool_t b) +{ + + b ? SETBIT(fFlags,kFourierSpectrumOK) : CLRBIT(fFlags,kFourierSpectrumOK); +} + + +// ------------------------------------------------------------------- +// +// Create the fourier spectrum using the class MFFT. +// The result is projected into a histogram and fitted by an exponential +// +void MHGausEvents::CreateFourierSpectrum() +{ + + if (fFExpFit) + return; + + // + // The number of entries HAS to be a potence of 2, + // so we can only cut out from the last potence of 2 to the rest. + // Another possibility would be to fill everything with + // zeros, but that gives a low frequency peak, which we would + // have to cut out later again. + // + // So, we have to live with the possibility that at the end + // of the calibration run, something has happened without noticing + // it... + // + + // This cuts only the non-used zero's, but MFFT will later cut the rest + MArray::StripZeros(fEvents); + + MFFT fourier; + + fPowerSpectrum = fourier.PowerSpectrumDensity(&fEvents); + fHPowerProbability = ProjectArray(*fPowerSpectrum, fPowerProbabilityBins, + "PowerProbability", + "Probability of Power occurrance"); + fHPowerProbability->SetXTitle("P(f)"); + fHPowerProbability->SetDirectory(NULL); + // + // First guesses for the fit (should be as close to reality as possible, + // + const Double_t xmax = fHPowerProbability->GetXaxis()->GetXmax(); + + fFExpFit = new TF1("FExpFit","exp([0]-[1]*x)",0.,xmax); + + const Double_t slope_guess = (TMath::Log(fHPowerProbability->GetEntries())+1.)/xmax; + const Double_t offset_guess = slope_guess*xmax; + + fFExpFit->SetParameters(offset_guess, slope_guess); + fFExpFit->SetParNames("Offset","Slope"); + fFExpFit->SetParLimits(0,offset_guess/2.,2.*offset_guess); + fFExpFit->SetParLimits(1,slope_guess/1.5,1.5*slope_guess); + fFExpFit->SetRange(0.,xmax); + + fHPowerProbability->Fit(fFExpFit,"RQL0"); + + if (GetExpProb() < fProbLimit) + SetExpFitOK(kFALSE); + + // + // For the moment, this is the only check, later we can add more... + // + SetFourierSpectrumOK(IsExpFitOK()); + + return; +} + +// ------------------------------------------------------------------- +// +// Fit fGausHist with a Gaussian after stripping zeros from both ends +// and rebinned to the number of bins specified in fBinsAfterStripping +// +// 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::FitGaus(Option_t *option) +{ + + if (IsGausFitOK()) + return kTRUE; + + // + // First, strip the zeros from the edges which contain only zeros and rebin + // to about fBinsAfterStripping bins. + // + // (ATTENTION: The Chisquare method is more sensitive, + // the _less_ bins, you have!) + // + StripZeros(&fHGausHist,fBinsAfterStripping); + + // + // Get the fitting ranges + // + Axis_t rmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst()); + Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast()); + + // + // First guesses for the fit (should be as close to reality as possible, + // + const Stat_t entries = fHGausHist.Integral("width"); + const Double_t mu_guess = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin()); + const Double_t sigma_guess = (rmax-rmin)/2.; + const Double_t area_guess = entries/TMath::Sqrt(TMath::TwoPi())/sigma_guess; + + fFGausFit = new TF1("GausFit","gaus",rmin,rmax); + + if (!fFGausFit) + { + *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit" << endl; + return kFALSE; + } + + fFGausFit->SetParameters(area_guess,mu_guess,sigma_guess); + fFGausFit->SetParNames("Area","#mu","#sigma"); + fFGausFit->SetParLimits(0,0.,entries); + fFGausFit->SetParLimits(1,rmin,rmax); + fFGausFit->SetParLimits(2,0.,rmax-rmin); + fFGausFit->SetRange(rmin,rmax); + + fHGausHist.Fit(fFGausFit,option); + + fMean = fFGausFit->GetParameter(1); + fSigma = fFGausFit->GetParameter(2); + fMeanErr = 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 (5) + // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%) + // + 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; +} + +// ----------------------------------------------------------------------------------- +// +// A default print +// +void MHGausEvents::Print(const Option_t *o) const +{ + + *fLog << all << endl; + *fLog << all << "Results of the Gauss Fit: " << endl; + *fLog << all << "Mean: " << GetMean() << endl; + *fLog << all << "Sigma: " << GetSigma() << endl; + *fLog << all << "Chisquare: " << GetChiSquare() << endl; + *fLog << all << "DoF: " << GetNdf() << endl; + *fLog << all << "Probability: " << GetProb() << endl; + *fLog << all << endl; + +} + +// ---------------------------------------------------------------------------------- +// +// Create a graph to display the array fEvents +// If the variable fEventFrequency is set, the x-axis is transformed into real time. +// +void MHGausEvents::CreateGraphEvents() +{ + + MArray::StripZeros(fEvents); + + const Int_t n = fEvents.GetSize(); + + fGraphEvents = new TGraph(n,CreateXaxis(n),fEvents.GetArray()); + fGraphEvents->SetTitle("Evolution of Events with time"); + fGraphEvents->GetXaxis()->SetTitle((fEventFrequency) ? "Time [s]" : "Event Nr."); +} + +// ---------------------------------------------------------------------------------- +// +// Create a graph to display the array fPowerSpectrum +// If the variable fEventFrequency is set, the x-axis is transformed into real frequency. +// +void MHGausEvents::CreateGraphPowerSpectrum() +{ + + MArray::StripZeros(*fPowerSpectrum); + + const Int_t n = fPowerSpectrum->GetSize(); + + fGraphPowerSpectrum = new TGraph(n,CreateXaxis(n),fPowerSpectrum->GetArray()); + fGraphPowerSpectrum->SetTitle("Power Spectrum Density"); + fGraphPowerSpectrum->GetXaxis()->SetTitle((fEventFrequency) ? "Frequency [Hz]" : "Frequency"); + fGraphPowerSpectrum->GetYaxis()->SetTitle("P(f)"); +} + +// ----------------------------------------------------------------------------- +// +// Create the x-axis for the graph +// +Float_t *MHGausEvents::CreateXaxis(Int_t n) +{ + + Float_t *xaxis = new Float_t[n]; + + if (fEventFrequency) + for (Int_t i=0;iDivide(1,win); + pad->cd(1); + + if (!IsEmpty()) + pad->SetLogy(); + + fHGausHist.Draw(opt); + + if (fFGausFit) + { + fFGausFit->SetLineColor(IsGausFitOK() ? kGreen : kRed); + fFGausFit->Draw("same"); + } + + switch (win) + { + case 2: + pad->cd(2); + DrawEvents(); + break; + case 3: + pad->cd(2); + DrawPowerSpectrum(*pad,3); + break; + case 4: + pad->cd(2); + DrawEvents(); + pad->cd(3); + DrawPowerSpectrum(*pad,4); + break; + } +} + +void MHGausEvents::DrawEvents() +{ + + if (!fGraphEvents) + CreateGraphEvents(); + + fGraphEvents->SetBit(kCanDelete); + fGraphEvents->SetTitle("Events with time"); + fGraphEvents->Draw("AL"); + +} + + +void MHGausEvents::DrawPowerSpectrum(TVirtualPad &pad, Int_t i) +{ + + if (fPowerSpectrum) + { + if (!fGraphPowerSpectrum) + CreateGraphPowerSpectrum(); + + fGraphPowerSpectrum->Draw("AL"); + fGraphPowerSpectrum->SetBit(kCanDelete); + } + + pad.cd(i); + + if (fHPowerProbability && fHPowerProbability->GetEntries() > 0) + { + gPad->SetLogy(); + fHPowerProbability->Draw(); + if (fFExpFit) + { + fFExpFit->SetLineColor(IsExpFitOK() ? kGreen : kRed); + fFExpFit->Draw("same"); + } + } +} + + Index: /trunk/MagicSoft/Mars/mcalib/MHGausEvents.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHGausEvents.h (revision 3154) +++ /trunk/MagicSoft/Mars/mcalib/MHGausEvents.h (revision 3154) @@ -0,0 +1,158 @@ +#ifndef MARS_MHGausEvents +#define MARS_MHGausEvents + +#ifndef ROOT_TH1 +#include +#endif + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TVirtualPad; +class TGraph; +class TArrayF; +class TH1F; +class TH1I; +class TF1; +class MHGausEvents : public MH +{ +private: + + const static Float_t fgProbLimit; // Default probability limit for judgement if fit is OK + const static Int_t fgNDFLimit; // Default NDF limit for judgement if fit is OK + const static Int_t fgPowerProbabilityBins; // Default number of bins for the projected power spectrum + const static Int_t fgBinsAfterStripping; // Default number of bins for the Gauss Histogram after stripping off the zeros at both end + + Float_t fProbLimit; // Probability limit for judgement if fit is OK + Int_t fNDFLimit; // NDF limit for judgement if fit is OK + Int_t fPowerProbabilityBins; // number of bins for the projected power spectrum + Int_t fBinsAfterStripping; // number of bins for the Gauss Histogram after stripping off the zeros at both end + Float_t fEventFrequency; // The event frequency in Hertz (to be set) + + TF1 *fFGausFit; //-> Gauss fit for fHGausHist + TF1 *fFExpFit; //-> Exponential fit for FHPowerProbability + + TH1I *fHPowerProbability; //-> Fourier transform of fEvents projected on y-axis + + TArrayF *fPowerSpectrum; //-> Fourier transform of fEvents + + TGraph *fGraphEvents; //! TGraph to display the event array (will not be cloned!!) + TGraph *fGraphPowerSpectrum; //! TGraph to display the power spectrum array (will not be cloned!!) + + Float_t *CreateXaxis(Int_t n); // Create an x-axis for the TGraphs + + Double_t fMean; // Mean of the Gauss fit + Double_t fSigma; // Sigma of the Gauss fit + Double_t fMeanErr; // Error of the mean of the Gauss fit + Double_t fSigmaErr; // Error of the sigma of the Gauss fit + 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 + + Byte_t fFlags; // Byte to hold the bits fit result bits + + ULong_t fCurrentSize; // Current size of the array fEvents + +protected: + + TH1F fHGausHist; // Histogram which should hold the Gaussian distribution + TArrayF fEvents; // Array which holds the entries of GausHist + + // Setters + void SetPowerProbabilityBins(const Int_t nbins=fgPowerProbabilityBins) { fPowerProbabilityBins = nbins; } + void SetBinsAfterStripping(const Int_t nbins=fgBinsAfterStripping) { fBinsAfterStripping = nbins; } + + void DrawEvents(); // Draw a graph of the array fEvents + void DrawPowerSpectrum(TVirtualPad &pad, Int_t i); // Draw a graph of the array fPowerSpectrum and the hist fHPowerProbability + +public: + + MHGausEvents(const char* name=NULL, const char* title=NULL); + ~MHGausEvents(); + + virtual void Clear(Option_t *o=""); + virtual void Reset(); + + // Setters + void SetEventFrequency(const Float_t f=0) { fEventFrequency = f; } + void SetMean(const Double_t d) { fMean = d; } + void SetMeanErr(const Double_t d) { fMeanErr = d; } + void SetSigma(const Double_t d) { fSigma = d; } + void SetSigmaErr(const Double_t d) { fSigmaErr = d; } + + 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); + + // Getters + const Double_t GetMean() const { return fMean; } + const Double_t GetMeanErr() const { return fMeanErr; } + const Double_t GetSigma() const { return fSigma; } + const Double_t GetSigmaErr() const { return fSigmaErr; } + const Double_t GetChiSquare() const; + const Double_t GetProb() const { return fProb; } + const Int_t GetNdf() const; + + const Double_t GetSlope() const; + const Double_t GetOffset() const; + const Double_t GetExpChiSquare() const; + const Double_t GetExpProb() const; + const Int_t GetExpNdf() const; + + TH1F *GetHGausHist() { return &fHGausHist; } + const TH1F *GetHGausHist() const { return &fHGausHist; } + + TArrayF *GetEvents() { return &fEvents; } + const TArrayF *GetEvents() const { return &fEvents; } + + TArrayF *GetPowerSpectrum() { return fPowerSpectrum; } + const TArrayF *GetPowerSpectrum() const { return fPowerSpectrum; } + + TF1 *GetFGausFit() { return fFGausFit; } + const TF1 *GetFGausFit() const { return fFGausFit; } + + TH1I *GetHPowerProbability() { return fHPowerProbability; } + const TH1I *GetHPowerProbability() const { return fHPowerProbability; } + + TF1 *GetFExpFit() { return fFExpFit; } + const TF1 *GetFExpFit() const { return fFExpFit; } + + TGraph *GetGraphEvents() { return fGraphEvents; } + const TGraph *GetGraphEvents() const { return fGraphEvents; } + + TGraph *GetGraphPowerSpectrum() { return fGraphPowerSpectrum; } + const TGraph *GetGraphPowerSpectrum() const { return fGraphPowerSpectrum; } + + const Bool_t IsGausFitOK() const; + const Bool_t IsExpFitOK() const; + const Bool_t IsEmpty() const; + const Bool_t IsFourierSpectrumOK() const; + + // Fill + void FillArray(Float_t f); // Fill only the array fEvents + Bool_t FillHist(Float_t f); // Fill only the histogram HGausHist + Bool_t FillHistAndArray(Float_t f); // Fill bothe the array fEvents and the histogram HGausHist + + // Fits + Bool_t FitGaus(Option_t *option="RQ0"); // Fit the histogram HGausHist with a Gaussian + + // Draws + virtual void Draw(Option_t *option=""); // Default Draw + + // Prints + virtual void Print(const Option_t *o="") const; // Default Print + + // Miscelleaneous + void CreateFourierSpectrum(); // Create the fourier spectrum out of fEvents + void CreateGraphEvents(); // Create the TGraph fGraphEvents of fEvents + void CreateGraphPowerSpectrum(); // Create the TGraph fGraphPowerSpectrum out of fPowerSpectrum + + ClassDef(MHGausEvents, 1) // Histograms for events with Gaussian distributed values +}; + +#endif Index: /trunk/MagicSoft/Mars/mcalib/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mcalib/Makefile (revision 3153) +++ /trunk/MagicSoft/Mars/mcalib/Makefile (revision 3154) @@ -48,5 +48,5 @@ MHCalibrationPixel.cc \ MMcCalibrationCalc.cc \ - MHGausEvent.cc + MHGausEvents.cc SRCS = $(SRCFILES)