/* ======================================================================== *\ ! ! * ! * 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 #include "MFFT.h" #include "MLog.h" #include "MLogManip.h" ClassImp(MHCalibrationBlindPixel); using namespace std; const Int_t MHCalibrationBlindPixel::fgBlindPixelChargeNbins = 1000; const Int_t MHCalibrationBlindPixel::fgBlindPixelTimeNbins = 22; 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; const Int_t MHCalibrationBlindPixel::fPSDNbins = 50; const Int_t MHCalibrationBlindPixel::fPulserFrequency = 500; // -------------------------------------------------------------------------- // // Default Constructor. // MHCalibrationBlindPixel::MHCalibrationBlindPixel(const char *name, const char *title) : fHBlindPixelPSD(NULL), fSinglePheFit(NULL), fTimeGausFit(NULL), fSinglePhePedFit(NULL), fPSDHiGain(NULL), fPSDLoGain(NULL), fHPSD(NULL), fPSDExpFit(NULL), fChargeXaxis(NULL), fPSDXaxis(NULL), fCurrentSize(1024), 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 = -400.; fBlindPixelChargelast = 1600.; fHBlindPixelCharge = new TH1I("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); fHiGains = new TArrayF(fCurrentSize); fLoGains = new TArrayF(fCurrentSize); fHSinglePheFADCSlices = new TH1I("HSinglePheFADCSlices","FADC slices Single Phe events",30,0.5,30.5); fHSinglePheFADCSlices->SetXTitle("FADC slice"); fHSinglePheFADCSlices->SetYTitle("Counts"); fHSinglePheFADCSlices->SetDirectory(NULL); fHPedestalFADCSlices = new TH1I("HPedestalFADCSlices", "FADC slices Pedestal events",30,0.5,30.5); fHPedestalFADCSlices->SetXTitle("FADC slice"); fHPedestalFADCSlices->SetYTitle("Counts"); fHPedestalFADCSlices->SetDirectory(NULL); Clear(); } MHCalibrationBlindPixel::~MHCalibrationBlindPixel() { if (fFitLegend) delete fFitLegend; delete fHBlindPixelCharge; delete fHBlindPixelTime; delete fHiGains; delete fLoGains; delete fHSinglePheFADCSlices; delete fHPedestalFADCSlices; if (fHBlindPixelPSD) delete fHBlindPixelPSD; if (fSinglePheFit) delete fSinglePheFit; if (fTimeGausFit) delete fTimeGausFit; if(fSinglePhePedFit) delete fSinglePhePedFit; if (fPSDExpFit) delete fPSDExpFit; if (fHPSD) delete fHPSD; if (fChargeXaxis) delete fChargeXaxis; if (fPSDXaxis) delete fPSDXaxis; } void MHCalibrationBlindPixel::Clear(Option_t *o) { fTotalEntries = 0; fCurrentSize = 1024; fBlindPixelChargefirst = -400.; fBlindPixelChargelast = 1600.; 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 = kEPoisson5; if (fFitLegend) delete fFitLegend; if (fHBlindPixelPSD) delete fHBlindPixelPSD; if (fSinglePheFit) delete fSinglePheFit; if (fTimeGausFit) delete fTimeGausFit; if(fSinglePhePedFit) delete fSinglePhePedFit; if (fPSDExpFit) delete fPSDExpFit; if (fHPSD) delete fHPSD; if (fChargeXaxis) delete fChargeXaxis; if (fPSDXaxis) delete fPSDXaxis; if (fPSDHiGain) delete fPSDHiGain; if (fPSDLoGain) delete fPSDLoGain; CLRBIT(fFlags,kFitOK); CLRBIT(fFlags,kOscillating); return; } void MHCalibrationBlindPixel::Reset() { Clear(); fHBlindPixelCharge->Reset(); fHBlindPixelTime->Reset(); fHiGains->Set(1024); fLoGains->Set(1024); fHiGains->Reset(0); fLoGains->Reset(0); } Bool_t MHCalibrationBlindPixel::CheckOscillations() { if (fPSDExpFit) return IsOscillating(); // This cuts only the non-used zero's, but MFFT will later cut the rest CutArrayBorder(fHiGains); CutArrayBorder(fLoGains); MFFT fourier; fPSDLoGain = fourier.PowerSpectrumDensity(fLoGains); fPSDHiGain = fourier.PowerSpectrumDensity(fHiGains); fHPSD = ProjectArray(*fPSDHiGain, fPSDNbins, "PSDProjectionBlindPixel", "Power Spectrum Density Projection HiGain Blind Pixel"); // // First guesses for the fit (should be as close to reality as possible, // const Double_t xmax = fHPSD->GetXaxis()->GetXmax(); fPSDExpFit = new TF1("PSDExpFit","exp([0]-[1]*x)",0.,xmax); const Double_t slope_guess = (TMath::Log(fHPSD->GetEntries())+1.)/xmax; const Double_t offset_guess = slope_guess*xmax; fPSDExpFit->SetParameters(offset_guess, slope_guess); fPSDExpFit->SetParNames("Offset","Slope"); fPSDExpFit->SetParLimits(0,offset_guess/2.,2.*offset_guess); fPSDExpFit->SetParLimits(1,slope_guess/1.5,1.5*slope_guess); fPSDExpFit->SetRange(0.,xmax); fHPSD->Fit(fPSDExpFit,"RQL0"); fPSDProb = fPSDExpFit->GetProb(); if (fPSDProb < gkProbLimit) { SETBIT(fFlags,kOscillating); return kTRUE; } CLRBIT(fFlags,kOscillating); return kFALSE; } void MHCalibrationBlindPixel::CreatePSDXaxis(Int_t n) { if (fPSDXaxis) return; fPSDXaxis = new TArrayF(n); for (Int_t i=0;iAddAt((Float_t)(fPulserFrequency*i)/2./n,i); } void MHCalibrationBlindPixel::CreateChargeXaxis(Int_t n) { if (!fChargeXaxis) { fChargeXaxis = new TArrayF(n); for (Int_t i=0;iAddAt((Float_t)i/fPulserFrequency,i); return; } if (fChargeXaxis->GetSize() == n) return; const Int_t diff = fChargeXaxis->GetSize()-n; fChargeXaxis->Set(n); if (diff < 0) for (Int_t i=n;iAddAt((Float_t)i/fPulserFrequency,i); } void MHCalibrationBlindPixel::CutArrayBorder(TArrayF *array) const { Int_t i; for (i=array->GetSize()-1;i>=0;i--) if (array->At(i) != 0) { array->Set(i+1); break; } } void MHCalibrationBlindPixel::CutArrayBorder(TArrayI *array) const { Int_t i; for (i=array->GetSize()-1;i>=0;i--) if (array->At(i) != 0) { array->Set(i+1); break; } } const Bool_t MHCalibrationBlindPixel::IsFitOK() const { return TESTBIT(fFlags,kFitOK); } const Bool_t MHCalibrationBlindPixel::IsOscillating() { if (fPSDExpFit) return TESTBIT(fFlags,kOscillating); return CheckOscillations(); } Bool_t MHCalibrationBlindPixel::FillGraphs(const Int_t qhi, const Int_t qlo) { if (fTotalEntries >= fCurrentSize) { fCurrentSize *= 2; fHiGains->Set(fCurrentSize); fLoGains->Set(fCurrentSize); } fHiGains->AddAt(qhi,fTotalEntries); fLoGains->AddAt(qlo,fTotalEntries); fTotalEntries++; return kTRUE; } Bool_t MHCalibrationBlindPixel::FillBlindPixelCharge(const Int_t q) { return fHBlindPixelCharge->Fill(q) > -1; } Bool_t MHCalibrationBlindPixel::FillBlindPixelTime(const Float_t t) { return fHBlindPixelTime->Fill(t) > -1; } // ------------------------------------------------------------------------- // // Draw a legend with the fit results // void MHCalibrationBlindPixel::DrawLegend() { fFitLegend = new TPaveText(0.05,0.05,0.95,0.95); if (IsFitOK()) 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 (IsFitOK()) { 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,5); gROOT->SetSelectedPad(NULL); c->cd(1); gPad->SetLogx(0); gPad->SetLogy(1); gPad->SetTicks(); fHBlindPixelCharge->Draw(opt); if (IsFitOK()) 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->cd(3); gPad->SetLogy(0); gPad->SetBorderMode(0); // fHBlindPixelTime->Draw(opt); fHSinglePheFADCSlices->Draw(opt); c->cd(4); gPad->SetLogy(0); gPad->SetBorderMode(0); // fHBlindPixelTime->Draw(opt); fHPedestalFADCSlices->Draw(opt); c->Modified(); c->Update(); CutArrayBorder(fHiGains); CreateChargeXaxis(fHiGains->GetSize()); c->cd(5); gPad->SetTicks(); TGraph *gr1 = new TGraph(fChargeXaxis->GetSize(), fChargeXaxis->GetArray(), fHiGains->GetArray()); gr1->SetTitle("Evolution of HiGain charges with time"); gr1->GetXaxis()->SetTitle("Time [s]"); gr1->GetYaxis()->SetTitle("Sum FADC slices"); gr1->SetBit(kCanDelete); gr1->Draw("AL"); CutArrayBorder(fLoGains); CreateChargeXaxis(fHiGains->GetSize()); c->cd(6); gPad->SetTicks(); TGraph *gr2 = new TGraph(fChargeXaxis->GetSize(), fChargeXaxis->GetArray(), fLoGains->GetArray()); gr2->SetTitle("Evolution of LoGain charges with time"); gr2->GetXaxis()->SetTitle("Time [s]"); gr2->GetYaxis()->SetTitle("Sum FADC slices"); gr2->SetBit(kCanDelete); gr2->Draw("AL"); c->Modified(); c->Update(); c->cd(7); TArrayF *array; if (!fPSDHiGain) return; array = fPSDHiGain; if (!fPSDXaxis) CreatePSDXaxis(array->GetSize()); TGraph *gr3 = new TGraph(fPSDXaxis->GetSize(),fPSDXaxis->GetArray(),array->GetArray()); gr3->SetTitle("Power Spectrum Density"); gr3->GetXaxis()->SetTitle("Freq. [Hz]"); gr3->GetYaxis()->SetTitle("P(f)"); gr3->SetBit(kCanDelete); gr3->Draw("AL"); c->Modified(); c->Update(); c->cd(8); gStyle->SetOptStat(111111); gPad->SetTicks(); if (fHPSD->Integral() > 0) gPad->SetLogy(); fHPSD->Draw(opt); fHPSD->GetXaxis()->SetTitle("P(f)"); fHPSD->GetYaxis()->SetTitle("Counts"); if (fPSDExpFit) { fPSDExpFit->SetLineColor(IsOscillating() ? kRed : kGreen); fPSDExpFit->Draw("same"); } c->cd(9); array = fPSDLoGain; if (!fPSDXaxis) CreatePSDXaxis(array->GetSize()); TGraph *gr4 = new TGraph(fPSDXaxis->GetSize(),fPSDXaxis->GetArray(),array->GetArray()); gr4->SetTitle("Power Spectrum Density"); gr4->GetXaxis()->SetTitle("Freq. [Hz]"); gr4->GetYaxis()->SetTitle("P(f)"); gr4->SetBit(kCanDelete); gr4->Draw("AL"); c->Modified(); c->Update(); c->cd(10); } 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) { *fLog << inf << "min: " << min << endl; *fLog << "max: " << max << endl; // // 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.1; 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; *fLog << inf << "pedestal: " << fMeanPedestal << endl; *fLog << "sigma: " << fSigmaPedestal << endl; // // 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.,0.5); if ((fMeanPedestal) && (fSigmaPedestal)) fSinglePheFit->SetParLimits(1, fMeanPedestal-5.*fMeanPedestalErr, fMeanPedestal+5.*fMeanPedestalErr); else fSinglePheFit->SetParLimits(1,-3.,0.); fSinglePheFit->SetParLimits(2,min,max); if ((fMeanPedestal) && (fSigmaPedestal)) fSinglePheFit->SetParLimits(3, fSigmaPedestal-5.*fSigmaPedestalErr, fSigmaPedestal+5.*fSigmaPedestalErr); else fSinglePheFit->SetParLimits(3,1.0,(max-min)/2.0); fSinglePheFit->SetParLimits(4,0.,(max-min)); fSinglePheFit->SetParLimits(5,norm-(0.5*norm),norm+(0.5*norm)); 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 << "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 50 events are in the single Photo-electron peak // if (fProb < gkProbLimit) { *fLog << warn << "WARNING - Fit Probability " << fProb << " is smaller than the allowed value: " << gkProbLimit << endl; CLRBIT(fFlags,kFitOK); return kFALSE; } Float_t contSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries; if (contSinglePhe < 50.) { *fLog << warn << "WARNING - Statistics is too low: Only " << contSinglePhe << " in the Single Photo-Electron peak " << endl; CLRBIT(fFlags,kFitOK); return kFALSE; } else *fLog << inf << contSinglePhe << " in Single Photo-Electron peak " << endl; SETBIT(fFlags,kFitOK); return kTRUE; } void MHCalibrationBlindPixel::CutAllEdges() { Int_t nbins = 25; StripZeros(fHBlindPixelCharge,nbins); fBlindPixelChargefirst = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetFirst()); fBlindPixelChargelast = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetLast())+fHBlindPixelCharge->GetBinWidth(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 << "Chisquare: " << fTimeGausFit->GetChisquare() << endl; *fLog << "Ndf: " << fTimeGausFit->GetNDF() << endl; return kTRUE; }