Changeset 9349

Timestamp:

02/19/09 16:27:33 (16 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars

Files:

• Changelog (modified) (1 diff)
• NEWS (modified) (1 diff)
• msim/MPhotonEvent.cc (modified) (8 diffs)
• msim/MPhotonEvent.h (modified) (3 diffs)
• msimcamera/MSimAPD.cc (modified) (2 diffs)
• msimreflector/MSimReflector.cc (modified) (7 diffs)

trunk/MagicSoft/Mars/Changelog

@@ -35 +35 @@
    * msim/MSimAbsorption.cc:
      - don't check wavelength range when theta should be used
+
+   * msimcamera/MSimAPD.cc:
+     - check for uninitialized indices
+
+   * msim/MPhotonEvent.[h,cc]:
+     - moved the code for MyClonesArray to the source file
+     - improved a lot the reading speed by reading larger blocks
+       of data from the file at once
+     - improved memory handling. This ensures that even the largest
+       events don't fill the memory forever and the allocated memory
+       is free'd again after some time
+
+   * msimreflector/MSimReflector.cc:
+     - Use the new Resize function of MPhotonEvent to make sure
+       that the memory is not allocated forever.
 
 

trunk/MagicSoft/Mars/NEWS

@@ -69 +69 @@
      layout
 
+   * Improved reading speed for corsika files
+
+   * Improved memory handling (if a large ampount of memory was needed
+     for a single event all further events were stored in the same 
+     memory and it was never freed, so the program took this memory
+     until the end)
+ 
  ;star
 

trunk/MagicSoft/Mars/msim/MPhotonEvent.cc

@@ -126 +126 @@
 using namespace std;
 
+// ==========================================================================
+
+class MyClonesArray : public TClonesArray
+{
+public:
+    TObject **FirstRef() { return fCont; }
+
+    // --------------------------------------------------------------------------
+    //
+    // This is an extremly optimized version of ExpandCreateFast. It only resets
+    // the marker for the last element (fLast) to n-1 and doen't change anything
+    // else. This implicitly assumes that the stored objects don't allocate
+    // memory. It does not necessarily mean that the slots after fLast
+    // are empty (set to 0). This is what is assumed in the TClonesArray.
+    // We also don't call Changed() because it would reset Sorted. If the
+    // array was sorted before it is also sorted now. You MUST make sure
+    // that you only set n in a range for which valid entries have been
+    // created before (e.g. by ExpandCreateFast).
+    //
+    void FastShrink(Int_t n)
+    {
+        fLast = n - 1;
+    }
+
+    // --------------------------------------------------------------------------
+    //
+    // This is a optimized (faster) version of Delete which deletes consequtive
+    // entries from index idx1 to idx2 (both included) and calls their
+    // destructor. Note that there is no range checking done!
+    //
+    void FastRemove(Int_t idx1, Int_t idx2)
+    {
+        // Remove object at index idx.
+
+        //if (!BoundsOk("RemoveAt", idx1)) return 0;
+        //if (!BoundsOk("RemoveAt", idx2)) return 0;
+
+        Long_t dtoronly = TObject::GetDtorOnly();
+
+        idx1 -= fLowerBound;
+        idx2 -= fLowerBound;
+
+        for (TObject **obj=fCont+idx1; obj<=fCont+idx2; obj++)
+        {
+            if (!*obj)
+                continue;
+
+            if ((*obj)->TestBit(kNotDeleted)) {
+                // Tell custom operator delete() not to delete space when
+                // object fCont[i] is deleted. Only destructors are called
+                // for this object.
+                TObject::SetDtorOnly(*obj);
+                delete *obj;
+            }
+
+            *obj = 0;
+            // recalculate array size
+        }
+        TObject::SetDtorOnly((void*)dtoronly);
+
+        if (idx1<=fLast && fLast<=idx2)
+        {
+            do {
+                fLast--;
+            } while (fLast >= 0 && fCont[fLast] == 0);
+        }
+
+        Changed();
+    }
+
+
+    //void SetSorted() { fSorted = kTRUE; }
+
+    // --------------------------------------------------------------------------
+    //
+    // This is an optimized version of Sort which doesn't check the
+    // IsSortable flag before. It only sorts the entries from 0
+    // to GetEntriesFast().
+    //
+    void UncheckedSort()
+    {
+        if (fSorted)
+            return;
+
+        const Int_t nentries = GetEntriesFast();
+        if (nentries <= 0)
+            return;
+
+        QSort(GetObjectRef(First()), fKeep->GetObjectRef(fKeep->First()),
+              0, TMath::Min(nentries, kMaxInt-fLowerBound));
+
+        fSorted = kTRUE;
+    }
+};
+
+// ==========================================================================
+
 // --------------------------------------------------------------------------
 //
@@ -140 +237 @@
 }
 
-/*
-const char *MPhotonEvent::GetClassName() const
-{
-    return static_cast<TObject*>(fData.GetClass())->GetName();
-}
-*/
+// --------------------------------------------------------------------------
+//
+// This is an extremly optimized version of ExpandCreateFast. It only resets
+// the marker for the last element (fLast) to n-1 and doen't change anything
+// else. This has the advantage that the allocated memory is kept but only
+// valid entries are written to a file.
+//
+// If the array was sorted before it is also sorted now. You MUST make sure
+// that you only set n in a range for which valid entries have been
+// created before (e.g. by ExpandCreateFast).
+//
+Int_t MPhotonEvent::Shrink(Int_t n)
+{
+    /*
+    // The number of object written by the streamer is defined by
+    // GetEntriesFast(), i.e. this number must be shrinked to
+    // the real array size. We use RemoveAt instead of ExpandCreate
+    // because RemoveAt doesn't free memory. Thus in the next
+    // iteration it can be reused and doesn't need to be reallocated.
+    // Do not change this. It is optimized for execution speed
+    //        for (int i=fData.GetSize()-1; i>=n; i--)
+    //          fData.RemoveAt(i);
+    const Bool_t sorted = fData.IsSorted();
+
+    MyClonesArray &loc = static_cast<MyClonesArray&>(fData);
+
+    loc.FastRemove(n, fData.GetSize()-1);
+
+    // If it was sorted before it is also sorted now.
+    if (sorted)
+        loc.SetSorted();
+    */
+
+    // fData.ExpandCreateFast(n);  // Just set fLast = n -1
+
+    // Just set fLast = n -1
+    static_cast<MyClonesArray&>(fData).FastShrink(n);
+    return fData.GetEntriesFast();
+}
+
+// --------------------------------------------------------------------------
+//
+// The resized the array. If it has to be increased in size it is done
+// with ExpandCreateFast. If it should be shrinked it is done with
+// ExpandCreateFast if n>fData.GetSize()/100 or n==0. This keeps the allocated
+// memory and just sets the marker for the last element in the array (fLast).
+//
+// If the allocated memory grew to huge we reset the allocated memory
+// by calling ExpandCreate(n) (frees the allocated storage for the
+// objects) and Expand(n) (frees the allocated memory for the list
+// of pointers to the objects)
+//
+// 100 hundred is an arbitrary number taking into account that todays
+// computers have a lot of memory and we don't want to free and allocate
+// new memory too often.
+//
+// In priciple there might be more clever methods to handle the memory.
+//
+void MPhotonEvent::Resize(Int_t n)
+{
+    if (n==0 || n*100>fData.GetSize())
+        fData.ExpandCreateFast(n);  // Just set fLast = n -1
+    else
+    {
+        fData.ExpandCreate(n);      // Free memory of allocated MPhotonData
+        fData.Expand(n);            // Free memory of allocated pointers
+    }
+}
 
 // --------------------------------------------------------------------------
@@ -159 +318 @@
     // Do not modify this. It is optimized for execution
     // speed and flexibility!
-    TObject *o = n<fData.GetSize() && fData.UncheckedAt(n) ? fData.UncheckedAt(n) : fData.New(n);
-
-    // Now we read a single cherenkov bunch
+    TObject *o = 0;
+    if (n<fData.GetSize() && fData.UncheckedAt(n))
+    {
+        o=fData.UncheckedAt(n);
+        static_cast<MyClonesArray&>(fData).FastShrink(n+1);
+    }
+    else
+    {
+        o=fData.New(n);
+    }
     return static_cast<MPhotonData&>(*o);
 }
@@ -175 +341 @@
 {
     return Add(GetNumPhotons());
+}
+
+void MPhotonEvent::Sort(Bool_t force)
+{
+    if (force)
+        fData.UnSort();
+
+    static_cast<MyClonesArray&>(fData).UncheckedSort(); /*Sort(GetEntriesFast())*/
 }
 
@@ -205 +379 @@
 MPhotonData *MPhotonEvent::GetFirst() const
 {
-    return static_cast<MPhotonData*>(fData.First());
+    return fData.GetEntriesFast()==0 ? 0 : static_cast<MPhotonData*>(fData.First());
 }
 
@@ -214 +388 @@
 MPhotonData *MPhotonEvent::GetLast() const
 {
-    return static_cast<MPhotonData*>(fData.Last());
+    return fData.GetEntriesFast()==0 ? 0 : static_cast<MPhotonData*>(fData.Last());
 }
 
@@ -282 +456 @@
     Int_t n = 0;
 
+    // --- old I/O ---
+    // Read only + Reflector (no absorption)
+    // Muons:   1.06GB/115s =  9.2MB/s (100kEvs)
+    // Gammas:  1.57GB/275s =  5.7MB/s (  1kEvs)
+
+    // Read only:
+    // Gammas:  1.57GB/158s =  9.9MB/s (  1kEvs)
+    // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)
+
+    // --- new I/O ---
+    // Read only (don't allocate storage space):
+    // Gammas:  1.57GB/143s = 11.0MB/s (  1kEvs)
+    // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)
+
+    // Read only in blocks (with storage allocation):
+    // Gammas:  1.57GB/28s  =  56MB/s (  1kEvs)
+    // Muons:   1.06GB/5.2s = 204MB/s (100kEvs)
+
+    // Read only in blocks (without storage allocation):
+    // similar to just copy
+
+    // Copy with cp
+    // 1.57GB/ 5s   CPU
+    // 1.57GB/28s   REAL
+    // 1.06GB/ 3s   CPU
+    // 1.06GB/22s   REAL
+
+    while (1)
+    {
+        // Stprage for one block
+        char c[273*4];
+
+        // Read the first four byte to check whether the next block
+        // doen't belong to the event anymore
+        fin.read(c, 4);
+        if (!fin)
+            return kFALSE;
+
+        // Check if the event is finished
+        if (!memcmp(c, "EVTE", 4))
+            break;
+
+        // Now read the rest of the data
+        fin.read(c+4, 272*4);
+
+        Float_t *ptr = reinterpret_cast<Float_t*>(c);
+        Float_t *end = ptr + 273;
+
+        // Loop over all sub-blocks (photons) in the block and if
+        // they contain valid data add them to the array
+        while (ptr<end)
+        {
+            // Get/Add the n-th entry from the array and
+            // fill it with the current 7 floats
+            const Int_t rc = Add(n).FillCorsika(ptr);
+            ptr += 7;
+
+            switch (rc)
+            {
+            case kCONTINUE:  continue;        // No data in this bunch... skip it.
+            case kERROR:     return kERROR;   // Error occured
+            //case kFALSE:     return kFALSE;   // End of stream
+            }
+
+            // This is a photon we would like to keep later.
+            // Increase the counter by one
+            n++;
+        }
+    }
+/*
     while (1)
     {
@@ -326 +570 @@
         n++;
     }
-
-    Shrink(n);
+  */
+
+    Resize(n);
     fData.UnSort();
 

trunk/MagicSoft/Mars/msim/MPhotonEvent.h

@@ -12 +12 @@
 #include <iosfwd>
 
+using namespace std;
+
 class MPhotonData;
 class MCorsikaRunHeader;
-
-class MyClonesArray : public TClonesArray
-{
-public:
-    TObject **FirstRef() { return fCont; }
-
-    void FastRemove(Int_t idx1, Int_t idx2)
-    {
-        // Remove object at index idx.
-
-        //if (!BoundsOk("RemoveAt", idx1)) return 0;
-        //if (!BoundsOk("RemoveAt", idx2)) return 0;
-
-        Long_t dtoronly = TObject::GetDtorOnly();
-
-        idx1 -= fLowerBound;
-        idx2 -= fLowerBound;
-
-        for (TObject **obj=fCont+idx1; obj<=fCont+idx2; obj++)
-        {
-            if (!*obj)
-                continue;
-
-            if ((*obj)->TestBit(kNotDeleted)) {
-                // Tell custom operator delete() not to delete space when
-                // object fCont[i] is deleted. Only destructors are called
-                // for this object.
-                TObject::SetDtorOnly(*obj);
-                delete *obj;
-            }
-
-            *obj = 0;
-            // recalculate array size
-        }
-        TObject::SetDtorOnly((void*)dtoronly);
-
-        if (idx1<=fLast && fLast<=idx2)
-        {
-            do {
-                fLast--;
-            } while (fLast >= 0 && fCont[fLast] == 0);
-        }
-
-        Changed();
-    }
-
-    void SetSorted() { fSorted = kTRUE; }
-};
 
 class MPhotonEvent : public MParContainer
@@ -71 +25 @@
     MPhotonEvent(const char *name=NULL, const char *title=NULL);
 
-    // TObjArray, FIXME: This could be improved if checking for IsSortable is omitted
-    void Sort(Bool_t force=kFALSE) { if (force) fData.UnSort(); fData.Sort(); }
+    void Sort(Bool_t force=kFALSE);
     Bool_t IsSorted() const { return fData.IsSorted(); }
-
 
     // Getter/Setter
@@ -96 +48 @@
     const MPhotonData &operator[](UInt_t idx) const;
 
-    Int_t Shrink(Int_t n)
-    {
-        // The number of object written by the streamer is defined by
-        // GetEntriesFast(), i.e. this number must be shrinked to
-        // the real array size. We use RemoveAt instead of ExpandCreate
-        // because RemoveAt doesn't free memory. Thus in the next
-        // iteration it can be reused and doesn't need to be reallocated.
-        // Do not change this. It is optimized for execution speed
-        //        for (int i=fData.GetSize()-1; i>=n; i--)
-        //          fData.RemoveAt(i);
-        const Bool_t sorted = fData.IsSorted();
-
-        MyClonesArray &loc = static_cast<MyClonesArray&>(fData);
-
-        loc.FastRemove(n, fData.GetSize()-1);
-
-        // If it was sorted before it is also sorted now.
-        if (sorted)
-            loc.SetSorted();
-
-        return fData.GetEntriesFast();
-    }
+    Int_t Shrink(Int_t n);
+    void Resize(Int_t n);
 
     // I/O

trunk/MagicSoft/Mars/msimcamera/MSimAPD.cc

@@ -182 +182 @@
 
         // Get arrival time of photon and idx
-        const Double_t t    = ph.GetTime();
-        const UInt_t   idx  = ph.GetTag();
+        const Double_t t = ph.GetTime();
+        const Int_t  idx = ph.GetTag();
+        if (idx<0)
+        {
+            *fLog << err << "ERROR - MSimAPD: Invalid index -1." << endl;
+            return kERROR;
+        }
 
         if (ph.GetWeight()!=1)
@@ -191 +196 @@
             return kERROR;
         }
-
         // Simulate hitting the APD (the signal height in effective
         // "number of photons" is returned)

trunk/MagicSoft/Mars/msimreflector/MSimReflector.cc

@@ -429 +429 @@
     TClonesArray &arr  = fEvt->GetArray();
 
-    TClonesArray &cpy0 = fMirror0->GetArray();
-    //TClonesArray &cpy1 = fMirror1->GetArray();
-    TClonesArray &cpy2 = fMirror2->GetArray();
-    TClonesArray &cpy3 = fMirror3->GetArray();
-    TClonesArray &cpy4 = fMirror4->GetArray();
-    cpy0.ExpandCreateFast(arr.GetEntriesFast());
-    //cpy1.ExpandCreateFast(arr.GetEntriesFast());
-    cpy2.ExpandCreateFast(arr.GetEntriesFast());
-    cpy3.ExpandCreateFast(arr.GetEntriesFast());
-    cpy4.ExpandCreateFast(arr.GetEntriesFast());
+    // Because we knwo in advance what the maximum storage space could
+    // be we allocated it in advance (or shrink it if it was extremely
+    // huge before)
+    // Note, that the drawback is that an extremly large event
+    //       will take about five times its storage space
+    //       for a moment even if a lot from it is unused.
+    //       It will be freed in the next step.
+    fMirror0->Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData
+    fMirror2->Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData
+    fMirror3->Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData
+    fMirror4->Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData
 
     // Initialize mirror properties
@@ -448 +449 @@
 
     // Rotation matrix to derotate sky
+    // For the new coordinate system see the Wiki
     TRotation rot;    // The signs are positive because we align the incident point on ground to the telescope axis
     rot.RotateZ( az); // Rotate point on ground to align it with the telescope axis
     rot.RotateX(-zd); // tilt the point from ground to make it parallel to the mirror plane
 
-    // Now: viewed from the backside of the mirror
-    //  x is pointing downwards
-    //  y is pointing left
-
-    // Rotate around z counterclockwise
-    // rot.RotateZ(TMath::Pi()/2); // Rotate x to point right / y downwards / z from mirror to camera
-    // rot.RotateX(TMath::Pi());   // Flip -y and y (also flips Z :( )
-
-    // Now: viewed from the backside of the mirror
-    //  x is pointing upwards
-    //  y is pointing right
-
+    // Now get the impact point from Corsikas output
     const TVector3 impact(fEvtHeader->GetX(), fEvtHeader->GetY(), 0);
 
@@ -500 +491 @@
         dat->SetDirection(w);
 
-        *static_cast<MPhotonData*>(cpy0.UncheckedAt(cnt[0]++)) = *dat;
+        (*fMirror0)[cnt[0]++] = *dat;
+        //*static_cast<MPhotonData*>(cpy0.UncheckedAt(cnt[0]++)) = *dat;
 
         // Check if the photon has hit the camera housing and holding
@@ -507 +499 @@
 
         // FIXME: Do we really need this one??
+        //(*fMirror1)[cnt[1]++] = *dat;
         //*static_cast<MPhotonData*>(cpy1.UncheckedAt(cnt[1]++)) = *dat;
 
@@ -513 +506 @@
             continue;
 
-        *static_cast<MPhotonData*>(cpy2.UncheckedAt(cnt[2]++)) = *dat;
+        (*fMirror2)[cnt[2]++] = *dat;
+        //*static_cast<MPhotonData*>(cpy2.UncheckedAt(cnt[2]++)) = *dat;
 
         // Now execute the reflection of the photon on the mirrors' surfaces
@@ -526 +520 @@
         dat->SetDirection(w);
 
-        *static_cast<MPhotonData*>(cpy3.UncheckedAt(cnt[3]++)) = *dat;
+        (*fMirror3)[cnt[3]++] = *dat;
+        //*static_cast<MPhotonData*>(cpy3.UncheckedAt(cnt[3]++)) = *dat;
 
         // Propagate the photon along its trajectory to the focal plane z=F
@@ -534 +529 @@
         dat->SetPosition(p);
 
-        *static_cast<MPhotonData*>(cpy4.UncheckedAt(cnt[4]++)) = *dat;
+        (*fMirror4)[cnt[4]++] = *dat;
+        //*static_cast<MPhotonData*>(cpy4.UncheckedAt(cnt[4]++)) = *dat;
 
         // FIXME: It make make sense to move this out of this class