Index: trunk/MagicSoft/Mars/mbase/MMath.cc
===================================================================
--- trunk/MagicSoft/Mars/mbase/MMath.cc	(revision 9225)
+++ trunk/MagicSoft/Mars/mbase/MMath.cc	(revision 9226)
@@ -1,4 +1,4 @@
 /* ======================================================================== *\
-! $Name: not supported by cvs2svn $:$Id: MMath.cc,v 1.43 2009-01-14 12:31:36 tbretz Exp $
+! $Name: not supported by cvs2svn $:$Id: MMath.cc,v 1.44 2009-01-17 14:52:37 tbretz Exp $
 ! --------------------------------------------------------------------------
 !
@@ -20,5 +20,5 @@
 !   Author(s): Thomas Bretz  3/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
 !
-!   Copyright: MAGIC Software Development, 2000-2005
+!   Copyright: MAGIC Software Development, 2000-2009
 !
 !
@@ -50,4 +50,8 @@
 #ifndef ROOT_TComplex
 #include <TComplex.h>
+#endif
+
+#ifndef ROOT_TRandom
+#include <TRandom.h>  // gRandom in RndmExp
 #endif
 
@@ -875,2 +879,13 @@
     e = error.Atof();
 }
+
+Double_t MMath::RndmExp(Double_t tau)
+{
+    // returns an exponential deviate.
+    //
+    //          exp( -t/tau )
+
+    const Double_t x = gRandom->Rndm(); // uniform on ] 0, 1 ]
+
+    return -tau * TMath::Log(x);       // convert to exponential distribution
+}
Index: trunk/MagicSoft/Mars/mbase/MMath.h
===================================================================
--- trunk/MagicSoft/Mars/mbase/MMath.h	(revision 9225)
+++ trunk/MagicSoft/Mars/mbase/MMath.h	(revision 9226)
@@ -87,4 +87,6 @@
 
     void Format(Double_t &v, Double_t &e);
+
+    Double_t RndmExp(Double_t tau);
 }
 
Index: trunk/MagicSoft/Mars/mextralgo/ExtralgoIncl.h
===================================================================
--- trunk/MagicSoft/Mars/mextralgo/ExtralgoIncl.h	(revision 9225)
+++ trunk/MagicSoft/Mars/mextralgo/ExtralgoIncl.h	(revision 9226)
@@ -1,3 +1,5 @@
 #ifndef __CINT__
 
+#include "../mbase/MArrayF.h"
+
 #endif // __CINT__
Index: trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc
===================================================================
--- trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc	(revision 9225)
+++ trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.cc	(revision 9226)
@@ -18,5 +18,5 @@
 !   Author(s): Markus Gaug 09/2004 <mailto:markus@ifae.es>
 !
-!   Copyright: MAGIC Software Development, 2002-2008
+!   Copyright: MAGIC Software Development, 2002-2009
 !
 !
@@ -120,13 +120,13 @@
 // --------------------------------------------------------------------------
 //
-// Returns the highest x value in [min;max[ at which the spline in
-// the bin i is equal to y
-//
-// min and max are defined to be [0;1]
-//
-// The default for min is 0, the default for max is 1
-// The defaule for y is 0
-//
-Double_t MExtralgoSpline::FindY(Int_t i, Bool_t downwards, Double_t y, Double_t min, Double_t max) const
+// Solve the polynomial
+//
+//    y = a*x^3 + b*x^2 + c*x + d'
+//    0 = a*x^3 + b*x^2 + c*x + d' - y
+//
+// to find y in the i-th bin. Return the result as x1, x2, x3 and the return
+// code from MMath::SolvPol3.
+//
+Int_t MExtralgoSpline::SolvePol3(Int_t i, Double_t y, Double_t &x1, Double_t &x2, Double_t &x3) const
 {
     // y = a*x^3 + b*x^2 + c*x + d'
@@ -145,35 +145,59 @@
     //     return -2;
 
+    return MMath::SolvePol3(a, b, c, d, x1, x2, x3);
+}
+
+// --------------------------------------------------------------------------
+//
+// Returns the highest x value in [min;max[ at which the spline in
+// the bin i is equal to y
+//
+// min and max must be in the range [0;1]
+//
+// The default for min is 0, the default for max is 1
+// The default for y is 0
+//
+Double_t MExtralgoSpline::FindYdn(Int_t i, Double_t y, Double_t min, Double_t max) const
+{
     Double_t x1, x2, x3;
-    const Int_t rc = MMath::SolvePol3(a, b, c, d, x1, x2, x3);
-
-    if (downwards==kTRUE)
-    {
-        Double_t x = -1;
-
-        if (rc>0 && x1>=min && x1<max && x1>x)
-            x = x1;
-        if (rc>1 && x2>=min && x2<max && x2>x)
-            x = x2;
-        if (rc>2 && x3>=min && x3<max && x3>x)
-            x = x3;
-
-        return x<0 ? -2 : x+i;
-    }
-    else
-    {
-        Double_t x = 2;
-
-        if (rc>0 && x1>min && x1<=max && x1<x)
-            x = x1;
-        if (rc>1 && x2>min && x2<=max && x2<x)
-            x = x2;
-        if (rc>2 && x3>min && x3<=max && x3<x)
-            x = x3;
-
-        return x>1 ? -2 : x+i;
-    }
-
-    return -2;
+    const Int_t rc = SolvePol3(i, y, x1, x2, x3);
+
+    Double_t x = -1;
+
+    if (rc>0 && x1>=min && x1<max && x1>x)
+        x = x1;
+    if (rc>1 && x2>=min && x2<max && x2>x)
+        x = x2;
+    if (rc>2 && x3>=min && x3<max && x3>x)
+        x = x3;
+
+    return x<0 ? -2 : x+i;
+}
+
+// --------------------------------------------------------------------------
+//
+// Returns the lowest x value in [min;max[ at which the spline in
+// the bin i is equal to y
+//
+// min and max must be in the range [0;1]
+//
+// The default for min is 0, the default for max is 1
+// The default for y is 0
+//
+Double_t MExtralgoSpline::FindYup(Int_t i, Double_t y, Double_t min, Double_t max) const
+{
+    Double_t x1, x2, x3;
+    const Int_t rc = SolvePol3(i, y, x1, x2, x3);
+
+    Double_t x = 2;
+
+    if (rc>0 && x1>min && x1<=max && x1<x)
+        x = x1;
+    if (rc>1 && x2>min && x2<=max && x2<x)
+        x = x2;
+    if (rc>2 && x3>min && x3<=max && x3<x)
+        x = x3;
+
+    return x>1 ? -2 : x+i;
 }
 
@@ -181,5 +205,5 @@
 //
 // Search analytically downward for the value y of the spline, starting
-// at x, until x==0. If y is not found -2 is returned.
+// at x, until x==0. If y is not found or out of range -2 is returned.
 //
 Double_t MExtralgoSpline::SearchYdn(Float_t x, Float_t y) const
@@ -192,11 +216,16 @@
         return -2;
 
-    Double_t rc = FindY(i, kTRUE, y, 0, x-i);
+    Double_t rc = FindYdn(i, y, 0, x-i);
     while (--i>=0 && rc<0)
-        rc = FindY(i, kTRUE, y);
+        rc = FindYdn(i, y);
 
     return rc;
 }
 
+// --------------------------------------------------------------------------
+//
+// Search analytically upwards for the value y of the spline, starting
+// at x, until x==fNum-1. If y is not found or out of range -2 is returned.
+//
 Double_t MExtralgoSpline::SearchYup(Float_t x, Float_t y) const
 {
@@ -208,7 +237,7 @@
         return -2;
 
-    Double_t rc = FindY(i, kFALSE, y, x-i, 1.);
+    Double_t rc = FindYup(i, y, x-i, 1.);
     while (++i<fNum-1 && rc<0)
-        rc = FindY(i, kFALSE, y);
+        rc = FindYup(i, y);
 
     return rc;
@@ -222,8 +251,4 @@
 Float_t MExtralgoSpline::CalcIntegral(Float_t pos) const
 {
-    // In the future we will calculate the intgeral analytically.
-    // It has been tested that it gives identical results within
-    // acceptable differences.
-
     // We allow extrapolation of 1/2 slice.
     const Float_t min = fRiseTime;        //-0.5+fRiseTime;
@@ -236,4 +261,25 @@
 
     return EvalInteg(pos-fRiseTime, pos+fFallTime);
+}
+
+MArrayF MExtralgoSpline::GetIntegral(bool norm) const
+{
+    MArrayF val(fNum);
+
+    //val[0] = 0;
+
+    Double_t integ = 0;
+    for (int i=0; i<fNum-1; i++)
+    {
+        integ += EvalInteg(i);
+
+        val[i+1] = integ;
+    }
+
+    if (norm)
+        for (int i=0; i<fNum-1; i++)
+            val[i+1] /= val[fNum-1];
+
+    return val;
 }
 
Index: trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h
===================================================================
--- trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h	(revision 9225)
+++ trunk/MagicSoft/Mars/mextralgo/MExtralgoSpline.h	(revision 9226)
@@ -6,4 +6,5 @@
 #endif
 
+class MArrayF;
 class TComplex;
 
@@ -75,5 +76,8 @@
     }
 
-    Double_t FindY(Int_t i, Bool_t downwards, Double_t y=0, Double_t min=0, Double_t max=1) const;
+    Int_t SolvePol3(Int_t i, Double_t y, Double_t &x1, Double_t &x2, Double_t &x3) const;
+    Double_t FindYdn(Int_t i, Double_t y=0, Double_t min=0, Double_t max=1) const;
+    Double_t FindYup(Int_t i, Double_t y=0, Double_t min=0, Double_t max=1) const;
+    //Double_t FindY(Int_t i, Bool_t downwards, Double_t y=0, Double_t min=0, Double_t max=1) const;
 
     Int_t EvalDerivEq0(const Int_t i, Double_t &x1, Double_t &x2) const;
@@ -153,5 +157,7 @@
     }
 
-    // Identical to sum EvalInteg(i, 0, 1) for i=0 to i<b but much faster
+    // Identical to sum of EvalInteg(i, 0, 1) for i=a to i=b-1,
+    // but much faster
+    // It is identical to EvalInteg(fVal[a], fVal[b])
     // Be carefull: NO RANGECHECK!
     inline Double_t EvalInteg(Int_t a, Int_t b) const
@@ -329,4 +335,6 @@
     Double_t SearchYdn(Float_t y) const { return SearchYdn(fNum, y); }
     Double_t SearchYup(Float_t y) const { return SearchYup(0,    y); }
+
+    MArrayF GetIntegral(bool norm=true) const;
 };
 
Index: trunk/MagicSoft/Mars/mpedestal/MPedestalSubtractedEvt.cc
===================================================================
--- trunk/MagicSoft/Mars/mpedestal/MPedestalSubtractedEvt.cc	(revision 9225)
+++ trunk/MagicSoft/Mars/mpedestal/MPedestalSubtractedEvt.cc	(revision 9226)
@@ -1,4 +1,4 @@
 /* ======================================================================== *\
-! $Name: not supported by cvs2svn $:$Id: MPedestalSubtractedEvt.cc,v 1.6 2007-06-16 22:08:00 tbretz Exp $
+! $Name: not supported by cvs2svn $:$Id: MPedestalSubtractedEvt.cc,v 1.7 2009-01-17 14:52:41 tbretz Exp $
 ! --------------------------------------------------------------------------
 !
@@ -32,6 +32,7 @@
 //
 /////////////////////////////////////////////////////////////////////////////
-
 #include "MPedestalSubtractedEvt.h"
+
+#include "MLogManip.h"
 
 ClassImp(MPedestalSubtractedEvt);
@@ -185,4 +186,28 @@
 }
 
+void MPedestalSubtractedEvt::Print(Option_t *o) const
+{
+    *fLog << all << GetDescriptor() << endl;
+    *fLog << " Num Pixels:  " << fNumPixels << " (" << fNumSamples << " samples)" << endl;
+    *fLog << " Samples raw:" << hex << endl;;
+    for (UInt_t idx=0; idx<fNumPixels; idx++)
+    {
+        *fLog << setw(4) << dec << idx << hex << ":";
+        for (UInt_t i=0; i<fNumSamples; i++)
+            *fLog << " " << fSamplesRaw[idx*fNumSamples+i];
+        *fLog << endl;
+    }
+    *fLog << dec << endl;
+    *fLog << " Samples:" << endl;;
+    for (UInt_t idx=0; idx<fNumPixels; idx++)
+    {
+        *fLog << setw(4) << idx << ":";
+        for (UInt_t i=0; i<fNumSamples; i++)
+            *fLog << " " << fSamples[idx*fNumSamples+i];
+        *fLog << endl;
+    }
+    *fLog << endl;
+}
+
 /*
 #include <TSpline.h>
Index: trunk/MagicSoft/Mars/mpedestal/MPedestalSubtractedEvt.h
===================================================================
--- trunk/MagicSoft/Mars/mpedestal/MPedestalSubtractedEvt.h	(revision 9225)
+++ trunk/MagicSoft/Mars/mpedestal/MPedestalSubtractedEvt.h	(revision 9226)
@@ -20,9 +20,9 @@
 {
 private:
-    MArrayF fSamples;        // list of all samples with pedestal subtracted
-    MArrayS fSamplesRaw;     // list of all samples (raw)
+    MArrayF fSamples;         // list of all samples with pedestal subtracted
+    MArrayS fSamplesRaw;      // list of all samples (raw)
 
-    UInt_t fNumSamples;      // number of samples per pixel
-    UInt_t fNumPixels;       // number of pixels
+    UInt_t fNumSamples;       // number of samples per pixel
+    UInt_t fNumPixels;        // number of pixels
 
 public:
@@ -114,4 +114,6 @@
     }
 
+    void Print(Option_t *o="") const;
+
     Bool_t GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type=0) const { return kTRUE; }
     void   DrawPixelContent(Int_t num) const { }
Index: trunk/MagicSoft/Mars/mraw/MRawEvtData.cc
===================================================================
--- trunk/MagicSoft/Mars/mraw/MRawEvtData.cc	(revision 9225)
+++ trunk/MagicSoft/Mars/mraw/MRawEvtData.cc	(revision 9226)
@@ -284,4 +284,5 @@
     name += pix.GetPixelId();
 
+
     Bool_t same = str.Contains("same");
 
@@ -392,4 +393,14 @@
 }
 
+void MRawEvtData::ResetPixels()
+{
+    if (!fRunHeader)
+        return;
+
+    // FIXME: Better give NumNormalPixels if numconnected==0 ?
+
+    ResetPixels(fRunHeader->GetNumNormalPixels(), fRunHeader->GetNumNormalPixels()-1);
+}
+
 // --------------------------------------------------------------------------
 //
@@ -435,5 +446,15 @@
 void MRawEvtData::Set(const MArrayI &data)
 {
-    const UInt_t n = fRunHeader->GetNumSamples()*fConnectedPixels;
+    fConnectedPixels = fHiGainPixId->GetSize();
+
+    UInt_t n = fConnectedPixels*fRunHeader->GetNumSamplesHiGain();
+    if (n!=data.GetSize())
+    {
+        *fLog << err << "MRawEvtData::Set: Size mismatch." << endl;
+        *fLog << " fConnectedPixels="   << fConnectedPixels << endl;
+        *fLog << " NumHiGainSamples="   << fRunHeader->GetNumSamplesHiGain() << endl;
+        *fLog << " data.GetSize()="     << data.GetSize() << endl;
+        return;
+    }
 
     Byte_t *dest = fHiGainFadcSamples->GetArray();
@@ -448,5 +469,5 @@
             Byte_t *ptr = reinterpret_cast<Byte_t*>(dest);
             while (src<end)
-                *ptr++ = *src++;
+                *ptr++ = Byte_t(*src++);
         }
         return;
@@ -456,5 +477,5 @@
             UShort_t *ptr = reinterpret_cast<UShort_t*>(dest);
             while (src<end)
-                *ptr++ = *src++;
+                *ptr++ = UShort_t(*src++);
         }
         return;
@@ -464,4 +485,18 @@
         return;
     }
+}
+
+void MRawEvtData::SetIndices(const MArrayS &idx)
+{
+    if (idx.GetSize()!=fHiGainPixId->GetSize())
+        return;
+
+    memcpy(fHiGainPixId->GetArray(), idx.GetArray(), idx.GetSize()*sizeof(UShort_t));
+}
+
+void MRawEvtData::SetIndices()
+{
+    for (UInt_t i=0; i<fHiGainPixId->GetSize(); i++)
+        (*fHiGainPixId)[i] = i;
 }
 
Index: trunk/MagicSoft/Mars/mraw/MRawEvtData.h
===================================================================
--- trunk/MagicSoft/Mars/mraw/MRawEvtData.h	(revision 9225)
+++ trunk/MagicSoft/Mars/mraw/MRawEvtData.h	(revision 9226)
@@ -57,4 +57,5 @@
     }
 
+
 public:
     MRawEvtData(const char *name=NULL, const char *title=NULL);
@@ -74,7 +75,10 @@
     void Draw (Option_t * = NULL);
 
+    void ResetPixels();
     void ResetPixels(UShort_t npix, UShort_t maxid);
     void AddPixel(UShort_t nOfPixel, const TArrayC &data);
     void Set(const MArrayI &data);
+    void SetIndices(const MArrayS &idx);
+    void SetIndices();
 
     UShort_t GetNumHiGainSamples() const;