Changeset 9234 for trunk/MagicSoft

Timestamp:

01/21/09 14:26:15 (16 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• mbase/BaseLinkDef.h (modified) (1 diff)
• mbase/MArrayI.cc (modified) (1 diff)
• mbase/MArrayI.h (modified) (2 diffs)
• mbase/MMath.cc (modified) (2 diffs)
• mbase/MQuaternion.cc (added)
• mbase/MQuaternion.h (added)
• mbase/Makefile (modified) (1 diff)
• mgeom/MGeomCam.h (modified) (2 diffs)
• mgeom/MGeomCamDwarf.cc (modified) (1 diff)
• mgeom/MGeomCamDwarf.h (modified) (1 diff)
• mgeom/MGeomCamMagic.cc (modified) (1 diff)
• mgeom/MGeomCamMagic.h (modified) (1 diff)
• msim/MHPhotonEvent.cc (modified) (4 diffs)

trunk/MagicSoft/Mars/mbase/BaseLinkDef.h

@@ -29 +29 @@
 
 #pragma link C++ class MSpline3+;
+#pragma link C++ class MQuaternion+;
 
 #pragma link C++ class MString+;

trunk/MagicSoft/Mars/mbase/MArrayI.cc

@@ -36 +36 @@
 #include "MArrayI.h"
 
+#include "MMath.h"
+
 ClassImp(MArrayI);
 
+void MArrayI::ReSort(Bool_t down)
+{
+    MMath::ReSort(fN, fArray, down);
+}

trunk/MagicSoft/Mars/mbase/MArrayI.h

@@ -73 +73 @@
         // 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 +178 @@
     }
 
+    void ReSort(Bool_t down=kFALSE);
+
     ClassDef(MArrayI, 1)  //Array of Int_t
 };

trunk/MagicSoft/Mars/mbase/MMath.cc

@@ -1 +1 @@
 /* ======================================================================== *\
-! $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 +347 @@
         *ptr = cpy[*idx++];
 
-    delete cpy;
-    delete pos;
+    delete [] cpy;
+    delete [] pos;
 }
 

trunk/MagicSoft/Mars/mbase/Makefile

@@ -24 +24 @@
            MMath.cc \
            MSpline3.cc \
+           MQuaternion.cc \
            MEnv.cc \
            MLog.cc \

trunk/MagicSoft/Mars/mgeom/MGeomCam.h

@@ -13 +13 @@
 #ifndef ROOT_TArrayS
 #include <TArrayS.h>
+#endif
+#ifndef MARS_MQuaternion
+#include "MQuaternion.h"
 #endif
 
@@ -109 +112 @@
     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;
 

trunk/MagicSoft/Mars/mgeom/MGeomCamDwarf.cc

@@ -112 +112 @@
 // --------------------------------------------------------------------------
 //
+// 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,

trunk/MagicSoft/Mars/mgeom/MGeomCamDwarf.h

@@ -28 +28 @@
     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
-

trunk/MagicSoft/Mars/mgeom/MGeomCamMagic.cc

@@ -49 +49 @@
     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;
 }
 

trunk/MagicSoft/Mars/mgeom/MGeomCamMagic.h

@@ -15 +15 @@
     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
 };

trunk/MagicSoft/Mars/msim/MHPhotonEvent.cc

@@ -44 +44 @@
 //
 // 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 +72 @@
 #include "MCorsikaRunHeader.h"
 
-#include "../msimreflector/MSimReflector.h" // MCamera, MReflector
+#include "MReflector.h"
 
 ClassImp(MHPhotonEvent);
@@ -198 +192 @@
     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 +198 @@
                 return kFALSE;
             }
-            xmax = fType==5 ? c->GetMaxRadius() : c->GetMaxRadius()/10;
+            xmax = fType==3 ? c->GetMaxRadius()/10 : c->GetMaxRadius()/100;
             num  = 50;