source: trunk/MagicSoft/Mars/manalysis/MMcPedestalNSBAdd.cc@ 2132

/* ======================================================================== *\
!
! *
! * 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): Oscar Blanch  11/2001 < mailto:blanch@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
//  MMcPedestalNSBAdd                                                      //
//  -----------------                                                      //
//                                                                         //
//  This Task adds the contribution of the diffuse NSB to the FADC         //
//  pedestals. We assume that NSB introduces larger fluctuation but does   //
//  not change the mean value.                                             //
//  To be precise we add quadratically to the rms, which is already in     //
//  MPedestalCam, the NSB contribution.                                    //
//  The number of photons from the diffuse NSB follows a poisson           // 
//  distribution with expected mean value X, then its standard deviation   //
//  is sqrt(X).                                                            //
//  X is the number of phe per ns so we have to convert in FADC counts per // 
//  slice:                                                                 //
//                                                                         //
//  Y = sqrt(X * FADC_time / Number_of_slices * pixel_size)                //
//      * Amp_single_phe_response                                          //
//                                                                         //
//  Input Containers:                                                      //
//   MMcFadcHeader                                                         //
//   MMcRunHeader                                                          //
//   MRawRunHeader                                                         //
//                                                                         //
//  Output Containers:                                                     //
//   MPedestalCam                                                          //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MMcPedestalNSBAdd.h"

#include "MParList.h"

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

#include "MGeomCam.h"
#include "MGeomPix.h"

#include "MPedestalPix.h"
#include "MPedestalCam.h"

#include "MRawRunHeader.h"
#include "MMcRunHeader.hxx"
#include "MMcFadcHeader.hxx"

ClassImp(MMcPedestalNSBAdd);

// --------------------------------------------------------------------------
//
// Default constructor.
//
MMcPedestalNSBAdd::MMcPedestalNSBAdd(const Float_t difnsb,
                                     const char *name, const char *title)
    : fDnsbPixel(difnsb)
{
    fName  = name  ? name  : "MMcPedestalNSBAdd";
    fTitle = title ? title : "Add diffuse NSB to the pedestal signal";

    //
    // This is not needed here using MReadMarsFile because for the
    // RunHeader tree the auto scheme is disabled by default
    //
    AddToBranchList("MMcFadcHeader.fPedesMean");
    AddToBranchList("MMcFadcHeader.fElecNoise");
}

// --------------------------------------------------------------------------
//
// Check for the run type. Return kTRUE if it is a MC run or if there
// is no MC run header (old camera files) kFALSE in case of a different
// run type
//
Bool_t MMcPedestalNSBAdd::CheckRunType(MParList *pList) const
{
    const MRawRunHeader *runheader = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    if (!runheader)
    {
        *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
        return kTRUE;
    }

    return runheader->GetRunType() == kRTMonteCarlo;
}

// --------------------------------------------------------------------------
//
// - check whether we have a monte carlo file. if not skip this task
// - try to get MMcFadcHeader, MGeomCam and MPedestalCam from the parameter
//   list
// - try to find a MMcRunHeader, too
//
Bool_t MMcPedestalNSBAdd::PreProcess(MParList *pList)
{
    if (!CheckRunType(pList))
    {
        *fLog << warn << dbginf << "Warning - MMcPedestalNSB is for Monte Carlo files only... removing this task from list." << endl;
        return kSKIP;
    }

    fFadc = (MMcFadcHeader*)pList->FindObject("MMcFadcHeader");
    if (!fFadc)
    {
        *fLog << err << dbginf << "MMcFadcHeader not found... aborting." << endl;
        return kFALSE;
    }

    fGeom = (MGeomCam*)pList->FindObject("MGeomCam");
    if (!fGeom)
    {
        *fLog << err << dbginf << "MGeomCam not found... aborting." << endl;
        return kFALSE;
    }

    fPedCam = (MPedestalCam*)pList->FindCreateObj("MPedestalCam");
    if (!fPedCam)
        return kFALSE;

    const MMcRunHeader *mcrunheader = (MMcRunHeader*)pList->FindObject("MMcRunHeader");
    if (!mcrunheader && fDnsbPixel<0)
    {
        *fLog << err << dbginf << "Using the default argument only ";
        *fLog << "allowed if MMcRunHeader is available... aborting." << endl;
        return kFALSE;
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// If a MMcRunHeader is available the DNSB MMcRunHeader::GetNumPheFromDNSB
// is returned. Otherwise the user given number is used.
//
Float_t MMcPedestalNSBAdd::GetDnsb(MParList *pList) const
{
    const MMcRunHeader *mcrunheader = (MMcRunHeader*)pList->FindObject("MMcRunHeader");
    if (!mcrunheader && fDnsbPixel<0)
    {
        *fLog << err << dbginf << "Using the default argument only ";
        *fLog << "allowed if MMcRunHeader is available... aborting." << endl;
        return -1;
    }

    if (!mcrunheader)
        return fDnsbPixel;

    if (fDnsbPixel >= 0 && fDnsbPixel != mcrunheader->GetNumPheFromDNSB())
    {
        *fLog << warn << dbginf << "The MC file has been generated with diffuse nsb " << mcrunheader->GetNumPheFromDNSB();
        *fLog <<" but you set up the diffuse NSB to " << fDnsbPixel << endl;

        return fDnsbPixel;
    }

    return mcrunheader->GetNumPheFromDNSB();
}

// --------------------------------------------------------------------------
//
//  This function is called each time MReadTree::Notify is called, which 
//  happens when it  changes the file to read from.
//  Here we add the contribution from NSB to the pedestals.
//  The ReInit searches for the following input containers:
//   - MRawRunHeader
//   - MMcRunHeader
//   - MMcFacdHeader
//   - MGeomCam
//
//   The following output containers are also searched and created if
//   they were not found:
//   - MPedestalCam
//
Bool_t MMcPedestalNSBAdd::ReInit(MParList *pList)
{
    if (!CheckRunType(pList))
        return kFALSE;

    const Float_t dnsbpix = GetDnsb(pList) * 50.0/15.0;

    if (dnsbpix < 0)
        return kFALSE;

    const int num = fFadc->GetNumPixel();
    
    fPedCam->InitSize(num);

    const Float_t size0 = (*fGeom)[0].GetR() * (*fGeom)[0].GetR();

    for (int i=0; i<num; i++)
    {
        MPedestalPix   &pix     = (*fPedCam)[i];
        const MGeomPix &pixgeom = (*fGeom)[i];

        const Float_t pedrms = pix.GetSigma();
        const Float_t size   = pixgeom.GetR()*pixgeom.GetR()/size0;

        const Float_t ampl   = fFadc->GetAmplitud();

        pix.SetSigma(sqrt(pedrms*pedrms + dnsbpix*ampl*ampl*size));
    }

    return kTRUE;
}