source: trunk/MagicSoft/Mars/manalysis/MCerPhotAnal2.cc@ 2325

/* ======================================================================== *\
!
! *
! * 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)
{
    fSkip = 0;

    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]==255)
                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)
                {
                    fSkip++;
                    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*=10; // tgb has measured 9, but Florian said it's 10.

        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;
}

Int_t MCerPhotAnal2::PostProcess()
{
    if (GetNumExecutions()==0 || fSkip==0)
        return kTRUE;

    *fLog << inf << endl;
    *fLog << GetDescriptor() << " execution statistics:" << endl;
    *fLog << dec << setfill(' ');
    *fLog << " " << setw(7) << fSkip << " (" << setw(3) << (int)(fSkip*100/GetNumExecutions()) << "%) Evts skipped due to: lo gain saturated." << endl;

    return kTRUE;
}