Changeset 3682 for trunk/MagicSoft/Mars

Timestamp:

04/08/04 19:52:13 (21 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (2 diffs)
• macros/readCT1.C (modified) (1 diff)
• macros/readMagic.C (modified) (1 diff)
• mfbase/MF.cc (modified) (7 diffs)
• mfbase/MFilterList.cc (modified) (1 diff)
• mfileio/MReadMarsFile.cc (modified) (1 diff)
• mfileio/MReadMarsFile.h (modified) (1 diff)
• mfileio/MReadTree.cc (modified) (1 diff)
• mfileio/MReadTree.h (modified) (1 diff)
• mfileio/MWriteRootFile.cc (modified) (1 diff)
• mhist/MHFalseSource.cc (modified) (19 diffs)
• mhist/MHFalseSource.h (modified) (3 diffs)
• mimage/MHHillasExt.cc (modified) (11 diffs)
• mimage/MHHillasExt.h (modified) (1 diff)
• mimage/MHNewImagePar.cc (modified) (7 diffs)
• mimage/MHNewImagePar.h (modified) (2 diffs)
• mimage/MHillasCalc.cc (modified) (2 diffs)
• mimage/MHillasExt.cc (modified) (5 diffs)
• mimage/MNewImagePar.cc (modified) (5 diffs)
• mimage/MNewImagePar.h (modified) (2 diffs)
• mmain/MEventDisplay.cc (modified) (2 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -18 +18 @@
 
                                                  -*-*- END OF LINE -*-*-
+ 2004/04/08: Thomas Bretz
+
+   * mfbase/MF.cc:
+     - optimized output
+
+   * mfbase/MFilterList.cc:
+     - added some comments
+
+   * mfileio/MReadMarsFile.[h,cc], mfileio/MReadTree.[h,cc]:
+     - added SortFiled-option -- RAQUEL, this is for you!
+
+   * mfileio/MWriteRootFile.cc:
+     - class AddSerialNumber in AddContainer
+
+   * mhist/MHFalseSource.[h,cc]:
+     - optimized setting of hist maximum
+     - added more source dependant cuts
+     - changed display layout
+     - scale number of entries to correct for different acceptance
+       (not perfect, but the best I can currently do)
+
+   * mimage/MHHillasExt.[h,cc]:
+     - added new histogram for fMaxDist
+
+   * mimage/MHNewImagePar.[h,cc]:
+     - added new histograms for Used/CoreArea
+
+   * mimage/MHillasCalc.cc:
+     - optimized output
+
+   * mimage/MHillasExt.cc:
+     - fMaxDist got a sign
+
+   * mimage/MNewImagePar.[h,cc]:
+     - enhanced comments
+     - added new Print() member function
+
+   * macros/readCT1.C, macros/readMagic.C, mmain/MEventDisplay.cc:
+     - forward geomcam to newimgepar.Print()
+
+
 
  2004/04/08: Markus Gaug
@@ -40 +81 @@
      - removed warning about low-gain saturation of Blind pixel
 
+
+
  2004/04/07: Markus Gaug
+
    * mcalib/MHGausEvents.[h,cc]
      - added fBlackout events
+
 
 

trunk/MagicSoft/Mars/macros/readCT1.C

@@ -125 +125 @@
         hillas.Print(*geomcam);
         hillasext.Print(*geomcam);
-        newimgpar.Print();
+        newimgpar.Print(*geomcam);
 
         if (!HandleInput())

trunk/MagicSoft/Mars/macros/readMagic.C

@@ -100 +100 @@
     ((MHillas*)fParList->FindObject("MHillas"))->Print(*geom);
     ((MHillasExt*)fParList->FindObject("MHillasExt"))->Print(*geom);
-    fParList->FindObject("MNewImagePar")->Print();
+    ((MNewImagePar*)fParList->FindObject("MNewImagePar"))->Print(*geom);
 
     return HandleInput();

trunk/MagicSoft/Mars/mfbase/MF.cc

@@ -41 +41 @@
 //   "MHillas.fWidth<0.5 && MHillas.fLength<0.6"
 //
-// You can also use brackets:
+// You can also use parantheses:
 //   "MHillas.fSize>200 || (MHillas.fWidth<0.5 && MHillas.fLength<0.6)"
 //
 // If you want to use mathematic expressions (as defined in MDataChain)
-// you must encapsulate it in {}-Brackets, eg:
+// you must encapsulate it in {}-parantheses, eg:
 //   "{log10(MHillas.fSize)}>3"
 //
@@ -60 +60 @@
 //
 //
-// Warning: There is no priority rule build in. So better use brackets
+// Warning: There is no priority rule build in. So better use parantheses
 //   to get correct results. The rule is parsed/evaluated from the left
 //   to the right, which means:
@@ -120 +120 @@
     fTitle = title ? title : gsDefTitle.Data();
 
-    *fLog << inf << "Trying to resolve filter rule..." << endl;
+    *fLog << inf << "Trying to resolve filter rule... " << flush;
     if (!(fF=ParseString(text, 1)))
     {
@@ -127 +127 @@
     }
 
-    *fLog << inf << endl;
-    *fLog << "Using Filter rule " << fF->GetName();
-    *fLog << " for " << fName << ":" << endl;
+    *fLog << inf << "found: ";
     fF->Print();
-    *fLog << endl << endl;
+    *fLog << endl;
 }
 
@@ -178 +176 @@
     {
         *fLog << err << dbginf << "Syntax Error: Two coditional signs found in " << txt << endl;
-        *fLog << "Currently you have to enclose all conditions in brackets, like: \"(x<y) && (z<5)\"" << endl;
+        *fLog << "Currently you have to enclose all conditions in parantheses, like: \"(x<y) && (z<5)\"" << endl;
         return NULL;
     }
@@ -226 +224 @@
             {
                 //
-                // Search for the corresponding bracket
+                // Search for the corresponding parantheses
                 //
                 Int_t first=1;
@@ -250 +248 @@
                 //
                 // Make a copy of the 'interieur' and delete the substringä
-                // including the brackets
+                // including the parantheses
                 //
                 TString sub = txt(1, first-1);

trunk/MagicSoft/Mars/mfbase/MFilterList.cc

@@ -79 +79 @@
 //
 //   Options:
-//      and, &   : is a bitwise and
-//      or, |    : is a bitwise or
-//      xor, ^   : is a bitwise exclusive or
+//      and,  &  : is a bitwise and
+//      or,   |  : is a bitwise or
+//      xor,  ^  : is a bitwise exclusive or
 //      land, && : is a logical and
-//      lor, ||  : is a logical or
+//      lor,  || : is a logical or
+//
+//   The bitwise operators are implemented for convinience. They may not
+//   make much sense. Maybe IsExpressionTrue should change its return
+//   type from Bool_t to Int_t.
 //
 MFilterList::MFilterList(const char *type, const char *name, const char *title)

trunk/MagicSoft/Mars/mfileio/MReadMarsFile.cc

@@ -116 +116 @@
 // --------------------------------------------------------------------------
 //
+//  Sort the files by their file-names
+//
+void MReadMarsFile::SortFiles()
+{
+    fRun->SortFiles();
+    MReadTree::SortFiles();
+}
+
+// --------------------------------------------------------------------------
+//
 //  This overload MReadTree::Notify. Before the MReadTree Notify
 //  TObjects are called the RunHeaders of the next files are read.

trunk/MagicSoft/Mars/mfileio/MReadMarsFile.h

@@ -24 +24 @@
     ~MReadMarsFile();
 
+    void SortFiles();
+
     Int_t AddFile(const char *fname, Int_t entries=-1);
 

trunk/MagicSoft/Mars/mfileio/MReadTree.cc

@@ -358 +358 @@
 // --------------------------------------------------------------------------
 //
+//  Sort the files by their file-names
+//
+void MReadTree::SortFiles()
+{
+    fChain->GetListOfFiles()->Sort();
+}
+
+// --------------------------------------------------------------------------
+//
 //  This function is called if Branch choosing method should get enabled.
 //  Branch choosing means, that only the enabled branches are read into

trunk/MagicSoft/Mars/mfileio/MReadTree.h

@@ -48 +48 @@
     ~MReadTree();
 
+    virtual void SortFiles();
+
     void   DisableAutoScheme() { fAutoEnable = kFALSE; }
     void   EnableBranch(const char *name);

trunk/MagicSoft/Mars/mfileio/MWriteRootFile.cc

@@ -198 +198 @@
     // add the entry to the list.
     //
-    MRootFileBranch *entry = new MRootFileBranch(cname, tname, must);
+    MRootFileBranch *entry = new MRootFileBranch(AddSerialNumber(cname), tname, must);
     fBranches.AddLast(entry);
 
     if (tname && tname[0])
-        AddToBranchList(Form("%s.%s", cname, tname));
+        AddToBranchList(Form("%s.%s", AddSerialNumber(cname), tname));
 }
 

trunk/MagicSoft/Mars/mhist/MHFalseSource.cc

@@ -46 +46 @@
 // pads) you can change the AlphaCut 'online'
 //
+// Each Z-Projection (Alpha-histogram) is scaled such, that the number
+// of entries fits the maximum number of entries in all Z-Projections.
+// This should correct for the different acceptance in the center
+// and at the vorder of the camera. This, however, produces more noise
+// at the border.
+//
 // Here is a slightly simplified version of the algorithm:
 // ------------------------------------------------------
@@ -139 +145 @@
 //
 MHFalseSource::MHFalseSource(const char *name, const char *title)
-    : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1),
-      fAlphaCut(12.5), fBgMean(55), fDistMin(-1), fDistMax(-1)
+    : fTime(0), fPointPos(0), fObservatory(0), fMm2Deg(-1), fAlphaCut(12.5),
+    fBgMean(55), fMinDist(-1), fMaxDist(-1), fMinLD(-1), fMaxLD(-1)
 {
     //
@@ -155 +161 @@
     fHist.SetYTitle("y [\\circ]");
     fHist.SetZTitle("\\alpha [\\circ]");
+}
+
+void MHFalseSource::MakeSymmetric(TH1 *h)
+{
+    const Float_t min = TMath::Abs(h->GetMinimum());
+    const Float_t max = TMath::Abs(h->GetMaximum());
+
+    const Float_t absmax = TMath::Max(min, max)*1.002;
+
+    h->SetMaximum( absmax);
+    h->SetMinimum(-absmax);
 }
 
@@ -337 +354 @@
             //        and/or MTaskList
             // Source dependant distance cut
-            if (fDistMin>0 && hsrc.GetDist()*fMm2Deg<fDistMin)
+            if (fMinDist>0 && hsrc.GetDist()*fMm2Deg<fMinDist)
                 continue;
-
-            if (fDistMax>0 && hil->GetLength()>fDistMax*hsrc.GetDist())
+            if (fMaxDist>0 && hsrc.GetDist()*fMm2Deg>fMaxDist)
+                continue;
+
+            if (fMaxLD>0 && hil->GetLength()>fMaxLD*hsrc.GetDist())
+                continue;
+            if (fMinLD>0 && hil->GetLength()<fMinLD*hsrc.GetDist())
                 continue;
 
@@ -358 +379 @@
 // the same number of bins than for on-data
 //
-void MHFalseSource::ProjectOff(TH2D *h2)
+void MHFalseSource::ProjectOff(TH2D *h2, TH2D *all)
 {
     TAxis &axe = *fHist.GetZaxis();
@@ -379 +400 @@
     // Move contents from projection to h2
     h2->Reset();
-    h2->Add(p);
+    h2->Add(p, all->GetMaximum());
+    h2->Divide(all);
 
     // Delete p
@@ -393 +415 @@
 // range (0, fAlphaCut)
 //
-void MHFalseSource::ProjectOn(TH2D *h3)
+void MHFalseSource::ProjectOn(TH2D *h3, TH2D *all)
 {
     TAxis &axe = *fHist.GetZaxis();
@@ -410 +432 @@
     // Move contents from projection to h3
     h3->Reset();
-    h3->Add(p);
+    h3->Add(p, all->GetMaximum());
+    h3->Divide(all);
+
+    // Delete p
     delete p;
 
@@ -454 +479 @@
     TH2D* h3;
     TH2D* h4;
+    TH2D* h5;
 
     // Update projection of all-events
-    padsave->GetPad(1)->cd(3);
+    padsave->GetPad(2)->cd(3);
     if ((h0 = (TH2D*)gPad->FindObject("Alpha_yx_all")))
         ProjectAll(h0);
 
+    // Update projection of on-events
+    padsave->GetPad(1)->cd(1);
+    if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on")))
+        ProjectOn(h3, h0);
+
     // Update projection of off-events
-    padsave->GetPad(2)->cd(1);
+    padsave->GetPad(1)->cd(3);
     if ((h2 = (TH2D*)gPad->FindObject("Alpha_yx_off")))
-        ProjectOff(h2);
-
-    // Update projection of on-events
+        ProjectOff(h2, h0);
+
     padsave->GetPad(2)->cd(2);
-    if ((h3 = (TH2D*)gPad->FindObject("Alpha_yx_on")))
-        ProjectOn(h3);
+    if ((h5 = (TH2D*)gPad->FindObject("Alpha_yx_diff")))
+    {
+        h5->Add(h3, h2, -1);
+        MakeSymmetric(h5);
+    }
 
     // Update projection of significance
-    padsave->GetPad(2)->cd(3);
+    padsave->GetPad(1)->cd(2);
     if (h2 && h3 && (h4 = (TH2D*)gPad->FindObject("Alpha_yx_sig")))
     {
@@ -483 +516 @@
         Int_t max = h4->GetBin(nx, ny);
 
-        for (int ix=0; ix<nx; ix++)
-            for (int iy=0; iy<ny; iy++)
+        for (int ix=1; ix<=nx; ix++)
+            for (int iy=1; iy<=ny; iy++)
             {
-                const Int_t n = h4->GetBin(ix+1, iy+1);
+                const Int_t n = h4->GetBin(ix, iy);
 
                 const Double_t s = h3->GetBinContent(n);
@@ -503 +536 @@
             }
 
+        MakeSymmetric(h4);
+
         // Update projection of 'the best alpha-plot'
-        padsave->GetPad(1)->cd(1);
+        padsave->GetPad(2)->cd(1);
         if ((h1 = (TH1D*)gPad->FindObject("Alpha")) && max>0)
         {
@@ -560 +595 @@
     TObject *catalog = GetCatalog();
 
-    // draw the 2D histogram Sigmabar versus Theta
+    // Initialize upper part
     pad->cd(1);
     gPad->SetBorderMode(0);
     gPad->Divide(3, 1);
-    delete pad->GetPad(1)->GetPad(2);
-
-    pad->GetPad(1)->cd(3);
-    gPad->SetBorderMode(0);
-    TH1 *h0 = fHist.Project3D("yx_all");
-    h0->SetDirectory(NULL);
-    h0->SetXTitle(fHist.GetXaxis()->GetTitle());
-    h0->SetYTitle(fHist.GetYaxis()->GetTitle());
-    h0->Draw("colz");
-    h0->SetBit(kCanDelete);
-    catalog->Draw("mirror same");
-
+
+    // PAD #1
     pad->GetPad(1)->cd(1);
-    gPad->SetBorderMode(0);
-
-    TH1 *h1 = fHist.ProjectionZ("Alpha");
-    h1->SetDirectory(NULL);
-    h1->SetTitle("Distribution of \\alpha");
-    h1->SetXTitle(fHist.GetZaxis()->GetTitle());
-    h1->SetYTitle("Counts");
-    h1->Draw(opt);
-    h1->SetBit(kCanDelete);
-
-    pad->cd(2);
-    gPad->SetBorderMode(0);
-    gPad->Divide(3, 1);
-
-    pad = gPad;
-
-    pad->cd(1);
-    gPad->SetBorderMode(0);
-    fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
-    TH1 *h2 = fHist.Project3D("yx_off");
-    h2->SetDirectory(NULL);
-    h2->SetXTitle(fHist.GetXaxis()->GetTitle());
-    h2->SetYTitle(fHist.GetYaxis()->GetTitle());
-    h2->Draw("colz");
-    h2->SetBit(kCanDelete);
-    catalog->Draw("mirror same");
-
-    pad->cd(2);
     gPad->SetBorderMode(0);
     fHist.GetZaxis()->SetRangeUser(0,fAlphaCut);
@@ -616 +613 @@
     catalog->Draw("mirror same");
 
-    pad->cd(3);
+    // PAD #2
+    pad->GetPad(1)->cd(2);
     gPad->SetBorderMode(0);
     fHist.GetZaxis()->SetRange(0,0);
@@ -627 +625 @@
     h4->Draw("colz");
     h4->SetBit(kCanDelete);
+    catalog->Draw("mirror same");
+
+    // PAD #3
+    pad->GetPad(1)->cd(3);
+    gPad->SetBorderMode(0);
+    fHist.GetZaxis()->SetRangeUser(fBgMean-fAlphaCut/2, fBgMean+fAlphaCut/2);
+    TH1 *h2 = fHist.Project3D("yx_off");
+    h2->SetDirectory(NULL);
+    h2->SetXTitle(fHist.GetXaxis()->GetTitle());
+    h2->SetYTitle(fHist.GetYaxis()->GetTitle());
+    h2->Draw("colz");
+    h2->SetBit(kCanDelete);
+    catalog->Draw("mirror same");
+
+    // Initialize lower part
+    pad->cd(2);
+    gPad->SetBorderMode(0);
+    gPad->Divide(3, 1);
+
+    // PAD #4
+    pad->GetPad(2)->cd(1);
+    gPad->SetBorderMode(0);
+    TH1 *h1 = fHist.ProjectionZ("Alpha");
+    h1->SetDirectory(NULL);
+    h1->SetTitle("Distribution of \\alpha");
+    h1->SetXTitle(fHist.GetZaxis()->GetTitle());
+    h1->SetYTitle("Counts");
+    h1->Draw(opt);
+    h1->SetBit(kCanDelete);
+
+    // PAD #5
+    pad->GetPad(2)->cd(2);
+    gPad->SetBorderMode(0);
+    TH1 *h5 = (TH1*)h3->Clone("Alpha_yx_diff");
+    h5->Add(h2, -1);
+    h5->SetTitle("Difference of on- and off-distribution");
+    h5->SetDirectory(NULL);
+    h5->Draw("colz");
+    h5->SetBit(kCanDelete);
+    catalog->Draw("mirror same");
+
+    // PAD #6
+    pad->GetPad(2)->cd(3);
+    gPad->SetBorderMode(0);
+    TH1 *h0 = fHist.Project3D("yx_all");
+    h0->SetDirectory(NULL);
+    h0->SetXTitle(fHist.GetXaxis()->GetTitle());
+    h0->SetYTitle(fHist.GetYaxis()->GetTitle());
+    h0->Draw("colz");
+    h0->SetBit(kCanDelete);
     catalog->Draw("mirror same");
 }
@@ -729 +777 @@
                        Form("Fit Region: Signal<%.1f\\circ, %.1f\\circ<Bg<%.1f\\circ",
                             sigmax, bgmin, bgmax));
-    hists->SetNameTitle("Excess",     Form("Number of excess events for \\alpha<%.0f\\circ", sigint));
-    histb->SetNameTitle("Background", Form("Number of background events for \\alpha<%.0f\\circ", sigint));
+    hists->SetName("Excess");
+    histb->SetName("Background");
     hist->SetXTitle(fHist.GetXaxis()->GetTitle());
     hists->SetXTitle(fHist.GetXaxis()->GetTitle());
@@ -776 +824 @@
     *fLog << "Polynom order:  " << (int)polynom << endl;
     *fLog << "Fitting False Source Plot..." << flush;
+
+    TH1 *h0 = fHist.Project3D("yx_entries");
+    Float_t entries = h0->GetMaximum();
+    delete h0;
 
     TH1 *h=0;
@@ -790 +842 @@
             if (alpha0==0)
                 continue;
+
+            h->Scale(entries/h->GetEntries());
 
             if (alpha0>maxalpha0)
@@ -887 +941 @@
     h0b.GetXaxis()->SetRangeUser(0, maxalpha0*1.5);
 
+    hists->SetTitle(Form("Excess events for \\alpha<%.0f\\circ (N_{max}=%d)", sigint, (int)hists->GetMaximum()));
+    histb->SetTitle(Form("Background events for \\alpha<%.0f\\circ", sigint));
+
     //hists->SetMinimum(GetMinimumGT(*hists));
     histb->SetMinimum(GetMinimumGT(*histb));
+
+    MakeSymmetric(hists);
+    MakeSymmetric(hist);
 
     clk.Stop();
@@ -950 +1010 @@
 
         TH1 *result = fHist.ProjectionZ("AlphaFit", maxx, maxx, maxy, maxy);
+        result->Scale(entries/h->GetEntries());
+
         result->SetDirectory(NULL);
         result->SetNameTitle("Result \\alpha", title);

trunk/MagicSoft/Mars/mhist/MHFalseSource.h

@@ -29 +29 @@
     Float_t fBgMean;             // Background mean
 
-    Float_t fDistMin;            // Min dist
-    Float_t fDistMax;            // Maximum distance in percent of dist
+    Float_t fMinDist;            // Min dist
+    Float_t fMaxDist;            // Max dist
+
+    Float_t fMinLD;              // Minimum distance in percent of dist
+    Float_t fMaxLD;              // Maximum distance in percent of dist
 
     TH3D    fHist;               // Alpha vs. x and y
@@ -40 +43 @@
     void Modified();
 
-    void ProjectOff(TH2D *h);
-    void ProjectOn(TH2D *h);
     void ProjectAll(TH2D *h);
+    void ProjectOff(TH2D *h, TH2D *all);
+    void ProjectOn(TH2D *h, TH2D *all);
 
     TObject *GetCatalog();
+
+    void MakeSymmetric(TH1 *h);
 
 public:
@@ -57 +62 @@
     void FitSignificanceStd() { FitSignificance(); } //*MENU*
 
-    void SetDistMin(Float_t dist)  { fDistMin = dist;  } // Absolute minimum distance
-    void SetDistMax(Float_t ratio) { fDistMax = ratio; } // Maximum ratio between length/dist
+    void SetMinDist(Float_t dist) { fMinDist = dist; } // Absolute minimum distance
+    void SetMaxDist(Float_t dist) { fMaxDist = dist; } // Absolute maximum distance
+    void SetMinLD(Float_t ratio)  { fMinLD = ratio; }  // Minimum ratio between length/dist
+    void SetMaxLD(Float_t ratio)  { fMaxLD = ratio; }  // Maximum ratio between length/dist
 
     void SetAlphaCut(Float_t alpha); //*MENU*

trunk/MagicSoft/Mars/mimage/MHHillasExt.cc

@@ -16 +16 @@
 !
 !
-!   Author(s): Thomas Bretz  2001 <mailto:tbretz@uni-sw.gwdg.de>
-!
-!   Copyright: MAGIC Software Development, 2000-2002
+!   Author(s): Thomas Bretz, 2001 <mailto:tbretz@astro.uni-wuerzburg.de>
+!
+!   Copyright: MAGIC Software Development, 2000-2004
 !
 !
@@ -75 +75 @@
     fHM3Long.SetDirectory(NULL);
     fHM3Trans.SetDirectory(NULL);
+    fHMaxDist.SetDirectory(NULL);
 
     fHAsym.UseCurrentStyle();
     fHM3Long.UseCurrentStyle();
     fHM3Trans.UseCurrentStyle();
+    fHMaxDist.UseCurrentStyle();
 
     fHAsym.SetName("Asymmetry");
     fHM3Long.SetName("3rd Mom Long");
     fHM3Trans.SetName("3rd Mom Trans");
+    fHMaxDist.SetName("Max Dist");
 
     fHAsym.SetTitle("Asymmetry");
     fHM3Long.SetTitle("3^{rd} Moment Longitudinal");
     fHM3Trans.SetTitle("3^{rd} Moment Transverse");
+    fHMaxDist.SetTitle("Distance of max distant pixel");
 
     fHAsym.SetXTitle("Asym [mm]");
     fHM3Long.SetXTitle("3^{rd} M_{l} [mm]");
     fHM3Trans.SetXTitle("3^{rd} M_{t} [mm]");
+    fHMaxDist.SetXTitle("D_{max} [mm]");
 
     fHAsym.SetYTitle("Counts");
     fHM3Long.SetYTitle("Counts");
     fHM3Trans.SetYTitle("Counts");
+    fHMaxDist.SetYTitle("Counts");
 
     fHAsym.SetFillStyle(4000);
     fHM3Long.SetFillStyle(4000);
     fHM3Trans.SetFillStyle(4000);
+    fHMaxDist.SetFillStyle(4000);
 
     fHM3Trans.SetLineColor(kBlue);
@@ -104 +111 @@
     MBinning bins;
 
-    bins.SetEdges(101, -326, 326);
+    bins.SetEdges(51, -326, 326);
     bins.Apply(fHM3Long);
     bins.Apply(fHM3Trans);
 
-    bins.SetEdges(101, -593, 593);
+    bins.SetEdges(51, -593, 593);
     bins.Apply(fHAsym);
+
+    bins.SetEdges(101, 0, 593);
+    bins.Apply(fHMaxDist);
 }
 
@@ -143 +153 @@
     ApplyBinning(*plist, "M3Long",  &fHM3Long);
     ApplyBinning(*plist, "M3Trans", &fHM3Trans);
+    ApplyBinning(*plist, "MaxDist", &fHMaxDist);
 
     return kTRUE;
@@ -161 +172 @@
     fHM3Long.Fill(scale*fHillasExt->GetM3Long(), w);
     fHM3Trans.Fill(scale*fHillasExt->GetM3Trans(), w);
-    //fHAsymna.Fill(scale*ext.GetAsymna());
-    //fHAsym0.Fill(scale*ext.GetAsym0());
+    fHMaxDist.Fill(TMath::Abs(scale*fHillasExt->GetMaxDist()), w);
 
     return kTRUE;
@@ -187 +197 @@
     MH::ScaleAxis(&fHM3Long,  scale);
     MH::ScaleAxis(&fHM3Trans, scale);
+    MH::ScaleAxis(&fHMaxDist, scale);
 
     if (mmscale)
@@ -193 +204 @@
         fHM3Long.SetXTitle("3^{rd} M_{l} [mm]");
         fHM3Trans.SetXTitle("3^{rd} M_{t} [mm]");
+        fHMaxDist.SetXTitle("D_{max} [mm]");
     }
     else
@@ -199 +211 @@
         fHM3Long.SetXTitle("3^{rd} M_{l} [\\circ]");
         fHM3Trans.SetXTitle("3^{rd} M_{t} [\\circ]");
+        fHMaxDist.SetXTitle("D_{max} [\\circ]");
     }
 
@@ -238 +251 @@
     AppendPad("");
 
-    pad->Divide(2, 1);
+    pad->Divide(2, 2);
 
     pad->cd(1);
@@ -244 +257 @@
     MH::DrawSame(fHM3Long, fHM3Trans, "3^{rd} Moments");
 
+    pad->cd(3);
+    gPad->SetBorderMode(0);
+    fHAsym.Draw();
+
     pad->cd(2);
     gPad->SetBorderMode(0);
-    fHAsym.Draw();
+    fHMaxDist.Draw();
+
+    delete pad->GetPad(4);
 
     pad->Modified();
@@ -260 +279 @@
     if (name.Contains("M3Trans", TString::kIgnoreCase))
         return &fHM3Trans;
+    if (name.Contains("MaxDist", TString::kIgnoreCase))
+        return &fHMaxDist;
 
     return NULL;

trunk/MagicSoft/Mars/mimage/MHHillasExt.h

@@ -19 +19 @@
     TH1F fHM3Long;  // [mm]    3rd moment (e-weighted) along major axis
     TH1F fHM3Trans; // [mm]    3rd moment (e-weighted) along minor axis
+    TH1F fHMaxDist; // [mm]    Distance between shower center maximum distant pixel
 
     Float_t fMm2Deg;

trunk/MagicSoft/Mars/mimage/MHNewImagePar.cc

@@ -19 +19 @@
 !   Author(s): Thomas Bretz, 04/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
 !
-!   Copyright: MAGIC Software Development, 2000-2003
+!   Copyright: MAGIC Software Development, 2000-2004
 !
 !
@@ -95 +95 @@
     fHistCorePix.SetFillStyle(4000);
 
+    fHistUsedArea.SetName("UsedArea");
+    fHistUsedArea.SetTitle("Area of used pixels");
+    fHistUsedArea.SetXTitle("Area [m^2]");
+    fHistUsedArea.SetYTitle("Counts");
+    fHistUsedArea.SetDirectory(NULL);
+    fHistUsedArea.UseCurrentStyle();
+    fHistUsedArea.SetLineColor(kBlue);
+    fHistUsedArea.SetFillStyle(4000);
+
+    fHistCoreArea.SetName("CoreArea");
+    fHistCoreArea.SetTitle("Area of core pixels");
+    fHistCoreArea.SetXTitle("Area [m^2]");
+    fHistCoreArea.SetYTitle("Counts");
+    fHistCoreArea.SetDirectory(NULL);
+    fHistCoreArea.UseCurrentStyle();
+    fHistCoreArea.SetLineColor(kBlack);
+    fHistCoreArea.SetFillStyle(4000);
+
     fHistConc.SetDirectory(NULL);
     fHistConc1.SetDirectory(NULL);
@@ -121 +139 @@
     bins.Apply(fHistConc1);
 
-    bins.SetEdges(150, 0, 150);
+    bins.SetEdges(75, 0, 150);
     bins.Apply(fHistUsedPix);
     bins.Apply(fHistCorePix);
+
+    bins.SetEdges(75, 0, 0.249);
+    bins.Apply(fHistUsedArea);
+    bins.Apply(fHistCoreArea);
 }
 
@@ -136 +158 @@
     ApplyBinning(*plist, "Leakage", &fHistLeakage2);
 
-    ApplyBinning(*plist, "Pixels", &fHistUsedPix);
-    ApplyBinning(*plist, "Pixels", &fHistCorePix);
-
-    ApplyBinning(*plist, "Conc",   &fHistConc);
-    ApplyBinning(*plist, "Conc1",  &fHistConc1);
+    ApplyBinning(*plist, "Pixels",  &fHistUsedPix);
+    ApplyBinning(*plist, "Pixels",  &fHistCorePix);
+
+    ApplyBinning(*plist, "Area",    &fHistUsedArea);
+    ApplyBinning(*plist, "Area",    &fHistCoreArea);
+
+    ApplyBinning(*plist, "Conc",    &fHistConc);
+    ApplyBinning(*plist, "Conc1",   &fHistConc1);
 
     return kTRUE;
@@ -166 +191 @@
     fHistCorePix.Fill(h.GetNumCorePixels(), w);
 
+    fHistUsedArea.Fill(h.GetUsedArea()/1000000, w);
+    fHistCoreArea.Fill(h.GetCoreArea()/1000000, w);
+
     fHistConc.Fill(h.GetConc(), w);
     fHistConc1.Fill(h.GetConc1(), w);
@@ -203 +231 @@
     pad->cd(4);
     gPad->SetBorderMode(0);
+    MH::DrawSame(fHistCoreArea, fHistUsedArea, "Area of core/used Pixels");
 
     pad->Modified();
@@ -222 +251 @@
     if (name.Contains("CorePix", TString::kIgnoreCase))
         return &fHistCorePix;
+    if (name.Contains("UsedArea", TString::kIgnoreCase))
+        return &fHistUsedArea;
+    if (name.Contains("CoreArea", TString::kIgnoreCase))
+        return &fHistCoreArea;
 
     return NULL;

trunk/MagicSoft/Mars/mimage/MHNewImagePar.h

@@ -14 +14 @@
 {
 private:
-    TH1F fHistLeakage1; //
-    TH1F fHistLeakage2; //
+    TH1F fHistLeakage1;  //
+    TH1F fHistLeakage2;  //
 
-    TH1F fHistUsedPix;  // Number of used pixels
-    TH1F fHistCorePix;  // Number of core pixels
+    TH1F fHistUsedPix;   // Number of used pixels
+    TH1F fHistCorePix;   // Number of core pixels
 
-    TH1F fHistConc;     // [ratio] concentration ratio: sum of the two highest pixels / fSize
-    TH1F fHistConc1;    // [ratio] concentration ratio: sum of the highest pixel / fSize
+    TH1F fHistUsedArea;  // Area of used pixels
+    TH1F fHistCoreArea;  // Area of core pixels
+
+    TH1F fHistConc;      // [ratio] concentration ratio: sum of the two highest pixels / fSize
+    TH1F fHistConc1;     // [ratio] concentration ratio: sum of the highest pixel / fSize
 
 public:
@@ -31 +34 @@
     TH1 *GetHistByName(const TString name);
 
-    TH1F &GetHistLeakage1() { return fHistLeakage1; }
-    TH1F &GetHistLeakage2() { return fHistLeakage2; }
+    TH1F &GetHistLeakage1()  { return fHistLeakage1; }
+    TH1F &GetHistLeakage2()  { return fHistLeakage2; }
 
-    TH1F &GetHistUsedPix()  { return fHistUsedPix; }
-    TH1F &GetHistCorePix()  { return fHistCorePix; }
+    TH1F &GetHistUsedPix()   { return fHistUsedPix; }
+    TH1F &GetHistCorePix()   { return fHistCorePix; }
 
-    TH1F &GetHistConc()     { return fHistConc; }
-    TH1F &GetHistConc1()    { return fHistConc1; }
+    TH1F &GetHistUsedArea()  { return fHistUsedArea; }
+    TH1F &GetHistCoreArea()  { return fHistCoreArea; }
+
+    TH1F &GetHistConc()      { return fHistConc; }
+    TH1F &GetHistConc1()     { return fHistConc1; }
 
     void Draw(Option_t *opt=NULL);

trunk/MagicSoft/Mars/mimage/MHillasCalc.cc

@@ -196 +196 @@
     *fLog << " " << setw(7) << fErrors[i] << " (";
     *fLog << setw(3) << (int)(100.*fErrors[i]/GetNumExecutions());
-    *fLog << "%) Evts skipped due to: " << str << endl;
+    *fLog << "%) Evts skipped: " << str << endl;
 }
 
@@ -212 +212 @@
     *fLog << GetDescriptor() << " execution statistics:" << endl;
     *fLog << dec << setfill(' ');
-    PrintSkipped(1, "Event has less than 3 pixels\n                                     (before image cleaning)");
-    PrintSkipped(2, "Calculated Size == 0\n                                     (no pixels survived image cleaning)");
+    PrintSkipped(1, "Less than 3 pixels (before cleaning)");
+    PrintSkipped(2, "Calculated Size == 0 (after cleaning)");
     PrintSkipped(3, "Number of used pixels < 3");
     PrintSkipped(4, "CorrXY==0");

trunk/MagicSoft/Mars/mimage/MHillasExt.cc

@@ -16 +16 @@
 !
 !
-!   Author(s): Thomas Bretz    12/2000 <mailto:tbretz@uni-sw.gwdg.de>
+!   Author(s): Thomas Bretz    12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
 !   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
+!   Copyright: MAGIC Software Development, 2000-2004
 !
 !
@@ -101 +101 @@
     fM3Trans =  0;
 
-    fMaxDist = -1;
+    fMaxDist =  0;
 }
 
@@ -170 +170 @@
         const Double_t dy = gpix.GetY() - hil.GetMeanY();      // [mm]
 
-        const Double_t dist = dx*dx+dy*dy;
-        if (dist>maxdist)
-            maxdist=dist;                                      // [mm^2]
-
         Double_t nphot = pix.GetNumPhotons();                  // [1]
 
         const Double_t dzx =  hil.GetCosDelta()*dx + hil.GetSinDelta()*dy; // [mm]
         const Double_t dzy = -hil.GetSinDelta()*dx + hil.GetCosDelta()*dy; // [mm]
+
+        const Double_t dist = dx*dx+dy*dy;
+        if (TMath::Abs(dist)>TMath::Abs(maxdist))
+            maxdist = dzx<0 ? -dist : dist;                    // [mm^2]
 
         m3x += nphot * dzx*dzx*dzx;                            // [mm^3]
@@ -210 +210 @@
     fM3Trans = m3y<0 ? -pow(-m3y, 1./3) : pow(m3y, 1./3);      // [mm]
 
-    fMaxDist = TMath::Sqrt(maxdist);                           // [mm]
+    const Double_t md = TMath::Sqrt(TMath::Abs(maxdist));
+    fMaxDist = maxdist<0 ? -md : md;                           // [mm]
 
     SetReadyToSave();
@@ -224 +225 @@
 void MHillasExt::Set(const TArrayF &arr)
 {
-    if (arr.GetSize() != 3)
+    if (arr.GetSize() != 4)
         return;
 
-    fAsym    = arr.At(0); // [mm] fDist minus dist: center of ellipse, highest pixel
-    fM3Long  = arr.At(1); // [mm] 3rd moment (e-weighted) along major axis
-    fM3Trans = arr.At(2); // [mm] 3rd moment (e-weighted) along minor axis
-}
+    fAsym    = arr.At(0);
+    fM3Long  = arr.At(1);
+    fM3Trans = arr.At(2);
+    fMaxDist = arr.At(3);
+}

trunk/MagicSoft/Mars/mimage/MNewImagePar.cc

@@ -17 +17 @@
 !
 !   Author(s): Wolfgang Wittek 03/2003 <mailto:wittek@mppmu.mpg.de>
+!   Author(s): Thomas Bretz            <mailto:tbretz@astro.uni-wuerzburg.de>
 !
 !   Copyright: MAGIC Software Development, 2000-2004
@@ -29 +30 @@
 // Storage Container for new image parameters
 //
-// fLeakage1 ratio: (photons in most outer ring of pixels) over fSize
-// fLeakage2 ratio: (photons in the 2 outer rings of pixels) over fSize
-// fNumSaturatedPixels: number of pixels in which at least one slice
-//                      of the low gain FADC was saturated.
+//    Float_t fLeakage1;             // (photons in most outer ring of pixels) over fSize
+//    Float_t fLeakage2;             // (photons in the 2 outer rings of pixels) over fSize
+//    Float_t fInnerLeakage1;        // (photons in most outer rings of inner pixels) over fInnerSize
+//    Float_t fInnerLeakage2;        // (photons in the 2 outer rings of inner pixels) over fInnerSize
+//    Float_t fInnerSize;            //
+//
+//    Float_t fConc;                 // [ratio] concentration ratio: sum of the two highest pixels / fSize
+//    Float_t fConc1;                // [ratio] concentration ratio: sum of the highest pixel / fSize
+//
+//    Float_t fUsedArea;             // Area of pixels which survived the image cleaning
+//    Float_t fCoreArea;             // Area of core pixels
+//    Short_t fNumUsedPixels;        // Number of pixels which survived the image cleaning
+//    Short_t fNumCorePixels;        // number of core pixels
+//    Short_t fNumHGSaturatedPixels; // number of pixels with saturating hi-gains
+//    Short_t fNumSaturatedPixels;   // number of pixels with saturating lo-gains
 //
 // Version 2:
@@ -46 +58 @@
 //  - added fUsedArea
 //  - added fCoreArea
-// 
-// 
+//
+//
 /////////////////////////////////////////////////////////////////////////////
 #include "MNewImagePar.h"
@@ -187 +199 @@
         }
 
-        // Compute Concetration 1 -2 
-
+        // Compute Concentration 1-2
         if (nphot>maxpix1)
         {
@@ -217 +228 @@
     *fLog << all;
     *fLog << "New Image Parameters (" << GetName() << ")" << endl;
-    *fLog << " - Leakage1       [1]   = " << fLeakage1      << endl;
-    *fLog << " - Leakage2       [1]   = " << fLeakage2      << endl;
-    *fLog << " - Conc           [1]   = " << fConc          << " (ratio)" << endl;
-    *fLog << " - Conc1          [1]   = " << fConc1         << " (ratio)" << endl;
-    *fLog << " - Used Pixels    [#]   = " << fNumUsedPixels << " Pixels" << endl;
-    *fLog << " - Core Pixels    [#]   = " << fNumCorePixels << " Pixels" << endl;
-    *fLog << " - Sat. Pixels (HG) [#]  = " << fNumHGSaturatedPixels << " Pixels" << endl;
-    *fLog << " - Sat. Pixels (LG) [#]  = " << fNumSaturatedPixels << " Pixels" << endl;
-}
+    *fLog << " - Leakage1       [1]   = " << fLeakage1             << endl;
+    *fLog << " - Leakage2       [1]   = " << fLeakage2             << endl;
+    *fLog << " - Conc           [1]   = " << fConc                 << " (ratio)" << endl;
+    *fLog << " - Conc1          [1]   = " << fConc1                << " (ratio)" << endl;
+    *fLog << " - Used Pixels    [#]   = " << fNumUsedPixels        << " Pixels" << endl;
+    *fLog << " - Core Pixels    [#]   = " << fNumCorePixels        << " Pixels" << endl;
+    *fLog << " - Used Area     [mm^2] = " << fUsedArea             << endl;
+    *fLog << " - Core Area     [mm^2] = " << fCoreArea             << endl;
+    *fLog << " - Sat.Pixels/HG  [#]   = " << fNumHGSaturatedPixels << " Pixels" << endl;
+    *fLog << " - Sat.Pixels/LG  [#]   = " << fNumSaturatedPixels   << " Pixels" << endl;
+}
+
+// -------------------------------------------------------------------------
+//
+// Print contents of MNewImagePar to *fLog, depending on the geometry in
+// units of deg.
+//
+void MNewImagePar::Print(const MGeomCam &geom) const
+{
+    *fLog << all;
+    *fLog << "New Image Parameters (" << GetName() << ")" << endl;
+    *fLog << " - Leakage1       [1]   = " << fLeakage1             << endl;
+    *fLog << " - Leakage2       [1]   = " << fLeakage2             << endl;
+    *fLog << " - Conc           [1]   = " << fConc                 << " (ratio)" << endl;
+    *fLog << " - Conc1          [1]   = " << fConc1                << " (ratio)" << endl;
+    *fLog << " - Used Pixels    [#]   = " << fNumUsedPixels        << " Pixels" << endl;
+    *fLog << " - Core Pixels    [#]   = " << fNumCorePixels        << " Pixels" << endl;
+    *fLog << " - Used Area    [deg^2] = " << fUsedArea*geom.GetConvMm2Deg()*geom.GetConvMm2Deg() << endl;
+    *fLog << " - Core Area    [deg^2] = " << fCoreArea*geom.GetConvMm2Deg()*geom.GetConvMm2Deg() << endl;
+    *fLog << " - Sat.Pixels/HG  [#]   = " << fNumHGSaturatedPixels << " Pixels" << endl;
+    *fLog << " - Sat.Pixels/LG  [#]   = " << fNumSaturatedPixels   << " Pixels" << endl;
+}

trunk/MagicSoft/Mars/mimage/MNewImagePar.h

@@ -43 +43 @@
     Short_t GetNumCorePixels() const { return fNumCorePixels; }
 
-    Float_t GetNumUsedArea() const { return fUsedArea; }
-    Float_t GetNumCoreArea() const { return fCoreArea; }
+    Float_t GetUsedArea() const { return fUsedArea; }
+    Float_t GetCoreArea() const { return fCoreArea; }
 
     Short_t GetNumSaturatedPixels() const { return fNumSaturatedPixels; }
@@ -50 +50 @@
 
     void Print(Option_t *opt=NULL) const;
+    void Print(const MGeomCam &geom) const;
 
     void Calc(const MGeomCam &geom, const MCerPhotEvt &evt,

trunk/MagicSoft/Mars/mmain/MEventDisplay.cc

@@ -82 +82 @@
 #include "MHillasExt.h"            // MHillasExt::Print(const MGeomCam&)
 #include "MHillasSrc.h"            // MHillasSrc::Print(const MGeomCam&)
+#include "MNewImagePar.h"          // MNewImagePar::Print(const MGeomCam&)
 #include "MHEvent.h"               // MHEvent
 #include "MHCamera.h"              // MHCamera
@@ -511 +512 @@
     // Print parameters
     //
-    ((MHillas*)   plist->FindObject("MHillas"))->Print(*geom);
-    ((MHillasExt*)plist->FindObject("MHillasExt"))->Print(*geom);
-    ((MHillasSrc*)plist->FindObject("MHillasSrc"))->Print(*geom);
-    plist->FindObject("MNewImagePar")->Print();
+    ((MHillas*)     plist->FindObject("MHillas"))->Print(*geom);
+    ((MHillasExt*)  plist->FindObject("MHillasExt"))->Print(*geom);
+    ((MHillasSrc*)  plist->FindObject("MHillasSrc"))->Print(*geom);
+    ((MNewImagePar*)plist->FindObject("MNewImagePar"))->Print(*geom);
 
     //