Changeset 1434

Timestamp:

07/25/02 09:14:36 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (2 diffs)
• NEWS (modified) (3 diffs)
• manalysis/MHillas.cc (modified) (14 diffs)
• manalysis/MHillas.h (modified) (2 diffs)
• manalysis/MHillasSrc.cc (modified) (6 diffs)
• manalysis/MHillasSrc.h (modified) (2 diffs)
• manalysis/MHillasSrcCalc.cc (modified) (1 diff)
• mgui/MineSweeper.h (modified) (1 diff)

trunk/MagicSoft/Mars/Changelog

@@ -1 +1 @@
                                                                   -*-*- END -*-*-
+
+ 2002/07/23: Thomas Bretz
+
+   * mgui/MineSweeper.cc, mgui/MCamDisplay.cc:
+     - fixed a bug which could caused the destructor of the created canvas
+       not to be called if the object was deleted.
+
+   * mgui/MCamDisplay.h:
+     - some reordering
+
+   * mgui/MHexagon.cc:
+     - return the correct distance to primitive
+     - do a sanity check in DistancetoPrimitive
+
+   * mhist/HistLinkDef.h:
+     - added missing classes
+
+
+
+ 2002/07/23: Wolfgang Wittek, Thomas Bretz
+
+   * manalysis/MHillas.cc, manalysis/MHillasSrc.cc:
+     - replaced algorithms by the ones from the corresponding TDAS note
+     - made results of calculations (hopefully) more accurate
+     - added some output warnings and sanity checks
+     - changed the default parameters in the Reset function
+
+   * manalysis/MHillasSrc.[h,cc]:
+     - added fCosDeltaAlpha
+     - made Calc returning a Bool_t
+
+   * manalysis/MHillasSrcCalc.cc:
+     - make use of the return value of MHillasSrc::Calc 
+
+
 
  2002/07/22: Abelardo Moralejo
@@ -8 +43 @@
    * macros/MagicHillas.C:
      -added example on how to use MCerPhotCalc2
+
+
 
  2002/07/22: Thomas Bretz

trunk/MagicSoft/Mars/NEWS

@@ -2 +2 @@
 
  *** Version 0.8
+
+   - added all subdirectories as Include-Paths for ACLIC (s. rootlogon.C)
 
    - changed trigrate.C: now it calculates the pure NSB triggers from
@@ -21 +23 @@
      (MMultiDimDistCalc, multidimdist.C, MHMatrix)
 
-   - Changed the Alpha range according to TDAS 02-03 to [-90°,90°]
-
    - Added a camera display displaying the pixel numbers to the
      event display gui.
@@ -34 +34 @@
    - Added the number of used and core pixels after image cleaning
      to MHillas
+
+   - Changed the algorithm used in MHillas according to TDAS 02-03
+
+   - Changed the Alpha range according to TDAS 02-03 to [-90°,90°]
+
+   - Changed the algorithm used in MHillasSrc according to TDAS 02-03
+
+   - Added fCosDeltaAlpha to MHillasSrc
 
    - Fixed a bug in the table describing the Next Neighbours of the

trunk/MagicSoft/Mars/manalysis/MHillas.cc

@@ -19 +19 @@
 !   Author(s): Harald Kornmayer 1/2001
 !   Author(s): Rudolf Bock     10/2001 <mailto:Rudolf.Bock@cern.ch>
+!   Author(s): Wolfgang Wittek  6/2002 <mailto:wittek@mppmu.mpg.de>
 !
 !   Copyright: MAGIC Software Development, 2000-2002
@@ -35 +36 @@
 // Version 1:
 // ----------
-// fLength   major axis of ellipse
-// fWidth    minor axis
-// fDelta    angle of major axis wrt x-axis
-// fSize     total sum of pixels
-// fMeanX    x of center of ellipse
-// fMeanY    y of center of ellipse
+// fLength   [mm]       major axis of ellipse
+// fWidth    [mm]       minor axis
+// fDelta    [rad]      angle of major axis with x-axis
+//                      by definition the major axis is pointing into
+//                      the hemisphere x>0, thus -pi/2 < delta < pi/2
+// fSize     [#CerPhot] total sum of pixels
+// fMeanX    [mm]       x of center of ellipse
+// fMeanY    [mm]       y of center of ellipse
 //
 // Version 2:
@@ -91 +94 @@
 // --------------------------------------------------------------------------
 //
+// Initializes the values with defaults. For the default values see the
+// source code.
+//
 void MHillas::Reset()
 {
-    fLength = 0;
-    fWidth  = 0;
-    fDelta  = 0;
+    fLength = -1;
+    fWidth  = -1;
+    fDelta  =  0;
 
     fSize   = -1;
@@ -122 +128 @@
     *fLog << " - Length   [mm]  = " << fLength << endl;
     *fLog << " - Width    [mm]  = " << fWidth  << endl;
+    *fLog << " - Delta    [deg] = " << fDelta*kRad2Deg << endl;
+    *fLog << " - Size     [1]   = " << fSize   << " #CherPhot"   << endl;
     *fLog << " - Meanx    [mm]  = " << fMeanX  << endl;
     *fLog << " - Meany    [mm]  = " << fMeanY  << endl;
-    *fLog << " - Delta    [deg] = " << fDelta*kRad2Deg << endl;
     *fLog << " - atg(y/x) [deg] = " << atg     << endl;
-    *fLog << " - Size     [1]   = " << fSize   << " #CherPhot"   << endl;
 }
 
@@ -207 +213 @@
 // Calculate the image parameters from a Cherenkov photon event
 // assuming Cher.photons/pixel and pixel coordinates are given
+// In case you don't call Calc from within an eventloop make sure, that
+// you call the Reset member function before.
 //
 Bool_t MHillas::Calc(const MGeomCam &geom, const MCerPhotEvt &evt)
 {
-    // FIXME: MHillas::Calc is rather slow at the moment because it loops
-    //    unnecessarily over all pixels in all its loops (s.MImgCleanStd)
-    //    The speed could be improved very much by using Hash-Tables
-    //    (linked lists, see THashTable and THashList, too)
-    //
-
     const UInt_t npixevt = evt.GetNumPixels();
 
@@ -222 +224 @@
     //
     if (npixevt <= 2)
+    {
+        *fLog << warn << "MHillas::Calc: Event has less than two pixels... skipped." << endl;
         return kFALSE;
+    }
 
     //
     // calculate mean value of pixel coordinates and fSize
     // -----------------------------------------------------
+    //
+    // Because this are only simple sums of roughly 600 values
+    // with an accuracy less than three digits there should not
+    // be any need to use double precision (Double_t) for the
+    // calculation of fMeanX, fMeanY and fSize.
     //
     fMeanX = 0;
@@ -235 +245 @@
     fNumCorePixels = 0;
 
-    //
-    // FIXME! npix should be retrieved from MCerPhotEvt
-    //
     for (UInt_t i=0; i<npixevt; i++)
     {
@@ -247 +254 @@
         const MGeomPix &gpix = geom[pix.GetPixId()];
 
-        const float nphot = pix.GetNumPhotons();
+        const Float_t nphot = pix.GetNumPhotons();
+
+        if (nphot==0)
+            *fLog << warn << GetDescriptor() << ": Pixel #" << pix.GetPixId() << " has no photons." << endl;
 
         fSize  += nphot;                             // [counter]
@@ -260 +270 @@
 
     //
-    // sanity check
-    //
-    if (fSize==0 || fNumUsedPixels<=2)
+    // sanity checks
+    //
+    if (fSize==0)
+    {
+        *fLog << warn << GetDescriptor() << ": Event has zero cerenkov photons... skipped." << endl;
         return kFALSE;
+    }
+
+    if (fNumUsedPixels<3)
+    {
+        *fLog << warn << GetDescriptor() << ": Event has less than 3 used pixels... skipped." << endl;
+        return kFALSE;
+    }
 
     fMeanX /= fSize;                                 // [mm]
@@ -270 +289 @@
     //
     // calculate 2nd moments
-    // -------------------
-    //
-    float corrxx=0;                                  // [m^2]
-    float corrxy=0;                                  // [m^2]
-    float corryy=0;                                  // [m^2]
+    // ---------------------
+    //
+    // To make sure that the more complicated sum is evaluated as
+    // accurate as possible (because it is needed for more
+    // complicated calculations (see below) we calculate it using
+    // double prcision.
+    //
+    Double_t corrxx=0;                               // [m^2]
+    Double_t corrxy=0;                               // [m^2]
+    Double_t corryy=0;                               // [m^2]
 
     for (UInt_t i=0; i<npixevt; i++)
@@ -284 +308 @@
 
         const MGeomPix &gpix = geom[pix.GetPixId()];
-        const float dx = gpix.GetX() - fMeanX;       // [mm]
-        const float dy = gpix.GetY() - fMeanY;       // [mm]
-
-        const float nphot = pix.GetNumPhotons();     // [#phot]
+
+        const Float_t dx = gpix.GetX() - fMeanX;     // [mm]
+        const Float_t dy = gpix.GetY() - fMeanY;     // [mm]
+
+        const Float_t nphot = pix.GetNumPhotons();   // [#phot]
 
         corrxx += nphot * dx*dx;                     // [mm^2]
@@ -295 +320 @@
 
     //
+    // If corrxy=0 (which should never happen, because fSize>0) we
+    // cannot calculate Length and Width. The calculation failed
+    // and returnes kFALSE
+    //
+    if (corrxy==0)
+    {
+        *fLog << warn << GetDescriptor() << ": Event has CorrXY==0... skipped." << endl;
+        return kFALSE;
+    }
+
+    //
     // calculate the basic Hillas parameters: orientation and size of axes
     // -------------------------------------------------------------------
     //
-    const float d = (corryy - corrxx)/(corrxy*2);
-
-    fDelta = atan(d + sqrt(d*d + 1));
-
-    fCosDelta = cos(fDelta);   // need these in derived classes
-    fSinDelta = sin(fDelta);   // like MHillasExt
-
-    float axis1 = ( fCosDelta*fSinDelta*corrxy*2 + fCosDelta*fCosDelta*corrxx + fSinDelta*fSinDelta*corryy) / fSize; // [mm^2]
-    float axis2 = (-fCosDelta*fSinDelta*corrxy*2 + fSinDelta*fSinDelta*corrxx + fCosDelta*fCosDelta*corryy) / fSize; // [mm^2]
- 
+    // fDelta is the angle between the major axis of the ellipse and
+    //  the x axis. It is independent of the position of the ellipse
+    //  in the camera it has values between -pi/2 and pi/2 degrees
+    //
+    const Double_t d0   = corryy - corrxx;
+    const Double_t d1   = corrxy*2;
+    const Double_t d2   = d0 + sqrt(d0*d0 + d1*d1);
+    const Double_t tand = d2 / d1;
+
+    fDelta = atan(tand);
+
+    const Double_t s2 = tand*tand+1;
+    const Double_t s  = sqrt(s2);
+
+    fCosDelta =  1.0/s;   // need these in derived classes
+    fSinDelta = tand/s;   // like MHillasExt
+
+    Double_t axis1 = (tand*tand*corryy + d2 + corrxx)/s2/fSize;
+    Double_t axis2 = (tand*tand*corrxx - d2 + corryy)/s2/fSize;
+
+    //
+    // fLength^2 is the second moment along the major axis of the ellipse
+    // fWidth^2  is the second moment along the minor axis of the ellipse
+    //
+    // From the algorithm we get: fWidth <= fLength is always true
+    //
     // very small numbers can get negative by rounding
-    if (axis1 < 0) axis1=0;
-    if (axis2 < 0) axis2=0; 
-
-    fLength = sqrt(axis1);  // [mm]
-    fWidth  = sqrt(axis2);  // [mm]
+    //
+    fLength = axis1<0 ? 0 : sqrt(axis1);  // [mm]
+    fWidth  = axis2<0 ? 0 : sqrt(axis2);  // [mm]
 
     SetReadyToSave();
@@ -322 +372 @@
 // --------------------------------------------------------------------------
 //
+/*
 void MHillas::AsciiRead(ifstream &fin)
 {
@@ -331 +382 @@
     fin >> fMeanY;
 }
-
+*/
 // --------------------------------------------------------------------------
 /*

trunk/MagicSoft/Mars/manalysis/MHillas.h

@@ -15 +15 @@
 private:
     // for description see MHillas.cc
-    Float_t   fLength;   // [mm]        major axis of ellipse
-    Float_t   fWidth;    // [mm]        minor axis of ellipse
-    Float_t   fDelta;    // [rad]       angle of major axis with x-axis
-    Float_t   fSize;     // [#CerPhot]  sum of content of all pixels (number of Cherenkov photons)
-    Float_t   fMeanX;    // [mm]        x-coordinate of center of ellipse
-    Float_t   fMeanY;    // [mm]        y-coordinate of center of ellipse
+    Float_t fLength;        // [mm]        major axis of ellipse
+    Float_t fWidth;         // [mm]        minor axis of ellipse
+    Float_t fDelta;         // [rad]       angle of major axis with x-axis
+    Float_t fSize;          // [#CerPhot]  sum of content of all pixels (number of Cherenkov photons)
+    Float_t fMeanX;         // [mm]        x-coordinate of center of ellipse
+    Float_t fMeanY;         // [mm]        y-coordinate of center of ellipse
 
-    Float_t   fSinDelta; //! [1] sin of Delta (to be used in derived classes)
-    Float_t   fCosDelta; //! [1] cos of Delta (to be used in derived classes)
+    Float_t fSinDelta;      //! [1] sin of Delta (to be used in derived classes)
+    Float_t fCosDelta;      //! [1] cos of Delta (to be used in derived classes)
 
     Short_t fNumUsedPixels; // Number of pixels which survived the image cleaning
     Short_t fNumCorePixels; // number of core pixels
 
-    TEllipse *fEllipse;  //! Graphical Object to Display Ellipse
+    TEllipse *fEllipse;     //! Graphical Object to Display Ellipse
 
 protected:
@@ -62 +62 @@
     Int_t GetNumCorePixels() const { return fNumCorePixels; }
 
-    virtual void AsciiRead(ifstream &fin);
+    //virtual void AsciiRead(ifstream &fin);
     //virtual void AsciiWrite(ofstream &fout) const;
 

trunk/MagicSoft/Mars/manalysis/MHillasSrc.cc

@@ -19 +19 @@
 !   Author(s): Harald Kornmayer 1/2001
 !   Author(s): Rudolf Bock     10/2001 <mailto:Rudolf.Bock@cern.ch>
+!   Author(s): Wolfgang Wittek 06/2002 <mailto:wittek@mppmu.mpg.de>
 !
 !   Copyright: MAGIC Software Development, 2000-2002
@@ -32 +33 @@
 //
 //    source-dependent image parameters
-// fAlpha    angle between major axis and line source-to-center
-// fDist     distance from source to center of ellipse
+//
+// Version 1:
+// ----------
+//  fAlpha          angle between major axis and line source-to-center
+//  fDist           distance from source to center of ellipse
+//
+// Version 2:
+// ----------
+//  fHeadTail
+//
+// Version 3:
+// ----------
+//  fCosDeltaLenth  cosine of angle between d and a, where
+//                   - d is the vector from the source position to the
+//                     center of the ellipse
+//                   - a is a vector along the main axis of the ellipse
 //
 /////////////////////////////////////////////////////////////////////////////
@@ -57 +72 @@
 }
 
+void MHillasSrc::Reset()
+{
+    fDist          = -1;
+    fAlpha         =  0;
+    fHeadTail      =  0;
+    fCosDeltaAlpha =  0;
+}
+
 // --------------------------------------------------------------------------
 //
-//  calculation of source-dependent parameters
+//  Calculation of source-dependent parameters
+//  In case you don't call Calc from within an eventloop make sure, that
+//  you call the Reset member function before.
 //
-void MHillasSrc::Calc(const MHillas *hillas)
+Bool_t MHillasSrc::Calc(const MHillas *hillas)
 {
-    const Double_t mx = hillas->GetMeanX();        // [mm]
-    const Double_t my = hillas->GetMeanY();        // [mm]
+    fHillas = hillas;
 
-    const Double_t sx = mx - fSrcPos->GetX();      // [mm]
-    const Double_t sy = my - fSrcPos->GetY();      // [mm]
+    const Double_t mx   = GetMeanX();            // [mm]
+    const Double_t my   = GetMeanY();            // [mm]
 
-    const Double_t sd = sin(hillas->GetDelta());   // [1]
-    const Double_t cd = cos(hillas->GetDelta());   // [1]
+    const Double_t sx   = mx - fSrcPos->GetX();  // [mm]
+    const Double_t sy   = my - fSrcPos->GetY();  // [mm]
 
+    const Double_t sd   = sin(GetDelta());       // [1]
+    const Double_t cd   = cos(GetDelta());       // [1]
 
-    fHeadTail = cd*sx + sd*sy;                     // [mm]
-    fDist     = sqrt(sx*sx + sy*sy);               // [mm]
-    fAlpha    = atan((cd*sy - sd*sx)/fHeadTail);   // [rad]
-    fAlpha   *= kRad2Deg;                          // [deg]
+    const Double_t tand = tan(GetDelta());       // [1]
 
-    fHillas   = hillas;
+    fHeadTail = cd*sx + sd*sy;                   // [mm]
+
+    //
+    // Distance from source position to center of ellipse.
+    // If the distance is 0 distance, Alpha is not specified.
+    // The calculation has failed and returnes kFALSE.
+    //
+    Double_t dist = sqrt(sx*sx + sy*sy);         // [mm]
+
+    if (dist==0)
+    {
+        *fLog << warn << GetDescriptor() << ": Event has Dist==0... skipped." << endl;
+        return kFALSE;
+    }
+
+    fDist = dist;
+                                                 // [mm]
+    //
+    // Calculate Alpha and Cosda = cos(d,a)
+    // The sign of Cosda will be used for quantities containing 
+    // a head-tail information
+    //
+    const Double_t arg = (sy-tand*sx) / (fDist*sqrt(tand*tand+1));
+
+    fAlpha         = asin(arg)*kRad2Deg;        // [deg]
+    fCosDeltaAlpha = fHeadTail/fDist;           // [1]
 
     SetReadyToSave();
+
+    return kTRUE;
 } 
 
@@ -87 +137 @@
     *fLog << all;
     *fLog << "Source dependant Image Parameters (" << GetName() << ")" << endl;
-    *fLog << " - Dist     = " << fDist     << " mm"  << endl;
-    *fLog << " - Alpha    = " << fAlpha    << " deg" << endl;
-    *fLog << " - HeadTail = " << fHeadTail << " mm"  << endl;
+    *fLog << " - Dist          [mm]  = " << fDist << endl;
+    *fLog << " - Alpha         [deg] = " << fAlpha << endl;
+    *fLog << " - HeadTail      [mm]  = " << fHeadTail << endl;
+    *fLog << " - CosDeltaAlpha       = " << fCosDeltaAlpha << endl;
 }
 
@@ -96 +147 @@
 // overloaded MParContainer to read MHillasSrc from an ascii file
 //
+/*
 void MHillasSrc::AsciiRead(ifstream &fin)
 {
@@ -102 +154 @@
     fin >> fHeadTail;
 }
-
+*/
 // -----------------------------------------------------------------------
 //

trunk/MagicSoft/Mars/manalysis/MHillasSrc.h

@@ -14 +14 @@
     const MSrcPosCam *fSrcPos; //! Source position in the camere
 
-    Float_t fAlpha;    // [deg]  angle of major axis with vector to src
-    Float_t fDist;     // [mm]   distance from src to center of ellipse
-    Float_t fHeadTail; // [mm]
+    Float_t fAlpha;         // [deg]  angle of major axis with vector to src
+    Float_t fDist;          // [mm]   distance between src and center of ellipse
+    Float_t fHeadTail;      // [mm]
+    Float_t fCosDeltaAlpha; // [1]    cosine of angle between d and a
 
 public:
@@ -24 +25 @@
     const MSrcPosCam *GetSrcPos() const   { return fSrcPos; }
 
-    Float_t GetLength()   const  { return fHillas->GetLength(); }
-    Float_t GetWidth()    const  { return fHillas->GetWidth(); }
-    Float_t GetDelta()    const  { return fHillas->GetDelta(); }
-    Float_t GetSize()     const  { return fHillas->GetSize(); }
-    Float_t GetMeanX()    const  { return fHillas->GetMeanX(); }
-    Float_t GetMeanY()    const  { return fHillas->GetMeanY(); }
-    Float_t GetAlpha()    const  { return fAlpha; }
-    Float_t GetDist()     const  { return fDist; }
-    Float_t GetHeadTail() const  { return fHeadTail; }
+    void Reset();
+
+    Float_t GetLength()        const { return fHillas->GetLength(); }
+    Float_t GetWidth()         const { return fHillas->GetWidth(); }
+    Float_t GetDelta()         const { return fHillas->GetDelta(); }
+    Float_t GetSize()          const { return fHillas->GetSize(); }
+    Float_t GetMeanX()         const { return fHillas->GetMeanX(); }
+    Float_t GetMeanY()         const { return fHillas->GetMeanY(); }
+    Float_t GetAlpha()         const { return fAlpha; }
+    Float_t GetDist()          const { return fDist; }
+    Float_t GetHeadTail()      const { return fHeadTail; }
+    Float_t GetCosDeltaAlpha() const { return fCosDeltaAlpha; }
 
     void Print(Option_t *opt=NULL) const;
 
-    virtual void Calc(const MHillas *hillas);
+    virtual Bool_t Calc(const MHillas *hillas);
 
-    virtual void AsciiRead(ifstream &fin);
+    //virtual void AsciiRead(ifstream &fin);
     //virtual void AsciiWrite(ofstream &fout) const;
 
-    ClassDef(MHillasSrc, 2) // Container to hold source position dependant parameters
+    ClassDef(MHillasSrc, 3) // Container to hold source position dependant parameters
 };
 
-/*
- Version 1: fAlpha, fDist
- Version 2: fAlpha, fDist, fSign
- */
-
 #endif

trunk/MagicSoft/Mars/manalysis/MHillasSrcCalc.cc

@@ -91 +91 @@
 Bool_t MHillasSrcCalc::Process()
 {
-    fHillasSrc->Calc(fHillas);
-
-    return kTRUE;
+    return fHillasSrc->Calc(fHillas) ? kTRUE : kCONTINUE;
 }
 

trunk/MagicSoft/Mars/mgui/MineSweeper.h

@@ -41 +41 @@
     enum
     {
-//        kBlack     =  1,      // schwarz
-//        kWhite     = 10,
         kHidden    = 50,
         kIsVisible = BIT(15),
         kHasBomb   = BIT(16),
         kHasFlag   = BIT(17),
-        kUserBits  = 0x0000ff00, // 14-23 are allowed
-
+        kUserBits  = 0x0000ff00 // 14-23 are allowed
     };