Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 3060) +++ trunk/MagicSoft/Mars/Changelog (revision 3061) @@ -28,6 +28,9 @@ * mcalib/MHCalibrationPixel.[h,cc] + * mcalib/MHCalibrationBlindPixel.[h,cc] * mcalib/MCalibrationPix.[h,cc] + * mcalib/MCalibrationBlindPix.[h,cc] * mcalib/MCalibrationPINDiode.[h,cc] + * mcalib/MHCalibrationPINDiode.[h,cc] - remove histograms MHChargevsN..., now keep TArrays directly - check for oscillations for all pixels (and you will not trust Index: trunk/MagicSoft/Mars/mcalib/MCalibrationBlindPix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationBlindPix.cc (revision 3060) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationBlindPix.cc (revision 3061) @@ -101,9 +101,4 @@ } -Bool_t MCalibrationBlindPix::FillRChargevsTime(const Float_t rq, const Int_t t) -{ - return fHist->FillBlindPixelChargevsN(rq,t); -} - Bool_t MCalibrationBlindPix::FitCharge() { @@ -141,2 +136,7 @@ } + +void MCalibrationBlindPix::CheckOscillations() +{ + fHist->CheckOscillations(); +} Index: trunk/MagicSoft/Mars/mcalib/MCalibrationBlindPix.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationBlindPix.h (revision 3060) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationBlindPix.h (revision 3061) @@ -59,5 +59,5 @@ Bool_t FillCharge(const Float_t q); Bool_t FillTime(const Float_t t); - Bool_t FillRChargevsTime(const Float_t rq, const Int_t t); + Bool_t FillGraphs(Float_t qhi,Float_t qlo) const { return fHist->FillGraphs(qhi,qlo); } // Fits @@ -69,4 +69,6 @@ void Draw(Option_t *opt="") { fHist->Draw(opt); } TObject *DrawClone(Option_t *opt="") const { return fHist->DrawClone(opt); } + + void CheckOscillations(); ClassDef(MCalibrationBlindPix, 1) // Storage Container for Calibration information of one pixel Index: trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.cc (revision 3060) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.cc (revision 3061) @@ -117,5 +117,5 @@ const Byte_t MCalibrationCalc::fgSaturationLimit = 254; const Byte_t MCalibrationCalc::fgBlindPixelFirst = 3; -const Byte_t MCalibrationCalc::fgBlindPixelLast = 12; +const Byte_t MCalibrationCalc::fgBlindPixelLast = 14; // -------------------------------------------------------------------------- @@ -507,19 +507,21 @@ { - Float_t blindpixelsum = 0.; + Float_t blindpixelsumhi = 0.; + Float_t blindpixelsumlo = 0.; // // We need a dedicated signal extractor for the blind pixel // MPedestalPix &ped = (*fPedestals)[pixid]; - if (!CalcSignalBlindPixel(pixel.GetHiGainSamples(), blindpixelsum, ped.GetPedestal())) + if (!CalcSignalBlindPixel(pixel.GetHiGainSamples(), blindpixelsumhi, ped.GetPedestal())) return kFALSE; + + CalcSignalBlindPixel(pixel.GetLoGainSamples(), blindpixelsumlo, ped.GetPedestal()); - if (!blindpixel.FillCharge(blindpixelsum)) + blindpixel.FillGraphs(blindpixelsumhi,blindpixelsumlo); + + if (!blindpixel.FillCharge(blindpixelsumhi)) *fLog << warn << - "Overflow or Underflow occurred filling Blind Pixel sum = " << blindpixelsum << endl; + "Overflow or Underflow occurred filling Blind Pixel sum = " << blindpixelsumhi << endl; - // if (!blindpixel.FillRChargevsTime(blindpixelsum,fEvents)) - // *fLog << warn << - // "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl; } /* if use blind pixel */ @@ -684,4 +686,6 @@ } + blindpixel.CheckOscillations(); + fCalibrations->SetBlindPixelMethodValid(kTRUE); blindpixel.DrawClone(); Index: trunk/MagicSoft/Mars/mcalib/MCalibrationPix.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MCalibrationPix.cc (revision 3060) +++ trunk/MagicSoft/Mars/mcalib/MCalibrationPix.cc (revision 3061) @@ -350,5 +350,8 @@ CalcFFactorMethod(); - return (fRSigmaCharge/fCharge)*TMath::Sqrt(fPheFFactorMethod); + if (fPheFFactorMethod > 0) + return (fRSigmaCharge/fCharge)*TMath::Sqrt(fPheFFactorMethod); + else + return -1.; } @@ -891,3 +894,2 @@ } - Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc (revision 3060) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc (revision 3061) @@ -63,12 +63,15 @@ #include "MHCalibrationConfig.h" + #include #include #include +#include #include #include #include +#include "MFFT.h" #include "MLog.h" #include "MLogManip.h" @@ -80,9 +83,12 @@ const Int_t MHCalibrationBlindPixel::fgBlindPixelChargeNbins = 1000; const Int_t MHCalibrationBlindPixel::fgBlindPixelTimeNbins = 22; -const Int_t MHCalibrationBlindPixel::fgBlindPixelChargevsNbins = 10000; const Axis_t MHCalibrationBlindPixel::fgBlindPixelTimeFirst = -9.00; const Axis_t MHCalibrationBlindPixel::fgBlindPixelTimeLast = 12.00; const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmp = 0.008; const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmpError = 0.002; + +const Axis_t MHCalibrationBlindPixel::fNyquistFreq = 1.0; +const Axis_t MHCalibrationBlindPixel::fMinFreq = 0.; +const Int_t MHCalibrationBlindPixel::fPSDNbins = 30; // -------------------------------------------------------------------------- @@ -95,4 +101,11 @@ fTimeGausFit(NULL), fSinglePhePedFit(NULL), + fPSDHiGain(NULL), + fPSDLoGain(NULL), + fHPSD(NULL), + fPSDExpFit(NULL), + fChargeXaxis(NULL), + fPSDXaxis(NULL), + fCurrentSize(1024), fFitLegend(NULL) { @@ -119,9 +132,6 @@ fHBlindPixelTime->SetDirectory(NULL); - fHBlindPixelChargevsN = new TH1I("HBlindPixelChargevsN","Sum of Charges vs. Event Number", - fgBlindPixelChargevsNbins,-0.5,(Axis_t)fgBlindPixelChargevsNbins-0.5); - fHBlindPixelChargevsN->SetXTitle("Event Nr."); - fHBlindPixelChargevsN->SetYTitle("Sum of FADC slices"); - fHBlindPixelChargevsN->SetDirectory(NULL); + fHiGains = new TArrayF(fCurrentSize); + fLoGains = new TArrayF(fCurrentSize); Clear(); @@ -136,5 +146,7 @@ delete fHBlindPixelCharge; delete fHBlindPixelTime; - delete fHBlindPixelChargevsN; + + delete fHiGains; + delete fLoGains; if (fHBlindPixelPSD) @@ -146,4 +158,12 @@ if(fSinglePhePedFit) delete fSinglePhePedFit; + if (fPSDExpFit) + delete fPSDExpFit; + if (fHPSD) + delete fHPSD; + if (fChargeXaxis) + delete fChargeXaxis; + if (fPSDXaxis) + delete fPSDXaxis; } @@ -151,4 +171,7 @@ void MHCalibrationBlindPixel::Clear(Option_t *o) { + + fTotalEntries = 0; + fCurrentSize = 1024; fBlindPixelChargefirst = -200.; @@ -196,4 +219,19 @@ 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; @@ -207,7 +245,148 @@ fHBlindPixelCharge->Reset(); fHBlindPixelTime->Reset(); - fHBlindPixelChargevsN->Reset(); - -} + + fHiGains->Set(1024); + fLoGains->Set(1024); + + fHiGains->Reset(0.); + fLoGains->Reset(0.); + + +} + +Bool_t MHCalibrationBlindPixel::CheckOscillations() +{ + + if (fPSDExpFit) + return IsOscillating(); + + MFFT fourier; + + fPSDLoGain = fourier.PowerSpectrumDensity(fLoGains); + fPSDHiGain = fourier.PowerSpectrumDensity(fHiGains); + + Int_t entries; + TArrayF *array; + + fHPSD = new TH1F("HBlindPixelPSD", + "Power Spectrum Density Projection Blind Pixel", + fPSDNbins,fMinFreq,fNyquistFreq); + + array = fPSDHiGain; + entries = array->GetSize(); + + for (Int_t i=0;iFill(array->At(i)); + + // + // First guesses for the fit (should be as close to reality as possible, + // + const Double_t area_guess = entries*10.; + + fPSDExpFit = new TF1("PSDExpFit","[0]*exp(-[1]*x)",0.,1.); + + fPSDExpFit->SetParameters(entries,10.); + fPSDExpFit->SetParNames("Area","slope"); + fPSDExpFit->SetParLimits(0,0.,3.*area_guess); + fPSDExpFit->SetParLimits(1,5.,20.); + fPSDExpFit->SetRange(fMinFreq,fNyquistFreq); + + 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)i,i); +} + +void MHCalibrationBlindPixel::CreateChargeXaxis(Int_t n) +{ + + if (!fChargeXaxis) + { + fChargeXaxis = new TArrayF(n); + for (Int_t i=0;iAddAt((Float_t)i,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,i); +} + +void MHCalibrationBlindPixel::CutArrayBorder(TArrayF *array) +{ + + 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(Float_t qhi,Float_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(Float_t q) @@ -221,9 +400,4 @@ } -Bool_t MHCalibrationBlindPixel::FillBlindPixelChargevsN(Stat_t rq, Int_t t) -{ - return fHBlindPixelChargevsN->Fill(t,rq) > -1; -} - // ------------------------------------------------------------------------- @@ -236,5 +410,5 @@ fFitLegend = new TPaveText(0.05,0.05,0.95,0.95); - if (fFitOK) + if (IsFitOK()) fFitLegend->SetFillColor(80); else @@ -284,5 +458,5 @@ t8->SetBit(kCanDelete); - if (fFitOK) + if (IsFitOK()) { TText *t = fFitLegend->AddText(0.,0.,"Result of the Fit: OK"); @@ -333,5 +507,5 @@ TCanvas *c = MakeDefCanvas(this,550,700); - c->Divide(2,3); + c->Divide(2,4); gROOT->SetSelectedPad(NULL); @@ -344,5 +518,5 @@ fHBlindPixelCharge->Draw(opt); - if (fFitOK) + if (IsFitOK()) fSinglePheFit->SetLineColor(kGreen); else @@ -366,20 +540,69 @@ fHBlindPixelTime->Draw(opt); - c->cd(4); - - fHBlindPixelChargevsN->Draw(opt); - + 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 event number"); + 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 HiGain charges with event number"); + gr2->SetBit(kCanDelete); + gr2->Draw("AL"); + c->Modified(); c->Update(); - - c->cd(5); - - // fHBlindPixelPSD = (TH1F*)MFFT::PowerSpectrumDensity(fHBlindPixelChargevsN); - // TH1F *hist = MFFT::PowerSpectrumDensity((*this)->fHBlindPixelChargevsN); - - // fHBlindPixelPSD->Draw(opt); + + 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->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); + + if (fPSDExpFit) + { + fPSDExpFit->SetLineColor(IsOscillating() ? kRed : kGreen); + fPSDExpFit->Draw("same"); + } + + c->Modified(); + c->Update(); + } @@ -653,5 +876,5 @@ *fLog << err << "ERROR: Fit Probability " << fProb << " is smaller than the allowed value: " << gkProbLimit << endl; - fFitOK = kFALSE; + CLRBIT(fFlags,kFitOK); return kFALSE; } @@ -663,5 +886,5 @@ *fLog << err << "ERROR: Statistics is too low: Only " << contSinglePhe << " in the Single Photo-Electron peak " << endl; - fFitOK = kFALSE; + CLRBIT(fFlags,kFitOK); return kFALSE; } @@ -669,5 +892,5 @@ *fLog << inf << contSinglePhe << " in Single Photo-Electron peak " << endl; - fFitOK = kTRUE; + SETBIT(fFlags,kFitOK); return kTRUE; @@ -678,5 +901,5 @@ { - Int_t nbins = 30; + Int_t nbins = 20; CutEdges(fHBlindPixelCharge,nbins); @@ -684,6 +907,4 @@ fBlindPixelChargefirst = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetFirst()); fBlindPixelChargelast = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetLast())+fHBlindPixelCharge->GetBinWidth(0); - - CutEdges(fHBlindPixelChargevsN,0); } Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h (revision 3060) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h (revision 3061) @@ -6,4 +6,5 @@ #endif +class TArrayF; class TH1F; class TH1I; @@ -19,13 +20,15 @@ static const Int_t fgBlindPixelChargeNbins; static const Int_t fgBlindPixelTimeNbins; - static const Int_t fgBlindPixelChargevsNbins; static const Axis_t fgBlindPixelTimeFirst; static const Axis_t fgBlindPixelTimeLast; static const Double_t fgBlindPixelElectronicAmp; static const Double_t fgBlindPixelElectronicAmpError; + + static const Axis_t fNyquistFreq; + static const Axis_t fMinFreq; + static const Int_t fPSDNbins; TH1F* fHBlindPixelCharge; // Histogram with the single Phe spectrum TH1F* fHBlindPixelTime; // Variance of summed FADC slices - TH1I* fHBlindPixelChargevsN; // Summed Charge vs. Event Nr. TH1F* fHBlindPixelPSD; // Power spectrum density of fHBlindPixelChargevsN @@ -34,11 +37,29 @@ TF1 *fSinglePhePedFit; + TArrayF* fPSDHiGain; //-> Power spectrum density of fHiGains + TArrayF* fPSDLoGain; //-> Power spectrum density of fLoGains + + TH1F* 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); TPaveText *fFitLegend; - Bool_t fFitOK; Double_t fLambda; @@ -70,4 +91,9 @@ Double_t fSigmaPedestal; Double_t fSigmaPedestalErr; + + Byte_t fFlags; + + enum { kFitOK, kOscillating }; + public: @@ -81,5 +107,5 @@ Bool_t FillBlindPixelCharge(Float_t q); Bool_t FillBlindPixelTime(Float_t t); - Bool_t FillBlindPixelChargevsN(Stat_t rq, Int_t t); + Bool_t FillGraphs(Float_t qhi, Float_t qlo); // Setters @@ -113,7 +139,7 @@ const Double_t GetSigmaTimeErr() const { return fSigmaTimeErr; } - const Bool_t IsFitOK() const { return fFitOK; } - - const TH1I *GetHBlindPixelChargevsN() const { return fHBlindPixelChargevsN; } + const Bool_t IsFitOK() const; + const Bool_t IsOscillating(); + const TH1F *GetHBlindPixelPSD() const { return fHBlindPixelPSD; } @@ -122,4 +148,5 @@ void Draw(Option_t *option=""); + // Fits enum FitFunc_t { kEPoisson4, kEPoisson5, kEPoisson6, kEPoisson7, kEPolya, kEMichele }; @@ -140,5 +167,7 @@ // Others void CutAllEdges(); - + Bool_t CheckOscillations(); + + private: Index: trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc =================================================================== --- trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc (revision 3060) +++ trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc (revision 3061) @@ -79,6 +79,6 @@ fChargeXaxis(NULL), fPSDXaxis(NULL), - fFitLegend(NULL), - fCurrentSize(1024) + fCurrentSize(1024), + fFitLegend(NULL) { @@ -271,5 +271,5 @@ fPSDNbins,fMinFreq,fNyquistFreq); - array = fPSDLoGain; + array = fPSDHiGain; }