source: trunk/MagicSoft/Mars/manalysis/MPedPhotCalc.cc@ 2841

/* ======================================================================== *\
!
! *
! * 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): Josep Flix  <mailto:jflix@ifae.es>
!              Javier Rico <mailto:jrico@ifae.es> 
!              (01/04)
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//   MPedPhotCalc                                                          //
//                                                                         //  
//  This is a Task class to compute, for each pixel, the signal mean and   //
//  rms from a pedestal run file that has undergone the standard signal    //
//  extration  and calibration procedure. The signal rms can be used as    //
//  reference to compute the significance of the measured signals in the   //
//  following data runs (e.g. during the image cleaning).                  //
//                                                                         //
//  Input Containers:                                                      //
//   MCerPhotEvt                                                           //
//                                                                         //
//  Output Containers:                                                     //
//   MPedPhotCam 
//                                                                         //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////

#include "MLog.h"
#include "MLogManip.h"

#include "MPedPhotCalc.h"
#include "MParList.h"
#include "MRawRunHeader.h"  
#include "MCerPhotPix.h"
#include "MCerPhotEvt.h"
#include "MPedPhotPix.h"
#include "MPedPhotCam.h"

ClassImp(MPedPhotCalc);

using namespace std;

MPedPhotCalc::MPedPhotCalc(const char *name, const char *title)
{
  fName  = name  ? name  : "MPedPhotCalc";
  fTitle = title ? title : "Task to calculate pedestals from the charge computed from pedestal runs (in units of photons)";
}

Int_t MPedPhotCalc::PreProcess( MParList *pList )
{      
  // Look for input container
  fCerPhot = (MCerPhotEvt*)pList->FindObject("MCerPhotEvt");
  if (!fCerPhot)
    {
      *fLog << err << "MPedPhotCalc::PreProcess Error: MCerPhotEvt not found... aborting." << endl;
      return kFALSE;
    }

  // Create output container
  fPedestals = (MPedPhotCam*)pList->FindCreateObj("MPedPhotCam");
  if (!fPedestals)
    return kFALSE;
  
  return kTRUE;
}

Bool_t MPedPhotCalc::ReInit(MParList *pList)
{
  const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
  if (!runheader)
    {
      *fLog << warn << dbginf;
      *fLog << "Warning - cannot check file type, MRawRunHeader not found." << endl;
    }
  else
    if (runheader->GetRunType() == kRTMonteCarlo)
      return kTRUE;
  

  fNumPixels = fPedestals->GetSize();

  // Initialize arrays
  if(fSumx.GetSize()==0)
    {
      fSumx.Set(fNumPixels);
      fSumx2.Set(fNumPixels);
  
      memset(fSumx.GetArray(),  0, sizeof(Float_t)*fNumPixels);
      memset(fSumx2.GetArray(), 0, sizeof(Float_t)*fNumPixels);
    }

    return kTRUE;
}

Int_t MPedPhotCalc::Process()
{
  for(UInt_t i=0;i<fCerPhot->GetNumPixels();i++)
    {
      MCerPhotPix &pix = (*fCerPhot)[i];
      Float_t nphot = pix.GetNumPhotons();
      Int_t idx     = pix.GetPixId();
      
      fSumx[idx]  += nphot;
      fSumx2[idx] += nphot*nphot;
    }
  
  fPedestals->SetReadyToSave();
  
  return kTRUE;
}

Int_t MPedPhotCalc::PostProcess()
  {
    // Compute pedestals and rms from fSumx and fSumx2 arrays
    const Int_t n  = GetNumExecutions();

    for(Int_t i=0;i<fNumPixels;i++)
      {
        const Float_t sum  = fSumx.At(i);
        const Float_t sum2 = fSumx2.At(i);
        
        const Float_t photped = sum/n;
        const Float_t photrms = TMath::Sqrt((sum2-sum*sum/n)/(n-1.));

        (*fPedestals)[i].Set(photped,photrms);  
      }
    
    return kTRUE;
}