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

/* ======================================================================== *\
!
! *
! * 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.
//
// The signal length of the digital signal emitted by the discriminator
// can either be bound to the time the signal is above the threshold
// defined by fDiscriminatorThreshold if fDigitalSignalLength<0 or set to a
// fixed length (fDigitalSignalLength>0).
//
// With a second look-up table fCoincidenceMap the analog channels are
// checked for coincidences. The coincidense must at least be of the length
// defined by fCoincidenceTime. The earliest coincide is then stored as
// trigger position.
//
//
// For MAGIC1:
//  - fDigitalSignalLength between 6ns and 12ns
//  - fCoincidenceTime between 0.25ns to 1ns
//
//
//  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"

#include "MTriggerPattern.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), fDigitalSignalLength(8), fCoincidenceTime(0.5)
{
    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);

    // ================== 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++)
        {
            // FIXME: Simulate clipping
            const UInt_t idx = row[j];
            patches[i].AddSignal((*fCamera)[idx]);
        }
    }

    // FIXME: Write patches

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

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

    for (UInt_t i=0; i<npatch; i++)
        ttls.AddAt(patches[i].Discriminate(fDiscriminatorThreshold, fDigitalSignalLength), 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;

    // Valid range in units of bins
    const Float_t min = fCamera->GetValidRangeMin()+pulspos;
    const Float_t max = fCamera->GetValidRangeMax()-(nsamp-pulspos);

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

    Int_t cnt  = 0;
    Int_t rmlo = 0;
    Int_t rmhi = 0;

    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]]),
                                             fCoincidenceTime);

            // 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)
            {
                delete arr->Remove(ttl);
                rmlo++;
            }
            if (ttl->GetStart()>max)
            {
                delete arr->Remove(ttl);
                rmhi++;
            }
        }

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

        cnt += arr->GetEntriesFast();

        // 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)
    {
        if (rmlo>0 || rmhi>0)
            *fLog << inf2 << rmlo << "/" << rmhi << " trigger out of valid range. No trigger raised." << endl;
        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();

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

    // Flag this event as triggered by the lvl1 trigger (FIXME?)
    fEvtHeader->SetTriggerPattern((UInt_t)~(MTriggerPattern::kTriggerLvl1 | (MTriggerPattern::kTriggerLvl1<<8)));
    fEvtHeader->SetReadyToSave();

    // inf2?
    *fLog << inf;
    *fLog << cnt << " triggers left in " << triggers.GetEntriesFast() << " patches (" << rmlo << "/" << rmhi << " trigger out of valid range), T=" << fTrigger->GetVal();
    *fLog << endl;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// FileNameRouteac:         routeac.txt
// FileNameCoincidenceMap:  coincidence.txt
// DiscriminatorTheshold:   3.5
// DigitalSignalLength:     8
// CoincidenceTime:         0.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);
    }

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

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

    return rc;
}