source: trunk/MagicSoft/Mars/manalysis/MCerPhotEvt.cc@ 4904

/* ======================================================================== *\
!
! *
! * 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): Thomas Bretz, 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
!   Author(s): Harald Kornmayer, 1/2001
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//            
// MCerPhotEvt
//
// NOTE: This container is NOT ment for I/O. Write it to a file on your own
//       risk!
//
// Class Version 3:
// ----------------
//  - added fNumIslands
//            
// Class Version 2:
// ----------------
//  - added fLut to accelerate GetPixById a lot
//            
// Class Version 1:
// ----------------
//  - first version
//
/////////////////////////////////////////////////////////////////////////////
#include "MCerPhotEvt.h"

#include <math.h>
#include <limits.h>
#include <fstream>

#include <TArrayD.h>
#include <TCanvas.h>

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

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

ClassImp(MCerPhotEvt);
ClassImp(MCerPhotEvtIter);

using namespace std;

// --------------------------------------------------------------------------
//
// Creates a MCerPhotPix object for each pixel in the event
//
MCerPhotEvt::MCerPhotEvt(const char *name, const char *title) : fNumPixels(0)
{
    fName  = name  ? name  : "MCerPhotEvt";
    fTitle = title ? title : "(Number of Photon)-Event Information";

    fPixels = new TClonesArray("MCerPhotPix", 0);
}

// --------------------------------------------------------------------------
//
// This is not yet implemented like it should.
//
/*
void MCerPhotEvt::Draw(Option_t* option) 
{
    //
    //   FIXME!!! Here the Draw function of the CamDisplay
    //   should be called to add the CamDisplay to the Pad.
    //   The drawing should be done in MCamDisplay::Paint
    //

    //    MGeomCam *geom = fType ? new MGeomCamMagic : new MGeomCamCT1;
    //    MCamDisplay *disp = new MCamDisplay(geom);
    //    delete geom;
    //    disp->DrawPhotNum(this);
}
*/

// --------------------------------------------------------------------------
//
// reset counter and delete netries in list.
//
void MCerPhotEvt::Reset()
{
    fNumPixels  =  0;
    fMaxIndex   = -1;
    fNumIslands = -1;
    fLut.Set(0);
    // fPixels->Delete();
}

void MCerPhotEvt::FixSize()
{
    fLut.Set(fMaxIndex+1);

    if (fPixels->GetEntriesFast() == (Int_t)fNumPixels)
        return;

    fPixels->ExpandCreateFast(fNumPixels);
}

// --------------------------------------------------------------------------
//
//  Dump the cerenkov photon event to *fLog
//
void MCerPhotEvt::Print(Option_t *) const
{
    const Int_t entries = fPixels->GetEntries();

    *fLog << GetDescriptor() << dec << endl;
    *fLog << " Number of Pixels: " << fNumPixels << "(" << entries << ")" << endl;

    for (Int_t i=0; i<entries; i++ )
        (*this)[i].Print();
}

// --------------------------------------------------------------------------
//
// Checks if in the pixel list is an entry with pixel id
//
Bool_t MCerPhotEvt::IsPixelExisting(Int_t id) const
{
    return GetPixById(id) ? kTRUE : kFALSE;
} 

// --------------------------------------------------------------------------
//
//   Checks if in the pixel list is an entry with pixel id
//
Bool_t MCerPhotEvt::IsPixelUsed(Int_t id) const
{
    const MCerPhotPix *pix = GetPixById(id);
    return pix ? pix->IsPixelUsed() : kFALSE;
} 

// --------------------------------------------------------------------------
//
//   Checks if in the pixel list is an entry with pixel id
//
Bool_t MCerPhotEvt::IsPixelCore(Int_t id) const
{
    const MCerPhotPix *pix = GetPixById(id);
    return pix ? pix->IsPixelCore() : kFALSE;
} 

// --------------------------------------------------------------------------
//
// get the minimum number of photons  of all valid pixels in the list
// If you specify a geometry the number of photons is weighted with the
// area of the pixel
//
Float_t MCerPhotEvt::GetNumPhotonsMin(const MGeomCam *geom) const
{
    if (fNumPixels <= 0)
        return -5.;

    const UInt_t n = geom->GetNumPixels();

    Float_t minval = FLT_MAX;

    for (UInt_t i=0; i<fNumPixels; i++)
    {
        const MCerPhotPix &pix = (*this)[i];
        if (!pix.IsPixelUsed())
            continue;

        const UInt_t id = pix.GetPixId();
        if (id<0 || id>=n)
            continue;

        Float_t testval = pix.GetNumPhotons();

        if (geom)
            testval *= geom->GetPixRatio(id);

        if (testval < minval)
            minval = testval;
    }

    return minval;
}

// --------------------------------------------------------------------------
//
// get the maximum number of photons of all valid pixels in the list
// If you specify a geometry the number of photons is weighted with the
// area of the pixel
//
Float_t MCerPhotEvt::GetNumPhotonsMax(const MGeomCam *geom) const
{
    if (fNumPixels <= 0)
        return 50.;

    const UInt_t n = geom->GetNumPixels();

    Float_t maxval = -FLT_MAX;

    for (UInt_t i=0; i<fNumPixels; i++)
    {
        const MCerPhotPix &pix = (*this)[i];
        if (!pix.IsPixelUsed())
            continue;

        const UInt_t id = pix.GetPixId();
        if (id<0 || id>=n)
            continue;

        Float_t testval = pix.GetNumPhotons();
        if (geom)
            testval *= geom->GetPixRatio(id);

        if (testval > maxval)
            maxval = testval;
    }
    return maxval;
}

// --------------------------------------------------------------------------
//
// get the minimum ratio of photons/error
//
Float_t MCerPhotEvt::GetRatioMin(const MGeomCam *geom) const
{
    if (fNumPixels <= 0)
        return -5.;

    Float_t minval = FLT_MAX;

    for (UInt_t i=0; i<fNumPixels; i++)
    {
        const MCerPhotPix &pix = (*this)[i];
        if (!pix.IsPixelUsed())
            continue;

        Float_t testval = pix.GetNumPhotons()/pix.GetErrorPhot();
        if (geom)
            testval *= TMath::Sqrt(geom->GetPixRatio(pix.GetPixId()));

        if (testval < minval)
            minval = testval;
    }

    return minval;
}

// --------------------------------------------------------------------------
//
// get the maximum ratio of photons/error
//
Float_t MCerPhotEvt::GetRatioMax(const MGeomCam *geom) const
{
    if (fNumPixels <= 0)
        return -5.;

    Float_t maxval = -FLT_MAX;

    for (UInt_t i=0; i<fNumPixels; i++)
    {
        const MCerPhotPix &pix = (*this)[i];
        if (!pix.IsPixelUsed())
            continue;

        Float_t testval = pix.GetNumPhotons()/pix.GetErrorPhot();
        if (geom)
            testval *= TMath::Sqrt(geom->GetPixRatio(pix.GetPixId()));

        if (testval > maxval)
            maxval = testval;
    }

    return maxval;
}

// --------------------------------------------------------------------------
//
// get the minimum of error
// If you specify a geometry the number of photons is weighted with the
// area of the pixel
//
Float_t MCerPhotEvt::GetErrorPhotMin(const MGeomCam *geom) const
{
    if (fNumPixels <= 0)
        return 50.;

    Float_t minval = FLT_MAX;

    for (UInt_t i=0; i<fNumPixels; i++)
    {
        const MCerPhotPix &pix = (*this)[i];
        if (!pix.IsPixelUsed())
            continue;

        Float_t testval = pix.GetErrorPhot();

        if (geom)
            testval *= geom->GetPixRatio(pix.GetPixId());

        if (testval < minval)
            minval = testval;
    }
    return minval;
}

// --------------------------------------------------------------------------
//
// get the maximum ratio of photons/error
// If you specify a geometry the number of photons is weighted with the
// area of the pixel
//
Float_t MCerPhotEvt::GetErrorPhotMax(const MGeomCam *geom) const
{
    if (fNumPixels <= 0)
        return 50.;

    Float_t maxval = -FLT_MAX;

    for (UInt_t i=0; i<fNumPixels; i++)
    {
        const MCerPhotPix &pix = (*this)[i];
        if (!pix.IsPixelUsed())
            continue;

        Float_t testval = pix.GetErrorPhot();

        if (geom)
            testval *= geom->GetPixRatio(pix.GetPixId());

        if (testval > maxval)
            maxval = testval;
    }
    return maxval;
}

void MCerPhotEvt::RemoveUnusedPixels()
{
    // Create iterator
    TIter Next(fPixels);
    MCerPhotPix *pix = NULL;

    fMaxIndex = -1;

    // Remove all unused pixels from list, calculate new fMaxIndex
    while ((pix=(MCerPhotPix*)Next()))
    {
        if (!pix->IsPixelUsed())
            fPixels->Remove(pix);
        else
            fMaxIndex = TMath::Max(fMaxIndex, pix->GetPixId());
    }

    // Crompress array
    fPixels->Compress();

    // Get new number of entries in array
    fNumPixels=fPixels->GetEntriesFast();

    // Rebuild lookup table
    RebuildLut();
}

// --------------------------------------------------------------------------
//
// Return a pointer to the pixel with the requested idx. NULL if it doesn't
// exist. The Look-up-table fLut is used. If its size is zero (according
// to Rene this will happen if an old class object is loaded) we still
// try to search in the array.
//
MCerPhotPix *MCerPhotEvt::GetPixById(Int_t idx) const
{
    if (idx<0)
        return 0;

    if (fLut.GetSize()>0)
    {
        if (idx>=fLut.GetSize())
            return 0;
        return fLut[idx]<0 ? 0 : (MCerPhotPix*)(fPixels->UncheckedAt(fLut[idx]));
    }

    TIter Next(fPixels);
    MCerPhotPix *pix = NULL;

    while ((pix=(MCerPhotPix*)Next()))
        if (pix->GetPixId()==idx)
            return pix;

    return NULL;
}

MCerPhotPix *MCerPhotEvt::AddPixel(Int_t idx, Float_t nph, Float_t er)
    {
        //
        // If this is too slow or takes to much space we might use
        // MGeomApply and an InitSize member function instead.
        //
        if (idx>=fLut.GetSize())
        {
            const Int_t n = fLut.GetSize();
            fLut.Set(idx*2+1); //idx+1 is slower than idx*2+1
            for (int i=n; i<idx*2+1; i++)
                fLut[i] = -1;
        }

        fLut[idx] = fNumPixels;
        if (idx>fMaxIndex)
            fMaxIndex=idx;

        return new ((*fPixels)[fNumPixels++]) MCerPhotPix(idx, nph, er);
    }

// --------------------------------------------------------------------------
//
// This function recursively finds all pixels of one island and assigns
// the number num as island number to the pixel.
//
//  1) Check whether a pixel with the index idx exists, is unused
//     and has not yet a island number assigned.
//  2) Assign the island number num to the pixel
//  3) Loop over all its neighbors taken from the geometry geom. For all
//     neighbors recursively call this function (CalcIsland)
//  4) Sum the size of the pixel and all neighbors newly assigned
//     (by CalcIsland) to this island
//
// Returns the sum of the pixel size.
//
Double_t MCerPhotEvt::CalcIsland(const MGeomCam &geom, Int_t idx, Int_t num)
{
    // Try to get the pixel information of a pixel with this index
    MCerPhotPix *pix = GetPixById(idx);

    // If a pixel with this index is not existing... do nothing.
    if (!pix)
        return 0;

    // If an island number was already assigned to this pixel... do nothing.
    if (pix->GetIdxIsland()>=0)
        return 0;

    // If the pixel is an unused pixel... do nothing.
    if (!pix->IsPixelUsed())
        return 0;

    // Assign the new island number num to this used pixel
    pix->SetIdxIsland(num);

    // Get the geometry information (neighbors) of this pixel
    const MGeomPix &gpix = geom[idx];

    // Get the size of this pixel
    Double_t size = pix->GetNumPhotons();

    // Now do the same with all its neighbors and sum the
    // sizes which they correspond to
    const Int_t n = gpix.GetNumNeighbors();
    for (int i=0; i<n; i++)
        size += CalcIsland(geom, gpix.GetNeighbor(i), num);

    // return size of this (sub)cluster
    return size;
}

// --------------------------------------------------------------------------
//
// Each pixel which is maked as used is assigned an island number
// (starting from 0). A pixel without an island number assigned
// has island number -1.
//
// The index 0 corresponds to the island with the highest size (sum
// of GetNumPhotons() in island). The size is decreasing with
// increasing indices.
//
// The information about pixel neighbory is taken from the geometry
// MGeomCam geom;
//
// You can access this island number of a pixel with a call to
// MCerPhotPix->GetIdxIsland. The total number of islands available
// can be accessed with MCerPhotEvt->GetNumIslands.
//
// CalcIslands returns the number of islands found. If an error occurs,
// eg the geometry has less pixels than the highest index stored, -1 is
// returned.
//
Int_t MCerPhotEvt::CalcIslands(const MGeomCam &geom)
{
    if (fMaxIndex<0 || fNumPixels<=0)
    {
        *fLog << err << "ERROR - MCerPhotEvt doesn't contain pixels!" << endl;
        fNumIslands = 0;
        return -1;
    }

    if ((UInt_t)fMaxIndex>=geom.GetNumPixels())
    {
        *fLog << err << "ERROR - MCerPhotEvt::CalcIslands: Size mismatch - geometry too small!" << endl;
        return -1;
    }

    // Create a list to hold the sizes of the islands (The maximum
    // number of islands possible is rougly fNumPixels/4)
    TArrayD size(fNumPixels/3);

    // Calculate Islands
    Int_t n=0;

    // We could loop over all indices which looks more straight
    // forward but should be a lot slower (assuming zero supression)
    MCerPhotPix *pix=0;

    TIter Next(*this);
    while ((pix=static_cast<MCerPhotPix*>(Next())))
    {
        // This 'if' is necessary (although it is done in GetIsland, too)
        // because otherwise the counter (n) would be wrong.
        // So only 'start' a new island for used pixels (selected by
        // using the Iterator) which do not yet belong to another island.
        if (pix->GetIdxIsland()<0)
        {
            Double_t sz = CalcIsland(geom, pix->GetPixId(), n);
            size[n++] = sz;
        }
    }

    // Create an array holding the indices
    TArrayI idx(n);

    // Sort the sizes descending
    TMath::Sort(n, size.GetArray(), idx.GetArray(), kTRUE);

    // Replace island numbers by size indices -- After this
    // islands indices are sorted by the island size
    Next.Reset();
    while ((pix=static_cast<MCerPhotPix*>(Next())))
    {
        const Short_t i = pix->GetIdxIsland();

        // Find new index
        Short_t j;
        for (j=0; j<n; j++)
            if (idx[j]==i)
                break;

        pix->SetIdxIsland(j==n ? -1 : j);
    }

    // Now assign number of islands found
    fNumIslands = n;

    // return number of island
    return fNumIslands;
}

// --------------------------------------------------------------------------
//
// Returns, depending on the type flag:
//
//  0: Number of Photons*PixRatio
//  1: Error*sqrt(PixRatio)
//  2: Cleaning level = Num Photons*sqrt(PixRatio)/Error
//  3: Number of Photons
//  4: Error
//  5: Island index
//
Bool_t MCerPhotEvt::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
{
    MCerPhotPix *pix = GetPixById(idx);

    if (!pix || !pix->IsPixelUsed())
        return kFALSE;

    const Double_t ratio = cam.GetPixRatio(idx);

    switch (type)
    {
    case 1:
        val = pix->GetErrorPhot()*TMath::Sqrt(ratio);
        return kTRUE;
    case 2:
        if (pix->GetErrorPhot()<=0)
            return kFALSE;
        val = pix->GetNumPhotons()*TMath::Sqrt(ratio)/pix->GetErrorPhot();
        return kTRUE;
    case 3:
        val = pix->GetNumPhotons();
        break;
    case 4:
        val = pix->GetErrorPhot();
        break;
    case 5:
        val = pix->GetIdxIsland();
        break;
    default:
        val = pix->GetNumPhotons()*ratio;
        return kTRUE;
    }
    return kTRUE;
}

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

TObject *MCerPhotEvtIter::Next()
{
    if (!fUsedOnly)
        return TObjArrayIter::Next();

    MCerPhotPix *pix;
    while ((pix = (MCerPhotPix*)TObjArrayIter::Next()))
        if (pix->IsPixelUsed())
            return pix;
    return pix;
}