source: trunk/MagicSoft/Mars/msignal/MExtractPINDiode.cc@ 3947

/* ======================================================================== *\
!
! *
! * 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
//
//  Extracts the signal from a fixed window in a given range.
//
//  Call: SetRange(fHiGainFirst, fHiGainLast, fLoGainFirst, fLoGainLast) 
//  to modify the ranges. Ranges have to be an even number. In case of odd 
//  ranges, the last slice will be reduced by one.
//  Defaults are: 
// 
//   fHiGainFirst =  fgHiGainFirst =  3 
//   fHiGainLast  =  fgHiGainLast  =  14
//   fLoGainFirst =  fgLoGainFirst =  3 
//   fLoGainLast  =  fgLoGainLast  =  14
//
//////////////////////////////////////////////////////////////////////////////
#include "MExtractPINDiode.h"
#include "MExtractor.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     = 100;
const Byte_t MExtractPINDiode::fgHiGainFirst =  0;
const Byte_t MExtractPINDiode::fgHiGainLast  =  14;
const Byte_t MExtractPINDiode::fgLoGainFirst =  0;
const Byte_t MExtractPINDiode::fgLoGainLast  =  14;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Calls: 
// - SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast)
// - SetPINDiodeIdx()
//
MExtractPINDiode::MExtractPINDiode(const char *name, const char *title)
{
  
  fName  = name  ? name  : "MExtractPINDiode";
  fTitle = title ? title : "Task to extract the signal from the FADC slices";
  
  SetRange(fgHiGainFirst, fgHiGainLast, fgLoGainFirst, fgLoGainLast);
  SetPINDiodeIdx();
}

// --------------------------------------------------------------------------
//
// SetRange: 
//
// Checks: 
// - if the window defined by (fHiGainLast-fHiGainFirst-1) are odd, subtract one
// - if the window defined by (fLoGainLast-fLoGainFirst-1) are odd, subtract one
// - if the Hi Gain window is smaller than 2, set fHiGainLast to fHiGainFirst+1
// - if the Lo Gain window is smaller than 2, set fLoGainLast to fLoGainFirst+1
// 
// Calls:
// - MExtractor::SetRange(hifirst,hilast,lofirst,lolast);
// 
// Sets:
// - fNumHiGainSamples to: (Float_t)(fHiGainLast-fHiGainFirst+1)
// - fNumLoGainSamples to: (Float_t)(fLoGainLast-fLoGainFirst+1)
// - fSqrtHiGainSamples to: TMath::Sqrt(fNumHiGainSamples)
// - fSqrtLoGainSamples to: TMath::Sqrt(fNumLoGainSamples)  
//  
void MExtractPINDiode::SetRange(Byte_t hifirst, Byte_t hilast, Byte_t lofirst, Byte_t lolast)
{

  const Byte_t window = hilast-hifirst+1+lolast-lofirst+1;
  const Byte_t weven  = window & ~1;

  if (weven != window)
    {
      *fLog << warn << GetDescriptor() 
            << Form("%s%2i%s%2i",": Total window size has to be even, set last slice from "
                    ,(int)lolast," to ",(int)(lolast-1)) << endl;
      lolast -= 1;
    }
  
  if (window<2) 
    {
      *fLog << warn << GetDescriptor() 
            << Form("%s%2i%s%2i",": Total window is smaller than 2 FADC sampes, set last slice from" 
                    ,(int)lolast," to ",(int)(lofirst+1)) << endl;
      hilast = hifirst+1;
    }
  

  MExtractor::SetRange(hifirst,hilast,lofirst,lolast);


  fNumHiGainSamples = (Float_t)(fHiGainLast-fHiGainFirst+1);
  fNumLoGainSamples = (Float_t)(fLoGainLast-fLoGainFirst+1);  

  fSqrtHiGainSamples = TMath::Sqrt(fNumHiGainSamples);
  fSqrtLoGainSamples = TMath::Sqrt(fNumLoGainSamples);  
  
  fNumSamples = fHiGainLast-fHiGainFirst+1+fLoGainLast-fLoGainFirst+1;
  fSqrtSamples = TMath::Sqrt((Float_t)fNumSamples);
  
}

// --------------------------------------------------------------------------
//
// Calls: 
// - MExtractor::PreProcess
//
// The following output containers are also searched and created if
// they were not found:
//
//  - MExtractedPINDiode
//
Int_t MExtractPINDiode::PreProcess(MParList *pList)
{

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

  fPINDiode->SetUsedFADCSlices(fHiGainFirst, fLoGainLast);

  const MPedestalPix &ped   = (*fPedestals)[fPINDiodeIdx]; 

    if (&ped)
      {
        fPedestal = ped.GetPedestal();
        fPedRms   = ped.GetPedestalRms();
      }
    else
      {
        *fLog << err << " Cannot find MPedestalPix of the PIN Diode (idx=" 
              << fPINDiodeIdx << ")" << endl;
        return kFALSE;
      }

    return kTRUE;
}

void MExtractPINDiode::FindSignalandVarianceHiGain(Byte_t *ptr, Int_t &sum, Int_t &sum2, Byte_t &sat) const
{

  Byte_t *end = ptr + fHiGainLast - fHiGainFirst + 1;

  while (ptr<end)
    {
      sum  += *ptr;
      sum2 += *ptr * *ptr;

      if (*ptr++ >= fSaturationLimit)
        sat++;
    }
}

void MExtractPINDiode::FindSignalandVarianceLoGain(Byte_t *ptr, Int_t &sum, Int_t &sum2, Byte_t &sat) const
{

  Byte_t *end = ptr +  fLoGainLast - fLoGainFirst  + 1;

  while (ptr<end)
    {
      sum  += *ptr;
      sum2 += *ptr * *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);
  
  Int_t sum   = 0;
  Int_t sum2  = 0;
  Byte_t sat  = 0;
  Int_t max   = 0;
  
  // 
  // Calculate first the time:
  //
  const Int_t maxhi = pixel.GetIdxMaxHiGainSample();
  const Int_t maxlo = pixel.GetIdxMaxLoGainSample();
  
  if (maxhi > maxlo)
    max = maxhi;
  else
    max = maxlo + pixel.GetNumHiGainSamples();
  
  FindSignalandVarianceHiGain(pixel.GetHiGainSamples()+fHiGainFirst,sum,sum2,sat);
  FindSignalandVarianceLoGain(pixel.GetLoGainSamples()+fLoGainFirst,sum,sum2,sat);

  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 - fPedestal*fNumSamples, fPedRms*fSqrtSamples);
  fPINDiode->SetExtractedRms   (rms, rms/2./fSqrtSamples);
  fPINDiode->SetExtractedTime  (max, rms/fSqrtSamples);
  fPINDiode->SetSaturation(sat);
  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+fHiGainFirst-1 != fgLoGainLast;
    const Bool_t arg1 = arg2 || fHiGainFirst != fgHiGainFirst;

    if (!arg1)
        return;

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