Changeset 5593 for trunk/MagicSoft/Mars/msignal

Timestamp:

12/14/04 17:15:28 (20 years ago)

Author:

gaug

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/msignal

Files:

• MExtractTimeAndChargeSpline.cc (modified) (27 diffs)
• MExtractTimeAndChargeSpline.h (modified) (3 diffs)

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc

@@ -325 +325 @@
   fAbMax        = 0.;
   fAbMaxPos     = 0.;
-  Byte_t maxpos = 0;
+  Int_t  maxpos = 0;
 
   //
@@ -340 +340 @@
         {
           fAbMax = signal;
-          maxpos =  p-first;
+          maxpos = count;
         }
 
       if (*p++ >= fSaturationLimit)
-        sat++;
-            
+            sat++;
+      
       count++;                                                    
     }
@@ -360 +360 @@
           const Float_t signal = (Float_t)*logain - pedmean[(ids+abflag) & 0x1];
           fHiGainSignal[range] = signal;
-          range++;
           
           if (signal > fAbMax)
             { 
               fAbMax = signal;
-              maxpos =  logain-first;
+              maxpos = range;
             }
           
@@ -371 +370 @@
             sat++;
 
+          range++;
           logain++;
         }
@@ -395 +395 @@
     fHiGainSecondDeriv[k] /= 6.;
   
-  if (IsNoiseCalculation() && IsExtractionType(kAmplitude))
-    {
-      if (fRandomIter == (TMath::Floor(1./fResolution)))
+  if (IsNoiseCalculation())
+    {
+      if (fRandomIter == int(1./fResolution))
         fRandomIter = 0;
-
-      const Float_t b = fRandomIter * fResolution;
-      const Float_t a = 1. - b;
-
+      
+      const Float_t nsx = fRandomIter * fResolution;
+
+      if (IsExtractionType(kAmplitude))
+        {
+          const Float_t b = nsx;
+          const Float_t a = 1. - nsx;
+          
+          sum = a*fHiGainSignal[1]
+            + b*fHiGainSignal[2]
+            + (a*a*a-a)*fHiGainSecondDeriv[1] 
+            + (b*b*b-b)*fHiGainSecondDeriv[2];
+        }
+      else
+        {
+          Float_t start = 2. + nsx;
+          Float_t last  = start + fRiseTime + fFallTime;
+      
+          if (int(last) > range)
+            {
+              const Int_t diff = range - int(last);
+              last  -= diff;
+              start -= diff;
+            }
+          
+          CalcIntegralHiGain(sum, start, last);
+        }
       fRandomIter++;
-
-      sum = a*fHiGainSignal[1]
-          + b*fHiGainSignal[2]
-          + (a*a*a-a)*fHiGainSecondDeriv[1] 
-          + (b*b*b-b)*fHiGainSecondDeriv[2];
       return;
     }
 
-  if (IsNoiseCalculation() && IsExtractionType(kIntegral))
-    {
-      //
-      // Take the spline value at the middle of the third slice (to avoid egde effects)
-      // 
-      Int_t first = 2;
-      Int_t last  = first + (Int_t)(fRiseTime+fFallTime);
-      CalcIntegralHiGain(sum,first,last);
-      return;
-    }
-
-  //
-  // Allow no saturated slice 
+  //
+  // Allow one saturated slice 
   // and 
-  // Don't start if the maxpos is too close to the left limit.
-  //
-  if ((sat || maxpos < 2))
+  // Don't start if the maxpos is too close to the limits.
+  //
+  if (sat > 1 || maxpos < 2 || maxpos > range-2)
     {
       time =  IsExtractionType(kMaximum) 
         ? (Float_t)(fHiGainFirst + maxpos) 
         : (Float_t)(fHiGainFirst + maxpos - 1);
-      sum  =  IsExtractionType(kAmplitude) 
-        ? fAbMax : 0.;
+
+      if (IsExtractionType(kAmplitude))
+        {
+          sum = fAbMax;
+          return;
+        }
+      
+      if (maxpos > range - 2)
+        CalcIntegralHiGain(sum, (Float_t)range - fRiseTime - fFallTime, (Float_t)range - 0.001);
+      else
+        CalcIntegralHiGain(sum, 0.001, fRiseTime + fFallTime);
+
       return;
     }
@@ -442 +459 @@
   //
   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 lower   = -1. + maxpos;
   Float_t upper   = (Float_t)maxpos;
   fAbMaxPos       = upper;
@@ -474 +491 @@
         }
 
-      //      *fLog << err << x << "  " << y << " " << fAbMaxPos<< endl;
     }
 
@@ -483 +499 @@
     {
 
-      upper   = (Float_t)maxpos+1.;
+      upper   = 1. + maxpos;
       lower   = (Float_t)maxpos;
       x       = lower;
@@ -508 +524 @@
               fAbMaxPos = x;
             }
-          //          *fLog << inf << x << "  " << y << " " << fAbMaxPos << endl;
-
         }
   }
-
 
   //
@@ -545 +558 @@
           fAbMaxPos = x;
         }
-      //      *fLog << inf << x << "  " << y << " " << fAbMaxPos << endl;      
     }
 
@@ -558 +570 @@
   // which requires new setting of klocont and khicont
   //
-  if (x < klo + fResolution/2.)
+  if (x < lower + fResolution/2.)
     {
       klo--;
       khi--;
-      upper--;
-      lower--;
+      upper += 1.;
+      lower -= 1.;
     }
 
@@ -586 +598 @@
           fAbMaxPos = x;
         }
-      //      *fLog << warn << x << "  " << y << "  " << fAbMaxPos << endl;
     }
 
@@ -674 +685 @@
       // Now integrate the whole thing!
       // 
-      Int_t startslice = (Int_t)(fAbMaxPos - fRiseTime);
-      Int_t lastslice  = (Int_t)(fAbMaxPos + fFallTime);
-      
-      if (lastslice > range)
-        {
-          lastslice = range;
-          startslice += (lastslice - range);
-        }
-      
-      CalcIntegralHiGain(sum, startslice, lastslice);
-    }
-  
+
+      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);
+    }
 }
 
@@ -703 +716 @@
   Int_t  count     = 0;
 
-  Float_t pedes        = ped.GetPedestal();
+  const Float_t pedes  = ped.GetPedestal();
   const Float_t ABoffs = ped.GetPedestalABoffset();
 
@@ -712 +725 @@
   fAbMax        = 0.;
   fAbMaxPos     = 0.;
-  Byte_t maxpos = 0;
+  Int_t  maxpos = 0;
 
   //
@@ -720 +733 @@
     {
 
-      const Int_t ids      = fLoGainFirst + count ;
+      const Int_t ids      =  count + fLoGainFirst;
       const Float_t signal = (Float_t)*p - pedmean[(ids+abflag) & 0x1];
       fLoGainSignal[count] = signal;
 
-      if (signal > fAbMax)
+      if (signal > fAbMax + 0.1)
         {
           fAbMax = signal;
-          maxpos =  p-first;
-        }
-
-      if (*p >= fSaturationLimit)
+          maxpos = count;
+        }
+
+      if (*p++ >= fSaturationLimit)
         sat++;
-            
-      p++;
+      
       count++;                                                    
     }
@@ -751 +763 @@
 
   fLoGainSecondDeriv[range-1] = 0.;
+
   for (Int_t k=range-2;k>=0;k--)
     fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k];
@@ -756 +769 @@
     fLoGainSecondDeriv[k] /= 6.;
   
-  if (IsNoiseCalculation() && IsExtractionType(kAmplitude))
-    {
-      //
-      // Take the spline value at the middle of the third slice (to avoid egde effects)
-      // 
-      sum = 0.5*fLoGainSignal[2]
-          + 0.5*fLoGainSignal[3]
-          + (-0.375)*fLoGainSecondDeriv[2] 
-          + (-0.375)*fLoGainSecondDeriv[3];
+  if (IsNoiseCalculation())
+    {
+      if (fRandomIter == int(1./fResolution))
+        fRandomIter = 0;
+      
+      const Float_t nsx = fRandomIter * fResolution;
+
+      if (IsExtractionType(kAmplitude))
+        {
+          const Float_t b = nsx;
+          const Float_t a = 1. - nsx;
+          
+          sum = a*fLoGainSignal[1]
+            + b*fLoGainSignal[2]
+            + (a*a*a-a)*fLoGainSecondDeriv[1] 
+            + (b*b*b-b)*fLoGainSecondDeriv[2];
+        }
+      else
+        {
+          Float_t start = 2. + nsx;
+          Float_t last  = start + fRiseTime + fFallTime;
+      
+          if (int(last) > range)
+            {
+              const Int_t diff = range - int(last);
+              last  -= diff;
+              start -= diff;
+            }
+          
+          CalcIntegralLoGain(sum, start, last);
+        }
+      fRandomIter++;
       return;
     }
-
-  if (IsNoiseCalculation() && IsExtractionType(kIntegral))
-    {
-      //
-      // Take the spline value at the middle of the third slice (to avoid egde effects)
-      // 
-      Int_t first = 2;
-      Int_t last  = first + (Int_t)(fRiseTime+fFallTime);
-      CalcIntegralLoGain(sum,first,last);
-      return;
-    }
-
   //
   // Allow no saturated slice 
   // and 
-  // Don't start if the maxpos is too close to the left limit.
-  //
-  if (sat || maxpos < 1)
+  // 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 (IsExtractionType(kAmplitude))
+        {
+          sum = fAbMax;
+          return;
+        }
+      
+      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;
     }
-      
-  if (maxpos < 2 && IsExtractionType(kHalfMaximum))
-    {
-      time = (Float_t)(fLoGainFirst + maxpos - 1);
-      return;
-    }
 
   //
@@ -802 +832 @@
   //
   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 lower   = -1. + maxpos;
   Float_t upper   = (Float_t)maxpos;
   fAbMaxPos       = upper;
@@ -834 +864 @@
         }
 
-      //      *fLog << err << x << "  " << y << " " << fAbMaxPos<< endl;
     }
 
@@ -844 +873 @@
     {
 
-      upper   = (Float_t)maxpos+1.;
+      upper   = 1. + maxpos;
       lower   = (Float_t)maxpos;
       x       = lower;
@@ -869 +898 @@
               fAbMaxPos = x;
             }
-          //          *fLog << inf << x << "  " << y << " " << fAbMaxPos << endl;
-
         }
   }
@@ -906 +933 @@
           fAbMaxPos = x;
         }
-      //      *fLog << inf << x << "  " << y << " " << fAbMaxPos << endl;      
     }
 
@@ -919 +945 @@
   // which requires new setting of klocont and khicont
   //
-  if (x < klo + fResolution/2.)
+  if (x < lower + fResolution/2.)
     {
       klo--;
       khi--;
-      upper--;
-      lower--;
+      upper += 1.;
+      lower -= 1.;
     }
 
@@ -947 +973 @@
           fAbMaxPos = x;
         }
-      //      *fLog << warn << x << "  " << y << "  " << fAbMaxPos << endl;
     }
 
   if (IsExtractionType(kMaximum))
     {
-      time = (Float_t)fLoGainFirst + fAbMaxPos;
+      time = fAbMaxPos + (Int_t)fLoGainFirst;
       dtime = fResolution;
     }
@@ -1020 +1045 @@
         }
 
-      time  = (Float_t)fLoGainFirst + x;
+      time  = x + (Int_t)fLoGainFirst;
       dtime = fResolution;
     }
@@ -1035 +1060 @@
       // Now integrate the whole thing!
       // 
-      Int_t startslice = (Int_t)(fAbMaxPos - fRiseTime);
-      Int_t lastslice  = (Int_t)(fAbMaxPos + fFallTime + 1);
-      
-      if (lastslice > range)
-        {
-          lastslice = range;
-          startslice += (lastslice - range);
-        }
-      CalcIntegralLoGain(sum, startslice, lastslice);
+      Float_t start = fAbMaxPos - fRiseTime;
+      Float_t last  = fAbMaxPos + fFallTime + 1.;
+      
+      const Int_t diff = int(last) - range;
+
+      if (diff > 0)
+        {
+          last  -= diff;
+          start -= diff;
+        }
+      CalcIntegralLoGain(sum, start, last);
     }
 }
 
-void MExtractTimeAndChargeSpline::CalcIntegralHiGain(Float_t &sum, Int_t startslice, Int_t lastslice)
+void MExtractTimeAndChargeSpline::CalcIntegralHiGain(Float_t &sum, Float_t start, Float_t last)
 {
 
-  if (startslice < 0)
-    {
-      lastslice -= startslice;
-      startslice = 0;
-    }
-  
-  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];
-  
-  sum += 0.5*fHiGainSignal[lastslice];
-
+  const Float_t step = 0.1;
+
+  if (start < 0)
+    {
+      last -= start;
+      start = 0.;
+    }
+  
+  Int_t klo = int(start);
+  Int_t khi = klo+1;
+
+  Float_t up = TMath::Ceil(start);
+  Float_t lo = TMath::Floor(start);
+
+  const Int_t m = int((start-klo)/step);
+  start = step*m + klo; // Correct start for the digitization due to resolution
+
+  Float_t x = start;
+  Float_t a = up-start;
+  Float_t b = start-lo;
+
+  while (1)
+    {
+      
+      while (x<up)
+        {
+          x += step;
+
+          if (x > last)
+            {
+              sum *= step;
+              return;
+            }
+          
+          a -= step;
+          b += step;
+          
+          sum += a*fHiGainSignal[klo]
+          + b*fHiGainSignal[khi]
+            + (a*a*a-a)*fHiGainSecondDeriv[klo] 
+            + (b*b*b-b)*fHiGainSecondDeriv[khi];
+        }
+
+      up += 1.;
+      lo += 1.;
+      klo++;
+      khi++;
+      start += 1.;
+      a = 1.;
+      b = 0.;
+    }
+  
 }
-
-void MExtractTimeAndChargeSpline::CalcIntegralLoGain(Float_t &sum, Int_t startslice, Int_t lastslice)
+void MExtractTimeAndChargeSpline::CalcIntegralLoGain(Float_t &sum, Float_t start, Float_t last)
 {
 
-  if (startslice < 0)
-    {
-      lastslice -= startslice;
-      startslice = 0;
-    }
-  
-  Int_t i = startslice;
-  sum = 0.5*fLoGainSignal[i];
-  
-  //
-  // We sum 1.5 times the second deriv. coefficients because these had been 
-  // divided by 6. already. Usually, 0.25*fLoGainSecondDeriv should be added.
-  //
-  for (i=startslice+1; i<lastslice; i++)
-    sum += fLoGainSignal[i] + 1.5*fLoGainSecondDeriv[i];
-  
-  sum += 0.5*fLoGainSignal[lastslice];
-
+  const Float_t step = 0.1;
+
+  if (start < 0)
+    {
+      last -= start;
+      start = 0.;
+    }
+  
+  Int_t klo = int(start);
+  Int_t khi = klo+1;
+
+  Float_t up = TMath::Ceil(start);
+  Float_t lo = TMath::Floor(start);
+
+  const Int_t m = int((start-klo)/step);
+  start = step*m + klo; // Correct start for the digitization due to resolution
+
+  Float_t x = start;
+  Float_t a = up-start;
+  Float_t b = start-lo;
+
+  while (1)
+    {
+      
+      while (x<up)
+        {
+          x += step;
+          
+          if (x > last)
+            {
+              sum *= step;
+              return;
+            }
+          
+          a -= step;
+          b += step;
+          
+          sum += a*fHiGainSignal[klo]
+          + b*fHiGainSignal[khi]
+            + (a*a*a-a)*fHiGainSecondDeriv[klo] 
+            + (b*b*b-b)*fHiGainSecondDeriv[khi];
+        }
+
+      up += 1.;
+      lo += 1.;
+      klo++;
+      khi++;
+      start += 1.;
+      a = 1.;
+      b = 0.;
+    }
 }
+
 
 

trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h

@@ -39 +39 @@
   Float_t fFallTime;                     // The usual fall time of the pulse
 
-  UInt_t  fRandomIter;                   // Counter used to randomize weights for noise calculation
+  Int_t   fRandomIter;                   // Counter used to randomize weights for noise calculation
   
   Byte_t  fFlags;                        // Bit-field to hold the time extraction types
@@ -47 +47 @@
   Bool_t InitArrays();
   
-  void CalcIntegralHiGain(Float_t &sum, Int_t startslice, Int_t lastslice);
-  void CalcIntegralLoGain(Float_t &sum, Int_t startslice, Int_t lastslice);
+  void CalcIntegralHiGain(Float_t &sum, Float_t start, Float_t last);
+  void CalcIntegralLoGain(Float_t &sum, Float_t start, Float_t last);
 
 public:
@@ -82 +82 @@
                                Byte_t &sat, const MPedestalPix &ped, const Bool_t abflag);
 
-  ClassDef(MExtractTimeAndChargeSpline, 0)   // Task to Extract Arrival Times and Charges using a Fast Cubic Spline
+  ClassDef(MExtractTimeAndChargeSpline, 1)   // Task to Extract Arrival Times and Charges using a Fast Cubic Spline
 };