source: trunk/MagicSoft/Mars/msimreflector/MReflector.cc@ 9307

/* ======================================================================== *\
!
! *
! * 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,  1/2009 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: CheObs Software Development, 2000-2009
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
//  MReflector
//
//////////////////////////////////////////////////////////////////////////////
#include "MReflector.h"

#include <fstream>
#include <errno.h>

#include <stdlib.h> // atof (Ubuntu 8.10)

#include <TClass.h>
#include <TSystem.h>

#include "MQuaternion.h"
#include "MMirror.h"

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

ClassImp(MReflector);

using namespace std;

// --------------------------------------------------------------------------
//
// Default constructor
//
MReflector::MReflector(const char *name, const char *title)
{
    fName  = name  ? name  : "MReflector";
    fTitle = title ? title : "Parameter container storing a collection of several mirrors (reflector)";

    fMirrors.SetOwner();
}

// --------------------------------------------------------------------------
//
// Calculate the maximum radius of th ereflector. This is not meant as
// a precise number but as a rough estimate e.g. to bin a histogram.
//
void MReflector::InitMaxR()
{
    fMaxR = 0;

    TIter Next(&fMirrors);
    MMirror *m = 0;
    while ((m=static_cast<MMirror*>(Next())))
    {
        // Take into account the maximum incident angle 8eg 10deg) and
        // the theta-angle of the mirror and the z-distance.
        const Double_t r = m->GetDist()+1.5*m->GetMaxR();
        if (r > fMaxR)
            fMaxR = r;
    }
}

// --------------------------------------------------------------------------
//
// Get the pointer to the first mirror. This is a very dangerous way of
// access, but the fastest possible. because it is the most often called
// function in ExecuteReflector we have to have a very fast access.
//
const MMirror **MReflector::GetFirstPtr() const
{
    return (const MMirror**)fMirrors.GetObjectRef(0);
}

// --------------------------------------------------------------------------
//
// Get number of mirrors. There should be no holes in the array!
//
const UInt_t MReflector::GetNumMirrors() const
{
    return fMirrors.GetEntriesFast();
}

// --------------------------------------------------------------------------
//
// Check with a rough estimate whether a photon can hit the reflector.
//
Bool_t MReflector::CanHit(const MQuaternion &p) const
{
    // p is given in the reflectory coordinate frame. This is meant as a
    // fast check without lengthy calculations to omit all photons which
    // cannot hit the reflector at all
    return p.R2()<fMaxR*fMaxR;
}

// --------------------------------------------------------------------------
//
// Read a reflector setup from a file. This needs improvemtn.
// FIXME: Documentation missing!
//
Bool_t MReflector::ReadFile(TString fname)
{
    SetTitle(fname);
    fMirrors.Delete();

    gSystem->ExpandPathName(fname);

    ifstream fin(fname);
    if (!fin)
    {
        *fLog << err << "Cannot open file " << fname << ": ";
        *fLog << (errno!=0?strerror(errno):"Insufficient memory for decompression") << endl;
        return kFALSE;
    }

/*
    Int_t idx[964];
    Int_t srt[964];
    for (int i=0; i<964; i++)
        srt[i] = gRandom->Integer(964);

    TMath::Sort(964, srt, idx);
    */

    while (1)
    {
        TString line;
        line.ReadLine(fin);
        if (!fin)
            break;

        if (line.BeginsWith("#"))
        {
            //cout << line << endl;
            continue;
        }

        line=line.Strip(TString::kBoth);

        if (line.IsNull())
            continue;

        TObjArray *arr = line.Tokenize(' ');

        if (arr->GetSize()<8)
        {
            cout << "Skip3: " <<line << endl;
            delete arr;
            continue;
        }

        const TVector3 pos(atof((*arr)[0]->GetName()),
                           atof((*arr)[1]->GetName()),
                           atof((*arr)[2]->GetName()));

        const TVector3 norm(atof((*arr)[3]->GetName()),
                            atof((*arr)[4]->GetName()),
                            atof((*arr)[5]->GetName()));

        const Double_t F = atof((*arr)[6]->GetName());

        TString type = (*arr)[7]->GetName();
        type.Prepend("MMirror");

        TString msg;
        TClass *cls = MParContainer::GetClass(type);
        if (!cls)
        {
            *fLog << err << dbginf << "ERROR - Class " << type << " not in dictionary." << endl;
            return kFALSE;
        }

        if (!cls->InheritsFrom(MMirror::Class()))
        {
            *fLog << err << dbginf << "Cannot create new instance of class " << type << ": " << endl;
            *fLog << "Class doesn't inherit from MMirror." << endl;
            return kFALSE;
        }

        MMirror *m = (MMirror*)cls->New();
        if (!m)
        {
            *fLog << err << dbginf << "Cannot create new instance of class " << type << ": " << endl;
            *fLog << " - Class has no default constructor." << endl;
            *fLog << " - An abstract member functions of a base class is not overwritten." << endl;
            return kFALSE;
        }

        m->SetFocalLength(F);
        m->SetPosition(pos);
        m->SetNorm(norm);

        Int_t n = m->ReadM(*arr);
        if (n<=0)
        {
            *fLog << err << dbginf << "ERROR - ReadM failed." << endl;
            return kFALSE;
        }

        fMirrors.Add(m);

        //maxr = TMath::Max(maxr, TMath::Hypot(pos[i].X()+24.75, pos[i].Y()+24.75));
        //maxr = TMath::Max(maxr, TMath::Hypot(pos.X()+24.75, pos.Y()+24.75));

        delete arr;
    }

    InitMaxR();

    return kTRUE;

//    fMirrors.Sort();
/*
    for (int i=0; i<964; i++)
    {
        MMirror &ref = (MMirror&)*fMirrors[i];

        TArrayD dist(964);
        for (int j=0; j<964; j++)
        {
            const MMirror &mir = (MMirror&)*fMirrors[j];
            dist[j] = (ref-mir).Mod();
        }

        TArrayI idx(964);
        TMath::Sort(964, dist.GetArray(), idx.GetArray(), kFALSE);

        for (int j=0; j<964; j++)
        {
            ref.fNeighbors.Add(fMirrors[idx[j]]);
        }
    }*/
}

// --------------------------------------------------------------------------
//
// Paint the collection of mirrors
//
void MReflector::Paint(Option_t *o)
{
    fMirrors.Paint(o);
    /*
    TIter Next(&fMirrors);
    MMirror *m = 0;

    while ((m=static_cast<MMirror*>(Next())))
        m->Paint(o);*/
}

// --------------------------------------------------------------------------
//
// FileName: reflector.txt
//
Int_t MReflector::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{
    Bool_t rc = kFALSE;
    if (IsEnvDefined(env, prefix, "FileName", print))
    {
        rc = kTRUE;
        if (!ReadFile(GetEnvValue(env, prefix, "FileName", "")))
            return kERROR;
    }

    return rc;
}