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

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

#include "MPedestalCam.h"
#include "MPedestalPix.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), fElectronicNoise(0), fGain(0),
    fDiscriminatorThreshold(-1), fDigitalSignalLength(8), fCoincidenceTime(0.5),
    fShiftBaseline(kTRUE), fUngainSignal(kTRUE)
{
    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;
    }

    fElectronicNoise = 0;
    if (fShiftBaseline)
    {
        fElectronicNoise = (MPedestalCam*)pList->FindObject("ElectronicNoise", "MPedestalCam");
        if (!fElectronicNoise)
        {
            *fLog << err << "ElectronicNoise [MPedestalCam] not found... aborting." << endl;
            return kFALSE;
        }
        *fLog << inf << "Baseline will be shifted back to 0 for discriminator." << endl;
    }

    fGain = 0;
    if (fUngainSignal)
    {
        fGain = (MPedestalCam*)pList->FindObject("Gain", "MPedestalCam");
        if (!fGain)
        {
            *fLog << err << "Gain [MPedestalCam] not found... aborting." << endl;
            return kFALSE;
        }
        *fLog << inf << "Discriminator will be multiplied by applied gain." << endl;
    }

    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;

    fRouteAC.Delete();
    if (!fNameRouteAC.IsNull() && fRouteAC.ReadFile(fNameRouteAC)<0)
        return kFALSE;

    fCoincidenceMap.Delete();
    if (!fNameCoincidenceMap.IsNull() && fCoincidenceMap.ReadFile(fNameCoincidenceMap)<0)
        return kFALSE;

    // ---------------- Consistency checks ----------------------

    if (!fRouteAC.IsEmpty() && !fCoincidenceMap.IsEmpty() &&
        fCoincidenceMap.GetMaxIndex()>fRouteAC.GetNumRows()-1)
    {
        *fLog << err;
        *fLog << "ERROR - AC routing produces " << fRouteAC.GetNumRows() << " analog channels," << endl;
        *fLog << "        but the coincidence map expects at least " << fCoincidenceMap.GetMaxIndex()+1 << " channels." << endl;
        return kERROR;
    }

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

    if (fElectronicNoise && !fRouteAC.IsEmpty() && !fRouteAC.IsDefault())
    {
        // FIXME: Apply to analog channels when summing
        *fLog << warn << "WARNING - A baseline shift doesn't make sense for sum-channels... reset." << endl;
        fElectronicNoise = 0;
    }

    if (fGain && !fRouteAC.IsEmpty() && !fRouteAC.IsDefault())
    {
        // FIXME: Apply to analog channels when summing
        *fLog << warn << "WARNING - Ungain doesn't make sense for sum-channels... reset." << endl;
        fGain = 0;
    }


    // ---------------- Information output ----------------------

    *fLog << inf;

    if (fRouteAC.IsEmpty())
        *fLog << "Re-routing/summing of analog channels before discriminator switched off." << endl;
    else
        *fLog << "Using " << fNameRouteAC << " for re-routing/summing of analog channels before discriminator." << endl;

    if (fCoincidenceMap.IsEmpty())
        *fLog << "No coincidences of digital channels will be checked. Signal-above-threshold trigger applied." << endl;
    else
        *fLog << "Using " << fNameCoincidenceMap << " to check for coincidences of the digital channels." << endl;

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

    return kTRUE;
}

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

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

    // FIXME: Before we can bundle the channels we have to make a copy
    //        and simulate clipping

    // Check if routing should be done
    const Bool_t empty = fRouteAC.IsEmpty();

    // If no channels are summed the number of patches stays the same
    const UInt_t npatch = empty ? fCamera->GetNumChannels() : fRouteAC.GetEntriesFast();

    // Use the given analog channels as default out. If channels are
    // summed overwrite with a newly allocated set of analog channels
    MAnalogChannels *patches = fCamera;
    if (!empty)
    {
        // FIXME: Can we add gain and offset here into a new container?

        patches = new MAnalogChannels(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();

    for (UInt_t i=0; i<npatch; i++)
    {
        // FIXME: What if the gain was also allpied to the baseline?
        const Double_t offset = fElectronicNoise ? (*fElectronicNoise)[i].GetPedestal() : 0;
        const Double_t gain   = fGain            ? (*fGain)[i].GetPedestal()            : 1;
        ttls.AddAt((*patches)[i].Discriminate(fDiscriminatorThreshold*gain+offset, fDigitalSignalLength), i);
    }

    // FIXME: Write TTLs!

    // If analog channels had been newly allocated free memmory
    if (patches!=fCamera)
        delete patches;

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

    // If the map is empty we create a one-pixel-coincidence map
    // FIMXE: This could maybe be accelerated if the Clone can be
    //        omitted in the loop
    if (fCoincidenceMap.IsEmpty())
        fCoincidenceMap.SetDefault(npatch);

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