Changeset 9349 for trunk


Ignore:
Timestamp:
02/19/09 16:27:33 (16 years ago)
Author:
tbretz
Message:
*** empty log message ***
Location:
trunk/MagicSoft/Mars
Files:
6 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Mars/Changelog

    r9348 r9349  
    3535   * msim/MSimAbsorption.cc:
    3636     - don't check wavelength range when theta should be used
     37
     38   * msimcamera/MSimAPD.cc:
     39     - check for uninitialized indices
     40
     41   * msim/MPhotonEvent.[h,cc]:
     42     - moved the code for MyClonesArray to the source file
     43     - improved a lot the reading speed by reading larger blocks
     44       of data from the file at once
     45     - improved memory handling. This ensures that even the largest
     46       events don't fill the memory forever and the allocated memory
     47       is free'd again after some time
     48
     49   * msimreflector/MSimReflector.cc:
     50     - Use the new Resize function of MPhotonEvent to make sure
     51       that the memory is not allocated forever.
    3752
    3853
  • trunk/MagicSoft/Mars/NEWS

    r9347 r9349  
    6969     layout
    7070
     71   * Improved reading speed for corsika files
     72
     73   * Improved memory handling (if a large ampount of memory was needed
     74     for a single event all further events were stored in the same
     75     memory and it was never freed, so the program took this memory
     76     until the end)
     77 
    7178 ;star
    7279
  • trunk/MagicSoft/Mars/msim/MPhotonEvent.cc

    r9342 r9349  
    126126using namespace std;
    127127
     128// ==========================================================================
     129
     130class MyClonesArray : public TClonesArray
     131{
     132public:
     133    TObject **FirstRef() { return fCont; }
     134
     135    // --------------------------------------------------------------------------
     136    //
     137    // This is an extremly optimized version of ExpandCreateFast. It only resets
     138    // the marker for the last element (fLast) to n-1 and doen't change anything
     139    // else. This implicitly assumes that the stored objects don't allocate
     140    // memory. It does not necessarily mean that the slots after fLast
     141    // are empty (set to 0). This is what is assumed in the TClonesArray.
     142    // We also don't call Changed() because it would reset Sorted. If the
     143    // array was sorted before it is also sorted now. You MUST make sure
     144    // that you only set n in a range for which valid entries have been
     145    // created before (e.g. by ExpandCreateFast).
     146    //
     147    void FastShrink(Int_t n)
     148    {
     149        fLast = n - 1;
     150    }
     151
     152    // --------------------------------------------------------------------------
     153    //
     154    // This is a optimized (faster) version of Delete which deletes consequtive
     155    // entries from index idx1 to idx2 (both included) and calls their
     156    // destructor. Note that there is no range checking done!
     157    //
     158    void FastRemove(Int_t idx1, Int_t idx2)
     159    {
     160        // Remove object at index idx.
     161
     162        //if (!BoundsOk("RemoveAt", idx1)) return 0;
     163        //if (!BoundsOk("RemoveAt", idx2)) return 0;
     164
     165        Long_t dtoronly = TObject::GetDtorOnly();
     166
     167        idx1 -= fLowerBound;
     168        idx2 -= fLowerBound;
     169
     170        for (TObject **obj=fCont+idx1; obj<=fCont+idx2; obj++)
     171        {
     172            if (!*obj)
     173                continue;
     174
     175            if ((*obj)->TestBit(kNotDeleted)) {
     176                // Tell custom operator delete() not to delete space when
     177                // object fCont[i] is deleted. Only destructors are called
     178                // for this object.
     179                TObject::SetDtorOnly(*obj);
     180                delete *obj;
     181            }
     182
     183            *obj = 0;
     184            // recalculate array size
     185        }
     186        TObject::SetDtorOnly((void*)dtoronly);
     187
     188        if (idx1<=fLast && fLast<=idx2)
     189        {
     190            do {
     191                fLast--;
     192            } while (fLast >= 0 && fCont[fLast] == 0);
     193        }
     194
     195        Changed();
     196    }
     197
     198
     199    //void SetSorted() { fSorted = kTRUE; }
     200
     201    // --------------------------------------------------------------------------
     202    //
     203    // This is an optimized version of Sort which doesn't check the
     204    // IsSortable flag before. It only sorts the entries from 0
     205    // to GetEntriesFast().
     206    //
     207    void UncheckedSort()
     208    {
     209        if (fSorted)
     210            return;
     211
     212        const Int_t nentries = GetEntriesFast();
     213        if (nentries <= 0)
     214            return;
     215
     216        QSort(GetObjectRef(First()), fKeep->GetObjectRef(fKeep->First()),
     217              0, TMath::Min(nentries, kMaxInt-fLowerBound));
     218
     219        fSorted = kTRUE;
     220    }
     221};
     222
     223// ==========================================================================
     224
    128225// --------------------------------------------------------------------------
    129226//
     
    140237}
    141238
    142 /*
    143 const char *MPhotonEvent::GetClassName() const
    144 {
    145     return static_cast<TObject*>(fData.GetClass())->GetName();
    146 }
    147 */
     239// --------------------------------------------------------------------------
     240//
     241// This is an extremly optimized version of ExpandCreateFast. It only resets
     242// the marker for the last element (fLast) to n-1 and doen't change anything
     243// else. This has the advantage that the allocated memory is kept but only
     244// valid entries are written to a file.
     245//
     246// If the array was sorted before it is also sorted now. You MUST make sure
     247// that you only set n in a range for which valid entries have been
     248// created before (e.g. by ExpandCreateFast).
     249//
     250Int_t MPhotonEvent::Shrink(Int_t n)
     251{
     252    /*
     253    // The number of object written by the streamer is defined by
     254    // GetEntriesFast(), i.e. this number must be shrinked to
     255    // the real array size. We use RemoveAt instead of ExpandCreate
     256    // because RemoveAt doesn't free memory. Thus in the next
     257    // iteration it can be reused and doesn't need to be reallocated.
     258    // Do not change this. It is optimized for execution speed
     259    //        for (int i=fData.GetSize()-1; i>=n; i--)
     260    //          fData.RemoveAt(i);
     261    const Bool_t sorted = fData.IsSorted();
     262
     263    MyClonesArray &loc = static_cast<MyClonesArray&>(fData);
     264
     265    loc.FastRemove(n, fData.GetSize()-1);
     266
     267    // If it was sorted before it is also sorted now.
     268    if (sorted)
     269        loc.SetSorted();
     270    */
     271
     272    // fData.ExpandCreateFast(n);  // Just set fLast = n -1
     273
     274    // Just set fLast = n -1
     275    static_cast<MyClonesArray&>(fData).FastShrink(n);
     276    return fData.GetEntriesFast();
     277}
     278
     279// --------------------------------------------------------------------------
     280//
     281// The resized the array. If it has to be increased in size it is done
     282// with ExpandCreateFast. If it should be shrinked it is done with
     283// ExpandCreateFast if n>fData.GetSize()/100 or n==0. This keeps the allocated
     284// memory and just sets the marker for the last element in the array (fLast).
     285//
     286// If the allocated memory grew to huge we reset the allocated memory
     287// by calling ExpandCreate(n) (frees the allocated storage for the
     288// objects) and Expand(n) (frees the allocated memory for the list
     289// of pointers to the objects)
     290//
     291// 100 hundred is an arbitrary number taking into account that todays
     292// computers have a lot of memory and we don't want to free and allocate
     293// new memory too often.
     294//
     295// In priciple there might be more clever methods to handle the memory.
     296//
     297void MPhotonEvent::Resize(Int_t n)
     298{
     299    if (n==0 || n*100>fData.GetSize())
     300        fData.ExpandCreateFast(n);  // Just set fLast = n -1
     301    else
     302    {
     303        fData.ExpandCreate(n);      // Free memory of allocated MPhotonData
     304        fData.Expand(n);            // Free memory of allocated pointers
     305    }
     306}
    148307
    149308// --------------------------------------------------------------------------
     
    159318    // Do not modify this. It is optimized for execution
    160319    // speed and flexibility!
    161     TObject *o = n<fData.GetSize() && fData.UncheckedAt(n) ? fData.UncheckedAt(n) : fData.New(n);
    162 
    163     // Now we read a single cherenkov bunch
     320    TObject *o = 0;
     321    if (n<fData.GetSize() && fData.UncheckedAt(n))
     322    {
     323        o=fData.UncheckedAt(n);
     324        static_cast<MyClonesArray&>(fData).FastShrink(n+1);
     325    }
     326    else
     327    {
     328        o=fData.New(n);
     329    }
    164330    return static_cast<MPhotonData&>(*o);
    165331}
     
    175341{
    176342    return Add(GetNumPhotons());
     343}
     344
     345void MPhotonEvent::Sort(Bool_t force)
     346{
     347    if (force)
     348        fData.UnSort();
     349
     350    static_cast<MyClonesArray&>(fData).UncheckedSort(); /*Sort(GetEntriesFast())*/
    177351}
    178352
     
    205379MPhotonData *MPhotonEvent::GetFirst() const
    206380{
    207     return static_cast<MPhotonData*>(fData.First());
     381    return fData.GetEntriesFast()==0 ? 0 : static_cast<MPhotonData*>(fData.First());
    208382}
    209383
     
    214388MPhotonData *MPhotonEvent::GetLast() const
    215389{
    216     return static_cast<MPhotonData*>(fData.Last());
     390    return fData.GetEntriesFast()==0 ? 0 : static_cast<MPhotonData*>(fData.Last());
    217391}
    218392
     
    282456    Int_t n = 0;
    283457
     458    // --- old I/O ---
     459    // Read only + Reflector (no absorption)
     460    // Muons:   1.06GB/115s =  9.2MB/s (100kEvs)
     461    // Gammas:  1.57GB/275s =  5.7MB/s (  1kEvs)
     462
     463    // Read only:
     464    // Gammas:  1.57GB/158s =  9.9MB/s (  1kEvs)
     465    // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)
     466
     467    // --- new I/O ---
     468    // Read only (don't allocate storage space):
     469    // Gammas:  1.57GB/143s = 11.0MB/s (  1kEvs)
     470    // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)
     471
     472    // Read only in blocks (with storage allocation):
     473    // Gammas:  1.57GB/28s  =  56MB/s (  1kEvs)
     474    // Muons:   1.06GB/5.2s = 204MB/s (100kEvs)
     475
     476    // Read only in blocks (without storage allocation):
     477    // similar to just copy
     478
     479    // Copy with cp
     480    // 1.57GB/ 5s   CPU
     481    // 1.57GB/28s   REAL
     482    // 1.06GB/ 3s   CPU
     483    // 1.06GB/22s   REAL
     484
     485    while (1)
     486    {
     487        // Stprage for one block
     488        char c[273*4];
     489
     490        // Read the first four byte to check whether the next block
     491        // doen't belong to the event anymore
     492        fin.read(c, 4);
     493        if (!fin)
     494            return kFALSE;
     495
     496        // Check if the event is finished
     497        if (!memcmp(c, "EVTE", 4))
     498            break;
     499
     500        // Now read the rest of the data
     501        fin.read(c+4, 272*4);
     502
     503        Float_t *ptr = reinterpret_cast<Float_t*>(c);
     504        Float_t *end = ptr + 273;
     505
     506        // Loop over all sub-blocks (photons) in the block and if
     507        // they contain valid data add them to the array
     508        while (ptr<end)
     509        {
     510            // Get/Add the n-th entry from the array and
     511            // fill it with the current 7 floats
     512            const Int_t rc = Add(n).FillCorsika(ptr);
     513            ptr += 7;
     514
     515            switch (rc)
     516            {
     517            case kCONTINUE:  continue;        // No data in this bunch... skip it.
     518            case kERROR:     return kERROR;   // Error occured
     519            //case kFALSE:     return kFALSE;   // End of stream
     520            }
     521
     522            // This is a photon we would like to keep later.
     523            // Increase the counter by one
     524            n++;
     525        }
     526    }
     527/*
    284528    while (1)
    285529    {
     
    326570        n++;
    327571    }
    328 
    329     Shrink(n);
     572  */
     573
     574    Resize(n);
    330575    fData.UnSort();
    331576
  • trunk/MagicSoft/Mars/msim/MPhotonEvent.h

    r9342 r9349  
    1212#include <iosfwd>
    1313
     14using namespace std;
     15
    1416class MPhotonData;
    1517class MCorsikaRunHeader;
    16 
    17 class MyClonesArray : public TClonesArray
    18 {
    19 public:
    20     TObject **FirstRef() { return fCont; }
    21 
    22     void FastRemove(Int_t idx1, Int_t idx2)
    23     {
    24         // Remove object at index idx.
    25 
    26         //if (!BoundsOk("RemoveAt", idx1)) return 0;
    27         //if (!BoundsOk("RemoveAt", idx2)) return 0;
    28 
    29         Long_t dtoronly = TObject::GetDtorOnly();
    30 
    31         idx1 -= fLowerBound;
    32         idx2 -= fLowerBound;
    33 
    34         for (TObject **obj=fCont+idx1; obj<=fCont+idx2; obj++)
    35         {
    36             if (!*obj)
    37                 continue;
    38 
    39             if ((*obj)->TestBit(kNotDeleted)) {
    40                 // Tell custom operator delete() not to delete space when
    41                 // object fCont[i] is deleted. Only destructors are called
    42                 // for this object.
    43                 TObject::SetDtorOnly(*obj);
    44                 delete *obj;
    45             }
    46 
    47             *obj = 0;
    48             // recalculate array size
    49         }
    50         TObject::SetDtorOnly((void*)dtoronly);
    51 
    52         if (idx1<=fLast && fLast<=idx2)
    53         {
    54             do {
    55                 fLast--;
    56             } while (fLast >= 0 && fCont[fLast] == 0);
    57         }
    58 
    59         Changed();
    60     }
    61 
    62     void SetSorted() { fSorted = kTRUE; }
    63 };
    6418
    6519class MPhotonEvent : public MParContainer
     
    7125    MPhotonEvent(const char *name=NULL, const char *title=NULL);
    7226
    73     // TObjArray, FIXME: This could be improved if checking for IsSortable is omitted
    74     void Sort(Bool_t force=kFALSE) { if (force) fData.UnSort(); fData.Sort(); }
     27    void Sort(Bool_t force=kFALSE);
    7528    Bool_t IsSorted() const { return fData.IsSorted(); }
    76 
    7729
    7830    // Getter/Setter
     
    9648    const MPhotonData &operator[](UInt_t idx) const;
    9749
    98     Int_t Shrink(Int_t n)
    99     {
    100         // The number of object written by the streamer is defined by
    101         // GetEntriesFast(), i.e. this number must be shrinked to
    102         // the real array size. We use RemoveAt instead of ExpandCreate
    103         // because RemoveAt doesn't free memory. Thus in the next
    104         // iteration it can be reused and doesn't need to be reallocated.
    105         // Do not change this. It is optimized for execution speed
    106         //        for (int i=fData.GetSize()-1; i>=n; i--)
    107         //          fData.RemoveAt(i);
    108         const Bool_t sorted = fData.IsSorted();
    109 
    110         MyClonesArray &loc = static_cast<MyClonesArray&>(fData);
    111 
    112         loc.FastRemove(n, fData.GetSize()-1);
    113 
    114         // If it was sorted before it is also sorted now.
    115         if (sorted)
    116             loc.SetSorted();
    117 
    118         return fData.GetEntriesFast();
    119     }
     50    Int_t Shrink(Int_t n);
     51    void Resize(Int_t n);
    12052
    12153    // I/O
  • trunk/MagicSoft/Mars/msimcamera/MSimAPD.cc

    r9328 r9349  
    182182
    183183        // Get arrival time of photon and idx
    184         const Double_t t    = ph.GetTime();
    185         const UInt_t   idx  = ph.GetTag();
     184        const Double_t t = ph.GetTime();
     185        const Int_t  idx = ph.GetTag();
     186        if (idx<0)
     187        {
     188            *fLog << err << "ERROR - MSimAPD: Invalid index -1." << endl;
     189            return kERROR;
     190        }
    186191
    187192        if (ph.GetWeight()!=1)
     
    191196            return kERROR;
    192197        }
    193 
    194198        // Simulate hitting the APD (the signal height in effective
    195199        // "number of photons" is returned)
  • trunk/MagicSoft/Mars/msimreflector/MSimReflector.cc

    r9336 r9349  
    429429    TClonesArray &arr  = fEvt->GetArray();
    430430
    431     TClonesArray &cpy0 = fMirror0->GetArray();
    432     //TClonesArray &cpy1 = fMirror1->GetArray();
    433     TClonesArray &cpy2 = fMirror2->GetArray();
    434     TClonesArray &cpy3 = fMirror3->GetArray();
    435     TClonesArray &cpy4 = fMirror4->GetArray();
    436     cpy0.ExpandCreateFast(arr.GetEntriesFast());
    437     //cpy1.ExpandCreateFast(arr.GetEntriesFast());
    438     cpy2.ExpandCreateFast(arr.GetEntriesFast());
    439     cpy3.ExpandCreateFast(arr.GetEntriesFast());
    440     cpy4.ExpandCreateFast(arr.GetEntriesFast());
     431    // Because we knwo in advance what the maximum storage space could
     432    // be we allocated it in advance (or shrink it if it was extremely
     433    // huge before)
     434    // Note, that the drawback is that an extremly large event
     435    //       will take about five times its storage space
     436    //       for a moment even if a lot from it is unused.
     437    //       It will be freed in the next step.
     438    fMirror0->Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData
     439    fMirror2->Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData
     440    fMirror3->Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData
     441    fMirror4->Resize(arr.GetEntriesFast()); // Free memory of allocated MPhotonData
    441442
    442443    // Initialize mirror properties
     
    448449
    449450    // Rotation matrix to derotate sky
     451    // For the new coordinate system see the Wiki
    450452    TRotation rot;    // The signs are positive because we align the incident point on ground to the telescope axis
    451453    rot.RotateZ( az); // Rotate point on ground to align it with the telescope axis
    452454    rot.RotateX(-zd); // tilt the point from ground to make it parallel to the mirror plane
    453455
    454     // Now: viewed from the backside of the mirror
    455     //  x is pointing downwards
    456     //  y is pointing left
    457 
    458     // Rotate around z counterclockwise
    459     // rot.RotateZ(TMath::Pi()/2); // Rotate x to point right / y downwards / z from mirror to camera
    460     // rot.RotateX(TMath::Pi());   // Flip -y and y (also flips Z :( )
    461 
    462     // Now: viewed from the backside of the mirror
    463     //  x is pointing upwards
    464     //  y is pointing right
    465 
     456    // Now get the impact point from Corsikas output
    466457    const TVector3 impact(fEvtHeader->GetX(), fEvtHeader->GetY(), 0);
    467458
     
    500491        dat->SetDirection(w);
    501492
    502         *static_cast<MPhotonData*>(cpy0.UncheckedAt(cnt[0]++)) = *dat;
     493        (*fMirror0)[cnt[0]++] = *dat;
     494        //*static_cast<MPhotonData*>(cpy0.UncheckedAt(cnt[0]++)) = *dat;
    503495
    504496        // Check if the photon has hit the camera housing and holding
     
    507499
    508500        // FIXME: Do we really need this one??
     501        //(*fMirror1)[cnt[1]++] = *dat;
    509502        //*static_cast<MPhotonData*>(cpy1.UncheckedAt(cnt[1]++)) = *dat;
    510503
     
    513506            continue;
    514507
    515         *static_cast<MPhotonData*>(cpy2.UncheckedAt(cnt[2]++)) = *dat;
     508        (*fMirror2)[cnt[2]++] = *dat;
     509        //*static_cast<MPhotonData*>(cpy2.UncheckedAt(cnt[2]++)) = *dat;
    516510
    517511        // Now execute the reflection of the photon on the mirrors' surfaces
     
    526520        dat->SetDirection(w);
    527521
    528         *static_cast<MPhotonData*>(cpy3.UncheckedAt(cnt[3]++)) = *dat;
     522        (*fMirror3)[cnt[3]++] = *dat;
     523        //*static_cast<MPhotonData*>(cpy3.UncheckedAt(cnt[3]++)) = *dat;
    529524
    530525        // Propagate the photon along its trajectory to the focal plane z=F
     
    534529        dat->SetPosition(p);
    535530
    536         *static_cast<MPhotonData*>(cpy4.UncheckedAt(cnt[4]++)) = *dat;
     531        (*fMirror4)[cnt[4]++] = *dat;
     532        //*static_cast<MPhotonData*>(cpy4.UncheckedAt(cnt[4]++)) = *dat;
    537533
    538534        // FIXME: It make make sense to move this out of this class
Note: See TracChangeset for help on using the changeset viewer.