/* ======================================================================== *\ ! ! * ! * 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 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // MHCalibrationPix // // A base class for events which are believed to follow a Gaussian distribution // with time, e.g. calibration events, observables containing white noise, ... // // MHCalibrationPix derives from MHGausEvents, thus all features of // MHGausEvents can be used by a class deriving from MHCalibrationPix // // See also: MHGausEvents // ////////////////////////////////////////////////////////////////////////////// #include "MHCalibrationPix.h" #include #include #include #include "MLog.h" #include "MLogManip.h" ClassImp(MHCalibrationPix); using namespace std; const Float_t MHCalibrationPix::fgBlackoutLimit = 5.; const Float_t MHCalibrationPix::fgPickupLimit = 5.; // -------------------------------------------------------------------------- // // Default Constructor. // Sets: // - the default number for fPickupLimit (fgPickupLimit) // - the default number for fBlackoutLimit (fgBlackoutLimit) // // Initializes: // - all variables to 0. // MHCalibrationPix::MHCalibrationPix(const char *name, const char *title) : fEventFrequency(0), fPixId(-1) { fName = name ? name : "MHCalibrationPix"; fTitle = title ? title : "Calibration histogram events"; Clear(); SetBlackoutLimit(); SetPickupLimit(); } // -------------------------------------------------------------------------- // // Default Clear(), can be overloaded. // // Sets: // - all other pointers to NULL // - all variables to 0., except fPixId to -1 and keep fEventFrequency // - all flags to kFALSE // // Deletes (if not NULL): // - all pointers // void MHCalibrationPix::Clear(Option_t *o) { MHGausEvents::Clear(); fSaturated = 0; } #if 0 // -------------------------------------------------------------------------- // // ATTENTION: This nasty Clone function is necessary since the ROOT cloning // lead to crashes on SOME machines (unfortunately not mine...). // This function is a workaround in order to achieve the correct // DrawClone() behaviour. // TObject *MHCalibrationPix::Clone(const char *name) const { MHCalibrationPix &pix = *new MHCalibrationPix(name ? name : fName.Data(),fTitle.Data()); // // Copy MHGausEvents data members // pix.fBinsAfterStripping = fBinsAfterStripping; pix.fCurrentSize = fCurrentSize; pix.fFlags = fFlags; pix.fPowerProbabilityBins = fPowerProbabilityBins; if (fHPowerProbability) pix.fHPowerProbability=(TH1I*)fHPowerProbability->Clone(); if (fPowerSpectrum) pix.fPowerSpectrum = new TArrayF(*fPowerSpectrum); pix.fEvents = fEvents; if (fFGausFit) pix.fFGausFit=(TF1*)fFGausFit->Clone(); if (fFExpFit) pix.fFExpFit=(TF1*)fFExpFit->Clone(); pix.fFirst = fFirst; if (fGraphEvents) pix.fGraphEvents=(TGraph*)fGraphEvents->Clone(); if (fGraphPowerSpectrum) pix.fGraphPowerSpectrum=(TGraph*)fGraphPowerSpectrum->Clone(); pix.fHGausHist = fHGausHist; pix.fLast = fLast; pix.fMean = fMean; pix.fMeanErr = fMeanErr; pix.fNbins = fNbins; pix.fNDFLimit = fNDFLimit; pix.fSigma = fSigma; pix.fSigmaErr = fSigmaErr; pix.fProb = fProb; pix.fProbLimit = fProbLimit; // // Copy data members // pix.fEventFrequency = fEventFrequency; pix.fBlackoutLimit = fBlackoutLimit; pix.fSaturated = fSaturated; pix.fPickupLimit = fPickupLimit; pix.fPixId = fPixId; return &pix; } #endif // ----------------------------------------------------------------------------- // // Bypasses the Gauss fit by taking mean and RMS from the histogram // // Errors are determined in the following way: // MeanErr = RMS / Sqrt(entries) // SigmaErr = RMS / (2.*Sqrt(entries) ) // void MHCalibrationPix::BypassFit() { const Stat_t entries = fHGausHist.GetEntries(); if (entries <= 0.) { *fLog << warn << GetDescriptor() << ": Cannot bypass fit. Number of entries smaller or equal 0 in pixel: " << fPixId << endl; return; } fMean = fHGausHist.GetMean(); fMeanErr = fHGausHist.GetRMS() / TMath::Sqrt(entries); fSigma = fHGausHist.GetRMS() ; fSigmaErr = fHGausHist.GetRMS() / TMath::Sqrt(entries) / 2.; } // -------------------------------------------------------------------------- // // - Set fPixId to id // // Add id to names and titles of: // - fHGausHist // void MHCalibrationPix::ChangeHistId(const 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); } // ----------------------------------------------------------------------------- // // Create the x-axis for the event graph // Float_t *MHCalibrationPix::CreateEventXaxis(Int_t n) { Float_t *xaxis = new Float_t[n]; if (fEventFrequency) for (Int_t i=0;iGetXaxis()->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 MHCalibrationPix::CreateGraphPowerSpectrum() { MHGausEvents::CreateGraphPowerSpectrum(); fGraphPowerSpectrum->GetXaxis()->SetTitle((fEventFrequency) ? "Frequency [Hz]" : "Frequency"); } // ------------------------------------------------------------------------------- // // Return the number of "blackout" events, which are events with values higher // than fBlackoutLimit sigmas from the mean // // const Double_t MHCalibrationPix::GetBlackout() const { if ((fMean == 0.) && (fSigma == 0.)) return -1.; const Int_t first = fHGausHist.GetXaxis()->GetFirst(); const Int_t last = fHGausHist.GetXaxis()->FindBin(fMean-fBlackoutLimit*fSigma); if (first >= last) return 0.; return fHGausHist.Integral(first, last, "width"); } // ------------------------------------------------------------------------------- // // Return the number of "pickup" events, which are events with values higher // than fPickupLimit sigmas from the mean // // const Double_t MHCalibrationPix::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"); } // -------------------------------------------------------------------------- // // Re-normalize the results, has to be overloaded // void MHCalibrationPix::Renorm() { } // ----------------------------------------------------------------------------- // // 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() - fBlackoutLimit * 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 MHCalibrationPix::RepeatFit(const Option_t *option) { if (IsGausFitOK()) return kTRUE; if ((fMean == 0.) && (fSigma == 0.)) return FitGaus(); // // Get new fitting ranges // Axis_t rmin = fMean - fBlackoutLimit * 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); 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; }