Index: /trunk/Mars/msimreflector/MLens.cc
===================================================================
--- /trunk/Mars/msimreflector/MLens.cc	(revision 19539)
+++ /trunk/Mars/msimreflector/MLens.cc	(revision 19539)
@@ -0,0 +1,531 @@
+/* ======================================================================== *\
+!
+! *
+! * 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,  6/2019 <mailto:tbretz@physik.rwth-aachen.de>
+!
+!   Copyright: CheObs Software Development, 2000-2009
+!
+!
+\* ======================================================================== */
+
+//////////////////////////////////////////////////////////////////////////////
+//
+//  MLens
+//
+//////////////////////////////////////////////////////////////////////////////
+#include "MLens.h"
+
+#include <fstream>
+#include <errno.h>
+
+#include <stdlib.h> // atof (Ubuntu 8.10)
+
+#include <TClass.h>
+#include <TSystem.h>
+#include <TEllipse.h>
+
+#include "MQuaternion.h"
+#include "MMirror.h"
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+#include "MH.h"
+
+ClassImp(MLens);
+
+using namespace std;
+
+// --------------------------------------------------------------------------
+//
+// Default constructor
+//
+MLens::MLens(const char *name, const char *title)
+{
+    fName  = name  ? name  : "MLens";
+    fTitle = title ? title : "Parameter container storing a collection of several mirrors (reflector)";
+
+    fMirrors.SetOwner();
+}
+
+// --------------------------------------------------------------------------
+//
+// Set the SigmaPSF of all mirrors currently stored.
+//
+void MLens::SetSigmaPSF(Double_t psf)
+{
+    fMirrors.R__FOR_EACH(MMirror, SetSigmaPSF)(psf);
+}
+
+// --------------------------------------------------------------------------
+//
+// 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 MLens::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;
+    }
+}
+
+// --------------------------------------------------------------------------
+//
+// Return the total Area of all mirrors. Note, that it is recalculated
+// with any call.
+//
+Double_t MLens::GetA() const
+{
+    Double_t A = 0;
+
+    TIter Next(&fMirrors);
+    MMirror *m = 0;
+    while ((m=static_cast<MMirror*>(Next())))
+        A += m->GetA();
+
+    return A;
+}
+
+// --------------------------------------------------------------------------
+//
+// 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 **MLens::GetFirstPtr() const
+{
+    return (const MMirror**)fMirrors.GetObjectRef(0);
+}
+
+// --------------------------------------------------------------------------
+//
+// Get number of mirrors. There should be no holes in the array!
+//
+const UInt_t MLens::GetNumMirrors() const
+{
+    return fMirrors.GetEntriesFast();
+}
+
+// --------------------------------------------------------------------------
+//
+// Check with a rough estimate whether a photon can hit the reflector.
+//
+Bool_t MLens::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;
+}
+
+// Jeder Spiegel sollte eine Liste aller andern Spiegel in der
+// reihenfolge Ihrer Entfernung enthalten. Wir starten mit der Suche
+// immer beim zuletzt getroffenen Spiegel!
+//
+// --------------------------------------------------------------------------
+//
+// Loops over all mirrors of the reflector. After doing a rough check
+// whether the mirror can be hit at all the reflection is executed
+// calling the ExecuteMirror function of the mirrors.
+//
+// If a mirror was hit its index is retuened, -1 otherwise.
+//
+// FIXME: Do to lopping over all mirrors for all photons this is the
+// most time consuming function in teh reflector simulation. By a more
+// intelligent way of finding the right mirror then just testing all
+// this could be accelerated a lot.
+//
+Int_t MLens::ExecuteOptics(MQuaternion &p, MQuaternion &u) const
+{
+    //static const TObjArray *arr = &((MMirror*)fMirrors[0])->fNeighbors;
+
+    // This way of access is somuch faster than the program is
+    // a few percent slower if accessed by UncheckedAt
+    const MMirror **s = GetFirstPtr();
+    const MMirror **e = s+GetNumMirrors();
+    //const MMirror **s = (const MMirror**)fMirrors.GetObjectRef(0);
+    //const MMirror **e = s+fMirrors.GetEntriesFast();
+    //const MMirror **s = (const MMirror**)arr->GetObjectRef(0);
+    //const MMirror **e = s+arr->GetEntriesFast();
+
+    // Loop over all mirrors
+    for (const MMirror **m=s; m<e; m++)
+    {
+        const MMirror &mirror = **m;
+
+        // FIXME: Can we speed up using lookup tables or
+        //        indexed tables?
+
+        // MirrorShape: Check if this mirror can be hit at all
+        // This is to avoid time consuming calculation if there is no
+        // chance of a coincidence.
+        // FIXME: Inmprove search algorithm (2D Binary search?)
+        if (!mirror.CanHit(p))
+            continue;
+
+        // Make a local copy of position and direction which can be
+        // changed by ExecuteMirror.
+        MQuaternion q(p);
+        MQuaternion v(u);
+
+        // Check if this mirror is hit, and if it is hit return
+        // the reflected position and direction vector.
+        // If the mirror is missed we go on with the next mirror.
+        if (!mirror.ExecuteMirror(q, v))
+            continue;
+
+        // We hit a mirror. Restore the local copy of position and
+        // direction back into p und u.
+        p = q;
+        u = v;
+
+        //arr = &mirror->fNeighbors;
+
+        return m-s;
+    }
+
+    return -1;
+}
+
+// --------------------------------------------------------------------------
+//
+// I/O helper for ReadFile to avoid calling "delete arr" before every return
+//
+MMirror *MLens::EvalTokens(TObjArray &arr, Double_t defpsf) const
+{
+    if (arr.GetEntries()<9)
+    {
+        *fLog << err << "ERROR - Not enough arguments..." << endl;
+        return 0;
+    }
+
+    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()));
+
+    UInt_t n = 6;
+
+    TString six = arr[n]->GetName();
+
+    Char_t shape = 'S';
+    if (!six.IsFloat())
+    {
+        shape = six[0];
+        n++;
+    }
+
+    const Double_t F = atof(arr[n++]->GetName());
+
+    const TString val = arr[n++]->GetName();
+
+    const Double_t psf = val.IsFloat() ? val.Atof() : -1;
+
+    n += val.IsFloat() ? 1 : 0;
+
+    TString type = arr[n-1]->GetName();
+    type.Prepend("MMirror");
+
+    for (UInt_t i=0; i<n; i++)
+        delete arr.RemoveAt(i);
+    arr.Compress();
+
+    TString msg;
+    TClass *cls = MParContainer::GetClass(type);
+    if (!cls)
+    {
+        *fLog << err << "ERROR - Class " << type << " not in dictionary." << endl;
+        return 0;
+    }
+
+    if (!cls->InheritsFrom(MMirror::Class()))
+    {
+        *fLog << err << "ERROR - Cannot create new instance of class " << type << ": " << endl;
+        *fLog << "Class doesn't inherit from MMirror." << endl;
+        return 0;
+    }
+
+    MMirror *m = static_cast<MMirror*>(cls->New());
+    if (!m)
+    {
+        *fLog << err << "ERROR - 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 0;
+    }
+
+    m->SetPosition(pos);
+    m->SetNorm(norm);
+    m->SetShape(shape);
+    m->SetFocalLength(F);
+    m->SetSigmaPSF(psf>=0 ? psf : defpsf);
+
+    if (m->ReadM(arr)>=0)
+        return m;
+
+    *fLog << err << "ERROR - " << type << "::ReadM failed." << endl;
+
+    delete m;
+    return 0;
+}
+
+// --------------------------------------------------------------------------
+//
+// Read a reflector setup from a file. This needs improvemtn.
+//
+// The file structur is like:
+//
+//     x y z nx ny nz F [psf] Type ...
+//
+//  x:         x-coordinate of the mirror's center
+//  y:         y-coordinate of the mirror's center
+//  z:         z-coordinate of the mirror's center
+//  nx:        x-component of the normal vecor in the mirror center
+//  ny:        y-component of the normal vecor in the mirror center
+//  nz:        z-component of the normal vecor in the mirror center
+//  [shape]:   S for spherical <default>, P for parabolic
+//  F:         Focal distance of a spherical mirror
+//  [psf]:     This number is the psf given in the units of x,y,z and
+//             defined at the focal distance F. It can be used to overwrite
+//             the second argument given in ReadFile for individual mirrors.
+//  Type:      A instance of a mirrot of the class Type MMirrorType is created
+//             (Type can be, for example, Hex for for MMirrorHex).
+//  ...:       Additional arguments as defined in MMirrorType::ReadM
+//
+//
+// Coordinate System:
+//  The coordinate system is local in the reflectors frame.
+//  It is defined viewing from the back side of the reflector
+//  towards the camera. (x "right", y "up", z from reflector to camera)
+//  Note, that it is left-handed!
+//
+Bool_t MLens::ReadFile(TString fname, Double_t defpsf)
+{
+    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);
+    */
+
+    Int_t cnt = 0;
+
+    while (1)
+    {
+        TString line;
+        line.ReadLine(fin);
+        if (!fin)
+            break;
+
+        // Count lines
+        cnt++;
+
+        // Skip comments
+        if (line.BeginsWith("#"))
+            continue;
+
+        // Remove leading and trailing whitespaces
+        line=line.Strip(TString::kBoth);
+
+        // Skip empty lines
+        if (line.IsNull())
+            continue;
+
+        // Tokenize line
+        TObjArray *arr = line.Tokenize(' ');
+
+        // Evaluate tokens
+        MMirror *m = EvalTokens(*arr, defpsf);
+
+        // Delete now obsolete array
+        delete arr;
+
+        // Check if a new mirror could be created successfully
+        if (!m)
+        {
+            *fLog << err << "Error in line " << cnt << ": " << line << endl;
+            return kFALSE;
+        }
+
+        // Add new mirror to array
+        fMirrors.Add(m);
+    }
+
+    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]]);
+        }
+    }*/
+}
+
+// ------------------------------------------------------------------------
+//
+Bool_t MLens::WriteFile(TString fname) const
+{
+    gSystem->ExpandPathName(fname);
+
+    ofstream fout(fname);
+    if (!fout)
+    {
+        *fLog << err << "Cannot open file " << fname << ": ";
+        *fLog << (errno!=0?strerror(errno):"Insufficient memory for decompression") << endl;
+        return kFALSE;
+    }
+
+    TIter Next(&fMirrors);
+    MMirror *m = 0;
+    while ((m=static_cast<MMirror*>(Next())))
+    {
+        //  x:         x-coordinate of the mirror's center
+        //  y:         y-coordinate of the mirror's center
+        //  z:         z-coordinate of the mirror's center
+        //  nx:        x-component of the normal vecor in the mirror center
+        //  ny:        y-component of the normal vecor in the mirror center
+        //  nz:        z-component of the normal vecor in the mirror center
+        //  [shape]:   S for spherical <default>, P for parabolic
+        //  F:         Focal distance of a spherical mirror
+        //  [psf]:     This number is the psf given in the units of x,y,z and
+        //             defined at the focal distance F. It can be used to overwrite
+        //             the second argument given in ReadFile for individual mirrors.
+        //  Type:      A instance of a mirrot of the class Type MMirrorType is created
+        //             (Type can be, for example, Hex for for MMirrorHex).
+        //  ...:       Additional arguments as defined in MMirrorType::ReadM
+
+        fout << setw(12) << m->X() << " ";
+        fout << setw(12) << m->Y() << " ";
+        fout << setw(12) << m->Z() << "  ";
+        fout << setw(12) << m->Nx() << " ";
+        fout << setw(12) << m->Ny() << " ";
+        fout << setw(12) << m->Nz() << " ";
+        if (m->GetShape()==1)
+            fout << "P ";
+        fout << m->GetFocalLength() << " ";
+        // cout << m->GetSigmaPSF() << " ";
+        fout << m->ClassName()+7 << " ";
+        m->WriteM(fout);
+        fout << endl;
+    }
+    return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// Print the collection of mirrors
+//
+void MLens::Print(Option_t *o) const
+{
+    *fLog << all << GetDescriptor() << " (" << GetNumMirrors() << "): " << endl;
+
+    fMirrors.Print(o);
+}
+
+// --------------------------------------------------------------------------
+//
+// Paint the collection of mirrors
+//
+void MLens::Paint(Option_t *o)
+{
+    if (!TString(o).Contains("same", TString::kIgnoreCase))
+        MH::SetPadRange(-fMaxR*1.01, -fMaxR*1.01, fMaxR*1.01, fMaxR*1.01);
+
+    fMirrors.Paint(o);
+
+    TEllipse e;
+    e.SetFillStyle(0);
+    e.SetLineColor(17);
+    e.PaintEllipse(0, 0, fMaxR, fMaxR, 0, 360, 0);
+}
+
+// --------------------------------------------------------------------------
+//
+// SigmaPSF: -1
+// FileName: reflector.txt
+//
+// SigmaPSF can be used to set a default for the psf of the mirrors
+// read from the file. Note, that this can be overwritten for individual
+// mirrors in the file. The SigmaPSF is the sigma of a 1D-Gauss fitted
+// to the radial profile of the light distribution.
+//
+// For details on the file structure see MLens::ReadFile
+//
+Int_t MLens::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
+{
+    Bool_t rc = kFALSE;
+
+    Double_t psf = -1;
+    if (IsEnvDefined(env, prefix, "SetSigmaPSF", print))
+    {
+        rc = kTRUE;
+        psf = GetEnvValue(env, prefix, "SetSigmaPSF", -1);
+    }
+
+    if (IsEnvDefined(env, prefix, "FileName", print))
+    {
+        rc = kTRUE;
+        if (!ReadFile(GetEnvValue(env, prefix, "FileName", ""), psf))
+            return kERROR;
+    }
+
+    return rc;
+}
Index: /trunk/Mars/msimreflector/MLens.h
===================================================================
--- /trunk/Mars/msimreflector/MLens.h	(revision 19539)
+++ /trunk/Mars/msimreflector/MLens.h	(revision 19539)
@@ -0,0 +1,66 @@
+#ifndef MARS_MLens
+#define MARS_MLens
+
+#ifndef MARS_MOptics
+#include "MOptics.h"
+#endif
+
+#ifndef ROOT_TObjArray
+#include <TObjArray.h>
+#endif
+
+class MQuaternion;
+class MMirror;
+
+class MLens : public MOptics
+{
+private:
+    //Simple Array               // 5.1s
+    TObjArray fMirrors;          // 6.1s   (Pointer)
+    //TObjArray fMirrors;        // 6.1s   (GetObjectRef)
+    //TObjArray fMirrors;        // 8.3s   (Next)
+    //TObjArray fMirrors;        // 10.1s  (UncheckedAt)
+    //TList fMirrors;            // 10.7s
+    //TOrdCollection fMirrors;   // 23.4s
+
+    Double_t fMaxR;
+
+    void InitMaxR();
+
+    // Helper for I/O
+    MMirror *EvalTokens(TObjArray &arr, Double_t defpsf) const;
+
+    // MParContainer
+    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
+
+public:
+    MLens(const char *name=NULL, const char *title=NULL);
+
+    const MMirror **GetFirstPtr() const;
+    const UInt_t GetNumMirrors() const;
+
+    const MMirror *GetMirror(UInt_t idx) const { return idx>=GetNumMirrors()?0:*(GetFirstPtr()+idx); }
+
+    Bool_t ReadFile(TString fname, Double_t defpsf=-1);
+    Bool_t WriteFile(TString fname) const;
+
+    Double_t GetMaxR() const { return fMaxR; }
+    Double_t GetA() const;
+
+    virtual Bool_t CanHit(const MQuaternion &p) const;
+
+    Int_t ExecuteOptics(MQuaternion &p, MQuaternion &u) const;
+
+    void SetSigmaPSF(Double_t psf);
+
+    virtual Bool_t IsValid() const { return fMirrors.GetEntries(); }
+
+
+    // TObject
+    void Paint(Option_t *o);
+    void Print(Option_t *o) const;
+
+    ClassDef(MLens, 1) // Parameter container storing the description of a lens
+};
+    
+#endif