source: trunk/MagicSoft/Mars/mpedestal/MPedCalcFromLoGain.cc@ 5105

/* ======================================================================== *\
!
! *
! * 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): Josep Flix 04/2001 <mailto:jflix@ifae.es>
!   Author(s): Thomas Bretz 05/2001 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Sebastian Commichau 12/2003 
!   Author(s): Javier Rico 01/2004 <mailto:jrico@ifae.es>
!   Author(s): Markus Gaug 01/2004 <mailto:markus@ifae.es>
!   Author(s): Florian Goebel 06/2004 <mailto:fgoebel@mppmu.mpg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//   MPedCalcLoGain
//
//
//  This task is derived form MPedCalcPedRun, described below. However, It
//  calculates the pedstals using the low gain slices, whenever the difference 
//  between the highest and the lowest slice in the high gain
//  slices is below a given threshold (SetMaxHiGainVar). In this case the receiver
//  boards do not switch to lo gain and the so called lo gain slices are actually
//  high gain slices. 
//
//  MPedCalcLoGain also fills the ABoffset in MPedestalPix which allows to correct 
//  the 150 MHz clock noise.
//
//  This task takes a pedestal run file and fills MPedestalCam during
//  the Process() with the pedestal and rms computed in an event basis.
//  In the PostProcess() MPedestalCam is finally filled with the pedestal
//  mean and rms computed in a run basis.
//  More than one run (file) can be merged
//
//  MPedCalcPedRun applies the following formula (1):
// 
//  Pedestal per slice = sum(x_i) / n / slices
//  PedRMS per slice   = Sqrt(  ( sum(x_i^2) - sum(x_i)^2/n ) / n-1 / slices )
// 
//  where x_i is the sum of "slices" FADC slices and sum means the sum over all
//  events. "n" is the number of events, "slices" is the number of summed FADC samples.
// 
//  Note that the slice-to-slice fluctuations are not Gaussian, but Poissonian, thus 
//  asymmetric and they are correlated.
// 
//  It is important to know that the Pedestal per slice and PedRMS per slice depend 
//  on the number of used FADC slices, as seen in the following plots:
//
//Begin_Html
/*
<img src="images/PedestalStudyInner.gif">
*/
//End_Html
//
//Begin_Html
/*
<img src="images/PedestalStudyOuter.gif">
*/
//
// The plots show the inner and outer pixels, respectivly and have the following meaning:
// 
// 1) The calculated mean pedestal per slice (from MPedCalcFromLoGain)
// 2) The fitted mean pedestal per slice     (from MHPedestalCam)
// 3) The calculated pedestal RMS per slice  (from MPedCalcFromLoGain)
// 4) The fitted sigma of the pedestal distribution per slice
//                                           (from MHPedestalCam)
// 5) The relative difference between calculation and histogram fit
//    for the mean
// 6) The relative difference between calculation and histogram fit
//    for the sigma or RMS, respectively.
// 
// The calculated means do not change significantly except for the case of 2 slices, 
// however the RMS changes from 5.7 per slice in the case of 2 extracted slices 
// to 8.3 per slice in the case of 26 extracted slices. This change is very significant.
// 
// The plots have been produced on run 20123. You can reproduce them using
// the macro pedestalstudies.C
// 
//  Usage of this class: 
//  ====================
// 
//  Call: SetRange(higainfirst, higainlast, logainfirst, logainlast) 
//  to modify the ranges in which the window is allowed to move. 
//  Defaults are: 
// 
//   fHiGainFirst =  fgHiGainFirst =  0 
//   fHiGainLast  =  fgHiGainLast  =  29
//   fLoGainFirst =  fgLoGainFirst =  0 
//   fLoGainLast  =  fgLoGainLast  =  14
//
//  Call: SetWindowSize(windowhigain, windowlogain) 
//  to modify the sliding window widths. Windows have to be an even number. 
//  In case of odd numbers, the window will be modified.
//
//  Defaults are:
//
//   fHiGainWindowSize = fgHiGainWindowSize = 14
//   fLoGainWindowSize = fgLoGainWindowSize = 0
//
//  Variables:
//   fgHiGainFirst;      First FADC slice Hi-Gain (currently set to: 3)
//   fgHiGainLast:       Last FADC slice Hi-Gain (currently set to: 14)
//   fgLoGainFirst:      First FADC slice Lo-Gain (currently set to: 3)
//   fgLoGainLast:       Last FADC slice Lo-Gain (currently set to: 14)
//   fgHiGainWindowSize: The extraction window Hi-Gain
//   fgLoGainWindowSize: The extraction window Lo-Gain
//   fgMaxHiGainVar:     The maximum difference between the highest and lowest slice
//                       in the high gain window allowed in order to use low gain
//
//  Input Containers:
//   MRawEvtData
//   MRawRunHeader
//   MGeomCam
//
//  Output Containers:
//   MPedestalCam
//
//  See also: MPedestalCam, MPedestalPix, MHPedestalCam, MExtractor
//
/////////////////////////////////////////////////////////////////////////////
#include "MPedCalcFromLoGain.h"
#include "MExtractor.h"

#include "MParList.h"

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

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

#include "MPedestalPix.h"
#include "MPedestalCam.h"

#include "MGeomPix.h"
#include "MGeomCam.h"

ClassImp(MPedCalcFromLoGain);

using namespace std;

const Byte_t MPedCalcFromLoGain::fgHiGainFirst      = 0;
const Byte_t MPedCalcFromLoGain::fgHiGainLast       = 11;
const Byte_t MPedCalcFromLoGain::fgLoGainFirst      = 1;
const Byte_t MPedCalcFromLoGain::fgLoGainLast       = 14;
const Byte_t MPedCalcFromLoGain::fgHiGainWindowSize = 12;
const Byte_t MPedCalcFromLoGain::fgLoGainWindowSize = 14;
const Byte_t MPedCalcFromLoGain::fgMaxHiGainVar     = 40;

// --------------------------------------------------------------------------
//
// Default constructor: 
//
// Sets:
// - all pointers to NULL
// - fWindowSizeHiGain to fgHiGainWindowSize
// - fWindowSizeLoGain to fgLoGainWindowSize
//
// Calls: 
// - AddToBranchList("fHiGainPixId");
// - AddToBranchList("fHiGainFadcSamples");
// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
// - Clear()
//
MPedCalcFromLoGain::MPedCalcFromLoGain(const char *name, const char *title)
    : fWindowSizeHiGain(fgHiGainWindowSize), 
      fWindowSizeLoGain(fgLoGainWindowSize), 
      fGeom(NULL), fPedContainerName("MPedestalCam")
{
    fName  = name  ? name  : "MPedCalcFromLoGain";
    fTitle = title ? title : "Task to calculate pedestals from pedestal runs raw data";

    AddToBranchList("fHiGainPixId");
    AddToBranchList("fLoGainPixId");
    AddToBranchList("fHiGainFadcSamples");
    AddToBranchList("fLoGainFadcSamples");

    SetRange();
    SetMaxHiGainVar();
    SetPedestalUpdate(kTRUE);
    SetNumEventsDump(500);

    Clear();
}

// --------------------------------------------------------------------------
//
// Sets:
// - fNumSamplesTot to 0
// - fRawEvt to NULL
// - fRunHeader to NULL
// - fPedestals to NULL
//
void MPedCalcFromLoGain::Clear(const Option_t *o)
{
    fRawEvt    = NULL;
    fRunHeader = NULL;
    fPedestals = NULL;

    // If the size is yet set, set the size
    if (fSumx.GetSize()>0)
    {
        // Reset contents of arrays.
        fSumx.Reset();
        fSumx2.Reset();
        fSumAB0.Reset();
        fSumAB1.Reset();
        fNumEventsUsed.Reset();
        fTotalCounter.Reset();
    }
}

// --------------------------------------------------------------------------
//
// SetRange: 
//
// Calls:
// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
// - SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
//
void MPedCalcFromLoGain::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
{
    MExtractor::SetRange(hifirst, hilast, lofirst, lolast);

    //
    // Redo the checks if the window is still inside the ranges
    //
    SetWindowSize(fWindowSizeHiGain,fWindowSizeLoGain);
}

// --------------------------------------------------------------------------
//
void MPedCalcFromLoGain::SetMaxHiGainVar(Byte_t maxvar)
{
    fMaxHiGainVar = maxvar;
}

// --------------------------------------------------------------------------
//
// Checks:
// - if a window is odd, subtract one
// - if a window is bigger than the one defined by the ranges, set it to the available range
// - if a window is smaller than 2, set it to 2
// 
// Sets:
// - fNumHiGainSamples to: (Float_t)fWindowSizeHiGain
// - fNumLoGainSamples to: (Float_t)fWindowSizeLoGain
// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)  
//  
void MPedCalcFromLoGain::SetWindowSize(Byte_t windowh, Byte_t windowl)
{
    fWindowSizeHiGain = windowh & ~1;
    fWindowSizeLoGain = windowl & ~1;
  
    if (fWindowSizeHiGain != windowh)
    {
        *fLog << warn;
        *fLog << GetDescriptor() << ": HiGain window has to be even, set to: ";
        *fLog << int(fWindowSizeHiGain) << " samples " << endl;
    }
    
    if (fWindowSizeLoGain != windowl)
    {
        *fLog << warn;
        *fLog << GetDescriptor() << ": Lo Gain window has to be even, set to: ";
        *fLog << int(fWindowSizeLoGain) << " samples " << endl;
    }
      
    const Byte_t availhirange = (fHiGainLast-fHiGainFirst+1) & ~1;
    const Byte_t availlorange = (fLoGainLast-fLoGainFirst+1) & ~1;
  
    if (fWindowSizeHiGain > availhirange)
    {
        *fLog << warn;
        *fLog << GetDescriptor() << ": HiGain window " << (int)fWindowSizeHiGain;
        *fLog << " out of range [" << (int)fHiGainFirst;
        *fLog << "," << (int)fHiGainLast << "]" << endl;
        *fLog << "Will set window size to " << (int)availhirange << endl;
        fWindowSizeHiGain = availhirange;
    }
    
    if (fWindowSizeLoGain > availlorange)
    {
        *fLog << warn;
        *fLog << GetDescriptor() << ": LoGain window " << (int)fWindowSizeLoGain;
        *fLog << " out of range [" << (int)fLoGainFirst;
        *fLog << "," << (int)fLoGainLast << "]" << endl;
        *fLog << "Will set window size to " << (int)availlorange << endl;
        fWindowSizeLoGain = availlorange;
    }
    /*
     fNumHiGainSamples = (Float_t)fWindowSizeHiGain;
     fNumLoGainSamples = (Float_t)fWindowSizeLoGain;

     fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
     fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);
    */
}

// --------------------------------------------------------------------------
//
// Look for the following input containers:
//
//  - MRawEvtData
//  - MRawRunHeader
//  - MGeomCam
// 
// The following output containers are also searched and created if
// they were not found:
//
//  - MPedestalCam
//
Int_t MPedCalcFromLoGain::PreProcess(MParList *pList)
{
    Clear();

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

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

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

    fPedestals = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", AddSerialNumber(fPedContainerName));
    if (!fPedestals)
        return kFALSE;

    if (fNumEventsDump<=0 && fPedestalUpdate)
    {
        *fLog << warn << "Pedestal Update switched on and Number of Events to dump <= 0... fNumEventsDump=1000" << endl;
        fNumEventsDump=1000;
    }

    *fLog << inf << "Continous update switched " << (fPedestalUpdate?"on":"off");
    if (fPedestalUpdate)
        *fLog << " (dump each " << fNumEventsDump << " events)" << endl;
    *fLog << endl;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// The ReInit searches for:
// -  MRawRunHeader::GetNumSamplesHiGain()
// -  MRawRunHeader::GetNumSamplesLoGain()
//
// In case that the variables fHiGainLast and fLoGainLast are smaller than 
// the even part of the number of samples obtained from the run header, a
// warning is given an the range is set back accordingly. A call to:  
// - SetRange(fHiGainFirst, fHiGainLast-diff, fLoGainFirst, fLoGainLast) or 
// - SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast-diff) 
// is performed in that case. The variable diff means here the difference 
// between the requested range (fHiGainLast) and the available one. Note that 
// the functions SetRange() are mostly overloaded and perform more checks, 
// modifying the ranges again, if necessary.
//
Bool_t MPedCalcFromLoGain::ReInit(MParList *pList)
{
    if (fRunHeader->GetNumSamplesHiGain()<1)
    {
        *fLog << err << "ERROR - Number of available HiGainSamples<1... abort." << endl;
        return kFALSE;
    }
    if (fRunHeader->GetNumSamplesLoGain()<1)
    {
        *fLog << err << "ERROR - Number of available LoGainSamples<1... abort." << endl;
        return kFALSE;
    }

    fHiGainFirst = TMath::Max(fHiGainFirst, 0);
    fHiGainLast  = TMath::Min(fHiGainLast,  fRunHeader->GetNumSamplesHiGain()-1);

    fLoGainFirst = TMath::Max(fLoGainFirst, 0);
    fLoGainLast  = TMath::Min(fLoGainLast,  fRunHeader->GetNumSamplesLoGain()-1);

    const Double_t wh = fWindowSizeHiGain;
    const Double_t wl = fWindowSizeLoGain;
    const Double_t fh = fHiGainFirst;
    const Double_t fl = fLoGainFirst;;

    fWindowSizeHiGain = TMath::Min(fWindowSizeHiGain, fHiGainLast-fHiGainFirst+1);
    fWindowSizeLoGain = TMath::Min(fWindowSizeLoGain, fLoGainLast-fLoGainFirst+1);

    SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast);

    if (wh!=fWindowSizeHiGain || fh!=fHiGainFirst || wl!=fWindowSizeLoGain || fl!=fLoGainFirst)
    {
        *fLog << inf << endl;
        *fLog << "Taking " << Form("%2d", (int)fWindowSizeHiGain) << " slices of Hi-Gain starting at slice " << (int)fHiGainFirst << endl;
        *fLog << "Taking " << Form("%2d", (int)fWindowSizeLoGain) << " slices of Lo-Gain starting at slice " << (int)fLoGainFirst << endl;
    }

    if (fWindowSizeHiGain==0)
    {
        *fLog << err << "ERROR - HiGain windows size == 0... abort." << endl;
        return kFALSE;
    }
    if (fWindowSizeLoGain==0)
    {
        *fLog << err << "ERROR - HiGain windows size == 0... abort." << endl;
        return kFALSE;
    }

    // If the size is not yet set, set the size
    if (fSumx.GetSize()==0)
    {
        const Int_t npixels = fPedestals->GetSize();

        fSumx. Set(npixels);
        fSumx2.Set(npixels);
        fSumAB0.Set(npixels);
        fSumAB1.Set(npixels);
        fNumEventsUsed.Set(npixels);
        fTotalCounter.Set(npixels);

        // Reset contents of arrays.
        fSumx.Reset();
        fSumx2.Reset();
        fSumAB0.Reset();
        fSumAB1.Reset();
        fNumEventsUsed.Reset();
        fTotalCounter.Reset();
    }

    return kTRUE;
}

void MPedCalcFromLoGain::Calc(ULong_t n, UInt_t idx)
{
    const ULong_t nsamplestot = n*fWindowSizeLoGain;

    const Float_t sum  = fSumx.At(idx);
    const Float_t sum2 = fSumx2.At(idx);
    const Float_t ped  = sum/nsamplestot;

    // 1. Calculate the Variance of the sums:
    Float_t var = (sum2-sum*sum/n)/(n-1.);

    // 2. Scale the variance to the number of slices:
    var /= (Float_t)(fWindowSizeLoGain);

    // 3. Calculate the RMS from the Variance:
    const Float_t rms = var<0 ? 0 : TMath::Sqrt(var);

    // 4. Calculate the amplitude of the 150MHz "AB" noise
    const Float_t abOffs = (fSumAB0[idx] - fSumAB1[idx]) / nsamplestot;

    (*fPedestals)[idx].Set(ped, rms, abOffs, n);

    fTotalCounter[idx]++;
}

// --------------------------------------------------------------------------
//
// Fill the MPedestalCam container with the signal mean and rms for the event.
// Store the measured signal in arrays fSumx and fSumx2 so that we can 
// calculate the overall mean and rms in the PostProcess()
//
Int_t MPedCalcFromLoGain::Process()
{
    MRawEvtPixelIter pixel(fRawEvt);

    while (pixel.Next())
    {
        const UInt_t idx = pixel.GetPixelId();

        Byte_t *ptr = pixel.GetHiGainSamples() + fHiGainFirst;
        Byte_t *end = ptr + fWindowSizeHiGain;

        UInt_t sum = 0;
        UInt_t sqr = 0;

        UInt_t max = 0;
        UInt_t min = 255;
    
        // Find the maximum and minimum signal per slice in the high gain window
        do {
            if (*ptr > max)
                max = *ptr;
            if (*ptr < min)
                min = *ptr;
        } while (++ptr != end);

        // If the maximum in the high gain window is smaller than
        if (max-min>=fMaxHiGainVar || max>=255)
            continue;

        ptr = pixel.GetLoGainSamples() + fLoGainFirst;
        end = ptr + fWindowSizeLoGain;

        Byte_t *firstSlice = ptr;

        do {
            sum += *ptr;
            sqr += *ptr * *ptr;
        } while (++ptr != end);

        const Float_t msum   = (Float_t)sum;
        const Float_t sqrsum = msum*msum;

        fSumx[idx]  += msum;
        fSumx2[idx] += sqrsum;
        fNumEventsUsed[idx]++;

        // Calculate the amplitude of the 150MHz "AB" noise

        if (pixel.IsABFlagValid())
        {
            const Int_t abFlag = (fRunHeader->GetNumSamplesHiGain()
                                  + fLoGainFirst + pixel.HasABFlag()) & 0x1;
            for (Int_t islice=0; islice<fWindowSizeLoGain; islice+=2)
            {
                const Int_t sliceAB0 = islice + abFlag;
                const Int_t sliceAB1 = islice - abFlag + 1;
                fSumAB0[idx] += firstSlice[sliceAB0];
                fSumAB1[idx] += firstSlice[sliceAB1];
            }
        }

        if (!fPedestalUpdate || fNumEventsUsed[idx]<fNumEventsDump)
            continue;

        Calc(fNumEventsDump, idx);

        fNumEventsUsed[idx]=0;
        fSumx[idx]=0;
        fSumx2[idx]=0;
        fSumAB0[idx]=0;
        fSumAB1[idx]=0;
    }

    if (fPedestalUpdate)
    {
        //fPedestals->ReCalc(*fGeom);
        fPedestals->SetReadyToSave();
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Compute signal mean and rms in the whole run and store it in MPedestalCam
//
Int_t MPedCalcFromLoGain::PostProcess()
{
    // Compute pedestals and rms from the whole run
    if (fPedestalUpdate || GetNumExecutions()<1)
        return kTRUE;

    *fLog << flush << inf << "Calculating pedestals..." << flush;

    Double_t sum = 0;

    const Int_t npix = fGeom->GetNumPixels();
    for (Int_t idx=0; idx<npix; idx++)
    {
        const ULong_t n = fNumEventsUsed[idx];
        if (n>1)
            Calc(n, idx);
        sum += n;
    }

    *fLog << flush << inf << "Calculating means..." << flush;

    fPedestals->SetTotalEntries((UInt_t)(sum/npix*(fWindowSizeLoGain+fWindowSizeHiGain)));
    fPedestals->ReCalc(*fGeom);
    fPedestals->SetReadyToSave();
    return kTRUE;
}

Int_t MPedCalcFromLoGain::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    if (MExtractor::ReadEnv(env, prefix, print)==kERROR)
        return kERROR;

    Bool_t rc=kFALSE;

    Byte_t hw = fWindowSizeHiGain;
    Byte_t lw = fWindowSizeLoGain;
    if (IsEnvDefined(env, prefix, "WindowSizeHiGain", print))
    {
        hw = GetEnvValue(env, prefix, "WindowSizeHiGain", hw);
        rc = kTRUE;
    }
    if (IsEnvDefined(env, prefix, "WindowSizeLoGain", print))
    {
        lw = GetEnvValue(env, prefix, "WindowSizeLoGain", lw);
        rc = kTRUE;
    }

    if (rc)
        SetWindowSize(hw, lw);

    if (IsEnvDefined(env, prefix, "NumEventsDump", print))
    {
        SetNumEventsDump(GetEnvValue(env, prefix, "NumEventsDump", fNumEventsDump));
        rc = kTRUE;
    }

    if (IsEnvDefined(env, prefix, "MaxHiGainVar", print))
    {
        SetMaxHiGainVar(GetEnvValue(env, prefix, "MaxHiGainVar", fMaxHiGainVar));
        rc = kTRUE;
    }

    if (IsEnvDefined(env, prefix, "PedestalUpdate", print))
    {
        SetPedestalUpdate(GetEnvValue(env, prefix, "PedestalUpdate", fPedestalUpdate));
        rc = kTRUE;
    }

    return rc;
}