Changeset 5221

Timestamp:

10/11/04 21:08:37 (20 years ago)

Author:

gaug

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• msignal/MExtractTimeAndChargeSpline.cc (modified) (21 diffs)
• msignal/MExtractTimeAndChargeSpline.h (modified) (3 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -26 +26 @@
        in MPedCalcFromLoGain. 
 
+   * msignal/MExtractTimeAndChargeSpline.[h,cc]
+     - fixed class documentation and some last bugs.
 
  2004/10/08: Markus Meyer and Keiichi Mase

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc

@@ -15 +15 @@
 ! *
 !
-!   Author(s): Markus Gaug       05/2004 <mailto:markus@ifae.es> 
+!   Author(s): Markus Gaug       09/2004 <mailto:markus@ifae.es> 
 !
 !   Copyright: MAGIC Software Development, 2002-2004
@@ -26 +26 @@
 //   MExtractTimeAndChargeSpline
 //
-//   Fast Spline extractor using a cubic spline algorithm of Numerical Recipes. 
-//   It returns the integral below the interpolating spline. 
+//   Fast Spline extractor using a cubic spline algorithm, adapted from 
+//   Numerical Recipes in C++, 2nd edition, pp. 116-119.
+//   
+//   The coefficients "ya" are here denoted as "fHiGainSignal" and "fLoGainSignal"
+//   which means the FADC value subtracted by the clock-noise corrected pedestal.
+//
+//   The coefficients "y2a" get immediately divided 6. and are called here 
+//   "fHiGainSecondDeriv" and "fLoGainSecondDeriv" although they are now not exactly 
+//   the second derivative coefficients any more. 
 // 
-//   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.
-//
-//  Defaults are: 
+//   The calculation of the cubic-spline interpolated value "y" on a point 
+//   "x" along the FADC-slices axis becomes:
 // 
-//   fHiGainFirst =  fgHiGainFirst =  3 
-//   fHiGainLast  =  fgHiGainLast  =  14
-//   fLoGainFirst =  fgLoGainFirst =  3 
-//   fLoGainLast  =  fgLoGainLast  =  14
-//
+//   y =    a*fHiGainSignal[klo] + b*fHiGainSignal[khi] 
+//       + (a*a*a-a)*fHiGainSecondDeriv[klo] + (b*b*b-b)*fHiGainSecondDeriv[khi]
+//
+//   with:
+//   a = (khi - x)
+//   b = (x - klo)
+//
+//   and "klo" being the lower bin edge FADC index and "khi" the upper bin edge FADC index.
+//   fHiGainSignal[klo] and fHiGainSignal[khi] are the FADC values at "klo" and "khi".
+//
+//   An analogues formula is used for the low-gain values.
+//
+//   The coefficients fHiGainSecondDeriv and fLoGainSecondDeriv are calculated with the 
+//   following simplified algorithm:
+//
+//   for (Int_t i=1;i<range-1;i++) {
+//       pp                   = fHiGainSecondDeriv[i-1] + 4.;
+//       fHiGainFirstDeriv[i] = fHiGainSignal[i+1] - 2.*fHiGainSignal[i] + fHiGainSignal[i-1]
+//       fHiGainFirstDeriv[i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
+//   }
+// 
+//   for (Int_t k=range-2;k>=0;k--)
+//       fHiGainSecondDeriv[k] = (fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k])/6.;
+// 
+//
+//   This algorithm takes advantage of the fact that the x-values are all separated by exactly 1
+//   which simplifies the Numerical Recipes algorithm.
+//   (Note that the variables "fHiGainFirstDeriv" are not real first derivative coefficients.)
+//
+//
+//   The algorithm to search the time proceeds as follows:
+//
+//   1) Calculate all fHiGainSignal from fHiGainFirst to fHiGainLast 
+//      (note that an "overlap" to the low-gain arrays is possible: i.e. fHiGainLast>14 in the case of 
+//      the MAGIC FADCs).
+//   2) Remember the position of the slice with the highest content "fAbMax" at "fAbMaxPos".
+//   3) If one or more slices are saturated or fAbMaxPos is less than 2 slices from fHiGainFirst, 
+//      return fAbMaxPos as time and fAbMax as charge (note that the pedestal is subtracted here).
+//   4) Calculate all fHiGainSecondDeriv from the fHiGainSignal array
+//   5) Search for the maximum, starting in interval fAbMaxPos-1 in steps of 0.2 till fAbMaxPos-0.2.
+//      If no maximum is found, go to interval fAbMaxPos+1. 
+//      --> 4 function evaluations
+//   6) Search for the absolute maximum from fAbMaxPos to fAbMaxPos+1 in steps of 0.2 
+//      --> 4 function  evaluations
+//   7) Try a better precision searching from new max. position fAbMaxPos-0.2 to fAbMaxPos+0.2
+//      in steps of 0.025 (83 psec. in the case of the MAGIC FADCs).
+//      --> 14 function evaluations
+//   8) If Time Extraction Type kMaximum has been chosen, the position of the found maximum is 
+//      returned, else:
+//   9) The Half Maximum is calculated. 
+//  10) fHiGainSignal is called beginning from fAbMaxPos-1 backwards until a value smaller than fHalfMax
+//      is found at "klo". 
+//  11) Then, the spline value between "klo" and "klo"+1 is halfed by means of bisection as long as 
+//      the difference between fHalfMax and spline evaluation is less than fResolution (default: 0.01).
+//      --> maximum 12 interations.
+//   
+//  The algorithm to search the charge proceeds as follows:
+//
+//  1) If Charge Type: kAmplitude was chosen, return the Maximum of the spline, found during the 
+//     time search.
+//  2) If Charge Type: kIntegral was chosen, sum the fHiGainSignal between:
+//     (Int_t)(fAbMaxPos - fRiseTime) and 
+//     (Int_t)(fAbMaxPos + fFallTime)
+//     (default: fRiseTime: 1.5, fFallTime: 4.5)
+//  3) Sum only half the values of the edge slices
+//  4) Sum 1.5*fHiGainSecondDeriv of the not-edge slices using the "natural cubic 
+//     spline with second derivatives set to 0. at the edges.
+//     (Remember that fHiGainSecondDeriv had been divided by 6.)
+//   
+//  The values: fNumHiGainSamples and fNumLoGainSamples are set to:
+//  1) If Charge Type: kAmplitude was chosen: 1.
+//  2) If Charge Type: kIntegral was chosen: TMath::Floor(fRiseTime + fFallTime)
+//
+//  Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast) 
+//        to modify the ranges.
+// 
+//        Defaults: 
+//        fHiGainFirst =  2 
+//        fHiGainLast  =  14
+//        fLoGainFirst =  3 
+//        fLoGainLast  =  14
+//
+//  Call: SetResolution() to define the resolution of the half-maximum search.
+//        Default: 0.01
+//
+//  Call: SetRiseTime() and SetFallTime() to define the integration ranges 
+//        for the case, the extraction type kIntegral has been chosen.
+//
+//  Call: - SetTimeType(MExtractTimeAndChargeSpline::kMaximum) for extraction 
+//          the position of the maximum (default)
+//          --> needs 22 function evaluations
+//        - SetTimeType(MExtractTimeAndChargeSpline::kHalfMaximum) for extraction 
+//          the position of the half maximum at the rising edge.
+//          --> needs max. 34 function evaluations
+//        - SetChargeType(MExtractTimeAndChargeSpline::kAmplitude) for the 
+//          computation of the amplitude at the maximum (default)
+//          --> no further function evaluation needed
+//        - SetChargeType(MExtractTimeAndChargeSpline::kIntegral) for the 
+//          computation of the integral beneith the spline between fRiseTime 
+//          from the position of the maximum to fFallTime after the position of 
+//          the maximum.
+//          --> needs one more simple summation loop over 7 slices.
+//        
 //////////////////////////////////////////////////////////////////////////////
 #include "MExtractTimeAndChargeSpline.h"
@@ -60 +162 @@
 const Byte_t  MExtractTimeAndChargeSpline::fgLoGainFirst  = 3;
 const Byte_t  MExtractTimeAndChargeSpline::fgLoGainLast   = 14;
-const Float_t MExtractTimeAndChargeSpline::fgResolution   = 0.001;
+const Float_t MExtractTimeAndChargeSpline::fgResolution   = 0.01;
 const Float_t MExtractTimeAndChargeSpline::fgRiseTime     = 2.0;
 const Float_t MExtractTimeAndChargeSpline::fgFallTime     = 4.0;
@@ -89 +191 @@
   SetFallTime();
   
+  SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
+
   SetTimeType();
   SetChargeType();
 
-  SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
-
-  fNumHiGainSamples = 1.;
-  fNumLoGainSamples = 1.;
-  fSqrtHiGainSamples = 1.;
-  fSqrtLoGainSamples = 1.;
-  
 }
 
@@ -119 +216 @@
 }
 
+//-------------------------------------------------------------------
+//
+// Set the ranges
+// In order to set the fNum...Samples variables correctly for the case, 
+// the integral is computed, have to overwrite this function and make an 
+// explicit call to SetChargeType().
+//
+void MExtractTimeAndChargeSpline::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
+{
+
+  MExtractor::SetRange(hifirst, hilast, lofirst, lolast);
+
+  if (IsExtractionType(kIntegral))
+    SetChargeType(kIntegral);
+  
+}
+
+
+//-------------------------------------------------------------------
+//
+// Set the Time Extraction type. Possible are:
+// - kMaximum:     Search the maximum of the spline and return its position
+// - kHalfMaximum: Search the half maximum left from the maximum and return
+//                 its position
+//
+void MExtractTimeAndChargeSpline::SetTimeType( ExtractionType_t typ )
+{
+
+  CLRBIT(fFlags,kMaximum);
+  CLRBIT(fFlags,kHalfMaximum);
+  SETBIT(fFlags,typ);
+
+}
+
+//-------------------------------------------------------------------
+//
+// Set the Charge Extraction type. Possible are:
+// - kAmplitude: Search the value of the spline at the maximum
+// - kIntegral:  Integral the spline from fHiGainFirst to fHiGainLast,   
+//               by counting the edge bins only half and setting the 
+//               second derivative to zero, there.
+//
+void MExtractTimeAndChargeSpline::SetChargeType( ExtractionType_t typ )
+{
+
+  CLRBIT(fFlags,kAmplitude);
+  CLRBIT(fFlags,kIntegral );
+
+  SETBIT(fFlags,typ);
+
+  if (IsExtractionType(kAmplitude))
+    {
+      fNumHiGainSamples = 1.;
+      fNumLoGainSamples = 1.;
+      fSqrtHiGainSamples = 1.;
+      fSqrtLoGainSamples = 1.;
+    }
+
+  if (IsExtractionType(kIntegral))
+    {
+      fNumHiGainSamples  = TMath::Floor(fRiseTime + fFallTime);
+      fNumLoGainSamples  = fNumHiGainSamples;
+      fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
+      fSqrtLoGainSamples = fSqrtHiGainSamples;
+    }
+}
+
+
 // --------------------------------------------------------------------------
 //
@@ -275 +440 @@
   if (sat || maxpos < 2)
     {
-      time =  IsTimeType(kMaximum) 
+      time =  IsExtractionType(kMaximum) 
         ? (Float_t)(fHiGainFirst + maxpos) 
         : (Float_t)(fHiGainFirst + maxpos - 1);
+      sum  =  IsExtractionType(kAmplitude) 
+        ? fAbMax : 0.;
       return;
     }
@@ -295 +462 @@
 
   fHiGainSecondDeriv[range-1] = 0.;
+
   for (Int_t k=range-2;k>=0;k--)
     fHiGainSecondDeriv[k] = (fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k])/6.;
@@ -338 +506 @@
 
   //
-  // Test the possibility that the absolute maximum has not been found before the 
-  // maxpos and search from maxpos to maxpos+1 in steps of 0.2
+  // Search for the absolute maximum from maxpos to maxpos+1 in steps of 0.2
   //
   if (fAbMaxPos > upper-0.1)
@@ -379 +546 @@
   // Try a better precision. 
   //
-  const Float_t up = fAbMaxPos+step-0.055;
-  const Float_t lo = fAbMaxPos-step+0.055;
+  const Float_t up = fAbMaxPos+step-0.035;
+  const Float_t lo = fAbMaxPos-step+0.035;
   const Float_t maxpossave = fAbMaxPos;
   
@@ -387 +554 @@
   b     = x - lower;
  
-  step  = 0.02; // step size of 42 ps 
+  step  = 0.025; // step size of 83 ps 
  
   while (x<up)
@@ -410 +577 @@
 
   //
-  // Second, try from time down to time-0.2 in steps of 0.04.
+  // Second, try from time down to time-0.2 in steps of 0.025.
   //
   x     = maxpossave;
@@ -416 +583 @@
   //
   // Test the possibility that the absolute maximum has not been found between
-  // maxpos and maxpos+0.02, then we have to look between maxpos-0.02 and maxpos
+  // maxpos and maxpos+0.025, then we have to look between maxpos-0.025 and maxpos
   // which requires new setting of klocont and khicont
   //
@@ -450 +617 @@
     }
 
-  if (IsTimeType(kMaximum))
+  if (IsExtractionType(kMaximum))
     {
       time = (Float_t)fHiGainFirst + fAbMaxPos;
-      dtime = 0.02;
+      dtime = 0.025;
     }
   else
@@ -482 +649 @@
       Bool_t back = kFALSE;
       
-      Int_t maxcnt = 1000;
+      Int_t maxcnt = 50;
       Int_t cnt    = 0;
 
@@ -524 +691 @@
     }
   
-  if (IsChargeType(kAmplitude))
+  if (IsExtractionType(kAmplitude))
     {
       sum = fAbMax;
@@ -530 +697 @@
     }
 
-  if (IsChargeType(kIntegral))
+  if (IsExtractionType(kIntegral))
     {
       //
@@ -544 +711 @@
         }
 
+      if (lastslice > range)
+        lastslice = range;
+
       Int_t i = startslice;
       sum = 0.5*fHiGainSignal[i];
       
+      //
+      // We sum 1.5 times the second deriv. coefficients because these had been 
+      // divided by 6. already. Usually, 0.25*fHiGainSecondDeriv should be added.
+      //
       for (i=startslice+1; i<lastslice; i++)
         sum += fHiGainSignal[i] + 1.5*fHiGainSecondDeriv[i];
@@ -611 +785 @@
   if (sat || maxpos < 2)
     {
-      time =  IsTimeType(kMaximum) 
+      time =  IsExtractionType(kMaximum) 
         ? (Float_t)(fLoGainFirst + maxpos) 
         : (Float_t)(fLoGainFirst + maxpos - 1);
@@ -786 +960 @@
     }
 
-  if (IsTimeType(kMaximum))
+  if (IsExtractionType(kMaximum))
     {
       time = (Float_t)fLoGainFirst + fAbMaxPos;
@@ -860 +1034 @@
     }
   
-  if (IsChargeType(kAmplitude))
+  if (IsExtractionType(kAmplitude))
     {
       sum = fAbMax;
@@ -866 +1040 @@
     }
 
-  if (IsChargeType(kIntegral))
+  if (IsExtractionType(kIntegral))
     {
       //

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h

@@ -35 +35 @@
   Float_t fHalfMax;                      // Current half maximum of the spline
 
-  Byte_t  fTimeFlags;                    // Flag to hold the time extraction type
-  Byte_t  fChargeFlags;                  // Flag to hold the charge extraction type
+  Byte_t  fFlags;                        // Bit-field to hold the time extraction types
   
   Int_t  PreProcess(MParList *pList);
@@ -50 +49 @@
 public:
 
-  enum TimeType_t   { kMaximum,   kHalfMaximum };    //! Possible time   extraction types
-  enum ChargeType_t { kAmplitude, kIntegral    };    //! Possible charge extraction types
+  enum ExtractionType_t { kMaximum,   kHalfMaximum,
+                          kAmplitude, kIntegral    };    //! Possible time and charge extraction types
 
 private:
 
-  Bool_t IsTimeType  ( TimeType_t   typ )  { return TESTBIT(fTimeFlags  , typ); }
-  Bool_t IsChargeType( ChargeType_t typ )  { return TESTBIT(fChargeFlags, typ); }  
+  Bool_t IsExtractionType ( ExtractionType_t typ )  { return TESTBIT(fFlags, typ); }
   
 public:
@@ -63 +61 @@
   ~MExtractTimeAndChargeSpline();  
 
+  void SetRange    ( Byte_t hifirst=0, Byte_t hilast=0, Byte_t lofirst=0, Byte_t lolast=0 );  
   void SetResolution   ( Float_t f=fgResolution  )  { fResolution  = f;  }
   void SetRiseTime     ( Float_t f=fgRiseTime    )  { fRiseTime    = f;  }
   void SetFallTime     ( Float_t f=fgFallTime    )  { fFallTime    = f;  }
 
-  void SetTimeType     ( TimeType_t typ=kMaximum    ) { fTimeFlags = 0;   SETBIT(fTimeFlags,typ);  }
-  void SetChargeType   ( ChargeType_t typ=kAmplitude) { fChargeFlags = 0; SETBIT(fChargeFlags,typ);}
+  void SetTimeType     ( ExtractionType_t typ=kMaximum ); 
+  void SetChargeType   ( ExtractionType_t typ=kAmplitude);
   
   ClassDef(MExtractTimeAndChargeSpline, 0)   // Task to Extract Arrival Times and Charges using a Fast Cubic Spline