source: trunk/MagicSoft/Mars/manalysis/MCerPhotAnal.cc@ 2893

/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//   MCerPhotAnal                                                           //
//                                                                          //
//   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 "MCerPhotAnal.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(MCerPhotAnal);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor.
//
MCerPhotAnal::MCerPhotAnal(const char *name, const char *title)
{
    fName  = name  ? name  : "MCerPhotAnal";
    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 MCerPhotAnal::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 MCerPhotAnal::Process()
{
    //fCerPhotEvt->InitSize(fRawEvt->GetNumPixels());

    MRawEvtPixelIter pixel(fRawEvt);

    while (pixel.Next())
    {
        Byte_t *ptr = pixel.GetHiGainSamples();
        const Byte_t *end = ptr + fRawEvt->GetNumHiGainSamples();
        const Byte_t *limit = end - 5;
        
        Int_t sum=0;
        Int_t sumpeak=0;
        Int_t sumlocal =0;
        Int_t slice=0;
        Int_t slicelocal=0;

        do
        {
            sumlocal = 0;
            for (Int_t i = 0; i<5;i++)
                sumlocal += *(ptr+i);

            if (sumpeak < sumlocal)
            {
                slice=slicelocal;
                sumpeak = sumlocal;
            }

            slicelocal++;
            sum += *ptr;
        } while (++ptr != limit);

        do sum += *ptr;
        while (++ptr != end);

        Float_t pedes = (Float_t)(sum-sumpeak)/(fRawEvt->GetNumHiGainSamples()-5);
        Float_t nphot = (Float_t)sumpeak - 5.0*pedes;

        Float_t sigmaped=0;

        slicelocal=0;
        sumlocal = 0;

        ptr = pixel.GetHiGainSamples();
        do
        {
            if (slicelocal==slice)
                ptr += 4;
            else
            {
                sumlocal = *ptr;

                const Float_t d = (Float_t)sumlocal-pedes;
                sigmaped += d*d;
            }
            slicelocal++;
        }
        while (++ptr != end);

        sigmaped /= (fRawEvt->GetNumHiGainSamples()-5);
        sigmaped = sqrt((float)sumlocal);

        const UInt_t pixid = pixel.GetPixelId();

        //
        // sanity check (old MC files sometimes have pixids>577)
        //
        if (fPedestals && !fPedestals->CheckBounds(pixid))
        {
            *fLog << inf << "Pixel ID larger than camera... skipping event." << endl;
            return kCONTINUE;
        }

        fCerPhotEvt->AddPixel(pixid, nphot, sigmaped/2.236);

        if (fPedestals)
            (*fPedestals)[pixid].Set(pedes, sigmaped);

        /*
         ped.SetPedestalRms(sigmaped/sqrt(fRawEvt->GetNumHiGainSamples()-5.),
                            sigmaped/sqrt(2.*(fRawEvt->GetNumHiGainSamples()-5)));
                            */
        // FIXME! Handling of Lo Gains is missing!
    }

    fCerPhotEvt->FixSize();
    fCerPhotEvt->SetReadyToSave();

    if (fPedestals)
        fPedestals->SetReadyToSave();

    return kTRUE;
}