Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 3010) +++ trunk/MagicSoft/Mars/Changelog (revision 3011) @@ -4,4 +4,11 @@ -*-*- END OF LINE -*-*- + 2004/02/03: Markus Gaug + + * manalysis/MExtractedSignalPix.[h,cc] + - added Thomas B. modified version of MHPedestalPixel. Later will + remove MHPedestalPixel + + 2004/02/03: Thomas Bretz Index: trunk/MagicSoft/Mars/manalysis/MHExtractedSignalPix.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHExtractedSignalPix.cc (revision 3011) +++ trunk/MagicSoft/Mars/manalysis/MHExtractedSignalPix.cc (revision 3011) @@ -0,0 +1,359 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! Thomas Bretz 02/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHExtractedSignalPix +// +////////////////////////////////////////////////////////////////////////////// +#include "MHExtractedSignalPix.h" + +#include +#include +#include + +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MH.h" + +ClassImp(MHExtractedSignalPix); + +using namespace std; + +const Float_t gkSq2Pi = TMath::Sqrt(TMath::TwoPi()); +const Float_t gkProbLimit = 0.01; +const Int_t MHExtractedSignalPix::gkChargeNbins = 500 ; +const Int_t MHExtractedSignalPix::gkChargevsNbins = 1000; +const Axis_t MHExtractedSignalPix::gkChargevsNFirst = -0.5; +const Axis_t MHExtractedSignalPix::gkChargevsNLast = 999.5; +const Axis_t MHExtractedSignalPix::gkChargeFirst = -0.5; +const Axis_t MHExtractedSignalPix::gkChargeLast = 99.5; + +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MHExtractedSignalPix::MHExtractedSignalPix() + : fPixId(-1), + fGausFit(NULL) +{ + + //fName = name ? name : "MHExtractedSignalPix"; + //fTitle = title ? title : "Fill the accumulated charges and times of all events and perform fits"; + + // 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(); +} + +MHExtractedSignalPix::~MHExtractedSignalPix() +{ + + delete fHPedestalCharge; + delete fHPedestalChargevsN; + + if (fGausFit) + delete fGausFit; +} + + +void MHExtractedSignalPix::Clear(Option_t *o) +{ + + fTotalEntries = 0; + + fChargeMean = -1.; + fChargeMeanErr = -1.; + fChargeSigma = -1; + fChargeSigmaErr = -1; + + fChargeChisquare = -1.; + fChargeProb = -1.; + fChargeNdf = -1; + + fChargeFirst = -0.5; + fChargeLast = 99.5; + + if (fGausFit) + delete fGausFit; + + CLRBIT(fFlags,kFitted); + CLRBIT(fFlags,kFitOK); + CLRBIT(fFlags,kOscillating); + +} + + +void MHExtractedSignalPix::Reset() +{ + Clear(); + + fHPedestalCharge->Reset(); + fHPedestalChargevsN->Reset(); +} + + +Bool_t MHExtractedSignalPix::IsEmpty() const +{ + return !fHPedestalCharge->GetEntries(); +} + +Bool_t MHExtractedSignalPix::IsFitted() const +{ + return TESTBIT(fFlags,kFitted); +} + +Bool_t MHExtractedSignalPix::IsFitOK() const +{ + return TESTBIT(fFlags,kFitOK); +} + + +Bool_t MHExtractedSignalPix::IsOscillating() const +{ + return TESTBIT(fFlags,kOscillating); +} + + +Bool_t MHExtractedSignalPix::FillCharge(Float_t q) +{ + return (fHPedestalCharge->Fill(q) > -1); +} + +Bool_t MHExtractedSignalPix::FillChargevsN(Float_t q) +{ + fTotalEntries++; + return (fHPedestalChargevsN->Fill(fTotalEntries,q) > -1); +} + +void MHExtractedSignalPix::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 *MHExtractedSignalPix::DrawClone(Option_t *option) const +{ + + gROOT->SetSelectedPad(NULL); + + MHExtractedSignalPix *newobj = (MHExtractedSignalPix*)Clone(); + + if (!newobj) + return 0; + newobj->SetBit(kCanDelete); + + + if (strlen(option)) + newobj->Draw(option); + else + newobj->Draw(GetDrawOption()); + + return newobj; +} + + +// ------------------------------------------------------------------------- +// +// Draw the histogram +// +void MHExtractedSignalPix::Draw(Option_t *opt) +{ + gStyle->SetOptFit(1); + gStyle->SetOptStat(111111); + + 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(); +} + + +Bool_t MHExtractedSignalPix::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("width"); + 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); + + // + // If we are not able to fit, try once again + // + if (fGausFit->GetProb() < gkProbLimit) + { + + Axis_t rtry = fGausFit->GetParameter(1) - 2.0*fGausFit->GetParameter(2); + rmin = (rtry < rmin ? rmin : rtry); + rmax = fGausFit->GetParameter(1) + 2.0*fGausFit->GetParameter(2); + 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)) + { + CLRBIT(fFlags,kFitOK); + return kFALSE; + } + + if (fChargeProbGetBinLowEdge(fHPedestalCharge->GetXaxis()->GetFirst()); + fChargeLast = fHPedestalCharge->GetBinLowEdge(fHPedestalCharge->GetXaxis()->GetLast()) + +fHPedestalCharge->GetBinWidth(0); + + MH::CutEdges(fHPedestalChargevsN,0); +} + +void MHExtractedSignalPix::Print(const Option_t *o) const +{ + + gLog << 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/MHExtractedSignalPix.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHExtractedSignalPix.h (revision 3011) +++ trunk/MagicSoft/Mars/manalysis/MHExtractedSignalPix.h (revision 3011) @@ -0,0 +1,94 @@ +#ifndef MARS_MHExtractedSignalPix +#define MARS_MHExtractedSignalPix + +#ifndef ROOT_TObject +#include +#endif + +class TH1F; +class TH1I; +class TF1; +#include "TH1.h" // Axis_t + +class MHExtractedSignalPix : public TObject +{ +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 + + 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: + + MHExtractedSignalPix(/*const char *name=NULL, const char *title=NULL*/); + ~MHExtractedSignalPix(); + + 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 IsFitted() const; + Bool_t IsEmpty() const; + Bool_t IsOscillating() const; + + // Fill histos + Bool_t FillCharge(Float_t q); + Bool_t FillChargevsN(Float_t q); + + // Fits + Bool_t FitCharge(Option_t *option="RQ0"); + + // Draws + void Draw(Option_t *option=""); + TObject *DrawClone(Option_t *option="") const; + + // Prints + void Print(const Option_t *o="") const; + + // Others + void CutAllEdges(); + + ClassDef(MHExtractedSignalPix, 1) // Histograms for each calibrated pixel +}; + +#endif