source: trunk/MagicSoft/Mars/manalysis/MArrivalTimeCalc.cc@ 3188

/* ======================================================================== *\
!
! *
! * 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): Sebastian Raducci 12/2003 <mailto:raducci@fisica.uniud.it>
!
!   Copyright: MAGIC Software Development, 2002-2004
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//   MArrivalTimeCalc
//
//   This is a task that calculates the arrival times of photons. 
//   It returns the absolute maximum of the spline that interpolates
//   the FADC slices 
//
// Input Containers:
//   MRawEvtData
//
// Output Containers:
//   MArrivalTime
//   MRawEvtData
//
//////////////////////////////////////////////////////////////////////////////
#include "MArrivalTimeCalc.h"

#include <TSpline.h>

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

#include "MParList.h"

#include "MGeomCam.h"
#include "MArrivalTime.h"
#include "MRawEvtData.h"      
#include "MRawEvtPixelIter.h"

ClassImp(MArrivalTimeCalc);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor.
//
// Initialize step size by default to 0.03 time slices == 100 ps.
// 
MArrivalTimeCalc::MArrivalTimeCalc(const char *name, const char *title)
    : fStepSize(0.03) 
{

    fName  = name  ? name  : "MArrivalTimeCalc";
    fTitle = title ? title : "Calculate photons arrival time";

}

// --------------------------------------------------------------------------
//
// The PreProcess searches for the following input containers:
//  - MRawEvtData
//  - MArrivalTime
//
// The following output containers are also searched and created if
// they were not found:
//  - MArrivalTime
//

Int_t MArrivalTimeCalc::PreProcess(MParList *pList)
{


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

    fArrTime = (MArrivalTime*)pList->FindCreateObj(AddSerialNumber("MArrivalTime"));
    if (!fArrTime)
        return kFALSE;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// The ReInit searches for the following input containers:
//  - MGeomCam
//
Bool_t MArrivalTimeCalc::ReInit(MParList *pList)
{
    MGeomCam *cam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));
    if (!cam)
    {
        *fLog << err << GetDescriptor() << ": No MGeomCam found... aborting." << endl;
        return kFALSE;
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Evaluation of the mean arrival times (spline interpolation)
// per pixel and store them in the MArrivalTime container.
//
Int_t MArrivalTimeCalc::Process()
{

  MRawEvtPixelIter pixel(fRawEvt);
  
  while (pixel.Next())
    {
      
      const UInt_t idx = pixel.GetPixelId();
      Float_t max = 0.;


      //
      // If pixel is saturated we use LoGains
      //
      if (pixel.GetMaxHiGainSample() == 0xff && pixel.HasLoGain())
        {
          
          const Short_t nslices = fRawEvt->GetNumLoGainSamples();
          max = Calc(pixel.GetLoGainSamples(),nslices);
        }
      

      //
      // Use HiGains
      //
      else if (pixel.HasLoGain())
        {
          
          const Short_t nslices = fRawEvt->GetNumHiGainSamples();
          max = Calc(pixel.GetHiGainSamples(),nslices);
        }
      
      //
      // If pixel is saturated and hasn't lo gains we do nothing, it's value remains -1
      //
      fArrTime->SetTime(idx,max);
      
    }
    
  fArrTime->SetReadyToSave();
  
  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Calculates the arrival time for each pixel 
// Possible Methods 
// Case 1: Spline5 (From TSpline5 Root Class)
//
Float_t MArrivalTimeCalc::Calc(const Byte_t *fadcSamples, const Short_t nslices)
{

  //
  // Initialize a double pointer with filled FADC slices
  //
  Double_t ptr[nslices];

  //
  // Initialize the spline
  //
  for (Int_t i=0; i<nslices; i++)
    ptr[i]=(Double_t)fadcSamples[i];
  
  TSpline5 spline("spline",0.,(Double_t)(nslices - 1),ptr,nslices);

  //
  // Now find the half maximum (!) 
  // evaluating the spline function at every fStepSize time slice
  //
  Double_t abscissa=0;
  Double_t maxAb   =0;
  Double_t maxOrd  =0;

  while (abscissa <= nslices - 1) 
    {
      const Double_t swap = spline.Eval(abscissa);

      if (swap > maxOrd)
        {
          maxOrd = swap;
          maxAb  = abscissa;
        }
      // make step size a bit bigger first
      abscissa += fStepSize;
      //      abscissa += fStepSize;
    }

  //
  // another (much smaller) loop to move back from the maximum 
  //
  Double_t halfMaxAb = 0;

  abscissa = maxAb;
  
  while (abscissa > 0)
    {
      
      const Double_t swap = spline.Eval(abscissa);
      
      if (swap < maxOrd/2.)
        {
          halfMaxAb = abscissa;
          break;
        }
      
      abscissa -= fStepSize;
    }
  

  //  return (Float_t)maxAb;
  return (Float_t)halfMaxAb;
}