Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 2950) +++ trunk/MagicSoft/Mars/Changelog (revision 2951) @@ -24,4 +24,13 @@ - put all missing InitSize(npixels) into MGeomApply + * manalysis/MPedCalcPedRun.[h,cc] + * manalysis/MPedestalCam.[h,cc] + * manalysis/MPedestalPix.[h,cc] + * manalysis/MHPedestalPixel.[h,cc] + * manalysis/Makefile + * manalysis/AnalysisLinkDef.h + - include possibility to fit the pedestals with Gaussian. + - old methods are not affected. + - add new method with MPedCalcPedRun::SetUseHists() Index: trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2950) +++ trunk/MagicSoft/Mars/manalysis/AnalysisLinkDef.h (revision 2951) @@ -17,4 +17,5 @@ #pragma link C++ class MPedestalPix+; +#pragma link C++ class MHPedestalPixel+; #pragma link C++ class MPedestalCam+; #pragma link C++ class MPedCalcPedRun+; Index: trunk/MagicSoft/Mars/manalysis/MHPedestalPixel.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHPedestalPixel.cc (revision 2951) +++ trunk/MagicSoft/Mars/manalysis/MHPedestalPixel.cc (revision 2951) @@ -0,0 +1,393 @@ +/* ======================================================================== *\ +! +! * +! * 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; + +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MHPedestalPixel::MHPedestalPixel(const char *name, const char *title) + : fPixId(-1), + fChargeNbins(500), + fChargevsNbins(1000), + fChargevsNFirst(-0.5), + fChargevsNLast(999.5), + fChargeFirst(-0.5), + fChargeLast(99.5), + fGausFit(NULL) +{ + + fName = name ? name : "MHPedestalPixel"; + 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 ", + fChargeNbins,fChargeFirst,fChargeLast); + 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 ", + fChargevsNbins,fChargevsNFirst,fChargevsNLast); + 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 = -0.5; + fChargeLast = 99.5; + + if (fGausFit) + delete fGausFit; + + 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; + + TString nameQ = TString(fHPedestalCharge->GetName()); + nameQ += id; + fHPedestalCharge->SetName(nameQ.Data()); + + TString nameQvsN = TString(fHPedestalChargevsN->GetName()); + nameQvsN += id; + fHPedestalChargevsN->SetName(nameQvsN.Data()); + + TString titleQ = TString(fHPedestalCharge->GetTitle()); + titleQ += id; + fHPedestalCharge->SetTitle(titleQ.Data()); + + TString titleQvsN = TString(fHPedestalChargevsN->GetTitle()); + titleQvsN += id; + fHPedestalChargevsN->SetTitle(titleQvsN.Data()); + +} + + +Bool_t MHPedestalPixel::SetupFill(const MParList *plist) +{ + + Reset(); + + return kTRUE; +} + + + + +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->SetOptFit(1); + gStyle->SetOptStat(111111); + + gROOT->SetSelectedPad(NULL); + + TCanvas *c = MakeDefCanvas(this,600,900); + + c->Divide(1,2); + + c->cd(1); + gPad->SetBorderMode(0); + gPad->SetTicks(); + + if (fHPedestalCharge->Integral() > 0) + gPad->SetLogy(1); + else + gPad->SetLogy(0); + + fHPedestalCharge->Draw(opt); + + c->Modified(); + c->Update(); + + if (fGausFit) + { + if (IsFitOK()) + fGausFit->SetLineColor(kGreen); + else + fGausFit->SetLineColor(kRed); + + fGausFit->Draw("same"); + } + + c->Modified(); + c->Update(); + + c->cd(2); + gPad->SetTicks(); + + fHPedestalChargevsN->Draw(opt); + c->Modified(); + c->Update(); + + return; +} + + +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 area_guess = entries/gkSq2Pi; + const Double_t mu_guess = fHPedestalCharge->GetBinCenter(fHPedestalCharge->GetMaximumBin()); + const Double_t sigma_guess = mu_guess/15.; + + TString name = TString("GausFit"); + name += fPixId; + + fGausFit = new TF1(name.Data(),"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) - 3.0*fGausFit->GetParameter(2); + rmin = (rtry < rmin ? rmin : rtry); + rmax = fGausFit->GetParameter(1) + 3.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 ((fChargeProb < gkProbLimit) || (TMath::IsNaN(fChargeProb))) + { + CLRBIT(fFlags,kFitOK); + return kFALSE; + } + + SETBIT(fFlags,kFitOK); + return kTRUE; +} + +void MHPedestalPixel::CutAllEdges() +{ + + Int_t nbins = 50; + + CutEdges(fHPedestalCharge,nbins); + + fChargeFirst = fHPedestalCharge->GetBinLowEdge(fHPedestalCharge->GetXaxis()->GetFirst()); + fChargeLast = fHPedestalCharge->GetBinLowEdge(fHPedestalCharge->GetXaxis()->GetLast()) + +fHPedestalCharge->GetBinWidth(0); + + CutEdges(fHPedestalChargevsN,0); + +} + +void MHPedestalPixel::Print(const Option_t *o) const +{ + + *fLog << all << "Pedestal Fits Pixel: " << fPixId << endl; + + if (TESTBIT(fFlags,kFitted)) + { + + *fLog << all << "Results of the Summed Charges Fit: " << endl; + *fLog << all << "Chisquare: " << fChargeChisquare << endl; + *fLog << all << "DoF: " << fChargeNdf << endl; + *fLog << all << "Probability: " << fChargeProb << endl; + *fLog << all << "Results of fit: " ; + if (TESTBIT(fFlags,kFitOK)) + *fLog << inf << "OK" << endl; + else + *fLog << err << "NOT OK" << endl; + *fLog << all << endl; + } + else + { + *fLog << all << "Pedestal Histogram has not yet been fitted" << endl; + *fLog << all << endl; + } + +} + Index: trunk/MagicSoft/Mars/manalysis/MHPedestalPixel.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MHPedestalPixel.h (revision 2951) +++ trunk/MagicSoft/Mars/manalysis/MHPedestalPixel.h (revision 2951) @@ -0,0 +1,93 @@ +#ifndef MARS_MHPedestalPixel +#define MARS_MHPedestalPixel + +#ifndef MARS_MH +#include "MH.h" +#endif + +class TH1F; +class TH1I; +class TF1; +class MHPedestalPixel : public MH +{ + +private: + + Int_t fPixId; // Pixel Nr + + const Int_t fChargeNbins; + const Int_t fChargevsNbins; + + const Axis_t fChargevsNFirst; + const Axis_t fChargevsNLast; + + Axis_t fChargeFirst; + Axis_t fChargeLast; + + 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 }; + +public: + + MHPedestalPixel(const char *name=NULL, const char *title=NULL); + ~MHPedestalPixel(); + + void Clear(Option_t *o=""); + void Reset(); + + void ChangeHistId(Int_t i); + + // Setters + 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; } + + Bool_t IsFitOK() const; + Bool_t IsEmpty() const; + + // Fill histos + Bool_t FillCharge(Float_t q); + Bool_t FillChargevsN(Float_t q); + + Bool_t SetupFill(const MParList *pList); + Bool_t Fill(const MParContainer *, const Stat_t w=1) { return kTRUE; } + + // 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(MHPedestalPixel, 1) // Histograms for each calibrated pixel +}; + +#endif Index: trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc (revision 2950) +++ trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.cc (revision 2951) @@ -58,4 +58,8 @@ #include "MPedestalCam.h" +#include "MExtractedSignalPix.h" +#include "MExtractedSignalCam.h" + + #include "MGeomCamMagic.h" @@ -69,4 +73,5 @@ // MPedCalcPedRun::MPedCalcPedRun(const char *name, const char *title) + : fRawEvt(NULL), fPedestals(NULL), fSignals(NULL) { fName = name ? name : "MPedCalcPedRun"; @@ -75,5 +80,24 @@ AddToBranchList("fHiGainPixId"); AddToBranchList("fHiGainFadcSamples"); -} + + Clear(); +} + +void MPedCalcPedRun::Clear(const Option_t *o) +{ + + fNumHiGainSamples = 0; + fNumPixels = 0; + fNumSamplesTot = 0; + fUseHists = kFALSE; + + fRawEvt = NULL; + fPedestals = NULL; + fSignals = NULL; + + return; + +} + // -------------------------------------------------------------------------- @@ -101,4 +125,12 @@ return kFALSE; + fSignals = (MExtractedSignalCam*)pList->FindObject("MExtractedSignalCam"); + + if (!fSignals && fUseHists) + { + *fLog << warn << "Cannot find MExtractedSignalCam... will not use histograms!" << endl; + fUseHists = kFALSE; + } + fNumSamplesTot=0; @@ -141,4 +173,7 @@ // the other one fNumHiGainSamples = runheader->GetNumSamplesHiGain() & ~1; + + if (fUseHists) + fPedestals->InitUseHists(); return kTRUE; @@ -170,11 +205,13 @@ const Float_t msum = (Float_t)sum; - const Float_t msqr = (Float_t)sqr; - - const Float_t higainped = msum/fNumHiGainSamples; - const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSamples)/(fNumHiGainSamples-1.)); - + const UInt_t idx = pixel.GetPixelId(); - (*fPedestals)[idx].Set(higainped, higainrms); + // + // These three lines have been uncommented by Markus Gaug + // If anybody needs them, please contact me!! + // + // const Float_t higainped = msum/fNumHiGainSamples; + // const Float_t higainrms = TMath::Sqrt((msqr-msum*msum/fNumHiGainSamples)/(fNumHiGainSamples-1.)); + // (*fPedestals)[idx].Set(higainped, higainrms); fSumx[idx] += msum; @@ -182,4 +219,5 @@ // The old version: // + // const Float_t msqr = (Float_t)sqr; // fSumx2[idx] += msqr; // @@ -187,8 +225,17 @@ // fSumx2[idx] += msum*msum; + + if (fUseHists) + { + MExtractedSignalPix &sig = (*fSignals)[idx]; + const Float_t signal = sig.GetExtractedSignalHiGain(); + const Float_t signalPerSlice = signal/(Float_t)fSignals->GetNumUsedFADCSlices(); + (*fPedestals)[idx].FillHists(signalPerSlice); + } } fPedestals->SetReadyToSave(); fNumSamplesTot += fNumHiGainSamples; + return kTRUE; @@ -225,4 +272,8 @@ (*fPedestals)[pixid].Set(higainped, higainrms); + + if (fUseHists) + (*fPedestals)[pixid].FitCharge(); + } Index: trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h (revision 2950) +++ trunk/MagicSoft/Mars/manalysis/MPedCalcPedRun.h (revision 2951) @@ -18,4 +18,5 @@ class MRawEvtData; class MPedestalCam; +class MExtractedSignalCam; class MPedCalcPedRun : public MTask @@ -25,10 +26,13 @@ ULong_t fNumSamplesTot; - MRawEvtData *fRawEvt; // raw event data (time slices) - MPedestalCam *fPedestals; // Pedestals of all pixels in the camera + MRawEvtData *fRawEvt; // raw event data (time slices) + MPedestalCam *fPedestals; // Pedestals of all pixels in the camera + MExtractedSignalCam *fSignals; // Signals of all pixels in the camera TArrayF fSumx; // sum of values TArrayF fSumx2; // sum of squared values + Bool_t fUseHists; + Bool_t ReInit(MParList *pList); @@ -40,4 +44,9 @@ MPedCalcPedRun(const char *name=NULL, const char *title=NULL); + void Clear(const Option_t *o=""); + + void SetUseHists(const Bool_t b = kTRUE) { fUseHists = b; } + Bool_t IsUseHists() const { return fUseHists; } + ClassDef(MPedCalcPedRun, 0) // Task to calculate pedestals from pedestal runs raw data }; Index: trunk/MagicSoft/Mars/manalysis/MPedestalCam.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalCam.cc (revision 2950) +++ trunk/MagicSoft/Mars/manalysis/MPedestalCam.cc (revision 2951) @@ -110,4 +110,9 @@ } +void MPedestalCam::InitUseHists() +{ + fArray->ForEach(MPedestalPix, InitUseHists)(); +} + void MPedestalCam::Print(Option_t *o) const { @@ -175,20 +180,69 @@ Bool_t MPedestalCam::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const { - switch (type) + + switch (type) { case 0: + if ((*this)[idx].IsValid()) val = (*this)[idx].GetPedestal(); - break; + else + return kFALSE; + break; case 1: + if ((*this)[idx].IsValid()) val = (*this)[idx].GetPedestalRms(); - break; + else + return kFALSE; + break; + case 2: + // if ((*this)[idx].IsFitValid()) + val = (*this)[idx].GetMean(); + // else + // return kFALSE; + break; + case 3: + // if ((*this)[idx].IsFitValid()) + val = (*this)[idx].GetMeanErr(); + // else + // return kFALSE; + break; + case 4: + // if ((*this)[idx].IsFitValid()) + val = (*this)[idx].GetSigma(); + // else + // return kFALSE; + break; + case 5: + // if ((*this)[idx].IsFitValid()) + val = (*this)[idx].GetSigmaErr(); + // else + // return kFALSE; + break; + case 6: + // if ((*this)[idx].IsFitValid()) + val = (*this)[idx].GetProb(); + // else + // return kFALSE; + break; + case 7: + // if ((*this)[idx].IsFitValid()) + val = ((*this)[idx].GetPedestal()-(*this)[idx].GetMean()); + // else + // return kFALSE; + break; + case 8: + // if ((*this)[idx].IsFitValid()) + val = ((*this)[idx].GetPedestalRms()-(*this)[idx].GetSigma()); + // else + // return kFALSE; + break; default: - return kFALSE; - } - return val>=0; -} - -void MPedestalCam::DrawPixelContent(Int_t num) const -{ - *fLog << warn << "MPedestalCam::DrawPixelContent - not available." << endl; -} + return kFALSE; + } + return kTRUE; +} + +void MPedestalCam::DrawPixelContent(Int_t idx) const +{ + (*this)[idx].Draw(); +} Index: trunk/MagicSoft/Mars/manalysis/MPedestalCam.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalCam.h (revision 2950) +++ trunk/MagicSoft/Mars/manalysis/MPedestalCam.h (revision 2951) @@ -21,5 +21,6 @@ void Clear(Option_t *o=""); - + void InitUseHists(); + void InitSize(const UInt_t i); Int_t GetSize() const; @@ -34,5 +35,5 @@ Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const; - void DrawPixelContent(Int_t num) const; + void DrawPixelContent(Int_t idx) const; ClassDef(MPedestalCam, 1) // Storage Container for all pedestal information of the camera Index: trunk/MagicSoft/Mars/manalysis/MPedestalPix.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalPix.cc (revision 2950) +++ trunk/MagicSoft/Mars/manalysis/MPedestalPix.cc (revision 2951) @@ -34,11 +34,21 @@ #include "MLog.h" +#include "MLogManip.h" ClassImp(MPedestalPix); -MPedestalPix::MPedestalPix() +using namespace std; + +MPedestalPix::MPedestalPix() : fHist(NULL), fFlags(0) { - Clear(); + Clear(); } + +MPedestalPix::~MPedestalPix() +{ + if (fHist) + delete fHist; +} + // ------------------------------------------------------------------------ @@ -48,6 +58,115 @@ void MPedestalPix::Clear(Option_t *o) { - fPedestal = 0.; - fPedestalRms = 0.; + fPedestal = -1.; + fPedestalRms = -1.; + + fMean = -1.; + fMeanErr = -1.; + fSigma = -1.; + fProb = -1.; + + CLRBIT(fFlags,kFitted); + CLRBIT(fFlags,kFitValid); + } +void MPedestalPix::InitUseHists() +{ + fPedestal = 0.; + fPedestalRms = 0.; + + fHist = new MHPedestalPixel("MHPedestalPixel","Pedestal Histograms Pixel "); + + if (!fHist) + *fLog << warn << dbginf << " Could not create MHCalibrationPixel " << endl; + +} + +void MPedestalPix::FillHists(const Float_t f) +{ + + if (!fHist) + return; + + fHist->FillCharge(f); + fHist->FillChargevsN(f); +} + + +void MPedestalPix::Set(Float_t m, Float_t r) +{ + fPedestal = m; + fPedestalRms = r; +} + +Bool_t MPedestalPix::IsValid() const +{ + return fPedestal>=0||fPedestalRms>=0; +} + +Bool_t MPedestalPix::IsFitted() const +{ + return TESTBIT(fFlags,kFitted); +} + +Bool_t MPedestalPix::IsFitValid() const +{ + return TESTBIT(fFlags,kFitValid); +} + +// -------------------------------------------------------------------------- +// +// 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 +// +void MPedestalPix::FitCharge() +{ + + // + // 1) Return if the charge distribution is already succesfully fitted + // or if the histogram is empty + // + if (fHist->IsFitOK() || fHist->IsEmpty()) + return; + + + fHist->CutAllEdges(); + + // + // 2) Fit the Lo Gain histograms with a Gaussian + // + if(fHist->FitCharge()) + SETBIT(fFlags,kFitted); + else + CLRBIT(fFlags,kFitted); + + // + // 6) Retrieve the results and store them in this class + // + fMean = fHist->GetChargeMean(); + fMeanErr = fHist->GetChargeMeanErr(); + fSigma = fHist->GetChargeSigma(); + fSigmaErr = fHist->GetChargeSigmaErr(); + fProb = fHist->GetChargeProb(); + + if (CheckFitValidity()) + SETBIT(fFlags,kFitValid); + else + CLRBIT(fFlags,kFitValid); + + return; +} + +Bool_t MPedestalPix::CheckFitValidity() +{ + + if (fProb < 0.001) + return kFALSE; + + return kTRUE; + +} Index: trunk/MagicSoft/Mars/manalysis/MPedestalPix.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MPedestalPix.h (revision 2950) +++ trunk/MagicSoft/Mars/manalysis/MPedestalPix.h (revision 2951) @@ -6,4 +6,9 @@ #endif +#ifndef MARS_MHPedestalPixel +#include "MHPedestalPixel.h" +#endif + +class MHPedestalPixel; class MPedestalPix : public MParContainer { @@ -12,23 +17,53 @@ Float_t fPedestalRms; // root mean square / sigma / standard deviation of pedestal + MHPedestalPixel *fHist; // Pointer to the histograms performing the fits, etc. + + Float_t fMean; + Float_t fMeanErr; + Float_t fSigma; + Float_t fSigmaErr; + Float_t fProb; + + Byte_t fFlags; + + enum { kFitted, kFitValid }; + public: MPedestalPix(); - + ~MPedestalPix(); + void Clear(Option_t *o=""); - Float_t GetPedestal() const { return fPedestal; } + // Using histograms + void InitUseHists(); + void FillHists(const Float_t f); + + // Setters + void SetPedestal(const Float_t f) { fPedestal = f; } + void SetPedestalRms(const Float_t f) { fPedestalRms = f; } + + void Set(const Float_t m, const Float_t r); + + // Getters + Float_t GetPedestal() const { return fPedestal; } Float_t GetPedestalRms() const { return fPedestalRms; } -/* - void SetMean(Float_t f) { fMean = f; } - void SetSigma(Float_t f) { fSigma = f; } - void SetMeanRms(Float_t f) { fMeanRms = f; } - void SetSigmaRms(Float_t f) { fSigmaRms = f; } - */ - void SetPedestal(Float_t f) { fPedestal = f; } - void SetPedestalRms(Float_t f) { fPedestalRms = f; } - void Set(Float_t m, Float_t r) { fPedestal = m; fPedestalRms = r; } - - Bool_t IsValid() const { return fPedestal>=0||fPedestalRms>=0; } + Float_t GetMean() const { return fMean; } + Float_t GetMeanErr() const { return fMeanErr; } + Float_t GetSigma() const { return fSigma; } + Float_t GetSigmaErr() const { return fSigmaErr; } + Float_t GetProb() const { return fProb; } + + Bool_t IsValid() const; + Bool_t IsFitted() const; + Bool_t IsFitValid() const; + + // Fits + void FitCharge(); + Bool_t CheckFitValidity(); + + // Draws + void Draw(Option_t *opt="") { fHist->Draw(opt); } + TObject *DrawClone(Option_t *opt="") const { return fHist->DrawClone(opt); } ClassDef(MPedestalPix, 1) // Storage Container for Pedestal information of one pixel Index: trunk/MagicSoft/Mars/manalysis/Makefile =================================================================== --- trunk/MagicSoft/Mars/manalysis/Makefile (revision 2950) +++ trunk/MagicSoft/Mars/manalysis/Makefile (revision 2951) @@ -33,4 +33,5 @@ SRCFILES = MPedestalPix.cc \ + MHPedestalPixel.cc \ MPedestalCam.cc \ MPedPhotPix.cc \