/* ======================================================================== *\ ! ! * ! * 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 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // // MHCalibrationBlindPixel // // // // Performs all the necessary fits to extract the mean number of photons // // out of the derived light flux // // // ////////////////////////////////////////////////////////////////////////////// #include "MHCalibrationBlindPixel.h" #include "MHCalibrationConfig.h" #include "MCalibrationFits.h" #include #include #include #include #include #include #include #include #include #include #include "MBinning.h" #include "MParList.h" #include "MLog.h" #include "MLogManip.h" ClassImp(MHCalibrationBlindPixel); using namespace std; // -------------------------------------------------------------------------- // // Default Constructor. // MHCalibrationBlindPixel::MHCalibrationBlindPixel(const char *name, const char *title) : fSinglePheFit(NULL), fTimeGausFit(NULL) { fName = name ? name : "MHCalibrationBlindPixel"; fTitle = title ? title : "Fill the accumulated charges and times all Blind Pixel events and perform fits"; // Create a large number of bins, later we will rebin fBPQfirst = 0; fBPQlast = gkStartBlindPixelBinNr; fBPQnbins = gkStartBlindPixelBinNr; fHBPQ = new TH1I("HBPQ","Distribution of Summed FADC Slices",fBPQnbins,fBPQfirst,fBPQlast); fHBPQ->SetXTitle("Sum FADC Slices"); fHBPQ->SetYTitle("Nr. of events"); fHBPQ->Sumw2(); fErrBPQfirst = 0.; fErrBPQlast = gkStartBlindPixelBinNr; fErrBPQnbins = gkStartBlindPixelBinNr; fHBPErrQ = new TH1F("HBPErrQ","Distribution of Variances of Summed FADC Slices", fErrBPQnbins,fErrBPQfirst,fErrBPQlast); fHBPErrQ->SetXTitle("Variance Summed FADC Slices"); fHBPErrQ->SetYTitle("Nr. of events"); fHBPErrQ->Sumw2(); Axis_t tfirst = -0.5; Axis_t tlast = 15.5; Int_t nbins = 16; fHBPT = new TH1I("HBPT","Distribution of Mean Arrival Times",nbins,tfirst,tlast); fHBPT->SetXTitle("Mean Arrival Times [FADC slice nr]"); fHBPT->SetYTitle("Nr. of events"); fHBPT->Sumw2(); // 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; fHBPQvsN = new TH1I("HBPQvsN","Sum of Charges vs. Event Number",nbins,nfirst,nlast); fHBPQvsN->SetXTitle("Event Nr."); fHBPQvsN->SetYTitle("Sum of FADC slices"); fgSinglePheFitFunc = &gfKto8; fgSinglePheFitNPar = 5; } MHCalibrationBlindPixel::~MHCalibrationBlindPixel() { delete fHBPQ; delete fHBPT; delete fHBPErrQ; if (fSinglePheFit) delete fSinglePheFit; if (fTimeGausFit) delete fTimeGausFit; } void MHCalibrationBlindPixel::ResetBin(Int_t i) { fHBPQ->SetBinContent (i, 1.e-20); fHBPErrQ->SetBinContent (i, 1.e-20); fHBPT->SetBinContent(i, 1.e-20); } // ------------------------------------------------------------------------- // // Draw a legend with the fit results // void MHCalibrationBlindPixel::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 single PhE Fit (to k=6):"); fFitLegend->SetTextSize(0.05); char line1[32]; sprintf(line1,"Mean: #lambda = %2.2f #pm %2.2f",GetLambda(),GetLambdaErr()); fFitLegend->AddText(line1); char line2[32]; sprintf(line2,"Pedestal: #mu_{0} = %2.2f #pm %2.2f",GetMu0(),GetMu0Err()); fFitLegend->AddText(line2); char line3[32]; sprintf(line3,"Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",GetSigma0(),GetSigma0Err()); fFitLegend->AddText(line3); char line4[32]; sprintf(line4,"1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",GetMu1(),GetMu1Err()); fFitLegend->AddText(line4); char line5[32]; sprintf(line5,"Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",GetSigma1(),GetSigma1Err()); fFitLegend->AddText(line5); char line7[32]; sprintf(line7,"#chi^{2} / N_{dof}: %4.2f / %3i",GetChiSquare(),GetNdf()); fFitLegend->AddText(line7); char line8[32]; sprintf(line8,"Probability: %4.2f ",GetProb()); 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 MHCalibrationBlindPixel::Draw(Option_t *opt) { gStyle->SetOptFit(0); gStyle->SetOptStat(1111111); TCanvas *c = MakeDefCanvas(this,550,700); c->Divide(2,2); gROOT->SetSelectedPad(NULL); c->cd(1); gPad->SetLogy(1); gPad->SetTicks(); fHBPQ->DrawCopy(opt); if (fSinglePheFit) { if (fFitOK) fSinglePheFit->SetLineColor(kGreen); else fSinglePheFit->SetLineColor(kRed); fSinglePheFit->DrawCopy("same"); c->Modified(); c->Update(); } c->cd(2); DrawLegend(); c->Modified(); c->Update(); c->cd(3); gPad->SetLogy(1); gPad->SetBorderMode(0); fHBPT->DrawCopy(opt); if (fHBPT->GetFunction("GausTime")) { TF1 *tfit = fHBPT->GetFunction("GausTime"); if (tfit->GetProb() < 0.01) tfit->SetLineColor(kRed); else tfit->SetLineColor(kGreen); tfit->DrawCopy("same"); c->Modified(); c->Update(); } c->cd(4); fHBPQvsN->DrawCopy(opt); c->Modified(); c->Update(); } Bool_t MHCalibrationBlindPixel::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1) { gRandom->SetSeed(); if (fHBPQ->GetEntries() != 0) { *fLog << err << "Histogram " << fHBPQ->GetTitle() << " is already filled. " << endl; *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl; return kFALSE; } TF1 *simulateSinglePhe = new TF1("simulateSinglePhe",fgSinglePheFitFunc, fBPQfirst,fBPQlast,fgSinglePheFitNPar); simulateSinglePhe->SetParameters(lambda,mu0,mu1,sigma0,sigma1); simulateSinglePhe->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1"); simulateSinglePhe->SetNpx(fBPQnbins); for (Int_t i=0;i<10000; i++) { fHBPQ->Fill(simulateSinglePhe->GetRandom()); } return kTRUE; } void MHCalibrationBlindPixel::ChangeFitFunc(BPFitFunc fitfunc, Int_t par) { fgSinglePheFitFunc = fitfunc; fgSinglePheFitNPar = par; } Bool_t MHCalibrationBlindPixel::FitSinglePhe(Axis_t rmin, Axis_t rmax, Option_t *opt) { if (fSinglePheFit) return kFALSE; // // Get the fitting ranges // rmin = (rmin != 0.) ? rmin : fBPQfirst; rmax = (rmax != 0.) ? rmax : fBPQlast; // // 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 = fHBPQ->GetSumOfWeights(); const Double_t lambda_guess = 0.2; const Double_t mu_0_guess = fHBPQ->GetBinCenter(fHBPQ->GetMaximumBin()); const Double_t si_0_guess = mu_0_guess/500.; const Double_t mu_1_guess = mu_0_guess + 2.; const Double_t si_1_guess = si_0_guess + si_0_guess; fSinglePheFit = new TF1("SinglePheFit",fgSinglePheFitFunc,rmin,rmax,fgSinglePheFitNPar+1); fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,entries); fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","area"); fSinglePheFit->SetParLimits(0,0.,5.); fSinglePheFit->SetParLimits(1,rmin,rmax); fSinglePheFit->SetParLimits(2,rmin,rmax); fSinglePheFit->SetParLimits(3,1.0,rmax-rmin); fSinglePheFit->SetParLimits(4,1.7,rmax-rmin); fSinglePheFit->SetParLimits(5,0.,2.*entries); // // Normalize the histogram to facilitate faster fitting of the area // For speed reasons, FKto8 is normalized to Sqrt(2 pi). // Therefore also normalize the histogram to that value // // ROOT gives us another nice example of user-unfriendly behavior: // Although the normalization of the function fSinglePhe and the // Histogram fHBPQ agree (!!), the fit does not normalize correctly INTERNALLY // in the fitting procedure !!! // // This has to do with the fact that the internal function histogramming // uses 100 bins and does not adapt to the binning of the fitted histogram, unlike PAW does // (very important if you use Sumw2, see e.g. ROOTTALK: Mon May 26 1997 - 09:56:03 MEST) // // So, WE have to adapt to that internal flaw of ROOT: // const Int_t npx = fSinglePheFit->GetNpx(); const Int_t bins = fHBPQ->GetXaxis()->GetLast()-fHBPQ->GetXaxis()->GetFirst(); // fHBPQ->Scale(gkSq2Pi*(float)bins/npx/entries); // // we need this, otherwise, ROOT does not calculate the area correctly // don't ask me why it does not behave correctly, it's one of the nasty // mysteries of ROOT which takes you a whole day to find out :-) // // fSinglePheFit->SetNpx(fQnbins); fHBPQ->Fit("SinglePheFit",opt); fLambda = fSinglePheFit->GetParameter(0); fMu0 = fSinglePheFit->GetParameter(1); fMu1 = fSinglePheFit->GetParameter(2); fSigma0 = fSinglePheFit->GetParameter(3); fSigma1 = fSinglePheFit->GetParameter(4); fLambdaErr = fSinglePheFit->GetParError(0); fMu0Err = fSinglePheFit->GetParError(1); fMu1Err = fSinglePheFit->GetParError(2); fSigma0Err = fSinglePheFit->GetParError(3); fSigma1Err = fSinglePheFit->GetParError(4); fProb = fSinglePheFit->GetProb(); fChisquare = fSinglePheFit->GetChisquare(); fNdf = fSinglePheFit->GetNDF(); *fLog << "Results of the Blind Pixel Fit: " << endl; *fLog << "Chisquare: " << fChisquare << endl; *fLog << "DoF: " << fNdf << endl; *fLog << "Probability: " << fProb << endl; *fLog << "Integral: " << fSinglePheFit->Integral(rmin,rmax); // // The fit result is accepted under condition // The Probability is greater than gkProbLimit (default 0.01 == 99%) // if (fProb < gkProbLimit) { *fLog << warn << "Prob: " << fProb << " is smaller than the allowed value: " << gkProbLimit << endl; fFitOK = kFALSE; return kFALSE; } fFitOK = kTRUE; return kTRUE; } void MHCalibrationBlindPixel::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; *fLog << "New number of bins in HSinQ: " << CutEdges(fHBPQ,nbins) << endl; fBPQfirst = fHBPQ->GetBinLowEdge(fHBPQ->GetXaxis()->GetFirst()); fBPQlast = fHBPQ->GetBinLowEdge(fHBPQ->GetXaxis()->GetLast())+fHBPQ->GetBinWidth(0); fBPQnbins = nbins; *fLog << "New number of bins in HErrQ: " << CutEdges(fHBPErrQ,30) << endl; fErrBPQfirst = fHBPErrQ->GetBinLowEdge(fHBPErrQ->GetXaxis()->GetFirst()); fErrBPQlast = fHBPErrQ->GetBinLowEdge(fHBPErrQ->GetXaxis()->GetLast())+fHBPErrQ->GetBinWidth(0); fErrBPQnbins = nbins; CutEdges(fHBPQvsN,0); } Bool_t MHCalibrationBlindPixel::FitT(Axis_t rmin, Axis_t rmax, Option_t *opt) { if (fTimeGausFit) return kFALSE; rmin = (rmin != 0.) ? rmin : 0.; rmax = (rmax != 0.) ? rmax : 16.; const Stat_t entries = fHBPT->GetEntries(); const Double_t mu_guess = fHBPT->GetBinCenter(fHBPT->GetMaximumBin()); const Double_t sigma_guess = (rmax - rmin)/2.; const Double_t area_guess = entries/gkSq2Pi; fTimeGausFit = new TF1("GausTime","gaus",rmin,rmax); fTimeGausFit->SetParameters(area_guess,mu_guess,sigma_guess); fTimeGausFit->SetParNames("Area","#mu","#sigma"); fTimeGausFit->SetParLimits(0,0.,entries); fTimeGausFit->SetParLimits(1,rmin,rmax); fTimeGausFit->SetParLimits(2,0.,rmax-rmin); fHBPT->Fit("GausTime",opt); fMeanT = fTimeGausFit->GetParameter(2); fSigmaT = fTimeGausFit->GetParameter(3); fMeanTErr = fTimeGausFit->GetParError(2); fSigmaTErr = fTimeGausFit->GetParError(3); Float_t prob = fTimeGausFit->GetProb(); *fLog << "Results of the Times Fit: " << endl; *fLog << "Chisquare: " << fTimeGausFit->GetChisquare() << endl; *fLog << "Ndf: " << fTimeGausFit->GetNDF() << endl; *fLog << "Probability: " << prob << endl; if (prob < gkProbLimit) { *fLog << warn << "Fit of the Arrival times failed ! " << endl; return kFALSE; } return kTRUE; }