/* ======================================================================== *\ ! ! * ! * 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 Single Photo-Electron Fit to extract // // the mean number of photons and to derive the light flux // // // // The fit result is accepted under condition that: // // 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%) // // 2) at least 100 events are in the single Photo-electron peak // // // ////////////////////////////////////////////////////////////////////////////// #include "MHCalibrationBlindPixel.h" #include "MHCalibrationConfig.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), fSinglePhePedFit(NULL), fLambda(0.), fMu0(0.), fMu1(0.), fSigma0(0.), fSigma1(0.), fLambdaErr(0.), fMu0Err(0.), fMu1Err(0.), fSigma0Err(0.), fSigma1Err(0.), fChisquare(0.), fProb(0.), fNdf(0), fMeanTime(0.), fMeanTimeErr(0.), fSigmaTime(0.), fSigmaTimeErr(0.), fLambdaCheck(0.), fLambdaCheckErr(0.), fgSinglePheFitFunc(&gfKto4),fgSinglePheFitNPar(6) { 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 fBlindPixelChargefirst = -1000.; fBlindPixelChargelast = gkStartBlindPixelBinNr; fBlindPixelChargenbins = gkStartBlindPixelBinNr+(int)fBlindPixelChargefirst; fHBlindPixelCharge = new TH1F("HBlindPixelCharge","Distribution of Summed FADC Slices", fBlindPixelChargenbins,fBlindPixelChargefirst,fBlindPixelChargelast); fHBlindPixelCharge->SetXTitle("Sum FADC Slices"); fHBlindPixelCharge->SetYTitle("Nr. of events"); fHBlindPixelCharge->Sumw2(); fHBlindPixelCharge->SetDirectory(NULL); Axis_t tfirst = -0.5; Axis_t tlast = 15.5; Int_t nbins = 16; fHBlindPixelTime = new TH1I("HBlindPixelTime","Distribution of Mean Arrival Times",nbins,tfirst,tlast); fHBlindPixelTime->SetXTitle("Mean Arrival Times [FADC slice nr]"); fHBlindPixelTime->SetYTitle("Nr. of events"); fHBlindPixelTime->Sumw2(); fHBlindPixelTime->SetDirectory(NULL); // 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; fHBlindPixelChargevsN = new TH1I("HBlindPixelChargevsN","Sum of Charges vs. Event Number",nbins,nfirst,nlast); fHBlindPixelChargevsN->SetXTitle("Event Nr."); fHBlindPixelChargevsN->SetYTitle("Sum of FADC slices"); fHBlindPixelChargevsN->SetDirectory(NULL); } MHCalibrationBlindPixel::~MHCalibrationBlindPixel() { delete fHBlindPixelCharge; delete fHBlindPixelTime; if (fSinglePheFit) delete fSinglePheFit; if (fSinglePhePedFit) delete fSinglePhePedFit; if (fTimeGausFit) delete fTimeGausFit; } void MHCalibrationBlindPixel::ResetBin(Int_t i) { fHBlindPixelCharge->SetBinContent (i, 1.e-20); fHBlindPixelTime->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); const TString line1 = Form("Mean: #lambda = %2.2f #pm %2.2f",GetLambda(),GetLambdaErr()); fFitLegend->AddText(line1); const TString line6 = Form("Mean #lambda (check) = %2.2f #pm %2.2f",GetLambdaCheck(),GetLambdaCheckErr()); fFitLegend->AddText(line6); const TString line2 = Form("Pedestal: #mu_{0} = %2.2f #pm %2.2f",GetMu0(),GetMu0Err()); fFitLegend->AddText(line2); const TString line3 = Form("Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",GetSigma0(),GetSigma0Err()); fFitLegend->AddText(line3); const TString line4 = Form("1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",GetMu1(),GetMu1Err()); fFitLegend->AddText(line4); const TString line5 = Form("Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",GetSigma1(),GetSigma1Err()); fFitLegend->AddText(line5); const TString line7 = Form("#chi^{2} / N_{dof}: %4.2f / %3i",GetChiSquare(),GetNdf()); fFitLegend->AddText(line7); const TString line8 = Form("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(); } TObject *MHCalibrationBlindPixel::DrawClone(Option_t *opt) const { gROOT->SetSelectedPad(NULL); TObject *newobj = Clone(); newobj->SetBit(kCanDelete); if (!newobj) return 0; newobj->Draw(opt); return newobj; } // ------------------------------------------------------------------------- // // 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(); fHBlindPixelCharge->DrawCopy(opt); if (fSinglePheFit) { if (fFitOK) fSinglePheFit->SetLineColor(kGreen); else fSinglePheFit->SetLineColor(kRed); fSinglePheFit->DrawCopy("same"); c->Modified(); c->Update(); if (fSinglePhePedFit) { fSinglePhePedFit->SetLineColor(kBlue); fSinglePhePedFit->DrawCopy("same"); } } c->cd(2); DrawLegend(); c->Modified(); c->Update(); c->cd(3); gPad->SetLogy(1); gPad->SetBorderMode(0); fHBlindPixelTime->DrawCopy(opt); if (fHBlindPixelTime->GetFunction("GausTime")) { TF1 *tfit = fHBlindPixelTime->GetFunction("GausTime"); if (tfit->GetProb() < 0.01) tfit->SetLineColor(kRed); else tfit->SetLineColor(kGreen); tfit->DrawCopy("same"); c->Modified(); c->Update(); } c->cd(4); fHBlindPixelChargevsN->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 (fHBlindPixelCharge->GetIntegral() != 0) { *fLog << err << "Histogram " << fHBlindPixelCharge->GetTitle() << " is already filled. " << endl; *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl; return kFALSE; } TF1 *simulateSinglePhe = new TF1("simulateSinglePhe",fgSinglePheFitFunc, fBlindPixelChargefirst,fBlindPixelChargelast,fgSinglePheFitNPar); simulateSinglePhe->SetParameters(lambda,mu0,mu1,sigma0,sigma1); simulateSinglePhe->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1"); simulateSinglePhe->SetNpx(fBlindPixelChargenbins); for (Int_t i=0;i<10000; i++) { fHBlindPixelCharge->Fill(simulateSinglePhe->GetRandom()); } return kTRUE; } void MHCalibrationBlindPixel::ChangeFitFunc(BlindPixelFitFunc 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 : fBlindPixelChargefirst; rmax = (rmax != 0.) ? rmax : fBlindPixelChargelast; // // 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 = fHBlindPixelCharge->Integral("width"); const Double_t lambda_guess = 0.5; const Double_t mu_0_guess = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin()); const Double_t si_0_guess = 20.; const Double_t mu_1_guess = mu_0_guess + 50.; const Double_t si_1_guess = si_0_guess + si_0_guess; const Double_t norm = entries/gkSq2Pi; fSinglePheFit = new TF1("SinglePheFit",fgSinglePheFitFunc,rmin,rmax,fgSinglePheFitNPar); // 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); fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm); // fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1"); fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","Area"); fSinglePheFit->SetParLimits(0,0.,1.); fSinglePheFit->SetParLimits(1,rmin,(rmax-rmin)/1.5); fSinglePheFit->SetParLimits(2,(rmax-rmin)/2.,(rmax-0.05*(rmax-rmin))); fSinglePheFit->SetParLimits(3,1.0,(rmax-rmin)/2.0); fSinglePheFit->SetParLimits(4,1.0,(rmax-rmin)/2.5); fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1); // fSinglePheFit->SetParLimits(5,entries/gkSq2Pi,entries/gkSq2Pi); // fSinglePheFit->SetParLimits(5,0.,1.5*entries); // // Normalize the histogram to facilitate faster fitting of the area // For speed reasons, gfKto4 is normalized to Sqrt(2 pi). // Therefore also normalize the histogram to that value // // fHBlindPixelCharge->Scale(gkSq2Pi*(float)bins/npx/entries); // fHBlindPixelCharge->Scale(gkSq2Pi/entries); // norm *= (Float_t)fSinglePheFit->GetNpx()/(Float_t)fBlindPixelChargenbins; fHBlindPixelCharge->Fit(fSinglePheFit,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(); // Perform the cross-check fitting only the pedestal: fSinglePhePedFit = new TF1("GausPed","gaus",rmin,fMu0); fHBlindPixelCharge->Fit(fSinglePhePedFit,opt); Double_t pedarea = fSinglePhePedFit->GetParameter(0)*gkSq2Pi*fSinglePhePedFit->GetParameter(2); fLambdaCheck = TMath::Log(entries/pedarea); fLambdaCheckErr = fSinglePhePedFit->GetParError(0)/fSinglePhePedFit->GetParameter(0) + fSinglePhePedFit->GetParError(2)/fSinglePhePedFit->GetParameter(2); *fLog << "Results of the Blind Pixel Fit: " << endl; *fLog << "Chisquare: " << fChisquare << endl; *fLog << "DoF: " << fNdf << endl; *fLog << "Probability: " << fProb << endl; // // The fit result is accepted under condition that: // 1) the Probability is greater than gkProbLimit (default 0.001 == 99.7%) // 2) at least 100 events are in the single Photo-electron peak // if (fProb < gkProbLimit) { *fLog << err << "Prob: " << fProb << " is smaller than the allowed value: " << gkProbLimit << endl; fFitOK = kFALSE; return kFALSE; } Float_t contSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries; if (contSinglePhe < 100.) { *fLog << err << "Statistics is too low: Only " << contSinglePhe << " in the first electron peak " << endl; fFitOK = kFALSE; return kFALSE; } else *fLog << GetDescriptor() << ": " << contSinglePhe << " in first electron peak " << endl; fFitOK = kTRUE; return kTRUE; } void MHCalibrationBlindPixel::CutAllEdges() { Int_t nbins = 30; CutEdges(fHBlindPixelCharge,nbins); fBlindPixelChargefirst = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetFirst()); fBlindPixelChargelast = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetLast())+fHBlindPixelCharge->GetBinWidth(0); fBlindPixelChargenbins = nbins; CutEdges(fHBlindPixelChargevsN,0); } Bool_t MHCalibrationBlindPixel::FitTime(Axis_t rmin, Axis_t rmax, Option_t *opt) { if (fTimeGausFit) return kFALSE; rmin = (rmin != 0.) ? rmin : 4.; rmax = (rmax != 0.) ? rmax : 9.; const Stat_t entries = fHBlindPixelTime->Integral(); const Double_t mu_guess = fHBlindPixelTime->GetBinCenter(fHBlindPixelTime->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); fHBlindPixelTime->Fit(fTimeGausFit,opt); rmin = fTimeGausFit->GetParameter(1) - 2.*fTimeGausFit->GetParameter(2); rmax = fTimeGausFit->GetParameter(1) + 2.*fTimeGausFit->GetParameter(2); fTimeGausFit->SetRange(rmin,rmax); fHBlindPixelTime->Fit(fTimeGausFit,opt); fMeanTime = fTimeGausFit->GetParameter(2); fSigmaTime = fTimeGausFit->GetParameter(3); fMeanTimeErr = fTimeGausFit->GetParError(2); fSigmaTimeErr = 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; }