Changeset 9941 for trunk

Timestamp:

09/24/10 15:35:55 (14 years ago)

Author:

tbretz

Message:

Implemented reading of a single telescope from an eventio file.

Location:

trunk/Mars

Files:

• Changelog (modified) (1 diff)
• mcorsika/MCorsikaEvtHeader.cc (modified) (1 diff)
• mcorsika/MCorsikaEvtHeader.h (modified) (2 diffs)
• mcorsika/MCorsikaFormat.cc (modified) (10 diffs)
• mcorsika/MCorsikaFormat.h (modified) (4 diffs)
• mcorsika/MCorsikaRead.cc (modified) (3 diffs)
• mcorsika/MCorsikaRead.h (modified) (2 diffs)
• msim/MPhotonData.cc (modified) (5 diffs)
• msim/MPhotonData.h (modified) (1 diff)
• msim/MPhotonEvent.h (modified) (2 diffs)

trunk/Mars/Changelog

@@ -18 +18 @@
 
                                                  -*-*- END OF LINE -*-*-
+
+ 2010/09/24 Thomas Bretz
+
+   * mcorsika/MCorsikaEvtHeader.cc:
+     - store impact parameters persistent for further use (because in
+       eventio format the header is not repeated this is needed)
+
+   * mcorsika/MCorsikaFormat.[h,cc]:
+     - adapted to be able to read all telescopes and all reused shower
+       events in eventio
+     - a lot of cosmetics
+     - many more comments
+     - PRELIMINARY
+
+   * mcorsika/MCorsikaRead.[h,cc]:
+     - implemented variable for selection of the telescope in eventio format
+     - implemented eventio format (PRELIMINARY)
+
+   * msim/MPhotonData.[h,cc], msim/MPhotonEvent.[h,cc]:
+     - implemented function to calculate mean event time
+     - implemented function to simulate the wavelength (lambda^-2)
+     - implemented function to shift all photons by a certain xy
+     - adapted ReadCorsikaEvt to changes in MCorsikaFormat
+
+
 
  2010/09/24 Daniela Dorner

trunk/Mars/mcorsika/MCorsikaEvtHeader.cc

@@ -160 +160 @@
     // f[96] // Number of reuse of event [max=20]
 
+    memcpy(fTempX, f +97, 20*sizeof(Float_t));
+    memcpy(fTempY, f+117, 20*sizeof(Float_t));
+
     // FIXME: Check fNumReuse<20
-    fX =  f[117 + fNumReuse];
-    fY = -f[ 97 + fNumReuse];
+//    fX =  f[117 + fNumReuse];
+//    fY = -f[ 97 + fNumReuse];
 
     fWeightedNumPhotons = 0;

trunk/Mars/mcorsika/MCorsikaEvtHeader.h

@@ -40 +40 @@
     Float_t  fWeightedNumPhotons;     // weighted number of photons arriving at observation level
 
+    Float_t fTempX[20];               //! Temporary storage for impact parameter
+    Float_t fTempY[20];               //! Temporary storage for impact parameter
+
 public:
     MCorsikaEvtHeader(const char *name=NULL, const char *title=NULL);
@@ -67 +70 @@
     void ResetNumReuse() { fNumReuse=(UInt_t)-1; }
 
+    void InitXY() { fX=fTempY[fNumReuse]; fY=-fTempX[fNumReuse]; }
+
     Int_t ReadEvt(MCorsikaFormat *informat);    // read in event header block
     Int_t ReadEvtEnd(MCorsikaFormat *informat); // read in event end block

trunk/Mars/mcorsika/MCorsikaFormat.cc

@@ -198 +198 @@
 // returns kFALSE;
 //
-Bool_t MCorsikaFormatRaw::GetNextEvent(Float_t ** buffer, Bool_t & readError)
+Int_t MCorsikaFormatRaw::GetNextEvent(Float_t **buffer, UInt_t)
 {
     static Float_t data[273];
@@ -204 +204 @@
 
     if (next == data + 273)
-        {
+    {
         // read next block of events
         if (!ReadData(273, data))
-            {
-            readError = kTRUE;
-            return kFALSE;
-            }
+            return kERROR;
 
         if (!memcmp(data, "EVTE", 4))
-            {
+        {
             // we found the end of list of events
             UnreadLastData();
             return kFALSE;
-            }
+        }
 
         next = data;
-        }
+    }
 
     *buffer = next;
@@ -226 +223 @@
 
     return kTRUE;
-
 }
 
@@ -259 +255 @@
 
         if (fIn->eof())
-            {
-            gLog << err << "ERROR - Missing identifier: " << id  <<
-                   " type: " << type << endl;
+        {
+            gLog << err << "MCorsikaFormatEventIO::SeekNextBlock: ERROR - unexpected end-of-file." << endl;
             return kFALSE;
-            }
-
-        unsigned short fileType = blockHeader[1] & 0xFFFF;
-        if (type == fileType)
+        }
+
+        const unsigned short objType = blockHeader[1] & 0xFFFF;
+        if (type == objType)
             break;
 
-         if (type == 1202 && fileType == 1210)
-            // we are looking for a event header, but found a runEnd. 
-            // This will happen at the end of a run. We stop here looking for
-            // the next event header
+        // we are looking for a event header, but found a runEnd.
+        // This will happen at the end of a run. We stop here looking for
+        // the next event header
+        if (type == 1202 && objType == 1210)
             return kFALSE;
 
         // a unknown block, we jump to the next one
-//gLog << "unknown: " <<  id << " type: " << fileType << " sub-blocks: " <<  (blockHeader[3]>>29);
-//gLog <<  " length: " << (blockHeader[3] & 0x3fffffff) <<   "  pos: " << fIn->tellg() << endl;
         const int length = blockHeader[3] & 0x3fffffff;
-        
-        fIn->seekg(length - 2 * (Int_t)(sizeof(int)), ios::cur);    
+
+        fIn->seekg(length - 2 * (Int_t)(sizeof(int)), ios::cur);
     }
 
@@ -307 +300 @@
 
     // is the RUNE block at the very end of the file?
-    std::streampos currentPos = fIn->tellg();
+    const streampos currentPos = fIn->tellg();
 
     fIn->seekg(-32, ios::end);
@@ -332 +325 @@
     UnreadLastHeader();
     fRunePos = fIn->tellg();
+
     return kTRUE;
 }
@@ -342 +336 @@
 // returns kFALSE;                                                           
 //
-Bool_t MCorsikaFormatEventIO::GetNextEvent(Float_t ** buffer,
-                                           Bool_t & readError) 
-{
-
-    static Float_t * data = NULL;
-    static Float_t * next;
-    static Int_t   topLevelLength = 0;
-    static Int_t   eventLength = 0;
-
-
-    while (eventLength == 0)
-        {
+Int_t MCorsikaFormatEventIO::GetNextEvent(Float_t **buffer, UInt_t telescope)
+{
+    static Float_t *data = NULL;
+    static Int_t    lengthEvent      = -1;
+    static Int_t    lengthPhotonData = 0;
+
+    while (lengthPhotonData == 0)
+    {
         delete [] data;
         data = NULL;
 
-        if (topLevelLength == 0)
-            {
-            if (!NextTopLevelBlock(topLevelLength, readError))
-                return kFALSE;
-            }
-
-        if (!NextEventBlock(eventLength, readError))
-            return kFALSE;
-        topLevelLength -= eventLength + 3 * sizeof(int);
-
-        // read next block of events
-        data = new Float_t [eventLength / sizeof(Float_t)];
-        if (!ReadData(eventLength / sizeof(Float_t), data, 0))
-            {
+        // If we are at the end of the event signal this and
+        // start next time with a new event
+        if (lengthEvent==0)
+        {
+            lengthEvent=-1;
+            return kFALSE; // Signal to start with a new event (process task list first)
+        }
+
+        // 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;
-            readError = kTRUE;
-            return kFALSE;
-            }
-        next = data + 3;        
-        eventLength -= 3 * sizeof(Float_t);
-        }
-
-    eventLength -= 8 * sizeof(Float_t);
-    *buffer = next;
-    next += 8;
-
+            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.
+    lengthPhotonData -= 8 * sizeof(Float_t);
+
+    // Return the pointer to the start of the current photon data
+    *buffer = data + 3;
 
     return kTRUE;
@@ -394 +417 @@
 // EventEnd block (1209). In this case kFALSE is returned
 //
-Bool_t MCorsikaFormatEventIO::NextTopLevelBlock(Int_t & length, 
-                                                Bool_t & readError) const
-{
-    Int_t blockHeader[4];
-
+Int_t MCorsikaFormatEventIO::NextEventObject(Int_t &length) const
+{
     while (1)
     {
@@ -405 +425 @@
         //         - identification field
         //         - length
+
+        UInt_t blockHeader[4];
         fIn->read((char*)blockHeader, 4 * sizeof(Int_t));
 
         if (fIn->eof())
-            {
-            gLog << err << "ERROR - Missing identifier: 1204 or 1209" << endl;
-            readError = kTRUE;
-            return kFALSE;
-            }
-
+        {
+            gLog << err << "MCorsikaFormatEventIO::NextEventObject: ERROR - Unexpected end-of-file." << endl;
+            return kERROR;
+        }
+
+        // For speed reasons we skip the check of the synchronization marker
+
+        // Decode the object type
+        const unsigned short objType = blockHeader[1] & 0xFFFF;
+
+        // Decode length of block
         length = blockHeader[3] & 0x3fffffff;
-        const unsigned short fileType = blockHeader[1] & 0xFFFF;
-        if (fileType == 1204)
+
+        // Check if we found the next array (reuse / impact parameter)
+        // blockHeader[2] == array number (reuse)
+        if (objType==1204)
             return kTRUE;
 
-        if (fileType == 1209)
-            {
-            // we found an eventEnd block, reset file pointer
-            fIn->seekg( -4 * (Int_t)(sizeof(Int_t)), ios::cur);
-            length = 0;
-            return kFALSE;
-            }
-
+        // we found an eventEnd block, reset file pointer
+        if (objType==1209)
+            break;
 
         // a unknown block, we jump to the next one
-        fIn->seekg(length, ios::cur);    
-    }
-
-    return kTRUE;
-}
-
-// --------------------------------------------------------------------------
-//
-Bool_t MCorsikaFormatEventIO::NextEventBlock(Int_t & length, 
-                                             Bool_t & readError) const
-{
-    Int_t blockHeader[3];
+        fIn->seekg(length, ios::cur);
+    }
+
+    fIn->seekg( -4 * (Int_t)(sizeof(Int_t)), ios::cur);
+    length = 0;
+
+    return kFALSE;
+}
+
+// --------------------------------------------------------------------------
+//
+Int_t MCorsikaFormatEventIO::NextPhotonObject(Int_t &length) const
+{
+    UInt_t blockHeader[3];
 
     // we read - synchronisation marker
@@ -446 +472 @@
     //         - identification field
     //         - length
-    fIn->read((char*)blockHeader, 3 * sizeof(Int_t));
+    fIn->read((char*)blockHeader, 3 * sizeof(UInt_t));
 
     if (fIn->eof())
-        {
-        gLog << err << "ERROR - Missing identifier: 1205" << endl;
-        readError = kTRUE;
-        return kFALSE;
-        }
-    
-    const unsigned short fileType = blockHeader[0] & 0xFFFF;
-    if (fileType != 1205)
-        {
-        gLog << err << "ERROR - Unexpected type: " << fileType << "expected 1205" << endl;
-        readError = kTRUE;
-        return kFALSE;
-        }
+    {
+        gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected end-of-file.." << endl;
+        return -1;
+    }
+
+    const unsigned short objType = blockHeader[0] & 0xFFFF;
+    if (objType != 1205)
+    {
+        gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected EventIO object type " << objType << " (expected 1205)." << endl;
+        return -1;
+    }
 
     const unsigned short version = (blockHeader[0] >> 20) & 0x0FFF;
     if (version != 0)
-        {
-        gLog << err << "ERROR - Unexpected version: " << version << "expected: 0" << endl;
-        readError = kTRUE;
-        return kFALSE;
-        }
+    {
+        gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected EventIO object version " << version << " (expected 0 for object type 1205)." << endl;
+        return -1;
+    }
 
     length = blockHeader[2] & 0x3fffffff;
 
-    return kTRUE;
-
-}
+    // blockHeader[1] == 1000 x array number (reuse) + telescope number
+    return blockHeader[1] % 1000;
+}

trunk/Mars/mcorsika/MCorsikaFormat.h

@@ -2 +2 @@
 #define MARS_MDataFormat
 
+#ifndef MARS_MAGIC
+#include "MAGIC.h"
+#endif
 
 #ifndef ROOT_Rtypes
@@ -35 +38 @@
    virtual void   ResetPos() const;
 
-   virtual Bool_t GetNextEvent(Float_t ** buffer, Bool_t & readError) = 0;
+   virtual Int_t  GetNextEvent(Float_t **buffer, UInt_t tel=0) = 0;
    virtual Bool_t IsEventioFormat() const = 0;
 
@@ -59 +62 @@
    Bool_t SeekEvtEnd();
 
-   Bool_t GetNextEvent(Float_t ** buffer, Bool_t & readError);
+   Int_t  GetNextEvent(Float_t **buffer, UInt_t);
    Bool_t IsEventioFormat() const {return kFALSE;}
 };
@@ -67 +70 @@
 {
 private:
-        std::streampos fRunePos; // file position of the RUNE block
+    std::streampos fRunePos; // file position of the RUNE block
 
 public:
-   MCorsikaFormatEventIO(std::istream * in)
+    MCorsikaFormatEventIO(std::istream *in)
         : MCorsikaFormat(in) {fRunePos = std::streampos(0);}
 
-   Bool_t SeekNextBlock(const char * id, unsigned short type) const;
-   void   UnreadLastHeader() const;
+    Bool_t SeekNextBlock(const char *id, unsigned short type) const;
+    void   UnreadLastHeader() const;
 
-   Bool_t SeekEvtEnd();
+    Bool_t SeekEvtEnd();
 
-   Bool_t GetNextEvent(Float_t ** buffer, Bool_t & readError);
-   Bool_t IsEventioFormat() const {return kTRUE;}
+    Int_t  GetNextEvent(Float_t **buffer, UInt_t tel);
+    Bool_t IsEventioFormat() const { return kTRUE; }
 
 private:
-   Bool_t NextTopLevelBlock(Int_t & length, Bool_t & readError) const;
-   Bool_t NextEventBlock(Int_t & length, Bool_t & readError) const;
-
+    Int_t  NextEventObject(Int_t &length) const;
+    Int_t  NextPhotonObject(Int_t &length) const;
 };
 

trunk/Mars/mcorsika/MCorsikaRead.cc

@@ -96 +96 @@
 //
 MCorsikaRead::MCorsikaRead(const char *fname, const char *name, const char *title)
-    : fRunHeader(0), fEvtHeader(0), fEvent(0), /*fEvtData(0),*/ fForceMode(kFALSE),
+    : fRunHeader(0), fEvtHeader(0), fEvent(0), fTelescopeIdx(0), fForceMode(kFALSE),
     fFileNames(0), fNumFile(0), fNumEvents(0), fNumTotalEvents(0),
     fIn(0),  fInFormat(0), fParList(0)
@@ -212 +212 @@
 
     const char *expname = gSystem->ExpandPathName(name);
-    fInFormat = CorsikaFormatFactory(fLog, expname);
+    fInFormat = MCorsikaFormat::CorsikaFormatFactory(expname);
     delete [] expname;
 
@@ -401 +401 @@
         fInFormat->StorePos();
         fEvtHeader->ResetNumReuse();
+
+        if (fInFormat->IsEventioFormat())
+        {
+            // Read the event header
+            const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat);
+            if (rc1!=kTRUE)
+                return rc1;
+        }
+
     }
     else
-        fInFormat->ResetPos();
+    {
+        if (!fInFormat->IsEventioFormat())
+            fInFormat->ResetPos();
+    }
+
+    if (!fInFormat->IsEventioFormat())
+    {
+        // Read the event header
+        const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat);
+        if (rc1!=kTRUE)
+            return rc1;
+
+    }
 
     fEvtHeader->IncNumReuse();
-
-    // Read the event header
-    const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat);
-    if (rc1!=kTRUE)
-        return rc1;
-
-    // Check if reusage number should be reset and increase it
-    //  Note, that the trick here is that after reset it is set to -1
+    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
-    const Int_t rc2 = fEvent->ReadCorsikaEvt(fInFormat, fEvtHeader->GetNumReuse()+1);
+    //  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

trunk/Mars/mcorsika/MCorsikaRead.h

@@ -22 +22 @@
     MPhotonEvent      *fEvent;      //! event information
 
+    UInt_t          fTelescopeIdx;  // Index of telescope to be extracted
     Bool_t          fForceMode;     // Force mode skipping defect RUNE
 
@@ -52 +53 @@
     ~MCorsikaRead();
 
-    //static Byte_t IsFileValid(const char *name);
-
-    //void SetInterleave(UInt_t i) { fInterleave = i; }
     void SetForceMode(Bool_t b=kTRUE) { fForceMode=b; }
+    void SetTelescopeIdx(UInt_t idx=0) { fTelescopeIdx = idx; }
 
     TString GetFullFileName() const;

trunk/Mars/msim/MPhotonData.cc

@@ -47 +47 @@
 
 #include <TMath.h>
+#include <TRandom.h>
 
 #include "MLog.h"
@@ -171 +172 @@
 // --------------------------------------------------------------------------
 //
+// Set the wavelength to a random lambda^-2 distributed value
+// between wmin and wmax.
+//
+void MPhotonData::SimWavelength(Float_t wmin, Float_t wmax)
+{
+    const Double_t w = gRandom->Uniform(wmin, wmax);
+
+    fWavelength = TMath::Nint(wmin*wmax / w);
+}
+
+
+// --------------------------------------------------------------------------
+//
 // Set the data member according to the 8 floats read from a reflector-file.
 // This function MUST reset all data-members, no matter whether these are
@@ -255 +269 @@
 Int_t MPhotonData::FillEventIO(Float_t f[8])
 {
-    if (TMath::Nint(f[6])!=1)
-    {
-        gLog << err << "ERROR - Bunch sizes != 1 are not supported." << endl;
-        return kFALSE;
-    }
+    // photons in this bunch
+    const UInt_t n = TMath::Nint(f[6]);
+    if (n==0)
+        return 0;
+
+    f[6] = n-1;
 
     fPosX             =  f[1];              // xpos relative to telescope [cm]
@@ -267 +282 @@
     fTime             =  f[4];              // a relative arival time [ns]
     fProductionHeight =  f[5];              // altitude of emission [cm]
-//    fNumPhotons       =  TMath::Nint(f[6]); // photons in this bunch
     fWavelength       =  TMath::Nint(f[7]); // so far always zeor = unspec. [nm]
 
@@ -275 +289 @@
     fWeight     =  1;
 
-    return kTRUE;
-}
+    return n;
+}
+
 /*
 // --------------------------------------------------------------------------

trunk/Mars/msim/MPhotonData.h

@@ -101 +101 @@
     void SetTime(Double_t t)  { fTime  = t; }
 
+    void SimWavelength(Float_t wmin, Float_t wmax);
+
     void Copy(TObject &obj) const;
 

trunk/Mars/msim/MPhotonEvent.h

@@ -39 +39 @@
     Double_t GetMeanX() const;
     Double_t GetMeanY() const;
+    Double_t GetMeanT() const;
 
     TClonesArray &GetArray() { return fData; }
@@ -54 +55 @@
     Int_t Shrink(Int_t n);
     void Resize(Int_t n);
+
+    void AddXY(Double_t x, Double_t y);
+    void SimWavelength(Float_t wmin, Float_t wmax);
 
     // I/O