Changeset 3081 for trunk/MagicSoft

Timestamp:

02/10/04 16:10:52 (21 years ago)

Author:

gaug

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• mcalib/MHCalibrationBlindPixel.cc (modified) (2 diffs)
• mcalib/MHCalibrationBlindPixel.h (modified) (2 diffs)
• mcalib/MHCalibrationPixel.cc (modified) (2 diffs)
• mcalib/MHCalibrationPixel.h (modified) (2 diffs)
• mhbase/MH.cc (modified) (1 diff)
• mhbase/MH.h (modified) (2 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -8 +8 @@
   * mcalib/MCalibrationCam.cc
     - fixed documentation
+
+  * mhbase/MH.[h,cc]
+    - new function ProjectArray
+
+  * mcalib/MHCalibrationPixel.[h,cc]
+  * mcalib/MHCalibrationBlindPixel.[h,cc]
+    - use ProjectArray from MH to plot the projection of the fourier 
+      spectrum
 
 

trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.cc

@@ -88 +88 @@
 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;
   
@@ -266 +264 @@
   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));
-
+  fHPSD = MH::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 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,0.,20.);
-  fPSDExpFit->SetRange(fMinFreq,fNyquistFreq);
+  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/2.,2.*slope_guess);
+  fPSDExpFit->SetRange(0.,xmax);
 
   fHPSD->Fit(fPSDExpFit,"RQL0");

trunk/MagicSoft/Mars/mcalib/MHCalibrationBlindPixel.h

@@ -25 +25 @@
   static const Double_t fgBlindPixelElectronicAmpError;
 
-  static const Axis_t  fNyquistFreq;
-  static const Axis_t  fMinFreq;
   static const Int_t   fPSDNbins;
   
@@ -40 +38 @@
   TArrayF* fPSDLoGain;           //-> Power spectrum density of fLoGains
 
-  TH1F* fHPSD;                   //-> 
+  TH1I* fHPSD;                   //-> 
   TF1*  fPSDExpFit;              //->
   

trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.cc

@@ -61 +61 @@
 const Int_t   MHCalibrationPixel::fNDFLimit          = 5;
 
-const Axis_t  MHCalibrationPixel::fMinFreq           = 0.;
-const Int_t   MHCalibrationPixel::fPSDNbins          = 50;
+const Int_t   MHCalibrationPixel::fPSDNbins          = 30;
   
 // --------------------------------------------------------------------------
@@ -268 +267 @@
   fPSDHiGain = fourier.PowerSpectrumDensity(fHiGains);
 
-  Int_t entries;
-  TArrayF *array;
-
-  const Double_t max_freq = (fChargeSigma > 0.) ? fChargeSigma*0.2 : 2.;
-  const Double_t min_freq = 0.;
-
   if (IsUseLoGain())
-    {
-      fHPSD = new TH1F(Form("%s%d","HPSD",fPixId),
-                       Form("%s%s","Power Spectrum Density Projection ","LoGain"),
-                       fPSDNbins,min_freq,max_freq);
-      
-      array = fPSDLoGain;
-    }
+    fHPSD = MH::ProjectArray(fPSDLoGain, fPSDNbins, 
+                             Form("%s%d","PSDProjection",fPixId),
+                             Form("%s%d","Power Spectrum Density Projection LoGain ",fPixId));
   else
-    {
-
-      fHPSD = new TH1F(Form("%s%d","HPSD",fPixId),
-                       Form("%s%s","Power Spectrum Density Projection ","HiGain"),
-                       fPSDNbins,min_freq,max_freq);
-      
-      array = fPSDHiGain;
-    }
-  
-  entries = array->GetSize();
-  
-  for (Int_t i=0;i<entries;i++)
-    fHPSD->Fill(array->At(i));
+    fHPSD = MH::ProjectArray(fPSDHiGain, fPSDNbins, 
+                             Form("%s%d","PSDProjection",fPixId),
+                             Form("%s%d","Power Spectrum Density Projection HiGain ",fPixId));
 
   //
   // First guesses for the fit (should be as close to reality as possible, 
   //
-  const Double_t area_guess  = entries*10.;
-
-  fPSDExpFit = new TF1(Form("%s%d","PSDExpFit",fPixId),"[0]*exp(-[1]*x)",0.,1.);
-
-  fPSDExpFit->SetParameters(area_guess,10.);
-  fPSDExpFit->SetParNames("Area","slope");
-  fPSDExpFit->SetParLimits(0,0.,20.*area_guess);
-  fPSDExpFit->SetParLimits(1,0.,30.);
-  fPSDExpFit->SetRange(min_freq,max_freq);
+  const Double_t xmax = fHPSD->GetXaxis()->GetXmax();
+
+  fPSDExpFit = new TF1(Form("%s%d","PSDExpFit",fPixId),"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/2.,2.*slope_guess);
+  fPSDExpFit->SetRange(0.,xmax);
 
   fHPSD->Fit(fPSDExpFit,"RQL0");

trunk/MagicSoft/Mars/mcalib/MHCalibrationPixel.h

@@ -26 +26 @@
   static const Int_t   fNDFLimit;  
 
-  static const Axis_t  fMinFreq;
   static const Int_t   fPSDNbins;
   
@@ -49 +48 @@
   TF1* fChargeGausFit;           //->
 
-  TH1F* fHPSD;                   //-> 
+  TH1I* fHPSD;                   //-> 
   TF1*  fPSDExpFit;              //->
   

trunk/MagicSoft/Mars/mhbase/MH.cc

@@ -1081 +1081 @@
     return NULL;
 }
+
+TH1I* MH::ProjectArray(const TArrayF *array, Int_t nbins, const char* name, const char* title)
+{
+
+  const Int_t size  = array->GetSize();
+  TH1I *h1=0;
+
+  //check if histogram with identical name exist
+  TObject *h1obj = gROOT->FindObject(name);
+  if (h1obj && h1obj->InheritsFrom("TH1I")) {
+    h1 = (TH1I*)h1obj;
+    h1->Reset();
+  }
+
+  Double_t max = 0.;
+  Double_t min = 0.;
+  Int_t newbins = 0;
+
+  // first loop over array to find the maximum:
+  for (Int_t i=0; i<size;i++)
+    if (array->At(i) > max)
+      max = array->At(i);
+
+  FindGoodLimits(nbins, newbins, min, max, kFALSE);
+
+  if (!h1)
+    {
+      h1 = new TH1I(name, title, newbins, min, max);
+      h1->SetXTitle("");
+      h1->SetYTitle("Counts");
+      h1->SetDirectory(NULL);
+    }
+  
+    // Second loop to fill the histogram
+    for (Int_t i=0;i<size;i++) 
+      h1->Fill(array->At(i));
+
+    h1->SetEntries(size);
+    
+    return h1;
+}
+
+TH1I* MH::ProjectArray(const TArrayD *array, Int_t nbins, const char* name, const char* title)
+{
+
+  const Int_t size  = array->GetSize();
+  
+  TH1I *h1=0;
+
+  //check if histogram with identical name exist
+  TObject *h1obj = gROOT->FindObject(name);
+  if (h1obj && h1obj->InheritsFrom("TH1I")) {
+    h1 = (TH1I*)h1obj;
+    h1->Reset();
+  }
+
+  Double_t max = 0.;
+  Double_t min = 0.;
+  Int_t newbins = 0;
+
+  // first loop over array to find the maximum:
+  for (Int_t i=0; i<size;i++)
+    if (array->At(i) > max)
+      max = array->At(i);
+
+  FindGoodLimits( nbins, newbins, min, max, kFALSE);
+
+  if (!h1)
+    {
+      h1 = new TH1I(name, title, newbins, min, max);
+      h1->SetXTitle("");
+      h1->SetYTitle("Counts");
+      h1->SetDirectory(NULL);
+    }
+  
+    // Second loop to fill the histogram
+    for (Int_t i=0;i<size;i++) 
+      h1->Fill(array->At(i));
+
+    h1->SetEntries(size);
+    
+    return h1;
+}
+

trunk/MagicSoft/Mars/mhbase/MH.h

@@ -7 +7 @@
 
 class TH1;
+class TH1I;
 class TH1D;
 class TH2;
 class TH3;
 class TAxis;
+class TArrayF;
 class TArrayD;
 class TCanvas;
@@ -89 +91 @@
     static Double_t GetMinimumGT(const TH1 &h, Double_t gt=0);
     static Int_t CutEdges(TH1 *h, Int_t nbins);
+
+    static TH1I* ProjectArray(const TArrayF *array, Int_t nbins=30, const char* name="ProjectArray",
+                                                              const char* title="Projected Array");
+    static TH1I* ProjectArray(const TArrayD *array, Int_t nbins=30, const char* name="ProjectArray",
+                                                              const char* title="Projected Array");
     
     static void ProjectionX(TH1D &dest, const TH2 &src, Int_t firstybin=-1, Int_t lastybin=9999);