Changeset 7055

Timestamp:

05/18/05 17:54:27 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

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

trunk/MagicSoft/Mars/Changelog

@@ -41 +41 @@
    * mjobs/MJStar.cc:
      - moved processing of CC-branch to the beginning of the tasklist
+
+   * msignal/MExtractTimeAndChargeSpline.[h,cc]:
+     - introduced some small changes to the validity range of
+       some variables
+     - determin the higher bound above which no search is done
+       analog to the lower bound using the fall-time
+     - CalcIntegral[Hi,Lo]Gain now returns sum. No need for a reference
+     - fixed calling Integral[HI,Lo]Gain in cases we are at the edge of
+       the valid range -- at a lot of position in the code random memory 
+       above or below the arrays have been accessed.
+     - improved the numercila stability of CalcIntegral[Hi,Lo]Gain
+       more by calculating the number of steps from the rise and fall time.
+       this should at least give consistent results on the same machine!
 
 

trunk/MagicSoft/Mars/NEWS

@@ -77 +77 @@
        adding the step-size (which results in numerical uncertanties
        exspecially if multiplied with large numbers)
+     + A lot of fixes have been introduced which effects integrating the
+       spline at the edges of the valid range. In this case any memory
+       was randomly accessed.
+     ! No result obtained with the Spline before can be trusted! Due to
+       random memory access it might by completely random!
 
    - callisto: set new defaults in MExtractTimeAndChargeDigitalFilter:

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc

@@ -432 +432 @@
         }
     }
-  
+
   fAbMax = fHiGainSignal[maxpos];
-
-  Float_t pp;
 
   fHiGainSecondDeriv[0] = 0.;
@@ -442 +440 @@
   for (Int_t i=1;i<range-1;i++)
     {
-      pp = fHiGainSecondDeriv[i-1] + 4.;
+      const Float_t pp = fHiGainSecondDeriv[i-1] + 4.;
       fHiGainSecondDeriv[i] = -1.0/pp;
-      fHiGainFirstDeriv [i] = fHiGainSignal[i+1] - fHiGainSignal[i] - fHiGainSignal[i] + fHiGainSignal[i-1];
+      fHiGainFirstDeriv [i] = 2*(fHiGainSignal[i+1] - fHiGainSignal[i]);
       fHiGainFirstDeriv [i] = (6.0*fHiGainFirstDeriv[i]-fHiGainFirstDeriv[i-1])/pp;
     }
@@ -474 +472 @@
         }
       else
-        {
-          Float_t start = 2. + nsx;
-          Float_t last  = start + fRiseTimeHiGain + fFallTimeHiGain;
-      
-          if (int(last) > range)
-            {
-              const Int_t diff = range - int(last);
-              last  -= diff;
-              start -= diff;
-            }
-          
-          CalcIntegralHiGain(sum, start, last);
-        }
+          sum = CalcIntegralHiGain(2. + nsx, range);
+
       fRandomIter++;
       return;
@@ -492 +479 @@
 
   //
-  // Allow no saturated slice 
-  // and 
+  // Allow no saturated slice and
   // Don't start if the maxpos is too close to the limits.
   //
-  if (sat || maxpos < TMath::Ceil(fRiseTimeHiGain) || maxpos > range-2)
-    {
-
+  const Bool_t limlo = maxpos <       TMath::Ceil(fRiseTimeHiGain);
+  const Bool_t limup = maxpos > range-TMath::Ceil(fFallTimeHiGain)-1;
+  if (sat || limlo || limup)
+    {
       dtime = 1.0;
       if (fExtractionType == kAmplitude)
@@ -506 +493 @@
           return;
         }
-      
-      if (maxpos > range - 2)
-        CalcIntegralHiGain(sum, (Float_t)range - fRiseTimeHiGain - fFallTimeHiGain, (Float_t)range - 0.001);
-      else
-        CalcIntegralHiGain(sum, 0.001, fRiseTimeHiGain + fFallTimeHiGain);
-
-      time =  (Float_t)(fHiGainFirst + maxpos - 1);
+
+      sum  = CalcIntegralHiGain(limlo ? 0 : range, range);
+      time = (Float_t)(fHiGainFirst + maxpos - 1);
       return;
     }
@@ -561 +544 @@
   if (fAbMaxPos > upper-0.1)
     {
-
       upper   = 1. + maxpos;
       lower   = (Float_t)maxpos;
@@ -678 +660 @@
   while (klo > 0)
     {
-      klo--;
-      if (fHiGainSignal[klo] < fHalfMax)
+      if (fHiGainSignal[--klo] < fHalfMax)
         break;
     }
@@ -719 +700 @@
         + (b*b*b-b)*fHiGainSecondDeriv[khi];
       
-      if (y > fHalfMax)
-        back = kTRUE;
-      else
-        back = kFALSE;
+      back = y > fHalfMax;
       
       if (++cnt > maxcnt)
@@ -730 +708 @@
     }
   
-  time  = (Float_t)fHiGainFirst + x;
   //
   // Now integrate the whole thing!
   // 
-  
-  Float_t start = fAbMaxPos - fRiseTimeHiGain;
-  Float_t last  = fAbMaxPos + fFallTimeHiGain;
-  
-  const Int_t diff = int(last) - range;
-  
-  if (diff > 0)
-    {
-      last  -= diff;
-      start -= diff;
-    }
-  
-  CalcIntegralHiGain(sum, start, last);
+  time = (Float_t)fHiGainFirst + x;
+  sum  = CalcIntegralHiGain(fAbMaxPos - fRiseTimeHiGain, range);
 }
 
@@ -795 +761 @@
   fAbMax = fLoGainSignal[maxpos];
   
-  Float_t pp;
-
   fLoGainSecondDeriv[0] = 0.;
   fLoGainFirstDeriv[0]  = 0.;
@@ -802 +766 @@
   for (Int_t i=1;i<range-1;i++)
     {
-      pp = fLoGainSecondDeriv[i-1] + 4.;
+      const Float_t pp = fLoGainSecondDeriv[i-1] + 4.;
       fLoGainSecondDeriv[i] = -1.0/pp;
-      fLoGainFirstDeriv [i] = fLoGainSignal[i+1] - fLoGainSignal[i] - fLoGainSignal[i] + fLoGainSignal[i-1];
+      fLoGainFirstDeriv [i] = 2*(fLoGainSignal[i+1] - fLoGainSignal[i]);
       fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
     }
@@ -833 +797 @@
         }
       else
-        {
-          Float_t start = 2. + nsx;
-          Float_t last  = start + fRiseTimeLoGain + fFallTimeLoGain;
-      
-          if (int(last) > range)
-            {
-              const Int_t diff = range - int(last);
-              last  -= diff;
-              start -= diff;
-            }
-          
-          CalcIntegralLoGain(sum, start, last);
-        }
+          sum = CalcIntegralLoGain(2. + nsx, range);
+
       fRandomIter++;
       return;
     }
   //
-  // Allow no saturated slice 
-  // and 
+  // Allow no saturated slice and
   // Don't start if the maxpos is too close to the limits.
   //
-  if (sat || maxpos < TMath::Ceil(fRiseTimeLoGain) || maxpos > range-2)
-    {
-
+  const Bool_t limlo = maxpos <       TMath::Ceil(fRiseTimeLoGain);
+  const Bool_t limup = maxpos > range-TMath::Ceil(fFallTimeLoGain)-1;
+  if (sat || limlo || limup)
+    {
       dtime = 1.0;
       if (fExtractionType == kAmplitude)
@@ -864 +817 @@
           return;
         }
-      
-      if (maxpos > range-2)
-        CalcIntegralLoGain(sum, (Float_t)range - fRiseTimeLoGain - fFallTimeLoGain -1., (Float_t)range - 0.001);
-      else
-        CalcIntegralLoGain(sum, 0.001, fRiseTimeLoGain + fFallTimeLoGain + 1.);
-
+
+      sum  = CalcIntegralLoGain(limlo ? 0 : range, range);
       time = (Float_t)(fLoGainFirst + maxpos - 1);
       return;
@@ -1080 +1029 @@
         + (b*b*b-b)*fLoGainSecondDeriv[khi];
       
-      if (y > fHalfMax)
-        back = kTRUE;
-      else
-        back = kFALSE;
+      back = y > fHalfMax;
       
       if (++cnt > maxcnt)
@@ -1091 +1037 @@
     }
   
-  time  = x + (Int_t)fLoGainFirst;
-
   //
   // Now integrate the whole thing!
   // 
-  Float_t start = fAbMaxPos - fRiseTimeLoGain;
-  Float_t last  = fAbMaxPos + fFallTimeLoGain;
-  
-  const Int_t diff = int(last) - range;
-  
-  if (diff > 0)
-    {
-      last  -= diff;
-      start -= diff;
-    }
-  CalcIntegralLoGain(sum, start, last);
+  time = x + (Int_t)fLoGainFirst;
+  sum  = CalcIntegralLoGain(fAbMaxPos - fRiseTimeLoGain, range);
 }
 
-void MExtractTimeAndChargeSpline::CalcIntegralHiGain(Float_t &sum, Float_t start, Float_t last)
+Float_t MExtractTimeAndChargeSpline::CalcIntegralHiGain(Float_t start, Float_t range) const
 {
-    const Float_t step = 0.2;
-
+    // The number of steps is calculated directly from the integration
+    // window. This is the only way to ensure we are not dealing with
+    // numerical rounding uncertanties, because we always get the same
+    // value under the same conditions -- it might still be different on
+    // other machines!
+    const Float_t step  = 0.2;
+    const Float_t width = fRiseTimeHiGain+fFallTimeHiGain;
+    const Float_t max   = range-1 - (width+step);
+    const Int_t   num   = TMath::Nint(width/step);
+
+    // The order is important. In some cases (limlo-/limup-check) it can
+    // happen than max<0. In this case we start at 0
+    if (start > max)
+        start = max;
     if (start < 0)
-    {
-        last -= start;
-        start = 0.;
-    }
-
-    const Int_t n = TMath::Nint((last-start)/step);
+        start = 0;
 
     start += step/2;
 
-    sum = 0.;
-    for (Int_t i=0; i<n; i++)
+    Double_t sum = 0.;
+    for (Int_t i=0; i<num; i++)
     {
         const Float_t x = start+i*step;
@@ -1149 +1091 @@
     }
     sum *= step; // Transform sum in integral
+    return sum;
 }
 
-void MExtractTimeAndChargeSpline::CalcIntegralLoGain(Float_t &sum, Float_t start, Float_t last)
+Float_t MExtractTimeAndChargeSpline::CalcIntegralLoGain(Float_t start, Float_t range) const
 {
-    const Float_t step = 0.2;
-
+    // The number of steps is calculated directly from the integration
+    // window. This is the only way to ensure we are not dealing with
+    // numerical rounding uncertanties, because we always get the same
+    // value under the same conditions -- it might still be different on
+    // other machines!
+    const Int_t   cnt   = 5;
+    const Float_t width = fRiseTimeLoGain+fFallTimeLoGain;
+    const Float_t max   = range-1 - (width+step);
+    const Int_t   num   = TMath::Nint(width/step);
+
+    // The order is important. In some cases (limlo-/limup-check) it can
+    // happen than max<0. In this case we start at 0
+    if (start > max)
+        start = max;
     if (start < 0)
-    {
-        last -= start;
-        start = 0.;
-    }
-
-    const Int_t n = TMath::Nint((last-start)/step);
+        start = 0;
 
     start += step/2;
 
-    sum = 0.;
-    for (Int_t i=0; i<n; i++)
+    Double_t sum = 0.;
+    for (Int_t i=0; i<num; i++)
     {
         const Float_t x = start+i*step;
@@ -1191 +1141 @@
     }
     sum *= step; // Transform sum in integral
+    return sum;
 }
 

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h

@@ -6 +6 @@
 #endif
 
-#ifndef MARS_MArrayF
-#include "MArrayF.h"
+#ifndef ROOT_TArrayF
+#include "TArrayF.h"
 #endif
 
@@ -27 +27 @@
   static const Float_t fgLoGainStartShift; //! Default for fLoGainStartShift (now set to -1.6)
   
-  MArrayF fHiGainSignal;                //! Need fast access to the signals in a float way
-  MArrayF fLoGainSignal;                //! Store them in separate arrays
-  MArrayF fHiGainFirstDeriv;            //! High-gain discretized first derivatives
-  MArrayF fLoGainFirstDeriv;            //! Low-gain discretized first derivatives
-  MArrayF fHiGainSecondDeriv;           //! High-gain discretized second derivatives
-  MArrayF fLoGainSecondDeriv;           //! Low-gain discretized second derivatives
+  TArrayF fHiGainSignal;                //! Need fast access to the signals in a float way
+  TArrayF fLoGainSignal;                //! Store them in separate arrays
+  TArrayF fHiGainFirstDeriv;            //! High-gain discretized first derivatives
+  TArrayF fLoGainFirstDeriv;            //! Low-gain discretized first derivatives
+  TArrayF fHiGainSecondDeriv;           //! High-gain discretized second derivatives
+  TArrayF fLoGainSecondDeriv;           //! Low-gain discretized second derivatives
 
   Float_t fAbMax;                       //! Current maximum of the spline
@@ -53 +53 @@
   Bool_t  InitArrays();
   
-  void CalcIntegralHiGain(Float_t &sum, Float_t start, Float_t last);
-  void CalcIntegralLoGain(Float_t &sum, Float_t start, Float_t last);
+  Float_t CalcIntegralHiGain(Float_t start, Float_t range) const;
+  Float_t CalcIntegralLoGain(Float_t start, Float_t range) const;
 
 public: