Changeset 5747 for trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc

Timestamp:

01/08/05 20:17:34 (20 years ago)

Author:

gaug

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc (modified) (23 diffs)

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc

@@ -66 +66 @@
 //   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:
@@ -102 +101 @@
 //     (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)
-//                 or: TMath::Floor(fRiseTime + fFallTime + 1.) in the case of the low-gain
+//  2) If Charge Type: kIntegral was chosen: fRiseTime + fFallTime
+//                 or: fRiseTime + fFallTime + 1. in the case of the low-gain
 //
 //  Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast) 
@@ -127 +122 @@
 //        for the case, the extraction type kIntegral has been chosen.
 //
-//  Call: - SetTimeType(MExtractTimeAndChargeSpline::kMaximum) for extraction 
+//  Call: - SetChargeType(MExtractTimeAndChargeSpline::kAmplitude) for the 
+//          computation of the amplitude at the maximum (default) and 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. The Low Gain is computed with one more slice at the falling
-//          edge.
-//          --> needs one more simple summation loop over 7 slices.
+//          the maximum. The Low Gain is computed with half a slice more at the rising
+//          edge and half a slice more at the falling edge.
+//          The time of the half maximum is returned.
+//          --> needs one function evaluations but is more precise
 //        
 //////////////////////////////////////////////////////////////////////////////
@@ -159 +150 @@
 const Byte_t  MExtractTimeAndChargeSpline::fgLoGainFirst  = 2;
 const Byte_t  MExtractTimeAndChargeSpline::fgLoGainLast   = 14;
-const Float_t MExtractTimeAndChargeSpline::fgResolution   = 0.025;
+const Float_t MExtractTimeAndChargeSpline::fgResolution   = 0.05;
 const Float_t MExtractTimeAndChargeSpline::fgRiseTime     = 1.5;
 const Float_t MExtractTimeAndChargeSpline::fgFallTime     = 4.5;
@@ -177 +168 @@
 //
 MExtractTimeAndChargeSpline::MExtractTimeAndChargeSpline(const char *name, const char *title) 
-    : fAbMax(0.), fAbMaxPos(0.), fHalfMax(0.), fRandomIter(0)
+    : fAbMax(0.), fAbMaxPos(0.), fHalfMax(0.), fHalfMaxPos(0.), fRandomIter(0)
 {
 
@@ -189 +180 @@
   SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
 
-  SetTimeType();
   SetChargeType();
 
@@ -212 +202 @@
     SetChargeType(kAmplitude);
   
-}
-
-
-//-------------------------------------------------------------------
-//
-// 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);
-
 }
 
@@ -316 +289 @@
   const Byte_t *end = first + range;
   Byte_t       *p  = first;
-  Int_t  count     = 0;
+
+  sat = 0;
 
   const Float_t pedes  = ped.GetPedestal();
   const Float_t ABoffs = ped.GetPedestalABoffset();
 
-  Float_t pedmean[2];
-  pedmean[0] = pedes + ABoffs;
-  pedmean[1] = pedes - ABoffs;
+  const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
 
   fAbMax        = 0.;
   fAbMaxPos     = 0.;
+  fHalfMax      = 0.;
+  fHalfMaxPos   = 0.;
   Int_t  maxpos = 0;
+  Int_t  max    = 0;
 
   //
   // Check for saturation in all other slices
   //
+  Int_t ids = fHiGainFirst;
+  Float_t *sample = fHiGainSignal.GetArray();
   while (p<end)
     {
 
-      const Int_t ids      = fHiGainFirst + count ;
-      const Float_t signal = (Float_t)*p - pedmean[(ids+abflag) & 0x1];
-      fHiGainSignal[count] = signal;
-
-      if (signal > fAbMax + 0.1) /* the 0.1 is necessary for the ultra-high enery events saturating many slices */
-        {
-          fAbMax = signal;
-          maxpos = count;
+      *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
+
+      if (*p > max)
+        {
+          maxpos = ids-fHiGainFirst-1;
+          max    = *p;
         }
 
       if (*p++ >= fSaturationLimit)
-            sat++;
-      
-      count++;                                                    
+        if (!sat)
+          sat = ids-3;
+      
     }
   
@@ -359 +334 @@
         {
 
-          const Int_t ids      = fHiGainFirst + range ;
-          const Float_t signal = (Float_t)*logain - pedmean[(ids+abflag) & 0x1];
-          fHiGainSignal[range] = signal;
-          
-          if (signal > fAbMax)
+          *sample++ = (Float_t)*logain - pedmean[(ids++ + abflag) & 0x1];
+
+          if (*logain > max)
             { 
-              fAbMax = signal;
-              maxpos = range;
+              maxpos = ids-fHiGainFirst-1;
+              max    = *logain;
             }
           
-          if (*logain >= fSaturationLimit)
-            sat++;
+          if (*logain++ >= fSaturationLimit)
+            if (!sat)
+              sat = ids-3;
 
           range++;
-          logain++;
-        }
-    }
-  
+        }
+    }
+  
+  fAbMax = fHiGainSignal[maxpos];
+
   Float_t pp;
 
@@ -434 +409 @@
 
   //
-  // Allow one saturated slice 
+  // Allow no saturated slice 
   // and 
   // Don't start if the maxpos is too close to the limits.
   //
-  if (sat > 1 || maxpos < 1 || maxpos > range-2)
-    {
-      time =  IsExtractionType(kMaximum) 
-        ? (Float_t)(fHiGainFirst + maxpos) 
-        : (Float_t)(fHiGainFirst + maxpos - 1);
-
+  if (sat || maxpos < 1 || maxpos > range-2)
+    {
+      dtime = 0.5;
       if (IsExtractionType(kAmplitude))
         {
-          sum = fAbMax;
+          sum  = fAbMax;
+          time = (Float_t)(fHiGainFirst + maxpos);
           return;
         }
@@ -455 +428 @@
         CalcIntegralHiGain(sum, 0.001, fRiseTime + fFallTime);
 
+      time =  (Float_t)(fHiGainFirst + maxpos - 1);
       return;
     }
       
+  dtime = fResolution;
+
   //
   // Now find the maximum  
@@ -533 +509 @@
   // Try a better precision. 
   //
-  const Float_t up = fAbMaxPos+step - 1.5*fResolution;
-  const Float_t lo = fAbMaxPos-step + 1.5*fResolution;
+  const Float_t up = fAbMaxPos+step - 3.0*fResolution;
+  const Float_t lo = fAbMaxPos-step + 3.0*fResolution;
   const Float_t maxpossave = fAbMaxPos;
   
@@ -541 +517 @@
   b     = x - lower;
  
-  step  = fResolution; // step size of 83 ps 
+  step  = 2.*fResolution; // step size of 0.1 FADC slices
  
   while (x<up)
@@ -569 +545 @@
   //
   // Test the possibility that the absolute maximum has not been found between
-  // maxpos and maxpos+0.025, then we have to look between maxpos-0.025 and maxpos
+  // maxpos and maxpos+0.05, then we have to look between maxpos-0.05 and maxpos
   // which requires new setting of klocont and khicont
   //
-  if (x < lower + fResolution/2.)
+  if (x < lower + fResolution)
     {
       klo--;
@@ -602 +578 @@
     }
 
-  if (IsExtractionType(kMaximum))
-    {
-      time = (Float_t)fHiGainFirst + fAbMaxPos;
-      dtime = fResolution;
-    }
-  else
-    {
-      fHalfMax = fAbMax/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 >= 0)
-        {
-          klo--;
-          if (fHiGainSignal[klo] < fHalfMax)
-            break;
-        }
-      
-      khi = klo+1;
-      //
-      // Loop from the beginning of the slice upwards to reach the fHalfMax:
-      // With means of bisection:
-      // 
-      x     = (Float_t)klo;
-      a     = 1.;
-      b     = 0.;
-      
-      step = 0.5;
-      Bool_t back = kFALSE;
-      
-      Int_t maxcnt = 20;
-      Int_t cnt    = 0;
-
-      while (TMath::Abs(y-fHalfMax) > fResolution)
-        {
-          
-          if (back)
-            {
-              x -= step;
-              a += step;
-              b -= step;
-            }
-          else
-            {
-              x += step;
-              a -= step;
-              b += step;
-            }
-          
-          y = a*fHiGainSignal[klo]
-            + b*fHiGainSignal[khi]
-            + (a*a*a-a)*fHiGainSecondDeriv[klo] 
-            + (b*b*b-b)*fHiGainSecondDeriv[khi];
-          
-          if (y > fHalfMax)
-            back = kTRUE;
-          else
-            back = kFALSE;
-
-          if (++cnt > maxcnt)
-            break;
-          
-          step /= 2.;
-        }
-
-      time  = (Float_t)fHiGainFirst + x;
-      dtime = fResolution;
-    }
-  
   if (IsExtractionType(kAmplitude))
     {
-      sum = fAbMax;
+      time  = fAbMaxPos + (Int_t)fHiGainFirst;
+      sum   = fAbMax;
       return;
     }
 
-  if (IsExtractionType(kIntegral))
-    {
-      //
-      // Now integrate the whole thing!
-      // 
-
-      Float_t start = fAbMaxPos - fRiseTime;
-      Float_t last  = fAbMaxPos + fFallTime;
-      
-      const Int_t diff = int(last) - range;
-
-      if (diff > 0)
-        {
-          last  -= diff;
-          start -= diff;
-        }
-
-      CalcIntegralHiGain(sum, start, last);
-    }
+  fHalfMax = fAbMax/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 >= 0)
+    {
+      klo--;
+      if (fHiGainSignal[klo] < fHalfMax)
+        break;
+    }
+  
+  khi = klo+1;
+  //
+  // Loop from the beginning of the slice upwards to reach the fHalfMax:
+  // With means of bisection:
+  // 
+  x     = (Float_t)klo;
+  a     = 1.;
+  b     = 0.;
+  
+  step = 0.5;
+  Bool_t back = kFALSE;
+  
+  Int_t maxcnt = 20;
+  Int_t cnt    = 0;
+  
+  while (TMath::Abs(y-fHalfMax) > fResolution)
+    {
+      
+      if (back)
+        {
+          x -= step;
+          a += step;
+          b -= step;
+        }
+      else
+        {
+          x += step;
+          a -= step;
+          b += step;
+        }
+      
+      y = a*fHiGainSignal[klo]
+        + b*fHiGainSignal[khi]
+        + (a*a*a-a)*fHiGainSecondDeriv[klo] 
+        + (b*b*b-b)*fHiGainSecondDeriv[khi];
+      
+      if (y > fHalfMax)
+        back = kTRUE;
+      else
+        back = kFALSE;
+      
+      if (++cnt > maxcnt)
+        break;
+      
+      step /= 2.;
+    }
+  
+  time  = (Float_t)fHiGainFirst + x;
+  //
+  // Now integrate the whole thing!
+  // 
+  
+  Float_t start = fAbMaxPos - fRiseTime;
+  Float_t last  = fAbMaxPos + fFallTime;
+  
+  const Int_t diff = int(last) - range;
+  
+  if (diff > 0)
+    {
+      last  -= diff;
+      start -= diff;
+    }
+  
+  CalcIntegralHiGain(sum, start, last);
 }
 
@@ -718 +682 @@
   const Float_t ABoffs = ped.GetPedestalABoffset();
 
-  Float_t pedmean[2];
-  pedmean[0] = pedes + ABoffs;
-  pedmean[1] = pedes - ABoffs;
+  const Float_t pedmean[2] = { pedes + ABoffs, pedes - ABoffs };
 
   fAbMax        = 0.;
   fAbMaxPos     = 0.;
   Int_t  maxpos = 0;
-  Int_t  count  = 0;
+  Int_t  max    = 0;
 
   //
   // Check for saturation in all other slices
   //
+  Int_t    ids    = fLoGainFirst;
+  Float_t *sample = fLoGainSignal.GetArray();
   while (p<end)
     {
 
-      const Int_t ids      =  count + fLoGainFirst;
-      const Float_t signal = (Float_t)*p - pedmean[(ids+abflag) & 0x1];
-      fLoGainSignal[count] = signal;
-
-      if (signal > fAbMax + 0.1)
-        {
-          fAbMax = signal;
-          maxpos = count;
+      *sample++ = (Float_t)*p - pedmean[(ids++ + abflag) & 0x1];
+
+      if (*p > max)
+        {
+          maxpos = ids-fLoGainFirst-1;
+          max    = *p;
         }
 
       if (*p++ >= fSaturationLimit)
         sat++;
-      
-      count++;                                                    
-    }
+    }
+  
+  fAbMax = fLoGainSignal[maxpos];
   
   Float_t pp;
@@ -808 +770 @@
   // Don't start if the maxpos is too close to the limits.
   //
-  if (sat || maxpos < 2 || maxpos > range-2)
-    {
-      time =  IsExtractionType(kMaximum) 
-        ? (Float_t)(fLoGainFirst + maxpos) 
-        : (Float_t)(fLoGainFirst + maxpos - 1);
-      
+  if (sat || maxpos < TMath::Ceil(fRiseTime+0.45) || maxpos > range-2)
+    {
+      dtime = 0.5;
       if (IsExtractionType(kAmplitude))
         {
+          time = (Float_t)(fLoGainFirst + maxpos);
           sum = fAbMax;
           return;
@@ -825 +785 @@
         CalcIntegralLoGain(sum, 0.001, fRiseTime + fFallTime + 1.);
 
+      time = (Float_t)(fLoGainFirst + maxpos - 1);
       return;
     }
 
-  if (maxpos < (Int_t)(fRiseTime+2.))
-    {
-      time =  IsExtractionType(kMaximum) 
-        ? (Float_t)(fLoGainFirst + maxpos) 
-        : (Float_t)(fLoGainFirst + maxpos - 1);
-
-      if (maxpos > range-2)
-        CalcIntegralLoGain(sum, (Float_t)range - fRiseTime - fFallTime-1., (Float_t)range - 0.001);
-      else
-        CalcIntegralLoGain(sum, 0.001, fRiseTime + fFallTime + 1.);
-
-      return;
-    }
+  dtime = fResolution;
 
   //
@@ -920 +869 @@
   // Try a better precision. 
   //
-  const Float_t up = fAbMaxPos+step - 1.5*fResolution;
-  const Float_t lo = fAbMaxPos-step + 1.5*fResolution;
+  const Float_t up = fAbMaxPos+step - 3.0*fResolution;
+  const Float_t lo = fAbMaxPos-step + 3.0*fResolution;
   const Float_t maxpossave = fAbMaxPos;
   
@@ -928 +877 @@
   b     = x - lower;
  
-  step  = fResolution; // step size of fResolution (33 ps )
+  step  = 2.*fResolution; // step size of 0.1 FADC slice
  
   while (x<up)
@@ -956 +905 @@
   //
   // 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.05, then we have to look between maxpos-0.05 and maxpos
   // which requires new setting of klocont and khicont
   //
-  if (x < lower + fResolution/2.)
+  if (x < lower + fResolution)
     {
       klo--;
@@ -989 +938 @@
     }
 
-  if (IsExtractionType(kMaximum))
+  if (IsExtractionType(kAmplitude))
     {
       time = fAbMaxPos + (Int_t)fLoGainFirst;
-      dtime = fResolution;
-    }
-  else
-    {
-      fHalfMax = fAbMax/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 > 0)
-        {
-          klo--;
-          if (fLoGainSignal[klo] < fHalfMax)
-            break;
-        }
-      
-      khi = klo+1;
-      //
-      // Loop from the beginning of the slice upwards to reach the fHalfMax:
-      // With means of bisection:
-      // 
-      x     = (Float_t)klo;
-      a     = 1.;
-      b     = 0.;
-      
-      step = 0.5;
-      Bool_t back = kFALSE;
-      
-      Int_t maxcnt = 20;
-      Int_t cnt    = 0;
-
-      while (TMath::Abs(y-fHalfMax) > fResolution)
-        {
-          
-          if (back)
-            {
-              x -= step;
-              a += step;
-              b -= step;
-            }
-          else
-            {
-              x += step;
-              a -= step;
-              b += step;
-            }
-          
-          y = a*fLoGainSignal[klo]
-            + b*fLoGainSignal[khi]
-            + (a*a*a-a)*fLoGainSecondDeriv[klo] 
-            + (b*b*b-b)*fLoGainSecondDeriv[khi];
-          
-          if (y > fHalfMax)
-            back = kTRUE;
-          else
-            back = kFALSE;
-
-          if (++cnt > maxcnt)
-            break;
-          
-          step /= 2.;
-        }
-
-      time  = x + (Int_t)fLoGainFirst;
-      dtime = fResolution;
-    }
-  
-  if (IsExtractionType(kAmplitude))
-    {
-      sum = fAbMax;
+      sum  = fAbMax;
       return;
     }
-
-  if (IsExtractionType(kIntegral))
-    {
-      //
-      // Now integrate the whole thing!
-      // 
-      Float_t start = fAbMaxPos - fRiseTime - 0.5;
-      Float_t last  = fAbMaxPos + fFallTime + 0.5;
-      
-      const Int_t diff = int(last) - range;
-
-      if (diff > 0)
-        {
-          last  -= diff;
-          start -= diff;
-        }
-      CalcIntegralLoGain(sum, start, last);
-      //      *fLog << inf << time << "  " << sum << "  " << start << "  " << last << endl;
-    }
-
+  
+  fHalfMax = fAbMax/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 > 0)
+    {
+      klo--;
+      if (fLoGainSignal[klo] < fHalfMax)
+        break;
+    }
+  
+  khi = klo+1;
+  //
+  // Loop from the beginning of the slice upwards to reach the fHalfMax:
+  // With means of bisection:
+  // 
+  x     = (Float_t)klo;
+  a     = 1.;
+  b     = 0.;
+  
+  step = 0.5;
+  Bool_t back = kFALSE;
+  
+  Int_t maxcnt = 20;
+  Int_t cnt    = 0;
+  
+  while (TMath::Abs(y-fHalfMax) > fResolution)
+    {
+      
+      if (back)
+        {
+          x -= step;
+          a += step;
+          b -= step;
+        }
+      else
+        {
+          x += step;
+          a -= step;
+          b += step;
+        }
+      
+      y = a*fLoGainSignal[klo]
+        + b*fLoGainSignal[khi]
+        + (a*a*a-a)*fLoGainSecondDeriv[klo] 
+        + (b*b*b-b)*fLoGainSecondDeriv[khi];
+      
+      if (y > fHalfMax)
+        back = kTRUE;
+      else
+        back = kFALSE;
+      
+      if (++cnt > maxcnt)
+        break;
+      
+      step /= 2.;
+    }
+  
+  time  = x + (Int_t)fLoGainFirst;
+
+  //
+  // Now integrate the whole thing!
+  // 
+  Float_t start = fAbMaxPos - fRiseTime - 0.5;
+  Float_t last  = fAbMaxPos + fFallTime + 0.5;
+  
+  const Int_t diff = int(last) - range;
+  
+  if (diff > 0)
+    {
+      last  -= diff;
+      start -= diff;
+    }
+  CalcIntegralLoGain(sum, start, last);
 }
 
@@ -1248 +1184 @@
       rc = kTRUE;
     }
-  if (IsEnvDefined(env, prefix, "Maximum", print))
-    {
-      b = GetEnvValue(env, prefix, "Maximum", IsExtractionType(kMaximum));
-      if (b) 
-        SetTimeType(kMaximum);
-      rc = kTRUE;
-    }
-  if (IsEnvDefined(env, prefix, "HalfMaximum", print))
-    {
-      b = GetEnvValue(env, prefix, "HalfMaximum", IsExtractionType(kHalfMaximum));
-      if (b) 
-        SetTimeType(kHalfMaximum);
-      rc = kTRUE;
-    }
 
   return MExtractTimeAndCharge::ReadEnv(env, prefix, print) ? kTRUE : rc;