Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 6328) +++ trunk/MagicSoft/Mars/Changelog (revision 6329) @@ -64,4 +64,7 @@ - initialize the bad events counters to 0 also in constructor - release the limit for bad events a bit + + * mcalib/MCalibrationTestCam.h + - take out include of TArrayI (already in base class) Index: trunk/MagicSoft/Mars/mcalib/MCalibrationTestCam.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationTestCam.h (revision 6328) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationTestCam.h (revision 6329) @@ -4,8 +4,4 @@ #ifndef MARS_MCalibrationCam #include "MCalibrationCam.h" -#endif - -#ifndef ROOT_TArrayI -#include #endif Index: trunk/MagicSoft/Mars/mhcalib/HCalibLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/HCalibLinkDef.h (revision 6328) +++ trunk/MagicSoft/Mars/mhcalib/HCalibLinkDef.h (revision 6329) @@ -16,4 +16,5 @@ #pragma link C++ class MHCalibrationHiLoCam+; #pragma link C++ class MHPedestalCam+; +#pragma link C++ class MHPedestalPix+; #pragma link C++ class MHCalibrationTestTimeCam+; Index: trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestTimeCam.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestTimeCam.h (revision 6328) +++ trunk/MagicSoft/Mars/mhcalib/MHCalibrationTestTimeCam.h (revision 6329) @@ -14,5 +14,5 @@ static const Axis_t fgFirst; //! Default for fFirst (now set to: -0.5 ) static const Axis_t fgLast; //! Default for fLast (now set to: 29.5 ) - static const Float_t fgProbLimit; //! The default for fProbLimit (now set to: 0.000001) + static const Float_t fgProbLimit; //! Default for fProbLimit (now set to: 0.000001) static const TString gsHistName; //! Default Histogram names Index: trunk/MagicSoft/Mars/mhcalib/MHPedestalCam.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHPedestalCam.cc (revision 6328) +++ trunk/MagicSoft/Mars/mhcalib/MHPedestalCam.cc (revision 6329) @@ -16,5 +16,4 @@ ! ! -! Author(s): Thomas Bretz 12/2000 ! Author(s): Markus Gaug 02/2004 ! @@ -91,6 +90,4 @@ #include "MHPedestalPix.h" -#include "MBadPixelsCam.h" - #include "MLog.h" #include "MLogManip.h" @@ -103,5 +100,10 @@ #include "MPedestalCam.h" #include "MPedestalPix.h" + +#include "MCalibrationIntensityCam.h" #include "MCalibrationPix.h" + +#include "MBadPixelsIntensityCam.h" +#include "MBadPixelsCam.h" #include "MGeomCam.h" @@ -118,11 +120,12 @@ using namespace std; -const Int_t MHPedestalCam::fgNbins = 100; -const Axis_t MHPedestalCam::fgFirst = -50.; -const Axis_t MHPedestalCam::fgLast = 50.; +const Int_t MHPedestalCam::fgNbins = 100; +const Axis_t MHPedestalCam::fgFirst = -49.; +const Axis_t MHPedestalCam::fgLast = 151.; const TString MHPedestalCam::gsHistName = "Pedestal"; const TString MHPedestalCam::gsHistTitle = "Pedestal"; const TString MHPedestalCam::gsHistXTitle = "Charge [FADC slices]"; const TString MHPedestalCam::gsHistYTitle = "Nr. events"; +const TString MHPedestalCam::fgNamePedestalCam = "MPedestalCam"; // -------------------------------------------------------------------------- // @@ -132,5 +135,5 @@ // - fExtractHiGainSlices to 0. // - the event frequency to 1200 Hz. -// - the fRenorm flag to kFALSE +// - the fRenormflag to kFALSE // // - fNbins to fgNbins @@ -144,5 +147,6 @@ // MHPedestalCam::MHPedestalCam(const char *name, const char *title) - : fNumEvents(0), fExtractHiGainSlices(0.) + : fNamePedestalCamOut(fgNamePedestalCam), fNumEvents(0), + fExtractHiGainSlices(0.), fPedestalsOut(NULL) { @@ -151,5 +155,5 @@ SetPulserFrequency(1200); - SetRenorm(kTRUE); + SetRenorm(kFALSE); SetLoGain(kFALSE); @@ -162,4 +166,5 @@ SetHistXTitle(gsHistXTitle.Data()); SetHistYTitle(gsHistYTitle.Data()); + } @@ -170,11 +175,11 @@ // Resets: // - fSum -// - fSum2 +// - fSumSquare // - fAreaSum -// - fAreaSum2 -// - fAreaVal +// - fAreaSumSquare +// - fAreaNum // - fSectorSum -// - fSectorSum2 -// - fSectorVal +// - fSectorSumSquare +// - fSectorNum // void MHPedestalCam::ResetHists() @@ -190,16 +195,15 @@ // Reset contents of arrays. fSum.Reset(); - fSum2.Reset(); + fSumSquare.Reset(); fAreaSum. Reset(); - fAreaSum2.Reset(); - fAreaVal.Reset(); + fAreaSumSquare.Reset(); + fAreaNum.Reset(); fSectorSum. Reset(); - fSectorSum2.Reset(); - fSectorVal.Reset(); - } -} - + fSectorSumSquare.Reset(); + fSectorNum.Reset(); + } +} // -------------------------------------------------------------------------- // @@ -269,15 +273,19 @@ MExtractedSignalCam *signal = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam"); - if (!signal) - { - *fLog << err << "Cannot find MExtractedSignalCam... abort." << endl; + if (!signal && fRenorm) + { + *fLog << err << "Cannot find MExtractedSignalCam, but re-normalization switched on..." + << " Cannot find number of used slices ...abort." << endl; return kFALSE; } - fPedestals = (MPedestalCam*)pList->FindObject("MPedestalCam"); - - if (!fPedestals) - *fLog << warn << "Cannot find MPedestalCam ..." << endl; + fPedestalsOut = (MPedestalCam*)pList->FindObject(AddSerialNumber(fNamePedestalCamOut),"MPedestalCam"); + + if (!fPedestalsOut) + { + *fLog << err << "Cannot find nor create outgoing MPedestalCam ..." << endl; + return kFALSE; + } const Int_t npixels = fGeom->GetNumPixels(); @@ -295,45 +303,125 @@ } + InitHiGainArrays(npixels,nareas,nsectors); + + if (fSum.GetSize()==0) + { + fSum. Set(npixels); + fSumSquare. Set(npixels); + fAreaSum. Set(nareas); + fAreaSumSquare. Set(nareas); + fAreaNum. Set(nareas); + fSectorSum. Set(nsectors); + fSectorSumSquare.Set(nsectors); + fSectorNum.Set(nsectors); + } + + fNumEvents = 0; + + if (!signal) + return kTRUE; + + Float_t sliceshi = signal->GetNumUsedHiGainFADCSlices(); - + if (sliceshi == 0.) { *fLog << err << "Number of used signal slices in MExtractedSignalCam is zero ... abort." - << endl; + << endl; return kFALSE; } - + if (fExtractHiGainSlices != 0. && sliceshi != fExtractHiGainSlices ) { *fLog << err << "Number of used High Gain signal slices changed in MExtractedSignalCam ... abort." - << endl; + << endl; return kFALSE; } - + fExtractHiGainSlices = sliceshi; - *fLog << endl; *fLog << inf << GetDescriptor() - << ": Number of found High Gain signal slices: " << fExtractHiGainSlices << endl; - - InitHiGainArrays(npixels,nareas,nsectors); - - if (fSum.GetSize()==0) - { - fSum. Set(npixels); - fSum2. Set(npixels); - fAreaSum. Set(nareas); - fAreaSum2. Set(nareas); - fAreaVal. Set(nareas); - fSectorSum. Set(nsectors); - fSectorSum2.Set(nsectors); - fSectorVal.Set(nsectors); - } - - fNumEvents = 0; + << ": Number of found High Gain signal slices: " << fExtractHiGainSlices << endl; return kTRUE; } +// -------------------------------------------------------------------------- +// +// Initializes the High Gain Arrays: +// +// - For every entry in the expanded arrays: +// * Initialize an MHCalibrationChargePix +// * Set Binning from fNbins, fFirst and fLast +// * Set Binning of Abs Times histogram from fAbsNbins, fAbsFirst and fAbsLast +// * Set Histgram names and titles from fHistName and fHistTitle +// * Set Abs Times Histgram names and titles from fAbsHistName and fAbsHistTitle +// * Set X-axis and Y-axis titles from fHistXTitle and fHistYTitle +// * Set X-axis and Y-axis titles of Abs Times Histogram from fAbsHistXTitle and fAbsHistYTitle +// * Call InitHists +// +// +void MHPedestalCam::InitHiGainArrays(const Int_t npixels, const Int_t nareas, const Int_t nsectors) +{ + + if (fHiGainArray->GetSize()==0) + { + for (Int_t i=0; iAddAt(new MHPedestalPix(Form("%sHiGainPix%04d",fHistName.Data(),i), + Form("%s High Gain Pixel%04d",fHistTitle.Data(),i)),i); + + MHPedestalPix &pix = (MHPedestalPix&)(*this)[i]; + + pix.SetNbins(fNbins); + pix.SetFirst(fFirst); + pix.SetLast (fLast); + + pix.SetProbLimit(fProbLimit); + + MBadPixelsPix &bad = fIntensBad ? (*fIntensBad)[i] : (*fBadPixels)[i]; + InitHists(pix,bad,i); + } + } + + if (!IsAverageing()) + return; + + if (fAverageHiGainAreas->GetSize()==0) + { + for (Int_t j=0; jAddAt(new MHPedestalPix(Form("%sHiGainArea%d",fHistName.Data(),j), + Form("%s High Gain Area Idx %d",fHistTitle.Data(),j)),j); + + MHPedestalPix &pix = (MHPedestalPix&)GetAverageHiGainArea(j); + + pix.SetNbins(fNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas)); + pix.SetFirst(fFirst); + pix.SetLast (fLast); + + InitHists(pix,fIntensCam ? fIntensCam->GetAverageBadArea(j) : fCam->GetAverageBadArea(j),j); + + } + } + + if (fAverageHiGainSectors->GetSize()==0) + { + for (Int_t j=0; jAddAt(new MHPedestalPix(Form("%sHiGainSector%02d",fHistName.Data(),j), + Form("%s High Gain Sector %02d",fHistTitle.Data(),j)),j); + + MHPedestalPix &pix = (MHPedestalPix&)GetAverageHiGainSector(j); + + pix.SetNbins(fNbins*(Int_t)TMath::Sqrt((Float_t)npixels/nareas)); + pix.SetFirst(fFirst); + pix.SetLast (fLast); + + InitHists(pix,fIntensCam ? fIntensCam->GetAverageBadSector(j) : fCam->GetAverageBadSector(j),j); + + } + } +} // ------------------------------------------------------------------------------- @@ -352,8 +440,8 @@ { - MExtractedSignalCam *signal = (MExtractedSignalCam*)par; - if (!signal) - { - *fLog << err << "No argument in MExtractedSignalCam::Fill... abort." << endl; + MPedestalCam *pedestal = (MPedestalCam*)par; + if (!pedestal) + { + *fLog << err << "No argument in MHPedestalCam::Fill... abort." << endl; return kFALSE; } @@ -376,25 +464,25 @@ continue; - const MExtractedSignalPix &pix = (*signal)[i]; - - const Float_t pedhi = pix.GetExtractedSignalHiGain(); + const MPedestalPix &ped = (*pedestal)[i]; + + const Float_t pedes = ped.GetPedestal(); const Int_t aidx = (*fGeom)[i].GetAidx(); const Int_t sector = (*fGeom)[i].GetSector(); - histhi.FillHistAndArray(pedhi) ; - - fSum[i] += pedhi; - fAreaSum[aidx] += pedhi; - fSectorSum[sector] += pedhi; - - const Float_t sqrsum = pedhi*pedhi; - fSum2[i] += sqrsum; - fAreaSum2[aidx] += sqrsum; - fSectorSum2[sector] += sqrsum; - - sumareahi [aidx] += pedhi; + histhi.FillHistAndArray(pedes) ; + + fSum[i] += pedes; + fAreaSum[aidx] += pedes; + fSectorSum[sector] += pedes; + + const Float_t sqrsum = pedes*pedes; + fSumSquare[i] += sqrsum; + fAreaSumSquare[aidx] += sqrsum; + fSectorSumSquare[sector] += sqrsum; + + sumareahi [aidx] += pedes; numareahi [aidx] ++; - sumsectorhi[sector] += pedhi; + sumsectorhi[sector] += pedes; numsectorhi[sector] ++; @@ -449,5 +537,5 @@ { - MHCalibrationPix &hist = (*this)[i]; + MHPedestalPix &hist = (MHPedestalPix&)(*this)[i]; MCalibrationPix &pix = (*fCam)[i]; @@ -456,5 +544,5 @@ // pix.SetLoGainMean ( fSum[i] / fNumEvents ); - const Double_t diff = fSum2[i] - fSum[i]*fSum[i]/fNumEvents; + const Double_t diff = fSumSquare[i] - fSum[i]*fSum[i]/fNumEvents; pix.SetLoGainSigma ( diff > 0. ? TMath::Sqrt( diff / (fNumEvents-1) ) : 0.); pix.SetLoGainMeanVar ( pix.GetLoGainSigma() * pix.GetLoGainSigma() / fNumEvents ); @@ -467,4 +555,7 @@ // 2) Fit the Hi Gain histograms with a Gaussian // + // if (i%100) + // hist.FitTripleGaus(); + // else hist.FitGaus(); // @@ -529,10 +620,10 @@ // 6) Store calculated means and variances in the low-gain slices // - const ULong_t aevts = fNumEvents * fAreaVal[j]; + const ULong_t aevts = fNumEvents * fAreaNum[j]; if (aevts <= 1) continue; pix.SetLoGainMean ( fAreaSum[j] / aevts ); - const Double_t diff = fAreaSum2[j] - fAreaSum[j]*fAreaSum[j]/aevts ; + const Double_t diff = fAreaSumSquare[j] - fAreaSum[j]*fAreaSum[j]/aevts ; pix.SetLoGainSigma( diff > 0. ? TMath::Sqrt( diff / (aevts-1) ) : 0.); } @@ -573,10 +664,10 @@ // 6) Store calculated means and variances in the low-gain slices // - const ULong_t sevts = fNumEvents * fSectorVal[j]; + const ULong_t sevts = fNumEvents * fSectorNum[j]; if (sevts <= 1) continue; pix.SetLoGainMean ( fSectorSum[j] / sevts ); - const Double_t diff = fSectorSum2[j] - fSectorSum[j]*fSectorSum[j]/sevts ; + const Double_t diff = fSectorSumSquare[j] - fSectorSum[j]*fSectorSum[j]/sevts ; pix.SetLoGainSigma( diff > 0. ? TMath::Sqrt( diff / (sevts-1) ) : 0.); } @@ -694,8 +785,8 @@ return kFALSE; - const Float_t ped = (*fPedestals)[idx].GetPedestal(); - const Float_t rms = (*fPedestals)[idx].GetPedestalRms(); - - const Float_t entsqr = TMath::Sqrt((Float_t)fPedestals->GetTotalEntries()); + const Float_t ped = (*fPedestalsOut)[idx].GetPedestal(); + const Float_t rms = (*fPedestalsOut)[idx].GetPedestalRms(); + + const Float_t entsqr = TMath::Sqrt((Float_t)fPedestalsOut->GetTotalEntries()); const Float_t pederr = rms/entsqr; Index: trunk/MagicSoft/Mars/mhcalib/MHPedestalCam.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHPedestalCam.h (revision 6328) +++ trunk/MagicSoft/Mars/mhcalib/MHPedestalCam.h (revision 6329) @@ -17,12 +17,16 @@ private: - static const Int_t fgNbins; //! Default number of bins (now set to: 100 ) - static const Axis_t fgFirst; //! Default lower histogram limit (now set to: -50. ) - static const Axis_t fgLast; //! Default upper histogram limit (now set to: 50. ) - - static const TString gsHistName; //! Default Histogram names - static const TString gsHistTitle; //! Default Histogram titles - static const TString gsHistXTitle; //! Default Histogram x-axis titles - static const TString gsHistYTitle; //! Default Histogram y-axis titles + static const Int_t fgNbins; //! Default number of bins (now set to: 100 ) + static const Axis_t fgFirst; //! Default lower histogram limit (now set to: -50. ) + static const Axis_t fgLast; //! Default upper histogram limit (now set to: 50. ) + + static const TString gsHistName; //! Default Histogram names + static const TString gsHistTitle; //! Default Histogram titles + static const TString gsHistXTitle; //! Default Histogram x-axis titles + static const TString gsHistYTitle; //! Default Histogram y-axis titles + + static const TString fgNamePedestalCam; //! "MPedestalCam" + + TString fNamePedestalCamOut; // Name of the 'MPedestalCam' container with results ULong_t fNumEvents; // Number of processed events @@ -31,20 +35,22 @@ Float_t fExtractLoGainSlices; // Number of FADC slices used for low gain signal extraction - MPedestalCam *fPedestals; //! Pedestal Cam filled by MPedCalcPedRun - + MPedestalCam *fPedestalsOut; //! Pedestal Cam with results + Bool_t fRenorm; // Flag if the results will be re-normalized - MArrayD fSum; //! sum of values - MArrayD fSum2; //! sum of squared values - MArrayD fAreaSum; //! averaged sum of values per area idx - MArrayD fAreaSum2; //! averaged sum of squared values per area idx - MArrayI fAreaVal; //! number of valid pixel with area idx - MArrayD fSectorSum; //! averaged sum of values per sector - MArrayD fSectorSum2; //! averaged sum of squared values per sector - MArrayI fSectorVal; //! number of valid pixel with sector idx + MArrayD fSum; //! sum of values + MArrayD fSumSquare; //! sum of squared values + MArrayD fAreaSum; //! averaged sum of values per area idx + MArrayD fAreaSumSquare; //! averaged sum of squared values per area idx + MArrayI fAreaNum; //! number of valid pixel with area idx + MArrayI fSectorNum; //! number of valid pixel with sector idx + MArrayD fSectorSum; //! averaged sum of values per sector + MArrayD fSectorSumSquare; //! averaged sum of squared values per sector Bool_t ReInitHists(MParList *pList); Bool_t FillHists(const MParContainer *par, const Stat_t w=1); Bool_t FinalizeHists(); + + void InitHiGainArrays( const Int_t npix, const Int_t nareas, const Int_t nsectors ); void RenormResults(); @@ -55,5 +61,5 @@ MHPedestalCam(const char *name=NULL, const char *title=NULL); - // Clone + // Clone TObject *Clone(const char *name="") const; @@ -62,5 +68,7 @@ void ResetHists(); - void SetRenorm( const Bool_t b=kTRUE ) { fRenorm = b; } + + void SetNamePedestalCamOut(const char *name) { fNamePedestalCamOut = name; } + void SetRenorm (const Bool_t b=kTRUE ) { fRenorm = b; } ClassDef(MHPedestalCam, 1) // Histogram class for Charge Camera Pedestals Index: trunk/MagicSoft/Mars/mhcalib/MHPedestalPix.cc =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHPedestalPix.cc (revision 6329) +++ trunk/MagicSoft/Mars/mhcalib/MHPedestalPix.cc (revision 6329) @@ -0,0 +1,236 @@ +/* ======================================================================== *\ +! +! * +! * 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/2005 +! +! Copyright: MAGIC Software Development, 2000-2005 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHPedestalPix +// +// A base class for events which are believed to follow a Gaussian distribution +// with time, e.g. calibration events, observables containing white noise, ... +// +// MHPedestalPix derives from MHGausEvents, thus all features of +// MHGausEvents can be used by a class deriving from MHPedestalPix +// +// As an additional feature to MHGausEvents, this class offers to skip the fitting +// to set mean, sigma and its errors directly from the histograms with the function +// BypassFit() +// +// See also: MHGausEvents +// +////////////////////////////////////////////////////////////////////////////// +#include "MHPedestalPix.h" + +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MHPedestalPix); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MHPedestalPix::MHPedestalPix(const char *name, const char *title) +{ + + fName = name ? name : "MHPedestalPix"; + fTitle = title ? title : "Pedestal histogram events"; + +} + +// ------------------------------------------------------------------------------- +// +// Fits the histogram to a double Gauss. +// +// +Bool_t MHPedestalPix::FitDoubleGaus(const Double_t xmin, const Double_t xmax, Option_t *option) +{ + + if (IsGausFitOK()) + return kTRUE; + + StripZeros(&fHGausHist,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 sigma_guess = fHGausHist.GetRMS(); + const Double_t area_guess = entries/TMath::Sqrt(TMath::TwoPi())/sigma_guess; + + fFGausFit = new TF1("GausFit","gaus(0)+gaus(3)",rmin,rmax); + + if (!fFGausFit) + { + *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit " + << "in: " << fName << endl; + return kFALSE; + } + + // + // For the fits, we have to take special care since ROOT + // has stored the function pointer in a global list which + // lead to removing the object twice. We have to take out + // the following functions of the global list of functions + // as well: + // + gROOT->GetListOfFunctions()->Remove(fFGausFit); + + fFGausFit->SetParameters(area_guess/2.,0.,sigma_guess/2.,area_guess/2.,25.,sigma_guess/2.); + fFGausFit->SetParNames("Area_{0}","#mu_{0}","#sigma_{0}","Area_{1}","#mu_{1}","#sigma_{1}"); + fFGausFit->SetParLimits(0,0.,area_guess*5.); + fFGausFit->SetParLimits(1,rmin,0.); + fFGausFit->SetParLimits(2,0.,rmax-rmin); + fFGausFit->SetParLimits(3,0.,area_guess*10.); + fFGausFit->SetParLimits(4,0.,rmax/2.); + fFGausFit->SetParLimits(5,0.,rmax-rmin); + fFGausFit->SetRange(rmin,rmax); + + fHGausHist.Fit(fFGausFit,option); + + SetMean (fFGausFit->GetParameter(4)-fFGausFit->GetParameter(1)); + SetSigma (TMath::Sqrt(fFGausFit->GetParameter(5)*fFGausFit->GetParameter(5) + +fFGausFit->GetParameter(2)*fFGausFit->GetParameter(2))); + SetMeanErr (TMath::Sqrt(fFGausFit->GetParError(4)*fFGausFit->GetParError(4) + +fFGausFit->GetParError(1)*fFGausFit->GetParError(1))); + SetSigmaErr (TMath::Sqrt(fFGausFit->GetParError(5)*fFGausFit->GetParError(5) + +fFGausFit->GetParError(2)*fFGausFit->GetParError(2))); + SetProb (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(GetMean()) + || TMath::IsNaN(GetMeanErr()) + || TMath::IsNaN(GetProb()) + || TMath::IsNaN(GetSigma()) + || TMath::IsNaN(GetSigmaErr()) + || fProb < fProbLimit ) + return kFALSE; + + SetGausFitOK(kTRUE); + return kTRUE; +} + + +// ------------------------------------------------------------------------------- +// +// Fits the histogram to a triple Gauss. +// +Bool_t MHPedestalPix::FitTripleGaus(const Double_t xmin, const Double_t xmax, Option_t *option) +{ + + if (IsGausFitOK()) + return kTRUE; + + StripZeros(&fHGausHist,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 sigma_guess = fHGausHist.GetRMS(); + const Double_t area_guess = entries/TMath::Sqrt(TMath::TwoPi())/sigma_guess; + + fFGausFit = new TF1("GausFit","gaus(0)+gaus(3)+gaus(6)",rmin,rmax); + + if (!fFGausFit) + { + *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit " + << "in: " << fName << endl; + return kFALSE; + } + + // + // For the fits, we have to take special care since ROOT + // has stored the function pointer in a global list which + // lead to removing the object twice. We have to take out + // the following functions of the global list of functions + // as well: + // + gROOT->GetListOfFunctions()->Remove(fFGausFit); + + fFGausFit->SetParameters(10.,-4.0,1.5,70.,1.5,6.,5.,7.,7.); + fFGausFit->SetParNames("Area_{0}","#mu_{0}","#sigma_{0}","Area_{1}","#mu_{1}","#sigma_{1}","Area_{2}","#mu_{2}","#sigma_{2}"); + fFGausFit->SetParLimits(0,0.,area_guess*2.5); + fFGausFit->SetParLimits(1,-9.0,-2.2); + fFGausFit->SetParLimits(2,-1.0,15.); + fFGausFit->SetParLimits(3,0.,area_guess*10.); + fFGausFit->SetParLimits(4,-4.5,2.); + fFGausFit->SetParLimits(5,0.,(rmax-rmin)/3.); + fFGausFit->SetParLimits(6,0.,area_guess*5.); + fFGausFit->SetParLimits(7,6.,20.); + fFGausFit->SetParLimits(8,5.,40.); + fFGausFit->SetRange(rmin,rmax); + + fHGausHist.Fit(fFGausFit,option); + + SetMean (fFGausFit->GetParameter(4)-fFGausFit->GetParameter(1)); + SetSigma (TMath::Sqrt(fFGausFit->GetParameter(5)*fFGausFit->GetParameter(5) + +fFGausFit->GetParameter(2)*fFGausFit->GetParameter(2))); + SetMeanErr (TMath::Sqrt(fFGausFit->GetParError(4)*fFGausFit->GetParError(4) + +fFGausFit->GetParError(1)*fFGausFit->GetParError(1))); + SetSigmaErr (TMath::Sqrt(fFGausFit->GetParError(5)*fFGausFit->GetParError(5) + +fFGausFit->GetParError(2)*fFGausFit->GetParError(2))); + SetProb (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(GetMean()) + || TMath::IsNaN(GetMeanErr()) + || TMath::IsNaN(GetProb()) + || TMath::IsNaN(GetSigma()) + || TMath::IsNaN(GetSigmaErr()) + || fProb < fProbLimit ) + return kFALSE; + + SetGausFitOK(kTRUE); + return kTRUE; +} + + Index: trunk/MagicSoft/Mars/mhcalib/MHPedestalPix.h =================================================================== --- trunk/MagicSoft/Mars/mhcalib/MHPedestalPix.h (revision 6329) +++ trunk/MagicSoft/Mars/mhcalib/MHPedestalPix.h (revision 6329) @@ -0,0 +1,25 @@ +#ifndef MARS_MHPedestalPix +#define MARS_MHPedestalPix + +#ifndef MARS_MHCalibrationPix +#include "MHCalibrationPix.h" +#endif + +class MHPedestalPix : public MHCalibrationPix +{ +public: + + MHPedestalPix(const char* name=NULL, const char* title=NULL); + + // Fits + Bool_t FitDoubleGaus( const Double_t xmin=0., + const Double_t xmax=0., + Option_t *option="RQ0"); + Bool_t FitTripleGaus( const Double_t xmin=0., + const Double_t xmax=0., + Option_t *option="RQ0"); + + ClassDef(MHPedestalPix, 1) // Base class for histogrammed pedestal events +}; + +#endif Index: trunk/MagicSoft/Mars/mhcalib/Makefile =================================================================== --- trunk/MagicSoft/Mars/mhcalib/Makefile (revision 6328) +++ trunk/MagicSoft/Mars/mhcalib/Makefile (revision 6329) @@ -43,4 +43,5 @@ MHCalibrationTestCam.cc \ MHPedestalCam.cc \ + MHPedestalPix.cc \ MHCalibrationTestTimeCam.cc \ MHGausEvents.cc