Changeset 4702 for trunk/MagicSoft/Mars/manalysis

Timestamp:

08/23/04 11:41:30 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/manalysis

Files:

• MCerPhotEvt.cc (modified) (8 diffs)
• MCerPhotEvt.h (modified) (4 diffs)
• MCerPhotPix.cc (modified) (2 diffs)
• MCerPhotPix.h (modified) (2 diffs)

trunk/MagicSoft/Mars/manalysis/MCerPhotEvt.cc

@@ -31 +31 @@
 //       risk!
 //
+// Class Version 3:
+// ----------------
+//  - added fNumIslands
+//            
 // Class Version 2:
 // ----------------
@@ -46 +50 @@
 #include <fstream>
 
+#include <TArrayD.h>
 #include <TCanvas.h>
 
@@ -52 +57 @@
 
 #include "MGeomCam.h"
+#include "MGeomPix.h"
 
 ClassImp(MCerPhotEvt);
@@ -96 +102 @@
 void MCerPhotEvt::Reset()
 {
-    fNumPixels =  0;
-    fMaxIndex  = -1;
+    fNumPixels  =  0;
+    fMaxIndex   = -1;
+    fNumIslands = -1;
     fLut.Set(0);
     // fPixels->Delete();
@@ -346 +353 @@
 void MCerPhotEvt::RemoveUnusedPixels()
 {
+    // Create iterator
     TIter Next(fPixels);
     MCerPhotPix *pix = NULL;
 
+    fMaxIndex = -1;
+
+    // Remove all unused pixels from list, calculate new fMaxIndex
     while ((pix=(MCerPhotPix*)Next()))
+    {
         if (!pix->IsPixelUsed())
             fPixels->Remove(pix);
-
+        else
+            fMaxIndex = TMath::Max(fMaxIndex, pix->GetPixId());
+    }
+
+    // Crompress array
     fPixels->Compress();
+
+    // Get new number of entries in array
     fNumPixels=fPixels->GetEntriesFast();
+
+    // Rebuild lookup table
+    RebuildLut();
 }
 
@@ -409 +430 @@
 // --------------------------------------------------------------------------
 //
+// This function recursively finds all pixels of one island and assigns
+// the number num as island number to the pixel.
+//
+//  1) Check whether a pixel with the index idx exists, is unused
+//     and has not yet a island number assigned.
+//  2) Assign the island number num to the pixel
+//  3) Loop over all its neighbors taken from the geometry geom. For all
+//     neighbors recursively call this function (CalcIsland)
+//  4) Sum the size of the pixel and all neighbors newly assigned
+//     (by CalcIsland) to this island
+//
+// Returns the sum of the pixel size.
+//
+Double_t MCerPhotEvt::CalcIsland(const MGeomCam &geom, Int_t idx, Int_t num)
+{
+    // Try to get the pixel information of a pixel with this index
+    MCerPhotPix *pix = GetPixById(idx);
+
+    // If a pixel with this index is not existing... do nothing.
+    if (!pix)
+        return 0;
+
+    // If an island number was already assigned to this pixel... do nothing.
+    if (pix->GetIdxIsland()>=0)
+        return 0;
+
+    // If the pixel is an unused pixel... do nothing.
+    if (!pix->IsPixelUsed())
+        return 0;
+
+    // Assign the new island number num to this used pixel
+    pix->SetIdxIsland(num);
+
+    // Get the geometry information (neighbors) of this pixel
+    const MGeomPix &gpix = geom[idx];
+
+    // Get the size of this pixel
+    Double_t size = pix->GetNumPhotons();
+
+    // Now do the same with all its neighbors and sum the
+    // sizes which they correspond to
+    const Int_t n = gpix.GetNumNeighbors();
+    for (int i=0; i<n; i++)
+        size += CalcIsland(geom, gpix.GetNeighbor(i), num);
+
+    // return size of this (sub)cluster
+    return size;
+}
+
+// --------------------------------------------------------------------------
+//
+// Each pixel which is maked as used is assigned an island number
+// (starting from 0). A pixel without an island number assigned
+// has island number -1.
+//
+// The index 0 corresponds to the island with the highest size (sum
+// of GetNumPhotons() in island). The size is decreasing with
+// increasing indices.
+//
+// The information about pixel neighbory is taken from the geometry
+// MGeomCam geom;
+//
+// You can access this island number of a pixel with a call to
+// MCerPhotPix->GetIdxIsland. The total number of islands available
+// can be accessed with MCerPhotEvt->GetNumIslands.
+//
+// CalcIslands returns the number of islands found. If an error occurs,
+// eg the geometry has less pixels than the highest index stored, -1 is
+// returned.
+//
+Int_t MCerPhotEvt::CalcIslands(const MGeomCam &geom)
+{
+    if (fMaxIndex<0 || fNumPixels<=0)
+    {
+        *fLog << err << "ERROR - MCerPhotEvt doesn't contain pixels!" << endl;
+        fNumIslands = 0;
+        return -1;
+    }
+
+    if ((UInt_t)fMaxIndex>=geom.GetNumPixels())
+    {
+        *fLog << err << "ERROR - MCerPhotEvt::CalcIslands: Size mismatch - geometry too small!" << endl;
+        return -1;
+    }
+
+    // Create a list to hold the sizes of the islands (The maximum
+    // number of islands possible is rougly fNumPixels/4)
+    TArrayD size(fNumPixels/3);
+
+    // Calculate Islands
+    Int_t n=0;
+
+    // We could loop over all indices which looks more straight
+    // forward but should be a lot slower (assuming zero supression)
+    MCerPhotPix *pix=0;
+
+    TIter Next(*this);
+    while ((pix=static_cast<MCerPhotPix*>(Next())))
+    {
+        // This 'if' is necessary (although it is done in GetIsland, too)
+        // because otherwise the counter (n) would be wrong.
+        // So only 'start' a new island for used pixels (selected by
+        // using the Iterator) which do not yet belong to another island.
+        if (pix->GetIdxIsland()<0)
+        {
+            Double_t sz = CalcIsland(geom, pix->GetPixId(), n);
+            size[n++] = sz;
+        }
+    }
+
+    // Create an array holding the indices
+    TArrayI idx(n);
+
+    // Sort the sizes descending
+    TMath::Sort(n, size.GetArray(), idx.GetArray(), kTRUE);
+
+    // Replace island numbers by size indices -- After this
+    // islands indices are sorted by the island size
+    Next.Reset();
+    while ((pix=static_cast<MCerPhotPix*>(Next())))
+    {
+        const Short_t i = pix->GetIdxIsland();
+
+        // Find new index
+        Short_t j;
+        for (j=0; j<n; j++)
+            if (idx[j]==i)
+                break;
+
+        pix->SetIdxIsland(j==n ? -1 : j);
+    }
+
+    // Now assign number of islands found
+    fNumIslands = n;
+
+    // return number of island
+    return fNumIslands;
+}
+
+// --------------------------------------------------------------------------
+//
 // Returns, depending on the type flag:
 //
@@ -416 +578 @@
 //  3: Number of Photons
 //  4: Error
+//  5: Island index
 //
 Bool_t MCerPhotEvt::GetPixelContent(Double_t &val, Int_t idx, const MGeomCam &cam, Int_t type) const
@@ -442 +605 @@
         val = pix->GetErrorPhot();
         break;
+    case 5:
+        val = pix->GetIdxIsland();
+        break;
     default:
         val = pix->GetNumPhotons()*ratio;

trunk/MagicSoft/Mars/manalysis/MCerPhotEvt.h

@@ -24 +24 @@
 private:
     UInt_t        fNumPixels;
+    Short_t       fNumIslands;
      Int_t        fMaxIndex;
     TArrayI       fLut;        // Lookup tabel to lookup pixel by index
     TClonesArray *fPixels;     //-> FIXME: Change TClonesArray away from a pointer?
+
+    void RebuildLut()
+    {
+        // Resize Lut
+        fLut.Set(fMaxIndex+1);
+
+        // Reset Lut
+        fLut.Reset(-1);
+
+        // Rebuild Lut
+        for (UInt_t i=0; i<GetNumPixels(); i++)
+        {
+            const MCerPhotPix &pix = (*this)[i];
+            fLut[pix.GetPixId()] = i;
+        }
+    }
+
+    Double_t CalcIsland(const MGeomCam &geom, Int_t idx, Int_t num);
 
 public:
@@ -32 +51 @@
     ~MCerPhotEvt() { delete fPixels; }
 
-    UInt_t GetNumPixels() const { return fNumPixels; }
-    //void   InitSize(UInt_t num) { fPixels->Expand(num); }
+    // Setter function to fill pixels
+    MCerPhotPix *AddPixel(Int_t idx, Float_t nph=0, Float_t er=0);
+    void FixSize();
 
-    MCerPhotPix *AddPixel(Int_t idx, Float_t nph=0, Float_t er=0);
-
-    void FixSize();
+    // Getter functions
+    UInt_t  GetNumPixels() const { return fNumPixels; }
+    Short_t GetNumIslands() const { return fNumIslands; };
 
     Bool_t  IsPixelExisting(Int_t id) const;
@@ -52 +72 @@
     Float_t GetErrorPhotMax(const MGeomCam *geom=NULL) const;
 
+    // Getter functions to access single pixels
     MCerPhotPix &operator[](int i)       { return *(MCerPhotPix*)(fPixels->UncheckedAt(i)); }
     MCerPhotPix &operator[](int i) const { return *(MCerPhotPix*)(fPixels->UncheckedAt(i)); }
 
+    MCerPhotPix *GetPixById(Int_t idx) const;
+
+    // Functions to change the contained data
     void Scale(Double_t f) { fPixels->ForEach(MCerPhotPix, Scale)(f); }
     void RemoveUnusedPixels();
+    Int_t CalcIslands(const MGeomCam &geom);
+    void Sort(Int_t upto = kMaxInt)
+    {
+        // Sort pixels by index
+        fPixels->Sort(upto);
+        RebuildLut();
+    } // For convinience: Sort pixels by index
 
-    MCerPhotPix *GetPixById(Int_t idx) const;
-
+    // class MParContainer
     void Reset();
 
+    // class TObject
     void Print(Option_t *opt=NULL) const;
     void Clear(Option_t *opt=NULL) { Reset(); }
@@ -69 +100 @@
     void DrawPixelContent(Int_t num) const;
 
+    // To build an iterator for this class
     operator TIterator*() const;
 

trunk/MagicSoft/Mars/manalysis/MCerPhotPix.cc

@@ -45 +45 @@
 //  - added fIsHGSaturated
 //
+// Version 5:
+// ----------
+//  - added fIdxIsland
+//
 ////////////////////////////////////////////////////////////////////////////
 #include "MCerPhotPix.h"
@@ -60 +64 @@
 //
 MCerPhotPix::MCerPhotPix(Int_t pix, Float_t phot, Float_t errphot) :
-    fPixId(pix), fRing(1), fIsCore(kFALSE), fPhot(phot), fErrPhot(errphot), 
+    fPixId(pix), fIsCore(kFALSE), fRing(1), fIdxIsland(-1),
+    fPhot(phot), fErrPhot(errphot),
     fIsSaturated(kFALSE), fIsHGSaturated(kFALSE)
 {
 } 
+
+// --------------------------------------------------------------------------
+//
+// From TObject:
+//  Compare abstract method. Must be overridden if a class wants to be able
+//  to compare itself with other objects. Must return -1 if this is smaller
+//  than obj, 0 if objects are equal and 1 if this is larger than obj.
+//
+// Here:
+//  Index numbers are compared
+//
+Int_t MCerPhotPix::Compare(const TObject *o) const
+{
+    const Int_t diff = fPixId - static_cast<const MCerPhotPix*>(o)->fPixId;
+    return diff==0 ? 0 : TMath::Sign(1, diff);
+}
 
 // --------------------------------------------------------------------------

trunk/MagicSoft/Mars/manalysis/MCerPhotPix.h

@@ -12 +12 @@
     Int_t    fPixId;     // the pixel Id
 
-   Short_t fRing;      // NT: number of analyzed rings around the core pixels, fRing>0 means: used, fRing= 0 means: unused, fRing= -1 means: unmapped (no possible to use in the calculation of the image parameters)
-    Bool_t   fIsCore;    // the pixel is a Core pixel -> kTRUE
+    Bool_t   fIsCore;     // the pixel is a Core pixel -> kTRUE
+    Short_t  fRing;       // NT: number of analyzed rings around the core pixels, fRing>0 means: used, fRing= 0 means: unused, fRing= -1 means: unmapped (no possible to use in the calculation of the image parameters)
+    Short_t  fIdxIsland;  // the pixel is a Core pixel -> kTRUE
 
     Float_t  fPhot;      // The number of Cerenkov photons
@@ -26 +27 @@
     MCerPhotPix(Int_t pix=-1, Float_t phot=0, Float_t errphot=0);
 
-    Int_t    GetPixId() const            { return fPixId;   }
-    Float_t  GetNumPhotons() const       { return fPhot;    }
-    Float_t  GetErrorPhot() const        { return fErrPhot; }
+    Int_t   GetPixId() const            { return fPixId;   }
+    Float_t GetNumPhotons() const       { return fPhot;    }
+    Float_t GetErrorPhot() const        { return fErrPhot; }
     
-    Bool_t   IsPixelUsed() const         { return fRing>0; }
-    Bool_t   IsPixelUnmapped() const     { return fRing==-1; }
-    void     SetPixelUnused()            { fRing=0; }
-    void     SetPixelUsed()              { fRing=1; }
-    void     SetPixelUnmapped()          { fRing=-1;}
+    Bool_t  IsPixelUsed() const         { return fRing>0; }
+    Bool_t  IsPixelUnmapped() const     { return fRing==-1; }
+    void    SetPixelUnused()            { fRing=0; }
+    void    SetPixelUsed()              { fRing=1; }
+    void    SetPixelUnmapped()          { fRing=-1;}
+    void    SetIdxIsland(Short_t num)   { fIdxIsland=num; }
+    Short_t GetIdxIsland() const        { return fIdxIsland; }
 
-    void     SetRing(UShort_t r)         { fRing = r;   }
-    Short_t  GetRing() const             { return fRing;}
+    void    SetRing(UShort_t r)         { fRing = r;   }
+    Short_t GetRing() const             { return fRing;}
 
-    void     SetPixelCore()              { fIsCore = kTRUE; }
-    Bool_t   IsPixelCore() const         { return fIsCore;  }
+    void    SetPixelCore()              { fIsCore = kTRUE; }
+    Bool_t  IsPixelCore() const         { return fIsCore;  }
 
-    void     SetNumPhotons(Float_t f)    { fPhot    = f; }
-    void     SetErrorPhot(Float_t f)     { fErrPhot = f; }
-    void     Set(Float_t np, Float_t ep) { fPhot = np; fErrPhot = ep; }
+    void    SetNumPhotons(Float_t f)    { fPhot    = f; }
+    void    SetErrorPhot(Float_t f)     { fErrPhot = f; }
+    void    Set(Float_t np, Float_t ep) { fPhot = np; fErrPhot = ep; }
 
-    void     SetPixelSaturated()         { fIsSaturated = kTRUE; }
-    Bool_t   IsPixelSaturated() const    { return fIsSaturated; }
+    void    SetPixelSaturated()         { fIsSaturated = kTRUE; }
+    Bool_t  IsPixelSaturated() const    { return fIsSaturated; }
 
-    void     SetPixelHGSaturated()         { fIsHGSaturated = kTRUE; }
-    Bool_t   IsPixelHGSaturated() const    { return fIsHGSaturated; }
+    void    SetPixelHGSaturated()       { fIsHGSaturated = kTRUE; }
+    Bool_t  IsPixelHGSaturated() const  { return fIsHGSaturated; }
 
-    void     AddNumPhotons(Float_t f)    { fPhot += f; }
+    void    AddNumPhotons(Float_t f)    { fPhot += f; }
 
-    void     Scale(Float_t f)            { fPhot/=f; }
+    void    Scale(Float_t f)            { fPhot/=f; }
 
-    void     Print(Option_t *opt = NULL) const;
+    void    Print(Option_t *opt = NULL) const;
+    Int_t   Compare(const TObject *obj) const;
+    Bool_t  IsSortable() const { return kTRUE; }
 
-    ClassDef(MCerPhotPix, 4)  // class containing information about the Cerenkov Photons in a pixel
+    ClassDef(MCerPhotPix, 5)  // class containing information about the Cerenkov Photons in a pixel
 };