source: trunk/MagicSoft/Mars/mcalib/MExtractPINDiode.cc@ 3193

/* ======================================================================== *\
!
! *
! * 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, 02/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MExtractPINDiode
//
//////////////////////////////////////////////////////////////////////////////
#include "MExtractPINDiode.h"

#include <fstream>

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

#include "MParList.h"

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

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

#include "MExtractedSignalPINDiode.h"

ClassImp(MExtractPINDiode);

using namespace std;

const UInt_t MExtractPINDiode::fgPINDiodeIdx     = 999;
const Byte_t MExtractPINDiode::fgSaturationLimit = 254;
const Byte_t MExtractPINDiode::fgFirst =  1;
const Byte_t MExtractPINDiode::fgLast  = 30;

// --------------------------------------------------------------------------
//
// Default constructor. 
//
MExtractPINDiode::MExtractPINDiode(const char *name, const char *title)
    : fSaturationLimit(fgSaturationLimit)
{
  
  fName  = name  ? name  : "MExtractPINDiode";
  fTitle = title ? title : "Task to extract the signal from the FADC slices";
  
  AddToBranchList("MRawEvtData.*");
  
  SetPINDiodeIdx();
  SetSaturationLimit();
  SetRange();
}

void MExtractPINDiode::SetRange(Byte_t first, Byte_t last)
{

    fNumSamples = last-first+1;
    fFirst      = first;
    fLast       = last;

    fSqrtSamples = TMath::Sqrt((Float_t)fNumSamples);
}

// --------------------------------------------------------------------------
//
// The PreProcess searches for the following input containers:
//  - MRawEvtData
//  - MPedestalCam
//
// The following output containers are also searched and created if
// they were not found:
//
//  - MExtractedPINDiode
//
Int_t MExtractPINDiode::PreProcess(MParList *pList)
{
    fRawEvt = (MRawEvtData*)pList->FindObject(AddSerialNumber("MRawEvtData"));
    if (!fRawEvt)
    {
        *fLog << err << AddSerialNumber("MRawEvtData") << " not found... aborting." << endl;
        return kFALSE;
    }


    fPINDiode = (MExtractedSignalPINDiode*)pList->FindCreateObj(AddSerialNumber("MExtractedSignalPINDiode"));
    if (!fPINDiode)
        return kFALSE;

    fPINDiode->SetUsedFADCSlices(fFirst, fFirst+fNumSamples-1);

    fPedestals = (MPedestalCam*)pList->FindObject(AddSerialNumber("MPedestalCam"));

    if (!fPedestals)
    {
        *fLog << err << AddSerialNumber("MPedestalCam") << " not found... aborting" << endl;
        return kFALSE;
    }

    return kTRUE;
}

void MExtractPINDiode::FindSignal(Byte_t *ptr, Int_t size, UInt_t &sum, UInt_t &sum2, UInt_t &sat, UInt_t &max) const
{

  Byte_t *end = ptr + size;
  
  while (ptr<end)
    {
      sum  += *ptr;
      sum2 += *ptr * *ptr;
      if (*ptr > max)
        max = *ptr;
      
      if (*ptr++ >= fSaturationLimit)
        sat++;
    }
}

// --------------------------------------------------------------------------
//
// Calculate the integral of the FADC time slices and store them as a new
// pixel in the MExtractedPINDiode container.
//
Int_t MExtractPINDiode::Process()
{

    MRawEvtPixelIter pixel(fRawEvt);

    fPINDiode->Clear();

    pixel.Jump(fPINDiodeIdx);
 
    UInt_t sum  = 0;
    UInt_t sum2 = 0;
    UInt_t sat  = 0;
    UInt_t max  = 0;

    FindSignal(pixel.GetHiGainSamples()+fFirst-1, 
               fRawEvt->GetNumHiGainSamples(), 
               sum, sum2, sat, max);
    FindSignal(pixel.GetLoGainSamples(), 
               fLast-fRawEvt->GetNumLoGainSamples(), 
               sum, sum2, sat, max);
    
    const MPedestalPix &ped   = (*fPedestals)[fPINDiodeIdx]; 

    Float_t pedes  = -999.;
    Float_t pedrms = -999.;

    if (&ped)
      {
        pedes  = ped.GetPedestal();
        pedrms = ped.GetPedestalRms();
      }
    else
      {
        *fLog << warn << " Cannot find MPedestalPix of the PIN Diode (idx=" 
              << fPINDiodeIdx << ")" << endl;

      }

    const Float_t var = ((Float_t)sum2 - (Float_t)sum*sum/fNumSamples)/(fNumSamples-1);
    const Float_t rms = TMath::Sqrt(var);
    
    // 
    // FIXME: The following formulae have to be revised!!
    //
    fPINDiode->SetExtractedSignal(sum - pedes*fNumSamples, pedrms*fSqrtSamples);
    fPINDiode->SetExtractedRms   (rms, rms/2./fSqrtSamples);
    fPINDiode->SetExtractedTime  (max, rms/fSqrtSamples);
    fPINDiode->SetSaturation(sat);

    if (sat)
      *fLog << warn << "WARNING - saturation occurred in the PIN Diode " << endl;

    fPINDiode->SetReadyToSave();

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Implementation of SavePrimitive. Used to write the call to a constructor
// to a macro. In the original root implementation it is used to write
// gui elements to a macro-file.
//
void MExtractPINDiode::StreamPrimitive(ofstream &out) const
{
    out << "   " << ClassName() << " " << GetUniqueName() << "(\"";
    out << "\"" << fName << "\", \"" << fTitle << "\");" << endl;

    if (fSaturationLimit!=fgSaturationLimit)
    {
        out << "   " << GetUniqueName() << ".SetSaturationLimit(";
        out << (int)fSaturationLimit << ");" << endl;
    }

    const Bool_t arg2 = fNumSamples+fFirst-1 != fgLast;
    const Bool_t arg1 = arg2 || fFirst != fgFirst;

    if (!arg1)
        return;

    out << "   " << GetUniqueName() << ".SetRange(";
    out << (int)fFirst;
    if (arg2)
      out << ", " << (int)(fNumSamples+fFirst-1);
    out << ");" << endl;
}