source: trunk/Mars/msimcamera/MSimReadout.cc@ 16485

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

//////////////////////////////////////////////////////////////////////////////
//
//  MSimReadout
//
// Task to convert the analog channels into a digital signal. This should
// simulate the conversion and saturation bahaviour of the FADC/readout
// system.
//
// You can give a conversion factor from the unitx of your analog signal
// to the units of your adc. This is a fixed factor because it is just
// a matter of what the meaning of an adc count is, nothing which could
// jitter or is a real part of the electronics. Such effects should
// be simulated somewhere else.
//
//
//  Input Containers:
//   MGeomCam
//   MAnalogChannels
//   TriggerPos [MParameterD]
//   IntendedPulsePos [MParameterD]
//   MRawRunHeader
//
//  Output Containers:
//   MRawEvtData
//   MRawEvtHeader
//
//////////////////////////////////////////////////////////////////////////////
#include "MSimReadout.h"

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

#include "MArrayI.h"

#include "MParList.h"
#include "MParameters.h"

#include "MGeomCam.h"

#include "MRawRunHeader.h"
#include "MRawEvtHeader.h"
#include "MRawEvtData.h"

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

ClassImp(MSimReadout);

using namespace std;


// ------------------------------------------------------------------------
//
// Default constructor
//
MSimReadout::MSimReadout(const char* name, const char *title)
    : fRunHeader(0), fEvtHeader(0), fCamera(0), fPulsePos(0), fTrigger(0), fData(0),
    fConversionFactor(1)
{
    fName  = name  ? name  : "MSimReadout";
    fTitle = title ? title : "Task to simulate the analog readout (FADCs)";
}

// ------------------------------------------------------------------------
//
// Look for the needed parameter containers.
// Initialize MRawEvtData from MRawEvtRunHeader.
//
Int_t MSimReadout::PreProcess(MParList *pList)
{
    fCamera = (MAnalogChannels*)pList->FindObject("MAnalogChannels");
    if (!fCamera)
    {
        *fLog << err << "MAnalogChannels not found... aborting." << endl;
        return kFALSE;
    }

    fTrigger = (MParameterD*)pList->FindObject("TriggerPos", "MParameterD");
    if (!fTrigger)
    {
        *fLog << err << "TriggerPos [MParameterD] not found... aborting." << endl;
        return kFALSE;
    }

    fPulsePos = (MParameterD*)pList->FindObject("IntendedPulsePos", "MParameterD");
    if (!fPulsePos)
    {
        *fLog << err << "IntendedPulsePos [MParameterD] not found... aborting." << endl;
        return kFALSE;
    }

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

    fEvtHeader = (MRawEvtHeader*)pList->FindCreateObj("MRawEvtHeader");
    if (!fEvtHeader)
        return kFALSE;

    fData = (MRawEvtData*)pList->FindCreateObj("MRawEvtData");
    if (!fData)
        return kFALSE;

    return kTRUE;
}

Bool_t MSimReadout::ReInit(MParList *plist)
{
    MGeomCam *cam = (MGeomCam*)plist->FindObject("MGeomCam");
    if (!cam)
    {
        *fLog << err << "MGeomCam not found... aborting." << endl;
        return kFALSE;
    }

    fRunHeader->InitPixels(cam->GetNumPixels());

    fData->InitRead(fRunHeader);
    fData->ResetPixels();
    fData->SetIndices();

    return kTRUE;
}

// ------------------------------------------------------------------------
//
// Convert (digitize) the analog channels into digital (FADC) data.
//
Int_t MSimReadout::Process()
{
    // Sanity checks
    if (fData->GetNumLoGainSamples()>0)
    {
        *fLog << err << "ERROR - MSimReadout: Lo-gains not implemented yet." << endl;
        return kERROR;
    }

    // Make sure that we have not more analog channels than pixels
    // FIXME: Is this really necessary?
    if (fCamera->GetNumChannels()>fData->GetNumPixels())
    {
        *fLog << err;
        *fLog << "ERROR - Number of analog channels " << fCamera->GetNumChannels();
        *fLog << " exceeds number of pixels " << fData->GetNumPixels() << endl;
        return kERROR;
    }

    if (fTrigger->GetVal()<0)
    {
        *fLog << err << "ERROR - MSimReadout: MSimReadout executed for an event which has no trigger." << endl;
        return kERROR;
    }

    // Get the intended pulse position and convert it to slices
    const Float_t pulpos = fPulsePos->GetVal()*fRunHeader->GetFreqSampling()/1000.;

    // Get trigger position and correct for intended pulse position
    const Int_t trig = TMath::CeilNint(fTrigger->GetVal()-pulpos);

    // Check if the position is valid
    if (trig<0)
    {
        *fLog << err;
        *fLog << "ERROR - Trigger position before analog signal." << endl;
        *fLog << "        Trigger:  " << fTrigger->GetVal() << endl;
        *fLog << "        PulsePos: " << pulpos << endl;
        return kERROR;
    }

    // Get Number of samples in analog channels
    const Int_t nsamp = fCamera->GetNumSamples();

    // Get number of samples to be digitized
    const Int_t nslices = fData->GetNumSamples();

    // Check if the whole requested signal can be digitized
    if (trig+nslices>=nsamp)
    {
        *fLog << err << "ERROR - Trigger position beyond valid analog signal range." << endl;
        *fLog << "        Trigger:    " << fTrigger->GetVal() << endl;
        *fLog << "        PulsePos:   " << pulpos << endl;
        *fLog << "        SamplesIn:  " << nsamp << endl;
        *fLog << "        SamplesOut: " << nslices << endl;
        return kERROR;
    }

    const Float_t offset    = 0;//128;
    const UInt_t  max       = fData->GetMax();

    // FIXME: Take this into account
//    const UInt_t scale      = 16;
//    const UInt_t resolution = 12;

    // Digitize into a buffer
    MArrayI buffer(nslices*fData->GetNumPixels());

    // Loop over all channels/pixels
    for (UInt_t i=0; i<fCamera->GetNumChannels(); i++)
    {
        // Get i-th canalog hannel
        const MAnalogSignal &sig = (*fCamera)[i];

        // Digitize all slices
        for (Int_t j=0; j<nslices; j++)
        {
            Float_t slice = j+trig>=(Int_t)sig.GetSize() ? offset :
                sig[j+trig] * fConversionFactor + offset;

            // FIXME: Handle/Implement saturation!
            if (slice>max)
                slice = max;
            if (slice<0)
                slice = 0;

            // We have a 16 bit FADC. The resolution of the DRS4 is just about 11.5 bit.
            // Therefore the last 4.5 bits (~22.5) are gaussian noise (is this right?)
            buffer[nslices*i + j] = TMath::Nint(slice);
        }
    }

    // Set samples as raw-data
    fData->Set(buffer);
    fData->SetReadyToSave();

    // Set the trigger/daq event number
    fEvtHeader->SetDAQEvtNumber(GetNumExecutions());
    fEvtHeader->SetReadyToSave();

    // FIMXE: This will never be stored correctly :(
    fRunHeader->SetNumEvents(fRunHeader->GetNumEvents()+1);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Read the parameters from the resource file.
//
//  ConversionFactor: 1
//
Int_t MSimReadout::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc = kFALSE;
    if (IsEnvDefined(env, prefix, "ConversionFactor", print))
    {
        rc = kTRUE;
        fConversionFactor = GetEnvValue(env, prefix, "ConversionFactor", fConversionFactor);
    }

    return rc;
}