Index: /trunk/MagicSoft/Mars/Changelog
===================================================================
--- /trunk/MagicSoft/Mars/Changelog	(revision 5592)
+++ /trunk/MagicSoft/Mars/Changelog	(revision 5593)
@@ -25,4 +25,8 @@
   * msignal/MExtractor.[h,cc]
     - add possibility to set pedestal pointer from outside
+
+  * msignal/MExtractTimeAndChargeSpline.cc
+    - fixed some smaller bugs for special cases
+
 
  2004/12/14: Markus Gaug and Abelardo Moralejo
Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc
===================================================================
--- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc	(revision 5592)
+++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.cc	(revision 5593)
@@ -325,5 +325,5 @@
   fAbMax        = 0.;
   fAbMaxPos     = 0.;
-  Byte_t maxpos = 0;
+  Int_t  maxpos = 0;
 
   //
@@ -340,10 +340,10 @@
         {
           fAbMax = signal;
-          maxpos =  p-first;
+          maxpos = count;
         }
 
       if (*p++ >= fSaturationLimit)
-        sat++;
-            
+            sat++;
+      
       count++;                                                    
     }
@@ -360,10 +360,9 @@
           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,4 +370,5 @@
             sat++;
 
+          range++;
           logain++;
         }
@@ -395,44 +395,61 @@
     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,5 +459,5 @@
   //
   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,5 +491,4 @@
         }
 
-      //      *fLog << err << x << "  " << y << " " << fAbMaxPos<< endl;
     }
 
@@ -483,5 +499,5 @@
     {
 
-      upper   = (Float_t)maxpos+1.;
+      upper   = 1. + maxpos;
       lower   = (Float_t)maxpos;
       x       = lower;
@@ -508,9 +524,6 @@
               fAbMaxPos = x;
             }
-          //          *fLog << inf << x << "  " << y << " " << fAbMaxPos << endl;
-
         }
   }
-
 
   //
@@ -545,5 +558,4 @@
           fAbMaxPos = x;
         }
-      //      *fLog << inf << x << "  " << y << " " << fAbMaxPos << endl;      
     }
 
@@ -558,10 +570,10 @@
   // 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,5 +598,4 @@
           fAbMaxPos = x;
         }
-      //      *fLog << warn << x << "  " << y << "  " << fAbMaxPos << endl;
     }
 
@@ -674,16 +685,18 @@
       // 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,5 +716,5 @@
   Int_t  count     = 0;
 
-  Float_t pedes        = ped.GetPedestal();
+  const Float_t pedes  = ped.GetPedestal();
   const Float_t ABoffs = ped.GetPedestalABoffset();
 
@@ -712,5 +725,5 @@
   fAbMax        = 0.;
   fAbMaxPos     = 0.;
-  Byte_t maxpos = 0;
+  Int_t  maxpos = 0;
 
   //
@@ -720,18 +733,17 @@
     {
 
-      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,4 +763,5 @@
 
   fLoGainSecondDeriv[range-1] = 0.;
+
   for (Int_t k=range-2;k>=0;k--)
     fLoGainSecondDeriv[k] = fLoGainSecondDeriv[k]*fLoGainSecondDeriv[k+1] + fLoGainFirstDeriv[k];
@@ -756,45 +769,62 @@
     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,5 +832,5 @@
   //
   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,5 +864,4 @@
         }
 
-      //      *fLog << err << x << "  " << y << " " << fAbMaxPos<< endl;
     }
 
@@ -844,5 +873,5 @@
     {
 
-      upper   = (Float_t)maxpos+1.;
+      upper   = 1. + maxpos;
       lower   = (Float_t)maxpos;
       x       = lower;
@@ -869,6 +898,4 @@
               fAbMaxPos = x;
             }
-          //          *fLog << inf << x << "  " << y << " " << fAbMaxPos << endl;
-
         }
   }
@@ -906,5 +933,4 @@
           fAbMaxPos = x;
         }
-      //      *fLog << inf << x << "  " << y << " " << fAbMaxPos << endl;      
     }
 
@@ -919,10 +945,10 @@
   // 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,10 +973,9 @@
           fAbMaxPos = x;
         }
-      //      *fLog << warn << x << "  " << y << "  " << fAbMaxPos << endl;
     }
 
   if (IsExtractionType(kMaximum))
     {
-      time = (Float_t)fLoGainFirst + fAbMaxPos;
+      time = fAbMaxPos + (Int_t)fLoGainFirst;
       dtime = fResolution;
     }
@@ -1020,5 +1045,5 @@
         }
 
-      time  = (Float_t)fLoGainFirst + x;
+      time  = x + (Int_t)fLoGainFirst;
       dtime = fResolution;
     }
@@ -1035,61 +1060,130 @@
       // 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.;
+    }
 }
+
 
 
Index: /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h
===================================================================
--- /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h	(revision 5592)
+++ /trunk/MagicSoft/Mars/msignal/MExtractTimeAndChargeSpline.h	(revision 5593)
@@ -39,5 +39,5 @@
   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,6 +47,6 @@
   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,5 +82,5 @@
                                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
 };