source: trunk/MagicSoft/Mars/manalysis/MArrivalTime.cc@ 2903

/* ======================================================================== *\
!
! *
! * 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, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MArrivalTime
//
// Times are calculated using the TSpline5 Root Class
// 
/////////////////////////////////////////////////////////////////////////////
#include "MArrivalTime.h"

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

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

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

ClassImp(MArrivalTime);

using namespace std;

// --------------------------------------------------------------------------
//
// Creates an object containing the arrival time for each pixel in the event
//
MArrivalTime::MArrivalTime(const char *name, const char *title)
{
    fName  = name  ? name  : "MArrivalTime";
    fTitle = title ? title : "Photons arrival times Information";
}

// Sets every pixel arrival time to -1
void MArrivalTime::CleanArray(const MGeomCam &geom)
{
    fData.Set(geom.GetNumPixels());
    fData.Reset(-1.0);
}

void MArrivalTime::InitSize(Int_t i)
{
  
  if (fData.GetSize() >= i)
    return;
  else
    fData.Set(i);
}

// 
// Set the arrival time in one pixel
//
void MArrivalTime::SetTime(Int_t i, Float_t t)
{

  if (i >= fData.GetSize())
      fData.Set(i+1);

  fData.AddAt(t,i);
}

//
// Calculates the arrival time for each pixel 
// Possible Methods 
// Case 1: Spline5 (From TSpline5 Root Class)
// 
//

void MArrivalTime::Calc(const Byte_t *fadcSamples, const Short_t nSlice,
              const Short_t idx, const MGeomCam &geom)
{

//Some casts are needed because TSpline5 constructor accepts only Double_t values
    Double_t ptr[nSlice];

    for (Int_t i = 0; i < nSlice; i++)
        ptr[i]=(Double_t)fadcSamples[i];

    TSpline5 *spline = new TSpline5("spline",(Double_t) 0,
                                    (Double_t)(nSlice - 1),ptr,nSlice);
//Now find the maximum evaluating the spline function at every 1/10 time slice
    Double_t abscissa=0.0;
    Double_t maxAb=0.0;
    Double_t maxOrd=0.0;
    Double_t swap;
    while (abscissa <= nSlice - 1) 
    {
        swap=spline->Eval(abscissa);
        if (swap > maxOrd)
        {
            maxOrd = swap;
            maxAb  = abscissa;
        }
        abscissa += 0.1;
    }
    fData[idx]=maxAb;
}

// -------------------------------------------------------------------------
//
// Finds the clusters with similar Arrival Time
//
// PRELIMINARY!! For now the Arrival Time is not the one calculated with
//               the splines, but it's the Max Slice Idx
//
// TEST!!        This method is used only for test. Do not use it until
//               it becomes stable
//

void MArrivalTime::EvalClusters(const MRawEvtData &evt, const MGeomCam &geom)
{
    const Int_t n = geom.GetNumPixels();

    fData2.Set(n);
    fData2.Reset(-1.0);
    fData3.Set(n);
    fData3.Reset(-1.0);
    fData4.Set(n);
    fData4.Reset(-1.0);
    fData5.Set(n);
    fData5.Reset(-1.0);

    MRawEvtPixelIter pixel((MRawEvtData*)&evt);

// Array that says if a pixel has been already checked or not

    for (Int_t i = 0; i < n; i++)
        fPixelChecked[i] = kFALSE;

// This fakeData array will be subsituted with the fData Array.
    fakeData.Set(n);
    fakeData.Reset();
   
    while ( pixel.Next() )
    {
        fakeData[pixel.GetPixelId()]=pixel.GetIdxMaxHiLoGainSample();
    }
// End of fakeData preparation
    *fLog << warn << "fin qui bene" << endl;
// Max dimension of cluster
    Short_t dimCluster;

    while ( pixel.Next() )
    {
        if (!fPixelChecked[pixel.GetPixelId()])
        {
            dimCluster = 0;
            fCluster.Set(n);
            fCluster.Reset(-1);
            fCluster[dimCluster]=pixel.GetPixelId();
            dimCluster++;

            CheckNeighbours(evt,geom,
                            pixel.GetPixelId(), &dimCluster);

            if (dimCluster > 4)
            {
                for (Int_t i = 0; i < dimCluster; i++)
                    fData5[fCluster[i]]=fakeData[fCluster[i]];
            }
            if (dimCluster > 3)
            {
                for (Int_t i = 0; i < dimCluster; i++)
                    fData4[fCluster[i]]=fakeData[fCluster[i]];
            }
            if (dimCluster > 2)
            {
                for (Int_t i = 0; i < dimCluster; i++)
                    fData3[fCluster[i]]=fakeData[fCluster[i]];
            }
            if (dimCluster > 1)
            {
                for (Int_t i = 0; i < dimCluster; i++)
                    fData2[fCluster[i]]=fakeData[fCluster[i]];
            }            
        }
    }

}

// --------------------------------------------------------------------------
//
// Function to check the nearest neighbour pixels arrivaltime
//

void MArrivalTime::CheckNeighbours(const MRawEvtData &evt, const MGeomCam &geom,
                                   Short_t idx, Short_t *dimCluster)
{
    Byte_t numNeighbors=geom[idx].GetNumNeighbors();
    fPixelChecked[idx] = kTRUE;

    for (Byte_t i=0x00; i < numNeighbors; i++)
    {
        
        if (!fPixelChecked[geom[idx].GetNeighbor(i)] &&
            fakeData[idx] == fakeData[geom[idx].GetNeighbor(i)])
        {    
            fCluster[*dimCluster]=geom[idx].GetNeighbor(i);    
            *dimCluster++;
            CheckNeighbours(evt, geom,
                            geom[idx].GetNeighbor(i), dimCluster);
        }
    }
}
// --------------------------------------------------------------------------
//
// Returns the arrival time value
//

Bool_t MArrivalTime::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
{
    if (idx<0 || idx>=fData.GetSize())
        return kFALSE;

    switch (type)
    {
    case 0:
    case 1:
        val = fData[idx];
        break;
    case 2:
        val = fData2[idx];
        break;
    case 3:
        val = fData3[idx];
        break;
    case 4:
        val = fData4[idx];
        break;
    case 5:
        val = fData5[idx];
        break;
    }

    return kTRUE;
}

void MArrivalTime::DrawPixelContent(Int_t num) const
{
    *fLog << warn << "MArrivalTime::DrawPixelContent - not available." << endl;
}