Changeset 1434 for trunk/MagicSoft/Mars/manalysis/MHillas.cc

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/manalysis/MHillas.cc (modified) (14 diffs)

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;
 }
-
+*/
 // --------------------------------------------------------------------------
 /*