Index: /trunk/Mars/msignal/MTreatSaturation.cc
===================================================================
--- /trunk/Mars/msignal/MTreatSaturation.cc	(revision 17858)
+++ /trunk/Mars/msignal/MTreatSaturation.cc	(revision 17858)
@@ -0,0 +1,151 @@
+/* ======================================================================== *\
+!
+! *
+! * This file is part of MARS, the MAGIC Analysis and Reconstruction
+! * Software. It is distributed to you in the hope that it can be a useful
+! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
+! * It is distributed WITHOUT ANY WARRANTY.
+! *
+! * Permission to use, copy, modify and distribute this software and its
+! * documentation for any purpose is hereby granted without fee,
+! * provided that the above copyright notice appear in all copies and
+! * that both that copyright notice and this permission notice appear
+! * in supporting documentation. It is provided "as is" without express
+! * or implied warranty.
+! *
+!
+!
+!   Author(s): Thomas Bretz, 2014 <mailto:tbretz@astro.uni-wuerzburg.de>
+!
+!   Copyright: MAGIC Software Development, 2014
+!
+!
+\* ======================================================================== */
+
+//////////////////////////////////////////////////////////////////////////////
+//
+//   MTreatSaturation
+//
+//////////////////////////////////////////////////////////////////////////////
+#include "MTreatSaturation.h"
+
+#include <algorithm>
+
+#include <TRandom.h>
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+#include "MParList.h"
+
+#include "MRawRunHeader.h"
+#include "MExtralgoSpline.h"
+
+#include "MArrivalTimeCam.h"
+#include "MArrivalTimePix.h"
+
+#include "MExtractedSignalCam.h"
+#include "MExtractedSignalPix.h"
+
+#include "MPedestalSubtractedEvt.h"
+
+ClassImp(MTreatSaturation);
+
+using namespace std;
+
+// --------------------------------------------------------------------------
+//
+// Default constructor. 
+//
+MTreatSaturation::MTreatSaturation(const char *name, const char *title)
+    : fRaw(0), fEvt(0)
+{
+    fName  = name  ? name  : "MTreatSaturation";
+    fTitle = title ? title : "Base class for replacing saturation with pulses";
+}
+
+// --------------------------------------------------------------------------
+//
+Int_t MTreatSaturation::PreProcess(MParList *pList)
+{
+    fRaw = (MRawEvtData*)pList->FindObject("MRawEvtData");
+    if (!fRaw)
+    {
+        *fLog << err << "MRawEvtData not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fEvt = (MPedestalSubtractedEvt*)pList->FindObject("MPedestalSubtractedEvt");
+    if (!fEvt)
+    {
+        *fLog << err << "MPedestalSubtractedEvt not found... aborting." << endl;
+        return kFALSE;
+    }
+    return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+Int_t MTreatSaturation::Process()
+{
+    const UInt_t fHiGainFirst =  50;
+    const UInt_t fHiGainLast  = 225;
+
+    const UInt_t npix = fEvt->GetNumPixels();
+
+    const UInt_t rangehi = fHiGainLast - fHiGainFirst - 1;
+
+    if (fDev1.GetSize()!=rangehi)
+    {
+        fDev1.Set(rangehi);
+        fDev2.Set(rangehi);
+    }
+
+    for (UInt_t ch=0; ch<npix; ch++)
+    {
+        Float_t *ptr = fEvt->GetSamples(ch);
+
+        const Float_t *pbeg = ptr+fHiGainFirst+1;
+        const Float_t *pend = ptr+fHiGainLast-1;
+
+        // Find maximum
+        const Float_t *pmax = pbeg;
+
+        for (const Float_t *p=pbeg; p<pend; p++)
+            if (*p>*pmax)
+                pmax = p;
+
+        if (*pmax<1800)
+            continue;
+
+        MExtralgoSpline s(pbeg, rangehi, fDev1.GetArray(), fDev2.GetArray());
+        s.SetExtractionType(MExtralgoSpline::kAmplitudeRel);
+
+        const double leading   = s.SearchYdn(1800., pmax-pbeg);
+        const double falling   = s.SearchYup(1800., pmax-pbeg);
+        const double width     = falling-leading;
+        const double amplitude = 1800/(0.898417 - 0.0187633*width + 0.000163919*width*width - 6.87369e-7*width*width*width + 1.13264e-9*width*width*width*width);
+        const double time      = leading-(4.46944 - 0.277272*width + 0.0035479*width*width);
+
+        double baseline = 0;
+        for (const Float_t *p=ptr+10; p<ptr+50; p++)
+            baseline += *ptr;
+        baseline /= 40;
+
+        const UInt_t beg = TMath::CeilNint(leading);
+        const UInt_t end = TMath::FloorNint(falling);
+
+        ptr += fHiGainFirst+1;
+        for (UInt_t i=0; i<end-beg; i++)
+        {
+            cout << i << " " << ptr[beg-i] << " ";
+            ptr[beg+i] = amplitude*(1-1/(1+exp((i-time)/2.14)))*exp((time-i)/38.8)+baseline;
+            cout << ptr[beg+i] << endl;
+        }
+
+        cout << endl;
+    }
+
+    return kTRUE;
+}
+
Index: /trunk/Mars/msignal/MTreatSaturation.h
===================================================================
--- /trunk/Mars/msignal/MTreatSaturation.h	(revision 17858)
+++ /trunk/Mars/msignal/MTreatSaturation.h	(revision 17858)
@@ -0,0 +1,35 @@
+#ifndef MARS_MTreatSaturation
+#define MARS_MTreatSaturation
+
+#ifndef MARS_MArrayF
+#include "MArrayF.h"
+#endif
+#ifndef MARS_MTask
+#include "MTask.h"
+#endif
+
+class MRawEvtData;
+class MPedestalSubtractedEvt;
+
+class MTreatSaturation : public MTask
+{
+private:
+    MRawEvtData            *fRaw;
+    MPedestalSubtractedEvt *fEvt;
+
+    MArrayF fDev1;
+    MArrayF fDev2;
+
+    Int_t  PreProcess(MParList *pList);
+    Int_t  Process();
+
+public:
+    MTreatSaturation(const char *name=NULL, const char *title=NULL);
+
+    // For MExtractPedestal
+    Bool_t ReInit(MParList *pList);
+
+    ClassDef(MTreatSaturation, 0)   // Replace saturating samples with pulse
+};
+
+#endif