Changeset 9943 for trunk/Mars/mcorsika

Timestamp:

09/28/10 10:11:31 (14 years ago)

Author:

tbretz

Message:

Implemented a working version of eventio which reads all telescope's data and all array's data.

Location:

trunk/Mars/mcorsika

Files:

• MCorsikaEvtHeader.h (modified) (1 diff)
• MCorsikaFormat.cc (modified) (9 diffs)
• MCorsikaFormat.h (modified) (4 diffs)
• MCorsikaRead.cc (modified) (10 diffs)
• MCorsikaRead.h (modified) (5 diffs)

trunk/Mars/mcorsika/MCorsikaEvtHeader.h

@@ -71 +71 @@
 
     void InitXY() { fX=fTempY[fNumReuse]; fY=-fTempX[fNumReuse]; }
+    void AddXY(Float_t x, Float_t y) { fX+=x; fY+=y; }
 
     Int_t ReadEvt(MCorsikaFormat *informat);    // read in event header block

trunk/Mars/mcorsika/MCorsikaFormat.cc

@@ -72 +72 @@
 }
 
+void MCorsikaFormat::Read(void *ptr, Int_t i) const
+{
+    fIn->read((char*)ptr, i);
+}
+
 // --------------------------------------------------------------------------
 //
@@ -122 +127 @@
 {
     fIn->seekg(fPos, ios::beg);
+}
+
+void MCorsikaFormat::Rewind() const
+{
+    fIn->seekg(0, ios::beg);
 }
 
@@ -198 +208 @@
 // returns kFALSE;
 //
-Int_t MCorsikaFormatRaw::GetNextEvent(Float_t **buffer, UInt_t)
+Int_t MCorsikaFormatRaw::GetNextEvent(Float_t **buffer, Int_t)
 {
     static Float_t data[273];
@@ -242 +252 @@
 Bool_t MCorsikaFormatEventIO::SeekNextBlock(const char * id, unsigned short type) const
 {
-    int blockHeader[7];
-
+    cout << "Seek " << type << endl;
     while (1)
     {
@@ -252 +261 @@
         //         - unknown field
         //         - id (first 4 bytes of data field)
-        fIn->read((char*)blockHeader, 6 * sizeof(int));
+        int blockHeader[4];
+        fIn->read((char*)blockHeader, 4 * sizeof(int));
 
         if (fIn->eof())
         {
-            gLog << err << "MCorsikaFormatEventIO::SeekNextBlock: ERROR - unexpected end-of-file." << endl;
+            gLog << err << "MCorsikaFormatEventIO::SeekNextBlock: ERROR - Unexpected end-of-file seeking " << id << " (" << type << ")." << endl;
             return kFALSE;
         }
 
         const unsigned short objType = blockHeader[1] & 0xFFFF;
+        cout << "objType=" << objType << endl;
         if (type == objType)
-            break;
+        {
+            if (!id)
+                break;
+
+            char c[9];
+            fIn->read(c, 8);
+
+            if (memcmp(c+4, id, 4)==0)
+                break;
+
+            c[8] = 0;
+            gLog << err << "MCorsikaFormatEventIO::SeekNextBlock: ERROR - '" << c+4 << "' doesn't match expected one '" << id << "' in object " << type << endl;
+            return kFALSE;
+        }
 
         // we are looking for a event header, but found a runEnd.
@@ -272 +296 @@
         // a unknown block, we jump to the next one
         const int length = blockHeader[3] & 0x3fffffff;
-
-        fIn->seekg(length - 2 * (Int_t)(sizeof(int)), ios::cur);
+        fIn->seekg(length, ios::cur);
     }
 
@@ -336 +359 @@
 // returns kFALSE;                                                           
 //
-Int_t MCorsikaFormatEventIO::GetNextEvent(Float_t **buffer, UInt_t telescope)
+Int_t MCorsikaFormatEventIO::GetNextEvent(Float_t **buffer, Int_t telescope)
 {
     static Float_t *data = NULL;
-    static Int_t    lengthEvent      = -1;
-    static Int_t    lengthPhotonData = 0;
-
-    while (lengthPhotonData == 0)
+
+    static int lengthPhotonData = 0;
+    static int lengthEvent      = 0;
+
+    if (lengthPhotonData>0)
+    {
+        cout << "Return Bunch l2=" << lengthPhotonData << endl;
+
+        lengthPhotonData -= 8 * sizeof(Float_t);
+        *buffer += 8; // Return the pointer to the start of the current photon data
+        return kTRUE;
+    }
+
+    if (data)
     {
         delete [] data;
         data = NULL;
-
-        // If we are at the end of the event signal this and
-        // start next time with a new event
-        if (lengthEvent==0)
+        cout << "Return end-of-tel LE=" << lengthEvent << endl;
+        return kFALSE;
+    }
+
+    if (lengthEvent==0)
+    {
+        cout << "Searching 1204... " << flush;
+        // The length of the block is stored and we use it to determine
+        // when the next top level block is expected
+        const Int_t rc = NextEventObject(lengthEvent);
+        if (rc==kERROR)
+            return kERROR;
+        if (!rc)
         {
-            lengthEvent=-1;
-            return kFALSE; // Signal to start with a new event (process task list first)
+            cout << "skip EVTE." << endl;
+            return kFALSE;
         }
 
-        // Look for the next top level block (an "event", object type 1204)
-        // A top level block contains the data of a full array simulation
-        // (a single re-usage of the shower)
-        if (lengthEvent<0)
-        {
-            // The length of the block is stored and we use it to determine
-            // when the next top level block is expected
-            const Int_t rc = NextEventObject(lengthEvent);
-            if (rc!=kTRUE)
-                return rc;
-        }
-
-        // Look for the next event photon bunch (object type 1205)
-        const Int_t tel = NextPhotonObject(lengthPhotonData);
-        if (tel<0)
-            return kERROR;
-
-        // Decrease the distance to the next "event" by the current photon bunches
-        lengthEvent -= lengthPhotonData + 3 * sizeof(int);
-
-        // If the telescope is not the one we want to read skip the photon data
-        if (UInt_t(tel)!=telescope)
-        {
-            fPrevPos = fIn->tellg();
-            fIn->seekg(lengthPhotonData, ios::cur);
-            lengthPhotonData=0;
-            continue;
-        }
-
-        const UInt_t cnt = lengthPhotonData / sizeof(Float_t);
-
-        // Read next object (photon bunch)
-        data = new Float_t[cnt];
-        if (!ReadData(cnt, data, 0))
-        {
-            delete [] data;
-            data = NULL;
-            return kERROR;
-        }
-
-        // The object containing the photon bunches starts with:
-        // Array[short], Telescope[Short], NumPhotons[Real], NumBunches[Long]
-        // eventLength contains now the number of bytes with photon data not yet evaluated
-        lengthPhotonData -= 3 * sizeof(Float_t);
-    }
-
-    // The photon data itself is 8 floats long. When we return we expect this to be
-    // evaluated already when the function is called the next time.
+        cout << "found new array." << endl;
+    }
+
+    cout << "Searching 1205..." << flush;
+
+    // Look for the next event photon bunch (object type 1205)
+    const Int_t tel = NextPhotonObject(lengthPhotonData);
+    if (tel<0)
+        return kERROR;
+
+    lengthEvent -= lengthPhotonData+12; // Length of data + Length of header
+
+    lengthPhotonData -= 12;     // Length of data before the photon bunches
+    fIn->seekg(12, ios::cur);   // Skip this data
+
+    cout << " found (tel=" << tel << ",LEN=" << lengthEvent << ",len=" << lengthPhotonData << ")... " << flush;
+
+    if (lengthPhotonData==0)
+    {
+        cout << "empty!" << endl;
+        return kFALSE;
+    }
+
+    const UInt_t cnt = lengthPhotonData / sizeof(Float_t);
+    // Read next object (photon bunch)
+    data = new Float_t[cnt];
+    if (!ReadData(cnt, data, 0))
+    {
+        delete [] data;
+        data = NULL;
+        return kERROR;
+    }
+
+    cout << "return!" << endl;
+
     lengthPhotonData -= 8 * sizeof(Float_t);
-
-    // Return the pointer to the start of the current photon data
-    *buffer = data + 3;
+    *buffer = data; // Return the pointer to the start of the current photon data
 
     return kTRUE;
+}
 }
 
@@ -425 +453 @@
         //         - identification field
         //         - length
-
         UInt_t blockHeader[4];
         fIn->read((char*)blockHeader, 4 * sizeof(Int_t));
@@ -440 +467 @@
         const unsigned short objType = blockHeader[1] & 0xFFFF;
 
+        cout << "objType=" << objType << endl;
+
         // Decode length of block
         length = blockHeader[3] & 0x3fffffff;

trunk/Mars/mcorsika/MCorsikaFormat.h

@@ -37 +37 @@
    virtual void   StorePos();
    virtual void   ResetPos() const;
+   virtual void   Rewind() const;
 
-   virtual Int_t  GetNextEvent(Float_t **buffer, UInt_t tel=0) = 0;
+   virtual Int_t  GetNextEvent(Float_t **buffer, Int_t tel=0) = 0;
    virtual Bool_t IsEventioFormat() const = 0;
 
@@ -44 +45 @@
 
    std::streampos GetCurrPos() const;
+
+   void Read(void *ptr, Int_t i=4) const;
 
    static MCorsikaFormat *CorsikaFormatFactory(const char *fileName);
@@ -62 +65 @@
    Bool_t SeekEvtEnd();
 
-   Int_t  GetNextEvent(Float_t **buffer, UInt_t);
+   Int_t  GetNextEvent(Float_t **buffer, Int_t);
    Bool_t IsEventioFormat() const {return kFALSE;}
 };
@@ -81 +84 @@
     Bool_t SeekEvtEnd();
 
-    Int_t  GetNextEvent(Float_t **buffer, UInt_t tel);
+    Int_t  GetNextEvent(Float_t **buffer, Int_t tel);
     Bool_t IsEventioFormat() const { return kTRUE; }
 

trunk/Mars/mcorsika/MCorsikaRead.cc

@@ -96 +96 @@
 //
 MCorsikaRead::MCorsikaRead(const char *fname, const char *name, const char *title)
-    : fRunHeader(0), fEvtHeader(0), fEvent(0), fTelescopeIdx(0), fForceMode(kFALSE),
+    : fRunHeader(0), fEvtHeader(0), fEvent(0), fTelescopeIdx(-1), fForceMode(kFALSE),
     fFileNames(0), fNumFile(0), fNumEvents(0), fNumTotalEvents(0),
-    fIn(0),  fInFormat(0), fParList(0)
+    /*fIn(0),*/  fInFormat(0), fParList(0), fNumTelescopes(1)
 {
     fName  = name  ? name  : "MRead";
@@ -117 +117 @@
 {
     delete fFileNames;
-    if (fIn)
-        delete fIn;
+//    if (fIn)
+//        delete fIn;
     if (fInFormat)
         delete fInFormat;
@@ -191 +191 @@
     // open the input stream and check if it is really open (file exists?)
     //
-    if (fIn)
+/*    if (fIn)
         delete fIn;
     fIn = NULL;
-
+*/
     if (fInFormat)
        delete fInFormat;
@@ -239 +239 @@
 //    if (!fEvtHeader->ReadRunHeader(*fIn, *fRunHeader))
 //        return kERROR;
+
+    fNumTelescopes = 1;
+    if (fInFormat->IsEventioFormat())
+    {
+        fInFormat->StorePos();
+        fInFormat->Rewind();
+
+        if (!fInFormat->SeekNextBlock(0, 1201))
+        {
+            *fLog << err << "ERROR - Object 1201 not found in file." << endl;
+            return kERROR;
+        }
+
+        fInFormat->Read(&fNumTelescopes);
+
+        if (fTelescopeIdx>=fNumTelescopes)
+        {
+            *fLog << err << "ERROR - Requested telescope index " << fTelescopeIdx << " exceeds number of telescopes " << fNumTelescopes << " in file." << endl;
+            return kERROR;
+        }
+
+        fTelescopeX.Set(fNumTelescopes);
+        fTelescopeY.Set(fNumTelescopes);
+        fTelescopeZ.Set(fNumTelescopes);
+        fTelescopeR.Set(fNumTelescopes);
+
+        fInFormat->Read(fTelescopeX.GetArray(), 4*fNumTelescopes);
+        fInFormat->Read(fTelescopeY.GetArray(), 4*fNumTelescopes);
+        fInFormat->Read(fTelescopeZ.GetArray(), 4*fNumTelescopes);
+        fInFormat->Read(fTelescopeR.GetArray(), 4*fNumTelescopes);
+
+        *fLog << all;
+        *fLog << "Number of telescopes: " << fNumTelescopes << " (using ";
+        if (fTelescopeIdx>=0)
+            *fLog << "telescope " << fTelescopeIdx;
+        else
+            *fLog << "all telescopes";
+        *fLog << ")" << endl;
+
+        const Int_t lo = fTelescopeIdx<0 ? 0              : fTelescopeIdx;
+        const Int_t up = fTelescopeIdx<0 ? fNumTelescopes : fTelescopeIdx+1;
+
+        for (int i=lo; i<up; i++)
+        {
+            *fLog << inf << "Telescope #" << setw(4) << i << " [X/Y/Z (R)]: ";
+            *fLog << setw(7) << fTelescopeX[i] << "/";
+            *fLog << setw(7) << fTelescopeY[i] << "/";
+            *fLog << setw(7) << fTelescopeZ[i] << "  (R=" << setw(4) << fTelescopeR[i] << ")" << endl;
+        }
+
+        fNumTelescope = 0;
+
+        fInFormat->ResetPos();
+    }
+
 
     fInFormat->StorePos();
@@ -319 +374 @@
             }
 
-            fNumTotalEvents += fRunHeader->GetNumEvents()*fRunHeader->GetNumReuse();
+            fNumTotalEvents += fRunHeader->GetNumEvents()*fRunHeader->GetNumReuse()*
+                (fTelescopeIdx<0 ? fNumTelescopes : 1);
             continue;
         }
         break;
     }
-
+/*
     if (fIn)
         delete fIn;
     fIn = NULL;
-
+*/
     return rc;
 }
@@ -356 +412 @@
     //
     fParList = 0;
-
+/*
     if (!OpenStream())
         return kFALSE;
-
+*/
     //
     //  check if all necessary containers exist in the Parameter list.
@@ -395 +451 @@
 Int_t MCorsikaRead::ReadEvent()
 {
-    // Store the position to re-read a single event several times
-    //  FIXME: Does this work with EventIO, too?
-    if (fEvtHeader->GetNumReuse()>=fRunHeader->GetNumReuse()-1)
-    {
-        fInFormat->StorePos();
-        fEvtHeader->ResetNumReuse();
-
-        if (fInFormat->IsEventioFormat())
+    if (fInFormat->IsEventioFormat())
+    {
+        if (fNumTelescope==0)
         {
-            // Read the event header
             const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat);
             if (rc1!=kTRUE)
                 return rc1;
+
+            if (fEvtHeader->GetNumReuse()==fRunHeader->GetNumReuse()-1)
+                fEvtHeader->ResetNumReuse();
+            fEvtHeader->IncNumReuse();
+
+            cout << "===== Array " << fEvtHeader->GetNumReuse() << " =====" << endl;
         }
 
+        const Int_t rc2 = fEvent->ReadCorsikaEvt(fInFormat, fTelescopeIdx);
+        if (rc2!=kTRUE)
+            return rc2;
+
+        fEvtHeader->InitXY();
+
+        const Double_t x = fTelescopeX[fNumTelescope];
+        const Double_t y = fTelescopeY[fNumTelescope];
+
+        fEvtHeader->AddXY(y, -x);
+        fEvent->AddXY(fEvtHeader->GetX(), fEvtHeader->GetY());
+        fEvent->SimWavelength(fRunHeader->GetWavelengthMin(), fRunHeader->GetWavelengthMax());
+
+        fNumTelescope++;
+        fNumTelescope%=fNumTelescopes;
+
+        return kTRUE;
     }
     else
     {
-        if (!fInFormat->IsEventioFormat())
+        // Store the position to re-read a single event several times
+        if (fEvtHeader->GetNumReuse()<fRunHeader->GetNumReuse()-1)
             fInFormat->ResetPos();
-    }
-
-    if (!fInFormat->IsEventioFormat())
-    {
+        else
+        {
+            fInFormat->StorePos();
+            fEvtHeader->ResetNumReuse();
+        }
+
         // Read the event header
         const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat);
@@ -424 +500 @@
             return rc1;
 
-    }
-
-    fEvtHeader->IncNumReuse();
-    fEvtHeader->InitXY();
-
-    // In the case of corsika raw data (MMCS only) we have too loop over one event
-    // several times to read all impact parameters (reusage of events)
-    // In case of the eventio format we have to decide, the data of which telescope
-    // we want to extract
-    const Int_t id = fInFormat->IsEventioFormat() ? fTelescopeIdx : fEvtHeader->GetNumReuse()+1;
-
-    // Read the photons corresponding to the i-th core location
-    //  EventIO: Number of telescope
-    //  Raw:     Number of re-use
-    const Int_t rc2 = fEvent->ReadCorsikaEvt(fInFormat, id);
-    if (rc2!=kTRUE)
-        return rc2;
-
-    //----
-    if (fInFormat->IsEventioFormat())
-    {
-        fEvent->AddXY(fEvtHeader->GetX(), fEvtHeader->GetY());
-        fEvent->SimWavelength(fRunHeader->GetWavelengthMin(), fRunHeader->GetWavelengthMax());
-        return kTRUE;
-    }
-    //----
-
-
-    // read event end
-    return fEvtHeader->ReadEvtEnd(fInFormat);
+        fEvtHeader->IncNumReuse();
+        fEvtHeader->InitXY();
+
+        // In the case of corsika raw data (MMCS only) we have too loop over one event
+        // several times to read all impact parameters (reusage of events)
+        // In case of the eventio format we have to decide, the data of which telescope
+        // we want to extract
+
+        // Read the photons corresponding to the i-th core location
+        //  EventIO: Number of telescope
+        //  Raw:     Number of re-use
+        const Int_t rc2 = fEvent->ReadCorsikaEvt(fInFormat, fEvtHeader->GetNumReuse()+1);
+        if (rc2!=kTRUE)
+            return rc2;
+
+        // read event end
+        return fEvtHeader->ReadEvtEnd(fInFormat);
+    }
 }
 
@@ -502 +566 @@
 Int_t MCorsikaRead::PostProcess()
 {
-    const UInt_t n = fNumEvents*fRunHeader->GetNumReuse();
+    const UInt_t n = fNumEvents*fRunHeader->GetNumReuse()*(fTelescopeIdx<0?fNumTelescopes:1);
 
     //
@@ -517 +581 @@
     return kTRUE;
 }
+
+// --------------------------------------------------------------------------
+//
+//    Telescope: 1
+//
+// In the case of eventio-format select which telescope to extract from the
+// file. -1 will extract all telescopes
+//
+Int_t MCorsikaRead::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
+{
+    Bool_t rc = kFALSE;
+    if (IsEnvDefined(env, prefix, "Telescope", print))
+    {
+        rc = kTRUE;
+        fTelescopeIdx = GetEnvValue(env, prefix, "Telescope", fTelescopeIdx);
+    }
+
+    return rc;
+}

trunk/Mars/mcorsika/MCorsikaRead.h

@@ -4 +4 @@
 #ifndef MARS_MRead
 #include "MRead.h"
+#endif
+
+#ifndef ROOT_TArrayF
+#include <TArrayF.h>
 #endif
 
@@ -22 +26 @@
     MPhotonEvent      *fEvent;      //! event information
 
-    UInt_t          fTelescopeIdx;  // Index of telescope to be extracted
+    Int_t           fTelescopeIdx;  // Index of telescope to be extracted
     Bool_t          fForceMode;     // Force mode skipping defect RUNE
 
@@ -30 +34 @@
     UInt_t    fNumTotalEvents; //! total number of events in all files
 
-    ifstream       *fIn;       //! input stream (file to read from)
+//    ifstream       *fIn;       //! input stream (file to read from)
     MCorsikaFormat *fInFormat; //! access to input corsika data
 
     MParList *fParList;        //! tasklist to call ReInit from
 
+    Int_t  fNumTelescopes;     //! Number of telescopes in array
+    Int_t  fNumTelescope;      //! Number of telescope currently being read
+    TArrayF fTelescopeX;       //! x pos (measured towards north, unit: cm)
+    TArrayF fTelescopeY;       //! y pos (measured towards west,  unit: cm)
+    TArrayF fTelescopeZ;       //! z pos (from detection level,   unit: cm)
+    TArrayF fTelescopeR;       //! Radii of spheres around tel. (unit: cm)
+
     //UInt_t    fInterleave;
     //Bool_t    fForce;
 
-    virtual Bool_t OpenStream() { return kTRUE; }
+//    virtual Bool_t OpenStream() { return kTRUE; }
 
     Bool_t ReadEvtEnd();
@@ -45 +56 @@
     Int_t  ReadEvent();
 
+    // MTask
     Int_t PreProcess(MParList *pList);
     Int_t Process();
     Int_t PostProcess();
+
+    // MParContainer
+    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
 
 public:
@@ -54 +69 @@
 
     void SetForceMode(Bool_t b=kTRUE) { fForceMode=b; }
-    void SetTelescopeIdx(UInt_t idx=0) { fTelescopeIdx = idx; }
+    void SetTelescopeIdx(Int_t idx=-1) { fTelescopeIdx = idx; }
 
     TString GetFullFileName() const;