source: trunk/MagicSoft/Mars/msim/MPhotonEvent.cc@ 9347

/* ======================================================================== *\
!
! *
! * This file is part of CheObs, the Modular 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 appears 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): Qi Zhe,       06/2007 <mailto:qizhe@astro.uni-wuerzburg.de>
!
!   Copyright: CheObs Software Development, 2000-2009
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MPhotonEvent
//
// Storage container to store photon collections
//
// The class is designed to be extremely fast which is important taking into
// account the extremely high number of photons. This has some impacts on
// its handling.
//
// The list has to be kept consistent, i.e. without holes.
//
// There are two ways to achieve this:
//
//   a) Use RemoveAt to remove an entry somewhere
//   b) Compress() the TClonesArray afterwards
//
// Compress is not the fastes, so there is an easier way.
//
//   a) When you loop over the list and want to remove an entry copy all
//      following entry backward in the list, so that the hole will
//      be created at its end.
//   b) Call Shrink(n) with n the number of valid entries in your list.
//
// To loop over the TClonesArray you can use a TIter which has some
// unnecessary overhead and therefore is slower than necessary.
//
// Since the list is kept consistent you can use a simple loop saving
// a lot of CPU time taking into account the high number of calls to
// TObjArrayIter::Next which you would create.
//
// Here is an example (how to remove every second entry)
//
//  ---------------------------------------------------------------------
//
//    Int_t cnt = 0;
//
//    const Int_t num = event->GetNumPhotons();
//    for (Int_t idx=0; idx<num; idx++)
//    {
//        if (idx%2==0)
//           continue;
//
//        MPhotonData *dat = (*event)[idx];
//
//        (*event)[cnt++] = *dat;
//     }
//
//     event->Shrink(cnt);
//
//  ---------------------------------- or -------------------------------
//
//    TClonesArray &arr = MPhotonEvent->GetArray();
//
//    Int_t cnt = 0;
//
//    const Int_t num = arr.GetEntriesFast();
//    for (Int_t idx=0; idx<num; idx++)
//    {
//        if (idx%2==0)
//           continue;
//
//        MPhotonData *dat = static_cast<MPhotonData*>(arr.UncheckedAt(idx));
//
//        *static_cast<MPhotonData*>(arr.UncheckedAt(cnt++)) = *dat;
//     }
//
//     MPhotonEvent->Shrink(cnt);
//
//  ---------------------------------------------------------------------
//
// The flag for a sorted array is for speed reasons not in all conditions
// maintained automatically. Especially Add() doesn't reset it.
//
// So be sure that if you want to sort your array it is really sorted.
//
//
//   Version 1:
//   ----------
//    - First implementation
//
/////////////////////////////////////////////////////////////////////////////
#include "MPhotonEvent.h"

#include <fstream>
#include <iostream>

#include <TMarker.h>

#include <MMath.h>

#include "MArrayF.h"

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

#include "MPhotonData.h"

ClassImp(MPhotonEvent);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor. It initializes all arrays with zero size.
//
MPhotonEvent::MPhotonEvent(const char *name, const char *title)
    : fData("MPhotonData", 1)
{
    fName  = name  ? name  : "MPhotonEvent";
    fTitle = title ? title : "Corsika Event Data Information";

    fData.SetBit(TClonesArray::kForgetBits);
    fData.BypassStreamer(kFALSE);
}

/*
const char *MPhotonEvent::GetClassName() const
{
    return static_cast<TObject*>(fData.GetClass())->GetName();
}
*/

// --------------------------------------------------------------------------
//
// If n is smaller than the current allocated array size a reference to
// the n-th entry is returned, otherwise an entry at n is created
// calling the default constructor. Note, that there is no range check
// but it is not recommended to call this function with
// n>fData.GetSize()
//
MPhotonData &MPhotonEvent::Add(Int_t n)
{
    // Do not modify this. It is optimized for execution
    // speed and flexibility!
    TObject *o = n<fData.GetSize() && fData.UncheckedAt(n) ? fData.UncheckedAt(n) : fData.New(n);

    // Now we read a single cherenkov bunch
    return static_cast<MPhotonData&>(*o);
}

// --------------------------------------------------------------------------
//
// Add a new photon (MPhtonData) at the end of the array.
// In this case the default constructor of MPhotonData is called.
//
// A reference to the new object is returned.
//
MPhotonData &MPhotonEvent::Add()
{
    return Add(GetNumPhotons());
}

// --------------------------------------------------------------------------
//
// Get the i-th photon from the array. Not, for speed reasons there is no
// range check so you are responsible that you do not excess the number
// of photons (GetNumPhotons)
//
MPhotonData &MPhotonEvent::operator[](UInt_t idx)
{
    return *static_cast<MPhotonData*>(fData.UncheckedAt(idx));
}

// --------------------------------------------------------------------------
//
// Get the i-th photon from the array. Not, for speed reasons there is no
// range check so you are responsible that you do not excess the number
// of photons (GetNumPhotons)
//
const MPhotonData &MPhotonEvent::operator[](UInt_t idx) const
{
    return *static_cast<MPhotonData*>(fData.UncheckedAt(idx));
}

// --------------------------------------------------------------------------
//
// Return a pointer to the first photon if available.
//
MPhotonData *MPhotonEvent::GetFirst() const
{
    return static_cast<MPhotonData*>(fData.First());
}

// --------------------------------------------------------------------------
//
// Return a pointer to the last photon if available.
//
MPhotonData *MPhotonEvent::GetLast() const
{
    return static_cast<MPhotonData*>(fData.Last());
}

// --------------------------------------------------------------------------
//
// Return the number of "external" photons, i.e. which are not NightSky
//
Int_t MPhotonEvent::GetNumExternal() const
{
    Int_t n=0;

    for (int i=0; i<GetNumPhotons(); i++)
        if ((*this)[i].GetPrimary()!=MMcEvtBasic::kNightSky)
            n++;

    return n;
}

// --------------------------------------------------------------------------
//
// Return time of first photon, 0 if none in array.
// Note: If you want this to be the earliest make sure that the array
// is properly sorted.
//
Float_t MPhotonEvent::GetTimeFirst() const
{
    const MPhotonData *dat=GetFirst();
    return dat ? dat->GetTime() : 0;
}

// --------------------------------------------------------------------------
//
// Return time of first photon, 0 if none in array.
// Note: If you want this to be the latest make sure that the array
// is properly sorted.
//
Float_t MPhotonEvent::GetTimeLast() const
{
    const MPhotonData *dat=GetLast();
    return dat ? dat->GetTime() : 0;
}

// --------------------------------------------------------------------------
//
// Return the median devian from the median of all arrival times.
// The median deviation is calculated using MMath::MedianDev.
// It is the half width in which one sigma (~68%) of all times are
// contained around the median.
//
Double_t MPhotonEvent::GetTimeMedianDev() const
{
    const UInt_t n = GetNumPhotons();

    MArrayF arr(n);
    for (UInt_t i=0; i<n; i++)
        arr[i] = operator[](i).GetTime();

    return MMath::MedianDev(n, arr.GetArray()/*, Double_t &med*/);
}

// --------------------------------------------------------------------------
//
// Read the Event section from the file
//
Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin)
{
    Int_t n = 0;

    while (1)
    {
        // Check the first four bytes
        char c[4];
        fin.read(c, 4);

        // End of stream
        if (!fin)
            return kFALSE;

        // Check if we found the end of the event
        if (!memcmp(c, "EVTE", 4))
            break;

        // The first for byte contained data already --> go back
        fin.seekg(-4, ios::cur);

        // Do not modify this. It is optimized for execution
        // speed and flexibility!
        MPhotonData &ph = Add(n);
        // It checks how many entries the lookup table has. If it has enough
        // entries and the entry was already allocated, we can re-use it,
        // otherwise we have to allocate it.

        // Now we read a single cherenkov bunch. Note that for speed reason we have not
        // called the constructor if the event was already constructed (virtual table
        // set), consequently we must make sure that ReadCorsikaEvent does reset
        // all data mebers no matter whether they are read or not.
        const Int_t rc = ph.ReadCorsikaEvt(fin);

        // Evaluate result from reading event
        switch (rc)
        {
        case kCONTINUE:  continue;        // No data in this bunch... skip it.
        case kFALSE:     return kFALSE;   // End of stream
        case kERROR:     return kERROR;   // Error occured
        }

        // FIXME: If fNumPhotons!=1 add the photon more than once

        // Now increase the number of entries which are kept,
        // i.e. keep this photon(s)
        n++;
    }

    Shrink(n);
    fData.UnSort();

    SetReadyToSave();

    //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
    return kTRUE;
}

// --------------------------------------------------------------------------
//
Int_t MPhotonEvent::ReadRflEvt(std::istream &fin)
{
    Int_t n = 0;

    while (1)
    {
        // Check the first four bytes
        char c[13];
        fin.read(c, 13);

        // End of stream
        if (!fin)
            return kFALSE;

        // Check if we found the end of the event
        if (!memcmp(c, "\nEND---EVENT\n", 13))
            break;

        // The first for byte contained data already --> go back
        fin.seekg(-13, ios::cur);

        // Do not modify this. It is optimized for execution
        // speed and flexibility!
        //TObject *o = n<fData.GetSize() && fData.UncheckedAt(n) ? fData.UncheckedAt(n) : fData.New(n);

        // Now we read a single cherenkov bunch
        //const Int_t rc = static_cast<MPhotonData*>(o)->ReadRflEvt(fin);
        const Int_t rc = Add(n).ReadRflEvt(fin);

        // Evaluate result from reading event
        switch (rc)
        {
        case kCONTINUE:  continue;        // No data in this bunch... skip it.
        case kFALSE:     return kFALSE;   // End of stream
        case kERROR:     return kERROR;   // Error occured
        }

        // Now increase the number of entries which are kept,
        // i.e. keep this photon(s)
        n++;
    }

    Shrink(n);

    SetReadyToSave();

    //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Print the array
//
void MPhotonEvent::Print(Option_t *) const
{
    fData.Print();
}

// ------------------------------------------------------------------------
//
// You can call Draw() to add the photons to the current pad.
// The photons are painted each tim ethe pad is updated.
// Make sure that you use the right (world) coordinate system,
// like created, eg. by the MHCamera histogram.
//
void MPhotonEvent::Paint(Option_t *)
{
    MPhotonData *ph=NULL;

    TMarker m;
    m.SetMarkerStyle(kFullDotMedium); // Gtypes.h

    TIter Next(&fData);
    while ((ph=(MPhotonData*)Next()))
    {
        m.SetX(ph->GetPosY()*10);  // north
        m.SetY(ph->GetPosX()*10);  // east
        m.Paint();
    }
}