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

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

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

#include "MSpline3.h"
#include "MPulseShape.h"

#include "MParList.h"
//#include "MParameters.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)//, fPulsePos(0)
{
    fName  = name  ? name  : "MSimCamera";
    fTitle = title ? title : "Task to simulate the electronic noise and to convert photons into pulses";
}

// --------------------------------------------------------------------------
//
// 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;
    }
/*
    fPulsePos = (MParameterD*)pList->FindObject("IntendedPulsePos", "MParameterD");
    if (!fPulsePos)
    {
        *fLog << err << "IntendedPulsePos [MParameterD] not found... aborting." << endl;
        return kFALSE;
    }
 */
    MPulseShape *pulse = (MPulseShape*)pList->FindObject("MPulseShape");
    if (!pulse)
    {
        *fLog << err << "MPulseShape not found... aborting." << endl;
        return kFALSE;
    }

    *fLog << warn << "FIXME - SCALE WITH THE SAMPLING FREQUENCY." << endl;

    fSpline = pulse->GetSpline();
    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;
    const Double_t end   = fStat->GetTimeLast() *freq;

    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;

    const Double_t pl = fSpline->GetXmin()*freq;
    const Double_t pr = fSpline->GetXmax()*freq;

    // Init the arrays and set the range which will contain valid data
    fCamera->Init(npix, nlen);
    fCamera->SetValidRange(TMath::FloorNint(pr), TMath::CeilNint(nlen+pl));

    // 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();

    // A random shift, uniformely distributed within one slice, to make sure that
    // the first photon is not always aligned identically with a sample edge.
    // FIXME: Make it switchable
    const Float_t rndm = gRandom->Uniform();

    // FIXME: Shell we add a random shift of [0,1] samples per channel?
    //        Or maybe per channel and run?

    // 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+rndm;// - fSpline->GetXmin();

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

        // === FIXME === FIXME === FIXME === Frequency!!!!
        (*fCamera)[idx].AddPulse(*fSpline, t, ph.GetWeight());
    }

    return kTRUE;
}

/*
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;
}
*/