Changeset 3061 for trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc

Timestamp:

02/08/04 22:17:27 (21 years ago)

Author:

gaug

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc (modified) (20 diffs)

trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc

@@ -63 +63 @@
 #include "MHCalibrationConfig.h"
 
+
 #include <TStyle.h>
 #include <TCanvas.h>
 #include <TPaveText.h>
 
+#include <TGraph.h>
 #include <TF1.h>
 #include <TH1.h>
 #include <TRandom.h>
 
+#include "MFFT.h"
 #include "MLog.h"
 #include "MLogManip.h"
@@ -80 +83 @@
 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 +101 @@
         fTimeGausFit(NULL),
         fSinglePhePedFit(NULL),
+        fPSDHiGain(NULL),
+        fPSDLoGain(NULL),
+        fHPSD(NULL),
+        fPSDExpFit(NULL),
+        fChargeXaxis(NULL),
+        fPSDXaxis(NULL),
+        fCurrentSize(1024),
         fFitLegend(NULL)
 {
@@ -119 +132 @@
     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 +146 @@
   delete fHBlindPixelCharge;
   delete fHBlindPixelTime;
-  delete fHBlindPixelChargevsN;
+
+  delete fHiGains;
+  delete fLoGains;
 
   if (fHBlindPixelPSD)
@@ -146 +158 @@
   if(fSinglePhePedFit)
     delete fSinglePhePedFit;
+  if (fPSDExpFit)
+    delete fPSDExpFit;
+  if (fHPSD)
+    delete fHPSD;
+  if (fChargeXaxis)
+    delete fChargeXaxis;
+  if (fPSDXaxis)
+    delete fPSDXaxis;
 
 }
@@ -151 +171 @@
 void MHCalibrationBlindPixel::Clear(Option_t *o)
 {
+
+  fTotalEntries            = 0;
+  fCurrentSize             = 1024;
   
   fBlindPixelChargefirst = -200.;
@@ -196 +219 @@
   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 +245 @@
   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;i<entries;i++)
+    fHPSD->Fill(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;i<n;i++)
+    fPSDXaxis->AddAt((Float_t)i,i);
+}
+
+void MHCalibrationBlindPixel::CreateChargeXaxis(Int_t n)
+{
+  
+  if (!fChargeXaxis)
+    {
+      fChargeXaxis = new TArrayF(n);
+      for (Int_t i=0;i<n;i++)
+        fChargeXaxis->AddAt((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;i<n+diff;i++)
+      fChargeXaxis->AddAt((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 +400 @@
 }
 
-Bool_t MHCalibrationBlindPixel::FillBlindPixelChargevsN(Stat_t rq, Int_t t)
-{
-    return fHBlindPixelChargevsN->Fill(t,rq) > -1;
-}
-
 
 // -------------------------------------------------------------------------
@@ -236 +410 @@
   fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
 
-  if (fFitOK) 
+  if (IsFitOK()) 
       fFitLegend->SetFillColor(80);
   else
@@ -284 +458 @@
   t8->SetBit(kCanDelete);
 
-  if (fFitOK)
+  if (IsFitOK())
     {
       TText *t = fFitLegend->AddText(0.,0.,"Result of the Fit: OK");
@@ -333 +507 @@
     TCanvas *c = MakeDefCanvas(this,550,700);
 
-    c->Divide(2,3);
+    c->Divide(2,4);
 
     gROOT->SetSelectedPad(NULL);
@@ -344 +518 @@
     fHBlindPixelCharge->Draw(opt);
     
-    if (fFitOK)
+    if (IsFitOK())
       fSinglePheFit->SetLineColor(kGreen);          
     else
@@ -366 +540 @@
     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 +876 @@
       *fLog << err << "ERROR: Fit Probability " << fProb 
             << " is smaller than the allowed value: " << gkProbLimit << endl;
-      fFitOK = kFALSE;
+      CLRBIT(fFlags,kFitOK);
       return kFALSE;
     }
@@ -663 +886 @@
       *fLog << err << "ERROR: Statistics is too low: Only " << contSinglePhe 
             << " in the Single Photo-Electron peak " << endl;
-      fFitOK = kFALSE;
+      CLRBIT(fFlags,kFitOK);      
       return kFALSE;
     } 
@@ -669 +892 @@
     *fLog << inf << contSinglePhe << " in Single Photo-Electron peak " << endl;
   
-  fFitOK = kTRUE;
+  SETBIT(fFlags,kFitOK);
     
   return kTRUE;
@@ -678 +901 @@
 {
 
-  Int_t nbins = 30;
+  Int_t nbins = 20;
 
   CutEdges(fHBlindPixelCharge,nbins);
@@ -684 +907 @@
   fBlindPixelChargefirst = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetFirst());
   fBlindPixelChargelast  = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetLast())+fHBlindPixelCharge->GetBinWidth(0);
-
-  CutEdges(fHBlindPixelChargevsN,0);
 
 }