source: trunk/MagicSoft/Mars/msimcamera/MSimAPD.cc@ 9328

/* ======================================================================== *\
!
! *
! * This file is part of CheObs, the Modular 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 appears 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,  1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: CheObs Software Development, 2000-2009
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  MSimAPD
//
// This tasks simulates the individual APDs. Before starting the APD is
// initialized randomly according to the photon rate hitting the APD. Such
// it is assumed that the initial condition of the APD is similar to the real
// one. In this context it is assumed that the events are independent, so
// that the APD is always in the same condition.
//
// For every photon and event the behaviour of the APD is simulated. The
// output is set as weight to the MPhotonData containers.
//
// Remark:
//   - The photons MUST be sorted increasing in time.
//   - The photon rate used to initialize the APD must match the one used
//     to "fill" the random photons. (FIXME: This should be stored somewhere)
//
//  Input Containers:
//   fNameGeomCam [MGeomCam]
//   MPhotonEvent
//   MPhotonStatistics
//
//  Output Containers:
//   MPhotonEvent
//
//////////////////////////////////////////////////////////////////////////////
#include "MSimAPD.h"

#include <TH2.h>
#include <TRandom.h>

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

#include "MMath.h"
#include "MParList.h"

#include "MGeomCam.h"

#include "MPhotonEvent.h"
#include "MPhotonData.h"

#include "MAvalanchePhotoDiode.h"

ClassImp(MSimAPD);

using namespace std;

// --------------------------------------------------------------------------
//
//  Default Constructor.
//
MSimAPD::MSimAPD(const char* name, const char *title)
: fGeom(0), fEvt(0), fStat(0)
{
    fName  = name  ? name  : "MSimAPD";
    fTitle = title ? title : " Task to simulate the detection behaviour of APDs";
}

// --------------------------------------------------------------------------
//
//  Get the necessary parameter containers
//
Int_t MSimAPD::PreProcess(MParList *pList)
{
    if (fNameGeomCam.IsNull())
    {
        *fLog << inf << "No geometry container... skipping." << endl;
        return kSKIP;
    }

    fGeom = (MGeomCam*)pList->FindObject(fNameGeomCam, "MGeomCam");
    if (!fGeom)
    {
        *fLog << inf << fNameGeomCam << " [MGeomCam] not found..." << endl;

        fGeom = (MGeomCam*)pList->FindCreateObj(fNameGeomCam);
        if (!fGeom)
            return kFALSE;
    }

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

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

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Initialize as many APDs as we have pixels in the fGeomCam
//
Bool_t MSimAPD::ReInit(MParList *plist)
{
    if (UInt_t(fAPDs.GetEntriesFast())!=fGeom->GetNumPixels())
    {
        fAPDs.Delete();

        // FIXME:
        //   * initialize an empty APD and read the APD setup from a file to
        //     allow different APDs.
        //   * Make the arguments a data member of MSimAPD

        for (UInt_t i=0; i<fGeom->GetNumPixels(); i++)
            //fAPDs.Add(new APD(60, 0.2, 3, 8.75));
            //fAPDs.Add(new APD(60/2, 0.2, 3e-9, 8.75e-9*4));
            fAPDs.Add(new APD(60/2, 0.2, 3, 8.75*4));
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Process all photons through the corresponding APD and set the output
// (weight) accordingly.
//
Int_t MSimAPD::Process()
{
    //const Double_t rate = 40e9;
    // FIXME: Where do we get this number from??
    // const Double_t rate = 0.04;

    // Make all APDs look neutral for the first hit by a photon according to the
    // average hit rate
    const UInt_t npix = fAPDs.GetEntriesFast();

    // Check if we can safely proceed (this can fail if we either haven't been
    // ReInit'ed or the max index in MPhotonStatistics is wrong)
    if ((Int_t)npix<fStat->GetMaxIndex())
    {
        *fLog << err << "ERROR - MSimAPD::Process: Only " << npix << " APDs initialized. At least " << fStat->GetMaxIndex() << " needed... abort." << endl;
        return kERROR;
    }

    for (UInt_t idx=0; idx<npix; idx++)
        static_cast<APD*>(fAPDs.UncheckedAt(idx))->FillRandom(fFreq, fStat->GetTimeFirst());

    // Get number of photons
    const Int_t num = fEvt->GetNumPhotons();

    // Loop over all photons
    for (Int_t i=0; i<num; i++)
    {
        // Get i-th photon
        MPhotonData &ph = (*fEvt)[i];

        // Get arrival time of photon and idx
        const Double_t t    = ph.GetTime();
        const UInt_t   idx  = ph.GetTag();

        if (ph.GetWeight()!=1)
        {
            *fLog << err << "ERROR - MSimAPD: Weight of " << i << "-th photon not 1, but " << ph.GetWeight() << endl;
            ph.Print();
            return kERROR;
        }

        // Simulate hitting the APD (the signal height in effective
        // "number of photons" is returned)
        const Double_t hits = static_cast<APD*>(fAPDs.UncheckedAt(idx))->HitRandomCell(t);

        // FIXME: Make a proper simulation of the excess noise!!!
        //const Double_t signal = gRandom->Gaus(hits, 0.2*TMath::Sqrt(hits));

        // Set the weight to the input
        ph.SetWeight(hits);
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// NameGeomCam
//
Int_t MSimAPD::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc = kFALSE;
    if (IsEnvDefined(env, prefix, "NameGeomCam", print))
    {
        rc = kTRUE;
        fNameGeomCam = GetEnvValue(env, prefix, "NameGeomCam", fNameGeomCam);
    }

    return rc;
}