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

/* ======================================================================== *\
!
! *
! * 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";
}

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

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

    fData.Set(n);
    fData.Reset();

    MRawEvtPixelIter pixel((MRawEvtData*)&evt);
   
    while ( pixel.Next() )
    {
        const UInt_t idx = pixel.GetPixelId();

//If pixel has saturated and hasn't lo gains we return -1

        if (pixel.GetMaxHiGainSample() == 0xff)
        {   
            if (!pixel.HasLoGain())        
            {
                fData[idx]=-1.0;
            }
            else //Calculate time using Lo Gains
            {
                Byte_t *ptr = pixel.GetLoGainSamples();

                Int_t nSlice = evt.GetNumLoGainSamples();
//Some casts are needed because TSpline5 constructor accepts only Double_t values
                Double_t ptr2[nSlice];
                for (Int_t i = 0; i < nSlice; i++)
                    ptr2[i]=(Double_t)ptr[i];
                TSpline5 *spline = new TSpline5("spline",(Double_t) 0,(Double_t)(nSlice - 1),ptr2,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;
            }
        }
        else //Calculate time using Hi Gains
        {
            Byte_t *ptr = pixel.GetHiGainSamples();

            Int_t nSlice = evt.GetNumHiGainSamples();
//Some casts are needed because TSpline5 constructor accepts only Double_t values
            Double_t  ptr2[nSlice];
            for (Int_t i = 0; i < nSlice; i++)
                ptr2[i]=(Double_t)ptr[i];
            TSpline5 *spline = new TSpline5("spline",(Double_t) 0,(Double_t)(nSlice - 1),ptr2,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);
    fData3.Set(n);
    fData3.Reset(-1);
    fData4.Set(n);
    fData4.Reset(-1);
    fData5.Set(n);
    fData5.Reset(-1);

    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

// 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;
}