Index: /trunk/MagicSoft/Mars/Changelog =================================================================== --- /trunk/MagicSoft/Mars/Changelog (revision 3107) +++ /trunk/MagicSoft/Mars/Changelog (revision 3108) @@ -40,4 +40,10 @@ of inverted events. + * mcalib/MHGausEvent.[h,cc] + * mcalib/Makefile + * mcalib/CalibLinkDef.h + - new base class for all kind of events with variables having a + Gaussian distribution. Will serve as base class for the + existing and new calibration histogram classes. Index: /trunk/MagicSoft/Mars/mcalib/CalibLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/CalibLinkDef.h (revision 3107) +++ /trunk/MagicSoft/Mars/mcalib/CalibLinkDef.h (revision 3108) @@ -19,3 +19,5 @@ #pragma link C++ class MMcCalibrationCalc++; +#pragma link C++ class MHGausEvent++; + #endif Index: /trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.cc (revision 3107) +++ /trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.cc (revision 3108) @@ -465,6 +465,4 @@ // We need a dedicated signal extractor for the blind pixel // - MPedestalPix &ped = (*fPedestals)[fBlindPixelId]; - Int_t diff = 0; Index: /trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc (revision 3107) +++ /trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc (revision 3108) @@ -88,6 +88,6 @@ const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmpError = 0.002; -const Int_t MHCalibrationBlindPixel::fPSDNbins = 30; - +const Int_t MHCalibrationBlindPixel::fPSDNbins = 50; +const Int_t MHCalibrationBlindPixel::fPulserFrequency = 500; // -------------------------------------------------------------------------- // @@ -113,6 +113,6 @@ // Create a large number of bins, later we will rebin - fBlindPixelChargefirst = -200.; - fBlindPixelChargelast = 800.; + fBlindPixelChargefirst = -400.; + fBlindPixelChargelast = 1600.; fHBlindPixelCharge = new TH1F("HBlindPixelCharge","Distribution of Summed FADC Slices", @@ -173,6 +173,6 @@ fCurrentSize = 1024; - fBlindPixelChargefirst = -200.; - fBlindPixelChargelast = 800.; + fBlindPixelChargefirst = -400.; + fBlindPixelChargelast = 1600.; fLambda = 0.; @@ -259,4 +259,8 @@ return IsOscillating(); + // This cuts only the non-used zero's, but MFFT will later cut the rest + CutArrayBorder(fHiGains); + CutArrayBorder(fLoGains); + MFFT fourier; @@ -267,4 +271,5 @@ "PSDProjectionBlindPixel", "Power Spectrum Density Projection HiGain Blind Pixel"); + // // First guesses for the fit (should be as close to reality as possible, @@ -280,5 +285,5 @@ fPSDExpFit->SetParNames("Offset","Slope"); fPSDExpFit->SetParLimits(0,offset_guess/2.,2.*offset_guess); - fPSDExpFit->SetParLimits(1,slope_guess/2.,2.*slope_guess); + fPSDExpFit->SetParLimits(1,slope_guess/1.5,1.5*slope_guess); fPSDExpFit->SetRange(0.,xmax); @@ -307,5 +312,6 @@ for (Int_t i=0;iAddAt((Float_t)i,i); + for (Int_t i=0;iAddAt((Float_t)(fPulserFrequency*i)/2./n,i); } @@ -317,5 +323,5 @@ fChargeXaxis = new TArrayF(n); for (Int_t i=0;iAddAt((Float_t)i,i); + fChargeXaxis->AddAt((Float_t)i/fPulserFrequency,i); return; } @@ -328,5 +334,5 @@ if (diff < 0) for (Int_t i=n;iAddAt((Float_t)i,i); + fChargeXaxis->AddAt((Float_t)i/fPulserFrequency,i); } @@ -498,5 +504,5 @@ TCanvas *c = MakeDefCanvas(this,550,700); - c->Divide(2,4); + c->Divide(2,5); gROOT->SetSelectedPad(NULL); @@ -539,5 +545,7 @@ fChargeXaxis->GetArray(), fHiGains->GetArray()); - gr1->SetTitle("Evolution of HiGain charges with event number"); + gr1->SetTitle("Evolution of HiGain charges with time"); + gr1->GetXaxis()->SetTitle("Time [s]"); + gr1->GetYaxis()->SetTitle("Sum FADC slices"); gr1->SetBit(kCanDelete); gr1->Draw("AL"); @@ -551,5 +559,7 @@ fChargeXaxis->GetArray(), fLoGains->GetArray()); - gr2->SetTitle("Evolution of HiGain charges with event number"); + gr2->SetTitle("Evolution of LoGain charges with time"); + gr2->GetXaxis()->SetTitle("Time [s]"); + gr2->GetYaxis()->SetTitle("Sum FADC slices"); gr2->SetBit(kCanDelete); gr2->Draw("AL"); @@ -571,4 +581,6 @@ TGraph *gr3 = new TGraph(fPSDXaxis->GetSize(),fPSDXaxis->GetArray(),array->GetArray()); gr3->SetTitle("Power Spectrum Density"); + gr3->GetXaxis()->SetTitle("Freq. [Hz]"); + gr3->GetYaxis()->SetTitle("P(f)"); gr3->SetBit(kCanDelete); gr3->Draw("AL"); @@ -586,4 +598,6 @@ fHPSD->Draw(opt); + fHPSD->GetXaxis()->SetTitle("P(f)"); + fHPSD->GetYaxis()->SetTitle("Counts"); if (fPSDExpFit) @@ -593,7 +607,20 @@ } + c->cd(9); + + array = fPSDLoGain; + + if (!fPSDXaxis) + CreatePSDXaxis(array->GetSize()); + + TGraph *gr4 = new TGraph(fPSDXaxis->GetSize(),fPSDXaxis->GetArray(),array->GetArray()); + gr4->SetTitle("Power Spectrum Density"); + gr4->GetXaxis()->SetTitle("Freq. [Hz]"); + gr4->GetYaxis()->SetTitle("P(f)"); + gr4->SetBit(kCanDelete); + gr4->Draw("AL"); + c->Modified(); c->Update(); - } @@ -892,5 +919,5 @@ { - Int_t nbins = 20; + Int_t nbins = 30; CutEdges(fHBlindPixelCharge,nbins); Index: /trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h (revision 3107) +++ /trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h (revision 3108) @@ -26,4 +26,5 @@ static const Int_t fPSDNbins; + static const Int_t fPulserFrequency; TH1F* fHBlindPixelCharge; // Histogram with the single Phe spectrum Index: /trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc (revision 3107) +++ /trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc (revision 3108) @@ -62,4 +62,5 @@ const Int_t MHCalibrationPixel::fPSDNbins = 30; +const Int_t MHCalibrationPixel::fPulserFrequency = 500; // -------------------------------------------------------------------------- @@ -316,5 +317,5 @@ for (Int_t i=0;iAddAt((Float_t)i,i); + fPSDXaxis->AddAt((Float_t)(fPulserFrequency*i)/2./n,i); } @@ -326,5 +327,5 @@ fChargeXaxis = new TArrayF(n); for (Int_t i=0;iAddAt((Float_t)i,i); + fChargeXaxis->AddAt((Float_t)i/fPulserFrequency,i); return; } @@ -646,5 +647,7 @@ fChargeXaxis->GetArray(), fHiGains->GetArray()); - gr1->SetTitle("Evolution of HiGain charges with event number"); + gr1->SetTitle("Evolution of HiGain charges with time"); + gr1->GetXaxis()->SetTitle("Time [s]"); + gr1->GetYaxis()->SetTitle("Sum FADC slices"); gr1->SetBit(kCanDelete); gr1->Draw("AL"); @@ -657,5 +660,7 @@ fChargeXaxis->GetArray(), fLoGains->GetArray()); - gr2->SetTitle("Evolution of LoGain charges with event number"); + gr2->SetTitle("Evolution of LoGain charges with time"); + gr2->GetXaxis()->SetTitle("Time [s]"); + gr2->GetYaxis()->SetTitle("Sum FADC slices"); gr2->SetBit(kCanDelete); gr2->Draw("AL"); @@ -689,4 +694,6 @@ TGraph *gr3 = new TGraph(fPSDXaxis->GetSize(),fPSDXaxis->GetArray(),array->GetArray()); gr3->SetTitle(title.Data()); + gr3->GetXaxis()->SetTitle("Freq. [Hz]"); + gr3->GetYaxis()->SetTitle("P(f)"); gr3->SetBit(kCanDelete); gr3->Draw("AL"); @@ -705,4 +712,6 @@ fHPSD->Draw(opt); + fHPSD->GetXaxis()->SetTitle("P(f)"); + fHPSD->GetYaxis()->SetTitle("Counts"); if (fPSDExpFit) Index: /trunk/MagicSoft/Mars/mcalib/MHGausEvent.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHGausEvent.cc (revision 3108) +++ /trunk/MagicSoft/Mars/mcalib/MHGausEvent.cc (revision 3108) @@ -0,0 +1,318 @@ +/* ======================================================================== *\ +! +! * +! * 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::CutEdges(fEvents); + + MFFT fourier; + + fPowerSpectrum = fourier.PowerSpectrumDensity(fEvents); + + fHPowerProbability = MH::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::CutEdges(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: /trunk/MagicSoft/Mars/mcalib/MHGausEvent.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHGausEvent.h (revision 3108) +++ /trunk/MagicSoft/Mars/mcalib/MHGausEvent.h (revision 3108) @@ -0,0 +1,97 @@ +#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/Makefile =================================================================== --- /trunk/MagicSoft/Mars/mcalib/Makefile (revision 3107) +++ /trunk/MagicSoft/Mars/mcalib/Makefile (revision 3108) @@ -45,5 +45,6 @@ MHCalibrationPINDiode.cc \ MHCalibrationPixel.cc \ - MMcCalibrationCalc.cc + MMcCalibrationCalc.cc \ + MHGausEvent.cc SRCS = $(SRCFILES)