Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc (revision 3289) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc (revision 3289) @@ -0,0 +1,850 @@ +/* ======================================================================== *\ +! +! * +! * 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 02/2004 +! +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MHCalibrationChargeBlindPix +// +// 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 "MHCalibrationChargeBlindPix.h" + +#include +#include +#include + +#include +#include +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MParList.h" + +#include "MRawEvtData.h" +#include "MRawEvtPixelIter.h" + +#include "MExtractedSignalBlindPixel.h" +#include "MCalibrationChargeBlindPix.h" + +ClassImp(MHCalibrationChargeBlindPix); + +using namespace std; + +const Double_t MHCalibrationChargeBlindPix::gkElectronicAmp = 0.008; +const Double_t MHCalibrationChargeBlindPix::gkElectronicAmpErr = 0.002; + +const Int_t MHCalibrationChargeBlindPix::fgChargeNbins = 1400; +const Axis_t MHCalibrationChargeBlindPix::fgChargeFirst = -200.5; +const Axis_t MHCalibrationChargeBlindPix::fgChargeLast = 1199.5; + +const Float_t MHCalibrationChargeBlindPix::fgSinglePheCut = 200.; +const Float_t MHCalibrationChargeBlindPix::fgNumSinglePheLimit = 50.; +// -------------------------------------------------------------------------- +// +// Default Constructor. +// +MHCalibrationChargeBlindPix::MHCalibrationChargeBlindPix(const char *name, const char *title) + : fBlindPix(NULL), fSignal(NULL), fRawEvt(NULL), + fASinglePheFADCSlices(30), fAPedestalFADCSlices(30), + fSinglePheFit(NULL), + fFitLegend(NULL), + fHSinglePheFADCSlices(NULL), fHPedestalFADCSlices(NULL) +{ + + fName = name ? name : "MHCalibrationChargeBlindPix"; + fTitle = title ? title : "Fill the accumulated charges and times of all Blind Pixel events and perform fits"; + + SetChargeNbins(); + SetChargeFirst(); + SetChargeLast(); + + SetSinglePheCut(); + SetNumSinglePheLimit(); + + fHGausHist.SetName("HCalibrationChargeBlindPix"); + fHGausHist.SetTitle("Distribution of Summed FADC slices Blind Pixel"); + fHGausHist.SetXTitle("Sum FADC Slices"); + fHGausHist.SetYTitle("Nr. of events"); + + Clear(); +} + +MHCalibrationChargeBlindPix::~MHCalibrationChargeBlindPix() +{ + + if (fFitLegend) + delete fFitLegend; + + if (fSinglePheFit) + delete fSinglePheFit; + + if (fHSinglePheFADCSlices) + delete fHSinglePheFADCSlices; + + if (fHPedestalFADCSlices) + delete fHPedestalFADCSlices; + +} +void MHCalibrationChargeBlindPix::Init() +{ + + fHGausHist.SetBins( fChargeNbins, fChargeFirst, fChargeLast); +} + +void MHCalibrationChargeBlindPix::Clear(Option_t *o) +{ + + fLambda = -999.; + fMu0 = -999.; + fMu1 = -999.; + fSigma0 = -999.; + fSigma1 = -999.; + + fLambdaErr = -999.; + fMu0Err = -999.; + fMu1Err = -999.; + fSigma0Err = -999.; + fSigma1Err = -999.; + + fLambdaCheck = -999.; + fLambdaCheckErr = -999.; + + fMeanPedestal = 0.; + fMeanPedestalErr = 0.; + fSigmaPedestal = 0.; + fSigmaPedestalErr = 0.; + + fFitFunc = kEPoisson5; + + fExtractSlices = 0; + fNumSinglePhes = 0; + fNumPedestals = 0; + + fNumSinglePhes = 0; + fNumPedestals = 0; + + fChisquare = 0.; + fNDF = 0 ; + fProb = 0.; + + SetSinglePheFitOK ( kFALSE ); + SetPedestalFitOK ( kFALSE ); + + if (fFitLegend) + delete fFitLegend; + + if (fSinglePheFit) + delete fSinglePheFit; + + MHCalibrationChargePix::Clear(); + return; +} + +void MHCalibrationChargeBlindPix::SetSinglePheFitOK ( const Bool_t b=kTRUE) +{ + b ? SETBIT(fFlags,kSinglePheFitOK) : CLRBIT(fFlags,kSinglePheFitOK); +} + +void MHCalibrationChargeBlindPix::SetPedestalFitOK( const Bool_t b=kTRUE) +{ + b ? SETBIT(fFlags,kPedestalFitOK) : CLRBIT(fFlags,kPedestalFitOK); +} + +const Bool_t MHCalibrationChargeBlindPix::IsSinglePheFitOK() const +{ + return TESTBIT(fFlags,kSinglePheFitOK); +} + +const Bool_t MHCalibrationChargeBlindPix::IsPedestalFitOK() const +{ + return TESTBIT(fFlags,kPedestalFitOK); +} + +Bool_t MHCalibrationChargeBlindPix::SetupFill(const MParList *pList) +{ + fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); + if (!fRawEvt) + { + *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl; + return kFALSE; + } + + fSignal = (MExtractedSignalBlindPixel*)pList->FindObject("MExtractedSignalBlindPixel"); + if (!fSignal) + { + *fLog << err << GetDescriptor() << ": ERROR: Could not find MExtractedSignalBlindPixel ... aborting " << endl; + return kFALSE; + } + + Init(); + + return kTRUE; +} + +Bool_t MHCalibrationChargeBlindPix::ReInit(MParList *pList) +{ + + fBlindPix = (MCalibrationChargeBlindPix*)pList->FindCreateObj("MCalibrationChargeBlindPix"); + if (!fBlindPix) + { + *fLog << err << GetDescriptor() << ": ERROR: Could not find MCalibrationChargeBlindPix ... aborting " << endl; + return kFALSE; + } + + return kTRUE; +} + +Bool_t MHCalibrationChargeBlindPix::Fill(const MParContainer *par, const Stat_t w) +{ + + Float_t slices = (Float_t)fSignal->GetNumFADCSamples(); + + if (slices == 0.) + { + *fLog << err << "Number of used signal slices in MExtractedSignalBlindPix is zero ... abort." + << endl; + return kFALSE; + } + + if (fExtractSlices != 0. && slices != fExtractSlices ) + { + *fLog << err << "Number of used signal slices changed in MExtractedSignalCam ... abort." + << endl; + return kFALSE; + } + fExtractSlices = slices; + + // + // Signal extraction and histogram filling + // + const Float_t signal = (Float_t)fSignal->GetExtractedSignal(); + FillHistAndArray(signal); + + // + // IN order to study the single-phe posistion, we extract the slices + // + MRawEvtPixelIter pixel(fRawEvt); + pixel.Jump(fSignal->GetBlindPixelIdx()); + + if (signal > fSinglePheCut) + FillSinglePheFADCSlices(pixel); + else + FillPedestalFADCSlices(pixel); + + return kTRUE; +} + +Bool_t MHCalibrationChargeBlindPix::Finalize() +{ + + if (IsEmpty()) + { + *fLog << err << GetDescriptor() << ": My histogram has not been filled !! " << endl; + return kFALSE; + } + + CreateFourierSpectrum(); + fBlindPix->SetOscillating ( !IsFourierSpectrumOK() ); + + fMeanPedestal = fSignal->GetPed(); + fMeanPedestalErr = fSignal->GetPedErr(); + fSigmaPedestal = fSignal->GetPedRms(); + fSigmaPedestalErr = fSignal->GetPedRmsErr(); + + if (fNumSinglePhes > 1) + for (Int_t i=0;i 1) + for (Int_t i=0;iSetFitted(kTRUE); + + fBlindPix->SetLambda ( fLambda ); + fBlindPix->SetMu0 ( fMu0 ); + fBlindPix->SetMu0Err ( fMu0Err ); + fBlindPix->SetMu1 ( fMu1 ); + fBlindPix->SetMu1Err ( fMu1Err ); + fBlindPix->SetSigma0 ( fSigma0 ); + fBlindPix->SetSigma0Err ( fSigma0Err ); + fBlindPix->SetSigma1 ( fSigma1 ); + fBlindPix->SetSigma1Err ( fSigma1Err ); + fBlindPix->SetProb ( fProb ); + + fBlindPix->SetLambdaCheck ( fLambdaCheck ); + fBlindPix->SetLambdaCheckErr ( fLambdaCheckErr ); + + return kTRUE; +} + +void MHCalibrationChargeBlindPix::FillSinglePheFADCSlices(const MRawEvtPixelIter &iter) +{ + + const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples(); + + if (fASinglePheFADCSlices.GetNrows() < n) + fASinglePheFADCSlices.ResizeTo(n); + + Int_t i=0; + + Byte_t *start = iter.GetHiGainSamples(); + Byte_t *end = start + iter.GetNumHiGainSamples(); + + for (Byte_t *ptr = start; ptr < end; ptr++, i++) + fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr; + + start = iter.GetLoGainSamples(); + end = start + iter.GetNumLoGainSamples(); + + for (Byte_t *ptr = start; ptr < end; ptr++, i++) + fASinglePheFADCSlices(i) = fASinglePheFADCSlices(i) + (Float_t)*ptr; + + fNumSinglePhes++; +} + +void MHCalibrationChargeBlindPix::FillPedestalFADCSlices(const MRawEvtPixelIter &iter) +{ + + const Int_t n = iter.GetNumHiGainSamples() + iter.GetNumLoGainSamples(); + + if (fAPedestalFADCSlices.GetNrows() < n) + fAPedestalFADCSlices.ResizeTo(n); + + Int_t i = 0; + Byte_t *start = iter.GetHiGainSamples(); + Byte_t *end = start + iter.GetNumHiGainSamples(); + + for (Byte_t *ptr = start; ptr < end; ptr++, i++) + fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr; + + start = iter.GetLoGainSamples(); + end = start + iter.GetNumLoGainSamples(); + + for (Byte_t *ptr = start; ptr < end; ptr++, i++) + fAPedestalFADCSlices(i) = fAPedestalFADCSlices(i)+ (Float_t)*ptr; + + fNumPedestals++; +} + + + +Bool_t MHCalibrationChargeBlindPix::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1) +{ + + gRandom->SetSeed(); + + if (fHGausHist.GetIntegral() != 0) + { + *fLog << err << "Histogram " << fHGausHist.GetTitle() << " is already filled. " << endl; + *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl; + return kFALSE; + } + + if (!InitFit()) + return kFALSE; + + for (Int_t i=0;i<10000; i++) + fHGausHist.Fill(fSinglePheFit->GetRandom()); + + return kTRUE; +} + +Bool_t MHCalibrationChargeBlindPix::InitFit() +{ + + // + // Get the fitting ranges + // + Axis_t rmin = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetFirst()); + Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetXaxis()->GetLast()); + + if (rmin < 0.) + rmin = 0.; + + // + // 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 = fHGausHist.Integral("width"); + const Double_t lambda_guess = 0.1; + const Double_t maximum_bin = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin()); + const Double_t norm = entries/TMath::Sqrt(TMath::TwoPi()); + + // + // Initialize the fit function + // + switch (fFitFunc) + { + case kEPoisson4: + fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto4,rmin,rmax,6); + break; + case kEPoisson5: + fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,rmin,rmax,6); + break; + case kEPoisson6: + fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto6,rmin,rmax,6); + break; + case kEPolya: + fSinglePheFit = new TF1("SinglePheFit",&fPolya,rmin,rmax,8); + break; + case kEMichele: + break; + + default: + *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl; + return kFALSE; + break; + } + + if (!fSinglePheFit) + { + *fLog << warn << dbginf << "WARNING: Could not create fit function for Single Phe fit" << endl; + return kFALSE; + } + + 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 = gkElectronicAmp; + const Double_t electronicAmp_limit = gkElectronicAmpErr; + + // + // Initialize boundaries and start parameters + // + switch (fFitFunc) + { + + case kEPoisson4: + fSinglePheFit->SetParNames( "#lambda", "#mu_{0}", "#mu_{1}", "#sigma_{0}", "#sigma_{1}","Area"); + fSinglePheFit->SetParameters(lambda_guess,fMeanPedestal,mu_1_guess,fSigmaPedestal,si_1_guess,norm); + + fSinglePheFit->SetParLimits(0,0.,0.5); + fSinglePheFit->SetParLimits(1, + fMeanPedestal-5.*fMeanPedestalErr, + fMeanPedestal+5.*fMeanPedestalErr); + fSinglePheFit->SetParLimits(2,rmin,rmax); + fSinglePheFit->SetParLimits(3, + fSigmaPedestal-5.*fSigmaPedestalErr, + fSigmaPedestal+5.*fSigmaPedestalErr); + fSinglePheFit->SetParLimits(4,0.,(rmax-rmin)); + 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,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); + break; + + case kEPolya: + fSinglePheFit->SetParameters(lambda_guess, excessPoisson_guess, + delta1_guess,delta2_guess, + electronicAmp_guess, + fSigmaPedestal, + norm, + fMeanPedestal); + 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); + fSinglePheFit->SetParLimits(5, + fSigmaPedestal-3.*fSigmaPedestalErr, + fSigmaPedestal+3.*fSigmaPedestalErr); + fSinglePheFit->SetParLimits(6,norm-0.1,norm+0.1); + fSinglePheFit->SetParLimits(7, + fMeanPedestal-3.*fMeanPedestalErr, + fMeanPedestal+3.*fMeanPedestalErr); + break; + case kEMichele: + break; + + default: + *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl; + return kFALSE; + break; + } + + fSinglePheFit->SetRange(rmin,rmax); + + return kTRUE; +} + +void MHCalibrationChargeBlindPix::ExitFit() +{ + + + // + // 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; + } + + fProb = fSinglePheFit->GetProb(); + fChisquare = fSinglePheFit->GetChisquare(); + fNDF = fSinglePheFit->GetNDF(); + + *fLog << all << "Results of the Blind Pixel Fit: " << endl; + *fLog << all << "Chisquare: " << fChisquare << endl; + *fLog << all << "DoF: " << fNDF << endl; + *fLog << all << "Probability: " << fProb << endl; + +} + + +Bool_t MHCalibrationChargeBlindPix::FitSinglePhe(Option_t *opt) +{ + + if (!InitFit()) + return kFALSE; + + fHGausHist.Fit(fSinglePheFit,opt); + + ExitFit(); + + // + // The fit result is accepted under condition: + // 1) The results are not nan's + // 2) The NDF is not smaller than fNDFLimit (5) + // 3) The Probability is greater than fProbLimit (default 0.001 == 99.9%) + // 4) at least fNumSinglePheLimit events are in the single Photo-electron peak + if ( TMath::IsNaN(fLambda) + || TMath::IsNaN(fLambdaErr) + || TMath::IsNaN(fProb) + || TMath::IsNaN(fMu0) + || TMath::IsNaN(fMu0Err) + || TMath::IsNaN(fMu1) + || TMath::IsNaN(fMu1Err) + || TMath::IsNaN(fSigma0) + || TMath::IsNaN(fSigma0Err) + || TMath::IsNaN(fSigma1) + || TMath::IsNaN(fSigma1Err) + || fNDF < fNDFLimit + || fProb < fProbLimit ) + return kFALSE; + + const Stat_t entries = fHGausHist.Integral("width"); + const Float_t numSinglePhe = TMath::Exp(-1.0*fLambda)*fLambda*entries; + + if (numSinglePhe < fNumSinglePheLimit) + { + *fLog << warn << "WARNING - Statistics is too low: Only " << numSinglePhe + << " in the Single Photo-Electron peak " << endl; + return kFALSE; + } + else + *fLog << all << numSinglePhe << " in Single Photo-Electron peak " << endl; + + SetSinglePheFitOK(); + return kTRUE; +} + +void MHCalibrationChargeBlindPix::FitPedestal (Option_t *opt) +{ + + // Perform the cross-check fitting only the pedestal: + const Axis_t rmin = 0.; + const Axis_t rmax = fHGausHist.GetBinCenter(fHGausHist.GetMaximumBin()); + + FitGaus(opt, rmin, rmax); + + const Stat_t entries = fHGausHist.Integral("width"); + const Double_t fitarea = fFGausFit->GetParameter(0); + const Double_t pedarea = fitarea * TMath::Sqrt(TMath::TwoPi()) * fFGausFit->GetParameter(2); + + fLambdaCheck = TMath::Log(entries/pedarea); + fLambdaCheckErr = fFGausFit->GetParError(0)/fFGausFit->GetParameter(0) + + fFGausFit->GetParError(2)/fFGausFit->GetParameter(2); + + + SetPedestalFitOK(); + return; +} + + +// ------------------------------------------------------------------------- +// +// Draw a legend with the fit results +// +void MHCalibrationChargeBlindPix::DrawLegend() +{ + + if (!fFitLegend) + { + fFitLegend = new TPaveText(0.05,0.05,0.95,0.95); + fFitLegend->SetLabel(Form("%s%s", "Results of the single PhE Fit (", + (fFitFunc = kEPoisson4) ? "Poisson(k=4))" : + (fFitFunc = kEPoisson5) ? "Poisson(k=5))" : + (fFitFunc = kEPoisson6) ? "Poisson(k=4))" : + (fFitFunc = kEPolya ) ? "Polya(k=4))" : + (fFitFunc = kEMichele ) ? "Michele)" : " none )" )); + fFitLegend->SetTextSize(0.05); + fFitLegend->SetBit(kCanDelete); + } + else + fFitLegend->Clear(); + + const TString line1 = + Form("Mean: #lambda = %2.2f #pm %2.2f",fLambda,fLambdaErr); + TText *t1 = fFitLegend->AddText(line1.Data()); + t1->SetBit(kCanDelete); + + const TString line6 = + Form("Mean #lambda (check) = %2.2f #pm %2.2f",fLambdaCheck,fLambdaCheckErr); + TText *t2 = fFitLegend->AddText(line6.Data()); + t2->SetBit(kCanDelete); + + const TString line2 = + Form("Pedestal: #mu_{0} = %2.2f #pm %2.2f",fMu0,fMu0Err); + TText *t3 = fFitLegend->AddText(line2.Data()); + t3->SetBit(kCanDelete); + + const TString line3 = + Form("Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",fSigma0,fSigma0Err); + TText *t4 = fFitLegend->AddText(line3.Data()); + t4->SetBit(kCanDelete); + + const TString line4 = + Form("1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",fMu1,fMu1Err); + 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",fSigma1,fSigma1Err); + TText *t6 = fFitLegend->AddText(line5.Data()); + t6->SetBit(kCanDelete); + + const TString line7 = + Form("#chi^{2} / N_{dof}: %4.2f / %3i",fChisquare,fNDF); + TText *t7 = fFitLegend->AddText(line7.Data()); + t7->SetBit(kCanDelete); + + const TString line8 = + Form("Probability: %4.2f ",fProb); + TText *t8 = fFitLegend->AddText(line8.Data()); + t8->SetBit(kCanDelete); + + if (IsSinglePheFitOK()) + { + 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->SetFillColor(IsSinglePheFitOK() ? 80 : 2); + fFitLegend->Draw(); + + return; +} + + +// ------------------------------------------------------------------------- +// +// Draw the histogram +// +void MHCalibrationChargeBlindPix::Draw(Option_t *opt) +{ + + TString option(opt); + option.ToLower(); + + Int_t win = 1; + + TVirtualPad *oldpad = gPad ? gPad : MH::MakeDefCanvas(this,900, 600); + TVirtualPad *pad = NULL; + + oldpad->SetBorderMode(0); + + if (option.Contains("all")) + { + option.ReplaceAll("all",""); + oldpad->Divide(2,1); + win = 2; + oldpad->cd(1); + TVirtualPad *newpad = gPad; + pad = newpad; + pad->Divide(2,2); + pad->cd(1); + } + else + { + pad = oldpad; + pad->Divide(2,2); + pad->cd(1); + } + + if (!IsEmpty()) + gPad->SetLogy(); + + gPad->SetTicks(); + + fHGausHist.Draw(opt); + if (fFGausFit) + { + fFGausFit->SetLineColor(kBlue); + fFGausFit->Draw("same"); + } + if (fSinglePheFit) + { + fSinglePheFit->SetLineColor(IsSinglePheFitOK() ? kGreen : kRed); + fSinglePheFit->Draw("same"); + } + + pad->cd(2); + DrawLegend(); + + pad->cd(3); + if (fHSinglePheFADCSlices) + delete fHSinglePheFADCSlices; + + fHSinglePheFADCSlices = new TH1F(fASinglePheFADCSlices); + fHSinglePheFADCSlices->SetName("SinglePheFADCSlices"); + fHSinglePheFADCSlices->SetTitle(Form("%s%f","Assumed Single Phe FADC Slices, Sum > ",fSinglePheCut)); + fHSinglePheFADCSlices->SetXTitle("FADC slice number"); + fHSinglePheFADCSlices->SetYTitle("FADC counts"); + fHSinglePheFADCSlices->Draw(opt); + + pad->cd(4); + if (fHPedestalFADCSlices) + delete fHPedestalFADCSlices; + + fHPedestalFADCSlices = new TH1F(fAPedestalFADCSlices); + fHPedestalFADCSlices->SetName("PedestalFADCSlices"); + fHPedestalFADCSlices->SetTitle(Form("%s%f","Assumed Single Phe FADC Slices, Sum > ",fSinglePheCut)); + fHPedestalFADCSlices->SetXTitle("FADC slice number"); + fHPedestalFADCSlices->SetYTitle("FADC counts"); + fHPedestalFADCSlices->Draw(opt); + + if (win < 2) + return; + + oldpad->cd(2); + MHGausEvents::Draw("fourierevents"); +} + + + + + + + + + + + Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.h (revision 3289) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.h (revision 3289) @@ -0,0 +1,520 @@ +#ifndef MARS_MHCalibrationChargeBlindPix +#define MARS_MHCalibrationChargeBlindPix + + +#ifndef MARS_MHCalibrationChargePix +#include "MHCalibrationChargePix.h" +#endif + +#ifndef ROOT_TMatrix +#include +#endif + +#ifndef ROOT_TF1 +#include +#endif + +class TH1F; +class TF1; +class TPaveText; +class TText; +class MRawEvtData; +class MRawEvtPixelIter; +class MCalibrationChargeBlindPix; +class MExtractedSignalBlindPixel; +class MHCalibrationChargeBlindPix : public MHCalibrationChargePix +{ +private: + + MCalibrationChargeBlindPix *fBlindPix; //! Storage container of the results + MExtractedSignalBlindPixel *fSignal; //! Storage container of the extracted signal + MRawEvtData *fRawEvt; //! Storage container of the raw data + + static const Int_t fgChargeNbins; + static const Axis_t fgChargeFirst; + static const Axis_t fgChargeLast; + + static const Float_t fgSinglePheCut; + + static const Float_t fgNumSinglePheLimit; + + static const Double_t gkElectronicAmp; + static const Double_t gkElectronicAmpErr; + + TVector fASinglePheFADCSlices; // Array containing the averaged FADC slice entries for the supposed single-phe events + TVector fAPedestalFADCSlices; // Array containing the averaged FADC slice entries for the supposed pedestal events + + TF1 *fSinglePheFit; //-> Single Phe Fit (Gaussians convoluted with Poisson) (needs to be a pointer, initialized later) + + UInt_t fNumSinglePhes; + UInt_t fNumPedestals; + + Float_t fSinglePheCut; + + Float_t fNumSinglePheLimit; + + Double_t fLambda; + Double_t fMu0; + Double_t fMu1; + Double_t fSigma0; + Double_t fSigma1; + + Double_t fLambdaErr; + Double_t fMu0Err; + Double_t fMu1Err; + Double_t fSigma0Err; + Double_t fSigma1Err; + + Double_t fLambdaCheck; + Double_t fLambdaCheckErr; + + Double_t fChisquare; + Int_t fNDF; + Double_t fProb; + + Double_t fMeanPedestal; + Double_t fMeanPedestalErr; + Double_t fSigmaPedestal; + Double_t fSigmaPedestalErr; + + Float_t fExtractSlices; + + Byte_t fFlags; + + enum { kSinglePheFitOK, kPedestalFitOK }; + + TPaveText *fFitLegend; //! + TH1F *fHSinglePheFADCSlices; //! + TH1F *fHPedestalFADCSlices; //! + + + // Fill histos + void FillSinglePheFADCSlices(const MRawEvtPixelIter &iter); + void FillPedestalFADCSlices( const MRawEvtPixelIter &iter); + + // Fit + Bool_t InitFit(); + void ExitFit(); + +public: + + MHCalibrationChargeBlindPix(const char *name=NULL, const char *title=NULL); + ~MHCalibrationChargeBlindPix(); + + void Clear(Option_t *o=""); + void Init(); + + Bool_t SetupFill(const MParList *pList); + Bool_t ReInit ( MParList *pList); + Bool_t Fill (const MParContainer *par, const Stat_t w=1); + Bool_t Finalize(); + + // Setters + void SetChargeNbins ( const Int_t bins =fgChargeNbins ) { fChargeNbins = bins; } + void SetChargeFirst ( const Axis_t first=fgChargeFirst ) { fChargeFirst = first; } + void SetChargeLast ( const Axis_t last =fgChargeLast ) { fChargeLast = last; } + void SetSinglePheCut ( const Float_t cut =fgSinglePheCut ) { fSinglePheCut = cut; } + void SetNumSinglePheLimit ( const Float_t lim =fgNumSinglePheLimit ) { fNumSinglePheLimit = lim; } + + void SetMeanPedestal ( const Float_t f ) { fMeanPedestal = f; } + void SetMeanPedestalErr ( const Float_t f ) { fMeanPedestalErr = f; } + void SetSigmaPedestal ( const Float_t f ) { fSigmaPedestal = f; } + void SetSigmaPedestalErr ( const Float_t f ) { fSigmaPedestalErr = f; } + + void SetSinglePheFitOK ( const Bool_t b=kTRUE); + void SetPedestalFitOK( const Bool_t b=kTRUE); + + // Getters + const Double_t GetLambda() const { return fLambda; } + const Double_t GetLambdaCheck() const { return fLambdaCheck; } + const Double_t GetMu0() const { return fMu0; } + const Double_t GetMu1() const { return fMu1; } + const Double_t GetSigma0() const { return fSigma0; } + const Double_t GetSigma1() const { return fSigma1; } + + const Double_t GetLambdaErr() const { return fLambdaErr; } + const Double_t GetLambdaCheckErr() const { return fLambdaCheckErr; } + const Double_t GetMu0Err() const { return fMu0Err; } + const Double_t GetMu1Err() const { return fMu1Err; } + const Double_t GetSigma0Err() const { return fSigma0Err; } + const Double_t GetSigma1Err() const { return fSigma1Err; } + + TVector &GetASinglePheFADCSlices() { return fASinglePheFADCSlices; } + const TVector &GetASinglePheFADCSlices() const { return fASinglePheFADCSlices; } + + TVector &GetAPedestalFADCSlices() { return fAPedestalFADCSlices; } + const TVector &GetAPedestalFADCSlices() const { return fAPedestalFADCSlices; } + + const Bool_t IsSinglePheFitOK() const; + const Bool_t IsPedestalFitOK() const; + + // Draws + void Draw(Option_t *opt=""); + +private: + void DrawLegend(); + + // Fits +public: + enum FitFunc_t { kEPoisson4, kEPoisson5, kEPoisson6, kEPoisson7, kEPolya, kEMichele }; + +private: + FitFunc_t fFitFunc; + +public: + Bool_t FitSinglePhe (Option_t *opt="RL0+Q"); + void FitPedestal (Option_t *opt="RL0+Q"); + + void ChangeFitFunc(FitFunc_t func) { fFitFunc = func; } + + // Simulation + Bool_t SimulateSinglePhe(Double_t lambda, + Double_t mu0,Double_t mu1, + Double_t sigma0,Double_t sigma1); + +private: + + const static Double_t fNoWay = 10000000000.0; + + inline static Double_t fFitFuncMichele(Double_t *x, Double_t *par) + { + + Double_t lambda1cat = par[0]; + Double_t lambda1dyn = par[1]; + Double_t mu0 = par[2]; + Double_t mu1cat = par[3]; + Double_t mu1dyn = par[4]; + Double_t sigma0 = par[5]; + Double_t sigma1cat = par[6]; + Double_t sigma1dyn = par[7]; + + Double_t sumcat = 0.; + Double_t sumdyn = 0.; + Double_t arg = 0.; + + if (mu1cat < mu0) + return fNoWay; + + if (sigma1cat < sigma0) + return fNoWay; + + // if (sigma1cat < sigma1dyn) + // return NoWay; + + //if (mu1cat < mu1dyn) + // return NoWay; + + // if (lambda1cat < lambda1dyn) + // return NoWay; + + Double_t mu2cat = (2.*mu1cat)-mu0; + Double_t mu2dyn = (2.*mu1dyn)-mu0; + Double_t mu3cat = (3.*mu1cat)-(2.*mu0); + Double_t mu3dyn = (3.*mu1dyn)-(2.*mu0); + + Double_t sigma2cat = TMath::Sqrt((2.*sigma1cat*sigma1cat) - (sigma0*sigma0)); + Double_t sigma2dyn = TMath::Sqrt((2.*sigma1dyn*sigma1dyn) - (sigma0*sigma0)); + Double_t sigma3cat = TMath::Sqrt((3.*sigma1cat*sigma1cat) - (2.*sigma0*sigma0)); + Double_t sigma3dyn = TMath::Sqrt((3.*sigma1dyn*sigma1dyn) - (2.*sigma0*sigma0)); + + Double_t lambda2cat = lambda1cat*lambda1cat; + Double_t lambda2dyn = lambda1dyn*lambda1dyn; + Double_t lambda3cat = lambda2cat*lambda1cat; + Double_t lambda3dyn = lambda2dyn*lambda1dyn; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sumcat = TMath::Exp(-0.5*arg*arg)/sigma0; + sumdyn =sumcat; + + // k=1cat: + arg = (x[0] - mu1cat)/sigma1cat; + sumcat += lambda1cat*TMath::Exp(-0.5*arg*arg)/sigma1cat; + // k=1dyn: + arg = (x[0] - mu1dyn)/sigma1dyn; + sumdyn += lambda1dyn*TMath::Exp(-0.5*arg*arg)/sigma1dyn; + + // k=2cat: + arg = (x[0] - mu2cat)/sigma2cat; + sumcat += 0.5*lambda2cat*TMath::Exp(-0.5*arg*arg)/sigma2cat; + // k=2dyn: + arg = (x[0] - mu2dyn)/sigma2dyn; + sumdyn += 0.5*lambda2dyn*TMath::Exp(-0.5*arg*arg)/sigma2dyn; + + + // k=3cat: + arg = (x[0] - mu3cat)/sigma3cat; + sumcat += 0.1666666667*lambda3cat*TMath::Exp(-0.5*arg*arg)/sigma3cat; + // k=3dyn: + arg = (x[0] - mu3dyn)/sigma3dyn; + sumdyn += 0.1666666667*lambda3dyn*TMath::Exp(-0.5*arg*arg)/sigma3dyn; + + sumcat = TMath::Exp(-1.*lambda1cat)*sumcat; + sumdyn = TMath::Exp(-1.*lambda1dyn)*sumdyn; + + return par[8]*(sumcat+sumdyn)/2.; + + } + + inline static Double_t fPoissonKto4(Double_t *x, Double_t *par) + { + + Double_t lambda = par[0]; + + Double_t sum = 0.; + Double_t arg = 0.; + + Double_t mu0 = par[1]; + Double_t mu1 = par[2]; + + if (mu1 < mu0) + return fNoWay; + + Double_t sigma0 = par[3]; + Double_t sigma1 = par[4]; + + if (sigma1 < sigma0) + return fNoWay; + + Double_t mu2 = (2.*mu1)-mu0; + Double_t mu3 = (3.*mu1)-(2.*mu0); + Double_t mu4 = (4.*mu1)-(3.*mu0); + + Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); + Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); + Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); + + Double_t lambda2 = lambda*lambda; + Double_t lambda3 = lambda2*lambda; + Double_t lambda4 = lambda3*lambda; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + return TMath::Exp(-1.*lambda)*par[5]*sum; + + } + + + inline static Double_t fPoissonKto5(Double_t *x, Double_t *par) + { + + Double_t lambda = par[0]; + + Double_t sum = 0.; + Double_t arg = 0.; + + Double_t mu0 = par[1]; + Double_t mu1 = par[2]; + + if (mu1 < mu0) + return fNoWay; + + Double_t sigma0 = par[3]; + Double_t sigma1 = par[4]; + + if (sigma1 < sigma0) + return fNoWay; + + + Double_t mu2 = (2.*mu1)-mu0; + Double_t mu3 = (3.*mu1)-(2.*mu0); + Double_t mu4 = (4.*mu1)-(3.*mu0); + Double_t mu5 = (5.*mu1)-(4.*mu0); + + Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); + Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); + Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); + Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0)); + + Double_t lambda2 = lambda*lambda; + Double_t lambda3 = lambda2*lambda; + Double_t lambda4 = lambda3*lambda; + Double_t lambda5 = lambda4*lambda; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + // k=5: + arg = (x[0] - mu5)/sigma5; + sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5; + + return TMath::Exp(-1.*lambda)*par[5]*sum; + + } + + + inline static Double_t fPoissonKto6(Double_t *x, Double_t *par) + { + + Double_t lambda = par[0]; + + Double_t sum = 0.; + Double_t arg = 0.; + + Double_t mu0 = par[1]; + Double_t mu1 = par[2]; + + if (mu1 < mu0) + return fNoWay; + + Double_t sigma0 = par[3]; + Double_t sigma1 = par[4]; + + if (sigma1 < sigma0) + return fNoWay; + + + Double_t mu2 = (2.*mu1)-mu0; + Double_t mu3 = (3.*mu1)-(2.*mu0); + Double_t mu4 = (4.*mu1)-(3.*mu0); + Double_t mu5 = (5.*mu1)-(4.*mu0); + Double_t mu6 = (6.*mu1)-(5.*mu0); + + Double_t sigma2 = TMath::Sqrt((2.*sigma1*sigma1) - (sigma0*sigma0)); + Double_t sigma3 = TMath::Sqrt((3.*sigma1*sigma1) - (2.*sigma0*sigma0)); + Double_t sigma4 = TMath::Sqrt((4.*sigma1*sigma1) - (3.*sigma0*sigma0)); + Double_t sigma5 = TMath::Sqrt((5.*sigma1*sigma1) - (4.*sigma0*sigma0)); + Double_t sigma6 = TMath::Sqrt((6.*sigma1*sigma1) - (5.*sigma0*sigma0)); + + Double_t lambda2 = lambda*lambda; + Double_t lambda3 = lambda2*lambda; + Double_t lambda4 = lambda3*lambda; + Double_t lambda5 = lambda4*lambda; + Double_t lambda6 = lambda5*lambda; + + // k=0: + arg = (x[0] - mu0)/sigma0; + sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + // k=5: + arg = (x[0] - mu5)/sigma5; + sum += 0.008333333333333*lambda5*TMath::Exp(-0.5*arg*arg)/sigma5; + + // k=6: + arg = (x[0] - mu6)/sigma6; + sum += 0.001388888888889*lambda6*TMath::Exp(-0.5*arg*arg)/sigma6; + + return TMath::Exp(-1.*lambda)*par[5]*sum; + + } + + inline static Double_t fPolya(Double_t *x, Double_t *par) + { + + const Double_t QEcat = 0.247; // mean quantum efficiency + const Double_t sqrt2 = 1.4142135623731; + const Double_t sqrt3 = 1.7320508075689; + const Double_t sqrt4 = 2.; + + const Double_t lambda = par[0]; // mean number of photons + + const Double_t excessPoisson = par[1]; // non-Poissonic noise contribution + const Double_t delta1 = par[2]; // amplification first dynode + const Double_t delta2 = par[3]; // amplification subsequent dynodes + + const Double_t electronicAmpl = par[4]; // electronic amplification and conversion to FADC charges + + const Double_t pmtAmpl = delta1*delta2*delta2*delta2*delta2*delta2; // total PMT gain + const Double_t A = 1. + excessPoisson - QEcat + + 1./delta1 + + 1./delta1/delta2 + + 1./delta1/delta2/delta2; // variance contributions from PMT and QE + + const Double_t totAmpl = QEcat*pmtAmpl*electronicAmpl; // Total gain and conversion + + const Double_t mu0 = par[7]; // pedestal + const Double_t mu1 = totAmpl; // single phe position + const Double_t mu2 = 2*totAmpl; // double phe position + const Double_t mu3 = 3*totAmpl; // triple phe position + const Double_t mu4 = 4*totAmpl; // quadruple phe position + + const Double_t sigma0 = par[5]; + const Double_t sigma1 = electronicAmpl*pmtAmpl*TMath::Sqrt(QEcat*A); + const Double_t sigma2 = sqrt2*sigma1; + const Double_t sigma3 = sqrt3*sigma1; + const Double_t sigma4 = sqrt4*sigma1; + + const Double_t lambda2 = lambda*lambda; + const Double_t lambda3 = lambda2*lambda; + const Double_t lambda4 = lambda3*lambda; + + //-- calculate the area---- + Double_t arg = (x[0] - mu0)/sigma0; + Double_t sum = TMath::Exp(-0.5*arg*arg)/sigma0; + + // k=1: + arg = (x[0] - mu1)/sigma1; + sum += lambda*TMath::Exp(-0.5*arg*arg)/sigma1; + + // k=2: + arg = (x[0] - mu2)/sigma2; + sum += 0.5*lambda2*TMath::Exp(-0.5*arg*arg)/sigma2; + + // k=3: + arg = (x[0] - mu3)/sigma3; + sum += 0.1666666667*lambda3*TMath::Exp(-0.5*arg*arg)/sigma3; + + // k=4: + arg = (x[0] - mu4)/sigma4; + sum += 0.041666666666667*lambda4*TMath::Exp(-0.5*arg*arg)/sigma4; + + return TMath::Exp(-1.*lambda)*par[6]*sum; + } + + + + ClassDef(MHCalibrationChargeBlindPix, 1) // Histograms from the Calibration Blind Pixel +}; + +#endif /* MARS_MHCalibrationChargeBlindPix */