Index: trunk/MagicSoft/Mars/Changelog =================================================================== --- trunk/MagicSoft/Mars/Changelog (revision 9348) +++ trunk/MagicSoft/Mars/Changelog (revision 9349) @@ -35,4 +35,19 @@ * 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. Index: trunk/MagicSoft/Mars/NEWS =================================================================== --- trunk/MagicSoft/Mars/NEWS (revision 9348) +++ trunk/MagicSoft/Mars/NEWS (revision 9349) @@ -69,4 +69,11 @@ 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 Index: trunk/MagicSoft/Mars/msim/MPhotonEvent.cc =================================================================== --- trunk/MagicSoft/Mars/msim/MPhotonEvent.cc (revision 9348) +++ trunk/MagicSoft/Mars/msim/MPhotonEvent.cc (revision 9349) @@ -126,4 +126,101 @@ 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,10 +237,72 @@ } -/* -const char *MPhotonEvent::GetClassName() const -{ - return static_cast(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(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(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,7 +318,14 @@ // Do not modify this. It is optimized for execution // speed and flexibility! - TObject *o = n(fData).FastShrink(n+1); + } + else + { + o=fData.New(n); + } return static_cast(*o); } @@ -175,4 +341,12 @@ { return Add(GetNumPhotons()); +} + +void MPhotonEvent::Sort(Bool_t force) +{ + if (force) + fData.UnSort(); + + static_cast(fData).UncheckedSort(); /*Sort(GetEntriesFast())*/ } @@ -205,5 +379,5 @@ MPhotonData *MPhotonEvent::GetFirst() const { - return static_cast(fData.First()); + return fData.GetEntriesFast()==0 ? 0 : static_cast(fData.First()); } @@ -214,5 +388,5 @@ MPhotonData *MPhotonEvent::GetLast() const { - return static_cast(fData.Last()); + return fData.GetEntriesFast()==0 ? 0 : static_cast(fData.Last()); } @@ -282,4 +456,74 @@ 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(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 +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,8 +25,6 @@ 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,26 +48,6 @@ 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(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 Index: trunk/MagicSoft/Mars/msimcamera/MSimAPD.cc =================================================================== --- trunk/MagicSoft/Mars/msimcamera/MSimAPD.cc (revision 9348) +++ trunk/MagicSoft/Mars/msimcamera/MSimAPD.cc (revision 9349) @@ -182,6 +182,11 @@ // 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,5 +196,4 @@ return kERROR; } - // Simulate hitting the APD (the signal height in effective // "number of photons" is returned) Index: trunk/MagicSoft/Mars/msimreflector/MSimReflector.cc =================================================================== --- trunk/MagicSoft/Mars/msimreflector/MSimReflector.cc (revision 9348) +++ trunk/MagicSoft/Mars/msimreflector/MSimReflector.cc (revision 9349) @@ -429,14 +429,15 @@ 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,20 +449,10 @@ // 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,5 +491,6 @@ dat->SetDirection(w); - *static_cast(cpy0.UncheckedAt(cnt[0]++)) = *dat; + (*fMirror0)[cnt[0]++] = *dat; + //*static_cast(cpy0.UncheckedAt(cnt[0]++)) = *dat; // Check if the photon has hit the camera housing and holding @@ -507,4 +499,5 @@ // FIXME: Do we really need this one?? + //(*fMirror1)[cnt[1]++] = *dat; //*static_cast(cpy1.UncheckedAt(cnt[1]++)) = *dat; @@ -513,5 +506,6 @@ continue; - *static_cast(cpy2.UncheckedAt(cnt[2]++)) = *dat; + (*fMirror2)[cnt[2]++] = *dat; + //*static_cast(cpy2.UncheckedAt(cnt[2]++)) = *dat; // Now execute the reflection of the photon on the mirrors' surfaces @@ -526,5 +520,6 @@ dat->SetDirection(w); - *static_cast(cpy3.UncheckedAt(cnt[3]++)) = *dat; + (*fMirror3)[cnt[3]++] = *dat; + //*static_cast(cpy3.UncheckedAt(cnt[3]++)) = *dat; // Propagate the photon along its trajectory to the focal plane z=F @@ -534,5 +529,6 @@ dat->SetPosition(p); - *static_cast(cpy4.UncheckedAt(cnt[4]++)) = *dat; + (*fMirror4)[cnt[4]++] = *dat; + //*static_cast(cpy4.UncheckedAt(cnt[4]++)) = *dat; // FIXME: It make make sense to move this out of this class