Changeset 4262

Timestamp:

06/02/04 02:34:43 (21 years ago)

Author:

gaug

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• msignal/MExtractTimeFastSpline.cc (modified) (14 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -52 +52 @@
        documentation
 
+   * msignal/MExtractTimeFastSpline.cc
+     - fixed some smaller bugs affecting a small part of the signals
 
 

trunk/MagicSoft/Mars/msignal/MExtractTimeFastSpline.cc

@@ -203 +203 @@
 
   for (Int_t k=range-2;k>=0;k--)
-    fHiGainSecondDeriv[k] = fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k];
+    fHiGainSecondDeriv[k] = (fHiGainSecondDeriv[k]*fHiGainSecondDeriv[k+1] + fHiGainFirstDeriv[k])/6.;
   
   //
@@ -211 +211 @@
   Float_t lower = (Float_t)maxpos-1.;
   Float_t upper = (Float_t)maxpos;
-  Float_t abmax = 0.;
+  Float_t x     = lower;
+  Float_t y     = 0.;
+  Float_t a     = 1.;
+  Float_t b     = 0.;
+  Int_t   klo = maxpos-1;
+  Int_t   khi = maxpos;
+  Float_t klocont = (Float_t)*(first+klo);
+  Float_t khicont = (Float_t)*(first+khi);
+  time          = lower;
+  Float_t abmax = klocont;
+
+  //
+  // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to 
+  // interval maxpos+1.
+  //
+  while (x<upper-0.1)
+    {
+
+      x += step;
+      a -= step;
+      b += step;
+
+      y = a*klocont
+        + b*khicont
+        + (a*a*a-a)*fHiGainSecondDeriv[klo] 
+        + (b*b*b-b)*fHiGainSecondDeriv[khi];
+
+      if (y > abmax)
+        {
+          abmax  = y;
+          time   = x;
+        }
+
+    }
+
+ if (time < lower+0.1)
+    {
+
+      upper = (Float_t)maxpos-1.;
+      lower = (Float_t)maxpos-2.;
+      x     = lower;
+      a     = 1.;
+      b     = 0.;
+      khi   = maxpos-1;
+      klo   = maxpos-2;
+      klocont = (Float_t)*(first+klo);
+      khicont = (Float_t)*(first+khi);
+
+      while (x<upper-0.1)
+        {
+
+          x += step;
+          a -= step;
+          b += step;
+          
+          y = a* klocont
+            + b* khicont
+            + (a*a*a-a)*fHiGainSecondDeriv[klo] 
+            + (b*b*b-b)*fHiGainSecondDeriv[khi];
+          
+          if (y > abmax)
+            {
+              abmax  = y;
+              time   = x;
+            }
+        }
+  }
+ 
+  const Float_t up = time+step-0.055;
+  const Float_t lo = time-step+0.055;
+  const Float_t maxpossave = time;
+
+  x     = time;
+  a     = upper - x;
+  b     = x - lower;
+
+  step  = 0.04; // step size of 83 ps 
+
+  while (x<up)
+    {
+
+      x += step;
+      a -= step;
+      b += step;
+      
+      y = a* klocont
+        + b* khicont
+        + (a*a*a-a)*fHiGainSecondDeriv[klo] 
+        + (b*b*b-b)*fHiGainSecondDeriv[khi];
+      
+      if (y > abmax)
+        {
+          abmax  = y;
+          time   = x;
+        }
+      
+    }
+
+  x     = maxpossave;
+  a     = upper - x;
+  b     = x - lower;
+
+  while (x>lo)
+    {
+
+      x -= step;
+      a += step;
+      b -= step;
+      
+      y = a* klocont
+        + b* khicont
+        + (a*a*a-a)*fHiGainSecondDeriv[klo] 
+        + (b*b*b-b)*fHiGainSecondDeriv[khi];
+      
+      if (y > abmax)
+        {
+          abmax  = y;
+          time   = x;
+        }
+      
+    }
+
+  const Float_t pedes   = ped.GetPedestal();
+  const Float_t halfmax = pedes + (abmax - pedes)/2.;
+
+  //
+  // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
+  // First, find the right FADC slice:
+  // 
+  klo   = maxpos - 1;
+  while (klo > maxpos-4)
+    {
+      if (*(first+klo) < (Byte_t)halfmax)
+        break;
+      klo--;
+    }
+
+  //
+  // Loop from the beginning of the slice upwards to reach the halfmax:
+  // With means of bisection:
+  // 
+  x     = (Float_t)klo;
+  a     = 1.;
+  b     = 0.;
+  klocont = (Float_t)*(first+klo);
+  khicont = (Float_t)*(first+klo+1);
+  time  = x;
+
+  step = 0.5;
+  Bool_t back = kFALSE;
+
+  while (step > fResolution)
+    {
+      
+      if (back)
+        {
+          x -= step;
+          a += step;
+          b -= step;
+        }
+      else
+        {
+          x += step;
+          a -= step;
+          b += step;
+        }
+      
+      y = a*klocont
+        + b*khicont
+        + (a*a*a-a)*fHiGainSecondDeriv[klo] 
+        + (b*b*b-b)*fHiGainSecondDeriv[khi];
+
+      if (y >= halfmax)
+        back = kTRUE;
+      else
+        back = kFALSE;
+
+      step /= 2.;
+      
+    }
+
+  time  = (Float_t)fHiGainFirst + x;
+  dtime = fResolution;
+}
+
+
+// --------------------------------------------------------------------------
+//
+// Calculates the arrival time for each pixel 
+//
+void MExtractTimeFastSpline::FindTimeLoGain(Byte_t *first, Float_t &time, Float_t &dtime, 
+                                        Byte_t &sat, const MPedestalPix &ped) const
+{
+  
+  const Int_t range = fLoGainLast - fLoGainFirst + 1;
+  const Byte_t *end = first + range;
+  Byte_t *p = first;
+  Byte_t max    = 0;
+  Byte_t maxpos = 0;
+
+  //
+  // Check for saturation in all other slices
+  //
+  while (++p<end)
+    {
+      if (*p > max)
+        {
+          max    = *p;
+          maxpos =  p-first;
+        }
+
+      if (*p >= fSaturationLimit)
+        {
+          sat++;
+          break;
+        }
+    }
+  
+  if (sat)
+    return;
+
+  if (maxpos < 2)
+    return;
+  
+  Float_t pp;
+
+  p = first;
+  fLoGainSecondDeriv[0] = 0.;
+  fLoGainFirstDeriv[0]  = 0.;
+
+  for (Int_t i=1;i<range-1;i++)
+    {
+      pp = fLoGainSecondDeriv[i-1] + 4.;
+      fLoGainSecondDeriv[i] = -1.0/pp;
+      fLoGainFirstDeriv [i] = *(p+2) - 2.* *(p+1) + *(p);
+      fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
+      p++;
+    }
+
+  fLoGainSecondDeriv[range-1] = 0.;
+
+  for (Int_t k=range-2;k>=0;k--)
+    fLoGainSecondDeriv[k] = (fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k])/6.;
+  
+  //
+  // Now find the maximum  
+  //
+  Float_t step  = 0.2; // start with step size of 1ns and loop again with the smaller one
+  Float_t lower = (Float_t)maxpos-1.;
+  Float_t upper = (Float_t)maxpos;
   Float_t x     = lower;
   Float_t y     = 0.;
@@ -221 +470 @@
   Float_t khicont = (Float_t)*(first+khi);
   time  = lower;
+  Float_t abmax = klocont;
 
   //
@@ -226 +476 @@
   // interval maxpos+1.
   //
-  while (x<upper+0.1)
+  while (x<upper-0.1)
     {
 
@@ -235 +485 @@
       y = a*klocont
         + b*khicont
-        + (  (a*a*a-a)*fHiGainSecondDeriv[klo] 
-             +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0;
+        + (a*a*a-a)*fLoGainSecondDeriv[klo] 
+        + (b*b*b-b)*fLoGainSecondDeriv[khi];
 
       if (y > abmax)
@@ -246 +496 @@
     }
 
- if (time > upper-0.1)
-    {
-
-      upper = (Float_t)maxpos+1.;
-      lower = (Float_t)maxpos;
+ if (time < lower+0.1)
+    {
+
+      upper = (Float_t)maxpos-1.;
+      lower = (Float_t)maxpos-2.;
       x     = lower;
       a     = 1.;
       b     = 0.;
-      khi   = maxpos+1;
-      klo   = maxpos;
+      khi   = maxpos-1;
+      klo   = maxpos-2;
       klocont = (Float_t)*(first+klo);
       khicont = (Float_t)*(first+khi);
 
-      while (x<upper+0.1)
+      while (x<upper-0.1)
         {
 
@@ -268 +518 @@
           y = a* klocont
             + b* khicont
-            + (  (a*a*a-a)*fHiGainSecondDeriv[klo] 
-                 +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0;
+            + (a*a*a-a)*fLoGainSecondDeriv[klo] 
+            + (b*b*b-b)*fLoGainSecondDeriv[khi];
           
           if (y > abmax)
@@ -279 +529 @@
   }
  
-  const Float_t up = time+step+0.015;
+  const Float_t up = time+step-0.055;
+  const Float_t lo = time-step+0.055;
+  const Float_t maxpossave = time;
 
   x     = time;
@@ -285 +537 @@
   b     = x - lower;
 
-  step  = 0.04; // step size of 83 ps 
+  step  = 0.025; // step size of 165 ps 
 
   while (x<up)
@@ -296 +548 @@
       y = a* klocont
         + b* khicont
-        + (  (a*a*a-a)*fHiGainSecondDeriv[klo] 
-             +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0;
+        + (a*a*a-a)*fLoGainSecondDeriv[klo] 
+        + (b*b*b-b)*fLoGainSecondDeriv[khi];
       
       if (y > abmax)
@@ -307 +559 @@
     }
 
-  Float_t lo = time-0.225;
-
-  x     = time;
+  x     = maxpossave;
   a     = upper - x;
   b     = x - lower;
-  khi--;
-  klo--;
-  klocont = (Float_t)*(first+klo);
-  khicont = (Float_t)*(first+khi);
 
   while (x>lo)
@@ -326 +572 @@
       y = a* klocont
         + b* khicont
-        + (  (a*a*a-a)*fHiGainSecondDeriv[klo] 
-             +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0;
+        + (a*a*a-a)*fLoGainSecondDeriv[klo] 
+        + (b*b*b-b)*fLoGainSecondDeriv[khi];
       
       if (y > abmax)
@@ -384 +630 @@
       y = a*klocont
         + b*khicont
-        + (  (a*a*a-a)*fHiGainSecondDeriv[klo] 
-             +(b*b*b-b)*fHiGainSecondDeriv[khi] )/6.0;
+        + (a*a*a-a)*fLoGainSecondDeriv[klo] 
+        + (b*b*b-b)*fLoGainSecondDeriv[khi];
 
       if (y >= halfmax)
@@ -393 +639 @@
 
       step /= 2.;
-    }
-
-  time  = (Float_t)fHiGainFirst + x;
-  dtime = fResolution;
-}
-
-
-// --------------------------------------------------------------------------
-//
-// Calculates the arrival time for each pixel 
-//
-void MExtractTimeFastSpline::FindTimeLoGain(Byte_t *first, Float_t &time, Float_t &dtime, 
-                                        Byte_t &sat, const MPedestalPix &ped) const
-{
-  
-  const Int_t range = fLoGainLast - fLoGainFirst + 1;
-  const Byte_t *end = first + range;
-  Byte_t *p = first;
-  Byte_t max    = 0;
-  Byte_t maxpos = 0;
-
-  //
-  // Check for saturation in all other slices
-  //
-  while (++p<end)
-    {
-      if (*p > max)
-        {
-          max    = *p;
-          maxpos =  p-first;
-        }
-
-      if (*p >= fSaturationLimit)
-        {
-          sat++;
-          break;
-        }
-    }
-  
-  if (sat)
-    return;
-
-  if (maxpos < 2)
-    return;
-  
-  Float_t pp;
-
-  p = first;
-  fLoGainSecondDeriv[0] = 0.;
-  fLoGainFirstDeriv[0]  = 0.;
-
-  for (Int_t i=1;i<range-1;i++)
-    {
-      pp = fLoGainSecondDeriv[i-1] + 4.;
-      fLoGainSecondDeriv[i] = -1.0/pp;
-      fLoGainFirstDeriv [i] = *(p+2) - 2.* *(p+1) + *(p);
-      fLoGainFirstDeriv [i] = (6.0*fLoGainFirstDeriv[i]-fLoGainFirstDeriv[i-1])/pp;
-      p++;
-    }
-
-  fLoGainSecondDeriv[range-1] = 0.;
-
-  for (Int_t k=range-2;k>=0;k--)
-    fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k];
-  
-  //
-  // Now find the maximum  
-  //
-  Float_t step  = 0.2; // start with step size of 1ns and loop again with the smaller one
-  Float_t lower = (Float_t)maxpos-1.;
-  Float_t upper = (Float_t)maxpos;
-  Float_t abmax = 0.;
-  Float_t x     = lower;
-  Float_t y     = 0.;
-  Float_t a     = 1.;
-  Float_t b     = 0.;
-  Int_t   klo = maxpos-1;
-  Int_t   khi = maxpos;
-  Float_t klocont = (Float_t)*(first+klo);
-  Float_t khicont = (Float_t)*(first+khi);
-  time  = lower;
-
-  //
-  // Search for the maximum, starting in interval maxpos-1. If no maximum is found, go to 
-  // interval maxpos+1.
-  //
-  while (x<upper+0.1)
-    {
-
-      x += step;
-      a -= step;
-      b += step;
-
-      y = a*klocont
-        + b*khicont
-        + (  (a*a*a-a)*fLoGainSecondDeriv[klo] 
-             +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0;
-
-      if (y > abmax)
-        {
-          abmax  = y;
-          time   = x;
-        }
-
-    }
-
- if (time > upper-0.1)
-    {
-
-      upper = (Float_t)maxpos+1.;
-      lower = (Float_t)maxpos;
-      x     = lower;
-      a     = 1.;
-      b     = 0.;
-      khi   = maxpos+1;
-      klo   = maxpos;
-      klocont = (Float_t)*(first+klo);
-      khicont = (Float_t)*(first+khi);
-
-      while (x<upper+0.1)
-        {
-
-          x += step;
-          a -= step;
-          b += step;
-          
-          y = a* klocont
-            + b* khicont
-            + (  (a*a*a-a)*fLoGainSecondDeriv[klo] 
-                 +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0;
-          
-          if (y > abmax)
-            {
-              abmax  = y;
-              time   = x;
-            }
-        }
-  }
- 
-  const Float_t up = time+step+0.015;
-
-  x     = time;
-  a     = upper - x;
-  b     = x - lower;
-
-  step  = 0.025; // step size of 165 ps 
-
-  while (x<up)
-    {
-
-      x += step;
-      a -= step;
-      b += step;
-      
-      y = a* klocont
-        + b* khicont
-        + (  (a*a*a-a)*fLoGainSecondDeriv[klo] 
-             +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0;
-      
-      if (y > abmax)
-        {
-          abmax  = y;
-          time   = x;
-        }
-      
-    }
-
-  Float_t lo = time-0.225;
-
-  x     = time;
-  a     = upper - x;
-  b     = x - lower;
-  khi--;
-  klo--;
-  klocont = (Float_t)*(first+klo);
-  khicont = (Float_t)*(first+khi);
-
-  while (x>lo)
-    {
-
-      x -= step;
-      a += step;
-      b -= step;
-      
-      y = a* klocont
-        + b* khicont
-        + (  (a*a*a-a)*fLoGainSecondDeriv[klo] 
-             +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0;
-      
-      if (y > abmax)
-        {
-          abmax  = y;
-          time   = x;
-        }
-      
-    }
-
-  const Float_t pedes   = ped.GetPedestal();
-  const Float_t halfmax = pedes + (abmax - pedes)/2.;
-
-  //
-  // Now, loop from the maximum bin leftward down in order to find the position of the half maximum.
-  // First, find the right FADC slice:
-  // 
-  klo   = maxpos;
-  while (klo > maxpos-4)
-    {
-      if (*(first+klo) < (Byte_t)halfmax)
-        break;
-      klo--;
-    }
-
-  //
-  // Loop from the beginning of the slice upwards to reach the halfmax:
-  // With means of bisection:
-  // 
-  x     = (Float_t)klo;
-  a     = 1.;
-  b     = 0.;
-  klocont = (Float_t)*(first+klo);
-  khicont = (Float_t)*(first+klo+1);
-  time  = x;
-
-  step = 0.5;
-  Bool_t back = kFALSE;
-
-  while (step > fResolution)
-    {
-      
-      if (back)
-        {
-          x -= step;
-          a += step;
-          b -= step;
-        }
-      else
-        {
-          x += step;
-          a -= step;
-          b += step;
-        }
-      
-      y = a*klocont
-        + b*khicont
-        + (  (a*a*a-a)*fLoGainSecondDeriv[klo] 
-             +(b*b*b-b)*fLoGainSecondDeriv[khi] )/6.0;
-
-      if (y >= halfmax)
-        back = kTRUE;
-      else
-        back = kFALSE;
-
-      step /= 2.;
+
     }