Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc (revision 4930) +++ (revision ) @@ -1,1178 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 -// -// Histogram class for the charge calibration of the Blind Pixel. -// Stores and fits the charges and stores the averaged assumed pedestal and -// single-phe FADC slice entries. Charges are taken from MExtractedSignalBlindPix. -// Performs the Single Photo-electron fit to extract the Poisson mean and its errors. -// -// Different fits can be chosen with the function ChangeFitFunc(). -// -// The fit result is accepted under the condition that: -// 1) the Probability is greater than fProbLimit (default 0.001 == 99.7%) -// 2) at least fNumSinglePheLimit events are found in the single Photo-electron peak -// -// The single FADC slice entries are averaged and stored in fASinglePheFADCSlices, if -// their sum exceeds fSinglePheCut, otherwise in fAPedestalFADCSlices. -// -// 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 "MExtractBlindPixel.h" - -#include -#include -#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 Float_t MHCalibrationChargeBlindPix::gkSignalInitializer = -9999.; - -const Int_t MHCalibrationChargeBlindPix::fgChargeNbins = 128; -const Axis_t MHCalibrationChargeBlindPix::fgChargeFirst = -0.5; -const Axis_t MHCalibrationChargeBlindPix::fgChargeLast = 511.5; -const Float_t MHCalibrationChargeBlindPix::fgSinglePheCut = 20.; -const Float_t MHCalibrationChargeBlindPix::fgNumSinglePheLimit = 50.; -// -------------------------------------------------------------------------- -// -// Default Constructor. -// -// Sets: -// - the default number for fNbins (fgChargeNbins) -// - the default number for fFirst (fgChargeFirst) -// - the default number for fLast (fgChargeLast) -// - the default number for fSinglePheCut (fgSingePheCut) -// - the default number for fNumSinglePheLimit (fgNumSinglePheLimit) -// - the default number of bins after stripping (30) -// -// - the default name of the fHGausHist ("HCalibrationChargeBlindPix") -// - the default title of the fHGausHist ("Distribution of Summed FADC slices Blind Pixel ") -// - the default x-axis title for fHGausHist ("Sum FADC Slices") -// - the default y-axis title for fHGausHist ("Nr. of events") -// -// Initializes: -// - all pointers to NULL -// - fASinglePheFADCSlices(0); -// - fAPedestalFADCSlices(0); -// - fPixId to 0 -// -// Calls: -// - Clear() -// -MHCalibrationChargeBlindPix::MHCalibrationChargeBlindPix(const char *name, const char *title) - : fBlindPix(NULL), fSignal(NULL), fRawEvt(NULL), - fSinglePheFit(NULL), - fFitLegend(NULL), - fHSinglePheFADCSlices(NULL), fHPedestalFADCSlices(NULL) -{ - - fName = name ? name : "MHCalibrationChargeBlindPix"; - fTitle = title ? title : "Statistics of the FADC sums of Blind Pixel calibration events"; - - SetNbins( fgChargeNbins ); - SetFirst( fgChargeFirst ); - SetLast ( fgChargeLast ); - - fASinglePheFADCSlices.ResizeTo(1); - fAPedestalFADCSlices.ResizeTo(1); - - SetSinglePheCut(); - SetNumSinglePheLimit(); - SetProbLimit(0.001); - SetBinsAfterStripping(0); - - fHGausHist.SetName("HCalibrationChargeBlindPix"); - fHGausHist.SetTitle("Distribution of Summed FADC slices Blind Pixel"); - fHGausHist.SetXTitle("Signal Amplitude"); - fHGausHist.SetYTitle("Nr. of events"); - - fPixId = 0; - - Clear(); -} - -// -------------------------------------------------------------------------- -// -// Default Destructor. -// -// Deletes (if Pointer is not NULL): -// -// - fSinglePheFit -// - fFitLegend -// - fHSinglePheFADCSlices -// - fHPedestalFADCSlices -// -MHCalibrationChargeBlindPix::~MHCalibrationChargeBlindPix() -{ - - if (fSinglePheFit) - delete fSinglePheFit; - - if (fFitLegend) - delete fFitLegend; - - if (fHSinglePheFADCSlices) - delete fHSinglePheFADCSlices; - - if (fHPedestalFADCSlices) - delete fHPedestalFADCSlices; - -} - -// -------------------------------------------------------------------------- -// -// Sets: -// - all variables to 0., except the fit result variables to gkSignalInitializer -// - all flags to kFALSE -// - all pointers to NULL -// - the default fit function (kEPoisson5) -// -// Deletes: -// - all pointers unequal NULL -// -// Calls: -// - MHCalibrationChargePix::Clear() -// -void MHCalibrationChargeBlindPix::Clear(Option_t *o) -{ - - fLambda = gkSignalInitializer; - fMu0 = gkSignalInitializer; - fMu1 = gkSignalInitializer; - fSigma0 = gkSignalInitializer; - fSigma1 = gkSignalInitializer; - fLambdaErr = gkSignalInitializer; - fMu0Err = gkSignalInitializer; - fMu1Err = gkSignalInitializer; - fSigma0Err = gkSignalInitializer; - fSigma1Err = gkSignalInitializer; - - fLambdaCheck = gkSignalInitializer; - fLambdaCheckErr = gkSignalInitializer; - - // fFitFunc = kEMichele; - fFitFunc = kEPoisson5; - - fNumSinglePhes = 0; - fNumPedestals = 0; - - fChisquare = 0.; - fNDF = 0 ; - fProb = 0.; - - SetSinglePheFitOK ( kFALSE ); - SetPedestalFitOK ( kFALSE ); - - if (fFitLegend) - { - delete fFitLegend; - fFitLegend = NULL; - } - - if (fSinglePheFit) - { - delete fSinglePheFit; - fSinglePheFit = NULL; - } - - if (fHSinglePheFADCSlices) - { - delete fHSinglePheFADCSlices; - fHSinglePheFADCSlices = NULL; - } - - if (fHPedestalFADCSlices) - { - delete fHPedestalFADCSlices; - fHPedestalFADCSlices = NULL; - } - - MHCalibrationPix::Clear(); - return; -} - -// -------------------------------------------------------------------------- -// -// Empty function to overload MHCalibrationChargePix::Reset() -// -void MHCalibrationChargeBlindPix::Reset() -{ -} - -/* -// -------------------------------------------------------------------------- -// -// Our own clone function is necessary since root 3.01/06 or Mars 0.4 -// I don't know the reason. -// -// Creates new MHCalibrationCam -// -TObject *MHCalibrationChargeBlindPix::Clone(const char *) const -{ - - MHCalibrationChargeBlindPix *pix = new MHCalibrationChargeBlindPix(); - this->Copy(*pix); - - this->fHGausHist.Copy(pix->fHGausHist); - this->fSinglePheFit->Copy(*(pix->fSinglePheFit)); - this->fHSinglePheFADCSlices->Copy(*(pix->fHSinglePheFADCSlices)); - this->fHPedestalFADCSlices->Copy(*(pix->fHPedestalFADCSlices)); - - - return pix; -} -*/ - -// -------------------------------------------------------------------------- -// -// Set bit kSinglePheFitOK from outside -// -void MHCalibrationChargeBlindPix::SetSinglePheFitOK (const Bool_t b ) -{ - b ? SETBIT(fFlags,kSinglePheFitOK) : CLRBIT(fFlags,kSinglePheFitOK); -} - -// -------------------------------------------------------------------------- -// -// Set bit kPedestalFitOK from outside -// -void MHCalibrationChargeBlindPix::SetPedestalFitOK(const Bool_t b) -{ - b ? SETBIT(fFlags,kPedestalFitOK) : CLRBIT(fFlags,kPedestalFitOK); -} - -// -------------------------------------------------------------------------- -// -// Ask for status of bit kSinglePheFitOK -// -const Bool_t MHCalibrationChargeBlindPix::IsSinglePheFitOK() const -{ - return TESTBIT(fFlags,kSinglePheFitOK); -} - -// -------------------------------------------------------------------------- -// -// Ask for status of bit kPedestalFitOK -// -const Bool_t MHCalibrationChargeBlindPix::IsPedestalFitOK() const -{ - return TESTBIT(fFlags,kPedestalFitOK); -} - -// -------------------------------------------------------------------------- -// -// Gets the pointers to: -// - MRawEvtData -// - MExtractedSignalBlindPixel -// -Bool_t MHCalibrationChargeBlindPix::SetupFill(const MParList *pList) -{ - - fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData"); - if (!fRawEvt) - { - *fLog << err << "MRawEvtData not found... aborting." << endl; - return kFALSE; - } - - fSignal = (MExtractedSignalBlindPixel*)pList->FindObject("MExtractedSignalBlindPixel"); - if (!fSignal) - { - *fLog << err << "MExtractedSignalBlindPixel not found... aborting " << endl; - return kFALSE; - } - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Gets or creates the pointers to: -// - MCalibrationChargeBlindPix -// -// Calls: -// - MHGausHist::InitBins() -// -Bool_t MHCalibrationChargeBlindPix::ReInit(MParList *pList) -{ - - fBlindPix = (MCalibrationChargeBlindPix*)pList->FindCreateObj("MCalibrationChargeBlindPix"); - if (!fBlindPix) - return kFALSE; - - const Int_t samples = fSignal->GetNumFADCSamples(); - const Int_t integ = fSignal->IsExtractionType( MExtractBlindPixel::kIntegral ); - - // - // Modify the histogram size in case, integrals have been used - // - if ( fLast < samples*integ*fgChargeLast ) - { - SetLast ( samples * (fgChargeLast+0.5) - 0.5 ); - SetSinglePheCut( samples * fgSinglePheCut ); - } - - MHCalibrationPix::InitBins(); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Retrieves from MExtractedSignalBlindPixel: -// - number of FADC samples -// - extracted signal -// - blind Pixel ID -// -// Resizes (if necessary): -// - fASinglePheFADCSlices to sum of HiGain and LoGain samples -// - fAPedestalFADCSlices to sum of HiGain and LoGain samples -// -// Fills the following histograms: -// - MHGausEvents::FillHistAndArray(signal) -// -// Creates MRawEvtPixelIter, jumps to blind pixel ID, -// fills the vectors fASinglePheFADCSlices and fAPedestalFADCSlices -// with the full FADC slices, depending on the size of the signal w.r.t. fSinglePheCut -// -Bool_t MHCalibrationChargeBlindPix::Fill(const MParContainer *par, const Stat_t w) -{ - - const Int_t samples = (Int_t)fRawEvt->GetNumHiGainSamples() - +(Int_t)fRawEvt->GetNumLoGainSamples(); - - if (!fASinglePheFADCSlices.IsValid()) - { - fASinglePheFADCSlices.ResizeTo(samples); - fAPedestalFADCSlices.ResizeTo(samples); - } - - if (fASinglePheFADCSlices.GetNrows() != samples) - { - fASinglePheFADCSlices.ResizeTo(samples); - fAPedestalFADCSlices.ResizeTo(samples); - } - - 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; - } - - // - // Signal extraction and histogram filling - // - const Float_t signal = fSignal->GetExtractedSignal(fPixId); - - if (signal > -0.5) - FillHist(signal); - else - return kTRUE; - - // - // In order to study the single-phe posistion, we extract the slices - // - const Int_t blindpixIdx = fSignal->GetBlindPixelIdx(fPixId); - - MRawEvtPixelIter pixel(fRawEvt); - pixel.Jump(blindpixIdx); - - if (signal > fSinglePheCut) - FillSinglePheFADCSlices(pixel); - else - FillPedestalFADCSlices(pixel); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// Returns kFALSE, if empty -// -// - Creates the fourier spectrum and sets bit MHGausEvents::IsFourierSpectrumOK() -// - Retrieves the pedestals from MExtractedSignalBlindPixel -// - Normalizes fASinglePheFADCSlices and fAPedestalFADCSlices -// - Executes FitPedestal() -// - Executes FitSinglePhe() -// - Retrieves fit results and stores them in MCalibrationChargeBlindPix -// -Bool_t MHCalibrationChargeBlindPix::Finalize() -{ - - if (IsEmpty()) - { - *fLog << err << GetDescriptor() << " ID: " << fPixId - << " My histogram has not been filled !! " << endl; - return kFALSE; - } - - 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;iSetValid(kTRUE); - - if (FitSinglePhe()) - fBlindPix->SetSinglePheFitOK(); - else - fBlindPix->SetValid(IsPedestalFitOK()); - - fBlindPix->SetLambda ( fLambda ); - fBlindPix->SetLambdaVar ( fLambdaErr*fLambdaErr ); - 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; -} - - -// -------------------------------------------------------------------------- -// -// Checks again for the size and fills fASinglePheFADCSlices with the FADC slice entries -// -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++; -} - -// -------------------------------------------------------------------------- -// -// Checks again for the size and fills fAPedestalFADCSlices with the FADC slice entries -// -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++; -} - - -// -------------------------------------------------------------------------- -// -// Task to simulate single phe spectrum with the given parameters -// -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; -} - -// -------------------------------------------------------------------------- -// -// - Get the ranges from the stripped histogram -// - choose reasonable start values for the fit -// - initialize the fit function depending on fFitFunc -// - initialize parameter names and limits depending on fFitFunc -// -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.5; - //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); - rmin += 6.5; - break; - case kEPoisson5: - fSinglePheFit = new TF1("SinglePheFit",&fPoissonKto5,rmin,rmax,6); - rmin = 0.; - 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: - fSinglePheFit = new TF1("SinglePheFit",&fFitFuncMichele,rmin,rmax,9); - 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 = 13.5; - const Double_t si_0_guess = 2.5; - const Double_t mu_1_guess = 30.; - const Double_t si_1_guess = si_0_guess + si_0_guess; - // Michele - const Double_t lambda_1cat_guess = 1.00; - const Double_t lambda_1dyn_guess = lambda_1cat_guess/10.; - const Double_t mu_1cat_guess = 50.; - const Double_t mu_1dyn_guess = 17.; - const Double_t si_1cat_guess = si_0_guess + si_0_guess; - const Double_t si_1dyn_guess = si_0_guess + si_0_guess/2.; - // 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,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,norm); - fSinglePheFit->SetParLimits(0,0.,2.); - fSinglePheFit->SetParLimits(1,10.,17.); - fSinglePheFit->SetParLimits(2,17.,50.); - fSinglePheFit->SetParLimits(3,1.,5.); - fSinglePheFit->SetParLimits(4,5.,30.); - fSinglePheFit->SetParLimits(5,norm-(0.5*norm),norm+(0.7*norm)); - break; - case kEPoisson5: - case kEPoisson6: - fSinglePheFit->SetParNames("#lambda","#mu_{0}","#mu_{1}","#sigma_{0}","#sigma_{1}","Area"); - fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,800.,si_0_guess,500.,norm); - fSinglePheFit->SetParLimits(0,0.,2.); - fSinglePheFit->SetParLimits(1,0.,100.); - fSinglePheFit->SetParLimits(2,300.,1500.); - fSinglePheFit->SetParLimits(3,30.,250.); - fSinglePheFit->SetParLimits(4,100.,1000.); - fSinglePheFit->SetParLimits(5,norm/1.5,norm*1.5); - 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: - fSinglePheFit->SetParNames("#lambda_{cat}","#lambda_{dyn}", - "#mu_{0}","#mu_{1cat}","#mu_{1dyn}", - "#sigma_{0}","#sigma_{1cat}","#sigma_{1dyn}", - "Area"); - fSinglePheFit->SetParameters(lambda_1cat_guess, lambda_1dyn_guess, - mu_0_guess, mu_1cat_guess,mu_1dyn_guess, - si_0_guess, si_1cat_guess,si_1dyn_guess, - norm); - fSinglePheFit->SetParLimits(0,0.01,2.0); - fSinglePheFit->SetParLimits(1,0.,0.5); - fSinglePheFit->SetParLimits(2,10.,16.); - fSinglePheFit->SetParLimits(3,25.,50.); - fSinglePheFit->SetParLimits(4,16.,18.5); - fSinglePheFit->SetParLimits(5,1.,5.); - fSinglePheFit->SetParLimits(6,10.,50.); - fSinglePheFit->SetParLimits(7,5.,10.); - fSinglePheFit->SetParLimits(8,norm/2.,norm*2.5); - break; - - default: - *fLog << warn << "WARNING: Could not find Fit Function for Blind Pixel " << endl; - return kFALSE; - break; - } - - fSinglePheFit->SetRange(rmin,rmax); - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -// - Retrieve the parameters depending on fFitFunc -// - Retrieve probability, Chisquare and NDF -// -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.; - case kEMichele: - fLambda = fSinglePheFit->GetParameter(0); - fMu0 = fSinglePheFit->GetParameter(2); - fMu1 = fSinglePheFit->GetParameter(3); - fSigma0 = fSinglePheFit->GetParameter(5); - fSigma1 = fSinglePheFit->GetParameter(6); - - fLambdaErr = fSinglePheFit->GetParError(0); - fMu0Err = fSinglePheFit->GetParError(2); - fMu1Err = fSinglePheFit->GetParError(3); - fSigma0Err = fSinglePheFit->GetParError(5); - fSigma1Err = fSinglePheFit->GetParError(6); - break; - 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; - -} - -// -------------------------------------------------------------------------- -// -// - Executes InitFit() -// - Fits the fHGausHist with fSinglePheFit -// - Executes 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 -// -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; -} - -// -------------------------------------------------------------------------- -// -// - Retrieves limits for the fit -// - Fits the fHGausHist with Gauss -// - Retrieves the results to fLambdaCheck and fLambdaCheckErr -// - Sets a flag IsPedestalFitOK() -// -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()); - const Axis_t rmax = fSinglePheCut; - - FitGaus(opt, rmin, rmax); - - const Stat_t entries = fHGausHist.Integral("width"); - const Double_t pedarea = fFGausFit->Integral(0.,fSinglePheCut); - - fLambdaCheck = TMath::Log(entries/pedarea); - // estimate the error by the error of the obtained area from the Gauss-function: - fLambdaCheckErr = fFGausFit->GetParError(0)/fFGausFit->GetParameter(0); - - SetPedestalFitOK(IsGausFitOK()); - return; -} - - -// ------------------------------------------------------------------------- -// -// Draw a legend with the fit results -// -void MHCalibrationChargeBlindPix::DrawLegend(Option_t *opt) -{ - - TString option(opt); - - 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=6))" : - (fFitFunc == kEPolya ) ? "Polya(k=4))" : - (fFitFunc == kEMichele ) ? "Michele)" : " none )" )); - fFitLegend->SetTextSize(0.05); - } - 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); - - if (option.Contains("datacheck")) - { - if (fLambda + 3.*fLambdaErr < fLambdaCheck - 3.*fLambdaCheckErr - || - fLambda - 3.*fLambdaErr > fLambdaCheck + 3.*fLambdaCheckErr ) - { - TText *t = fFitLegend->AddText("#lambda and #lambda_{check} more than 3#sigma apart!"); - t->SetBit(kCanDelete); - } - } - else - { - - 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: %6.4f ",fProb); - TText *t8 = fFitLegend->AddText(line8.Data()); - t8->SetBit(kCanDelete); - - if (IsSinglePheFitOK()) - { - TText *t = fFitLegend->AddText("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 -// -// The following options can be chosen: -// -// "": displays the fHGausHist, the fits, the legend and fASinglePheFADCSlices and fAPedestalFADCSlices -// "all": executes additionally MHGausEvents::Draw(), with option "fourierevents" -// "datacheck" display the fHGausHist, the fits and the legend -// -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; - - 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 if (option.Contains("datacheck")) - { - pad = oldpad; - pad->Divide(1,2); - pad->cd(1); - fHGausHist.SetStats(0); - } - else - { - pad = oldpad; - pad->Divide(2,2); - pad->cd(1); - } - - if (!IsEmpty() && !IsOnlyOverflow() && !IsOnlyUnderflow()) - gPad->SetLogy(); - - gPad->SetTicks(); - - fHGausHist.Draw(); - if (fFGausFit ) - { - fFGausFit->SetLineColor(kBlue); - fFGausFit->Draw("same"); - if (!option.Contains("datacheck")) - { - TLine *line = new TLine(fSinglePheCut, 0., fSinglePheCut, 10.); - line->SetBit(kCanDelete); - line->SetLineColor(kBlue); - line->SetLineWidth(3); - line->DrawLine(fSinglePheCut, 0., fSinglePheCut, 2.); - } - } - - if (fSinglePheFit) - { - fSinglePheFit->SetFillStyle(0); - fSinglePheFit->SetLineWidth(3); - fSinglePheFit->SetLineColor(IsSinglePheFitOK() ? kGreen : kRed); - fSinglePheFit->Draw("same"); - } - - pad->cd(2); - DrawLegend(option.Data()); - - if (option.Contains("datacheck")) - return; - - pad->cd(3); - - if (fASinglePheFADCSlices.GetNrows()!=1) - { - if (fHSinglePheFADCSlices) - delete fHSinglePheFADCSlices; - - fHSinglePheFADCSlices = new TH1F(fASinglePheFADCSlices); - fHSinglePheFADCSlices->SetName("SinglePheFADCSlices"); - fHSinglePheFADCSlices->SetTitle(Form("%s%4.1f","Assumed Single Phe FADC Slices, Sum > ",fSinglePheCut)); - fHSinglePheFADCSlices->SetXTitle("FADC slice number"); - fHSinglePheFADCSlices->SetYTitle("FADC counts"); - const Int_t nbins = fHSinglePheFADCSlices->GetNbinsX(); - TH2D *nulls = new TH2D("Nulls",fHSinglePheFADCSlices->GetTitle(),nbins,0., - fHSinglePheFADCSlices->GetXaxis()->GetBinCenter(nbins), - 100,0.,50.); - nulls->SetDirectory(NULL); - nulls->SetBit(kCanDelete); - nulls->GetXaxis()->SetTitle(fHSinglePheFADCSlices->GetXaxis()->GetTitle()); - nulls->GetYaxis()->SetTitle(fHSinglePheFADCSlices->GetYaxis()->GetTitle()); - nulls->GetXaxis()->CenterTitle(); - nulls->GetYaxis()->CenterTitle(); - nulls->SetStats(0); - nulls->Draw(); - fHSinglePheFADCSlices->Draw("same"); - } - - pad->cd(4); - if (fAPedestalFADCSlices.GetNrows()!=1) - { - - if (fHPedestalFADCSlices) - delete fHPedestalFADCSlices; - - fHPedestalFADCSlices = new TH1F(fAPedestalFADCSlices); - fHPedestalFADCSlices->SetName("PedestalFADCSlices"); - fHPedestalFADCSlices->SetTitle(Form("%s%4.1f","Pedestal FADC Slices, Sum < ",fSinglePheCut)); - fHPedestalFADCSlices->SetXTitle("FADC slice number"); - fHPedestalFADCSlices->SetYTitle("FADC counts"); - const Int_t nbins = fHPedestalFADCSlices->GetNbinsX(); - TH2D *nullp = new TH2D("Nullp",fHPedestalFADCSlices->GetTitle(),nbins,0., - fHPedestalFADCSlices->GetXaxis()->GetBinCenter(nbins), - 100,0.,50.); - nullp->SetDirectory(NULL); - nullp->SetBit(kCanDelete); - nullp->GetXaxis()->SetTitle(fHPedestalFADCSlices->GetXaxis()->GetTitle()); - nullp->GetYaxis()->SetTitle(fHPedestalFADCSlices->GetYaxis()->GetTitle()); - nullp->GetXaxis()->CenterTitle(); - nullp->GetYaxis()->CenterTitle(); - nullp->SetStats(0); - nullp->Draw(); - fHPedestalFADCSlices->Draw("same"); - } - - if (win < 2) - return; - - oldpad->cd(2); - MHCalibrationPix::Draw("fourierevents"); -} - - - - - - - - - - - Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationTestPix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationTestPix.cc (revision 4930) +++ (revision ) @@ -1,98 +1,0 @@ -/* ======================================================================== *\ -! -! * -! * 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 -! -! -\* ======================================================================== */ -////////////////////////////////////////////////////////////////////////////// -// -// MHCalibrationTestPix -// -// Histogram class for the charge calibration. -// Stores and fits the charges and stores the location of the maximum FADC -// slice. Tests are taken from MExtractedSignalPix. -// -////////////////////////////////////////////////////////////////////////////// -#include "MHCalibrationTestPix.h" - -#include -#include - -#include -#include -#include -#include - -#include "MH.h" - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MHCalibrationTestPix); - -using namespace std; - -const Int_t MHCalibrationTestPix::fgChargeNbins = 4000; -const Axis_t MHCalibrationTestPix::fgChargeFirst = -0.5; -const Axis_t MHCalibrationTestPix::fgChargeLast = 39999.5; -// -------------------------------------------------------------------------- -// -// Default Constructor. -// -// Sets: -// - the default number for fNbins (fgChargeNbins) -// - the default number for fFirst (fgChargeFirst) -// - the default number for fLast (fgChargeLast) -// -// - the default name of the fHGausHist ("HCalibrationTest") -// - the default title of the fHGausHist ("Distribution of calibrated FADC slices Pixel ") -// - the default x-axis title for fHGausHist ("Sum FADC Slices") -// - the default y-axis title for fHGausHist ("Nr. of events") -// -// Calls: -// - Clear(); -// -MHCalibrationTestPix::MHCalibrationTestPix(const char *name, const char *title) -{ - - fName = name ? name : "MHCalibrationTestPix"; - fTitle = title ? title : "Statistics of the calibrated FADC sums of calibration events"; - - SetNbins ( fgChargeNbins ); - SetFirst ( fgChargeFirst ); - SetLast ( fgChargeLast ); - - fHGausHist.SetName("HCalibrationTest"); - fHGausHist.SetTitle("Distribution of calibrated Photons Pixel "); - fHGausHist.SetXTitle("Nr. Photons"); - fHGausHist.SetYTitle("Nr. of events"); - -} - - -// -------------------------------------------------------------------------- -// -// returns fHGausHist.Integral("width") -// -const Float_t MHCalibrationTestPix::GetIntegral() const -{ - return fHGausHist.Integral("width"); -} - Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationTestPix.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationTestPix.h (revision 4930) +++ (revision ) @@ -1,27 +1,0 @@ -#ifndef MARS_MHCalibrationTestPix -#define MARS_MHCalibrationTestPix - -#ifndef MARS_MHCalibrationPix -#include "MHCalibrationPix.h" -#endif - -class MHCalibrationTestPix : public MHCalibrationPix -{ - -private: - - static const Int_t fgChargeNbins; // Default for fNBins (now set to: 2000 ) - static const Axis_t fgChargeFirst; // Default for fFirst (now set to: -0.5 ) - static const Axis_t fgChargeLast; // Default for fLast (now set to: 1999.5) - -public: - - MHCalibrationTestPix(const char *name=NULL, const char *title=NULL); - ~MHCalibrationTestPix() {} - - const Float_t GetIntegral() const; - - ClassDef(MHCalibrationTestPix, 1) // Base Histogram class for Test Pixel Calibration -}; - -#endif