/* ======================================================================== *\ ! ! * ! * This file is part of CheObs, the Modular Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appears in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): Thomas Bretz, 12/2000 ! Author(s): Qi Zhe, 06/2007 ! ! Copyright: CheObs Software Development, 2000-2009 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // MPhotonEvent // // Storage container to store photon collections // // The class is designed to be extremely fast which is important taking into // account the extremely high number of photons. This has some impacts on // its handling. // // The list has to be kept consistent, i.e. without holes. // // There are two ways to achieve this: // // a) Use RemoveAt to remove an entry somewhere // b) Compress() the TClonesArray afterwards // // Compress is not the fastes, so there is an easier way. // // a) When you loop over the list and want to remove an entry copy all // following entry backward in the list, so that the hole will // be created at its end. // b) Call Shrink(n) with n the number of valid entries in your list. // // To loop over the TClonesArray you can use a TIter which has some // unnecessary overhead and therefore is slower than necessary. // // Since the list is kept consistent you can use a simple loop saving // a lot of CPU time taking into account the high number of calls to // TObjArrayIter::Next which you would create. // // Here is an example (how to remove every second entry) // // --------------------------------------------------------------------- // // Int_t cnt = 0; // // const Int_t num = event->GetNumPhotons(); // for (Int_t idx=0; idxShrink(cnt); // // ---------------------------------- or ------------------------------- // // TClonesArray &arr = MPhotonEvent->GetArray(); // // Int_t cnt = 0; // // const Int_t num = arr.GetEntriesFast(); // for (Int_t idx=0; idx(arr.UncheckedAt(idx)); // // *static_cast(arr.UncheckedAt(cnt++)) = *dat; // } // // MPhotonEvent->Shrink(cnt); // // --------------------------------------------------------------------- // // The flag for a sorted array is for speed reasons not in all conditions // maintained automatically. Especially Add() doesn't reset it. // // So be sure that if you want to sort your array it is really sorted. // // // Version 1: // ---------- // - First implementation // ///////////////////////////////////////////////////////////////////////////// #include "MPhotonEvent.h" #include #include #include #include "MLog.h" #include "MLogManip.h" #include "MPhotonData.h" ClassImp(MPhotonEvent); using namespace std; // -------------------------------------------------------------------------- // // Default constructor. It initializes all arrays with zero size. // MPhotonEvent::MPhotonEvent(const char *name, const char *title) : fData("MPhotonData", 1) { fName = name ? name : "MPhotonEvent"; fTitle = title ? title : "Corsika Event Data Information"; fData.SetBit(TClonesArray::kForgetBits); fData.BypassStreamer(kFALSE); } /* const char *MPhotonEvent::GetClassName() const { return static_cast(fData.GetClass())->GetName(); } */ // -------------------------------------------------------------------------- // // If n is smaller than the current allocated array size a reference to // the n-th entry is returned, otherwise an entry at n is created // calling the default constructor. Note, that there is no range check // but it is not recommended to call this function with // n>fData.GetSize() // MPhotonData &MPhotonEvent::Add(Int_t n) { // Do not modify this. It is optimized for execution // speed and flexibility! TObject *o = n(*o); } // -------------------------------------------------------------------------- // // Add a new photon (MPhtonData) at the end of the array. // In this case the default constructor of MPhotonData is called. // // A reference to the new object is returned. // MPhotonData &MPhotonEvent::Add() { return Add(GetNumPhotons()); } // -------------------------------------------------------------------------- // // Get the i-th photon from the array. Not, for speed reasons there is no // range check so you are responsible that you do not excess the number // of photons (GetNumPhotons) // MPhotonData &MPhotonEvent::operator[](UInt_t idx) { return *static_cast(fData.UncheckedAt(idx)); } // -------------------------------------------------------------------------- // // Get the i-th photon from the array. Not, for speed reasons there is no // range check so you are responsible that you do not excess the number // of photons (GetNumPhotons) // const MPhotonData &MPhotonEvent::operator[](UInt_t idx) const { return *static_cast(fData.UncheckedAt(idx)); } // -------------------------------------------------------------------------- // // Return a pointer to the first photon if available. // MPhotonData *MPhotonEvent::GetFirst() const { return static_cast(fData.First()); } // -------------------------------------------------------------------------- // // Return a pointer to the last photon if available. // MPhotonData *MPhotonEvent::GetLast() const { return static_cast(fData.Last()); } // -------------------------------------------------------------------------- // // Read the Event section from the file // Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin) { Int_t n = 0; while (1) { // Check the first four bytes char c[4]; fin.read(c, 4); // End of stream if (!fin) return kFALSE; // Check if we found the end of the event if (!memcmp(c, "EVTE", 4)) break; // The first for byte contained data already --> go back fin.seekg(-4, ios::cur); // Do not modify this. It is optimized for execution // speed and flexibility! MPhotonData &ph = Add(n); // It checks how many entries the lookup table has. If it has enough // entries and the entry was already allocated, we can re-use it, // otherwise we have to allocate it. // Now we read a single cherenkov bunch. Note that for speed reason we have not // called the constructor if the event was already constructed (virtual table // set), consequently we must make sure that ReadCorsikaEvent does reset // all data mebers no matter whether they are read or not. const Int_t rc = ph.ReadCorsikaEvt(fin); // Evaluate result from reading event switch (rc) { case kCONTINUE: continue; // No data in this bunch... skip it. case kFALSE: return kFALSE; // End of stream case kERROR: return kERROR; // Error occured } // FIXME: If fNumPhotons!=1 add the photon more than once // Now increase the number of entries which are kept, // i.e. keep this photon(s) n++; } Shrink(n); fData.UnSort(); SetReadyToSave(); //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl; return kTRUE; } // -------------------------------------------------------------------------- // Int_t MPhotonEvent::ReadRflEvt(std::istream &fin) { Int_t n = 0; while (1) { // Check the first four bytes char c[13]; fin.read(c, 13); // End of stream if (!fin) return kFALSE; // Check if we found the end of the event if (!memcmp(c, "\nEND---EVENT\n", 13)) break; // The first for byte contained data already --> go back fin.seekg(-13, ios::cur); // Do not modify this. It is optimized for execution // speed and flexibility! //TObject *o = n(o)->ReadRflEvt(fin); const Int_t rc = Add(n).ReadRflEvt(fin); // Evaluate result from reading event switch (rc) { case kCONTINUE: continue; // No data in this bunch... skip it. case kFALSE: return kFALSE; // End of stream case kERROR: return kERROR; // Error occured } // Now increase the number of entries which are kept, // i.e. keep this photon(s) n++; } Shrink(n); SetReadyToSave(); //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl; return kTRUE; } // -------------------------------------------------------------------------- // // Print the array // void MPhotonEvent::Print(Option_t *) const { fData.Print(); } // ------------------------------------------------------------------------ // // You can call Draw() to add the photons to the current pad. // The photons are painted each tim ethe pad is updated. // Make sure that you use the right (world) coordinate system, // like created, eg. by the MHCamera histogram. // void MPhotonEvent::Paint(Option_t *) { MPhotonData *ph=NULL; TMarker m; m.SetMarkerStyle(kFullDotMedium); // Gtypes.h TIter Next(&fData); while ((ph=(MPhotonData*)Next())) { m.SetX(ph->GetPosY()*10); // north m.SetY(ph->GetPosX()*10); // east m.Paint(); } }