Changeset 7043 for trunk/MagicSoft/Mars/msignal/MExtractPINDiode.cc

Timestamp:

05/17/05 12:08:31 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/msignal/MExtractPINDiode.cc (modified) (10 diffs)

trunk/MagicSoft/Mars/msignal/MExtractPINDiode.cc

@@ -25 +25 @@
 //////////////////////////////////////////////////////////////////////////////
 //
-//   MExtractPINDiode
+//  MExtractPINDiode
 //
 //  Extracts the signal from a fixed window in a given range.
 //
 //  Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast) 
-//  to modify the ranges. Ranges have to be an even number. In case of odd 
-//  ranges, the last slice will be reduced by one.
+//  to modify the ranges.
+//
 //  Defaults are: 
 // 
-//   fHiGainFirst =  fgHiGainFirst =  3 
-//   fHiGainLast  =  fgHiGainLast  =  14
-//   fLoGainFirst =  fgLoGainFirst =  3 
-//   fLoGainLast  =  fgLoGainLast  =  14
+//   fHiGainFirst =  fgHiGainFirst =  0 
+//   fHiGainLast  =  fgHiGainLast  =  29
+//   fLoGainFirst =  fgLoGainFirst =  0
+//   fLoGainLast  =  fgLoGainLast  =  0
+//
+//  The FADC slices are fit by a Gaussian around the pulse maximum. 
+//  The following figures show two typical (high-intensity and low-intensity) 
+//  pulses together with the applied fit:
+//
+//Begin_Html
+/*
+<img src="images/PINDiode_pulse_high.gif">
+<img src="images/PINDiode_pulse_low.gif">
+*/
+//End_Html
+//
+// The fit ranges can be modified with the functions:
+// - SetLowerFitLimit() 
+// - SetUpperFitLimit()
+//
+// Defaults are:
+//   - fLowerFitLimit: 2
+//   - fUpperFitLimit: 5
 //
 //////////////////////////////////////////////////////////////////////////////
@@ -44 +63 @@
 #include <fstream>
 
+#include <TF1.h>
+#include <TH1.h>
+#include <TPad.h>
+
 #include "MLog.h"
 #include "MLogManip.h"
@@ -61 +84 @@
 using namespace std;
 
-const UInt_t MExtractPINDiode::fgPINDiodeIdx     = 100;
-const Byte_t MExtractPINDiode::fgHiGainFirst =  0;
-const Byte_t MExtractPINDiode::fgHiGainLast  =  14;
-const Byte_t MExtractPINDiode::fgLoGainFirst =  0;
-const Byte_t MExtractPINDiode::fgLoGainLast  =  14;
+const UInt_t MExtractPINDiode::fgPINDiodeIdx   =  3;
+const Byte_t MExtractPINDiode::fgHiGainFirst   =  0;
+const Byte_t MExtractPINDiode::fgHiGainLast    = 29;
+const Byte_t MExtractPINDiode::fgLoGainFirst   =  0;
+const Byte_t MExtractPINDiode::fgLoGainLast    =  0;
+const Byte_t MExtractPINDiode::fgLowerFitLimit =  2;
+const Byte_t MExtractPINDiode::fgUpperFitLimit =  5;
+
 // --------------------------------------------------------------------------
 //
@@ -72 +98 @@
 // Calls: 
 // - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
-// - SetPINDiodeIdx()
+//
+// - Set fPINDiodeIdx   to fgPINDiodeIdx
+// - Set fLowerFitLimit to fgLowerFitLimit
+// - Set fUpperFitLimit to fgUpperFitLimit
 //
 MExtractPINDiode::MExtractPINDiode(const char *name, const char *title)
-{
-  
+    : fSlices(NULL)
+{
+
   fName  = name  ? name  : "MExtractPINDiode";
   fTitle = title ? title : "Task to extract the signal from the FADC slices";
@@ -82 +112 @@
   SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
   SetPINDiodeIdx();
+
+  SetLowerFitLimit();
+  SetUpperFitLimit();  
+
+  fPedMean.Set(2);
+}
+
+// --------------------------------------------------------------------------
+//
+// - delete the histogram fSlices, if it exists
+//
+MExtractPINDiode::~MExtractPINDiode()
+{
+  if (fSlices)
+    delete fSlices;
 }
 
@@ -89 +134 @@
 //
 // Checks: 
-// - if the window defined by (fHiGainLast-fHiGainFirst-1) are odd, subtract one
-// - if the window defined by (fLoGainLast-fLoGainFirst-1) are odd, subtract one
-// - if the Hi Gain window is smaller than 2, set fHiGainLast to fHiGainFirst+1
-// - if the Lo Gain window is smaller than 2, set fLoGainLast to fLoGainFirst+1
+// - if the Hi Gain window is smaller than 4, set fHiGainLast to fHiGainFirst+3
 // 
 // Calls:
@@ -99 +141 @@
 // Sets:
 // - fNumHiGainSamples to: (Float_t)(fHiGainLast-fHiGainFirst+1)
-// - fNumLoGainSamples to: (Float_t)(fLoGainLast-fLoGainFirst+1)
+// - fNumLoGainSamples to: 0.
 // - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
-// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)  
+// - fSqrtLoGainSamples to: 0.
 //  
 void MExtractPINDiode::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
 {
 
-  const Byte_t window = hilast-hifirst+1+lolast-lofirst+1;
-  const Byte_t weven  = window & ~1;
-
-  if (weven != window)
+  lofirst = 0;
+  lolast  = 0;
+
+  const Byte_t window = hilast-hifirst+1;
+
+  if (window<4) 
     {
       *fLog << warn << GetDescriptor() 
-            << Form("%s%2i%s%2i",": Total window size has to be even, set last slice from "
-                    ,(int)lolast," to ",(int)(lolast-1)) << endl;
-      lolast -= 1;
-    }
-  
-  if (window<2) 
-    {
-      *fLog << warn << GetDescriptor() 
-            << Form("%s%2i%s%2i",": Total window is smaller than 2 FADC sampes, set last slice from" 
-                    ,(int)lolast," to ",(int)(lofirst+1)) << endl;
-      hilast = hifirst+1;
-    }
-  
+            << Form("%s%2i%s%2i",": Total window is smaller than 4 FADC sampes, set last slice from" 
+                    ,(int)lolast," to ",(int)(lofirst+3)) << endl;
+      hilast = hifirst+3;
+    }
 
   MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
 
   fNumHiGainSamples = (Float_t)(fHiGainLast-fHiGainFirst+1);
-  fNumLoGainSamples = (fLoGainLast == 0) ? 0. : (Float_t)(fLoGainLast-fLoGainFirst+1);  
+  fNumLoGainSamples = 0.;
 
   fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
-  fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);  
-  
-  fNumSamples = (fLoGainLast == 0) 
-    ? fHiGainLast-fHiGainFirst+1 
-    : fHiGainLast-fHiGainFirst+1+fLoGainLast-fLoGainFirst+1;
-  fSqrtSamples = TMath::Sqrt((Float_t)fNumSamples);
+  fSqrtLoGainSamples = 0.;
   
 }
@@ -149 +179 @@
 // they were not found:
 //
+//  - MRawEvtData2
 //  - MExtractedPINDiode
+//
+// Initializes fPedMean to:
+// - fPedMean[0]: pedestal + AB-offset
+// - fPedMean[1]: pedestal - AB-offset
+//
+// Initializes TH1F fSlices to [fHiGainFirst-0.5,fHiGainLast+0.5]
 //
 Int_t MExtractPINDiode::PreProcess(MParList *pList)
@@ -157 +194 @@
     return kFALSE;
   
+  fRawEvt = NULL;
+  fRawEvt = (MRawEvtData*)pList->FindObject(AddSerialNumber("MRawEvtData2"));
+  if (!fRawEvt)
+    {
+      *fLog << err << AddSerialNumber("MRawEvtData2") << " not found... aborting." << endl;
+      return kFALSE;
+    }
+
   fPINDiode = (MExtractedSignalPINDiode*)pList->FindCreateObj(AddSerialNumber("MExtractedSignalPINDiode"));
   if (!fPINDiode)
     return kFALSE;
 
+  fPedMean.Reset();
+
   const MPedestalPix &ped   = (*fPedestals)[fPINDiodeIdx]; 
 
-    if (&ped)
-      {
-        fPedestal = ped.GetPedestal();
-        fPedRms   = ped.GetPedestalRms();
-      }
-    else
-      {
-        *fLog << err << " Cannot find MPedestalPix of the PIN Diode (idx=" 
-              << fPINDiodeIdx << ")" << endl;
+  if (&ped)
+    {
+      fPedMean[0] = ped.GetPedestal() + ped.GetPedestalABoffset();
+      fPedMean[1] = ped.GetPedestal() - ped.GetPedestalABoffset();      
+    }
+  else
+    {
+      *fLog << err << " Cannot find MPedestalPix of the PIN Diode (idx=" 
+            << fPINDiodeIdx << ")" << endl;
         return kFALSE;
-      }
-
-    return kTRUE;
+    }
+  
+  if (fSlices)
+    delete fSlices;
+
+  fSlices = new TH1F("PINDiode","PIN Diode fitted slices",(Int_t)(fHiGainLast-fHiGainFirst+1),
+                     fHiGainFirst-0.5,fHiGainLast+0.5);
+  fSlices->SetDirectory(NULL);
+
+  return kTRUE;
 }
 
@@ -190 +244 @@
   
   fPINDiode->SetUsedFADCSlices(fHiGainFirst, fLoGainLast);
-  
+
+  *fLog << endl;
+  *fLog << inf << "Taking " << fNumHiGainSamples
+        << " HiGain samples from slice " << (Int_t)fHiGainFirst
+        << " to " << (Int_t)(fHiGainLast+fHiLoLast) << " incl" << endl;
+
   return kTRUE;
-  
-}
-
-
-
-void MExtractPINDiode::FindSignalandVarianceHiGain(Byte_t *ptr, Int_t &sum, Int_t &sum2, Byte_t &sat) const
-{
-
-  Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1;
-
+}
+
+
+
+
+// ----------------------------------------------------------------------------------
+//
+// Extracts the (pedestal-subtracted) FADC slices between fHiGainFirst and fHiGainLast 
+// and fills them into the histogram fSlices. Memorizes the position of maximum at 
+// maxpos.
+//
+// Checks for saturation
+// 
+// Fits fSlices to a Gaussian in the ranges: maxpos-fLowerFitLimit, maxpos+fUpperFitLimit
+//
+// Writes fit results into MExtractedSignalPINDiode
+//
+Int_t MExtractPINDiode::Process()
+{
+
+
+  MRawEvtPixelIter pixel(fRawEvt);
+  
+  fPINDiode->Clear();
+  fSlices->Reset();
+  
+  pixel.Jump(fPINDiodeIdx);
+  
+  Byte_t sat  = 0;
+
+  Int_t higainsamples = pixel.GetNumHiGainSamples();
+  Int_t logainsamples = pixel.GetNumLoGainSamples();
+  
+  const Bool_t higainabflag = pixel.HasABFlag();
+  Byte_t *ptr = pixel.GetHiGainSamples()+fHiGainFirst;
+  Byte_t *end = ptr+higainsamples;
+  
+  Int_t cnt=0;
+
+  Float_t max = 0.;
+  Int_t maxpos = 0;
+  
   while (ptr<end)
     {
-      sum  += *ptr;
-      sum2 += *ptr * *ptr;
-
-      if (*ptr++ >= fSaturationLimit)
-        sat++;
-    }
-}
-
-void MExtractPINDiode::FindSignalandVarianceLoGain(Byte_t *ptr, Int_t &sum, Int_t &sum2, Byte_t &sat) const
-{
-
-  Byte_t *end = ptr +  fLoGainLast - fLoGainFirst  + 1;
-
-  while (ptr<end)
-    {
-      sum  += *ptr;
-      sum2 += *ptr * *ptr;
       
-      if (*ptr++ >= fSaturationLimit)
-        sat++;
-    }
-}
-
-
-
-// --------------------------------------------------------------------------
-//
-// Calculate the integral of the FADC time slices and store them as a new
-// pixel in the MExtractedPINDiode container.
-//
-Int_t MExtractPINDiode::Process()
-{
-
-
-  MRawEvtPixelIter pixel(fRawEvt);
-  
-  fPINDiode->Clear();
-  
-  pixel.Jump(fPINDiodeIdx);
-  
-  Int_t sum   = 0;
-  Int_t sum2  = 0;
-  Byte_t sat  = 0;
-  Int_t max   = 0;
-  
-  // 
-  // Calculate first the time:
-  //
-  const Int_t maxhi = pixel.GetIdxMaxHiGainSample();
-  const Int_t maxlo = pixel.GetIdxMaxLoGainSample();
-  
-  if (maxhi > maxlo)
-    max = maxhi;
-  else
-    max = maxlo + pixel.GetNumHiGainSamples();
-  
-  FindSignalandVarianceHiGain(pixel.GetHiGainSamples()+fHiGainFirst,sum,sum2,sat);
-  if (pixel.HasLoGain())
-    FindSignalandVarianceLoGain(pixel.GetLoGainSamples()+fLoGainFirst,sum,sum2,sat);
-
-  const Float_t var = ((Float_t)sum2 - (Float_t)sum*sum/fNumSamples)/(fNumSamples-1);
-  const Float_t rms = TMath::Sqrt(var);
-  
-  // 
-  // FIXME: The following formulae have to be revised!!
-  //
-  fPINDiode->SetExtractedSignal(sum - fPedestal*fNumSamples, fPedRms*fSqrtSamples);
-  fPINDiode->SetExtractedRms   (rms, rms/2./fSqrtSamples);
-  fPINDiode->SetExtractedTime  (max, rms/fSqrtSamples);
-  fPINDiode->SetSaturation(sat);
+      if (*ptr >= fSaturationLimit)
+        {
+          sat++;
+          break;
+        }
+
+      const Float_t cont = (Float_t)*ptr - fPedMean[(cnt + higainabflag) & 0x1];
+      fSlices->Fill(cnt,cont);
+
+      if (cont > max)
+      {
+        max = cont;
+        maxpos = cnt;
+      }
+
+      ptr++;
+      cnt++;
+    }
+  
+  cnt = 0;
+  
+  if (pixel.HasLoGain() && !sat)
+    {
+      
+      ptr = pixel.GetLoGainSamples();
+      end = ptr+logainsamples;
+      
+      const Bool_t logainabflag = (higainabflag + pixel.GetNumHiGainSamples()) & 0x1;
+      
+      while (ptr<end)
+        {
+          
+          if (*ptr >= fSaturationLimit)
+            {
+              sat++;
+              break;
+            }
+
+          const Float_t cont = (Float_t)*ptr - fPedMean[(cnt + logainabflag) & 0x1];
+          
+          fSlices->Fill(cnt+ higainsamples,cont);
+          
+          if (cont > max) 
+            {
+              max    = cont;
+              maxpos = cnt+higainsamples;
+            }
+          ptr++;
+          cnt++;
+        }
+    }
+          
+  if (sat)
+    {
+      fPINDiode->SetSaturation(1);
+      return kTRUE;
+    }
+
+  fSlices->Fit("gaus", "RQ", "", maxpos-fLowerFitLimit,maxpos+fUpperFitLimit);
+  
+  TF1 *gausfunc = fSlices->GetFunction("gaus");
+
+  fPINDiode->SetExtractedSignal(gausfunc->GetParameter(0), gausfunc->GetParError(0));
+  fPINDiode->SetExtractedTime  (gausfunc->GetParameter(1), gausfunc->GetParError(1));
+  fPINDiode->SetExtractedSigma (gausfunc->GetParameter(2), gausfunc->GetParError(2));
+  fPINDiode->SetExtractedChi2  (gausfunc->GetChisquare());
   fPINDiode->SetReadyToSave();
-  
+
   return kTRUE;
 }
 
+// ----------------------------------------------------------------------------------
+//
+// deletes fSlices and sets pointer to NULL
+//
+Int_t MExtractPINDiode::PostProcess()
+{
+  
+  delete fSlices;
+  fSlices = NULL;
+  
+  return kTRUE;
+}