source: trunk/MagicSoft/Mars/msimcamera/MSimTrigger.cc@ 9272

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

//////////////////////////////////////////////////////////////////////////////
//
//  MSimTrigger
//
// This task takes the pure analog channels and simulates a trigger
// electronics.
//
// In a first step several channels can be summed together by a look-up table
// fRouteAC.
//
// In a second step from these analog channels the output of a discriminator
// is calculated using a threshold and optional a fixed digital signal length.
//
// With a second look-up table fCoincidenceMap the analog channels are
// checked for coincidences. The earliest coincide is then stored as
// trigger position.
//
//
//  Input Containers:
//   IntendedPulsePos [MParameterD]
//   MAnalogChannels
//   MRawRunHeader
//
//  Output Containers:
//   TriggerPos [MParameterD]
//   MRawEvtHeader
//
//////////////////////////////////////////////////////////////////////////////
#include "MSimTrigger.h"

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

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

#include "MLut.h"
#include "MArrayI.h"

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

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

ClassImp(MSimTrigger);

using namespace std;

// --------------------------------------------------------------------------
//
//  Default Constructor.
//
MSimTrigger::MSimTrigger(const char *name, const char *title)
    : fCamera(0), fPulsePos(0), fTrigger(0), fRunHeader(0), fEvtHeader(0),
    fDiscriminatorThreshold(-1)
{
    fName  = name  ? name  : "MSimTrigger";
    fTitle = title ? title : "Task to simulate trigger electronics";
}

// --------------------------------------------------------------------------
//
// Take two TObjArrays with a collection of digital signals.
// Every signal from one array is compared with any from the other array.
// For all signals whihc overlaps and which have an overlap time >gate
// a new digital signal is created storing start time and length of overlap.
// They are collected in a newly allocated TObjArray. A pointer to this array
// is returned.
//
// Th euser gains owenership of the object, ie.e., the user is responsible of
// deleting the memory.
//
TObjArray *MSimTrigger::CalcCoincidence(const TObjArray &arr1, const TObjArray &arr2, Float_t gate) const
{
    TObjArray *res = new TObjArray;

    if (arr1.GetEntriesFast()==0 || arr2.GetEntriesFast()==0)
        return res;

    TIter Next1(&arr1);
    MDigitalSignal *ttl1 = 0;
    while ((ttl1=static_cast<MDigitalSignal*>(Next1())))
    {
        TIter Next2(&arr2);
        MDigitalSignal *ttl2 = 0;
        while ((ttl2=static_cast<MDigitalSignal*>(Next2())))
        {
            MDigitalSignal *ttl = new MDigitalSignal(*ttl1, *ttl2);

            if (ttl->GetLength()<=gate)
            {
                delete ttl;
                continue;
            }

            res->Add(ttl);
        }
    }

    res->SetOwner();

    return res;
}

// --------------------------------------------------------------------------
//
// Check for the necessary parameter containers. Read the luts.
//
Int_t MSimTrigger::PreProcess(MParList *pList)
{
    fTrigger = (MParameterD*)pList->FindCreateObj("MParameterD", "TriggerPos");
    if (!fTrigger)
        return kFALSE;

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

    fCamera = (MAnalogChannels*)pList->FindObject("MAnalogChannels");
    if (!fCamera)
    {
        *fLog << err << "MAnalogChannels 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;


    if (fRouteAC.ReadFile(fNameRouteAC)<0)
        return kFALSE;

    if (fCoincidenceMap.ReadFile(fNameCoincidenceMap)<0)
        return kFALSE;

    if (fDiscriminatorThreshold<=0)
    {
        *fLog << err << "ERROR - Discriminator threshold " << fDiscriminatorThreshold << " invalid." << endl;
        return kFALSE;
    }

    *fLog << inf << "Using discriminator threshold of " << fDiscriminatorThreshold << endl;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
Int_t MSimTrigger::Process()
{
    // Invalidate trigger
    fTrigger->SetVal(-1);

    // FIXME: Set from somewhere else? (see also MSimCalibrationSignal)
    fEvtHeader->SetDAQEvtNumber(GetNumExecutions());
    fEvtHeader->SetReadyToSave();

    // ================== Simulate channel bundling ====================

    const UInt_t npatch = fRouteAC.GetEntriesFast();

    MAnalogChannels patches(npatch, fCamera->GetNumSamples());

    for (UInt_t i=0; i<npatch; i++)
    {
        const MArrayI &row = fRouteAC.GetRow(i);
        for (UInt_t j=0; j<row.GetSize(); j++)
        {
            const UInt_t idx = row[j];
            patches[i].AddSignal((*fCamera)[idx]);
        }
    }

    // FIXME: Write patches

    // ================== Simulate discriminators ====================

    TObjArray ttls(npatch);
    ttls.SetOwner();

        // MAGIC: Constant length of signals between 6ns and 12ns
    const Double_t siglen = 8; // Signal length (-1: As long as Signal-above-Threshold)

    for (UInt_t i=0; i<npatch; i++)
        ttls.AddAt(patches[i].Discriminate(fDiscriminatorThreshold, siglen), i);

    // FIXME: Write TTLs!

    // =================== Simulate coincidences ======================

    // Calculate the minimum and maximum time for a valid trigger
    const Double_t freq    = fRunHeader->GetFreqSampling()/1000.;
    const Float_t  nsamp   = fRunHeader->GetNumSamplesHiGain();
    const Float_t  pulspos = fPulsePos->GetVal()/freq;

    const Float_t min = fCamera->GetValidRangeMin()+pulspos;
    const Float_t max = fCamera->GetValidRangeMax()+pulspos-nsamp;

    // Define gate time (minimum coincidence time)
    const Double_t gate = 1; // MAGIC: minimum coincidence time 0.25ns to 1ns,

    // Create an array for the individual triggers
    TObjArray triggers;
    triggers.SetOwner();

    for (int j=0; j<fCoincidenceMap.GetEntries(); j++)
    {
        const MArrayI &idx = fCoincidenceMap.GetRow(j);

        // Start with a copy of the first coincidence channel
        TObjArray *arr = static_cast<TObjArray*>(ttls[idx[0]]->Clone());
        arr->SetOwner();

        // compare to all other channels in this coincidence patch, one by one
        for (UInt_t k=1; k<idx.GetSize() && arr->GetEntriesFast()>0; k++)
        {
            TObjArray *res = CalcCoincidence(*arr, *static_cast<TObjArray*>(ttls[idx[k]]),
                                             gate);

            // Delete the original array and keep the new one
            delete arr;
            arr = res;
        }

        // Remove all signals which are not in the valid digitization range
        // (This is not the digitization window, but the region in which
        //  the analog channels contain usefull data)
        TIter Next(arr);
        MDigitalSignal *ttl = 0;
        while ((ttl=static_cast<MDigitalSignal*>(Next())))
            if (ttl->GetStart()<min || ttl->GetStart()>max)
                delete arr->Remove(ttl);

        // Remove the empty slots
        arr->Compress();

        // If we have at least one trigger keep the earliest one.
        // FIXME: The triggers should be ordered in time automatically: To be checked!
        // FIXME: Simulate trigger dead-time!
        if (arr->GetEntriesFast()>0)
            triggers.Add(arr->RemoveAt(0));

        // delete the allocated space
        delete arr;
    }

    // No trigger issued. Go on.
    if (triggers.GetEntriesFast()==0)
        return kTRUE;

    // There are usually not enough entries that it is worth to search
    // for the earliest instead of just sorting and taking the first one
    //  FIXME: This could be improved if checking for IsSortable
    //         is omitted
    triggers.Sort();

    // inf2?
    *fLog << all;
    *fLog << triggers.GetEntriesFast() << " triggers in " << triggers.GetEntriesFast() << " patches.";
    *fLog << endl;

    // FIXME: Jitter! (Own class?)
    fTrigger->SetVal(static_cast<MDigitalSignal*>(triggers[0])->GetStart());


    return kTRUE;
}

// --------------------------------------------------------------------------
//
// FileNameRouteac:         routeac.txt
// FileNameCoincidenceMap:  coincidence.txt
// DiscriminatorTheshold:   3.5
//
Int_t MSimTrigger::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc = kFALSE;
    if (IsEnvDefined(env, prefix, "FileNameRouteAC", print))
    {
        rc = kTRUE;
        fNameRouteAC = GetEnvValue(env, prefix, "FileNameRouteAC", fNameRouteAC);
    }

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

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

    return rc;
}