Index: /trunk/MagicSoft/Mars/manalysis/MCerPhotAnal2.cc
===================================================================
--- /trunk/MagicSoft/Mars/manalysis/MCerPhotAnal2.cc	(revision 2249)
+++ /trunk/MagicSoft/Mars/manalysis/MCerPhotAnal2.cc	(revision 2249)
@@ -0,0 +1,237 @@
+/* ======================================================================== *\
+!
+! *
+! * 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  12/2000 <mailto:tbretz@uni-sw.gwdg.de>
+!
+!   Copyright: MAGIC Software Development, 2000-2001
+!
+!
+\* ======================================================================== */
+
+//////////////////////////////////////////////////////////////////////////////
+//                                                                          //
+//   MCerPhotAnal2                                                           //
+//                                                                          //
+//   This is a task which calculates the number of photons from the FADC    //
+//   time slices. At the moment it integrates simply the FADC values.       //
+//                                                                          //
+//  Input Containers:                                                       //
+//   MRawEvtData, MPedesdtalCam                                             //
+//                                                                          //
+//  Output Containers:                                                      //
+//   MCerPhotEvt                                                            //
+//                                                                          //
+//////////////////////////////////////////////////////////////////////////////
+
+#include "MCerPhotAnal2.h"
+
+#include "MParList.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+#include "MRawRunHeader.h"
+#include "MRawEvtData.h"       // MRawEvtData::GetNumPixels
+#include "MCerPhotEvt.h"
+#include "MPedestalPix.h"
+#include "MPedestalCam.h"
+#include "MRawEvtPixelIter.h"
+
+ClassImp(MCerPhotAnal2);
+
+using namespace std;
+
+// --------------------------------------------------------------------------
+//
+// Default constructor. b is the number of slices before the maximum slice,
+// a the number of slices behind the maximum slice which is taken as signal.
+//
+MCerPhotAnal2::MCerPhotAnal2(Byte_t b, Byte_t a, const char *name, const char *title)
+    : fBefore(b), fAfter(a)
+{
+    fName  = name  ? name  : "MCerPhotAnal2";
+    fTitle = title ? title : "Task to calculate Cerenkov photons from raw data";
+
+    AddToBranchList("MRawEvtData.fHiGainPixId");
+    AddToBranchList("MRawEvtData.fLoGainPixId");
+    AddToBranchList("MRawEvtData.fHiGainFadcSamples");
+    AddToBranchList("MRawEvtData.fLoGainFadcSamples");
+}
+
+// --------------------------------------------------------------------------
+//
+// The PreProcess searches for the following input containers:
+//  - MRawRunHeader
+//  - MRawEvtData
+//  - MPedestalCam
+//
+// The following output containers are also searched and created if
+// they were not found:
+//  - MCerPhotEvt
+//
+Int_t MCerPhotAnal2::PreProcess(MParList *pList)
+{
+    fRunHeader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
+    if (!fRunHeader)
+    {
+        *fLog << dbginf << "MRawRunHeader not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
+    if (!fRawEvt)
+    {
+        *fLog << dbginf << "MRawEvtData not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    fCerPhotEvt = (MCerPhotEvt*)pList->FindCreateObj("MCerPhotEvt");
+    if (!fCerPhotEvt)
+        return kFALSE;
+
+    const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
+    if (!runheader)
+        *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
+    else
+        if (runheader->GetRunType() == kRTMonteCarlo)
+        {
+            fPedestals=NULL;
+            return kTRUE;
+        }
+
+    fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam");
+    if (runheader && !fPedestals)
+        return kFALSE;
+
+    return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// Calculate the integral of the FADC time slices and store them as a new
+// pixel in the MCerPhotEvt container.
+//
+Int_t MCerPhotAnal2::Process()
+{
+    MRawEvtPixelIter pixel(fRawEvt);
+
+    while (pixel.Next())
+    {
+        Byte_t max = pixel.GetNumMaxHiGainSample();
+	Byte_t num = fRawEvt->GetNumHiGainSamples();
+
+        Byte_t *ptr = pixel.GetHiGainSamples();
+
+        ULong_t sumb  = 0;   // sum background
+        ULong_t sqb   = 0;   // sum sqares background
+        ULong_t sumsb = 0;   // sum signal+background
+        ULong_t sqsb  = 0;   // sum sqares signal+background
+
+        Int_t sat = 0; // saturates?
+        Int_t nb  = 0;
+        Int_t nsb = 0;
+        for (int i=0; i<num; i++)
+        {
+            if (ptr[i]>250)
+                sat++;
+            if (sat>1)
+                continue;
+
+            if (i<max-fBefore || i>max+fAfter)
+            {
+                sumb += ptr[i];
+                sqb  += ptr[i]*ptr[i];
+                nb++;
+            }
+            else
+            {
+                sumsb += ptr[i];
+                sqsb  += ptr[i]*ptr[i];
+                nsb++;
+            }
+        }
+
+        if (sat>1)
+        {
+            // Area: x9
+            max = pixel.GetNumMaxLoGainSample();
+            num = fRawEvt->GetNumLoGainSamples();
+
+            ptr = pixel.GetLoGainSamples();
+
+            sumsb = 0;   // sum signal+background
+            sqsb  = 0;   // sum sqares signal+background
+            sumb  = 0;   // sum background
+            sqb   = 0;   // sum sqares background
+
+            nb = 0;
+            nsb = 0;
+            for (int i=0; i<num; i++)
+            {
+                if (ptr[i]>250)
+                    return kCONTINUE;
+                if (i<max-fBefore || i>max+fAfter)
+                {
+                    sumb += ptr[i];
+                    sqb  += ptr[i]*ptr[i];
+                    nb++;
+                }
+                else
+                {
+                    sumsb += ptr[i];
+                    sqsb  += ptr[i]*ptr[i];
+                    nsb++;
+                }
+            }
+        }
+
+        Float_t b  = (float)sumb/nb;       // background
+        Float_t sb = (float)sumsb/nsb;     // signal+background
+
+        Float_t msb  = (float)sqb/nb;      // mean square background
+        //Float_t mssb = (float)sqsb/nsb;    // mean square signal+background
+
+        Float_t sigb  = sqrt(msb-b*b);     // sigma background
+        //Float_t sigsb = sqrt(mssb-sb*sb);  // sigma signal+background
+
+        Float_t s   = sb-b;                // signal
+        Float_t sqs = sqsb-nsb*b;          // sum sqaures signal
+
+        Float_t mss  = (float)sqs/nsb;     // mean quare signal
+        Float_t sigs = sqrt(mss-s*s);      // sigma signal
+
+        if (sat>1)
+            s*=9;
+
+        Int_t idx = pixel.GetPixelId();
+        fCerPhotEvt->AddPixel(idx, s, sigs);
+
+        // Preliminary: Do not overwrite pedestals calculated by
+        // MMcPedestalCopy and MMcPedestalNSBAdd
+        if (fPedestals)
+            (*fPedestals)[idx].Set(b, sigb);
+    }
+
+    fCerPhotEvt->FixSize();
+    fCerPhotEvt->SetReadyToSave();
+
+    if (fPedestals)
+        fPedestals->SetReadyToSave();
+
+    return kTRUE;
+}
+