Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 3158) +++ trunk/MagicSoft/Mars/Changelog (revision 3159) @@ -9,4 +9,7 @@ * manalysis/MHPedestalPix.[h,cc] - new version of MHPedestalPixel, deriving from MHGausEvents. + + * manalysis/MHPedestalPixel.[h,cc] + - old version removed, since obsolete (and erroneous) * manalysis/MHPedestalCam.[h,cc] Index: trunk/MagicSoft/Mars/manalysis/MHPedestalPixel.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHPedestalPixel.cc (revision 3158) +++ (revision ) @@ -1,392 +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 02/2004 -! -! Copyright: MAGIC Software Development, 2000-2004 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// -// MHPedestalPixel -// -// Performs all the necessary fits to extract the mean number of photons -// out of the derived light flux -// -////////////////////////////////////////////////////////////////////////////// -#include "MHPedestalPixel.h" - -#include -#include - -#include -#include - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHPedestalPixel); - -using namespace std; -// Square root of 2pi: -const Float_t gkSq2Pi = 2.506628274631; -const Float_t gkProbLimit = 0.01; -const Int_t MHPedestalPixel::gkChargeNbins = 450 ; -const Axis_t MHPedestalPixel::gkChargeFirst = -0.5; -const Axis_t MHPedestalPixel::gkChargeLast = 449.5; -const Int_t MHPedestalPixel::gkChargevsNbins = 1000; -const Axis_t MHPedestalPixel::gkChargevsNFirst = -0.5; -const Axis_t MHPedestalPixel::gkChargevsNLast = 999.5; - -// -------------------------------------------------------------------------- -// -// Default Constructor. -// -MHPedestalPixel::MHPedestalPixel() - : fPixId(-1), - fGausFit(NULL) -{ - - // Create a large number of bins, later we will rebin - fHPedestalCharge = new TH1F("HPedestalCharge","Distribution of Summed FADC Pedestal Slices Pixel ", - gkChargeNbins,gkChargeFirst,gkChargeLast); - fHPedestalCharge->SetXTitle("Sum FADC Slices"); - fHPedestalCharge->SetYTitle("Nr. of events"); - fHPedestalCharge->Sumw2(); - - // We define a reasonable number and later enlarge it if necessary - fHPedestalChargevsN = new TH1I("HChargevsN","Sum of Charges vs. Event Number Pixel ", - gkChargevsNbins,gkChargevsNFirst,gkChargevsNLast); - fHPedestalChargevsN->SetXTitle("Event Nr."); - fHPedestalChargevsN->SetYTitle("Sum of FADC slices"); - - fHPedestalCharge->SetDirectory(NULL); - fHPedestalChargevsN->SetDirectory(NULL); - - Clear(); -} - -MHPedestalPixel::~MHPedestalPixel() -{ - - delete fHPedestalCharge; - delete fHPedestalChargevsN; - - if (fGausFit) - delete fGausFit; -} - - -void MHPedestalPixel::Clear(Option_t *o) -{ - - fTotalEntries = 0; - - fChargeMean = -1.; - fChargeMeanErr = -1.; - fChargeSigma = -1; - fChargeSigmaErr = -1; - - fChargeChisquare = -1.; - fChargeProb = -1.; - fChargeNdf = -1; - - fChargeFirst = gkChargeFirst; - fChargeLast = gkChargeLast; - - if (fGausFit) - delete fGausFit; - - CLRBIT(fFlags,kFitted); - CLRBIT(fFlags,kFitOK); - CLRBIT(fFlags,kOscillating); - - return; -} - - -void MHPedestalPixel::Reset() -{ - - Clear(); - - fHPedestalCharge->Reset(); - fHPedestalChargevsN->Reset(); - -} - - - -Bool_t MHPedestalPixel::IsEmpty() const -{ - return !fHPedestalCharge->GetEntries(); -} - -Bool_t MHPedestalPixel::IsFitOK() const -{ - return TESTBIT(fFlags,kFitOK); -} - -Bool_t MHPedestalPixel::FillCharge(Float_t q) -{ - return (fHPedestalCharge->Fill(q) > -1); -} - -Bool_t MHPedestalPixel::FillChargevsN(Float_t q) -{ - - fTotalEntries++; - return (fHPedestalChargevsN->Fill(fTotalEntries,q) > -1); -} - -void MHPedestalPixel::ChangeHistId(Int_t id) -{ - - fPixId = id; - - fHPedestalCharge->SetName(Form("%s%d", fHPedestalCharge->GetName(), id)); - fHPedestalChargevsN->SetName(Form("%s%d", fHPedestalChargevsN->GetName(), id)); - fHPedestalCharge->SetTitle(Form("%s%d", fHPedestalCharge->GetTitle(), id)); - fHPedestalChargevsN->SetTitle(Form("%s%d", fHPedestalChargevsN->GetTitle(), id)); - -} - -TObject *MHPedestalPixel::DrawClone(Option_t *option) const -{ - - gROOT->SetSelectedPad(NULL); - - MHPedestalPixel *newobj = (MHPedestalPixel*)Clone(); - - if (!newobj) - return 0; - newobj->SetBit(kCanDelete); - - - if (strlen(option)) - newobj->Draw(option); - else - newobj->Draw(GetDrawOption()); - - return newobj; -} - - -// ------------------------------------------------------------------------- -// -// Draw the histogram -// -void MHPedestalPixel::Draw(Option_t *opt) -{ - - gStyle->SetOptStat(111111); - gStyle->SetOptFit(); - - gROOT->SetSelectedPad(NULL); - - TCanvas *c = MH::MakeDefCanvas(this,600,900); - - // c->Divide(1,2); - - c->cd(1); - gPad->SetBorderMode(0); - gPad->SetTicks(); - - if (fHPedestalCharge->Integral() > 0) - gPad->SetLogy(); - - fHPedestalCharge->Draw(opt); - - c->Modified(); - c->Update(); - - if (fGausFit) - { - fGausFit->SetLineColor(IsFitOK() ? kGreen : kRed); - fGausFit->Draw("same"); - } - - c->Modified(); - c->Update(); - - /* - c->cd(2); - gPad->SetTicks(); - - fHPedestalChargevsN->Draw(opt); - */ - - c->Modified(); - c->Update(); - - return; -} - - -// -------------------------------------------------------------------------- -// -// 1) Return if the charge distribution is already succesfully fitted -// or if the histogram is empty -// 2) Cut the histograms empty edges -// 3) Fit the histograms with a Gaussian -// 4) In case of failure print out the fit results -// 5) Retrieve the results and store them in this class -// -Bool_t MHPedestalPixel::FitCharge(Option_t *option) -{ - - if (fGausFit) - return kFALSE; - - // - // Get the fitting ranges - // - Axis_t rmin = fChargeFirst; - Axis_t rmax = fChargeLast; - - // - // First guesses for the fit (should be as close to reality as possible, - // otherwise the fit goes gaga because of high number of dimensions ... - // - const Stat_t entries = fHPedestalCharge->Integral(); - const Double_t mu_guess = fHPedestalCharge->GetBinCenter(fHPedestalCharge->GetMaximumBin()); - const Double_t sigma_guess = mu_guess/15.; - const Double_t area_guess = entries/gkSq2Pi/sigma_guess; - - fGausFit = new TF1(Form("%s%d","GausFit ",fPixId),"gaus",rmin,rmax); - - if (!fGausFit) - { - *fLog << warn << dbginf << "WARNING: Could not create fit function for Gauss fit" << endl; - return kFALSE; - } - - fGausFit->SetParameters(area_guess,mu_guess,sigma_guess); - fGausFit->SetParNames("Area","#mu","#sigma"); - fGausFit->SetParLimits(0,0.,entries); - fGausFit->SetParLimits(1,rmin,rmax); - fGausFit->SetParLimits(2,0.,rmax-rmin); - fGausFit->SetRange(rmin,rmax); - - fHPedestalCharge->Fit(fGausFit,option); - - // - // In order not to be affected by outliers, - // try once again with stricter borders - // - Axis_t rtry = fGausFit->GetParameter(1) - 2.5*fGausFit->GetParameter(2); - rmin = (rtry < rmin ? rmin : rtry); - rtry = fGausFit->GetParameter(1) + 2.5*fGausFit->GetParameter(2); - rmax = (rtry > rmax ? rmax : rtry); - fGausFit->SetRange(rmin,rmax); - - fHPedestalCharge->Fit(fGausFit,option); - - fChargeChisquare = fGausFit->GetChisquare(); - fChargeNdf = fGausFit->GetNDF(); - fChargeProb = fGausFit->GetProb(); - fChargeMean = fGausFit->GetParameter(1); - fChargeMeanErr = fGausFit->GetParError(1); - fChargeSigma = fGausFit->GetParameter(2); - fChargeSigmaErr = fGausFit->GetParError(2); - - SETBIT(fFlags,kFitted); - // - // The fit result is accepted under condition: - // The Results are not nan's - // The Probability is greater than gkProbLimit (default 0.001 == 99.9%) - // - if (TMath::IsNaN(fChargeMean) || - TMath::IsNaN(fChargeMeanErr) || - TMath::IsNaN(fChargeSigma) || - TMath::IsNaN(fChargeSigmaErr)) - { - Clear(); - // CLRBIT(fFlags,kFitOK); - return kFALSE; - } - - if ((fChargeProb < gkProbLimit) || (TMath::IsNaN(fChargeProb))) - { - CLRBIT(fFlags,kFitOK); - return kFALSE; - } - - SETBIT(fFlags,kFitOK); - return kTRUE; -} - -void MHPedestalPixel::Renorm(const Float_t nslices) -{ - - if (!TESTBIT(fFlags,kFitOK)) - return; - - Float_t sqslices = TMath::Sqrt(nslices); - - fChargeMean /= nslices; - // - // Mean error goes with PedestalRMS/Sqrt(entries) -> scale with slices - // - fChargeMeanErr /= nslices; - // - // Sigma goes like PedestalRMS -> scale with sqrt(slices) - // - fChargeSigma /= sqslices; - // - // Sigma error goes like PedestalRMS/2.(entries) -> scale with slices - // - fChargeSigmaErr /= nslices; - -} - -void MHPedestalPixel::CutAllEdges() -{ - - Int_t nbins = 15; - - MH::StripZeros(fHPedestalCharge,nbins); - - fChargeFirst = fHPedestalCharge->GetBinLowEdge(fHPedestalCharge->GetXaxis()->GetFirst()); - fChargeLast = fHPedestalCharge->GetBinLowEdge(fHPedestalCharge->GetXaxis()->GetLast()) - +fHPedestalCharge->GetBinWidth(0); - - MH::StripZeros(fHPedestalChargevsN,0); - -} - -void MHPedestalPixel::Print(const Option_t *o) const -{ - - *fLog << all << "Pedestal Fits Pixel: " << fPixId << endl; - - if (!TESTBIT(fFlags,kFitted)) - { - gLog << "Pedestal Histogram has not yet been fitted" << endl; - return; - } - - gLog << "Results of the Summed Charges Fit: " << endl; - gLog << "Chisquare: " << fChargeChisquare << endl; - gLog << "DoF: " << fChargeNdf << endl; - gLog << "Probability: " << fChargeProb << endl; - gLog << "Results of fit: " << (TESTBIT(fFlags,kFitOK) ? "Ok.": "Not OK!") << endl; - -} - Index: trunk/MagicSoft/Mars/manalysis/MHPedestalPixel.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHPedestalPixel.h (revision 3158) +++ (revision ) @@ -1,98 +1,0 @@ -#ifndef MARS_MHPedestalPixel -#define MARS_MHPedestalPixel - -#ifndef MARS_MH -#include "MH.h" -#endif - -class TH1F; -class TH1I; -class TF1; -class MHPedestalPixel : public MH -{ - -private: - - static const Int_t gkChargeNbins; - static const Int_t gkChargevsNbins; - static const Axis_t gkChargevsNFirst; - static const Axis_t gkChargevsNLast; - static const Axis_t gkChargeFirst; - static const Axis_t gkChargeLast; - - Int_t fPixId; // Pixel Nr - - Axis_t fChargeFirst; - Axis_t fChargeLast; - - const Int_t fChargevsNbins; - const Axis_t fChargevsNFirst; - const Axis_t fChargevsNLast; - - TH1F* fHPedestalCharge; //-> Summed FADC slices - TH1I* fHPedestalChargevsN; //-> Summed FADC slices vs Event nr. - TF1* fGausFit; // Fit the the Summed FADC slices - - Int_t fTotalEntries; // Number of entries - - Double_t fChargeChisquare; - Double_t fChargeProb; - Int_t fChargeNdf; - - Double_t fChargeMean; - Double_t fChargeMeanErr; - Double_t fChargeSigma; - Double_t fChargeSigmaErr; - - Byte_t fFlags; - - enum { kFitted, kFitOK, kOscillating }; - -public: - - MHPedestalPixel(); - ~MHPedestalPixel(); - - void Clear(Option_t *o=""); - void Reset(); - - void ChangeHistId(Int_t i); - - // Getters - const TH1F *GetHPedestalCharge() const { return fHPedestalCharge; } - - Double_t GetChargeMean() const { return fChargeMean; } - Double_t GetChargeMeanErr() const { return fChargeMeanErr; } - Double_t GetChargeSigma() const { return fChargeSigma; } - Double_t GetChargeSigmaErr() const { return fChargeSigmaErr; } - Double_t GetChargeChiSquare() const { return fChargeChisquare; } - Double_t GetChargeProb() const { return fChargeProb; } - Int_t GetChargeNdf() const { return fChargeNdf; } - - Int_t GetTotalEntries() const { return fTotalEntries; } - - Bool_t IsFitOK() const; - Bool_t IsEmpty() const; - - // Fill histos - Bool_t FillCharge(Float_t q); - Bool_t FillChargevsN(Float_t q); - - // Fits - Bool_t FitCharge(Option_t *option="RQ0"); - void Renorm(const Float_t nslices); - - // Draws - void Draw(Option_t *option=""); - TObject *DrawClone(Option_t *option="") const; - - // Prints - void Print(const Option_t *o="") const; - - // Others - void CutAllEdges(); - - ClassDef(MHPedestalPixel, 1) // Histograms for each calibrated pixel -}; - -#endif Index: trunk/MagicSoft/Mars/manalysis/MPedestalPix.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalPix.h (revision 3158) +++ trunk/MagicSoft/Mars/manalysis/MPedestalPix.h (revision 3159) @@ -6,9 +6,4 @@ #endif -#ifndef MARS_MHPedestalPixel -#include "MHPedestalPixel.h" -#endif - -class MHPedestalPixel; class MPedestalPix : public MParContainer { @@ -17,5 +12,4 @@ Float_t fPedestal; // mean value of pedestal (PMT offset) Float_t fPedestalRms; // root mean square / sigma / standard deviation of pedestal - public: