source: trunk/Mars/msim/MPhotonEvent.cc@ 9969

/* ======================================================================== *\
!
! *
! * 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"
#include "MCorsikaFormat.h"

ClassImp(MPhotonEvent);

using namespace std;

// ==========================================================================

class MyClonesArray : public TClonesArray
{
public:
    TObject **FirstRef() { return fCont; }

    // --------------------------------------------------------------------------
    //
    // This is an extremly optimized version of ExpandCreateFast. It only resets
    // the marker for the last element (fLast) to n-1 and doen't change anything
    // else. This implicitly assumes that the stored objects don't allocate
    // memory. It does not necessarily mean that the slots after fLast
    // are empty (set to 0). This is what is assumed in the TClonesArray.
    // We also don't call Changed() because it would reset Sorted. If the
    // array was sorted before it is also sorted now. You MUST make sure
    // that you only set n in a range for which valid entries have been
    // created before (e.g. by ExpandCreateFast).
    //
    void FastShrink(Int_t n)
    {
        fLast = n - 1;
    }

    // --------------------------------------------------------------------------
    //
    // This is a optimized (faster) version of Delete which deletes consequtive
    // entries from index idx1 to idx2 (both included) and calls their
    // destructor. Note that there is no range checking done!
    //
    void FastRemove(Int_t idx1, Int_t idx2)
    {
        // Remove object at index idx.

        //if (!BoundsOk("RemoveAt", idx1)) return 0;
        //if (!BoundsOk("RemoveAt", idx2)) return 0;

        Long_t dtoronly = TObject::GetDtorOnly();

        idx1 -= fLowerBound;
        idx2 -= fLowerBound;

        for (TObject **obj=fCont+idx1; obj<=fCont+idx2; obj++)
        {
            if (!*obj)
                continue;

            if ((*obj)->TestBit(kNotDeleted)) {
                // Tell custom operator delete() not to delete space when
                // object fCont[i] is deleted. Only destructors are called
                // for this object.
                TObject::SetDtorOnly(*obj);
                delete *obj;
            }

            *obj = 0;
            // recalculate array size
        }
        TObject::SetDtorOnly((void*)dtoronly);

        if (idx1<=fLast && fLast<=idx2)
        {
            do {
                fLast--;
            } while (fLast >= 0 && fCont[fLast] == 0);
        }

        Changed();
    }


    //void SetSorted() { fSorted = kTRUE; }

    // --------------------------------------------------------------------------
    //
    // This is an optimized version of Sort which doesn't check the
    // IsSortable flag before. It only sorts the entries from 0
    // to GetEntriesFast().
    //
    void UncheckedSort()
    {
        if (fSorted)
            return;

        const Int_t nentries = GetEntriesFast();
        if (nentries <= 0)
            return;

        QSort(GetObjectRef(First()), fKeep->GetObjectRef(fKeep->First()),
              0, TMath::Min(nentries, kMaxInt-fLowerBound));

        fSorted = kTRUE;
    }
};

// ==========================================================================

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

// --------------------------------------------------------------------------
//
// This is an extremly optimized version of ExpandCreateFast. It only resets
// the marker for the last element (fLast) to n-1 and doen't change anything
// else. This has the advantage that the allocated memory is kept but only
// valid entries are written to a file.
//
// If the array was sorted before it is also sorted now. You MUST make sure
// that you only set n in a range for which valid entries have been
// created before (e.g. by ExpandCreateFast).
//
Int_t MPhotonEvent::Shrink(Int_t n)
{
    /*
    // The number of object written by the streamer is defined by
    // GetEntriesFast(), i.e. this number must be shrinked to
    // the real array size. We use RemoveAt instead of ExpandCreate
    // because RemoveAt doesn't free memory. Thus in the next
    // iteration it can be reused and doesn't need to be reallocated.
    // Do not change this. It is optimized for execution speed
    //        for (int i=fData.GetSize()-1; i>=n; i--)
    //          fData.RemoveAt(i);
    const Bool_t sorted = fData.IsSorted();

    MyClonesArray &loc = static_cast<MyClonesArray&>(fData);

    loc.FastRemove(n, fData.GetSize()-1);

    // If it was sorted before it is also sorted now.
    if (sorted)
        loc.SetSorted();
    */

    // fData.ExpandCreateFast(n);  // Just set fLast = n -1

    // Just set fLast = n -1
    static_cast<MyClonesArray&>(fData).FastShrink(n);
    return fData.GetEntriesFast();
}

// --------------------------------------------------------------------------
//
// The resized the array. If it has to be increased in size it is done
// with ExpandCreateFast. If it should be shrinked it is done with
// ExpandCreateFast if n>fData.GetSize()/100 or n==0. This keeps the allocated
// memory and just sets the marker for the last element in the array (fLast).
//
// If the allocated memory grew to huge we reset the allocated memory
// by calling ExpandCreate(n) (frees the allocated storage for the
// objects) and Expand(n) (frees the allocated memory for the list
// of pointers to the objects)
//
// 100 hundred is an arbitrary number taking into account that todays
// computers have a lot of memory and we don't want to free and allocate
// new memory too often.
//
// In priciple there might be more clever methods to handle the memory.
//
void MPhotonEvent::Resize(Int_t n)
{
    if (n==0 || n*100>fData.GetSize())
        fData.ExpandCreateFast(n);  // Just set fLast = n -1
    else
    {
        fData.ExpandCreate(n);      // Free memory of allocated MPhotonData
        fData.Expand(n);            // Free memory of allocated pointers
    }
}

// --------------------------------------------------------------------------
//
// 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 = 0;
    if (n<fData.GetSize() && fData.UncheckedAt(n))
    {
        o=fData.UncheckedAt(n);
        static_cast<MyClonesArray&>(fData).FastShrink(n+1);
    }
    else
    {
        o=fData.New(n);
    }
    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());
}

void MPhotonEvent::Sort(Bool_t force)
{
    if (force)
        fData.UnSort();

    static_cast<MyClonesArray&>(fData).UncheckedSort(); /*Sort(GetEntriesFast())*/
}

// --------------------------------------------------------------------------
//
// 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 fData.GetEntriesFast()==0 ? 0 : static_cast<MPhotonData*>(fData.First());
}

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

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

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

    return rc;
}

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

Double_t MPhotonEvent::GetMeanX() const
{
    const UInt_t n = GetNumPhotons();
    if (n==0)
        return 0;

    Double_t mean = 0;
    for (UInt_t i=0; i<n; i++)
        mean += operator[](i).GetPosX();

    return mean / n;
}

Double_t MPhotonEvent::GetMeanY() const
{
    const UInt_t n = GetNumPhotons();
    if (n==0)
        return 0;

    Double_t mean = 0;
    for (UInt_t i=0; i<n; i++)
        mean += operator[](i).GetPosY();

    return mean / n;
}

Double_t MPhotonEvent::GetMeanT() const
{
    const UInt_t n = GetNumPhotons();
    if (n==0)
        return 0;

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

    return mean / n;
}

void MPhotonEvent::AddXY(Double_t x, Double_t y)
{
    const UInt_t n = GetNumPhotons();

    for (UInt_t i=0; i<n; i++)
    {
        MPhotonData &p = operator[](i);
        p.SetPosition(p.GetPosX()+x, p.GetPosY()+y);
    }
}

void MPhotonEvent::SimWavelength(Float_t wmin, Float_t wmax)
{
    const UInt_t n = GetNumPhotons();

    for (UInt_t i=0; i<n; i++)
        operator[](i).SimWavelength(wmin, wmax);
}

// --------------------------------------------------------------------------
//
// Read the Event section from the file
//
Int_t MPhotonEvent::ReadCorsikaEvt(MCorsikaFormat *fInFormat, Int_t id)
{
    Int_t n = 0;

    // --- old I/O ---
    // Read only + Reflector (no absorption)
    // Muons:   1.06GB/115s =  9.2MB/s (100kEvs)
    // Gammas:  1.57GB/275s =  5.7MB/s (  1kEvs)

    // Read only:
    // Gammas:  1.57GB/158s =  9.9MB/s (  1kEvs)
    // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)

    // --- new I/O ---
    // Read only (don't allocate storage space):
    // Gammas:  1.57GB/143s = 11.0MB/s (  1kEvs)
    // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)

    // Read only in blocks (with storage allocation):
    // Gammas:  1.57GB/28s  =  56MB/s (  1kEvs)
    // Muons:   1.06GB/5.2s = 204MB/s (100kEvs)

    // Read only in blocks (without storage allocation):
    // similar to just copy

    // Copy with cp
    // 1.57GB/ 5s   CPU
    // 1.57GB/28s   REAL
    // 1.06GB/ 3s   CPU
    // 1.06GB/22s   REAL
    Float_t *buffer = 0;

    if (fInFormat->IsEventioFormat())
    {
        while (1)
        {
            const Int_t rc = fInFormat->GetNextEvent(&buffer, id);
            if (rc==kERROR)
                return kERROR;
            if (rc==kFALSE)
                break;

            // Loop over number of photons in bunch
            while (Add(n).FillEventIO(buffer))
                n++;
        }
    }
    else
    {
        while (1)
        {
            const Int_t rc1 = fInFormat->GetNextEvent(&buffer);
            if (rc1==kERROR)
                return kERROR;
            if (rc1==kFALSE)
                break;

            const Int_t rc2 = Add(n).FillCorsika(buffer, id);
            switch (rc2)
            {
            case kCONTINUE:  continue;        // No data in this bunch... skip it.
            case kERROR:     return kERROR;   // Error occured
            //case kFALSE:     return kFALSE;   // End of stream
            }

            // This is a photon we would like to keep later.
            // Increase the counter by one
            n++;
        }
    }

/*

    while (1)
    {
        Float_t buffer[273];
        Float_t * ptr = buffer;


        if (!fInFormat->ReadData(273, buffer))
            return kFALSE;

        if (!memcmp(ptr, "EVTE", 4))
            {

            fInFormat->UnreadLastData();
            break;
            }

        Float_t *end = ptr + 273;

        // Loop over all sub-blocks (photons) in the block and if
        // they contain valid data add them to the array
        while (ptr<end)
        {
            // Get/Add the n-th entry from the array and
            // fill it with the current 7 floats
            const Int_t rc = Add(n).FillCorsika(ptr);
            ptr += 7;

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

            // This is a photon we would like to keep later.
            // Increase the counter by one
            n++;
        }
    }

*/
    Resize(n);
    fData.UnSort();

    SetReadyToSave();

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

}

Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin, Int_t i)
{
    Int_t n = 0;

    // --- old I/O ---
    // Read only + Reflector (no absorption)
    // Muons:   1.06GB/115s =  9.2MB/s (100kEvs)
    // Gammas:  1.57GB/275s =  5.7MB/s (  1kEvs)

    // Read only:
    // Gammas:  1.57GB/158s =  9.9MB/s (  1kEvs)
    // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)

    // --- new I/O ---
    // Read only (don't allocate storage space):
    // Gammas:  1.57GB/143s = 11.0MB/s (  1kEvs)
    // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)

    // Read only in blocks (with storage allocation):
    // Gammas:  1.57GB/28s  =  56MB/s (  1kEvs)
    // Muons:   1.06GB/5.2s = 204MB/s (100kEvs)

    // Read only in blocks (without storage allocation):
    // similar to just copy

    // Copy with cp
    // 1.57GB/ 5s   CPU
    // 1.57GB/28s   REAL
    // 1.06GB/ 3s   CPU
    // 1.06GB/22s   REAL

    while (1)
    {
        // Stprage for one block
        char c[273*4];

        // Read the first four byte to check whether the next block
        // doen't belong to the event anymore
        fin.read(c, 4);
        if (!fin)
            return kFALSE;

        // Check if the event is finished
        if (!memcmp(c, "EVTE", 4))
            break;

        // Now read the rest of the data
        fin.read(c+4, 272*4);

        Float_t *ptr = reinterpret_cast<Float_t*>(c);
        Float_t *end = ptr + 273;

        // Loop over all sub-blocks (photons) in the block and if
        // they contain valid data add them to the array
        while (ptr<end)
        {
            // Get/Add the n-th entry from the array and
            // fill it with the current 7 floats
            const Int_t rc = Add(n).FillCorsika(ptr, i);
            ptr += 7;

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

            // This is a photon we would like to keep later.
            // Increase the counter by one
            n++;
        }
    }
/*
    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++;
    }
  */

    Resize(n);
    fData.UnSort();

    SetReadyToSave();

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

// --------------------------------------------------------------------------
/*
Int_t MPhotonEvent::ReadRflEvt(std::istream &fin, Int_t i)
{
    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, i);

        // 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
{
    // This is much faster than looping over all entries and discarding
    // the empty ones
    const UInt_t n = GetNumPhotons();
    for (UInt_t i=0; i<n; i++)
        operator[](i).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();
    }
}