source: trunk/MagicSoft/Mars/msignal/MExtractTimeAndCharge.cc@ 8154

/* ======================================================================== *\
! $Name: not supported by cvs2svn $:$Id: MExtractTimeAndCharge.cc,v 1.55 2006-10-24 08:24:52 tbretz Exp $
! --------------------------------------------------------------------------
!
! *
! * This file is part of MARS, the MAGIC 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 appear 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): Markus Gaug, 05/2004 <mailto:markus@ifae.es>
!   Author(s): Thomas Bretz, 05/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2006
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MExtractTimeAndCharge
//
//   Base class for the signal extractors which extract the arrival time 
//   and the signal at the same time. Uses the functions 
//   FindTimeAndChargeHiGain() and FindTimeAndChargeLoGain() to extract the signal and 
//   substract the pedestal value.
//
//   The following figure gives and example of possible inheritance trees. 
//   An extractor class can inherit from each of the following base classes:
//    - MExtractor
//    - MExtractTime
//    - MExtractTimeAndCharge
//
//Begin_Html
/*
<img src="images/ExtractorClasses.gif">
*/
//End_Html
//
//   The following variables have to be set by the derived class and 
//   do not have defaults:
//   - fNumHiGainSamples
//   - fNumLoGainSamples
//   - fSqrtHiGainSamples
//   - fSqrtLoGainSamples
//
// Input Containers:
//   MRawEvtData
//   MRawRunHeader
//   MPedestalCam
//
// Output Containers:
//   MArrivalTimeCam
//   MExtractedSignalCam
//
// For weired events check: Run 94127 Event 672, 1028
//
//////////////////////////////////////////////////////////////////////////////
#include "MExtractTimeAndCharge.h"

#include <TRandom.h>

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

#include "MParList.h"

//#include "MArrayB.h"
//#include "MArrayF.h"

//#include "MRawEvtData.h"
#include "MRawEvtPixelIter.h"
//#include "MRawRunHeader.h"

//#include "MPedestalCam.h"
//#include "MPedestalPix.h"
#include "MPedestalSubtractedEvt.h"

#include "MArrivalTimeCam.h"
#include "MArrivalTimePix.h"

#include "MExtractedSignalCam.h"
#include "MExtractedSignalPix.h"

ClassImp(MExtractTimeAndCharge);

using namespace std;

const Float_t MExtractTimeAndCharge::fgLoGainStartShift = -6.0;
const Byte_t  MExtractTimeAndCharge::fgLoGainSwitch     =  120;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Sets: 
// - fWindowSizeHiGain and fWindowSizeLoGain to 0
// - fLoGainStartShift to fgLoGainStartShift+fgOffsetLoGain
// - fLoGainSwitch     to fgLoGainSwitch
//
MExtractTimeAndCharge::MExtractTimeAndCharge(const char *name, const char *title)
    : /*fLoGainFirstSave(0),*/ fWindowSizeHiGain(0), fWindowSizeLoGain(0)
{
    fName  = name  ? name  : "MExtractTimeAndCharge";
    fTitle = title ? title : "Base class for signal and time extractors";

    SetLoGainStartShift();
    SetLoGainSwitch();
}

// --------------------------------------------------------------------------
//
// The PreProcess searches for the following input containers:
//  - MRawEvtData
//  - MRawRunHeader
//  - MPedestalCam
//
// The following output containers are also searched and created if
// they were not found:
//
//  - MExtractedSignalCam
//  - MArrivalTimeCam    
//
Int_t MExtractTimeAndCharge::PreProcess(MParList *pList)
{

  if (!MExtractTime::PreProcess(pList))
    return kFALSE;
  
  fSignals = (MExtractedSignalCam*)pList->FindCreateObj("MExtractedSignalCam",AddSerialNumber(fNameSignalCam));
  if (!fSignals)
      return kFALSE;

  *fLog << flush << inf;
  return kTRUE;
}

// --------------------------------------------------------------------------
//
// The ReInit calls:
// -  MExtractor::ReInit()
//
// Call: 
// - MArrivalTimeCam::SetUsedFADCSlices(fHiGainFirst, fHiGainLast, fNumHiGainSamples,
//                                      fLoGainFirst, fLoGainLast, fNumLoGainSamples);
// - InitArrays();
//
Bool_t MExtractTimeAndCharge::ReInit(MParList *pList)
{

  if (!MExtractTime::ReInit(pList))
    return kFALSE;

  if (!InitArrays())
    return kFALSE;
  
  if (fSignals)
    fSignals->SetUsedFADCSlices(fHiGainFirst, fHiGainLast/*+fHiLoLast*/, fNumHiGainSamples,
                                fLoGainFirst, fLoGainLast, fNumLoGainSamples);

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Calculate the integral of the FADC time slices and store them as a new
// pixel in the MArrivalTimeCam container.
// Calculate the integral of the FADC time slices and store them as a new
// pixel in the MExtractedSignalCam container. 
// The functions FindTimeAndChargeHiGain() and FindTimeAndChargeLoGain() are 
// supposed to extract the signal themselves.
//
Int_t MExtractTimeAndCharge::Process()
{
    MRawEvtPixelIter pixel(fRawEvt);

    while (pixel.Next())
    {
        const Int_t pixidx = pixel.GetPixelId();

        //
        // Preparations for the pedestal subtraction (with AB-noise correction)
        //
        Float_t *sig = fSignal->GetSamples(pixidx);

        // ====================================================================
        //    MOVE THIS TO MExtractTimeAndCharge
        // ====================================================================
        // number of hi- and lo-gains
        const Int_t numh = pixel.GetNumHiGainSamples();
        const Int_t numl = pixel.GetNumLoGainSamples();
/*
        // COPY HERE PRODUCING ARRAY WITH SAMPLES
        static MArrayB sample;
        sample.Set(numh+numl);

        if (numl>0)
        {
            memcpy(sample.GetArray(),      pixel.GetHiGainSamples(), numh);
            memcpy(sample.GetArray()+numh, pixel.GetLoGainSamples(), numl);
        }

        // get pedestal information
        const MPedestalPix &pedpix = (*fPedestals)[pixidx];

        // pedestal information
        const Int_t   ab  = pixel.HasABFlag() ? 1 : 0;
        const Float_t ped = pedpix.GetPedestal();

        // destination array
        static MArrayF sig;
        sig.Set(numh+numl);

        // start of destination array, end of hi-gain destination array
        // and start of hi-gain samples
        Float_t *beg = sig.GetArray();
        Float_t *end = beg + sig.GetSize();

        // determine with which pedestal (+/- AB offset) to start
        const Bool_t  noswap  = (ab&1)==((beg-beg)&1);
        const Float_t offh    = noswap ? pedpix.GetPedestalABoffset() : -pedpix.GetPedestalABoffset();
        const Float_t mean[2] = { ped + offh, ped - offh };

        //const Int_t rangeh = fHiGainLast - fHiGainFirst + 1  + fHiLoLast;
        //const Int_t rangel = fLoGainLast - fLoGainFirst + 1;

        const Byte_t *src = numl>0 ? sample.GetArray() : pixel.GetHiGainSamples();
        //const Byte_t *src = sample.GetArray();

        // Copy hi-gains into array and substract pedestal
        // FIXME: Shell we really subtract the pedestal from saturating slices???
        for (Float_t *ptr=beg; ptr<end; ptr++)
            *ptr = (Float_t)*src++ - mean[(ptr-beg)&1];

        // Determin saturation of hi-gains
        Byte_t *p0 = fSignal->GetSamplesRaw(pixidx); //numl>0 ? sample.GetArray() : pixel.GetHiGainSamples();
//        Byte_t *p0 = sample.GetArray();

        Byte_t maxcont  =  0;
        Int_t  maxposhi = -1;

        Byte_t *sathi0 = 0; // first saturating hi-gain slice
        Byte_t *sathi1 = 0; // last saturating hi-gain slice
        for (Byte_t *ptr=p0+fHiGainFirst; ptr<p0+fHiGainLast+fHiLoLast+1; ptr++)
        {
            // Get hi-gain maxbincontent
            // Put this into its own function and loop?
            if (*ptr>maxcont)
            {
                // ORGINAL:
                //Int_t range = (fHiGainLast - fHiGainFirst + 1 + fHiLoLast) - fWindowSizeHiGain + 1;
                // maxpos>1 && maxpos<=range

                // This is a strange workaround to fit the previous
                // Spline setup: TO BE CHECKED!
                const Int_t pos = ptr-p0;
                if (pos<fHiGainLast+1)
                {
                    maxposhi = pos-fHiGainFirst;

                    if (maxposhi>1 && maxposhi<fHiGainLast-fHiGainFirst*//*+1 -fWindowSizeHiGain+1*//*+fHiLoLast*//*)
                        maxcont = *ptr;
                }
            }

            // FIXME: Do we also have to really COUNT the number
            // of saturating slices, eg. for crosschecks?
            if (*ptr>=fSaturationLimit)
            {
                sathi1 = ptr;
                if (!sathi0)
                    sathi0 = ptr;
            }
        }
*/
        Int_t sathi0 = fHiGainFirst;          // First slice to extract and first saturating slice
        Int_t sathi1 = fHiGainLast/*+fHiLoLast*/; // Last  slice to extract and last saturating slice

        Int_t maxcont;
        Int_t maxposhi = fSignal->GetMax(pixidx, sathi0, sathi1, maxcont);
        // Would it be better to take lastsat-firstsat?
        Int_t numsathi = fSignal->GetSaturation(pixidx, fSaturationLimit, sathi0, sathi1);
/*
        // sathi2 is the number (not the index) of first saturating slice
//        const Byte_t sathi2   = sathi0<0 ? 0 : sathi0+1;

        // Number (not index) of first saturating slice
 //       const Byte_t sathi2 = sathi0==0 ? 0 : sathi0-p0+1;
//        const Byte_t numsathi = sathi0==0 ? 0 : sathi1-sathi0+1;

//        Int_t maxb = maxcont;

        // ====================================================================
        // FIXME: Move range out of the extractors...

        // Initialize maxcont for the case, it gets not set by the derived class
        Float_t sumhi2 =0., deltasumhi2 =-1; // Set hi-gain of MExtractedSignalPix valid
        Float_t timehi2=0., deltatimehi2=-1; // Set hi-gain of MArrivalTimePix valid
*/
        Float_t sumhi =0., deltasumhi =-1; // Set hi-gain of MExtractedSignalPix valid
        Float_t timehi=0., deltatimehi=-1; // Set hi-gain of MArrivalTimePix valid
        if (numsathi<1)
        {

        const Int_t rangehi = fHiGainLast - fHiGainFirst + 1/* + fHiLoLast*/;
        FindTimeAndChargeHiGain2(sig/*.GetArray()*/+fHiGainFirst, rangehi,
                                 sumhi, deltasumhi, timehi, deltatimehi,
                                 numsathi, maxposhi);

        timehi += fHiGainFirst;
        }
/*
        Float_t sumhi = sumhi2;
        Float_t deltasumhi = deltasumhi2;
        Float_t timehi = timehi2;
        Float_t deltatimehi=deltatimehi2;
//        Byte_t sathi = sathi2;


        //
        // Find signal in hi- and lo-gain
        //
        Float_t sumhi =0., deltasumhi =-1; // Set hi-gain of MExtractedSignalPix valid
        Float_t timehi=0., deltatimehi=-1; // Set hi-gain of MArrivalTimePix valid
        Byte_t sathi=0;

        // Initialize maxcont for the case, it gets not set by the derived class
        maxcont = fLoGainSwitch + 1;

        //const Int_t pixidx = pixel.GetPixelId();
        //const MPedestalPix  &ped = (*fPedestals)[pixidx];
        const Bool_t higainabflag = pixel.HasABFlag();
        FindTimeAndChargeHiGain(pixel.GetHiGainSamples()+fHiGainFirst, pixel.GetLoGainSamples(),
                                sumhi, deltasumhi, timehi, deltatimehi,
                                sathi, pedpix, higainabflag);
        timehi += fHiGainFirst;

  */
        // Make sure that in cases the time couldn't be correctly determined
        // more meaningfull default values are assigned
        if (timehi<=fHiGainFirst || timehi>=fHiGainLast)
            timehi = gRandom->Uniform(fHiGainLast-fHiGainFirst)+fHiGainFirst;

        Float_t sumlo =0., deltasumlo =-1.; // invalidate logain of MExtractedSignalPix
        Float_t timelo=0., deltatimelo=-1;  // invalidate logain of MArrivalTimePix
        Byte_t satlo=0;
        Int_t numsatlo=0;

        //
        // Adapt the low-gain extraction range from the obtained high-gain time
        //

        // IN THIS CASE THE PIXEL SHOULD BE MARKED BAD!!!!
        // MEANS: Hi gain has saturated, but the signal is to dim
        // to extract the lo-gain properly
        // THIS produces pulse positions ~= -1
        // The signal might be handled in MCalibrateData, but hwat's about
        // the arrival times in MCalibrateRelTime
        if (numsathi>0 && maxcont<=fLoGainSwitch)
            deltasumlo=deltasumhi=deltatimelo=deltatimehi=-1;

        // If more than 8 hi-gain slices have saturated this is
        // no physical event. We just assume that something with
        // the extraction is wrong
        if (numsathi>8) // FIXME: Should be something like 2?
            deltasumhi=deltatimehi=-1;

        // FIXME: What to do with the pixel if it saturates too early???
        if (pixel.HasLoGain() && (maxcont > fLoGainSwitch /*|| sathi>0*/) )
        {
            // To determin the window in which the lo-gain is extracted
            // clearly more information about the relation between the
            // extraction window and the reslting time is necessary.
            //fLoGainStartShift = -6.0;

            const Byte_t fLoGainFirstSave = fLoGainFirst;

            // sathi is the number (not index!) of the first saturating slice
            // 0 indicates that no saturation was found
            // FIXME: Is 0.5 should be expressed by the rise time of
            //        the hi-gain signal!
            const Float_t pos = numsathi==0 ? timehi : sathi0-0.5;

            const Int_t lostart = TMath::FloorNint(pos+6);

            const Int_t newfirst = TMath::FloorNint(pos+fLoGainStartShift);
            fLoGainFirst = newfirst>fLoGainFirstSave ? newfirst : fLoGainFirstSave;

//            if (fLoGainLast-fLoGainFirst >= fWindowSizeLoGain)
            {
/*
                Float_t sumlo2 =0., deltasumlo2 =-1.; // invalidate logain of MExtractedSignalPix
                Float_t timelo2=0., deltatimelo2=-1.;  // invalidate logain of MArrivalTimePix

                // Determin saturation in lo-gains
                Byte_t *satlo0 = 0; // first saturating hi-gain slice
                Byte_t *satlo1 = 0; // last saturating hi-gain slice
                Int_t maxposlo = -1;
                Byte_t maxlo = 0;

                Byte_t *p0 = fSignal->GetSamplesRaw(pixidx);
                for (Byte_t *ptr=p0+numh+fLoGainFirst; ptr<p0+numh+fLoGainLast+1; ptr++)
                {
                    if (*ptr>maxlo)
                    {
                        // This is a strange workaround to fit the previous
                        // Spline setup: TO BE CHECKED!
                        const Int_t ipos = ptr-p0-numh;
                        if (ipos>=fLoGainFirst && ipos<=fLoGainLast)
                        {
                            maxposlo = ipos-fLoGainFirst;
                            maxlo = *ptr;
                        }
                    }
                    if (*ptr>=fSaturationLimit)
                    {
                        satlo1 = ptr;
                        if (!satlo0)
                            satlo0 = ptr;
                    }
                    // FIXME: Do we also have to really COUNT the number
                    // of saturating slices, eg. for crosschecks?
                }

                // Number of saturating lo-gain slices (this is not necessary
                // identical to a counter counting the number of saturating
                // slices!)
                const Byte_t numsatlo = satlo0==0 ? 0 : satlo1-satlo0+1;
*/

                Int_t satlo0 = numh+fLoGainFirst;   // First slice to extract and first saturating slice
                Int_t satlo1 = numh+fLoGainLast;    // Last  slice to extract and last saturating slice

                // Would it be better to take lastsat-firstsat?
                Int_t maxlo;
                Int_t maxposlo = fSignal->GetMax(pixidx, satlo0, satlo1, maxlo);
                numsatlo = fSignal->GetSaturation(pixidx, fSaturationLimit, satlo0, satlo1);

                //                const Byte_t satlo2 = numsatlo;

                const Int_t rangelo = fLoGainLast - fLoGainFirst + 1;
                FindTimeAndChargeLoGain2(sig/*.GetArray()*/+numh+fLoGainFirst, rangelo,
                                         sumlo, deltasumlo, timelo, deltatimelo,
                                         numsatlo, maxposlo);
                timelo += fLoGainFirst;

//                sumlo =sumlo2, deltasumlo=deltasumlo2; // invalidate logain of MExtractedSignalPix
//                timelo=timelo2, deltatimelo=deltatimelo2;  // invalidate logain of MArrivalTimePix
//                satlo = satlo2;

                /*
                // THIS IS NOW THE JOB OF THE EXTRACTOR!
                //deltasumlo  = 0; // make logain of MExtractedSignalPix valid
                //deltatimelo = 0; // make logain of MArrivalTimePix valid

                const Bool_t logainabflag = (higainabflag + pixel.GetNumHiGainSamples()) & 0x1;
                FindTimeAndChargeLoGain(pixel.GetLoGainSamples()+fLoGainFirst,
                                        sumlo, deltasumlo, timelo, deltatimelo,
                                        satlo, pedpix, logainabflag);
                timelo += fLoGainFirst;
                */


                // If more than 6 lo-gain slices have saturated this is
                // no physical event. We just assume that something with
                // the extraction is wrong
                if (numsatlo>6)
                    deltasumlo=deltatimelo=-1;
            }
            fLoGainFirst = fLoGainFirstSave;

            // Make sure that in cases the time couldn't be correctly determined
            // more meaningfull default values are assigned
            if (timelo<=fLoGainFirst || timelo>=fLoGainLast)
                timelo = gRandom->Uniform(fLoGainLast-fLoGainFirst)+fLoGainFirst;
       }

        MExtractedSignalPix &pix = (*fSignals)[pixidx];
        MArrivalTimePix     &tix = (*fArrTime)[pixidx];
        pix.SetExtractedSignal(sumhi, deltasumhi, sumlo, deltasumlo);
        pix.SetGainSaturation(numsathi, numsatlo);
//        pix.SetGainSaturation(sathi, satlo);

        tix.SetArrivalTime(timehi, deltatimehi, timelo-fOffsetLoGain, deltatimelo);
        tix.SetGainSaturation(numsathi, numsatlo);
//        tix.SetGainSaturation(sathi, satlo);
    } /* while (pixel.Next()) */

    fArrTime->SetReadyToSave();
    fSignals->SetReadyToSave();

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// In addition to the resources of the base-class MExtractor:
//   MJPedestal.MExtractor.LoGainStartShift: -2.8
//
Int_t MExtractTimeAndCharge::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{

    Bool_t rc = MExtractTime::ReadEnv(env, prefix, print);

    if (rc)
      SetLoGainStartShift();

    if (IsEnvDefined(env, prefix, "LoGainStartShift", print))
    {
        fLoGainStartShift = GetEnvValue(env, prefix, "LoGainStartShift", fgLoGainStartShift);
        SetLoGainStartShift(fLoGainStartShift);
        rc = kTRUE;
    }

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

    return rc;
}

void MExtractTimeAndCharge::Print(Option_t *o) const
{
  MExtractTime::Print(o);
//  if (IsA()==Class())
//    *fLog << GetDescriptor() << ":" << endl;

  *fLog << dec;
  //*fLog << " Taking " << fWindowSizeHiGain
  //      << " HiGain samples from slice " << (Int_t)fHiGainFirst
  //      << " to " << (Int_t)(fHiGainLast/*+fHiLoLast*/) << " incl" << endl;
  //*fLog << " Taking " << fWindowSizeLoGain
  //      << " LoGain samples from slice " << (Int_t)fLoGainFirst
  //      << " to " << (Int_t)fLoGainLast << " incl" << endl;
  //
  *fLog << " LoGainStartShift:   " << fLoGainStartShift << endl;
  *fLog << " LoGainSwitch:       " << (Int_t)fLoGainSwitch << endl;
}