Index: /trunk/MagicSoft/Mars/mbase/BaseLinkDef.h
===================================================================
--- /trunk/MagicSoft/Mars/mbase/BaseLinkDef.h	(revision 9233)
+++ /trunk/MagicSoft/Mars/mbase/BaseLinkDef.h	(revision 9234)
@@ -29,4 +29,5 @@
 
 #pragma link C++ class MSpline3+;
+#pragma link C++ class MQuaternion+;
 
 #pragma link C++ class MString+;
Index: /trunk/MagicSoft/Mars/mbase/MArrayI.cc
===================================================================
--- /trunk/MagicSoft/Mars/mbase/MArrayI.cc	(revision 9233)
+++ /trunk/MagicSoft/Mars/mbase/MArrayI.cc	(revision 9234)
@@ -36,4 +36,10 @@
 #include "MArrayI.h"
 
+#include "MMath.h"
+
 ClassImp(MArrayI);
 
+void MArrayI::ReSort(Bool_t down)
+{
+    MMath::ReSort(fN, fArray, down);
+}
Index: /trunk/MagicSoft/Mars/mbase/MArrayI.h
===================================================================
--- /trunk/MagicSoft/Mars/mbase/MArrayI.h	(revision 9233)
+++ /trunk/MagicSoft/Mars/mbase/MArrayI.h	(revision 9234)
@@ -73,4 +73,25 @@
         // Add char c at position i. Check for out of bounds.
         fArray[i] = c;
+    }
+
+    void Add(Int_t c)
+    {
+        Set(fN+1);
+        fArray[fN-1] = c;
+    }
+
+    Int_t Find(Int_t c) const
+    {
+        for (Int_t *ptr=fArray; ptr<fArray+fN; ptr++)
+            if (*ptr==c)
+                return ptr-fArray;
+
+        return -1;
+    }
+
+    void AddUniq(Int_t c)
+    {
+        if (Find(c)<0)
+            Add(c);
     }
 
@@ -157,4 +178,6 @@
     }
 
+    void ReSort(Bool_t down=kFALSE);
+
     ClassDef(MArrayI, 1)  //Array of Int_t
 };
Index: /trunk/MagicSoft/Mars/mbase/MMath.cc
===================================================================
--- /trunk/MagicSoft/Mars/mbase/MMath.cc	(revision 9233)
+++ /trunk/MagicSoft/Mars/mbase/MMath.cc	(revision 9234)
@@ -1,4 +1,4 @@
 /* ======================================================================== *\
-! $Name: not supported by cvs2svn $:$Id: MMath.cc,v 1.44 2009-01-17 14:52:37 tbretz Exp $
+! $Name: not supported by cvs2svn $:$Id: MMath.cc,v 1.45 2009-01-21 14:22:39 tbretz Exp $
 ! --------------------------------------------------------------------------
 !
@@ -347,6 +347,6 @@
         *ptr = cpy[*idx++];
 
-    delete cpy;
-    delete pos;
+    delete [] cpy;
+    delete [] pos;
 }
 
Index: /trunk/MagicSoft/Mars/mbase/MQuaternion.cc
===================================================================
--- /trunk/MagicSoft/Mars/mbase/MQuaternion.cc	(revision 9234)
+++ /trunk/MagicSoft/Mars/mbase/MQuaternion.cc	(revision 9234)
@@ -0,0 +1,74 @@
+/* ======================================================================== *\
+!
+! *
+! * 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  11/2008 <mailto:tbretz@astro.uni-wuerzburg.de>
+!
+!   Copyright: MAGIC Software Development, 2000-2009
+!
+!
+\* ======================================================================== */
+
+//////////////////////////////////////////////////////////////////////////////
+//
+//  MQuaternion
+//
+// The MQuaternion is derived from TQuaternion. A quaternion is a four vector
+// which can store space and time (like a lorentz vector).
+//
+// There are a few advantages of the TQuaternion class over the
+// TLorentzVector, namely the implementation of a direct algebra with
+// just the space part of the vector keeping the time as it is.
+// (This is useful, e.g, for rotations and shift just in space).
+//
+//  - You can construct the MQuaternion from a TQuaternion or a TVector3 and
+//    time.
+//  - Multiplying the MQuaternion with a TRotation with rotate just the
+//    space-part.
+//  - You can access the data members with X(), Y(), Z() and T()
+//  - To get the length or squared-length of the space-vector use R() and R2()
+//  - Access the 2D vector (x/y) with XYvector()
+//  - Thera are a few new function to propagate a MQuaternion along a trajectory
+//    in space and time (also expressed as an MQuaternion with a direction
+//    vector and a speed) Here we assume v>0.
+//
+//    + PropagateDz(MQuaternion &w, Double_t dz)
+//
+//      If dz is positive the position is propagated along the given trajectory
+//      in space (such that the z-component will increase by dz) and
+//      forward in time. If dz<0 the result is vice versa.
+//
+//    + PropagateZ0(MQuaternion &w)
+//
+//      This is an abbreviation for Propagate(w, -Z()). It propagates the
+//      position such that its z-component will vanish. If this is along
+//      the given trajectory time will increase if it is backward time
+//      will decrease.
+//
+//    + PropagateZ(MQuaternion &w, Double_t z)
+//
+//      This is an abbreviation for Propagate(w, z-Z()). It propagates the
+//      position such that its z-component will become z. If this is along
+//      the given trajectory time will increase if it is backward time
+//      will decrease.
+//
+//////////////////////////////////////////////////////////////////////////////
+#include "MQuaternion.h"
+
+ClassImp(MQuaternion);
+
+using namespace std;
+
Index: /trunk/MagicSoft/Mars/mbase/MQuaternion.h
===================================================================
--- /trunk/MagicSoft/Mars/mbase/MQuaternion.h	(revision 9234)
+++ /trunk/MagicSoft/Mars/mbase/MQuaternion.h	(revision 9234)
@@ -0,0 +1,133 @@
+#ifndef MARS_MQuaternion
+#define MARS_MQuaternion
+
+#if 1
+
+// We prefer to derive from TQuaternion instead of TLorantzVector
+// because TQuaternion implements vector algebra with just the 3D vector
+
+#ifndef ROOT_TQuaternion
+#include <TQuaternion.h>
+#endif
+
+class MQuaternion : public TQuaternion
+{
+public:
+    MQuaternion(const TQuaternion &q) : TQuaternion(q) { }
+    MQuaternion(const TVector3 &v, Double_t t=0) : TQuaternion(v, t) { }
+    void operator*=(const TRotation &r)
+    {
+        fVectorPart *= r;
+    }
+    Double_t X() const { return fVectorPart.X(); }
+    Double_t Y() const { return fVectorPart.Y(); }
+    Double_t Z() const { return fVectorPart.Z(); }
+    Double_t T() const { return fRealPart; }
+
+    // It seems to be a little bit faster than X*X+Y*Y
+    Double_t R2() const { return XYvector().Mod2(); }
+    Double_t R() const { return XYvector().Mod(); }
+
+    void PropagateDz(const MQuaternion &w, const Double_t dz)
+    {
+        *this += dz/w.Z()*w;
+    }
+
+    // Propagates the particle by a distance f in z along
+    // its trajectory w, if f is positive, in the opposite
+    // direction otherwise.
+    void PropagateZ(const MQuaternion &w, const Double_t z)
+    {
+        PropagateDz(w, z-Z());
+
+        // z=3400, Z= 1700, t=0, c=1    -=  3400/-5*-5     -= 3400       Z=0,    c>0
+        //                              +=  1700/-5*-5     += 1700       Z=1700, c>0
+        // z=3400, Z=-1700, t=0, c=1    -= -3400/-5*-5     -= -1700      Z=0, c<0
+
+        // z=3400, Z= 1700, t=0, c=1    -=  (3400-1700)/-5*-5     -= 3400       Z=0,    c>0
+    }
+
+    // Move the photon along its trajectory to the x/y plane
+    // so that z=0. Therefor stretch the vector until
+    // its z-component vanishes.
+    //p -= p.Z()/u.Z()*u;
+    void PropagateZ0(const MQuaternion &w)
+    {
+        // If z>0 we still have to move by a distance of z.
+        // If z<0 we have to move in the opposite direction.
+        //  --> z has the right sign for PropagateZ
+        PropagateDz(w, -Z());
+
+        // Z= 1700, t=0, c=1    -=  1700/-5*-5     -=  1700 +c     Z=0, c>0
+        // Z=-1700, t=0, c=1    -= -1700/-5*-5     -= -1700 -c     Z=0, c<0
+
+
+        // Z= 1700, t=0, c=1    -=  1700/ 5* 5     -=  1700 -c     Z=0, c<0
+        // Z=-1700, t=0, c=1    -= -1700/ 5* 5     -= -1700 +c     Z=0, c>0
+
+        //PropagateZ(w, Z());
+    }
+
+    TVector2 XYvector() const { return fVectorPart.XYvector(); }
+
+    //void Normalize() { fVectorPart *= TMath::Sqrt(1 - R2())/Z(); }
+
+    ClassDef(MQuaternion, 1)
+};
+
+#else
+
+#ifndef ROOT_TLorentzVector
+#include <TLorentzVector.h>
+#endif
+
+class MQuaternion : public TLorentzVector
+{
+public:
+    //MQuaternion(const TLorentzVector &q) : TLorentzVector(q) { }
+    MQuaternion(const TVector3 &v, Double_t t=0) : TLorentzVector(v, t) { }
+    /*
+    void operator*=(const TRotation &r)
+    {
+        fVectorPart *= r;
+    }
+    Double_t X() const { return fVectorPart.X(); }
+    Double_t Y() const { return fVectorPart.Y(); }
+    Double_t Z() const { return fVectorPart.Z(); }
+    Double_t T() const { return fRealPart; }
+    */
+
+    // It seems to be a little bit faster than X*X+Y*Y
+    Double_t R2() const { return Perp2(); }
+    Double_t R() const { return Perp(); }
+
+    // Propagates the particle by a distance f in z along
+    // its trajectory w, if f is positive, in the opposite
+    // direction otherwise.
+    void PropagateZ(const MQuaternion &w, const Double_t f)
+    {
+        *this += f/TMath::Abs(w.Z())*w;
+    }
+
+    // Move the photon along its trajectory to the x/y plane
+    // so that z=0. Therefor stretch the vector until
+    // its z-component vanishes.
+    //p -= p.Z()/u.Z()*u;
+    void PropagateZ0(const MQuaternion &w)
+    {
+        // If z>0 we still have to move by a distance of z.
+        // If z<0 we have to move in th eopposite direction.
+        //  --> z has the right sign for PropagateZ
+        PropagateZ(w, Z());
+    }
+
+    TVector2 XYvector() const { return Vect().XYvector(); }
+
+    //void Normalize() { fVectorPart *= TMath::Sqrt(1 - R2())/Z(); }
+
+    ClassDef(MQuaternion, 0)
+};
+
+#endif
+
+#endif
Index: /trunk/MagicSoft/Mars/mbase/Makefile
===================================================================
--- /trunk/MagicSoft/Mars/mbase/Makefile	(revision 9233)
+++ /trunk/MagicSoft/Mars/mbase/Makefile	(revision 9234)
@@ -24,4 +24,5 @@
            MMath.cc \
            MSpline3.cc \
+	   MQuaternion.cc \
            MEnv.cc \
 	   MLog.cc \
Index: /trunk/MagicSoft/Mars/mgeom/MGeomCam.h
===================================================================
--- /trunk/MagicSoft/Mars/mgeom/MGeomCam.h	(revision 9233)
+++ /trunk/MagicSoft/Mars/mgeom/MGeomCam.h	(revision 9234)
@@ -13,4 +13,7 @@
 #ifndef ROOT_TArrayS
 #include <TArrayS.h>
+#endif
+#ifndef MARS_MQuaternion
+#include "MQuaternion.h"
 #endif
 
@@ -109,4 +112,7 @@
     Int_t GetNeighbor(UInt_t idx, Int_t dir) const;
 
+    virtual Bool_t HitFrame(MQuaternion p, const MQuaternion &u) const { return kFALSE; }
+    virtual Bool_t HitDetector(const MQuaternion &p) const { return kFALSE; }
+
     void Print(Option_t *opt=NULL) const;
 
Index: /trunk/MagicSoft/Mars/mgeom/MGeomCamDwarf.cc
===================================================================
--- /trunk/MagicSoft/Mars/mgeom/MGeomCamDwarf.cc	(revision 9233)
+++ /trunk/MagicSoft/Mars/mgeom/MGeomCamDwarf.cc	(revision 9234)
@@ -112,4 +112,41 @@
 // --------------------------------------------------------------------------
 //
+// Check if the photon which is flying along the trajectory u has passed
+// (or will pass) the frame of the camera (and consequently get
+// absorbed). The position p and direction u must be in the
+// telescope coordinate frame, which is z parallel to the focal plane,
+// x to the right and y upwards, looking from the mirror towards the camera.
+//
+// The units are cm.
+//
+Bool_t MGeomCamDwarf::HitFrame(MQuaternion p, const MQuaternion &u) const
+{
+    // z is defined from the mirror (0) to the camera (z>0).
+    // Thus we just propagate to the focal plane (z=fDist)
+    //p -= 1700./u.Z()*u;
+    p.PropagateZ(u, GetCameraDist()*100); // m->cm
+
+    // Add 10% to the max radius and convert from mm to cm
+    return p.R()<GetMaxRadius()*0.11;//TMath::Abs(p.X())<65 && TMath::Abs(p.Y())<65;
+}
+
+// --------------------------------------------------------------------------
+//
+// Check if the position given in the focal plane (so z can be ignored)
+// is a position which might hit the detector. It is meant to be a rough
+// and fast estimate not a precise calculation. All positions dicarded
+// must not hit the detector. All positions accepted might still miss
+// the detector.
+//
+// The units are cm.
+//
+Bool_t MGeomCamDwarf::HitDetector(const MQuaternion &v) const
+{
+    // Add 10% to the max radius and convert from mm to cm
+    return v.R()<GetMaxRadius()*0.11;//60.2;
+}
+
+// --------------------------------------------------------------------------
+//
 // Calculate in the direction 0-5 (kind of sector) in the ring-th ring
 // the x and y coordinate of the i-th pixel. The unitx are unity,
Index: /trunk/MagicSoft/Mars/mgeom/MGeomCamDwarf.h
===================================================================
--- /trunk/MagicSoft/Mars/mgeom/MGeomCamDwarf.h	(revision 9233)
+++ /trunk/MagicSoft/Mars/mgeom/MGeomCamDwarf.h	(revision 9234)
@@ -28,7 +28,9 @@
     TObject *Clone(const char *newname) const;
 
-    ClassDef(MGeomCamDwarf, 1)		// Geometry class for the Dwarf camera
+    Bool_t HitFrame(MQuaternion p, const MQuaternion &u) const;
+    Bool_t HitDetector(const MQuaternion &p) const;
+
+    ClassDef(MGeomCamDwarf, 1) // Geometry class for the Dwarf camera
 };
 
 #endif
-
Index: /trunk/MagicSoft/Mars/mgeom/MGeomCamMagic.cc
===================================================================
--- /trunk/MagicSoft/Mars/mgeom/MGeomCamMagic.cc	(revision 9233)
+++ /trunk/MagicSoft/Mars/mgeom/MGeomCamMagic.cc	(revision 9234)
@@ -49,4 +49,35 @@
     CreateNN();
     InitGeometry();
+}
+
+// --------------------------------------------------------------------------
+//
+// Check if the photon which is flying along the trajectory u has passed
+// (or will pass) the frame of the camera (and consequently get
+// absorbed). The position p and direction u must be in the
+// telescope coordinate frame, which is z parallel to the focal plane,
+// x to the right and y upwards, looking from the mirror towards the camera.
+//
+Bool_t MGeomCamMagic::HitFrame(MQuaternion p, const MQuaternion &u) const
+{
+    // z is defined from the mirror (0) to the camera (z>0).
+    // Thus we just propagate to the focal plane (z=fDist)
+    //p -= 1700./u.Z()*u;
+    p.PropagateZ(u, GetCameraDist()*100);
+
+    return TMath::Abs(p.X())<65 && TMath::Abs(p.Y())<65;
+}
+
+// --------------------------------------------------------------------------
+//
+// Check if the position given in the focal plane (so z can be ignored)
+// is a position which might hit the detector. It is meant to be a rough
+// and fast estimate not a precise calculation. All positions dicarded
+// must not hit the detector. All positions accepted might still miss
+// the detector.
+//
+Bool_t MGeomCamMagic::HitDetector(const MQuaternion &v) const
+{
+    return v.R()<60.2;
 }
 
Index: /trunk/MagicSoft/Mars/mgeom/MGeomCamMagic.h
===================================================================
--- /trunk/MagicSoft/Mars/mgeom/MGeomCamMagic.h	(revision 9233)
+++ /trunk/MagicSoft/Mars/mgeom/MGeomCamMagic.h	(revision 9234)
@@ -15,4 +15,7 @@
     MGeomCamMagic(const char *name=NULL);
 
+    Bool_t HitFrame(MQuaternion p, const MQuaternion &u) const;
+    Bool_t HitDetector(const MQuaternion &p) const;
+
     ClassDef(MGeomCamMagic, 1)		// Geometry class for the Magic camera
 };
Index: /trunk/MagicSoft/Mars/msim/MHPhotonEvent.cc
===================================================================
--- /trunk/MagicSoft/Mars/msim/MHPhotonEvent.cc	(revision 9233)
+++ /trunk/MagicSoft/Mars/msim/MHPhotonEvent.cc	(revision 9234)
@@ -44,17 +44,11 @@
 //
 // Type 3:
-//   The maximum in x and y is determined from MCamera->GetMaxR();
+//   The maximum in x and y is determined from MGeomCam->GetMaxR();
 //
 // Type 4:
 //   As type 3 but divided by 10.
 //
-// Type 5:
-//   The maximum in x and y is determined from MGeomCam->GetMaxR();
-//
-// Type 6:
-//   As type 5 but divided by 10.
-//
 // The binning is optimized using MH::FindGoodLimits. The number of bins
-// in 100 in the default case and 50 for type 3-6.
+// in 100 in the default case and 50 for type 3-4.
 //
 // Fill expects a MPhotonEvent (the second argumnet in MFillH).
@@ -78,5 +72,5 @@
 #include "MCorsikaRunHeader.h"
 
-#include "../msimreflector/MSimReflector.h" // MCamera, MReflector
+#include "MReflector.h"
 
 ClassImp(MHPhotonEvent);
@@ -198,18 +192,4 @@
     case 4:
         {
-            MCamera *c = (MCamera*)pList->FindObject("MCamera");
-            if (!c)
-            {
-                *fLog << err << "MCamera not found... aborting." << endl;
-                return kFALSE;
-            }
-            xmax = fType==3 ? c->GetMaxR() : c->GetMaxR()/10;
-            num  = 50;
-
-            break;
-        }
-    case 5:
-    case 6:
-        {
             MGeomCam *c = (MGeomCam*)pList->FindObject("MGeomCam");
             if (!c)
@@ -218,5 +198,5 @@
                 return kFALSE;
             }
-            xmax = fType==5 ? c->GetMaxRadius() : c->GetMaxRadius()/10;
+            xmax = fType==3 ? c->GetMaxRadius()/10 : c->GetMaxRadius()/100;
             num  = 50;