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

/* ======================================================================== *\
!
! *
! * 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): Qi Zhe,       06/2007 <mailto:qizhe@astro.uni-wuerzburg.de>
!
!   Copyright: 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);
//
//
//   Version 1:
//   ----------
//    - First implementation
//
/////////////////////////////////////////////////////////////////////////////
#include "MPhotonEvent.h"

#include <fstream>
#include <iostream>

#include <TMarker.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));
}

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

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

void MPhotonEvent::Print(Option_t *) const
{
    fData.Print();
}

/*
// --------------------------------------------------------------------------
//
// Fills all photon distances scaled with scale (default=1) into a TH1
//
void MPhotonEvent::FillRad(TH1 &hist, Float_t scale) const
{
    MPhotonData *ph=NULL;

    TIter Next(&fData);
    while ((ph=(MPhotonData*)Next()))
        ph->FillRad(hist, scale);
}

// --------------------------------------------------------------------------
//
// Fills all photon distances scaled with scale (default=1) versus x
// into a TH2
//
void MPhotonEvent::FillRad(TH2 &hist, Double_t x, Float_t scale) const
{
    MPhotonData *ph=NULL;

    TIter Next(&fData);
    while ((ph=(MPhotonData*)Next()))
        ph->FillRad(hist, x, scale);
}

// --------------------------------------------------------------------------
//
// Fills all photons (x,y) scaled with scale (default=1) into a TH2.
//
void MPhotonEvent::Fill(TH2 &hist, Float_t scale) const
{
    MPhotonData *ph=NULL;

    TIter Next(&fData);
    while ((ph=(MPhotonData*)Next()))
        ph->Fill(hist, scale);
}

// --------------------------------------------------------------------------
//
// Fills all photons (x,y) scaled with scale (default=1) versus z into a TH3
//
void MPhotonEvent::Fill(TH3 &hist, Double_t z, Float_t scale) const
{
    MPhotonData *ph=NULL;

    TIter Next(&fData);
    while ((ph=(MPhotonData*)Next()))
        ph->Fill(hist, z, scale);
}

void MPhotonEvent::DrawPixelContent(Int_t num) const
{
}

#include "MHexagon.h"
#include "MGeomCam.h"
#include <TMarker.h>

// ------------------------------------------------------------------------
//
// Fill a reflector event. Sums all pixels in each pixel as the
// pixel contents.
//
// WARNING: Due to the estimation in DistanceToPrimitive and the
//          calculation in pixels instead of x, y this is only a
//          rough estimate.
//
Bool_t MPhotonEvent::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
{
    //
    // sum the photons content in each pixel
    //
    val = 0;

    MHexagon hex(cam[idx]);

    MPhotonData *ph=NULL;

    TIter Next(&fData);

    switch (type)
    {

    case 1: // time
        while ((ph=(MPhotonData*)Next()))
            if (hex.DistanceToPrimitive(ph->GetPosY(), ph->GetPosX())<=0)
                val += ph->GetTime();
        return val/fData.GetEntriesFast()>0;

    default: // signal
        while ((ph=(MPhotonData*)Next()))
            if (hex.DistanceToPrimitive(ph->GetPosY()*10, ph->GetPosX()*10)<=0)
                val += cam.GetPixRatio(idx);
        return val>0;
    }

    return kFALSE;

}
*/

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