Changeset 7177

Timestamp:

07/08/05 11:49:41 (19 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• NEWS (modified) (1 diff)
• mimage/MHillas.cc (modified) (10 diffs)
• mimage/MHillasCalc.cc (modified) (1 diff)

trunk/MagicSoft/Mars/Changelog

@@ -43 +43 @@
    * mmuon/MMuonCalibPar.cc, mmuon/MMuonSearchPar.cc:
      - fixed Print output
+
+   * mimage/MHillas.cc:
+     - implemented plotting a estimate of disp
+     - fixed some units in output
+     - fixed the special case in which the shower is symmetric to
+       the axis
+
+   * mimage/MHillasCalc.cc:
+     - skipped printing the removed case 4 from MHillas::Calc
 
 

trunk/MagicSoft/Mars/NEWS

@@ -14 +14 @@
 
    - general: MH::SetPalette offers a lot of new palettes
+
+   - general: MHillas - the case of CorrXY==0 is now handled properly
 
    - showplot: got a new option to start a root interpreter, too

trunk/MagicSoft/Mars/mimage/MHillas.cc

@@ -159 +159 @@
     *fLog << " - Width          [mm]  = " << fWidth  << endl;
     *fLog << " - Delta          [deg] = " << fDelta*kRad2Deg << endl;
-    *fLog << " - Size           [1]   = " << fSize   << " #CherPhot"   << endl;
+    *fLog << " - Size           [phe] = " << fSize   << endl;
     *fLog << " - Meanx          [mm]  = " << fMeanX  << endl;
     *fLog << " - Meany          [mm]  = " << fMeanY  << endl;
@@ -182 +182 @@
     *fLog << " - Width          [deg] = " << fWidth *geom.GetConvMm2Deg() << endl;
     *fLog << " - Delta          [deg] = " << fDelta*kRad2Deg << endl;
-    *fLog << " - Size           [1]   = " << fSize   << " #CherPhot"   << endl;
+    *fLog << " - Size           [phe] = " << fSize << endl;
     *fLog << " - Meanx          [deg] = " << fMeanX *geom.GetConvMm2Deg() << endl;
     *fLog << " - Meany          [deg] = " << fMeanY *geom.GetConvMm2Deg() << endl;
@@ -215 +215 @@
     line.PaintLine(+w*fSinDelta+fMeanX, -w*fCosDelta+fMeanY,
                    -w*fSinDelta+fMeanX, +w*fCosDelta+fMeanY);
+
+    // Rough estimate for disp
+    const Double_t p = 362*(1-fWidth/fLength);
+
+    // Disp line
+    line.PaintLine( p*fCosDelta+20*fSinDelta+fMeanX,  p*fSinDelta-20*fCosDelta+fMeanY,
+                    p*fCosDelta-20*fSinDelta+fMeanX,  p*fSinDelta+20*fCosDelta+fMeanY);
+    line.PaintLine(-p*fCosDelta+20*fSinDelta+fMeanX, -p*fSinDelta-20*fCosDelta+fMeanY,
+                   -p*fCosDelta-20*fSinDelta+fMeanX, -p*fSinDelta+20*fCosDelta+fMeanY);
 
     TEllipse e(fMeanX, fMeanY, fLength, fWidth, 0, 360, fDelta*kRad2Deg+180);
@@ -259 +268 @@
 
     UInt_t numused = 0;
-/*
-    MSignalPix *pix = 0;
-
-    TIter Next(evt);
-    while ((pix=(MSignalPix*)Next()))*/
+
     for (UInt_t i=0; i<numpix; i++)
     {
@@ -273 +278 @@
             continue;
 
-        const MGeomPix &gpix = geom[i/*pix->GetPixId()*/];
+        const MGeomPix &gpix = geom[i];
 
         const Float_t nphot = pix.GetNumPhotons();
@@ -312 +317 @@
     Double_t corryy=0;                               // [m^2]
 
-    //Next.Reset();
-    //while ((pix=(MSignalPix*)Next()))
-    //{
     for (UInt_t i=0; i<numpix; i++)
     {
@@ -324 +326 @@
             continue;
 
-        const MGeomPix &gpix = geom[i/*pix->GetPixId()*/];
+        const MGeomPix &gpix = geom[i];
 
         const Float_t dx = gpix.GetX() - fMeanX;     // [mm]
@@ -337 +339 @@
 
     //
-    // If corrxy=0 (which should happen not very often, because fSize>0)
-    // we cannot calculate Length and Width.
-    // In reallity it is almost impossible to have a distribution of
-    // cerenkov photons in the used pixels which is exactly symmetric
-    // along one of the axis.
-    // It seems to be less than 0.1% of all events.
-    //
-    if (corrxy==0)
-        return 4;
-
-    //
     // calculate the basic Hillas parameters: orientation and size of axes
     // -------------------------------------------------------------------
@@ -360 +351 @@
     const Double_t d0    = corryy - corrxx;
     const Double_t d1    = corrxy*2;
-    const Double_t d2    = d0 + TMath::Sqrt(d0*d0 + d1*d1);
-    const Double_t tand  = d2 / d1;
+    const Double_t d2    = TMath::Sqrt(d0*d0 + d1*d1) + d0; // [0
+
+    const Double_t tand  = d2==0 ? 0 : d2 / d1; // Force 0 also if d1==0
     const Double_t tand2 = tand*tand;
 
-    fDelta = TMath::ATan(tand);
-
-    const Double_t s2 = tand2+1;
-    const Double_t s  = TMath::Sqrt(s2);
-
-    fCosDelta =  1.0/s;   // need these in derived classes
-    fSinDelta = tand/s;   // like MHillasExt
-
-    const Double_t axis1 = (tand2*corryy + d2 + corrxx)/s2/fSize;
-    const Double_t axis2 = (tand2*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
+    const Double_t s2    = tand2+1;
+    const Double_t s     = TMath::Sqrt(s2);
+
+    // Set default for the case in which the image is symmetric on the y-axis
+    fDelta    = TMath::Pi()/2;
+    fCosDelta = 0;
+    fSinDelta = 1;
+
+    Double_t axis1 = corryy;
+    Double_t axis2 = corrxx;
+
+    // This are all cases in which the image is not symmetric on the y-axis
+    if (d1!=0 || d2==0)
+    {
+        fDelta    = TMath::ATan(tand);
+
+        fCosDelta = 1.0 /s;   // need these in derived classes
+        fSinDelta = tand/s;   // like MHillasExt
+
+        axis1 = (tand2*corryy + d2 + corrxx)/s2;
+        axis2 = (tand2*corrxx - d2 + corryy)/s2;
+    }
+
+    //
+    // fLength^2 is the second moment along the major axis of the distribution
+    // fWidth^2  is the second moment along the minor axis of the distribution
     //
     // From the algorithm we get: fWidth <= fLength is always true
@@ -383 +387 @@
     // very small numbers can get negative by rounding
     //
-    fLength = axis1<0 ? 0 : TMath::Sqrt(axis1);  // [mm]
-    fWidth  = axis2<0 ? 0 : TMath::Sqrt(axis2);  // [mm]
+    fLength = axis1<0 ? 0 : TMath::Sqrt(axis1/fSize);  // [mm]
+    fWidth  = axis2<0 ? 0 : TMath::Sqrt(axis2/fSize);  // [mm]
 
     SetReadyToSave();

trunk/MagicSoft/Mars/mimage/MHillasCalc.cc

@@ -376 +376 @@
         PrintSkipped(fErrors[2], "Calculated Size == 0 (after cleaning)");
         PrintSkipped(fErrors[3], "Number of used pixels < 3");
-        PrintSkipped(fErrors[4], "CorrXY==0");
+        //PrintSkipped(fErrors[4], "CorrXY==0");
     }
     if (TestFlag(kCalcHillasSrc))