/* ======================================================================== *\ ! ! * ! * 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 // // Used numbers are the following: // // Electronic conversion factor: // Assume, we have N_e electrons at the anode, // thus a charge of N_e*e (e = electron charge) Coulomb. // // This charge is AC coupled and runs into a R_pre = 50 Ohm resistency. // The corresponding current is amplified by a gain factor G_pre = 400 // (the precision of this value still has to be checked !!!) and again AC coupled to // the output. // The corresponding signal goes through the whole transmission and // amplification chain and is digitized in the FADCs. // The conversion Signal Area to FADC counts (Conv_trans) has been measured // by David and Oscar to be approx. 3.9 pVs^-1 // // Thus: Conversion FADC counts to Number of Electrons at Anode: // FADC counts = (1/Conv_tran) * G_pre * R_pre * e * N_e = 8 * 10^-4 N_e. // // Also: FADC counts = 8*10^-4 * GAIN * N_phe // // In the blind pixel, there is an additional pre-amplifier with an amplification of // about 10. Therefore, we have for the blind pixel: // // FADC counts (Blind Pixel) = 8*10^-3 * GAIN * N_phe // ////////////////////////////////////////////////////////////////////////////// #include "MHCalibrationBlindPixel.h" #include "MHCalibrationConfig.h" #include #include #include #include #include #include #include "MLog.h" #include "MLogManip.h" ClassImp(MHCalibrationBlindPixel); using namespace std; const Int_t MHCalibrationBlindPixel::fgBlindPixelChargeNbins = 1000; const Int_t MHCalibrationBlindPixel::fgBlindPixelTimeNbins = 22; const Int_t MHCalibrationBlindPixel::fgBlindPixelChargevsNbins = 10000; const Axis_t MHCalibrationBlindPixel::fgBlindPixelTimeFirst = -9.00; const Axis_t MHCalibrationBlindPixel::fgBlindPixelTimeLast = 12.00; const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmp = 0.008; const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmpError = 0.002; // -------------------------------------------------------------------------- // // Default Constructor. // MHCalibrationBlindPixel::MHCalibrationBlindPixel(const char *name, const char *title) : fHBlindPixelPSD(NULL), fSinglePheFit(NULL), fTimeGausFit(NULL), fSinglePhePedFit(NULL), fFitLegend(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 fBlindPixelChargefirst = -200.; fBlindPixelChargelast = 800.; fHBlindPixelCharge = new TH1F("HBlindPixelCharge","Distribution of Summed FADC Slices", fgBlindPixelChargeNbins,fBlindPixelChargefirst,fBlindPixelChargelast); fHBlindPixelCharge->SetXTitle("Sum FADC Slices"); fHBlindPixelCharge->SetYTitle("Nr. of events"); fHBlindPixelCharge->Sumw2(); fHBlindPixelCharge->SetDirectory(NULL); fHBlindPixelTime = new TH1F("HBlindPixelTime","Distribution of Mean Arrival Times", fgBlindPixelTimeNbins,fgBlindPixelTimeFirst,fgBlindPixelTimeLast); fHBlindPixelTime->SetXTitle("Mean Arrival Times [FADC slice nr]"); fHBlindPixelTime->SetYTitle("Nr. of events"); fHBlindPixelTime->Sumw2(); fHBlindPixelTime->SetDirectory(NULL); fHBlindPixelChargevsN = new TH1I("HBlindPixelChargevsN","Sum of Charges vs. Event Number", fgBlindPixelChargevsNbins,-0.5,(Axis_t)fgBlindPixelChargevsNbins-0.5); fHBlindPixelChargevsN->SetXTitle("Event Nr."); fHBlindPixelChargevsN->SetYTitle("Sum of FADC slices"); fHBlindPixelChargevsN->SetDirectory(NULL); Clear(); } 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; } void MHCalibrationBlindPixel::Clear(Option_t *o) { fBlindPixelChargefirst = -200.; fBlindPixelChargelast = 800.; fLambda = 0.; fMu0 = 0.; fMu1 = 0.; fSigma0 = 0.; fSigma1 = 0.; fLambdaErr = 0.; fMu0Err = 0.; fMu1Err = 0.; fSigma0Err = 0.; fSigma1Err = 0.; fChisquare = -1.; fProb = -1.; fNdf = -1; fMeanTime = -1.; fMeanTimeErr = -1.; fSigmaTime = -1.; fSigmaTimeErr = -1.; fLambdaCheck = -1.; fLambdaCheckErr = -1.; fMeanPedestal = 0.; fMeanPedestalErr = 0.; fSigmaPedestal = 0.; fSigmaPedestalErr = 0.; fFitFunc = kEPoisson4; if (fFitLegend) delete fFitLegend; if (fHBlindPixelPSD) delete fHBlindPixelPSD; if (fSinglePheFit) delete fSinglePheFit; if (fTimeGausFit) delete fTimeGausFit; if(fSinglePhePedFit) delete fSinglePhePedFit; return; } void MHCalibrationBlindPixel::Reset() { Clear(); fHBlindPixelCharge->Reset(); fHBlindPixelTime->Reset(); fHBlindPixelChargevsN->Reset(); } Bool_t MHCalibrationBlindPixel::FillBlindPixelCharge(Float_t q) { return fHBlindPixelCharge->Fill(q) > -1; } Bool_t MHCalibrationBlindPixel::FillBlindPixelTime(Float_t t) { return fHBlindPixelTime->Fill(t) > -1; } Bool_t MHCalibrationBlindPixel::FillBlindPixelChargevsN(Stat_t rq, Int_t t) { return fHBlindPixelChargevsN->Fill(t,rq) > -1; } // ------------------------------------------------------------------------- // // 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(1); 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 maximum_bin = fHBlindPixelCharge->GetBinCenter(fHBlindPixelCharge->GetMaximumBin()); const Double_t norm = entries/gkSq2Pi; // // Initialize the fit function // switch (fFitFunc) { case kEPoisson4: fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,min,max,6); break; case kEPoisson5: fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,min,max,6); break; case kEPoisson6: fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,min,max,6); break; case kEPolya: fSinglePheFit = new TF1("SinglePheFit",&fPolya,min,max,8); break; case kEMichele: break; default: *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl; return kFALSE; break; } const Double_t mu_0_guess = maximum_bin; const Double_t si_0_guess = 40.; const Double_t mu_1_guess = mu_0_guess + 100.; const Double_t si_1_guess = si_0_guess + si_0_guess; // Michele const Double_t lambda_1cat_guess = 0.5; const Double_t lambda_1dyn_guess = 0.5; const Double_t mu_1cat_guess = mu_0_guess + 50.; const Double_t mu_1dyn_guess = mu_0_guess + 20.; const Double_t si_1cat_guess = si_0_guess + si_0_guess; const Double_t si_1dyn_guess = si_0_guess; // Polya const Double_t excessPoisson_guess = 0.5; const Double_t delta1_guess = 8.; const Double_t delta2_guess = 5.; const Double_t electronicAmp_guess = fgBlindPixelElectronicAmp; const Double_t electronicAmp_limit = fgBlindPixelElectronicAmpError; // // Initialize boundaries and start parameters // switch (fFitFunc) { case kEPoisson4: if ((fMeanPedestal) && (fSigmaPedestal)) fSinglePheFit->SetParameters(lambda_guess,fMeanPedestal,mu_1_guess,fSigmaPedestal,si_1_guess,norm); else 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.); if ((fMeanPedestal) && (fSigmaPedestal)) fSinglePheFit->SetParLimits(1, fMeanPedestal-1.*fMeanPedestalErr, fMeanPedestal+1.*fMeanPedestalErr); else fSinglePheFit->SetParLimits(1,-3.,0.); fSinglePheFit->SetParLimits(2,(max-min)/2.,max); if ((fMeanPedestal) && (fSigmaPedestal)) fSinglePheFit->SetParLimits(3, fSigmaPedestal-3.*fSigmaPedestalErr, fSigmaPedestal+3.*fSigmaPedestalErr); else fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0); fSinglePheFit->SetParLimits(4,1.0,(max-min)); fSinglePheFit->SetParLimits(5,norm-0.5,norm+0.5); break; case kEPoisson5: case kEPoisson6: 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 kEPolya: if ((fMeanPedestal) && (fSigmaPedestal)) fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess, delta1_guess,delta2_guess, electronicAmp_guess, fSigmaPedestal, norm, fMeanPedestal); else fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess, delta1_guess,delta2_guess, electronicAmp_guess, si_0_guess, norm, mu_0_guess); fSinglePheFit->SetParNames("#lambda","b_{tot}", "#delta_{1}","#delta_{2}", "amp_{e}","#sigma_{0}", "Area", "#mu_{0}"); fSinglePheFit->SetParLimits(0,0.,1.); fSinglePheFit->SetParLimits(1,0.,1.); fSinglePheFit->SetParLimits(2,6.,12.); fSinglePheFit->SetParLimits(3,3.,8.); fSinglePheFit->SetParLimits(4,electronicAmp_guess-electronicAmp_limit, electronicAmp_guess+electronicAmp_limit); if ((fMeanPedestal) && (fSigmaPedestal)) fSinglePheFit->SetParLimits(5, fSigmaPedestal-3.*fSigmaPedestalErr, fSigmaPedestal+3.*fSigmaPedestalErr); else fSinglePheFit->SetParLimits(5,min,(max-min)/1.5); fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1); if ((fMeanPedestal) && (fSigmaPedestal)) fSinglePheFit->SetParLimits(7, fMeanPedestal-3.*fMeanPedestalErr, fMeanPedestal+3.*fMeanPedestalErr); else fSinglePheFit->SetParLimits(7,-35.,15.); 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; case kEPolya: fLambda = fSinglePheFit->GetParameter(0); fMu0 = fSinglePheFit->GetParameter(7); fMu1 = 0.; fSigma0 = fSinglePheFit->GetParameter(5); fSigma1 = 0.; fLambdaErr = fSinglePheFit->GetParError(0); fMu0Err = fSinglePheFit->GetParError(7); fMu1Err = 0.; fSigma0Err = fSinglePheFit->GetParError(5); fSigma1Err = 0.; 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,0.); 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 50 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; } Float_t contSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries; if (contSinglePhe < 50.) { *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); 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; }