/* ======================================================================== *\ ! ! * ! * 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 "MLog.h" #include "MLogManip.h" ClassImp(MHCalibrationBlindPixel); using namespace std; // -------------------------------------------------------------------------- // // Default Constructor. // MHCalibrationBlindPixel::MHCalibrationBlindPixel(const char *name, const char *title) : fHBlindPixelCharge(NULL), fHBlindPixelTime(NULL), fHBlindPixelChargevsN(NULL), fHBlindPixelPSD(NULL), fSinglePheFit(NULL), fTimeGausFit(NULL), fSinglePhePedFit(NULL), fBlindPixelChargefirst(0.), fBlindPixelChargelast(0.), fBlindPixelChargenbins(0), fFitLegend(NULL), fFitOK(kFALSE), 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.), fFitFunc(kEPoisson4) { 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() { if (fFitLegend) delete fFitLegend; delete fHBlindPixelCharge; delete fHBlindPixelTime; delete fHBlindPixelChargevsN; if (fHBlindPixelPSD) delete fHBlindPixelPSD; if (fSinglePheFit) delete fSinglePheFit; if (fTimeGausFit) delete fTimeGausFit; if(fSinglePhePedFit) delete fSinglePhePedFit; } Bool_t MHCalibrationBlindPixel::FillBlindPixelCharge(Float_t q) { return fHBlindPixelCharge->Fill(q) > -1; } Bool_t MHCalibrationBlindPixel::FillBlindPixelTime(Int_t t) { return fHBlindPixelTime->Fill(t) > -1; } Bool_t MHCalibrationBlindPixel::FillBlindPixelChargevsN(Stat_t rq, Int_t t) { return fHBlindPixelChargevsN->Fill(t,rq) > -1; } 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()); TText *t1 = fFitLegend->AddText(line1.Data()); t1->SetBit(kCanDelete); const TString line6 = Form("Mean #lambda (check) = %2.2f #pm %2.2f",GetLambdaCheck(),GetLambdaCheckErr()); TText *t2 = fFitLegend->AddText(line6.Data()); t2->SetBit(kCanDelete); const TString line2 = Form("Pedestal: #mu_{0} = %2.2f #pm %2.2f",GetMu0(),GetMu0Err()); TText *t3 = fFitLegend->AddText(line2.Data()); t3->SetBit(kCanDelete); const TString line3 = Form("Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",GetSigma0(),GetSigma0Err()); TText *t4 = fFitLegend->AddText(line3.Data()); t4->SetBit(kCanDelete); const TString line4 = Form("1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",GetMu1(),GetMu1Err()); TText *t5 = fFitLegend->AddText(line4.Data()); t5->SetBit(kCanDelete); const TString line5 = Form("Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",GetSigma1(),GetSigma1Err()); TText *t6 = fFitLegend->AddText(line5.Data()); t6->SetBit(kCanDelete); const TString line7 = Form("#chi^{2} / N_{dof}: %4.2f / %3i",GetChiSquare(),GetNdf()); TText *t7 = fFitLegend->AddText(line7.Data()); t7->SetBit(kCanDelete); const TString line8 = Form("Probability: %4.2f ",GetProb()); TText *t8 = fFitLegend->AddText(line8.Data()); t8->SetBit(kCanDelete); if (fFitOK) { TText *t = fFitLegend->AddText(0.,0.,"Result of the Fit: OK"); t->SetBit(kCanDelete); } else { TText *t = fFitLegend->AddText("Result of the Fit: NOT OK"); t->SetBit(kCanDelete); } fFitLegend->SetBit(kCanDelete); fFitLegend->Draw(); return; } TObject *MHCalibrationBlindPixel::DrawClone(Option_t *option) const { gROOT->SetSelectedPad(NULL); MHCalibrationBlindPixel *newobj = (MHCalibrationBlindPixel*)Clone(); if (!newobj) return 0; newobj->SetBit(kCanDelete); if (strlen(option)) newobj->Draw(option); else newobj->Draw(GetDrawOption()); return newobj; } // ------------------------------------------------------------------------- // // Draw the histogram // void MHCalibrationBlindPixel::Draw(Option_t *opt) { gStyle->SetOptFit(0); gStyle->SetOptStat(111111); TCanvas *c = MakeDefCanvas(this,550,700); c->Divide(2,3); gROOT->SetSelectedPad(NULL); c->cd(1); gPad->SetLogx(0); gPad->SetLogy(1); gPad->SetTicks(); fHBlindPixelCharge->Draw(opt); if (fFitOK) fSinglePheFit->SetLineColor(kGreen); else fSinglePheFit->SetLineColor(kRed); fSinglePheFit->Draw("same"); c->Modified(); c->Update(); fSinglePhePedFit->SetLineColor(kBlue); fSinglePhePedFit->Draw("same"); c->cd(2); DrawLegend(); c->Modified(); c->Update(); c->cd(3); gPad->SetLogy(1); gPad->SetBorderMode(0); fHBlindPixelTime->Draw(opt); c->cd(4); fHBlindPixelChargevsN->Draw(opt); c->Modified(); c->Update(); c->cd(5); // fHBlindPixelPSD = (TH1F*)MFFT::PowerSpectrumDensity(fHBlindPixelChargevsN); // TH1F *hist = MFFT::PowerSpectrumDensity((*this)->fHBlindPixelChargevsN); // fHBlindPixelPSD->Draw(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; } if (!InitFit(fBlindPixelChargefirst,fBlindPixelChargelast)) return kFALSE; for (Int_t i=0;i<10000; i++) fHBlindPixelCharge->Fill(fSinglePheFit->GetRandom()); return kTRUE; } Bool_t MHCalibrationBlindPixel::InitFit(Axis_t min, Axis_t max) { // // 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; // // Initialize // switch (fFitFunc) { case kEPoisson4: fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,min,max,6); 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","Area"); fSinglePheFit->SetParLimits(0,0.,1.); fSinglePheFit->SetParLimits(1,min,(max-min)/1.5); fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min))); fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0); fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5); fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1); break; case kEPoisson5: fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,min,max,6); 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","Area"); fSinglePheFit->SetParLimits(0,0.,1.); fSinglePheFit->SetParLimits(1,min,(max-min)/1.5); fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min))); fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0); fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5); fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1); break; case kEPoisson6: fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,min,max,6); 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","Area"); fSinglePheFit->SetParLimits(0,0.,1.); fSinglePheFit->SetParLimits(1,min,(max-min)/1.5); fSinglePheFit->SetParLimits(2,(max-min)/2.,(max-0.05*(max-min))); fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0); fSinglePheFit->SetParLimits(4,1.0,(max-min)/2.5); fSinglePheFit->SetParLimits(5,norm-0.1,norm+0.1); break; case kEPoisson7: break; case kEPolya: break; case kEMichele: break; default: *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl; return kFALSE; break; } return kTRUE; } void MHCalibrationBlindPixel::ExitFit(TF1 *f) { // // Finalize // switch (fFitFunc) { case kEPoisson4: case kEPoisson5: case kEPoisson6: case kEPoisson7: 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); break; default: break; } } Bool_t MHCalibrationBlindPixel::FitSinglePhe(Axis_t rmin, Axis_t rmax, Option_t *opt) { // // Get the fitting ranges // rmin = (rmin != 0.) ? rmin : fBlindPixelChargefirst; rmax = (rmax != 0.) ? rmax : fBlindPixelChargelast; if (!InitFit(rmin,rmax)) return kFALSE; fHBlindPixelCharge->Fit(fSinglePheFit,opt); ExitFit(fSinglePheFit); 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); const Stat_t entries = fHBlindPixelCharge->Integral("width"); 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 << inf << "Results of the Blind Pixel Fit: " << endl; *fLog << inf << "Chisquare: " << fChisquare << endl; *fLog << inf << "DoF: " << fNdf << endl; *fLog << inf << "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 << "ERROR: Fit Probability " << fProb << " is smaller than the allowed value: " << gkProbLimit << endl; fFitOK = kFALSE; return kFALSE; } if (fProb < 0.01) *fLog << warn << "WARNING: Fit Probability " << fProb << " is smaller than 1% " << endl; Float_t contSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries; if (contSinglePhe < 100.) { *fLog << err << "ERROR: Statistics is too low: Only " << contSinglePhe << " in the Single Photo-Electron peak " << endl; fFitOK = kFALSE; return kFALSE; } else *fLog << inf << contSinglePhe << " in Single Photo-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) { 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); *fLog << inf << "Results of the Times Fit: " << endl; *fLog << inf << "Chisquare: " << fTimeGausFit->GetChisquare() << endl; *fLog << inf << "Ndf: " << fTimeGausFit->GetNDF() << endl; return kTRUE; }