Index: trunk/MagicSoft/Mars/mhcalib/HCalibLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/HCalibLinkDef.h (revision 4930) +++ trunk/MagicSoft/Mars/mhcalib/HCalibLinkDef.h (revision 4930) @@ -0,0 +1,23 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MHCalibrationCam+; +#pragma link C++ class MHCalibrationPix+; +#pragma link C++ class MHCalibrationChargeCam+; +#pragma link C++ class MHCalibrationChargePix+; +#pragma link C++ class MHCalibrationChargeBlindCam+; +#pragma link C++ class MHCalibrationChargeBlindPix+; +#pragma link C++ class MHCalibrationChargePINDiode+; +#pragma link C++ class MHCalibrationRelTimePix+; +#pragma link C++ class MHCalibrationRelTimeCam+; +#pragma link C++ class MHCalibrationTestCam+; +#pragma link C++ class MHCalibrationTestPix+; +#pragma link C++ class MHCalibrationTestTimeCam+; +#pragma link C++ class MHCalibrationTestTimePix+; + +#pragma link C++ class MHGausEvents++; + +#endif Index: trunk/MagicSoft/Mars/mhcalib/MHGausEvents.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHGausEvents.cc (revision 4930) +++ trunk/MagicSoft/Mars/mhcalib/MHGausEvents.cc (revision 4930) @@ -0,0 +1,912 @@ +/* ======================================================================== *\ +! +! * +! * 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 to follow a Gaussian distribution +// with time, e.g. calibration events, observables containing white noise, ... +// +// 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. 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 +// 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, a flag IsGausFitOK() is set. +// - One can repeat the fit within a given amount of sigmas from the previous mean +// (specified by the variables fPickupLimit and fBlackoutLimit) 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: +// ========================== +// +// It is created with 0 entries and not expanded unless FillArray() or FillHistAndArray() +// are 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. +// Be aware that you might lose information at the end of your analysis. +// - Calling the function CreateFourierSpectrum() 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 IsExpFitOK() is set. +// - The flag IsFourierSpectrumOK() is set accordingly to IsExpFitOK(). +// 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. +// +// To see an example, have a look at: Draw() +// +////////////////////////////////////////////////////////////////////////////// +#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 Int_t MHGausEvents::fgNDFLimit = 2; +const Float_t MHGausEvents::fgProbLimit = 0.001; +const Int_t MHGausEvents::fgPowerProbabilityBins = 30; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// Sets: +// - the default Probability Bins for fPowerProbabilityBins (fgPowerProbabilityBins) +// - the default Probability Limit for fProbLimit (fgProbLimit) +// - the default NDF Limit for fNDFLimit (fgNDFLimit) +// - 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: +// - fEvents to 0 entries +// - fHGausHist() +// - all other pointers to NULL +// - all variables to 0. +// - all flags to kFALSE +// +MHGausEvents::MHGausEvents(const char *name, const char *title) + : fHPowerProbability(NULL), + fPowerSpectrum(NULL), + fFGausFit(NULL), fFExpFit(NULL), + fFirst(0.), + fGraphEvents(NULL), fGraphPowerSpectrum(NULL), + fLast(100.), fNbins(100) +{ + + fName = name ? name : "MHGausEvents"; + fTitle = title ? title : "Events with expected Gaussian distributions"; + + Clear(); + + SetBinsAfterStripping(); + SetNDFLimit(); + SetPowerProbabilityBins(); + SetProbLimit(); + + fHGausHist.SetName("HGausHist"); + fHGausHist.SetTitle("Histogram of Events with Gaussian Distribution"); + // important, other ROOT will not draw the axes: + fHGausHist.UseCurrentStyle(); + fHGausHist.SetDirectory(NULL); + fHGausHist.GetYaxis()->CenterTitle(); +} + + + +// -------------------------------------------------------------------------- +// +// Default Destructor. +// +// Deletes (if Pointer is not NULL): +// +// - fHPowerProbability +// - fFGausFit +// - fFExpFit +// - fPowerSpectrum +// - fGraphEvents +// - fGraphPowerSpectrum +// +MHGausEvents::~MHGausEvents() +{ + + // delete histograms +// if (fHPowerProbability) +// if (gROOT->FindObject(fHPowerProbability->GetName())) +// 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; +} + +// -------------------------------------------------------------------------- +// +// Default Clear(), can be overloaded. +// +// Sets: +// - all other pointers to NULL +// - all variables to 0. +// - all flags to kFALSE +// +// Deletes (if not NULL): +// - all pointers +// +void MHGausEvents::Clear(Option_t *o) +{ + + SetGausFitOK ( kFALSE ); + SetExpFitOK ( kFALSE ); + SetFourierSpectrumOK( kFALSE ); + SetExcluded ( kFALSE ); + + fMean = 0.; + fSigma = 0.; + fMeanErr = 0.; + fSigmaErr = 0.; + fProb = 0.; + + fCurrentSize = 0; + + if (fHPowerProbability) + { + if (gROOT->FindObject(fHPowerProbability->GetName())) + delete fHPowerProbability; + fHPowerProbability = NULL; + } + + // delete fits + if (fFGausFit) + { + delete fFGausFit; + fFGausFit = NULL; + } + + if (fFExpFit) + { + delete fFExpFit; + fFExpFit = NULL; + } + + // delete arrays + if (fPowerSpectrum) + { + delete fPowerSpectrum; + fPowerSpectrum = NULL; + } + + // delete graphs + if (fGraphEvents) + { + delete fGraphEvents; + fGraphEvents = NULL; + } + + if (fGraphPowerSpectrum) + { + delete fGraphPowerSpectrum; + fGraphPowerSpectrum = NULL; + } +} + +TObject *MHGausEvents::Clone(const char *name) const +{ + + MHGausEvents &gaus = *new MHGausEvents; + + gaus.fBinsAfterStripping = fBinsAfterStripping; + gaus.fCurrentSize = fCurrentSize; + gaus.fFlags = fFlags; + gaus.fPowerProbabilityBins = fPowerProbabilityBins; + /* + if (fHPowerProbability) + gaus.fHPowerProbability=(TH1I*)fHPowerProbability->Clone(); + + if (fPowerSpectrum) + gaus.fPowerSpectrum = new TArrayF(*fPowerSpectrum); + */ + + gaus.fEvents.Set(fEvents.GetSize()); + gaus.fEvents = fEvents; + /* + if (fFGausFit) + gaus.fFGausFit=(TF1*)fFGausFit->Clone(); + if (fFExpFit) + gaus.fFExpFit=(TF1*)fFExpFit->Clone(); + */ + + gaus.fFirst = fFirst; + + /* + if (fGraphEvents) + gaus.fGraphEvents=(TGraph*)fGraphEvents->Clone(); + if (fGraphPowerSpectrum) + gaus.fGraphPowerSpectrum=(TGraph*)fGraphPowerSpectrum->Clone(); + */ + + Copy(gaus.fHGausHist); + + gaus.fLast = fLast; + gaus.fMean = fMean; + gaus.fMeanErr = fMeanErr; + gaus.fNbins = fNbins; + gaus.fNDFLimit = fNDFLimit; + gaus.fSigma = fSigma; + gaus.fSigma = fSigmaErr; + gaus.fProb = fProb; + gaus.fProbLimit = fProbLimit; + + return &gaus; +} + + +// ----------------------------------------------------------------------------- +// +// Create the x-axis for the event graph +// +Float_t *MHGausEvents::CreateEventXaxis(Int_t n) +{ + + Float_t *xaxis = new Float_t[n]; + + for (Int_t i=0;iSetXTitle("P(f)"); + fHPowerProbability->SetYTitle("Counts"); + fHPowerProbability->GetYaxis()->CenterTitle(); + fHPowerProbability->SetDirectory(NULL); + fHPowerProbability->SetBit(kCanDelete); + // + // 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(kTRUE); + + // + // For the moment, this is the only check, later we can add more... + // + SetFourierSpectrumOK(IsExpFitOK()); + + return; +} + +// ---------------------------------------------------------------------------------- +// +// 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(); + + if (n==0) + return; + + fGraphEvents = new TGraph(n,CreateEventXaxis(n),fEvents.GetArray()); + fGraphEvents->SetTitle("Evolution of Events with time"); + fGraphEvents->GetYaxis()->SetTitle(fHGausHist.GetXaxis()->GetTitle()); + fGraphEvents->GetYaxis()->CenterTitle(); +} + +// ---------------------------------------------------------------------------------- +// +// 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,CreatePSDXaxis(n),fPowerSpectrum->GetArray()); + fGraphPowerSpectrum->SetTitle("Power Spectrum Density"); + fGraphPowerSpectrum->GetYaxis()->SetTitle("P(f)"); + fGraphPowerSpectrum->GetYaxis()->CenterTitle(); + +} + + +// ----------------------------------------------------------------------------- +// +// Create the x-axis for the event graph +// +Float_t *MHGausEvents::CreatePSDXaxis(Int_t n) +{ + + Float_t *xaxis = new Float_t[n]; + + for (Int_t i=0;i +*/ +//End_Html +// +void MHGausEvents::Draw(const Option_t *opt) +{ + + TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this,600, 900); + + TString option(opt); + option.ToLower(); + + Int_t win = 1; + + if (option.Contains("events")) + { + option.ReplaceAll("events",""); + win += 1; + } + if (option.Contains("fourier")) + { + option.ReplaceAll("fourier",""); + win += 2; + } + + pad->SetBorderMode(0); + pad->Divide(1,win); + pad->cd(1); + + if (!IsEmpty() && !IsOnlyOverflow() && !IsOnlyUnderflow()) + gPad->SetLogy(); + + gPad->SetTicks(); + + fHGausHist.Draw(option); + + 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(); + + if (!fGraphEvents) + return; + + 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"); + } + } +} + + +// -------------------------------------------------------------------------- +// +// Fill fEvents with f +// If size of fEvents is 0, initializes it to 512 +// If size of fEvents is smaller than fCurrentSize, double the size +// Increase fCurrentSize by 1 +// +void MHGausEvents::FillArray(const Float_t f) +{ + if (fEvents.GetSize() == 0) + fEvents.Set(512); + + if (fCurrentSize >= fEvents.GetSize()) + fEvents.Set(fEvents.GetSize()*2); + + fEvents.AddAt(f,fCurrentSize++); +} + + +// -------------------------------------------------------------------------- +// +// Fills fHGausHist with f +// Returns kFALSE, if overflow or underflow occurred, else kTRUE +// +Bool_t MHGausEvents::FillHist(const Float_t f) +{ + return fHGausHist.Fill(f) > -1; +} + +// -------------------------------------------------------------------------- +// +// Executes: +// - FillArray() +// - FillHist() +// +Bool_t MHGausEvents::FillHistAndArray(const Float_t f) +{ + + FillArray(f); + return FillHist(f); +} + +// ------------------------------------------------------------------- +// +// 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, const Double_t xmin, const Double_t xmax) +{ + + if (IsGausFitOK()) + return kTRUE; + + // + // First, strip the zeros from the edges which contain only zeros and rebin + // to fBinsAfterStripping bins. If fBinsAfterStripping is 0, reduce bins only by + // a factor 10. + // + // (ATTENTION: The Chisquare method is more sensitive, + // the _less_ bins, you have!) + // + StripZeros(&fHGausHist, + fBinsAfterStripping ? fBinsAfterStripping + : (fNbins > 1000 ? fNbins/10 : 0)); + + TAxis *axe = fHGausHist.GetXaxis(); + // + // Get the fitting ranges + // + Axis_t rmin = ((xmin==0.) && (xmax==0.)) ? fHGausHist.GetBinCenter(axe->GetFirst()) : xmin; + Axis_t rmax = ((xmin==0.) && (xmax==0.)) ? fHGausHist.GetBinCenter(axe->GetLast()) : xmax; + + // + // First guesses for the fit (should be as close to reality as possible, + // + const Stat_t entries = fHGausHist.Integral(axe->FindBin(rmin),axe->FindBin(rmax),"width"); + const Double_t mu_guess = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin()); + const Double_t sigma_guess = fHGausHist.GetRMS(); + 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 " + << "in: " << fName << endl; + return kFALSE; + } + + fFGausFit->SetParameters(area_guess,mu_guess,sigma_guess); + fFGausFit->SetParNames("Area","#mu","#sigma"); + fFGausFit->SetParLimits(0,0.,area_guess*25.); + 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 (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; +} + + +const Double_t MHGausEvents::GetChiSquare() const +{ + return ( fFGausFit ? fFGausFit->GetChisquare() : 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 Int_t MHGausEvents::GetNdf() const +{ + return ( fFGausFit ? fFGausFit->GetNDF() : 0); +} + +const Double_t MHGausEvents::GetOffset() const +{ + return ( fFExpFit ? fFExpFit->GetParameter(0) : 0.); +} + + + +// -------------------------------------------------------------------------- +// +// If fFExpFit exists, returns fit parameter 1 (Slope of Exponential fit), +// otherwise 0. +// +const Double_t MHGausEvents::GetSlope() const +{ + return ( fFExpFit ? fFExpFit->GetParameter(1) : 0.); +} + +// -------------------------------------------------------------------------- +// +// Default InitBins, can be overloaded. +// +// Executes: +// - fHGausHist.SetBins(fNbins,fFirst,fLast) +// +void MHGausEvents::InitBins() +{ + fHGausHist.SetBins(fNbins,fFirst,fLast); +} + +// -------------------------------------------------------------------------- +// +// Return kFALSE if number of entries is 0 +// +const Bool_t MHGausEvents::IsEmpty() const +{ + return !(fHGausHist.GetEntries()); +} + +// -------------------------------------------------------------------------- +// +// Return kTRUE if number of entries is bin content of fNbins+1 +// +const Bool_t MHGausEvents::IsOnlyOverflow() const +{ + return fHGausHist.GetEntries() == fHGausHist.GetBinContent(fNbins+1); +} + +// -------------------------------------------------------------------------- +// +// Return kTRUE if number of entries is bin content of 0 +// +const Bool_t MHGausEvents::IsOnlyUnderflow() const +{ + return fHGausHist.GetEntries() == fHGausHist.GetBinContent(0); +} + + +const Bool_t MHGausEvents::IsExcluded() const +{ + return TESTBIT(fFlags,kExcluded); +} + + +const Bool_t MHGausEvents::IsExpFitOK() const +{ + return TESTBIT(fFlags,kExpFitOK); +} + +const Bool_t MHGausEvents::IsFourierSpectrumOK() const +{ + return TESTBIT(fFlags,kFourierSpectrumOK); +} + + +const Bool_t MHGausEvents::IsGausFitOK() const +{ + return TESTBIT(fFlags,kGausFitOK); +} + + +// ----------------------------------------------------------------------------------- +// +// A default print +// +void MHGausEvents::Print(const Option_t *o) const +{ + + *fLog << all << endl; + *fLog << all << "Results of the Gauss Fit in: " << fName << 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; + +} + + +// -------------------------------------------------------------------------- +// +// Default Reset(), can be overloaded. +// +// Executes: +// - Clear() +// - fHGausHist.Reset() +// - fEvents.Set(0) +// +void MHGausEvents::Reset() +{ + + Clear(); + fHGausHist.Reset(); + fEvents.Set(0); + +} + +// -------------------------------------------------------------------------- +// +// Set Excluded bit from outside +// +void MHGausEvents::SetExcluded(const Bool_t b) +{ + b ? SETBIT(fFlags,kExcluded) : CLRBIT(fFlags,kExcluded); +} + + +// ------------------------------------------------------------------- +// +// 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); +} + + +// ------------------------------------------------------------------- +// +// 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); + +} + Index: trunk/MagicSoft/Mars/mhcalib/MHGausEvents.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHGausEvents.h (revision 4930) +++ trunk/MagicSoft/Mars/mhcalib/MHGausEvents.h (revision 4930) @@ -0,0 +1,161 @@ +#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 Int_t fgNDFLimit; //! Default for fNDFLimit (now set to: 2) + const static Float_t fgProbLimit; //! Default for fProbLimit (now set to: 0.001) + const static Int_t fgPowerProbabilityBins; //! Default for fPowerProbabilityBins (now set to: 20) + + Int_t fBinsAfterStripping; // Bins for the Gauss Histogram after stripping off the zeros at both ends + Int_t fCurrentSize; // Current size of the array fEvents + Byte_t fFlags; // Bit field for the fit result bits + Int_t fPowerProbabilityBins; // Bins for the projected power spectrum + + TH1I *fHPowerProbability; // Fourier transform of fEvents projected on y-axis + TArrayF *fPowerSpectrum; // Fourier transform of fEvents + + enum { kGausFitOK, + kExpFitOK, + kFourierSpectrumOK, + kExcluded }; // Bits for information about fit results + +protected: + + TArrayF fEvents; // Array which holds the entries of GausHist + TF1 *fFGausFit; // Gauss fit for fHGausHist + TF1 *fFExpFit; // Exponential fit for FHPowerProbability + Axis_t fFirst; // Lower histogram edge for fHGausHist (used by InitBins()) + 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!!) + TH1F fHGausHist; // Histogram to hold the Gaussian distribution + Axis_t fLast; // Upper histogram edge for fHGausHist (used by InitBins()) + Double_t fMean; // Mean of the Gauss fit + Double_t fMeanErr; // Error of the mean of the Gauss fit + Int_t fNbins; // Number histogram bins for fHGausHist (used by InitBins()) + Int_t fNDFLimit; // NDF limit for judgement if fit is OK + Double_t fSigma; // Sigma of the Gauss fit + Double_t fSigmaErr; // Error of the sigma of the Gauss fit + Double_t fProb; // Probability of the Gauss fit + Float_t fProbLimit; // Probability limit for judgement if fit is OK + + virtual Float_t *CreateEventXaxis(Int_t n); // Create an x-axis for the Event TGraphs + virtual Float_t *CreatePSDXaxis(Int_t n); // Create an x-axis for the PSD TGraphs + virtual void CreateGraphEvents(); // Create the TGraph fGraphEvents of fEvents + virtual void CreateGraphPowerSpectrum(); // Create the TGraph fGraphPowerSpectrum out of fPowerSpectrum + + void DrawPowerSpectrum(TVirtualPad &pad, Int_t i); // Draw graph of fPowerSpectrum and fHPowerProbability + + // Setters + void SetBinsAfterStripping ( const Int_t nbins=0 ) { fBinsAfterStripping =nbins; } + void SetPowerProbabilityBins ( const Int_t nbins=fgPowerProbabilityBins ) { fPowerProbabilityBins=nbins; } + +public: + + MHGausEvents(const char* name=NULL, const char* title=NULL); + ~MHGausEvents(); + + TObject *Clone(const char* name="") const; + + void Clear(Option_t *o=""); + void Reset(); + + // Draws + void Draw(Option_t *option=""); // Default Draw + void DrawEvents(); // Draw graph of fEvents + + // Inits + virtual void InitBins(); + + // Getters + const Double_t GetChiSquare() const; + const Double_t GetExpChiSquare() const; + const Int_t GetExpNdf() const; + const Double_t GetExpProb() const; + TArrayF *GetEvents() { return &fEvents; } + const TArrayF *GetEvents() const { return &fEvents; } + TF1 *GetFExpFit() { return fFExpFit; } + const TF1 *GetFExpFit() const { return fFExpFit; } + TF1 *GetFGausFit() { return fFGausFit; } + const TF1 *GetFGausFit() const { return fFGausFit; } + TGraph *GetGraphEvents() { return fGraphEvents; } + const Double_t GetFirst() const { return fFirst; } + const Double_t GetLast () const { return fLast ; } + const TGraph *GetGraphEvents() const { return fGraphEvents; } + TGraph *GetGraphPowerSpectrum() { return fGraphPowerSpectrum; } + const TGraph *GetGraphPowerSpectrum() const { return fGraphPowerSpectrum; } + TH1F *GetHGausHist() { return &fHGausHist; } + const TH1F *GetHGausHist() const { return &fHGausHist; } + TH1I *GetHPowerProbability() { return fHPowerProbability; } + const TH1I *GetHPowerProbability() const { return fHPowerProbability; } + const Double_t GetMean() const { return fMean; } + const Double_t GetMeanErr() const { return fMeanErr; } + const Int_t GetNdf() const; + const Double_t GetOffset() const; + TArrayF *GetPowerSpectrum() { return fPowerSpectrum; } + const TArrayF *GetPowerSpectrum() const { return fPowerSpectrum; } + const Double_t GetProb() const { return fProb; } + const Double_t GetSigma() const { return fSigma; } + const Double_t GetSigmaErr() const { return fSigmaErr; } + const Double_t GetSlope() const; + + const Bool_t IsExcluded() const; + const Bool_t IsExpFitOK() const; + const Bool_t IsEmpty() const; + const Bool_t IsFourierSpectrumOK() const; + const Bool_t IsGausFitOK() const; + const Bool_t IsOnlyOverflow() const; + const Bool_t IsOnlyUnderflow() const; + + // Fill + void FillArray ( const Float_t f ); // Fill only the array fEvents + Bool_t FillHist ( const Float_t f ); // Fill only the histogram HGausHist + Bool_t FillHistAndArray( const Float_t f ); // Fill bothe the array fEvents and the histogram HGausHist + + // Fits + Bool_t FitGaus( Option_t *option="RQ0", + const Double_t xmin=0., + const Double_t xmax=0.); // Fit the histogram HGausHist with a Gaussian + + // Prints + virtual void Print(const Option_t *o="") const; // Default Print + + // Setters + void SetExcluded ( const Bool_t b=kTRUE ); + void SetExpFitOK ( const Bool_t b=kTRUE ); + void SetFourierSpectrumOK( const Bool_t b=kTRUE ); + void SetGausFitOK ( const Bool_t b=kTRUE ); + void SetLast ( const Double_t d ) { fLast = d; } + void SetFirst ( const Double_t d ) { fFirst = d; } + void SetMean ( const Double_t d ) { fMean = d; } + void SetMeanErr ( const Double_t d ) { fMeanErr = d; } + void SetNbins ( const Int_t i ) { fNbins = i; } + void SetNDFLimit ( const Int_t lim=fgNDFLimit ) { fNDFLimit = lim; } + void SetProb ( const Double_t d ) { fProb = d; } + void SetProbLimit ( const Float_t lim=fgProbLimit ) { fProbLimit = lim; } + void SetSigma ( const Double_t d ) { fSigma = d; } + void SetSigmaErr ( const Double_t d ) { fSigmaErr = d; } + + void CreateFourierSpectrum(); // Create the fourier spectrum out of fEvents + + ClassDef(MHGausEvents, 1) // Base class for events with Gaussian distributed values +}; + +#endif Index: trunk/MagicSoft/Mars/mhcalib/Makefile =================================================================== --- trunk/MagicSoft/Mars/mhcalib/Makefile (revision 4930) +++ trunk/MagicSoft/Mars/mhcalib/Makefile (revision 4930) @@ -0,0 +1,52 @@ +################################################################## +# +# makefile +# +# for the MARS software +# +################################################################## +include ../Makefile.conf.$(OSTYPE) +include ../Makefile.conf.general + +#------------------------------------------------------------------------------ + +# +# Handling name of the Root Dictionary Files +# +CINT = HCalib + +# +# connect the include files defined in the config.mk file +# +INCLUDES = -I. -I../mbase -I../mhbase \ + -I../manalysis -I../mraw -I../mtools -I../mmc -I../mhist \ + -I../mimage -I../msignal -I../mbadpixels -I../mpedestal + +# mhbase: MBinning MH +# manalysis: MExtractedSignal, MCerPhotEvt (move to mcalib?) +# mraw: MRawRunHeader, MRawEvtHeader, MRawEvtPixelIter (3xMCalibrationCalc) +# mmc: MMcFadcHeader, MMcEvt +# mimage MHillas + +SRCFILES = MHCalibrationChargeBlindPix.cc \ + MHCalibrationChargeBlindCam.cc \ + MHCalibrationChargePix.cc \ + MHCalibrationCam.cc \ + MHCalibrationPix.cc \ + MHCalibrationChargeCam.cc \ + MHCalibrationChargePINDiode.cc \ + MHCalibrationRelTimeCam.cc \ + MHCalibrationRelTimePix.cc \ + MHCalibrationTestCam.cc \ + MHCalibrationTestPix.cc \ + MHCalibrationTestTimeCam.cc \ + MHCalibrationTestTimePix.cc \ + MHGausEvents.cc + +############################################################ + +all: $(OBJS) + +include ../Makefile.rules + +mrproper: clean rmbak