Changeset 3061 for trunk


Ignore:
Timestamp:
02/08/04 22:17:27 (21 years ago)
Author:
gaug
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/Changelog

    r3060 r3061  
    2828
    2929   * mcalib/MHCalibrationPixel.[h,cc]
     30   * mcalib/MHCalibrationBlindPixel.[h,cc]
    3031   * mcalib/MCalibrationPix.[h,cc]
     32   * mcalib/MCalibrationBlindPix.[h,cc]
    3133   * mcalib/MCalibrationPINDiode.[h,cc]
     34   * mcalib/MHCalibrationPINDiode.[h,cc]
    3235     - remove histograms MHChargevsN..., now keep TArrays directly
    3336     - check for oscillations for all pixels (and you will not trust

  • trunk/MagicSoft/Mars/mcalib/MCalibrationBlindPix.cc

    r3030 r3061  
    101101}
    102102
    103 Bool_t MCalibrationBlindPix::FillRChargevsTime(const Float_t rq, const Int_t t)
    104 {
    105   return fHist->FillBlindPixelChargevsN(rq,t);
    106 }
    107 
    108103Bool_t MCalibrationBlindPix::FitCharge()
    109104{
    … …  
    141136
    142137}
     138
     139void MCalibrationBlindPix::CheckOscillations()
     140{
     141  fHist->CheckOscillations();
     142}

  • trunk/MagicSoft/Mars/mcalib/MCalibrationBlindPix.h

    r3029 r3061  
    5959  Bool_t FillCharge(const Float_t q);
    6060  Bool_t FillTime(const Float_t t);
    61   Bool_t FillRChargevsTime(const Float_t rq, const Int_t t);
     61  Bool_t FillGraphs(Float_t qhi,Float_t qlo) const { return fHist->FillGraphs(qhi,qlo); }
    6262 
    6363  // Fits
    … …  
    6969  void Draw(Option_t *opt="")                  { fHist->Draw(opt); }
    7070  TObject *DrawClone(Option_t *opt="") const  {  return fHist->DrawClone(opt); } 
     71
     72  void  CheckOscillations(); 
    7173 
    7274  ClassDef(MCalibrationBlindPix, 1)     // Storage Container for Calibration information of one pixel

  • trunk/MagicSoft/Mars/mcalib/MCalibrationCalc.cc

    r3057 r3061  
    117117const Byte_t MCalibrationCalc::fgSaturationLimit = 254;
    118118const Byte_t MCalibrationCalc::fgBlindPixelFirst = 3;
    119 const Byte_t MCalibrationCalc::fgBlindPixelLast  = 12;
     119const Byte_t MCalibrationCalc::fgBlindPixelLast  = 14;
    120120
    121121// --------------------------------------------------------------------------
    … …  
    507507            {
    508508         
    509               Float_t blindpixelsum = 0.;
     509              Float_t blindpixelsumhi = 0.;
     510              Float_t blindpixelsumlo = 0.;
    510511              //
    511512              // We need a dedicated signal extractor for the blind pixel
    512513              //
    513514              MPedestalPix &ped  = (*fPedestals)[pixid];
    514               if (!CalcSignalBlindPixel(pixel.GetHiGainSamples(), blindpixelsum, ped.GetPedestal()))
     515              if (!CalcSignalBlindPixel(pixel.GetHiGainSamples(), blindpixelsumhi, ped.GetPedestal()))
    515516                return kFALSE;
     517
     518              CalcSignalBlindPixel(pixel.GetLoGainSamples(), blindpixelsumlo, ped.GetPedestal());
    516519             
    517               if (!blindpixel.FillCharge(blindpixelsum))
     520              blindpixel.FillGraphs(blindpixelsumhi,blindpixelsumlo);
     521         
     522              if (!blindpixel.FillCharge(blindpixelsumhi))
    518523                *fLog << warn <<
    519                   "Overflow or Underflow occurred filling Blind Pixel sum = " << blindpixelsum << endl;
     524                  "Overflow or Underflow occurred filling Blind Pixel sum = " << blindpixelsumhi << endl;
    520525             
    521               //              if (!blindpixel.FillRChargevsTime(blindpixelsum,fEvents))
    522                   //                  *fLog << warn <<
    523                   //                    "Overflow or Underflow occurred filling Blind Pixel eventnr = " << fEvents << endl;
    524526            } /* if use blind pixel */
    525527         
    … …  
    684686        }
    685687
     688      blindpixel.CheckOscillations();
     689
    686690      fCalibrations->SetBlindPixelMethodValid(kTRUE);
    687691      blindpixel.DrawClone();

  • trunk/MagicSoft/Mars/mcalib/MCalibrationPix.cc

    r3056 r3061  
    350350    CalcFFactorMethod();
    351351 
    352   return (fRSigmaCharge/fCharge)*TMath::Sqrt(fPheFFactorMethod);
     352  if (fPheFFactorMethod > 0)
     353    return (fRSigmaCharge/fCharge)*TMath::Sqrt(fPheFFactorMethod);
     354  else
     355    return -1.;
    353356}
    354357
    … …  
    891894 
    892895}
    893 

  • trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc

    r3026 r3061  
    6363#include "MHCalibrationConfig.h"
    6464
     65
    6566#include <TStyle.h>
    6667#include <TCanvas.h>
    6768#include <TPaveText.h>
    6869
     70#include <TGraph.h>
    6971#include <TF1.h>
    7072#include <TH1.h>
    7173#include <TRandom.h>
    7274
     75#include "MFFT.h"
    7376#include "MLog.h"
    7477#include "MLogManip.h"
    … …  
    8083const Int_t    MHCalibrationBlindPixel::fgBlindPixelChargeNbins        = 1000;
    8184const Int_t    MHCalibrationBlindPixel::fgBlindPixelTimeNbins          = 22;
    82 const Int_t    MHCalibrationBlindPixel::fgBlindPixelChargevsNbins      = 10000;
    8385const Axis_t   MHCalibrationBlindPixel::fgBlindPixelTimeFirst          = -9.00;
    8486const Axis_t   MHCalibrationBlindPixel::fgBlindPixelTimeLast           = 12.00;
    8587const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmp      = 0.008;
    8688const Double_t MHCalibrationBlindPixel::fgBlindPixelElectronicAmpError = 0.002;
     89 
     90const Axis_t  MHCalibrationBlindPixel::fNyquistFreq       = 1.0;
     91const Axis_t  MHCalibrationBlindPixel::fMinFreq           = 0.;
     92const Int_t   MHCalibrationBlindPixel::fPSDNbins          = 30;
    8793 
    8894// --------------------------------------------------------------------------
    … …  
    95101        fTimeGausFit(NULL),
    96102        fSinglePhePedFit(NULL),
     103        fPSDHiGain(NULL),
     104        fPSDLoGain(NULL),
     105        fHPSD(NULL),
     106        fPSDExpFit(NULL),
     107        fChargeXaxis(NULL),
     108        fPSDXaxis(NULL),
     109        fCurrentSize(1024),
    97110        fFitLegend(NULL)
    98111{
    … …  
    119132    fHBlindPixelTime->SetDirectory(NULL);
    120133
    121     fHBlindPixelChargevsN = new TH1I("HBlindPixelChargevsN","Sum of Charges vs. Event Number",
    122                                      fgBlindPixelChargevsNbins,-0.5,(Axis_t)fgBlindPixelChargevsNbins-0.5);
    123     fHBlindPixelChargevsN->SetXTitle("Event Nr.");
    124     fHBlindPixelChargevsN->SetYTitle("Sum of FADC slices");
    125     fHBlindPixelChargevsN->SetDirectory(NULL);
     134    fHiGains = new TArrayF(fCurrentSize);
     135    fLoGains = new TArrayF(fCurrentSize);
    126136
    127137    Clear();
    … …  
    136146  delete fHBlindPixelCharge;
    137147  delete fHBlindPixelTime;
    138   delete fHBlindPixelChargevsN;
     148
     149  delete fHiGains;
     150  delete fLoGains;
    139151
    140152  if (fHBlindPixelPSD)
    … …  
    146158  if(fSinglePhePedFit)
    147159    delete fSinglePhePedFit;
     160  if (fPSDExpFit)
     161    delete fPSDExpFit;
     162  if (fHPSD)
     163    delete fHPSD;
     164  if (fChargeXaxis)
     165    delete fChargeXaxis;
     166  if (fPSDXaxis)
     167    delete fPSDXaxis;
    148168
    149169}
    … …  
    151171void MHCalibrationBlindPixel::Clear(Option_t *o)
    152172{
     173
     174  fTotalEntries            = 0;
     175  fCurrentSize             = 1024;
    153176 
    154177  fBlindPixelChargefirst = -200.;
    … …  
    196219  if(fSinglePhePedFit)
    197220    delete fSinglePhePedFit;
     221  if (fPSDExpFit)
     222    delete fPSDExpFit;
     223  if (fHPSD)
     224    delete fHPSD;
     225  if (fChargeXaxis)
     226    delete fChargeXaxis;
     227  if (fPSDXaxis)
     228    delete fPSDXaxis;
     229  if (fPSDHiGain)
     230    delete fPSDHiGain;
     231  if (fPSDLoGain)
     232    delete fPSDLoGain;
     233
     234  CLRBIT(fFlags,kFitOK);
     235  CLRBIT(fFlags,kOscillating);
    198236
    199237  return;
    … …  
    207245  fHBlindPixelCharge->Reset();
    208246  fHBlindPixelTime->Reset();
    209   fHBlindPixelChargevsN->Reset();
    210  
    211 }
     247 
     248  fHiGains->Set(1024);
     249  fLoGains->Set(1024);
     250
     251  fHiGains->Reset(0.);
     252  fLoGains->Reset(0.);
     253
     254
     255}
     256
     257Bool_t MHCalibrationBlindPixel::CheckOscillations()
     258{
     259
     260  if (fPSDExpFit)
     261    return IsOscillating();
     262
     263  MFFT fourier;
     264
     265  fPSDLoGain = fourier.PowerSpectrumDensity(fLoGains);
     266  fPSDHiGain = fourier.PowerSpectrumDensity(fHiGains);
     267
     268  Int_t entries;
     269  TArrayF *array;
     270 
     271  fHPSD = new TH1F("HBlindPixelPSD",
     272                   "Power Spectrum Density Projection Blind Pixel",
     273                   fPSDNbins,fMinFreq,fNyquistFreq);
     274
     275  array   = fPSDHiGain;
     276  entries = array->GetSize();
     277 
     278  for (Int_t i=0;i<entries;i++)
     279    fHPSD->Fill(array->At(i));
     280
     281  //
     282  // First guesses for the fit (should be as close to reality as possible,
     283  //
     284  const Double_t area_guess  = entries*10.;
     285
     286  fPSDExpFit = new TF1("PSDExpFit","[0]*exp(-[1]*x)",0.,1.);
     287
     288  fPSDExpFit->SetParameters(entries,10.);
     289  fPSDExpFit->SetParNames("Area","slope");
     290  fPSDExpFit->SetParLimits(0,0.,3.*area_guess);
     291  fPSDExpFit->SetParLimits(1,5.,20.);
     292  fPSDExpFit->SetRange(fMinFreq,fNyquistFreq);
     293
     294  fHPSD->Fit(fPSDExpFit,"RQL0");
     295 
     296  fPSDProb  = fPSDExpFit->GetProb();
     297
     298  if (fPSDProb < gkProbLimit)
     299    {
     300      SETBIT(fFlags,kOscillating);
     301      return kTRUE;
     302    }
     303 
     304  CLRBIT(fFlags,kOscillating);
     305
     306  return kFALSE;
     307}
     308
     309void MHCalibrationBlindPixel::CreatePSDXaxis(Int_t n)
     310{
     311 
     312  if (fPSDXaxis)
     313    return;
     314
     315  fPSDXaxis = new TArrayF(n);
     316
     317  for (Int_t i=0;i<n;i++)
     318    fPSDXaxis->AddAt((Float_t)i,i);
     319}
     320
     321void MHCalibrationBlindPixel::CreateChargeXaxis(Int_t n)
     322{
     323 
     324  if (!fChargeXaxis)
     325    {
     326      fChargeXaxis = new TArrayF(n);
     327      for (Int_t i=0;i<n;i++)
     328        fChargeXaxis->AddAt((Float_t)i,i);
     329      return;
     330    }
     331
     332  if (fChargeXaxis->GetSize() == n)
     333    return;
     334
     335  const Int_t diff = fChargeXaxis->GetSize()-n;
     336  fChargeXaxis->Set(n);
     337  if (diff < 0)
     338    for (Int_t i=n;i<n+diff;i++)
     339      fChargeXaxis->AddAt((Float_t)i,i);
     340}
     341
     342void MHCalibrationBlindPixel::CutArrayBorder(TArrayF *array)
     343{
     344 
     345  Int_t i;
     346
     347  for (i=array->GetSize()-1;i>=0;i--)
     348    if (array->At(i) != 0)
     349      {
     350        array->Set(i+1);
     351        break;
     352      }
     353}
     354
     355const Bool_t MHCalibrationBlindPixel::IsFitOK() const
     356{
     357  return TESTBIT(fFlags,kFitOK);
     358}
     359
     360const Bool_t MHCalibrationBlindPixel::IsOscillating()
     361{
     362
     363  if (fPSDExpFit)
     364    return TESTBIT(fFlags,kOscillating);
     365
     366  return CheckOscillations();
     367
     368}
     369
     370Bool_t MHCalibrationBlindPixel::FillGraphs(Float_t qhi,Float_t qlo)
     371{
     372
     373  if (fTotalEntries >= fCurrentSize)
     374    {
     375      fCurrentSize *= 2;
     376
     377      fHiGains->Set(fCurrentSize);
     378      fLoGains->Set(fCurrentSize);
     379    }
     380 
     381  fHiGains->AddAt(qhi,fTotalEntries);
     382  fLoGains->AddAt(qlo,fTotalEntries);
     383
     384  fTotalEntries++;
     385
     386  return kTRUE;
     387
     388}
     389
     390
    212391
    213392Bool_t MHCalibrationBlindPixel::FillBlindPixelCharge(Float_t q)
    … …  
    221400}
    222401
    223 Bool_t MHCalibrationBlindPixel::FillBlindPixelChargevsN(Stat_t rq, Int_t t)
    224 {
    225     return fHBlindPixelChargevsN->Fill(t,rq) > -1;
    226 }
    227 
    228402
    229403// -------------------------------------------------------------------------
    … …  
    236410  fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);
    237411
    238   if (fFitOK)
     412  if (IsFitOK())
    239413      fFitLegend->SetFillColor(80);
    240414  else
    … …  
    284458  t8->SetBit(kCanDelete);
    285459
    286   if (fFitOK)
     460  if (IsFitOK())
    287461    {
    288462      TText *t = fFitLegend->AddText(0.,0.,"Result of the Fit: OK");
    … …  
    333507    TCanvas *c = MakeDefCanvas(this,550,700);
    334508
    335     c->Divide(2,3);
     509    c->Divide(2,4);
    336510
    337511    gROOT->SetSelectedPad(NULL);
    … …  
    344518    fHBlindPixelCharge->Draw(opt);
    345519   
    346     if (fFitOK)
     520    if (IsFitOK())
    347521      fSinglePheFit->SetLineColor(kGreen);         
    348522    else
    … …  
    366540    fHBlindPixelTime->Draw(opt);
    367541   
    368     c->cd(4);
    369 
    370     fHBlindPixelChargevsN->Draw(opt);
    371 
     542    CutArrayBorder(fHiGains);
     543    CreateChargeXaxis(fHiGains->GetSize());
     544   
     545    c->cd(5);
     546    gPad->SetTicks();
     547    TGraph *gr1 = new TGraph(fChargeXaxis->GetSize(),
     548                             fChargeXaxis->GetArray(),
     549                             fHiGains->GetArray()); 
     550    gr1->SetTitle("Evolution of HiGain charges with event number");
     551    gr1->SetBit(kCanDelete);
     552    gr1->Draw("AL");
     553   
     554    CutArrayBorder(fLoGains); 
     555    CreateChargeXaxis(fHiGains->GetSize());
     556   
     557    c->cd(6);
     558    gPad->SetTicks();
     559    TGraph *gr2 = new TGraph(fChargeXaxis->GetSize(),
     560                             fChargeXaxis->GetArray(),
     561                             fLoGains->GetArray()); 
     562    gr2->SetTitle("Evolution of HiGain charges with event number");
     563    gr2->SetBit(kCanDelete);
     564    gr2->Draw("AL");
     565 
    372566    c->Modified();
    373567    c->Update();
    374 
    375     c->cd(5);
    376    
    377     //    fHBlindPixelPSD = (TH1F*)MFFT::PowerSpectrumDensity(fHBlindPixelChargevsN);
    378     //    TH1F *hist = MFFT::PowerSpectrumDensity((*this)->fHBlindPixelChargevsN);
    379 
    380     //    fHBlindPixelPSD->Draw(opt);
     568   
     569    c->cd(7);
     570 
     571    TArrayF *array;
     572   
     573    if (!fPSDHiGain)
     574      return;
     575    array = fPSDHiGain;
     576 
     577    if (!fPSDXaxis)
     578      CreatePSDXaxis(array->GetSize());
     579
     580    TGraph *gr3 = new TGraph(fPSDXaxis->GetSize(),fPSDXaxis->GetArray(),array->GetArray());
     581    gr3->SetTitle("Power Spectrum Density");
     582    gr3->SetBit(kCanDelete);
     583    gr3->Draw("AL");
     584   
    381585    c->Modified();
    382586    c->Update();
    383 
     587   
     588    c->cd(8);
     589   
     590    gStyle->SetOptStat(111111);
     591    gPad->SetTicks(); 
     592   
     593    if (fHPSD->Integral() > 0)
     594      gPad->SetLogy();
     595   
     596    fHPSD->Draw(opt);
     597   
     598    if (fPSDExpFit)
     599      {
     600        fPSDExpFit->SetLineColor(IsOscillating() ? kRed : kGreen);         
     601        fPSDExpFit->Draw("same");
     602      }
     603   
     604    c->Modified();
     605    c->Update();
     606   
    384607}
    385608
    … …  
    653876      *fLog << err << "ERROR: Fit Probability " << fProb
    654877            << " is smaller than the allowed value: " << gkProbLimit << endl;
    655       fFitOK = kFALSE;
     878      CLRBIT(fFlags,kFitOK);
    656879      return kFALSE;
    657880    }
    … …  
    663886      *fLog << err << "ERROR: Statistics is too low: Only " << contSinglePhe
    664887            << " in the Single Photo-Electron peak " << endl;
    665       fFitOK = kFALSE;
     888      CLRBIT(fFlags,kFitOK);     
    666889      return kFALSE;
    667890    }
    … …  
    669892    *fLog << inf << contSinglePhe << " in Single Photo-Electron peak " << endl;
    670893 
    671   fFitOK = kTRUE;
     894  SETBIT(fFlags,kFitOK);
    672895   
    673896  return kTRUE;
    … …  
    678901{
    679902
    680   Int_t nbins = 30;
     903  Int_t nbins = 20;
    681904
    682905  CutEdges(fHBlindPixelCharge,nbins);
    … …  
    684907  fBlindPixelChargefirst = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetFirst());
    685908  fBlindPixelChargelast  = fHBlindPixelCharge->GetBinLowEdge(fHBlindPixelCharge->GetXaxis()->GetLast())+fHBlindPixelCharge->GetBinWidth(0);
    686 
    687   CutEdges(fHBlindPixelChargevsN,0);
    688909
    689910}

  • trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h

    r2997 r3061  
    66#endif
    77
     8class TArrayF;
    89class TH1F;
    910class TH1I;
    … …  
    1920  static const Int_t    fgBlindPixelChargeNbins;
    2021  static const Int_t    fgBlindPixelTimeNbins;
    21   static const Int_t    fgBlindPixelChargevsNbins;
    2222  static const Axis_t   fgBlindPixelTimeFirst;
    2323  static const Axis_t   fgBlindPixelTimeLast;
    2424  static const Double_t fgBlindPixelElectronicAmp;
    2525  static const Double_t fgBlindPixelElectronicAmpError;
     26
     27  static const Axis_t  fNyquistFreq;
     28  static const Axis_t  fMinFreq;
     29  static const Int_t   fPSDNbins;
    2630 
    2731  TH1F* fHBlindPixelCharge;        // Histogram with the single Phe spectrum
    2832  TH1F* fHBlindPixelTime;          // Variance of summed FADC slices
    29   TH1I* fHBlindPixelChargevsN;     // Summed Charge vs. Event Nr.
    3033  TH1F* fHBlindPixelPSD;           // Power spectrum density of fHBlindPixelChargevsN
    3134 
    … …  
    3437  TF1 *fSinglePhePedFit;
    3538
     39  TArrayF* fPSDHiGain;           //-> Power spectrum density of fHiGains
     40  TArrayF* fPSDLoGain;           //-> Power spectrum density of fLoGains
     41
     42  TH1F* fHPSD;                   //->
     43  TF1*  fPSDExpFit;              //->
     44 
     45  TArrayF *fHiGains;             //->
     46  TArrayF *fLoGains;             //->
     47  TArrayF *fChargeXaxis;         //
     48  TArrayF *fPSDXaxis;            //
     49
     50  Float_t  fPSDProb;
     51
     52  Int_t fTotalEntries;           // Number of entries
     53  Int_t fCurrentSize;
     54 
    3655  Axis_t  fBlindPixelChargefirst;
    3756  Axis_t  fBlindPixelChargelast;
    3857 
    3958  void DrawLegend();
     59  void CreateChargeXaxis(Int_t n);
     60  void CreatePSDXaxis(Int_t n);
     61  void CutArrayBorder(TArrayF *array);
    4062
    4163  TPaveText *fFitLegend;                 
    42   Bool_t fFitOK; 
    4364 
    4465  Double_t  fLambda;
    … …  
    7091  Double_t  fSigmaPedestal;
    7192  Double_t  fSigmaPedestalErr;
     93
     94  Byte_t   fFlags;
     95
     96  enum     { kFitOK, kOscillating };
     97
    7298 
    7399public:
    … …  
    81107  Bool_t FillBlindPixelCharge(Float_t q);
    82108  Bool_t FillBlindPixelTime(Float_t t);
    83   Bool_t FillBlindPixelChargevsN(Stat_t rq, Int_t t);
     109  Bool_t FillGraphs(Float_t qhi, Float_t qlo);
    84110 
    85111  // Setters
    … …  
    113139  const Double_t GetSigmaTimeErr()   const { return fSigmaTimeErr;   }
    114140
    115   const Bool_t IsFitOK()             const { return fFitOK;          }
    116 
    117   const TH1I *GetHBlindPixelChargevsN() const  { return fHBlindPixelChargevsN; }
     141  const Bool_t IsFitOK()        const;
     142  const Bool_t IsOscillating();
     143 
    118144  const TH1F *GetHBlindPixelPSD()       const  { return fHBlindPixelPSD;       }
    119145 
    … …  
    122148  void Draw(Option_t *option="");
    123149
     150 
    124151  // Fits
    125152  enum FitFunc_t  { kEPoisson4, kEPoisson5, kEPoisson6, kEPoisson7, kEPolya, kEMichele };
    … …  
    140167  // Others
    141168  void CutAllEdges();
    142 
     169  Bool_t CheckOscillations();
     170
     171 
    143172private:
    144173

  • trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc

    r3056 r3061  
    7979        fChargeXaxis(NULL),
    8080        fPSDXaxis(NULL),
    81         fFitLegend(NULL),
    82         fCurrentSize(1024)
     81        fCurrentSize(1024),
     82        fFitLegend(NULL)
    8383{
    8484
    … …  
    271271                       fPSDNbins,fMinFreq,fNyquistFreq);
    272272     
    273       array = fPSDLoGain;
     273      array = fPSDHiGain;
    274274    }
    275275 
Note: See TracChangeset for help on using the changeset viewer.