/* ======================================================================== *\ ! ! * ! * 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 11/2003 ! ! Copyright: MAGIC Software Development, 2000-2002 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // // MHCalibrationPixel // // // // Performs all the necessary fits to extract the mean number of photons // // out of the derived light flux // // // ////////////////////////////////////////////////////////////////////////////// #include "MHCalibrationPixel.h" #include "MHCalibrationConfig.h" #include "MCalibrationFits.h" #include #include #include #include #include #include #include #include #include "MParList.h" #include "MLog.h" #include "MLogManip.h" ClassImp(MHCalibrationPixel); using namespace std; // -------------------------------------------------------------------------- // // Default Constructor. // MHCalibrationPixel::MHCalibrationPixel(Int_t pix, const char *name, const char *title) : fFitOK(kFALSE), fPixId(pix), fTGausFit(NULL), fQGausFit(NULL), fFitLegend(NULL) { fName = name ? name : "MHCalibrationPixel"; fTitle = title ? title : "Fill the accumulated charges and times of all events and perform fits"; TString qname = "HQ"; TString qtitle = "Distribution of Summed FADC Slices Pixel "; qname += pix; qtitle += pix; // Create a large number of bins, later we will rebin fQfirst = -0.5; fQlast = gkStartQlast - 0.5; fQnbins = gkStartPixelBinNr; fHQ = new TH1I( qname.Data(),qtitle.Data(), fQnbins,fQfirst,fQlast); fHQ->SetXTitle("Sum FADC Slices"); fHQ->SetYTitle("Nr. of events"); fHQ->Sumw2(); TString tname = "HT"; TString ttitle = "Distribution of Mean Arrival Times Pixel "; tname += pix; ttitle += pix; Axis_t tfirst = -0.5; Axis_t tlast = 15.5; Int_t nbins = 16; fHT = new TH1I(tname.Data(),ttitle.Data(), nbins,tfirst,tlast); fHT->SetXTitle("Mean Arrival Times [FADC slice nr]"); fHT->SetYTitle("Nr. of events"); fHT->Sumw2(); TString qvsnname = "HQvsN"; TString qvsntitle = "Sum of Charges vs. Event Number Pixel "; qvsnname += pix; qvsntitle += pix; // We define a reasonable number and later enlarge it if necessary nbins = 20000; Axis_t nfirst = -0.5; Axis_t nlast = (Axis_t)nbins - 0.5; fHQvsN = new TH1I(qvsnname.Data(),qvsntitle.Data(), nbins,nfirst,nlast); fHQvsN->SetXTitle("Event Nr."); fHQvsN->SetYTitle("Sum of FADC slices"); fQChisquare = -1.; fQProb = -1.; fQNdf = -1; fTChisquare = -1.; fTProb = -1.; fTNdf = -1; } MHCalibrationPixel::~MHCalibrationPixel() { delete fHQ; delete fHT; delete fHQvsN; if (fQGausFit) delete fQGausFit; if (fTGausFit) delete fTGausFit; if (fFitLegend) delete fFitLegend; } void MHCalibrationPixel::Reset() { for (Int_t i = fHQ->FindBin(fQfirst); i <= fHQ->FindBin(fQlast); i++) fHQ->SetBinContent(i, 1.e-20); for (Int_t i = 0; i < 16; i++) fHT->SetBinContent(i, 1.e-20); fQlast = gkStartQlast; fHQ->GetXaxis()->SetRangeUser(0.,fQlast); return; } // ------------------------------------------------------------------------- // // Set the binnings and prepare the filling of the histograms // Bool_t MHCalibrationPixel::SetupFill(const MParList *plist) { fHQ->Reset(); fHT->Reset(); return kTRUE; } // ------------------------------------------------------------------------- // // Draw a legend with the fit results // void MHCalibrationPixel::DrawLegend() { fFitLegend = new TPaveText(0.05,0.05,0.95,0.95); if (fFitOK) fFitLegend->SetFillColor(80); else fFitLegend->SetFillColor(2); fFitLegend->SetLabel("Results of the Gauss Fit:"); fFitLegend->SetTextSize(0.05); char line1[32]; sprintf(line1,"Mean: Q_{#mu} = %2.2f #pm %2.2f",fQMean,fQMeanErr); fFitLegend->AddText(line1); char line4[32]; sprintf(line4,"Sigma: #sigma_{Q} = %2.2f #pm %2.2f",fQSigma,fQSigmaErr); fFitLegend->AddText(line4); char line7[32]; sprintf(line7,"#chi^{2} / N_{dof}: %4.2f / %3i",fQChisquare,fQNdf); fFitLegend->AddText(line7); char line8[32]; sprintf(line8,"Probability: %4.2f ",fQProb); fFitLegend->AddText(line8); if (fFitOK) fFitLegend->AddText("Result of the Fit: OK"); else fFitLegend->AddText("Result of the Fit: NOT OK"); fFitLegend->SetBit(kCanDelete); fFitLegend->Draw(); } // ------------------------------------------------------------------------- // // Draw the histogram // void MHCalibrationPixel::Draw(Option_t *opt) { gStyle->SetOptFit(0); gStyle->SetOptStat(1111111); TCanvas *c = MakeDefCanvas(this,600,900); gROOT->SetSelectedPad(NULL); c->Divide(2,2); c->cd(1); gPad->SetBorderMode(0); gPad->SetLogy(1); gPad->SetTicks(); fHQ->DrawCopy(opt); if (fQGausFit) { if (fFitOK) fQGausFit->SetLineColor(kGreen); else fQGausFit->SetLineColor(kRed); fQGausFit->DrawCopy("same"); c->Modified(); c->Update(); } c->cd(2); DrawLegend(); c->Update(); c->cd(3); gStyle->SetOptStat(1111111); gPad->SetLogy(1); fHT->DrawCopy(opt); if (fHT->GetFunction("GausTime")) { TF1 *tfit = fHT->GetFunction("GausTime"); if (tfit->GetProb() < 0.01) tfit->SetLineColor(kRed); else tfit->SetLineColor(kGreen); tfit->DrawCopy("same"); c->Modified(); c->Update(); } c->Modified(); c->Update(); c->cd(4); fHQvsN->DrawCopy(opt); } Bool_t MHCalibrationPixel::FitT(Axis_t rmin, Axis_t rmax, Option_t *option) { if (fTGausFit) return kFALSE; rmin = (rmin != 0.) ? rmin : -0.5; rmax = (rmax != 0.) ? rmax : 15.5; const Stat_t entries = fHT->GetEntries(); const Double_t mu_guess = fHT->GetBinCenter(fHT->GetMaximumBin()); const Double_t sigma_guess = (rmax - rmin)/2.; const Double_t area_guess = entries/gkSq2Pi; fTGausFit = new TF1("GausTime","gaus",rmin,rmax); if (!fTGausFit) { *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl; return kFALSE; } fTGausFit->SetParameters(area_guess,mu_guess,sigma_guess); fTGausFit->SetParNames("Area","#mu","#sigma"); fTGausFit->SetParLimits(0,0.,entries); fTGausFit->SetParLimits(1,rmin,rmax); fTGausFit->SetParLimits(2,0.,rmax-rmin); fHT->Fit("GausTime",option); fTChisquare = fTGausFit->GetChisquare(); fTNdf = fTGausFit->GetNDF(); fTProb = fTGausFit->GetProb(); fTMean = fTGausFit->GetParameter(1); fTSigma = fTGausFit->GetParameter(2); if (fTProb < gkProbLimit) { *fLog << warn << "Fit of the Arrival times failed ! " << endl; return kFALSE; } return kTRUE; } Bool_t MHCalibrationPixel::FitQ(Axis_t rmin, Axis_t rmax, Option_t *option) { if (fQGausFit) return kFALSE; // // Get the fitting ranges // rmin = (rmin != 0.) ? rmin : fQfirst; rmax = (rmax != 0.) ? rmax : fQlast; // // 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 = fHQ->GetEntries(); const Double_t ar_guess = entries/gkSq2Pi; const Double_t mu_guess = fHQ->GetBinCenter(fHQ->GetMaximumBin()); const Double_t si_guess = mu_guess/500.; fQGausFit = new TF1("QGausFit","gaus",rmin,rmax); if (!fQGausFit) { *fLog << err << dbginf << "Could not create fit function for Gauss fit" << endl; return kFALSE; } fQGausFit->SetParameters(ar_guess,mu_guess,si_guess); fQGausFit->SetParNames("Area","#mu","#sigma"); fQGausFit->SetParLimits(0,0.,entries); fQGausFit->SetParLimits(1,rmin,rmax); fQGausFit->SetParLimits(2,0.,rmax-rmin); fHQ->Fit("QGausFit",option); fQChisquare = fQGausFit->GetChisquare(); fQNdf = fQGausFit->GetNDF(); fQProb = fQGausFit->GetProb(); fQMean = fQGausFit->GetParameter(1); fQMeanErr = fQGausFit->GetParError(1); fQSigma = fQGausFit->GetParameter(2); fQSigmaErr = fQGausFit->GetParError(2); // // The fit result is accepted under condition // The Probability is greater than gkProbLimit (default 0.01 == 99%) // if (fQProb < gkProbLimit) { *fLog << warn << "Prob: " << fQProb << " is smaller than the allowed value: " << gkProbLimit << endl; fFitOK = kFALSE; return kFALSE; } fFitOK = kTRUE; return kTRUE; } void MHCalibrationPixel::CutAllEdges() { // // The number 100 is necessary because it is the internal binning // of ROOT functions. A call to SetNpx() does NOT help // If you find another solution which WORKS!!, please tell me!! // Int_t nbins = 100; CutEdges(fHQ,nbins); fQfirst = fHQ->GetBinLowEdge(fHQ->GetXaxis()->GetFirst()); fQlast = fHQ->GetBinLowEdge(fHQ->GetXaxis()->GetLast())+fHQ->GetBinWidth(0); fQnbins = nbins; CutEdges(fHQvsN,0); } void MHCalibrationPixel::PrintQFitResult() { *fLog << "Results of the Summed Charges Fit: " << endl; *fLog << "Chisquare: " << fQChisquare << endl; *fLog << "DoF: " << fQNdf << endl; *fLog << "Probability: " << fQProb << endl; *fLog << endl; } void MHCalibrationPixel::PrintTFitResult() { *fLog << "Results of the Arrival Time Slices Fit: " << endl; *fLog << "Chisquare: " << fTChisquare << endl; *fLog << "Ndf: " << fTNdf << endl; *fLog << "Probability: " << fTProb << endl; *fLog << endl; }