source: trunk/MagicSoft/Mars/msimcamera/MSimCamera.cc@ 9256

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

//////////////////////////////////////////////////////////////////////////////
//
//  MSimCamera
//
//  This task initializes the analog channels with analog noise and simulated
//  the analog pulses from the photon signal.
//
//  Input Containers:
//   MPhotonEvent
//   MPhotonStatistics
//   MRawRunHeader
//
//  Output Containers:
//   MAnalogChannels
//
//////////////////////////////////////////////////////////////////////////////
#include "MSimCamera.h"

#include <TF1.h>

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

#include "MSpline3.h"

#include "MParList.h"

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

#include "MAnalogSignal.h"
#include "MAnalogChannels.h"

#include "MRawRunHeader.h"

ClassImp(MSimCamera);

using namespace std;

// --------------------------------------------------------------------------
//
//  Default Constructor.
//
MSimCamera::MSimCamera(const char* name, const char *title)
: fEvt(0), fStat(0), fRunHeader(0), fCamera(0), fSpline(0),
  fFunction("exp(-(x/2)^2/2)"), fNpx(25), fXmin(-25), fXmax(25)
{
    fName  = name  ? name  : "MSimCamera";
    fTitle = title ? title : "Task to simulate the electronic noise and to convert photons into pulses";
}

// --------------------------------------------------------------------------
//
//  Call Clear()
//
MSimCamera::~MSimCamera()
{
    Clear();
}

// --------------------------------------------------------------------------
//
//  Delete fSpline if set and set it to 0
//
void MSimCamera::Clear(Option_t *)
{
    if (fSpline)
        delete fSpline;
    fSpline=0;
}

// --------------------------------------------------------------------------
//
//  Read the intended pulse shape from a file and initialize the spline
// accordingly
//
Bool_t MSimCamera::ReadFile(const char *fname)
{
    if (!fRunHeader)
        return kFALSE;

    if (fname)
        fFileName = fname;

    *fLog << inf << "Reading pulse shape from " << fFileName << endl;

    const TGraph g(fFileName);
    if (g.GetN()==0)
    {
        *fLog << err << "ERROR - No data points from " << fFileName << "." << endl;
        return kFALSE;
    }

    // option: b1/e1 b2/e2   (first second derivative?)
    // option: valbeg/valend (first second derivative?)

    Clear();
    fSpline = new MSpline3(g, fRunHeader->GetFreqSampling()/1000.);

    return kTRUE;
}

void MSimCamera::SetFunction(const TF1 &f)
{
    // FIXME: Use TF1 directly? (In most cases this seems to be slower)
    if (!fRunHeader)
        return;// kFALSE;

    // option: b1/e1 b2/e2   (first second derivative?)
    // option: valbeg/valend (first second derivative?)

    // if (f.GetNpar()==0)
    // No SUPPORT

    Clear();
    fSpline = new MSpline3(f, fRunHeader->GetFreqSampling()/1000.);

    fFunction = f.GetTitle();
}

Bool_t MSimCamera::SetFunction(const char *func, Int_t n, Double_t xmin, Double_t xmax)
{
    // FIXME: Use TF1 directly? (In most cases this seems to be slower)
    TF1 f("f", func, xmin, xmax);
    f.SetNpx(n);

    SetFunction(f);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// SetReadyToSave for fRunHeader
//
Bool_t MSimCamera::ReInit(MParList *pList)
{
    // make that the run-header gets written to the file
    fRunHeader->SetReadyToSave();

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Search for the necessayr parameter containers.
// Setup spline for pulse shape.
//
Int_t MSimCamera::PreProcess(MParList *pList)
{
    fCamera = (MAnalogChannels*)pList->FindCreateObj("MAnalogChannels");
    if (!fCamera)
        return kFALSE;

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

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

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


    // FIMXE: Move to ReInit in case fRunHeader is read form file?
    if (!fFileName.IsNull())
        return ReadFile(fFileName);

    if (!fFunction.IsNull())
        return SetFunction(fFunction, fNpx, fXmin, fXmax);

    if (!fSpline)
    {
        *fLog << err << "No spline initialized." << endl;
        return kFALSE;
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// fStat->GetMaxIndex must return the maximum index possible
// (equiv. number of pixels) not just the maximum index stored!
//
Int_t MSimCamera::Process()
{
    // Calculate start time, end time and corresponding number of samples
    const Double_t freq = fRunHeader->GetFreqSampling()/1000.;

    const Double_t start = fStat->GetTimeFirst()*freq + fSpline->GetXmin();
    const Double_t end   = fStat->GetTimeLast() *freq + fSpline->GetXmax();

    const UInt_t   nlen  = TMath::CeilNint(end-start);

    // FIXME: Jitter the start point of digitization by one sample [0;1]
    // FIXME: TAKE PULSE WIDTH out of the calculation and start at TimeFirst

    // Get number of pixels/channels
    const UInt_t npix = fStat->GetMaxIndex()+1;

    // Init the arrays and set the range which will contain valid data
    fCamera->Init(npix, nlen);
    fCamera->SetValidRange(0   -TMath::FloorNint(fSpline->GetXmin()),
                           nlen-TMath::CeilNint( fSpline->GetXmax()));

    // Add electronic noise to empty channels
    for (UInt_t i=0; i<npix; i++)
    {
        // FIXME: We might add the base line here already.
        (*fCamera)[i].AddGaussianNoise(22.5/64);
    }

    // FIXME: Simulate correlations with neighboring pixels

    const Int_t num = fEvt->GetNumPhotons();

    // Simulate pulses
    for (Int_t i=0; i<num; i++)
    {
        const MPhotonData &ph = (*fEvt)[i];

        const UInt_t   idx = ph.GetTag();
        const Double_t t   = (ph.GetTime()-fStat->GetTimeFirst())*freq - fSpline->GetXmin();

        // FIXME: Time jitter?
        // FIXME: Add additional routing here?

        (*fCamera)[idx].AddPulse(*fSpline, t, ph.GetWeight());
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// FileName: pulse-shape.txt
// Function.Name: gaus
// Function.Npx:    50
// Function.Xmin:  -5
// Function.Xmax:   5
//
Int_t MSimCamera::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc = kFALSE;
    if (IsEnvDefined(env, prefix, "FileName", print))
    {
        rc = kTRUE;
        SetFileName(GetEnvValue(env, prefix, "FileName", fFileName));
    }

    if (IsEnvDefined(env, prefix, "Function.Name", print))
    {
        rc = kTRUE;
        SetFunction(GetEnvValue(env, prefix, "Function.Name", fFunction));

        if (IsEnvDefined(env, prefix, "Function.Npx", print))
            fNpx = GetEnvValue(env, prefix, "Function.Npx", fNpx);
        if (IsEnvDefined(env, prefix, "Function.Xmin", print))
            fXmin = GetEnvValue(env, prefix, "Function.Xmin", fXmin);
        if (IsEnvDefined(env, prefix, "Function.Xmax", print))
            fXmax = GetEnvValue(env, prefix, "Function.Xmax", fXmax);
    }

    return rc;
}

/*
MSimCameraFiles::Process()
{
    // fCorsikaHeader->GetEvtNumber()   -->   fMcEvt->SetEvtNumber()
    // fCorsikaHeader->GetTotalEnergy() -->   fMcEvt->SetEnergy()
    // fCorsikaHeader->GetParticleID()  -->   fMcEvt->SetParticleID()
    // fCorsikaHeader->GetImpact()      -->   fMcEvt->SetImpact()
    // fCorsikaHeader->GetTheta()       -->   fMcEvt->SetTheta()
    // fCorsikaHeader->GetPhi()         -->   fMcEvt->SetPhi()
    // fPointingPos->GetTheta()         -->   fMcEvt->SetTelescopeTheta()
    // fPointingPos->GetPhi()           -->   fMcEvt->SetTelescopePhi()
    // fStats->GetTimeFirst()           -->   fMcEvt->SetTimeFirst()
    // fStats->GetTimeLast()            -->   fMcEvt->SetTimeLast()
    //                                        fMcEvt->SetReuse()
    // MMcCorsikaRunHeader: fSlopeSpec, fELowLim, fEUppLim;

    fMcEvt->Fill(*fCorsikaHeader, *fPointingPos, *fStats);

    return kTRUE;
}
*/