Index: /trunk/MagicSoft/Mars/mcalib/CalibLinkDef.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/CalibLinkDef.h (revision 3296) +++ /trunk/MagicSoft/Mars/mcalib/CalibLinkDef.h (revision 3297) @@ -14,5 +14,4 @@ #pragma link C++ class MCalibrationChargePINDiode+; -#pragma link C++ class MHCalibrationChargeBlindPix+; #pragma link C++ class MHCalibrationChargeBlindPix+; #pragma link C++ class MHCalibrationChargePix+; Index: unk/MagicSoft/Mars/mcalib/MCalibrationBlindPix.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MCalibrationBlindPix.cc (revision 3296) +++ (revision ) @@ -1,144 +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 11/2003 -! -! Copyright: MAGIC Software Development, 2000-2001 -! -! -\* ======================================================================== */ - -///////////////////////////////////////////////////////////////////////////// -// // -// MCalibrationBlindPix // -// // -// This is the storage container to hold informations about the calibration// -// blind pixel // -// // -///////////////////////////////////////////////////////////////////////////// -#include "MCalibrationBlindPix.h" -#include "MHCalibrationBlindPixel.h" - -#include - -#include "MLog.h" -#include "MLogManip.h" - -ClassImp(MCalibrationBlindPix); - -using namespace std; -// -------------------------------------------------------------------------- -// -// Default Constructor. -// -MCalibrationBlindPix::MCalibrationBlindPix(const char *name, const char *title) -{ - - fName = name ? name : "MCalibrationBlindPix"; - fTitle = title ? title : "Container of the MHCalibrationBlindPixel and the fit results"; - - fHist = new MHCalibrationBlindPixel(); - - if (!fHist) - *fLog << warn << dbginf << " Could not create MHCalibrationBlindPixel " << endl; - - Clear(); -} - -MCalibrationBlindPix::~MCalibrationBlindPix() -{ - delete fHist; -} - -// ------------------------------------------------------------------------ -// -// Invalidate values -// -void MCalibrationBlindPix::Clear(Option_t *o) -{ - fHist->Reset(); - - fLambda = -1.; - fMu0 = -1.; - fMu1 = -1.; - fSigma0 = -1.; - fSigma1 = -1.; - fErrLambda = -1.; - fErrMu0 = -1.; - fErrMu1 = -1.; - fErrSigma0 = -1.; - fErrSigma1 = -1.; - fTime = -1.; - fErrTime = -1; - -} - -Bool_t MCalibrationBlindPix::IsFitOK() const -{ - return fHist->IsFitOK(); -} - -Bool_t MCalibrationBlindPix::FillCharge(const Int_t q) -{ - return fHist->FillBlindPixelCharge(q); -} - -Bool_t MCalibrationBlindPix::FillTime(const Float_t t) -{ - return fHist->FillBlindPixelTime(t); -} - -Bool_t MCalibrationBlindPix::FitCharge() -{ - - if (!fHist->FitSinglePhe()) - return kFALSE; - - fLambda = TMath::IsNaN(fHist->GetLambda()) ? -1. : fHist->GetLambda(); - fMu0 = TMath::IsNaN(fHist->GetMu0()) ? -1. : fHist->GetMu0(); - fMu1 = TMath::IsNaN(fHist->GetMu1()) ? -1. : fHist->GetMu1(); - fSigma0 = TMath::IsNaN(fHist->GetSigma0()) ? -1. : fHist->GetSigma0(); - fSigma1 = TMath::IsNaN(fHist->GetSigma1()) ? -1. : fHist->GetSigma1(); - - fErrLambda = TMath::IsNaN(fHist->GetLambdaErr()) ? -1. : fHist->GetLambdaErr(); - fErrMu0 = TMath::IsNaN(fHist->GetMu0Err()) ? -1. : fHist->GetMu0Err(); - fErrMu1 = TMath::IsNaN(fHist->GetMu1Err()) ? -1. : fHist->GetMu1Err(); - fErrSigma0 = TMath::IsNaN(fHist->GetSigma0Err()) ? -1. : fHist->GetSigma0Err(); - fErrSigma1 = TMath::IsNaN(fHist->GetSigma1Err()) ? -1. : fHist->GetSigma1Err(); - - return kTRUE; -} - - - -Bool_t MCalibrationBlindPix::FitTime() -{ - - if(!fHist->FitTime()) - return kFALSE; - - fTime = fHist->GetMeanTime(); - fErrTime = fHist->GetMeanTimeErr(); - - return kTRUE; - -} - -Bool_t MCalibrationBlindPix::CheckOscillations() -{ - return fHist->CheckOscillations(); -} Index: unk/MagicSoft/Mars/mcalib/MCalibrationBlindPix.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MCalibrationBlindPix.h (revision 3296) +++ (revision ) @@ -1,83 +1,0 @@ -#ifndef MARS_MCalibrationBlindPix -#define MARS_MCalibrationBlindPix - -#ifndef MARS_MHCalibrationBlindPixel -#include "MHCalibrationBlindPixel.h" -#endif - -#ifndef MARS_MH -#include "MH.h" -#endif - -class TH1I; -class MCalibrationBlindPix : public MParContainer -{ -private: - - MHCalibrationBlindPixel *fHist; // Pointer to the histograms performing the fits, etc. - - Float_t fLambda; // The mean charge after the fit - Float_t fMu0; // The position of the pedestal-peak - Float_t fMu1; // The position of the first phe-peak - Float_t fSigma0; // The width of the pedestal-peak - Float_t fSigma1; // The width of the first phe-peak - - Float_t fErrLambda; // The error of the mean charge after the fit - Float_t fErrMu0; // The error of the position of the pedestal-peak - Float_t fErrMu1; // The error of the position of the first phe-peak - Float_t fErrSigma0; // The error of the width of the pedestal-peak - Float_t fErrSigma1; // The error of the width of the first phe-peak - - Float_t fTime; // The mean arrival time after the fit - Float_t fErrTime; // The error of the mean arrival time after the fit - - -public: - - MCalibrationBlindPix(const char *name=NULL, const char *title=NULL); - ~MCalibrationBlindPix(); - - void Clear(Option_t *o=""); - - // Getters - Float_t GetLambda() const { return fLambda; } - Float_t GetMu0() const { return fMu0; } - Float_t GetMu1() const { return fMu1; } - Float_t GetSigma0() const { return fSigma0; } - Float_t GetSigma1() const { return fSigma1; } - - Float_t GetErrLambda() const { return fErrLambda; } - Float_t GetErrMu0() const { return fErrMu0; } - Float_t GetErrMu1() const { return fErrMu1; } - Float_t GetErrSigma0() const { return fErrSigma0; } - Float_t GetErrSigma1() const { return fErrSigma1; } - - Float_t GetTime() const { return fTime; } - Float_t GetErrTime() const { return fErrTime; } - - MHCalibrationBlindPixel *GetHist() const { return fHist; } - MHCalibrationBlindPixel *GetHist() { return fHist; } - - Bool_t IsFitOK() const; - - // Fill histos - Bool_t FillCharge(const Int_t q); - Bool_t FillTime(const Float_t t); - Bool_t FillGraphs(const Int_t qhi, const Int_t qlo) const { return fHist->FillGraphs(qhi,qlo); } - - // Fits - Bool_t FitCharge(); - Bool_t FitTime(); - void ChangeFitFunc(MHCalibrationBlindPixel::FitFunc_t f) const { fHist->ChangeFitFunc(f); } - - // Draws - void Draw(Option_t *opt="") { fHist->Draw(opt); } - TObject *DrawClone(Option_t *opt="") const { return fHist->DrawClone(opt); } - - Bool_t CheckOscillations(); - - ClassDef(MCalibrationBlindPix, 1) // Container for Calibration Blind Pixel -}; - -#endif - Index: unk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc (revision 3296) +++ (revision ) @@ -1,1009 +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 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; - -} - Index: unk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h (revision 3296) +++ (revision ) @@ -1,525 +1,0 @@ -#ifndef MARS_MHCalibrationBlindPixel -#define MARS_MHCalibrationBlindPixel - -#ifndef MARS_MH -#include "MH.h" -#endif - -class TArrayF; -class TArrayI; -class TH1F; -class TH1I; -class TF1; -class TPaveText; - -class TMath; -class MParList; -class MHCalibrationBlindPixel : public MH -{ -private: - - static const Int_t fgBlindPixelChargeNbins; - static const Int_t fgBlindPixelTimeNbins; - static const Axis_t fgBlindPixelTimeFirst; - static const Axis_t fgBlindPixelTimeLast; - static const Double_t fgBlindPixelElectronicAmp; - static const Double_t fgBlindPixelElectronicAmpError; - - static const Int_t fPSDNbins; - static const Int_t fPulserFrequency; - - TH1I* fHBlindPixelCharge; // Histogram with the single Phe spectrum - TH1F* fHBlindPixelTime; // Variance of summed FADC slices - TH1F* fHBlindPixelPSD; // Power spectrum density of fHBlindPixelChargevsN - - TH1I* fHSinglePheFADCSlices; - TH1I* fHPedestalFADCSlices; - - TF1 *fSinglePheFit; - TF1 *fTimeGausFit; - TF1 *fSinglePhePedFit; - - TArrayF* fPSDHiGain; //-> Power spectrum density of fHiGains - TArrayF* fPSDLoGain; //-> Power spectrum density of fLoGains - - TH1I* fHPSD; //-> - TF1* fPSDExpFit; //-> - - TArrayF *fHiGains; //-> - TArrayF *fLoGains; //-> - TArrayF *fChargeXaxis; // - TArrayF *fPSDXaxis; // - - Float_t fPSDProb; - - Int_t fTotalEntries; // Number of entries - Int_t fCurrentSize; - - Axis_t fBlindPixelChargefirst; - Axis_t fBlindPixelChargelast; - - void DrawLegend(); - void CreateChargeXaxis(Int_t n); - void CreatePSDXaxis(Int_t n); - void CutArrayBorder(TArrayF *array) const; - void CutArrayBorder(TArrayI *array) const; - - TPaveText *fFitLegend; - - 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 fChisquare; - Double_t fProb; - Int_t fNdf; - - Double_t fMeanTime; - Double_t fMeanTimeErr; - Double_t fSigmaTime; - Double_t fSigmaTimeErr; - - Double_t fLambdaCheck; - Double_t fLambdaCheckErr; - - Double_t fMeanPedestal; - Double_t fMeanPedestalErr; - Double_t fSigmaPedestal; - Double_t fSigmaPedestalErr; - - Byte_t fFlags; - - enum { kFitOK, kOscillating }; - -public: - - MHCalibrationBlindPixel(const char *name=NULL, const char *title=NULL); - ~MHCalibrationBlindPixel(); - - void Clear(Option_t *o=""); - void Reset(); - - Bool_t FillBlindPixelCharge(const Int_t q); - Bool_t FillBlindPixelTime(const Float_t t); - Bool_t FillGraphs(const Int_t qhi, const Int_t qlo); - - // Setters - 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; } - - // 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; } - - const Double_t GetChiSquare() const { return fChisquare; } - const Double_t GetProb() const { return fProb; } - const Int_t GetNdf() const { return fNdf; } - - const Double_t GetMeanTime() const { return fMeanTime; } - const Double_t GetMeanTimeErr() const { return fMeanTimeErr; } - const Double_t GetSigmaTime() const { return fSigmaTime; } - const Double_t GetSigmaTimeErr() const { return fSigmaTimeErr; } - - TH1I *GetHSinglePheFADCSlices() { return fHSinglePheFADCSlices; } - TH1I *GetHPedestalFADCSlices() { return fHPedestalFADCSlices; } - - const Bool_t IsFitOK() const; - const Bool_t IsOscillating(); - - const TH1F *GetHBlindPixelPSD() const { return fHBlindPixelPSD; } - - // Draws - TObject *DrawClone(Option_t *option="") const; - void Draw(Option_t *option=""); - - - // Fits - enum FitFunc_t { kEPoisson4, kEPoisson5, kEPoisson6, kEPoisson7, kEPolya, kEMichele }; - -private: - FitFunc_t fFitFunc; - -public: - Bool_t FitSinglePhe(Axis_t rmin=0, Axis_t rmax=0, Option_t *opt="RL0+Q"); - Bool_t FitTime(Axis_t rmin=0., Axis_t rmax=0.,Option_t *opt="R0+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); - - // Others - void CutAllEdges(); - Bool_t CheckOscillations(); - -private: - - const static Double_t fNoWay = 10000000000.0; - - Bool_t InitFit(Axis_t min, Axis_t max); - void ExitFit(TF1 *f); - - 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(MHCalibrationBlindPixel, 1) // Histograms from the Calibration Blind Pixel -}; - -#endif /* MARS_MHCalibrationBlindPixel */ Index: /trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc (revision 3296) +++ /trunk/MagicSoft/Mars/mcalib/MHCalibrationChargeBlindPix.cc (revision 3297) @@ -827,5 +827,5 @@ fHPedestalFADCSlices = new TH1F(fAPedestalFADCSlices); fHPedestalFADCSlices->SetName("PedestalFADCSlices"); - fHPedestalFADCSlices->SetTitle(Form("%s%f","Assumed Single Phe FADC Slices, Sum > ",fSinglePheCut)); + fHPedestalFADCSlices->SetTitle(Form("%s%f","Pedestal FADC Slices, Sum > ",fSinglePheCut)); fHPedestalFADCSlices->SetXTitle("FADC slice number"); fHPedestalFADCSlices->SetYTitle("FADC counts"); Index: unk/MagicSoft/Mars/mcalib/MHCalibrationConfig.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationConfig.h (revision 3296) +++ (revision ) @@ -1,19 +1,0 @@ -#ifndef MARS_MHCalibrationConfig -#define MARS_MHCalibrationConfig - -///////////////////////////////////////////////////////////////////////////// -// // -// MHCalibrationConfig // -// // -// Contains all configuration data of the Calibration // -// i // -///////////////////////////////////////////////////////////////////////////// - - -// Global rejection criteria for the acceptance of a fit: Prob=0.001 == 99.7% Probability -const Float_t gkProbLimit = 0.001; - -// Square root of 2 pi: -const Float_t gkSq2Pi = 2.506628274631; - -#endif /* MARS_MHCalibrationBlindPixelConfig */ Index: unk/MagicSoft/Mars/mcalib/MHCalibrationPINDiode.cc =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationPINDiode.cc (revision 3296) +++ (revision ) @@ -1,128 +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 11/2003 -! -! Copyright: MAGIC Software Development, 2000-2002 -! -! -\* ======================================================================== */ - -////////////////////////////////////////////////////////////////////////////// -// -// MHCalibrationPINDiode -// -// Performs all the necessary fits to extract the mean number of photons -// out of the derived light flux -// -////////////////////////////////////////////////////////////////////////////// -#include "MHCalibrationPINDiode.h" -#include "MHCalibrationConfig.h" - -#include -#include -#include - -ClassImp(MHCalibrationPINDiode); - -using namespace std; - -// -------------------------------------------------------------------------- -// -// Default Constructor. -// -MHCalibrationPINDiode::MHCalibrationPINDiode(const char *name, const char *title) - : fChargeNbins(1000), - fChargevsNbins(10000) -{ - - fName = name ? name : "MHCalibrationPINDiode"; - fTitle = title ? title : "Fill the accumulated charges and times all PINDiode events and perform fits"; - - // Create a large number of bins, later we will rebin - fChargeFirstHiGain = -100.5; - fChargeLastHiGain = 1999.5; - - fHPINDiodeCharge = new TH1F("HPINDiodeCharge","Distribution of Summed FADC Slices PINDiode", - fChargeNbins,fChargeFirstHiGain,fChargeLastHiGain); - fHPINDiodeCharge->SetXTitle("Sum FADC Slices"); - fHPINDiodeCharge->SetYTitle("Nr. of events"); - fHPINDiodeCharge->Sumw2(); - fHPINDiodeCharge->SetDirectory(NULL); - - fHPINDiodeChargevsN = new TH1I("HPINDiodeChargevsN","Sum of Hi Gain Charges vs. Event Number Pixel ", - fChargevsNbins,-0.5,(Axis_t)fChargevsNbins - 0.5); - fHPINDiodeChargevsN->SetXTitle("Event Nr."); - fHPINDiodeChargevsN->SetYTitle("Sum of Hi Gain FADC slices"); - fHPINDiodeChargevsN->SetDirectory(NULL); - - Clear(); -} - -MHCalibrationPINDiode::~MHCalibrationPINDiode() -{ - delete fHPINDiodeCharge; - delete fHPINDiodeChargevsN; -} - -void MHCalibrationPINDiode::Clear(Option_t *o) -{ - fTotalEntries = 0; - - fChargeFirstHiGain = -100.5; - fChargeLastHiGain = 1999.5; - - fChargeChisquare = -1; - fChargeProb = -1; - fChargeNdf = -1; - - MHCalibrationPixel::Clear(); -} - -void MHCalibrationPINDiode::Reset() -{ - Clear(); - - fHPINDiodeCharge->Reset(); - fHPINDiodeChargevsN->Reset(); -} - - -Bool_t MHCalibrationPINDiode::FillCharge(Float_t q) -{ - return (fHPINDiodeCharge->Fill(q) > -1); -} - -Bool_t MHCalibrationPINDiode::FillChargevsN(Float_t q, Int_t n) -{ - return (fHPINDiodeChargevsN->Fill(n,q) > -1); -} - -void MHCalibrationPINDiode::CutAllEdges() -{ - - Int_t nbins = 20; - - StripZeros(fHPINDiodeCharge,nbins); - - fChargeFirstHiGain = fHPINDiodeCharge->GetBinLowEdge(fHPINDiodeCharge->GetXaxis()->GetFirst()); - fChargeLastHiGain = fHPINDiodeCharge->GetBinLowEdge(fHPINDiodeCharge->GetXaxis()->GetLast()) - +fHPINDiodeCharge->GetBinWidth(0); - StripZeros(fHPINDiodeChargevsN,0); - -} - Index: unk/MagicSoft/Mars/mcalib/MHCalibrationPINDiode.h =================================================================== --- /trunk/MagicSoft/Mars/mcalib/MHCalibrationPINDiode.h (revision 3296) +++ (revision ) @@ -1,50 +1,0 @@ -#ifndef MARS_MHCalibrationPINDiode -#define MARS_MHCalibrationPINDiode - -#ifndef MARS_MHCalibrationPixel -#include "MHCalibrationPixel.h" -#endif - -class TH1I; -class TH1F; -class TF1; -class TPaveText; - -class MHCalibrationPINDiode : public MHCalibrationPixel -{ -private: - - TH1F* fHPINDiodeCharge; // Histogram containing the summed 32 PINDiode slices - TH1F* fHPINDiodeErrCharge; // Variance of summed FADC slices - - TH1I* fHPINDiodeChargevsN; // Summed Charge vs. Event Nr. - - const Int_t fChargeNbins; - const Int_t fChargevsNbins; - - const Axis_t fTimeFirst; - const Axis_t fTimeLast; - -public: - - MHCalibrationPINDiode(const char *name=NULL, const char *title=NULL); - ~MHCalibrationPINDiode(); - - void Clear(Option_t *o=""); - void Reset(); - - // Fill histos - Bool_t FillCharge(Float_t q); - Bool_t FillChargevsN(Float_t q, Int_t n); - - // Fits -- not yet implemented - Bool_t FitCharge(Option_t *option="RQ0") { return kTRUE; } - - // Others - void CutAllEdges(); - - ClassDef(MHCalibrationPINDiode, 0) // Histograms from the Calibration PIN Diode -}; - -#endif - Index: /trunk/MagicSoft/Mars/mjobs/MJCalibration.cc =================================================================== --- /trunk/MagicSoft/Mars/mjobs/MJCalibration.cc (revision 3296) +++ /trunk/MagicSoft/Mars/mjobs/MJCalibration.cc (revision 3297) @@ -49,4 +49,6 @@ #include "MCalibrationChargeCam.h" #include "MCalibrationChargePINDiode.h" +#include "MCalibrationChargeBlindPix.h" +#include "MCalibrationChargeCalc.h" #include "MReadMarsFile.h" @@ -57,5 +59,4 @@ #include "MExtractBlindPixel.h" #include "MExtractSignal2.h" -#include "MCalibrationChargeCalc.h" #include "MFCosmics.h" #include "MContinue.h" @@ -477,12 +478,13 @@ static_cast(read).AddFiles(*fRuns); - MCalibrationChargePINDiode pindiode; - pindiode.SetColor(MCalibrationChargePINDiode::kECT1); - // // As long, as we don't have digital modules, // we have to set the color of the pulser LED by hand // - fCalibrationCam.SetColor(MCalibrationChargeCam::kECT1); + MCalibrationChargePINDiode pindiode; + pindiode.SetColor(MCalibrationChargePINDiode::kECT1); + + MCalibrationChargeBlindPix blindpix; + blindpix.SetColor(MCalibrationChargeBlindPix::kECT1); // Setup Tasklist @@ -491,4 +493,5 @@ plist.AddToList(&fCalibrationCam); plist.AddToList(&pindiode); + plist.AddToList(&blindpix); MTaskList tlist; @@ -504,4 +507,5 @@ MFillH fillpin("MHCalibrationChargePINDiode", "MExtractedSignalPINDiode"); + MFillH fillblind("MHCalibrationChargeBlindPix", "MExtractedSignalBlindPixel"); // // Apply a filter against cosmics @@ -522,4 +526,5 @@ tlist.AddToList(&cont); tlist.AddToList(&fillpin); + tlist.AddToList(&fillblind); tlist.AddToList(&calcalc);