Changeset 1489 for trunk/MagicSoft/Mars/mhist

Timestamp:

08/08/02 11:58:59 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/mhist

Files:

• MHHillas.cc (modified) (10 diffs)
• MHHillas.h (modified) (1 diff)
• MHHillasExt.cc (modified) (4 diffs)
• MHHillasSrc.cc (modified) (5 diffs)
• MHMatrix.cc (modified) (14 diffs)
• MHMatrix.h (modified) (3 diffs)

trunk/MagicSoft/Mars/mhist/MHHillas.cc

@@ -67 +67 @@
     fTitle = title ? title : "Source independent image parameters";
 
-    fLength = new TH1F("Length", "Length of Ellipse",               100,   0, 296.7);
-    fWidth  = new TH1F("Width",  "Width of Ellipse",                100,   0, 296.7);
-    fDistC  = new TH1F("DistC",  "Distance from center of camera",  100,   0, 445);
-    fDelta  = new TH1F("Delta",  "Angle (Main axis - x-axis)",      101, -90,  90);
+    fLength  = new TH1F("Length",  "Length of Ellipse",               100,   0, 296.7);
+    fWidth   = new TH1F("Width",   "Width of Ellipse",                100,   0, 296.7);
+    fDistC   = new TH1F("DistC",   "Distance from center of camera",  100,   0, 445);
+    fDelta   = new TH1F("Delta",   "Angle (Main axis - x-axis)",      101, -90,  90);
+    fUsedPix = new TH1F("UsedPix", "Number of used pixels",           150,   0, 150);
+    fCorePix = new TH1F("CorePix", "Number of core pixels",           150,   0, 150);
 
     fLength->SetDirectory(NULL);
@@ -76 +78 @@
     fDistC->SetDirectory(NULL);
     fDelta->SetDirectory(NULL);
+    fUsedPix->SetDirectory(NULL);
+    fCorePix->SetDirectory(NULL);
 
     fLength->SetXTitle("Length [mm]");
@@ -81 +85 @@
     fDistC->SetXTitle("Distance [mm]");
     fDelta->SetXTitle("Delta [\\circ]");
+    fUsedPix->SetXTitle("Number of Pixels");
+    fCorePix->SetXTitle("Number of Pixels");
 
     fLength->SetYTitle("Counts");
@@ -86 +92 @@
     fDistC->SetYTitle("Counts");
     fDelta->SetYTitle("Counts");
+    fUsedPix->SetYTitle("Counts");
+    fCorePix->SetYTitle("Counts");
 
     MBinning bins;
@@ -122 +130 @@
     delete fSize;
     delete fCenter;
+
+    delete fUsedPix;
+    delete fCorePix;
 }
 
@@ -150 +161 @@
     ApplyBinning(*plist, "Delta",  fDelta);
     ApplyBinning(*plist, "Size",   fSize);
+    ApplyBinning(*plist, "Pixels", fUsedPix);
+    ApplyBinning(*plist, "Pixels", fCorePix);
 
     const MBinning *bins = (MBinning*)plist->FindObject("BinningCamera");
@@ -235 +248 @@
     const Double_t scale = fUseMmScale ? 1 : fMm2Deg;
 
-    fLength->Fill(scale*h.GetLength());
-    fWidth ->Fill(scale*h.GetWidth());
-    fDistC ->Fill(scale*d);
-    fCenter->Fill(scale*h.GetMeanX(), scale*h.GetMeanY());
-    fDelta ->Fill(kRad2Deg*h.GetDelta());
-    fSize  ->Fill(h.GetSize());
+    fLength ->Fill(scale*h.GetLength());
+    fWidth  ->Fill(scale*h.GetWidth());
+    fDistC  ->Fill(scale*d);
+    fCenter ->Fill(scale*h.GetMeanX(), scale*h.GetMeanY());
+    fDelta  ->Fill(kRad2Deg*h.GetDelta());
+    fSize   ->Fill(h.GetSize());
+    fUsedPix->Fill(h.GetNumUsedPixels());
+    fCorePix->Fill(h.GetNumCorePixels());
 
     return kTRUE;
@@ -269 +284 @@
 TObject *MHHillas::DrawClone(Option_t *opt) const
 {
-    TCanvas *c = MakeDefCanvas("Hillas", fTitle, 720, 810);
+    TCanvas *c = MakeDefCanvas(this, 720, 810);
     c->Divide(2,3);
 
@@ -275 +290 @@
 
     c->cd(1);
-    fLength->DrawCopy();
+    fWidth->DrawCopy();
+    fLength->SetLineColor(kBlue);
+    fLength->DrawCopy("same");
 
     c->cd(2);
-    fWidth->DrawCopy();
-
-    c->cd(3);
     gPad->SetLogx();
     fSize->DrawCopy();
+
+    c->cd(3);
+    fCorePix->SetLineColor(kRed);
+    fCorePix->DrawCopy();
+    fUsedPix->SetLineColor(kGreen);
+    fUsedPix->DrawCopy("same");
 
     c->cd(4);
@@ -309 +329 @@
 {
     if (!gPad)
-        MakeDefCanvas("Hillas", fTitle, 720, 810);
+        MakeDefCanvas(this, 720, 810);
 
     gPad->Divide(2,3);
 
     gPad->cd(1);
-    fLength->Draw();
+    fWidth->Draw();
+    fLength->SetLineColor(kBlue);
+    fLength->Draw("same");
 
     gPad->cd(2);
-    fWidth->Draw();
-
-    gPad->cd(3);
     gPad->SetLogx();
     fSize->Draw();
+
+    gPad->cd(3);
+    fCorePix->SetLineColor(kRed);
+    fCorePix->Draw();
+    fUsedPix->SetLineColor(kGreen);
+    fUsedPix->Draw("same");
 
     gPad->cd(4);

trunk/MagicSoft/Mars/mhist/MHHillas.h

@@ -22 +22 @@
     TH1F *fSize;   //->
     TH2F *fCenter; //->
+
+    TH1F *fUsedPix; //->
+    TH1F *fCorePix; //->
 
     void SetColors() const;

trunk/MagicSoft/Mars/mhist/MHHillasExt.cc

@@ -62 +62 @@
     //
     fName  = name  ? name  : "MHHillasExt";
-    fTitle = title ? title : "Container for Hillas (ext) histograms";
+    fTitle = title ? title : "Container for extended Hillas histograms";
 
     //
@@ -69 +69 @@
     // connect all the histogram with the container fHist
     //
+    fHConc1.SetLineColor(kBlue);
+    fHConc.SetFillStyle(0);
+
     fHConc.SetDirectory(NULL);
     fHConc1.SetDirectory(NULL);
@@ -220 +223 @@
 TObject *MHHillasExt::DrawClone(Option_t *opt) const
 {
-    TCanvas &c = *MakeDefCanvas("Hillas", "Histograms of Hillas Parameters",
-                                720, 810);
-    c.Divide(2, 3);
+    TCanvas &c = *MakeDefCanvas(this, 720, 540);
+    c.Divide(2, 2);
 
     gROOT->SetSelectedPad(NULL);
 
-    //
-    // This is necessary to get the expected bahviour of DrawClone
-    //
+
     c.cd(1);
-    ((TH1F&)fHConc).DrawCopy();
+    ((TH1F)fHConc1).DrawCopy();
+    ((TH1F)fHConc).DrawCopy("same");
 
     c.cd(2);
-    ((TH1F&)fHConc1).DrawCopy();
-
-    c.cd(5);
-    ((TH1F&)fHAsym).DrawCopy();
+    ((TH1F)fHAsym).DrawCopy();
 
     c.cd(3);
-    ((TH1F&)fHM3Long).DrawCopy();
+    ((TH1F)fHM3Long).DrawCopy();
 
     c.cd(4);
-    ((TH1F&)fHM3Trans).DrawCopy();
+    ((TH1F)fHM3Trans).DrawCopy();
 
     c.Modified();
@@ -261 +259 @@
 {
     if (!gPad)
-        MakeDefCanvas("Hillas", "Histograms of Hillas Parameters",
-                      720, 810);
-
-    gPad->Divide(2, 3);
+        MakeDefCanvas(this, 720, 540);
+
+    gPad->Divide(2, 2);
 
     gPad->cd(1);
-    fHConc.DrawCopy();
+    fHConc1.SetLineColor(kBlue);
+    fHConc1.Draw();
+    fHConc.Draw("same");
 
     gPad->cd(2);
-    fHConc1.DrawCopy();
-
-    gPad->cd(5);
-    fHAsym.DrawCopy();
+    fHAsym.Draw();
 
     gPad->cd(3);
-    fHM3Long.DrawCopy();
+    fHM3Long.Draw();
 
     gPad->cd(4);
-    fHM3Trans.DrawCopy();
+    fHM3Trans.Draw();
 
     gPad->Modified();

trunk/MagicSoft/Mars/mhist/MHHillasSrc.cc

@@ -68 +68 @@
     // connect all the histogram with the container fHist
     //
-    fAlpha    = new TH1F("Alpha",    "Alpha of Ellipse",             90,    0,  90);
+    fAlpha    = new TH1F("Alpha",    "Alpha of Ellipse",            181,  -90,  90);
     fDist     = new TH1F("Dist",     "Dist of Ellipse",             100,    0, 445);
     fHeadTail = new TH1F("HeadTail", "HeadTail of Ellipse",         101, -445, 445);
@@ -81 +81 @@
     fDist->SetXTitle("Dist [mm]");
     fHeadTail->SetXTitle("Head-Tail [mm]");
-    fCosDA->SetXTitle("cos(\\delta,\\alpha) [mm]");
+    fCosDA->SetXTitle("cos(\\delta,\\alpha)");
 
     fAlpha->SetYTitle("Counts");
@@ -138 +138 @@
     const MHillasSrc &h = *(MHillasSrc*)par;
 
-    fAlpha   ->Fill(fabs(h.GetAlpha()));
+    fAlpha   ->Fill(h.GetAlpha());
     fDist    ->Fill(fUseMmScale ? h.GetDist() : fMm2Deg*h.GetDist());
     fHeadTail->Fill(fUseMmScale ? h.GetHeadTail() : fMm2Deg*h.GetHeadTail());
@@ -211 +211 @@
 TObject *MHHillasSrc::DrawClone(Option_t *opt) const
 {
-    TCanvas *c = MakeDefCanvas("HillasSrc", "Histograms of Source dependant Parameters",
-                               700, 500);
+    TCanvas *c = MakeDefCanvas(this, 700, 500);
     c->Divide(2, 2);
 
@@ -250 +249 @@
 {
     if (!gPad)
-        MakeDefCanvas("HillasSrc", "Histograms of Src dependant Parameters",
-                      700, 500);
+        MakeDefCanvas(this, 700, 500);
 
     // FIXME: Display Source position

trunk/MagicSoft/Mars/mhist/MHMatrix.cc

@@ -45 +45 @@
 #include "MHMatrix.h"
 
+#include <fstream.h>
+
 #include <TList.h>
 #include <TArrayD.h>
@@ -54 +56 @@
 #include "MParList.h"
 
-#include "MDataChain.h"
+#include "MData.h"
+#include "MDataArray.h"
 
 ClassImp(MHMatrix);
 
-// --------------------------------------------------------------------------
-//
-// Default Constructor
+static const TString gsDefName  = "MHMatrix";
+static const TString gsDefTitle = "Multidimensional Matrix";
+
+// --------------------------------------------------------------------------
+//
+//  Default Constructor
 //
 MHMatrix::MHMatrix(const char *name, const char *title)
-    : fNumRow(0)
-{
-    fName  = name  ? name  : "MHMatrix";
-    fTitle = title ? title : "Multidimensional Matrix";
-
-    fData = new TList;
-    fData->SetOwner();
-
-    fRules = new TList;
-    fRules->SetOwner();
-}
-
-// --------------------------------------------------------------------------
-//
-// Destructor
+    : fNumRow(0), fData(NULL)
+{
+    fName  = name  ? name  : gsDefName.Data();
+    fTitle = title ? title : gsDefTitle.Data();
+}
+
+// --------------------------------------------------------------------------
+//
+//  Constructor. Initializes the columns of the matrix with the entries
+//  from a MDataArray
+//
+MHMatrix::MHMatrix(MDataArray *mat, const char *name, const char *title)
+    : fNumRow(0), fData(mat)
+{
+    fName  = name  ? name  : gsDefName.Data();
+    fTitle = title ? title : gsDefTitle.Data();
+}
+
+// --------------------------------------------------------------------------
+//
+//  Destructor. Does not deleted a user given MDataArray, except IsOwner
+//  was called.
 //
 MHMatrix::~MHMatrix()
 {
-    delete fData;
-    delete fRules;
+    if (TestBit(kIsOwner) && fData)
+        delete fData;
 }
 
@@ -99 +112 @@
     }
 
-    MDataChain &chain = *new MDataChain(rule);
-    if (!chain.IsValid())
-    {
-        *fLog << err << "Error - Rule cannot be translated... ignored." << endl;
-        delete &chain;
+    if (!fData)
+    {
+        fData = new MDataArray;
+        SetBit(kIsOwner);
+    }
+
+    fData->AddEntry(rule);
+}
+
+void MHMatrix::AddColumns(MDataArray *matrix)
+{
+    if (fM.IsValid())
+    {
+        *fLog << warn << "Warning - matrix is already in use. Can't add new columns... skipped." << endl;
         return;
     }
-    fData->Add(&chain);
-
-    TNamed *name = new TNamed(rule, rule); // Fimxe, in 3.02/07 the title can't be "", why?
-    fRules->Add(name);
-}
-
-void MHMatrix::AddColumns(const MHMatrix *matrix)
-{
-    TIter Next(matrix->fRules);
-    TObject *rule=NULL;
-    while ((rule=Next()))
-        AddColumn(rule->GetName());
+
+    if (fData)
+        *fLog << warn << "Warning - columns already added... replacing." << endl;
+
+    if (fData && TestBit(kIsOwner))
+    {
+        delete fData;
+        ResetBit(kIsOwner);
+    }
+
+    fData = matrix;
 }
 
@@ -127 +148 @@
 Bool_t MHMatrix::SetupFill(const MParList *plist)
 {
-    if (fData->GetSize()==0)
-    {
-        *fLog << err << "Error - No Column specified... aborting." << endl;
+    if (!fData)
+    {
+        *fLog << err << "Error - No Columns initialized... aborting." << endl;
         return kFALSE;
     }
 
-    TIter Next(fData);
-    MData *data = NULL;
-    while ((data=(MData*)Next()))
-        if (!data->PreProcess(plist))
-            return kFALSE;
-
-    return kTRUE;
+    return fData->PreProcess(plist);
 }
 
@@ -155 +170 @@
     if (!fM.IsValid())
     {
-        fM.ResizeTo(1, fData->GetSize());
+        fM.ResizeTo(1, fData->GetNumEntries());
         return;
     }
@@ -161 +176 @@
     TMatrix m(fM);
 
-    fM.ResizeTo(fM.GetNrows()*2, fData->GetSize());
+    fM.ResizeTo(fM.GetNrows()*2, fData->GetNumEntries());
 
     for (int x=0; x<m.GetNcols(); x++)
@@ -176 +191 @@
     AddRow();
 
-    int col=0;
-    TIter Next(fData);
-    MData *data = NULL;
-    while ((data=(MData*)Next()))
-        fM(fNumRow-1, col++) = data->GetValue();
+    for (int col=0; col<fData->GetNumEntries(); col++)
+        fM(fNumRow-1, col) = (*fData)(col);
 
     return kTRUE;
@@ -195 +207 @@
     // It's not a fatal error so we don't need to stop PostProcessing...
     //
-    if (fData->GetSize()<1 || fNumRow<1)
+    if (fData->GetNumEntries()==0 || fNumRow<1)
         return kTRUE;
 
     TMatrix m(fM);
 
-    fM.ResizeTo(fNumRow, fData->GetSize());
+    fM.ResizeTo(fNumRow, fData->GetNumEntries());
 
     for (int x=0; x<fM.GetNcols(); x++)
@@ -287 +299 @@
 void MHMatrix::Print(Option_t *) const
 {
-    int n=0;
-
-    TIter Next(fData->GetSize() ? fData : fRules);
-    MData *data = NULL;
-    while ((data=(MData*)Next()))
-    {
-        *fLog << all << " Column " << setw(3) << n++ << ": " << flush;
-        data->Print();
-        *fLog << endl;
-    }
+    if (fData)
+        fData->Print();
+    else
+        *fLog << all << "Sorry, no column information available." << endl;
 
     fM.Print();
@@ -336 +342 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Calculated the distance of vector evt from the reference sample
+// represented by the covariance metrix m.
+//  - If n<0 the kernel method is applied and
+//    sum(epx(-d)/n is returned.
+//  - For n>=0 the n nearest neighbors are summed and
+//    sqrt(sum(d))/n is returned.
+//  - if n==0 all distances are summed
+//
 Double_t MHMatrix::CalcDist(const TMatrix &m, const TVector &evt, Int_t num) const
 {
+    if (num==0) // may later be used for another method
+    {
+        TVector d = evt;
+        d *= m;
+        return sqrt(d*evt); 
+    }
+
     const Int_t rows = fM.GetNrows();
     const Int_t cols = fM.GetNcols();
@@ -358 +381 @@
 
         dists[i] = d2*d; // square of distance
+
+        //
+        // This corrects for numerical uncertanties in cases of very
+        // small distances...
+        //
+        if (dists[i]<0)
+            dists[i]=0;
     }
 
@@ -363 +393 @@
     TMath::Sort(dists.GetSize(), dists.GetArray(), idx.GetArray(), kFALSE);
 
-    const Int_t n = num<rows ? num : rows;
-
-    Double_t res = 0;
-    for (int i=0; i<n; i++)
-        res += dists[idx[i]];
-
-    return sqrt(res)/n;
-}
-
+    const Int_t n = abs(num)<rows ? abs(num) : rows;
+
+    if (num<0)
+    {
+        //
+        // Kernel function sum
+        //
+        const Double_t h = 1;
+
+        Double_t res = 0;
+        for (int i=0; i<n; i++)
+            res += exp(-dists[idx[i]]/h);
+
+        return log(res/n);
+    }
+    else
+    {
+        //
+        // Nearest Neighbor sum
+        //
+        Double_t res = 0;
+        for (int i=0; i<n; i++)
+            res += dists[idx[i]];
+
+        return sqrt(res/n);
+    }
+}
+
+// --------------------------------------------------------------------------
+//
+// Calls calc dist. In the case of the first call the covariance matrix
+// fM2 is calculated.
+//  - If n<0 it is divided by (nrows-1)/h while h is the kernel factor.
+//
 Double_t MHMatrix::CalcDist(const TVector &evt, Int_t num)
 {
@@ -379 +434 @@
         fM2.ResizeTo(m);
         fM2 = m;
+        fM2 *= fM.GetNrows()-1;
         delete &m;
     }
@@ -391 +447 @@
     fM.ResizeTo(m);
     fM = m;
-
-    TIter Next(fRules);
-    TObject *obj = NULL;
-    while ((obj=Next()))
-        fData->Add(new MDataChain(obj->GetName()));
-}
+}
+
+void MHMatrix::StreamPrimitive(ofstream &out) const
+{
+    Bool_t data = fData && !TestBit(kIsOwner);
+
+    if (data)
+    {
+        fData->SavePrimitive(out);
+        out << endl;
+    }
+
+    out << "   MHMatrix " << GetUniqueName();
+
+    if (data || fName!=gsDefName || fTitle!=gsDefTitle)
+    {
+        out << "(";
+        if (data)
+            out << "&" << fData->GetUniqueName();
+        if (fName!=gsDefName || fTitle!=gsDefTitle)
+        {
+            if (data)
+                out << ", ";
+            out << "\"" << fName << "\"";
+            if (fTitle!=gsDefTitle)
+                out << ", \"" << fTitle << "\"";
+        }
+    }
+    out << ");" << endl;
+
+    if (fData && TestBit(kIsOwner))
+        for (int i=0; i<fData->GetNumEntries(); i++)
+            out << "   " << GetUniqueName() << ".AddColumn(\"" << (*fData)[i].GetRule() << "\");" << endl;
+}

trunk/MagicSoft/Mars/mhist/MHMatrix.h

@@ -9 +9 @@
 #endif
 
-class MDataChain;
+class MDataArray;
 
 class MHMatrix : public MH
 {
 protected:
-    Int_t   fNumRow; //! Number of dimensions of histogram
-    TMatrix fM;      // Matrix to be filled
+    Int_t   fNumRow;    //! Number of dimensions of histogram
+    TMatrix fM;         // Matrix to be filled
 
-    TMatrix fM2;     //!
+    TMatrix fM2;        //!
 
-    TList  *fData;   //! List of data members (columns)
-    TList  *fRules;  // List of data members as text for storage
+    MDataArray *fData;  // List of data members (columns)
+
+    enum {
+        kIsOwner = BIT(14)
+    };
 
     void AddRow();
@@ -28 +31 @@
     Bool_t Finalize();
 
+    void StreamPrimitive(ofstream &out) const;
+
 public:
     MHMatrix(const char *name=NULL, const char *title=NULL);
+    MHMatrix(MDataArray *mat, const char *name=NULL, const char *title=NULL);
     ~MHMatrix();
 
     void AddColumn(const char *name);
-    void AddColumns(const MHMatrix *matrix);
+    void AddColumns(MDataArray *mat);
 
-    //    TMatrix &GetM() { return fM; }
+    MDataArray *GetColumns() { return fData; }
+
     const TMatrix &GetM() const { return fM; }
-    const TList *GetRules() const { return fRules; }
 
     //void Draw(Option_t *opt=NULL);
@@ -49 +55 @@
     Double_t CalcDist(const TVector &v, Int_t num = 25);
 
+    void SetIOwner(Bool_t b=kTRUE) { b ? SetBit(kIsOwner) : ResetBit(kIsOwner); }
+
     void Reassign();