/* ======================================================================== *\ ! ! * ! * 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::CutEdges(fHPedestalCharge,nbins); fChargeFirst = fHPedestalCharge->GetBinLowEdge(fHPedestalCharge->GetXaxis()->GetFirst()); fChargeLast = fHPedestalCharge->GetBinLowEdge(fHPedestalCharge->GetXaxis()->GetLast()) +fHPedestalCharge->GetBinWidth(0); MH::CutEdges(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; }