source: trunk/Mars/msimcamera/MSimAPD.cc@ 10047

/* ======================================================================== *\
!
! *
! * 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 "MPedestalCam.h"
#include "MPedestalPix.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), fType(1)
{
    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->FindObject("MGeomCam");
        if (!fGeom)
        {
            *fLog << err << "MGeomCam not found... aborting." << endl;
            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;
    }

    fRates = (MPedestalCam*)pList->FindObject("AccidentalPhotonRates", "MPedestalCam");
    if (!fRates)
    {
        *fLog << inf;
        *fLog << "AccidentalPhotonRates [MPedestalCam] not found..." << endl;
        *fLog << " using " << fFreq << " as default for all G-APDs." << endl;
    }

    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())
        return kTRUE;

    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

    Int_t   ncells    = 0;
    Float_t crosstalk = 0;
    Float_t deadtime  = 0;
    Float_t recovery  = 0;

    switch (fType)
    {
    case 1:
        ncells    = 30;
        crosstalk = 0.2;
        deadtime  = 3;
        recovery  = 8.75*4;
        break;

    case 2:
        ncells    = 60;
        crosstalk = 0.2;
        deadtime  = 3;
        recovery  = 8.75;
        break;

    case 3:
        ncells    = 60;
        crosstalk = 0.15;
        deadtime  = 3;
        recovery  = 8.75;
        break;

    default:
        *fLog << err << "ERROR - APD type " << fType << " undefined." << endl;
        return kFALSE;
    }

    for (UInt_t i=0; i<fGeom->GetNumPixels(); i++)
        fAPDs.Add(new APD(ncells, crosstalk, deadtime, recovery));

    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++)
    {
        const Double_t freq = fRates ? (*fRates)[idx].GetPedestal() : fFreq;
        static_cast<APD*>(fAPDs.UncheckedAt(idx))->FillRandom(freq, fStat->GetTimeFirst());
    }
*/

    // This tries to initialize dead and relaxing cells properly. If
    // the APD has not been initialized before the chip is randomsly
    // filled, otherwise a time window of the default relaxing time
    // is simulated, so that the previous influence is less than a permille.
    for (UInt_t idx=0; idx<npix; idx++)
    {
        const Double_t freq = fRates ? (*fRates)[idx].GetPedestal() : fFreq;
        static_cast<APD*>(fAPDs.UncheckedAt(idx))->Init(freq);
    }

    // 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()-fStat->GetTimeFirst();
        const Int_t  idx = ph.GetTag();
        if (idx<0)
        {
            *fLog << err << "ERROR - MSimAPD: Invalid index -1." << endl;
            return kERROR;
        }

        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))->HitRandomCellRelative(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);
    }

    // Now we have to shift the evolved time of all APDs to the end of our
    // simulated time.
    for (UInt_t idx=0; idx<npix; idx++)
        static_cast<APD*>(fAPDs.UncheckedAt(idx))->IncreaseTime(fStat->GetTimeLast());

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// NameGeomCam
// Type: 1
//
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);
    }

    if (IsEnvDefined(env, prefix, "Type", print))
    {
        rc = kTRUE;
        fType = GetEnvValue(env, prefix, "Type", fType);
    }

    return rc;
}