Changeset 10060 for trunk


Ignore:
Timestamp:
11/29/10 14:24:23 (14 years ago)
Author:
rohlfs
Message:
two new command line arguments of readcorsika: -A=arrayNo and -T=telescopeNo. New design of program flow in MCorsikaRead: It is now determined by the order of the data blocks in the input file.
Location:
trunk/Mars
Files:
13 edited

Legend:

Unmodified
Added
Removed

  • trunk/Mars/Changelog

    r10056 r10060  
    1818
    1919                                                 -*-*- END OF LINE -*-*-
     20
     21 2010/11/29 Reiner Rohlfs
     22
     23   * msim/MPhotonData.cc:
     24     - Use all data if telescope array is not defined
     25
     26   * msim/MPhotonEvent.[h,cc]:
     27     - remove following methods:
     28         Int_t ReadCorsikaEvt(MCorsikaFormat *fInFormat, Int_t i);
     29         Int_t ReadCorsikaEvt(istream &fin, Int_t i);
     30     - replace those methods by
     31         Int_t ReadEventIoEvt(MCorsikaFormat *fInFormat);
     32         Int_t ReadCorsikaEvt(Float_t * data, Int_t numEvents, Int_t arrayIdx);
     33 
     34   * mcorsika/MCorsikaRunHeader.h mcorsika/MCorsikaRunHeader.cc
     35     - Split method ReadEvtEnd() into two functions:
     36        ReadEvtEnd() and ReadEventHeader()
     37
     38   * mcorsika/MCorsikaEvtHeader.[h,cc]:
     39     - method ReadEvt does not read data from file by itself, but gets
     40       buffer as input
     41     - test that number of reuse of same shour is not greater than 20
     42     - method ReadEvtEnd does not seek to the EVTE block and does not
     43       read the Bloch header any more.
     44
     45   * mcorsika/MCorsikaFormat.[h,cc]:
     46     - remove following methods:
     47       SeekNextBlock(), UnreadLastHeader(), UnreadLastData(), StorePos(),
     48       ResetPos(), Rewind(), GetNextEvent(), GetCurrPos(),
     49       NextEventObject() and NextPhotonObject()
     50     - new method: NextBlock()
     51
     52   * mcorsika/MCorsikaRead.[h,cc]:
     53     - new design of program flow. It is now determined by the order of
     54       the data blocks in the input file.
     55     - Depending on the variables fTelescopeIdx and fArrayIdx, only data
     56       of the requested telescope and array are stored in the output
     57       file.
     58       
     59   * readcorsika.cc:     
     60     - two new command line arguments: -A=arrayNo and -T=telescopeNo 
     61       The values of these parameters are store in the MCorsikaRead -
     62       class.
     63
    2064
    2165 2010/11/22 Thomas Bretz

  • trunk/Mars/mcorsika/MCorsikaEvtHeader.cc

    r9949 r10060  
    124124// --------------------------------------------------------------------------
    125125//
    126 // read the EVTH information from the input stream
    127 // return FALSE if there is no  header anymore, else TRUE
     126// get the EVTH information from the input parameter f
    128127//
    129 Int_t MCorsikaEvtHeader::ReadEvt(MCorsikaFormat *fInFormat)
     128Int_t MCorsikaEvtHeader::ReadEvt(Float_t * f)
    130129{
    131     const Int_t rc=fInFormat->SeekNextBlock("EVTH", 1202);
    132     if (rc!=kTRUE)
    133         return rc;
    134 
    135     Float_t f[273];
    136     if (!fInFormat->ReadData(272, f))
    137         return kERROR;
    138130
    139131    fEvtNumber  = TMath::Nint(f[0]);
    … …  
    160152
    161153    // f[96] // Number of reuse of event [max=20]
     154    fTotReuse = f[96];
     155    if (fTotReuse > 20)
     156      {
     157      *fLog << err << "Number of reuse of shower is " << fTotReuse
     158                   << " But maximum implemented: 20" << endl;
     159      return kFALSE;
     160      }
    162161
    163162    memcpy(fTempX, f +97, 20*sizeof(Float_t));
    164163    memcpy(fTempY, f+117, 20*sizeof(Float_t));
     164
    165165
    166166    // FIXME: Check fNumReuse<20
    … …  
    170170    fWeightedNumPhotons = 0;
    171171
    172     return fInFormat->Eof() ? kERROR : kTRUE;
     172    return kTRUE;
    173173}
    174174
    … …  
    178178Int_t MCorsikaEvtHeader::ReadEvtEnd(MCorsikaFormat * fInFormat)
    179179{
    180     if (fInFormat->SeekNextBlock("EVTE", 1209)!=kTRUE)
    181         return kERROR;
    182 
    183     //fin.seekg(-1088,ios::cur);
    184 
    185     Float_t f[2];
    186     if (!fInFormat->ReadData(2, f))
     180    Float_t f[272];
     181    if (!fInFormat->Read(f, 272 * sizeof(Float_t)))
    187182        return kERROR;
    188183

  • trunk/Mars/mcorsika/MCorsikaEvtHeader.h

    r9943 r10060  
    2121    UInt_t   fEvtNumber;              // Event number
    2222    UInt_t   fNumReuse;               // Counter of the reuse of the same shower
     23    UInt_t   fTotReuse;               //! number of reuse of the same shower     
    2324//    UInt_t   fParticleID;             // Particle ID (see MMcEvtBasic or CORSIKA manual)
    2425    Float_t  fTotalEnergy;            // [GeV]
    … …  
    5152    UInt_t GetEvtNumber() const { return fEvtNumber; }
    5253    UInt_t GetNumReuse() const { return fNumReuse; }
     54    UInt_t GetTotReuse() const { return fTotReuse; }
    5355//    UInt_t GetParticleID() const { return fParticleID; }
    5456
    … …  
    6769    Double_t GetImpact() const;
    6870
     71    void GetArrayOffset(Int_t arrayIdx, Float_t & xArrOff, Float_t & yArrOff)
     72                     {xArrOff = fTempY[arrayIdx]; yArrOff=-fTempX[arrayIdx]; }
     73    void SetTelescopeOffset(Int_t arrayIdx, Float_t xTelOff, Float_t yTelOff)
     74                     {fNumReuse = arrayIdx; fX = xTelOff; fY = yTelOff;}
     75
    6976    void IncNumReuse() { fNumReuse++; }
    70     void ResetNumReuse() { fNumReuse=(UInt_t)-1; }
     77    void ResetNumReuse() { fNumReuse=0; }
    7178
    7279    void InitXY() { fX=fTempY[fNumReuse]; fY=-fTempX[fNumReuse]; }
    7380    void AddXY(Float_t x, Float_t y) { fX+=x; fY+=y; }
    7481
    75     Int_t ReadEvt(MCorsikaFormat *informat);    // read in event header block
     82    Int_t ReadEvt(Float_t * f);                 // read in event header block
    7683    Int_t ReadEvtEnd(MCorsikaFormat *informat); // read in event end block
    7784

  • trunk/Mars/mcorsika/MCorsikaFormat.cc

    r9949 r10060  
    3636#include "MLogManip.h"
    3737
     38
    3839using namespace std;
     40
    3941
    4042const unsigned int MCorsikaFormat::kSyncMarker = 0xd41f8a37;
    … …  
    5557    }
    5658
    57     char buffer[5]="\0\0\0\0";
     59    char *buffer = new char[5];
     60    memset(buffer, 0, 5);
    5861    fileIn->read(buffer, 4);
    5962    fileIn->seekg(-4, ios::cur);
    6063
    6164    if (strcmp(buffer, "RUNH") == 0)
     65    {
     66        delete [] buffer;
    6267        return new MCorsikaFormatRaw(fileIn);
     68    }
    6369
    6470    if (*reinterpret_cast<unsigned int*>(buffer) == kSyncMarker)
     71    {
     72        delete [] buffer;
    6573        return new MCorsikaFormatEventIO(fileIn);
     74    }
    6675
    6776    gLog << err << "File " << fileName <<
    6877            " is neither a CORSIKA raw nor EventIO file" << endl;
    6978    delete fileIn;
     79    delete [] buffer;
    7080
    7181    return NULL;
    7282}
    7383
    74 void MCorsikaFormat::Read(void *ptr, Int_t i) const
    75 {
    76     fIn->read((char*)ptr, i);
    77 }
    78 
     84Bool_t MCorsikaFormat::Read(void *ptr, Int_t i) const
     85{
     86   fIn->read((char*)ptr, i);
     87   return !fIn->fail();
     88
     89}
    7990// --------------------------------------------------------------------------
    8091//
    … …  
    8697// --------------------------------------------------------------------------
    8798//
    88 streampos MCorsikaFormat::GetCurrPos() const
    89 {
    90     return fIn->tellg();
    91 }
    92 
    93 // --------------------------------------------------------------------------
    94 //
    95 Bool_t MCorsikaFormat::ReadData(Int_t numValues, Float_t * buffer,
    96                                 Int_t minSeekValues)
    97 {
    98     fPrevPos = fIn->tellg();
    99 
    100     fIn->read((char*)buffer, numValues * sizeof(Float_t));
    101 
    102     if (numValues < minSeekValues)
    103         // skip the remaining data of this block
    104         fIn->seekg( (minSeekValues - numValues) * 4, ios::cur);
    105 
    106     return !fIn->fail();
    107 
    108 }
    109 
    110 // --------------------------------------------------------------------------
    111 //
    112 void MCorsikaFormat::UnreadLastData() const
    113 {
    114     fIn->seekg(fPrevPos, ios::beg);
    115 }
    116 
    117 // --------------------------------------------------------------------------
    118 //
    119 void MCorsikaFormat::StorePos()
    120 {
    121     fPos = fIn->tellg();
    122 }
    123 
    124 // --------------------------------------------------------------------------
    125 //
    126 void MCorsikaFormat::ResetPos() const
    127 {
    128     fIn->seekg(fPos, ios::beg);
    129 }
    130 
    131 void MCorsikaFormat::Rewind() const
    132 {
    133     fIn->seekg(0, ios::beg);
    134 }
    135 
    136 // --------------------------------------------------------------------------
    137 //
    13899MCorsikaFormat::~MCorsikaFormat()
    139100{
    … …  
    141102}
    142103
    143 // --------------------------------------------------------------------------
    144 //
    145 // Reads the next 4 bytes which should be the id. (for example "RUNH", RUNE".
    146 // "EVTH")
    147 // Returns kTRUE if the functions finds the id
    148 //         kFALSE if the functions does not find the "id" as the next 4
    149 //                bytes in the file.
    150 // After this call the file position pointer points just after the 4 bytes
    151 // of the id.
    152 //
    153 Int_t MCorsikaFormatRaw::SeekNextBlock(const char * id, unsigned short type) const
    154 {
    155     char blockHeader[5]="\0\0\0\0";
     104
     105// --------------------------------------------------------------------------
     106//
     107// After a call to this function, the file pointer is located after the   
     108// header of the block. As the event block has no header it is located   
     109// at the beginning of this block, which is as the beginning of the data 
     110Bool_t MCorsikaFormatRaw::NextBlock(Bool_t  subBlock,
     111                                    Int_t & blockType,
     112                                    Int_t & blockVersion,
     113                                    Int_t & blockIdentifier,
     114                                    Int_t & blockLength) const
     115{
     116    char blockHeader[5];
     117    blockHeader[4] = 0;
    156118    fIn->read(blockHeader, 4);
    157 
    158     if (strcmp(blockHeader, id)==0)
    159         return kTRUE;
    160 
    161     // at the end of a file we are looking for the next Event header,
    162     // but find the end of a run. This is expected, therefor no error
    163     // message.
    164     if (strcmp(blockHeader, "RUNE")==0)
     119    if (fIn->eof())
    165120        return kFALSE;
    166121
    167     gLog << err << "ERROR - Wrong identifier: " << id << " expected.";
    168     gLog << " But read " << blockHeader << " from file." << endl;
    169 
    170     return kERROR;
    171 }
    172 
    173 // --------------------------------------------------------------------------
    174 //
    175 void MCorsikaFormatRaw::UnreadLastHeader() const
    176 {
    177     fIn->seekg(-4, ios::cur);
    178 }
    179 
     122    blockVersion = 0;
     123    blockIdentifier = 0;
     124    blockLength     = 272 * 4;
     125
     126    if (strcmp(blockHeader, "RUNH") == 0)
     127        blockType = 1200;
     128    else if (strcmp(blockHeader, "RUNE") == 0)
     129        blockType = 1210;
     130    else if (strcmp(blockHeader, "EVTH") == 0)
     131        blockType = 1202;
     132    else if (strcmp(blockHeader, "EVTE") == 0)
     133        blockType = 1209;
     134    else    // the events, they don't have a specific header
     135        {
     136        blockType = 1105;
     137        fIn->seekg(-4, ios::cur);
     138        blockLength += 4;
     139        }
     140
     141    return kTRUE;
     142}
    180143// --------------------------------------------------------------------------
    181144//
    … …  
    192155        if (!strcmp(runh, "RUNE"))
    193156        {
    194             fIn->seekg(-4, ios::cur);
     157//            fIn->seekg(-4, ios::cur);
    195158            return kTRUE;
    196159        }
    … …  
    199162    return kTRUE;
    200163}
    201 
    202 // --------------------------------------------------------------------------
    203 //                                                                           
    204 // Returns one event (7 Float values) after the other until the EventEnd     
    205 // block is found.                                                           
    206 // If a read error occurred, the readError is set to kTRUE and the function 
    207 // returns kFALSE;
    208 //
    209 Int_t MCorsikaFormatRaw::GetNextEvent(Float_t **buffer, Int_t)
    210 {
    211     static Float_t data[273];
    212     static Float_t * next = data + 273;
    213 
    214     if (next == data + 273)
    215     {
    216         // read next block of events
    217         if (!ReadData(273, data))
    218             return kERROR;
    219 
    220         if (!memcmp(data, "EVTE", 4))
    221         {
    222             // we found the end of list of events
    223             UnreadLastData();
    224             return kFALSE;
    225         }
    226 
    227         next = data;
    228     }
    229 
    230     *buffer = next;
    231     next += 7;
    232 
    233     return kTRUE;
    234 }
    235 
    236164///////////////////////////////////////////////////////////////////////////////
    237165///////////////////////////////////////////////////////////////////////////////
    238166
    239 // --------------------------------------------------------------------------
    240 //
    241 // Jumps from one top level object to the next until if finds the object with
    242 // correct type and correct id. The id is identifier of the raw corsika block,
    243 // for example "RUNH", RUNE", "EVTH"
    244 //
    245 // Returns kTRUE if the functions finds the type / id
    246 //         kFALSE if the functions does not find the type / id.
    247 //
    248 // After this call the file position pointer points just after the 4 bytes
    249 // of the id.
    250 //
    251 Int_t MCorsikaFormatEventIO::SeekNextBlock(const char * id, unsigned short type) const
    252 {
    253     cout << "Seek " << type << endl;
    254     while (1)
    255     {
    256         // we read - synchronisation marker
    257         //         - type / version field
    258         //         - identification field
    259         //         - length
    260         //         - unknown field
    261         //         - id (first 4 bytes of data field)
    262         int blockHeader[4];
    263         fIn->read((char*)blockHeader, 4 * sizeof(int));
    264 
    265         if (fIn->eof())
    266         {
    267             gLog << err << "MCorsikaFormatEventIO::SeekNextBlock: ERROR - Unexpected end-of-file seeking " << id << " (" << type << ")." << endl;
    268             return kERROR;
    269         }
    270 
    271         const unsigned short objType = blockHeader[1] & 0xFFFF;
    272         cout << "objType=" << objType << endl;
    273         if (type == objType)
    274         {
    275             if (!id)
    276                 break;
    277 
    278             char c[9];
    279             fIn->read(c, 8);
    280 
    281             if (memcmp(c+4, id, 4)==0)
    282                 break;
    283 
    284             c[8] = 0;
    285             gLog << err << "MCorsikaFormatEventIO::SeekNextBlock: ERROR - '" << c+4 << "' doesn't match expected one '" << id << "' in object " << type << endl;
    286             return kFALSE;
    287         }
    288 
    289         // we are looking for a event header, but found a runEnd.
    290         // This will happen at the end of a run. We stop here looking for
    291         // the next event header
    292         if (type == 1202 && objType == 1210)
    293             return kFALSE;
    294 
    295         // a unknown block, we jump to the next one
    296         const int length = blockHeader[3] & 0x3fffffff;
    297         fIn->seekg(length, ios::cur);
    298     }
    299 
    300 
     167Bool_t MCorsikaFormatEventIO::NextBlock(Bool_t  subBlock,
     168                                        Int_t & blockType,
     169                                        Int_t & blockVersion,
     170                                        Int_t & blockIdentifier,
     171                                        Int_t & blockLength) const
     172{
     173// we read - synchronisation markerif not subBlock
     174//         - type / version field
     175//         - identification field
     176//         - length
     177//         - unknown field
     178//         - id (first 4 bytes of data field)
     179
     180   if (subBlock)
     181      {
     182      int blockHeader[3];
     183      fIn->read((char*)blockHeader, 3 * sizeof(int));
     184
     185      if (fIn->eof())
     186         return kFALSE;
     187
     188
     189      blockType       = blockHeader[0] & 0xFFFF;
     190      blockVersion    = (blockHeader[0] & 0xFFF00000) >> 20;
     191      blockIdentifier = blockHeader[1];
     192      blockLength     = blockHeader[2] & 0x3FFFFFFF;
     193      }
     194   else
     195      {
     196       int blockHeader[4];
     197       fIn->read((char*)blockHeader, 4 * sizeof(int));
     198
     199       if (fIn->eof())
     200           return kFALSE;
     201
     202
     203       blockType       = blockHeader[1] & 0xFFFF;
     204       blockVersion    = (blockHeader[1] & 0xFFF00000) >> 20;
     205       blockIdentifier = blockHeader[2];
     206       blockLength     = blockHeader[3] & 0x3FFFFFFF;
     207
     208       if (blockType == 1200  || blockType == 1210 ||
     209           blockType == 1202  || blockType == 1209    )
     210           {
     211           // read the "old" corsika header plus additional 4 unknown byte
     212           char tmp[8];
     213           fIn->read(tmp, 8);
     214           blockLength -= 8;
     215           }
     216   
     217      }
    301218    return kTRUE;
    302219}
    … …  
    304221// --------------------------------------------------------------------------
    305222//
    306 void MCorsikaFormatEventIO::UnreadLastHeader() const
    307 {
    308     fIn->seekg( -6 * (Int_t)(sizeof(int)), ios::cur);
    309 }
    310 
    311 // --------------------------------------------------------------------------
    312 //
    313223Bool_t MCorsikaFormatEventIO::SeekEvtEnd()
    314224{
    315     if (fRunePos != streampos(0))
    316     {
    317         fIn->seekg(fRunePos, ios::beg);
    318         return kTRUE;
    319     }
    320 
    321     // it is the first time we are looking for the RUNE block
    322 
    323     // is the RUNE block at the very end of the file?
    324     const streampos currentPos = fIn->tellg();
    325 
     225
     226    // the RUNE block it at the very end of the file.
    326227    fIn->seekg(-32, ios::end);
    327228
    … …  
    334235    {
    335236        // this seams to be a RUNE (1210)  block
    336         fIn->seekg( -4 * (Int_t)(sizeof(int)), ios::cur);
    337         fRunePos = fIn->tellg();
     237        fIn->seekg( 8, ios::cur);
    338238        return kTRUE;
    339239    }
    340240
    341     // we do not find a RUNE block at the end of the file
    342     // we have to search in the file
    343     fIn->seekg(currentPos, ios::beg);
    344     if (SeekNextBlock("RUNE", 1210)!=kTRUE)
    345         return kFALSE;
    346 
    347     UnreadLastHeader();
    348     fRunePos = fIn->tellg();
    349 
    350     return kTRUE;
    351 }
    352 
    353 // --------------------------------------------------------------------------
    354 //                                                                           
    355 // Returns one event (7 Float values) after the other until the EventEnd     
    356 // block is found.                                                           
    357 // If a read error occurred, the readError is set to kTRUE and the function 
    358 // returns kFALSE;                                                           
    359 //
    360 Int_t MCorsikaFormatEventIO::GetNextEvent(Float_t **buffer, Int_t telescope)
    361 {
    362     static Float_t *data = NULL;
    363 
    364     static int lengthPhotonData = 0;
    365     static int lengthEvent      = 0;
    366 
    367     if (lengthPhotonData>0)
    368     {
    369         cout << "Return Bunch l2=" << lengthPhotonData << endl;
    370 
    371         lengthPhotonData -= 8 * sizeof(Float_t);
    372         *buffer += 8; // Return the pointer to the start of the current photon data
    373         return kTRUE;
    374     }
    375 
    376     if (data)
    377     {
    378         delete [] data;
    379         data = NULL;
    380         cout << "Return end-of-tel LE=" << lengthEvent << endl;
    381         return kFALSE;
    382     }
    383 
    384     if (lengthEvent==0)
    385     {
    386         cout << "Searching 1204... " << flush;
    387         // The length of the block is stored and we use it to determine
    388         // when the next top level block is expected
    389         const Int_t rc = NextEventObject(lengthEvent);
    390         if (rc==kERROR)
    391             return kERROR;
    392         if (!rc)
    393         {
    394             cout << "skip EVTE." << endl;
    395             return kFALSE;
    396         }
    397 
    398         cout << "found new array." << endl;
    399     }
    400 
    401     cout << "Searching 1205..." << flush;
    402 
    403     // Look for the next event photon bunch (object type 1205)
    404     const Int_t tel = NextPhotonObject(lengthPhotonData);
    405     if (tel<0)
    406         return kERROR;
    407 
    408     lengthEvent -= lengthPhotonData+12; // Length of data + Length of header
    409 
    410     lengthPhotonData -= 12;     // Length of data before the photon bunches
    411     fIn->seekg(12, ios::cur);   // Skip this data
    412 
    413     cout << " found (tel=" << tel << ",LEN=" << lengthEvent << ",len=" << lengthPhotonData << ")... " << flush;
    414 
    415     if (lengthPhotonData==0)
    416     {
    417         cout << "empty!" << endl;
    418         return kFALSE;
    419     }
    420 
    421     const UInt_t cnt = lengthPhotonData / sizeof(Float_t);
    422     // Read next object (photon bunch)
    423     data = new Float_t[cnt];
    424     if (!ReadData(cnt, data, 0))
    425     {
    426         delete [] data;
    427         data = NULL;
    428         return kERROR;
    429     }
    430 
    431     cout << "return!" << endl;
    432 
    433     lengthPhotonData -= 8 * sizeof(Float_t);
    434     *buffer = data; // Return the pointer to the start of the current photon data
    435 
    436     return kTRUE;
    437 }
    438 
    439 // --------------------------------------------------------------------------
    440 //                                                                           
    441 // Looks for the next Block with type 1204 and return kTRUE.                 
    442 // The function also stops moving forward in the file, if it finds a         
    443 // EventEnd block (1209). In this case kFALSE is returned
    444 //
    445 Int_t MCorsikaFormatEventIO::NextEventObject(Int_t &length) const
    446 {
    447     while (1)
    448     {
    449         // we read - synchronisation marker
    450         //         - type / version field
    451         //         - identification field
    452         //         - length
    453         UInt_t blockHeader[4];
    454         fIn->read((char*)blockHeader, 4 * sizeof(Int_t));
    455 
    456         if (fIn->eof())
    457         {
    458             gLog << err << "MCorsikaFormatEventIO::NextEventObject: ERROR - Unexpected end-of-file." << endl;
    459             return kERROR;
    460         }
    461 
    462         // For speed reasons we skip the check of the synchronization marker
    463 
    464         // Decode the object type
    465         const unsigned short objType = blockHeader[1] & 0xFFFF;
    466 
    467         cout << "objType=" << objType << endl;
    468 
    469         // Decode length of block
    470         length = blockHeader[3] & 0x3fffffff;
    471 
    472         // Check if we found the next array (reuse / impact parameter)
    473         // blockHeader[2] == array number (reuse)
    474         if (objType==1204)
    475             return kTRUE;
    476 
    477         // we found an eventEnd block, reset file pointer
    478         if (objType==1209)
    479             break;
    480 
    481         // a unknown block, we jump to the next one
    482         fIn->seekg(length, ios::cur);
    483     }
    484 
    485     fIn->seekg( -4 * (Int_t)(sizeof(Int_t)), ios::cur);
    486     length = 0;
    487 
    488241    return kFALSE;
    489242}
    490243
    491 // --------------------------------------------------------------------------
    492 //
    493 Int_t MCorsikaFormatEventIO::NextPhotonObject(Int_t &length) const
    494 {
    495     UInt_t blockHeader[3];
    496 
    497     // we read - synchronisation marker
    498     //         - type / version field
    499     //         - identification field
    500     //         - length
    501     fIn->read((char*)blockHeader, 3 * sizeof(UInt_t));
    502 
    503     if (fIn->eof())
    504     {
    505         gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected end-of-file.." << endl;
    506         return -1;
    507     }
    508 
    509     const unsigned short objType = blockHeader[0] & 0xFFFF;
    510     if (objType != 1205)
    511     {
    512         gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected EventIO object type " << objType << " (expected 1205)." << endl;
    513         return -1;
    514     }
    515 
    516     const unsigned short version = (blockHeader[0] >> 20) & 0x0FFF;
    517     if (version != 0)
    518     {
    519         gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected EventIO object version " << version << " (expected 0 for object type 1205)." << endl;
    520         return -1;
    521     }
    522 
    523     length = blockHeader[2] & 0x3fffffff;
    524 
    525     // blockHeader[1] == 1000 x array number (reuse) + telescope number
    526     return blockHeader[1] % 1000;
    527 }

  • trunk/Mars/mcorsika/MCorsikaFormat.h

    r9949 r10060  
    1717   std::istream  *fIn;
    1818
    19    std::streampos fPrevPos; // file position before previous read
    20    std::streampos fPos;
    21 
    2219public:
    2320    static const unsigned int kSyncMarker;
    … …  
    2724   virtual ~MCorsikaFormat();
    2825
    29    virtual Int_t  SeekNextBlock(const char * id, unsigned short type) const = 0;
    30    virtual void   UnreadLastHeader() const = 0;
     26   virtual Bool_t NextBlock(Bool_t  subBlock, Int_t & blockType, Int_t & blockVersion,
     27                            Int_t & blockIdentifier, Int_t & blockLength) const = 0;
    3128
    32    virtual Bool_t ReadData(Int_t numValues, Float_t * buffer,
    33                                 Int_t minSeekValues = 272);
    34    virtual void   UnreadLastData() const;
     29           void   Seek(std::streampos offset)    {fIn->seekg(offset, std::ios::cur);}
    3530
    3631   virtual Bool_t SeekEvtEnd() = 0;
    37    virtual void   StorePos();
    38    virtual void   ResetPos() const;
    39    virtual void   Rewind() const;
    4032
    41    virtual Int_t  GetNextEvent(Float_t **buffer, Int_t tel=0) = 0;
    4233   virtual Bool_t IsEventioFormat() const = 0;
    4334
    4435   virtual Bool_t Eof() const;
    4536
    46    std::streampos GetCurrPos() const;
    47 
    48    void Read(void *ptr, Int_t i=4) const;
     37           Bool_t Read(void *ptr, Int_t i) const;
    4938
    5039   static MCorsikaFormat *CorsikaFormatFactory(const char *fileName);
    … …  
    6049        : MCorsikaFormat(in) {}
    6150
    62    Int_t  SeekNextBlock(const char * id, unsigned short type) const;
    63    void   UnreadLastHeader() const;
     51   Bool_t NextBlock(Bool_t  subBlock, Int_t & blockType, Int_t & blockVersion,
     52                    Int_t & blockIdentifier, Int_t & blockLength) const;
    6453
    6554   Bool_t SeekEvtEnd();
    6655
    67    Int_t  GetNextEvent(Float_t **buffer, Int_t);
    6856   Bool_t IsEventioFormat() const {return kFALSE;}
    6957};
    … …  
    7260class MCorsikaFormatEventIO : public MCorsikaFormat
    7361{
    74 private:
    75     std::streampos fRunePos; // file position of the RUNE block
    7662
    7763public:
    7864    MCorsikaFormatEventIO(std::istream *in)
    79         : MCorsikaFormat(in) {fRunePos = std::streampos(0);}
     65        : MCorsikaFormat(in) {}
    8066
    81     Int_t  SeekNextBlock(const char *id, unsigned short type) const;
    82     void   UnreadLastHeader() const;
     67
     68    Bool_t NextBlock(Bool_t  subBlock, Int_t & blockType, Int_t & blockVersion,
     69                     Int_t & blockIdentifier, Int_t & blockLength) const;
    8370
    8471    Bool_t SeekEvtEnd();
    8572
    86     Int_t  GetNextEvent(Float_t **buffer, Int_t tel);
    8773    Bool_t IsEventioFormat() const { return kTRUE; }
    8874
    89 private:
    90     Int_t  NextEventObject(Int_t &length) const;
    91     Int_t  NextPhotonObject(Int_t &length) const;
    9275};
    9376

  • trunk/Mars/mcorsika/MCorsikaRead.cc

    r9943 r10060  
    1717!
    1818!   Author(s): Thomas Bretz  11/2008 <mailto:tbretz@astro.uni-wuerzburg.de>
     19               Reiner Rohlfs 11/2010 <mailto:Reiner.Rohlfs@unige.ch>
    1920!
    2021!   Copyright: Software Development, 2000-2008
    … …  
    5859
    5960using namespace std;
     61
    6062
    6163/*  ----------- please don't delete and don't care about (Thomas) ------------
    … …  
    9698//
    9799MCorsikaRead::MCorsikaRead(const char *fname, const char *name, const char *title)
    98     : fRunHeader(0), fEvtHeader(0), fEvent(0), fTelescopeIdx(-1), fForceMode(kFALSE),
    99     fFileNames(0), fNumFile(0), fNumEvents(0), fNumTotalEvents(0),
    100     /*fIn(0),*/  fInFormat(0), fParList(0), fNumTelescopes(1)
     100    : fRunHeader(0), fEvtHeader(0), fEvent(0), fTelescopeIdx(-1), fArrayIdx(-1),
     101    fForceMode(kFALSE), fFileNames(0), fNumFile(0), fNumEvents(0),
     102    fNumTotalEvents(0), fInFormat(0), fParList(0), fNumTelescopes(1), fReadState(0)
    101103{
    102104    fName  = name  ? name  : "MRead";
    … …  
    117119{
    118120    delete fFileNames;
    119 //    if (fIn)
    120 //        delete fIn;
    121121    if (fInFormat)
    122122        delete fInFormat;
    … …  
    161161
    162162}
    163 
    164 Bool_t MCorsikaRead::ReadEvtEnd()
    165 {
    166     if (!fInFormat->SeekEvtEnd())
    167     {
    168         *fLog << (fForceMode?warn:err) << "Error: RUNE section not found in file." << endl;
    169         if (!fForceMode)
    170             return kFALSE;
    171     }
    172 
    173     if (!fRunHeader->ReadEvtEnd(fInFormat))
    174     {
    175         *fLog << (fForceMode?warn:err) << "Error: Reading RUNE section failed." << endl;
    176         if (!fForceMode)
    177             return kFALSE;
    178     }
    179 
    180     return kTRUE;
    181 }
    182 
    183163// --------------------------------------------------------------------------
    184164//
    … …  
    191171    // open the input stream and check if it is really open (file exists?)
    192172    //
    193 /*    if (fIn)
    194         delete fIn;
    195     fIn = NULL;
    196 */
    197173    if (fInFormat)
    198174       delete fInFormat;
    … …  
    231207
    232208    fNumFile++;
    233 
    234     //
    235     // Read RUN HEADER (see specification) from input stream
    236     //
    237     if (!fRunHeader->ReadEvt(fInFormat))
    238         return kERROR;
    239 //    if (!fEvtHeader->ReadRunHeader(*fIn, *fRunHeader))
    240 //        return kERROR;
    241 
    242     fNumTelescopes = 1;
    243     if (fInFormat->IsEventioFormat())
    244     {
    245         fInFormat->StorePos();
    246         fInFormat->Rewind();
    247 
    248         if (!fInFormat->SeekNextBlock(0, 1201))
    249         {
    250             *fLog << err << "ERROR - Object 1201 not found in file." << endl;
    251             return kERROR;
    252         }
    253 
    254         fInFormat->Read(&fNumTelescopes);
    255 
    256         if (fTelescopeIdx>=fNumTelescopes)
    257         {
    258             *fLog << err << "ERROR - Requested telescope index " << fTelescopeIdx << " exceeds number of telescopes " << fNumTelescopes << " in file." << endl;
    259             return kERROR;
    260         }
    261 
    262         fTelescopeX.Set(fNumTelescopes);
    263         fTelescopeY.Set(fNumTelescopes);
    264         fTelescopeZ.Set(fNumTelescopes);
    265         fTelescopeR.Set(fNumTelescopes);
    266 
    267         fInFormat->Read(fTelescopeX.GetArray(), 4*fNumTelescopes);
    268         fInFormat->Read(fTelescopeY.GetArray(), 4*fNumTelescopes);
    269         fInFormat->Read(fTelescopeZ.GetArray(), 4*fNumTelescopes);
    270         fInFormat->Read(fTelescopeR.GetArray(), 4*fNumTelescopes);
    271 
    272         *fLog << all;
    273         *fLog << "Number of telescopes: " << fNumTelescopes << " (using ";
    274         if (fTelescopeIdx>=0)
    275             *fLog << "telescope " << fTelescopeIdx;
    276         else
    277             *fLog << "all telescopes";
    278         *fLog << ")" << endl;
    279 
    280         const Int_t lo = fTelescopeIdx<0 ? 0              : fTelescopeIdx;
    281         const Int_t up = fTelescopeIdx<0 ? fNumTelescopes : fTelescopeIdx+1;
    282 
    283         for (int i=lo; i<up; i++)
    284         {
    285             *fLog << inf << "Telescope #" << setw(4) << i << " [X/Y/Z (R)]: ";
    286             *fLog << setw(7) << fTelescopeX[i] << "/";
    287             *fLog << setw(7) << fTelescopeY[i] << "/";
    288             *fLog << setw(7) << fTelescopeZ[i] << "  (R=" << setw(4) << fTelescopeR[i] << ")" << endl;
    289         }
    290 
    291         fNumTelescope = 0;
    292 
    293         fInFormat->ResetPos();
    294     }
    295 
    296 
    297     fInFormat->StorePos();
    298     if (!ReadEvtEnd())
    299         return kERROR;
    300     fInFormat->ResetPos();
    301 
    302 
    303     fNumEvents += fRunHeader->GetNumEvents();
    304     fRunHeader->SetReadyToSave();
    305 
    306     //
    307     // Print Run Header
    308     //  We print it after the first event was read because
    309     //  we still miss information which is stored in the event header?!?
    310     if (print)
    311         fRunHeader->Print();
    312209
    313210    if (!fParList)
    … …  
    368265            break;
    369266        case kTRUE:
    370             if (!ReadEvtEnd())
     267
     268            // read run header(1200), telescope position(1201) and
     269            // first event header(1202)
     270            Bool_t status = kTRUE;
     271            Int_t blockType, blockVersion, blockIdentifier, blockLength;
     272            while (status &&
     273                   fInFormat->NextBlock(fReadState == 4, blockType, blockVersion,
     274                              blockIdentifier, blockLength))
     275               {
     276               if (blockType == 1200)
     277                  status = fRunHeader->ReadEvt(fInFormat);
     278
     279               else if(blockType == 1201)
     280                   status = ReadTelescopePosition();
     281
     282               else if (blockType == 1202)
     283                  {
     284                  Float_t buffer[272];
     285                  status = fInFormat->Read(buffer, 272 * sizeof(Float_t));
     286                  status = fRunHeader->ReadEventHeader(buffer);
     287                  break;
     288                  }
     289               else
     290                  fInFormat->Seek(blockLength);
     291               }
     292                 
     293            if (status != kTRUE)
     294               return status;
     295
     296            if (!fInFormat->SeekEvtEnd())
    371297            {
    372                 rc = kFALSE;
    373                 break;
     298               *fLog << (fForceMode?warn:err) << "Error: RUNE section not found in file." << endl;
     299               if (!fForceMode)
     300                  return fForceMode ? kTRUE : kFALSE;
     301            }
     302
     303            if (!fRunHeader->ReadEvtEnd(fInFormat, kTRUE))
     304            {
     305               *fLog << (fForceMode?warn:err) << "Error: Reading RUNE section failed." << endl;
     306               if (!fForceMode)
     307                  return kFALSE;
    374308            }
    375309
    … …  
    380314        break;
    381315    }
    382 /*
    383     if (fIn)
    384         delete fIn;
    385     fIn = NULL;
    386 */
     316
    387317    return rc;
     318}
     319
     320Int_t MCorsikaRead::ReadTelescopePosition()
     321{
     322   if (!fInFormat->Read(&fNumTelescopes, 4)) return kERROR;
     323
     324   if (fTelescopeIdx>=fNumTelescopes)
     325   {
     326      *fLog << err << "ERROR - Requested telescope index " << fTelescopeIdx <<
     327            " exceeds number of telescopes " << fNumTelescopes << " in file." << endl;
     328      return kERROR;
     329   }
     330
     331   fTelescopeX.Set(fNumTelescopes);
     332   fTelescopeY.Set(fNumTelescopes);
     333   fTelescopeZ.Set(fNumTelescopes);
     334   fTelescopeR.Set(fNumTelescopes);
     335
     336   if (!fInFormat->Read(fTelescopeX.GetArray(), 4*fNumTelescopes)) return kERROR;
     337   if (!fInFormat->Read(fTelescopeY.GetArray(), 4*fNumTelescopes)) return kERROR;
     338   if (!fInFormat->Read(fTelescopeZ.GetArray(), 4*fNumTelescopes)) return kERROR;
     339   if (!fInFormat->Read(fTelescopeR.GetArray(), 4*fNumTelescopes)) return kERROR;
     340
     341   *fLog << all;
     342   *fLog << "Number of telescopes: " << fNumTelescopes << " (using ";
     343   if (fTelescopeIdx>=0)
     344      *fLog << "telescope " << fTelescopeIdx;
     345   else
     346      *fLog << "all telescopes";
     347   *fLog << ")" << endl;
     348
     349   const Int_t lo = fTelescopeIdx<0 ? 0              : fTelescopeIdx;
     350   const Int_t up = fTelescopeIdx<0 ? fNumTelescopes : fTelescopeIdx+1;
     351
     352   for (int i=lo; i<up; i++)
     353   {
     354      *fLog << inf << "Telescope #" << setw(4) << i << " [X/Y/Z (R)]: ";
     355      *fLog << setw(7) << fTelescopeX[i] << "/";
     356      *fLog << setw(7) << fTelescopeY[i] << "/";
     357      *fLog << setw(7) << fTelescopeZ[i] << "  (R=" << setw(4) << fTelescopeR[i] << ")" << endl;
     358   }
     359
     360   fNumTelescope = 0;
     361
     362   return kTRUE;
    388363}
    389364
    … …  
    412387    //
    413388    fParList = 0;
    414 /*
    415     if (!OpenStream())
    416         return kFALSE;
    417 */
     389
    418390    //
    419391    //  check if all necessary containers exist in the Parameter list.
    … …  
    445417}
    446418
    447 // --------------------------------------------------------------------------
    448 //
    449 // Read a single event from the stream
    450 //
    451 Int_t MCorsikaRead::ReadEvent()
    452 {
    453     if (fInFormat->IsEventioFormat())
    454     {
    455         if (fNumTelescope==0)
    456         {
    457             const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat);
    458             if (rc1!=kTRUE)
    459                 return rc1;
    460 
    461             if (fEvtHeader->GetNumReuse()==fRunHeader->GetNumReuse()-1)
    462                 fEvtHeader->ResetNumReuse();
    463             fEvtHeader->IncNumReuse();
    464 
    465             cout << "===== Array " << fEvtHeader->GetNumReuse() << " =====" << endl;
    466         }
    467 
    468         const Int_t rc2 = fEvent->ReadCorsikaEvt(fInFormat, fTelescopeIdx);
    469         if (rc2!=kTRUE)
    470             return rc2;
    471 
    472         fEvtHeader->InitXY();
    473 
    474         const Double_t x = fTelescopeX[fNumTelescope];
    475         const Double_t y = fTelescopeY[fNumTelescope];
    476 
    477         fEvtHeader->AddXY(y, -x);
    478         fEvent->AddXY(fEvtHeader->GetX(), fEvtHeader->GetY());
    479         fEvent->SimWavelength(fRunHeader->GetWavelengthMin(), fRunHeader->GetWavelengthMax());
    480 
    481         fNumTelescope++;
    482         fNumTelescope%=fNumTelescopes;
    483 
    484         return kTRUE;
    485     }
    486     else
    487     {
    488         // Store the position to re-read a single event several times
    489         if (fEvtHeader->GetNumReuse()<fRunHeader->GetNumReuse()-1)
    490             fInFormat->ResetPos();
    491         else
    492         {
    493             fInFormat->StorePos();
    494             fEvtHeader->ResetNumReuse();
    495         }
    496 
    497         // Read the event header
    498         const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat);
    499         if (rc1!=kTRUE)
    500             return rc1;
    501 
    502         fEvtHeader->IncNumReuse();
    503         fEvtHeader->InitXY();
    504 
    505         // In the case of corsika raw data (MMCS only) we have too loop over one event
    506         // several times to read all impact parameters (reusage of events)
    507         // In case of the eventio format we have to decide, the data of which telescope
    508         // we want to extract
    509 
    510         // Read the photons corresponding to the i-th core location
    511         //  EventIO: Number of telescope
    512         //  Raw:     Number of re-use
    513         const Int_t rc2 = fEvent->ReadCorsikaEvt(fInFormat, fEvtHeader->GetNumReuse()+1);
    514         if (rc2!=kTRUE)
    515             return rc2;
    516 
    517         // read event end
    518         return fEvtHeader->ReadEvtEnd(fInFormat);
    519     }
    520 }
    521419
    522420// --------------------------------------------------------------------------
    … …  
    530428Int_t MCorsikaRead::Process()
    531429{
    532     while (1)
    533     {
    534         if (fInFormat)
    535         {
    536             // Read a single event from file
    537             const Int_t rc = ReadEvent();
    538 
    539             // kFALSE means: end of file (try next one)
    540             if (rc!=kFALSE)
    541                 return rc;
    542 
    543 
    544             fInFormat->UnreadLastHeader();
    545             if (!fRunHeader->ReadEvtEnd(fInFormat))
    546                 if (!fForceMode)
    547                     return kERROR;
    548         }
    549 
    550         //
    551         // If an event could not be read from file try to open new file
    552         //
    553         const Int_t rc = OpenNextFile();
    554         if (rc!=kTRUE)
    555             return rc;
    556     }
    557     return kTRUE;
     430
     431   if (fReadState == 11)
     432      {
     433      // we are currently saving the events of the raw format in the root file
     434      if (fEvtHeader->GetNumReuse() == fRunHeader->GetNumReuse())
     435         {
     436         // all data are saved
     437         fRawEvemtBuffer.resize(0);
     438         fReadState = 3;
     439         }
     440      else
     441         {
     442         fEvtHeader->InitXY();
     443         Int_t rc = fEvent->ReadCorsikaEvt(&fRawEvemtBuffer[0],
     444                                           fRawEvemtBuffer.size() / 7,
     445                                           fEvtHeader->GetNumReuse()+1);
     446         fEvtHeader->IncNumReuse();
     447         return rc;
     448         }
     449      }
     450
     451   while (1)
     452   {
     453      if (fInFormat)
     454      {
     455         Int_t blockType, blockVersion, blockIdentifier, blockLength;
     456
     457         while (fInFormat->NextBlock(fReadState == 4, blockType, blockVersion,
     458                                     blockIdentifier, blockLength))
     459         {
     460//            gLog << "Next block:  type=" << blockType << " version=" << blockVersion;
     461//            gLog << " identifier=" << blockIdentifier << " length=" << blockLength << endl;
     462
     463
     464            if (fReadState == 4)
     465               {
     466               fTopBlockLength -= blockLength + 12;
     467               if (fTopBlockLength <= 0)
     468                  // all data of a top block are read, go back to normal state
     469                  fReadState = 3;
     470               }
     471
     472            Int_t status = kTRUE;
     473            switch (blockType)
     474               {
     475               case 1200:       // the run header
     476                  status = fRunHeader->ReadEvt(fInFormat);
     477                  fReadState = 1;  // RUNH is read
     478                  break;
     479
     480               case 1201:       // telescope position
     481                  status = ReadTelescopePosition();
     482                  break;
     483
     484               case 1202:  // the event header
     485                  Float_t buffer[272];
     486                  if (!fInFormat->Read(buffer, 272 * sizeof(Float_t)))
     487                     return kFALSE;
     488
     489                  if (fReadState == 1)  // first event after RUN header
     490                     {
     491                     fRunHeader->ReadEventHeader(buffer);
     492                     fRunHeader->Print();
     493                     }
     494
     495                  status = fEvtHeader->ReadEvt(buffer);
     496                  if (fArrayIdx >= (Int_t)fEvtHeader->GetTotReuse())
     497                     {
     498                     *fLog << err << "ERROR - Requested array index " << fArrayIdx <<
     499                           " exceeds number of arrays " << fEvtHeader->GetTotReuse() <<
     500                           " in file." << endl;
     501                     return kERROR;
     502                     }
     503
     504
     505                  fReadState = 2;
     506                  break;
     507
     508               case 1204:
     509                  if (fArrayIdx < 0 || fArrayIdx == blockIdentifier)
     510                     {
     511                     fReadState = 4;
     512                     fTopBlockLength = blockLength;
     513                     }
     514                  else
     515                     // skip this array of telescopes
     516                     fInFormat->Seek(blockLength);
     517
     518                  break;
     519               
     520               case 1205:
     521                  {
     522                  Int_t telIdx   = blockIdentifier % 1000;
     523                  if (blockVersion == 0                               &&
     524                      (fTelescopeIdx < 0 || fTelescopeIdx ==  telIdx) &&
     525                       blockLength > 12)
     526                     {
     527                     status = fEvent->ReadEventIoEvt(fInFormat);
     528
     529                     Int_t arrayIdx = blockIdentifier / 1000;
     530                     Float_t xArrOff, yArrOff;
     531                     fEvtHeader->GetArrayOffset(arrayIdx, xArrOff, yArrOff);
     532                     fEvtHeader->SetTelescopeOffset(arrayIdx,
     533                                                    xArrOff + fTelescopeY[telIdx],
     534                                                    yArrOff - fTelescopeX[telIdx] );
     535                     fEvent->AddXY(xArrOff + fTelescopeY[telIdx],
     536                                   yArrOff - fTelescopeX[telIdx]);
     537                     fEvent->SimWavelength(fRunHeader->GetWavelengthMin(),
     538                                           fRunHeader->GetWavelengthMax());
     539
     540                     if (status == kTRUE)
     541                        // end of reading for one telescope in one array ==
     542                        // end of this Process - step
     543                        return status;
     544                     }
     545                  else
     546                     // skip this telescope
     547                     fInFormat->Seek(blockLength);
     548                  }
     549                  break;
     550
     551               case 1209:  // the event end
     552                  status = fEvtHeader->ReadEvtEnd(fInFormat);
     553                  if (fReadState == 10)
     554                     {
     555                     // all particles of this event are read, now save them
     556                     fReadState = 11;
     557                     fEvtHeader->ResetNumReuse();
     558
     559                     fEvtHeader->InitXY();
     560                     Int_t rc = fEvent->ReadCorsikaEvt(&fRawEvemtBuffer[0],
     561                                                       fRawEvemtBuffer.size() / 7,
     562                                                       fEvtHeader->GetNumReuse()+1);
     563                     fEvtHeader->IncNumReuse();
     564                     return rc;
     565                     }
     566                  else
     567                     fReadState = 3;  // event closed, run still open
     568                  break;
     569
     570               case 1210:  // the run end
     571                  status = fRunHeader->ReadEvtEnd(fInFormat, kTRUE);
     572                  fNumEvents += fRunHeader->GetNumEvents();
     573                  fRunHeader->SetReadyToSave();
     574                  fReadState = 5;           // back to  starting point
     575                  return status;
     576                  break;
     577
     578               case 1105:  // event block of raw format
     579                  if (fReadState == 2 || fReadState == 10)
     580                     {
     581                     Int_t oldSize = fRawEvemtBuffer.size();
     582                     fRawEvemtBuffer.resize(oldSize + blockLength / sizeof(Float_t));
     583                     status = fInFormat->Read(&fRawEvemtBuffer[oldSize], blockLength);
     584                     fReadState = 10;
     585                     }
     586                  else
     587                     fInFormat->Seek(blockLength);
     588                  break;
     589
     590               default:
     591                  // unknown block, we skip it
     592                  fInFormat->Seek(blockLength);
     593
     594               }
     595
     596            if (status != kTRUE)
     597                return status;
     598           
     599         }
     600
     601      }
     602
     603      //
     604      // If an event could not be read from file try to open new file
     605      //
     606      const Int_t rc = OpenNextFile();
     607      if (rc!=kTRUE)
     608          return rc;
     609   }
     610   return kTRUE;
    558611}
    559612
    … …  
    566619Int_t MCorsikaRead::PostProcess()
    567620{
     621
    568622    const UInt_t n = fNumEvents*fRunHeader->GetNumReuse()*(fTelescopeIdx<0?fNumTelescopes:1);
    569623
    … …  
    575629
    576630    *fLog << warn << dec;
    577     *fLog << "Warning - number of read events (" << GetNumExecutions()-1;
     631    *fLog << "Warning - number of read events (" << GetNumExecutions()-2;
    578632    *fLog << ") doesn't match number in run header(s) (";
    579633    *fLog << fRunHeader->GetNumReuse() << "*" << fNumEvents << "=" << n << ")." << endl;
    … …  
    582636}
    583637
    584 // --------------------------------------------------------------------------
    585 //
    586 //    Telescope: 1
    587 //
    588 // In the case of eventio-format select which telescope to extract from the
    589 // file. -1 will extract all telescopes
    590 //
    591 Int_t MCorsikaRead::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
    592 {
    593     Bool_t rc = kFALSE;
    594     if (IsEnvDefined(env, prefix, "Telescope", print))
    595     {
    596         rc = kTRUE;
    597         fTelescopeIdx = GetEnvValue(env, prefix, "Telescope", fTelescopeIdx);
    598     }
    599 
    600     return rc;
    601 }

  • trunk/Mars/mcorsika/MCorsikaRead.h

    r9943 r10060  
    11#ifndef MARS_MCorsikaRead
    22#define MARS_MCorsikaRead
     3
     4#include <vector>
    35
    46#ifndef MARS_MRead
    … …  
    2729
    2830    Int_t           fTelescopeIdx;  // Index of telescope to be extracted
     31    Int_t           fArrayIdx;      // Index of telescope-array to be extracted
    2932    Bool_t          fForceMode;     // Force mode skipping defect RUNE
    3033
    … …  
    3437    UInt_t    fNumTotalEvents; //! total number of events in all files
    3538
    36 //    ifstream       *fIn;       //! input stream (file to read from)
    3739    MCorsikaFormat *fInFormat; //! access to input corsika data
    3840
    … …  
    4648    TArrayF fTelescopeR;       //! Radii of spheres around tel. (unit: cm)
    4749
     50    Int_t   fReadState;        //! 0: nothing read yet
     51                               //  1: RUNH is already read
     52                               //  2: EVTH is already read
     53                               //  3: EVTE is already read, run still open
     54                               //  4: inside of a top level block (1205)
     55                               //  5: RUNE is already read
     56                               // 10: raw events are currently read
     57                               // 11: raw events are currently saved
     58    Int_t   fTopBlockLength;   // remaining length of the current top-level block 1204
     59
     60    std::vector<Float_t>  fRawEvemtBuffer;     //! buffer of raw event data
    4861    //UInt_t    fInterleave;
    4962    //Bool_t    fForce;
    5063
    51 //    virtual Bool_t OpenStream() { return kTRUE; }
    5264
    53     Bool_t ReadEvtEnd();
    5465    Int_t  OpenNextFile(Bool_t print=kTRUE);
    5566    Bool_t CalcNumTotalEvents();
    56     Int_t  ReadEvent();
     67    Int_t  ReadTelescopePosition();
    5768
    5869    // MTask
    … …  
    6172    Int_t PostProcess();
    6273
    63     // MParContainer
    64     Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print);
    6574
    6675public:
    … …  
    6978
    7079    void SetForceMode(Bool_t b=kTRUE) { fForceMode=b; }
    71     void SetTelescopeIdx(Int_t idx=-1) { fTelescopeIdx = idx; }
     80    void SetTelescopeIdx(Int_t idx)   { fTelescopeIdx = idx; }
     81    void SetArrayIdx(Int_t idx)       { fArrayIdx = idx; }
    7282
    7383    TString GetFullFileName() const;

  • trunk/Mars/mcorsika/MCorsikaRunHeader.cc

    r9949 r10060  
    8585Bool_t MCorsikaRunHeader::ReadEvt(MCorsikaFormat * fInFormat)
    8686{
    87     if (fInFormat->SeekNextBlock("RUNH", 1200)!=kTRUE)
    88         return kFALSE;
    89 
    9087    Float_t f[272];
    91     if (!fInFormat->ReadData(272, f))
     88    if (!fInFormat->Read(f, 272 * sizeof(Float_t)))
    9289        return kFALSE;
    9390
    … …  
    137134    memcpy(fAtmosphericCoeffB, f+258, 5*4);
    138135    memcpy(fAtmosphericCoeffC, f+263, 5*4);
     136
     137    return kTRUE;
     138}
     139
     140// --------------------------------------------------------------------------
     141//
     142// Read in one event header. It is called for the first event header after 
     143// a run header                                                             
     144//
     145Bool_t MCorsikaRunHeader::ReadEventHeader(Float_t * g)
     146{
    139147
    140148    // -------------------- Read first event header -------------------
    … …  
    158166    // f[145] Muon multiple scattering flag
    159167
    160     if (fInFormat->SeekNextBlock("EVTH", 1202)!=kTRUE)
    161         return kFALSE;
    162 
    163     Float_t g[273];
    164     if (!fInFormat->ReadData(272, g))
    165         return kFALSE;
    166 
    167     fInFormat->UnreadLastData();
    168     fInFormat->UnreadLastHeader();
    169168
    170169    fNumReuse = TMath::Nint(g[96]);  // Number i of uses of each cherenkov event
    171     if (fNumReuse!=1)
    172     {
    173         *fLog << err  << "ERROR   - Currently only one impact parameter per event is supported." << endl;
    174         *fLog << warn << "WARNING - This error is replaced by a warning." << endl;
    175     }
    176170
    177171    fParticleID = TMath::Nint(g[1]);
    … …  
    219213}
    220214
    221 Bool_t MCorsikaRunHeader::ReadEvtEnd(MCorsikaFormat * fInFormat)
    222 {
    223     if (fInFormat->SeekNextBlock("RUNE", 1210)!=kTRUE)
     215Bool_t MCorsikaRunHeader::ReadEvtEnd(MCorsikaFormat * fInFormat, Bool_t runNumberVerify)
     216{
     217    Float_t f[2];
     218    if (!fInFormat->Read(f, 2 * sizeof(Float_t)))
    224219        return kFALSE;
    225220
    226     Float_t f[2];
    227     if (!fInFormat->ReadData(2, f))
    228         return kFALSE;
    229 
    230     const UInt_t runnum = TMath::Nint(f[0]);
    231     if (runnum!=fRunNumber)
    232     {
    233         *fLog << err << "ERROR - Mismatch in stream: Run number in RUNE (";
    234         *fLog << runnum << ") doesn't match RUNH (" << fRunNumber << ")." << endl;
    235         return kFALSE;
    236     }
     221    if (runNumberVerify)
     222      {
     223       const UInt_t runnum = TMath::Nint(f[0]);
     224       if (runnum!=fRunNumber)
     225       {
     226           *fLog << err << "ERROR - Mismatch in stream: Run number in RUNE (";
     227           *fLog << runnum << ") doesn't match RUNH (" << fRunNumber << ")." << endl;
     228           return kFALSE;
     229       }
     230      }
    237231
    238232    fNumEvents = TMath::Nint(f[1]);
    … …  
    345339    *fLog << endl;
    346340
    347 }
    348 
     341
     342}
     343

  • trunk/Mars/mcorsika/MCorsikaRunHeader.h

    r9937 r10060  
    124124    // I/O
    125125    Bool_t ReadEvt(MCorsikaFormat * fInFormat);
    126     Bool_t ReadEvtEnd(MCorsikaFormat * fInFormat);
     126    Bool_t ReadEventHeader(Float_t * g);
     127    Bool_t ReadEvtEnd(MCorsikaFormat * fInFormat, Bool_t runNumberVerify);
    127128    Bool_t SeekEvtEnd(istream &fin);
    128129

  • trunk/Mars/msim/MPhotonData.cc

    r9941 r10060  
    222222{
    223223    const UInt_t n = TMath::Nint(f[0]);
     224
    224225    if (n==0)
    225226        // FIXME: Do we need to decode the rest anyway?
    … …  
    227228
    228229    // Check reuse
    229     const Int_t reuse = (n/1000)%100; // Force this to be 1!
    230     if (reuse!=i)
    231         return kCONTINUE;
     230    if (i >=0)
     231      {
     232       const Int_t reuse = (n/1000)%100; // Force this to be 1!
     233       if (reuse!=i)
     234           return kCONTINUE;
     235      }
    232236
    233237    // This seems to be special to mmcs

  • trunk/Mars/msim/MPhotonEvent.cc

    r9949 r10060  
    507507        operator[](i).SimWavelength(wmin, wmax);
    508508}
    509 
    510 // --------------------------------------------------------------------------
    511 //
    512 // Read the Event section from the file
    513 //
    514 Int_t MPhotonEvent::ReadCorsikaEvt(MCorsikaFormat *fInFormat, Int_t id)
    515 {
    516     Int_t n = 0;
    517 
    518     // --- old I/O ---
    519     // Read only + Reflector (no absorption)
    520     // Muons:   1.06GB/115s =  9.2MB/s (100kEvs)
    521     // Gammas:  1.57GB/275s =  5.7MB/s (  1kEvs)
    522 
    523     // Read only:
    524     // Gammas:  1.57GB/158s =  9.9MB/s (  1kEvs)
    525     // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)
    526 
    527     // --- new I/O ---
    528     // Read only (don't allocate storage space):
    529     // Gammas:  1.57GB/143s = 11.0MB/s (  1kEvs)
    530     // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)
    531 
    532     // Read only in blocks (with storage allocation):
    533     // Gammas:  1.57GB/28s  =  56MB/s (  1kEvs)
    534     // Muons:   1.06GB/5.2s = 204MB/s (100kEvs)
    535 
    536     // Read only in blocks (without storage allocation):
    537     // similar to just copy
    538 
    539     // Copy with cp
    540     // 1.57GB/ 5s   CPU
    541     // 1.57GB/28s   REAL
    542     // 1.06GB/ 3s   CPU
    543     // 1.06GB/22s   REAL
    544     Float_t *buffer = 0;
    545 
    546     if (fInFormat->IsEventioFormat())
    547     {
    548         while (1)
    549         {
    550             const Int_t rc = fInFormat->GetNextEvent(&buffer, id);
    551             if (rc==kERROR)
    552                 return kERROR;
    553             if (rc==kFALSE)
    554                 break;
    555 
    556             // Loop over number of photons in bunch
    557             while (Add(n).FillEventIO(buffer))
    558                 n++;
    559         }
    560     }
    561     else
    562     {
    563         while (1)
    564         {
    565             const Int_t rc1 = fInFormat->GetNextEvent(&buffer);
    566             if (rc1==kERROR)
    567                 return kERROR;
    568             if (rc1==kFALSE)
    569                 break;
    570 
    571             const Int_t rc2 = Add(n).FillCorsika(buffer, id);
    572             switch (rc2)
    573             {
    574             case kCONTINUE:  continue;        // No data in this bunch... skip it.
    575             case kERROR:     return kERROR;   // Error occured
    576             //case kFALSE:     return kFALSE;   // End of stream
    577             }
    578 
    579             // This is a photon we would like to keep later.
    580             // Increase the counter by one
    581             n++;
    582         }
    583     }
    584 
    585 /*
    586 
    587     while (1)
    588     {
    589         Float_t buffer[273];
    590         Float_t * ptr = buffer;
    591 
    592 
    593         if (!fInFormat->ReadData(273, buffer))
    594             return kFALSE;
    595 
    596         if (!memcmp(ptr, "EVTE", 4))
    597             {
    598 
    599             fInFormat->UnreadLastData();
    600             break;
    601             }
    602 
    603         Float_t *end = ptr + 273;
    604 
    605         // Loop over all sub-blocks (photons) in the block and if
    606         // they contain valid data add them to the array
    607         while (ptr<end)
    608         {
    609             // Get/Add the n-th entry from the array and
    610             // fill it with the current 7 floats
    611             const Int_t rc = Add(n).FillCorsika(ptr);
    612             ptr += 7;
    613 
    614             switch (rc)
    615             {
    616             case kCONTINUE:  continue;        // No data in this bunch... skip it.
    617             case kERROR:     return kERROR;   // Error occured
    618             //case kFALSE:     return kFALSE;   // End of stream
    619             }
    620 
    621             // This is a photon we would like to keep later.
    622             // Increase the counter by one
    623             n++;
    624         }
    625     }
    626 
    627 */
    628     Resize(n);
    629     fData.UnSort();
    630 
    631     SetReadyToSave();
    632 
    633     //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
    634     return kTRUE;
    635 
    636 }
    637 
    638 Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin, Int_t i)
    639 {
    640     Int_t n = 0;
    641 
    642     // --- old I/O ---
    643     // Read only + Reflector (no absorption)
    644     // Muons:   1.06GB/115s =  9.2MB/s (100kEvs)
    645     // Gammas:  1.57GB/275s =  5.7MB/s (  1kEvs)
    646 
    647     // Read only:
    648     // Gammas:  1.57GB/158s =  9.9MB/s (  1kEvs)
    649     // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)
    650 
    651     // --- new I/O ---
    652     // Read only (don't allocate storage space):
    653     // Gammas:  1.57GB/143s = 11.0MB/s (  1kEvs)
    654     // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)
    655 
    656     // Read only in blocks (with storage allocation):
    657     // Gammas:  1.57GB/28s  =  56MB/s (  1kEvs)
    658     // Muons:   1.06GB/5.2s = 204MB/s (100kEvs)
    659 
    660     // Read only in blocks (without storage allocation):
    661     // similar to just copy
    662 
    663     // Copy with cp
    664     // 1.57GB/ 5s   CPU
    665     // 1.57GB/28s   REAL
    666     // 1.06GB/ 3s   CPU
    667     // 1.06GB/22s   REAL
    668 
    669     while (1)
    670     {
    671         // Stprage for one block
    672         char c[273*4];
    673 
    674         // Read the first four byte to check whether the next block
    675         // doen't belong to the event anymore
    676         fin.read(c, 4);
    677         if (!fin)
    678             return kFALSE;
    679 
    680         // Check if the event is finished
    681         if (!memcmp(c, "EVTE", 4))
    682             break;
    683 
    684         // Now read the rest of the data
    685         fin.read(c+4, 272*4);
    686 
    687         Float_t *ptr = reinterpret_cast<Float_t*>(c);
    688         Float_t *end = ptr + 273;
    689 
    690         // Loop over all sub-blocks (photons) in the block and if
    691         // they contain valid data add them to the array
    692         while (ptr<end)
    693         {
    694             // Get/Add the n-th entry from the array and
    695             // fill it with the current 7 floats
    696             const Int_t rc = Add(n).FillCorsika(ptr, i);
    697             ptr += 7;
    698 
    699             switch (rc)
    700             {
    701             case kCONTINUE:  continue;        // No data in this bunch... skip it.
    702             case kERROR:     return kERROR;   // Error occured
    703             //case kFALSE:     return kFALSE;   // End of stream
    704             }
    705 
    706             // This is a photon we would like to keep later.
    707             // Increase the counter by one
    708             n++;
    709         }
    710     }
    711 /*
    712     while (1)
    713     {
    714         // Check the first four bytes
    715         char c[4];
    716         fin.read(c, 4);
    717 
    718         // End of stream
    719         if (!fin)
    720             return kFALSE;
    721 
    722         // Check if we found the end of the event
    723         if (!memcmp(c, "EVTE", 4))
    724             break;
    725 
    726         // The first for byte contained data already --> go back
    727         fin.seekg(-4, ios::cur);
    728 
    729         // Do not modify this. It is optimized for execution
    730         // speed and flexibility!
    731         MPhotonData &ph = Add(n);
    732         // It checks how many entries the lookup table has. If it has enough
    733         // entries and the entry was already allocated, we can re-use it,
    734         // otherwise we have to allocate it.
    735 
    736         // Now we read a single cherenkov bunch. Note that for speed reason we have not
    737         // called the constructor if the event was already constructed (virtual table
    738         // set), consequently we must make sure that ReadCorsikaEvent does reset
    739         // all data mebers no matter whether they are read or not.
    740         const Int_t rc = ph.ReadCorsikaEvt(fin);
    741 
    742         // Evaluate result from reading event
    743         switch (rc)
    744         {
    745         case kCONTINUE:  continue;        // No data in this bunch... skip it.
    746         case kFALSE:     return kFALSE;   // End of stream
    747         case kERROR:     return kERROR;   // Error occured
    748         }
    749 
    750         // FIXME: If fNumPhotons!=1 add the photon more than once
    751 
    752         // Now increase the number of entries which are kept,
    753         // i.e. keep this photon(s)
    754         n++;
    755     }
    756   */
    757 
    758     Resize(n);
    759     fData.UnSort();
    760 
    761     SetReadyToSave();
    762 
    763     //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
    764     return kTRUE;
    765 }
    766 
    767 // --------------------------------------------------------------------------
    768 /*
    769 Int_t MPhotonEvent::ReadRflEvt(std::istream &fin, Int_t i)
    770 {
    771     Int_t n = 0;
    772 
    773     while (1)
    774     {
    775         // Check the first four bytes
    776         char c[13];
    777         fin.read(c, 13);
    778 
    779         // End of stream
    780         if (!fin)
    781             return kFALSE;
    782 
    783         // Check if we found the end of the event
    784         if (!memcmp(c, "\nEND---EVENT\n", 13))
    785             break;
    786 
    787         // The first for byte contained data already --> go back
    788         fin.seekg(-13, ios::cur);
    789 
    790         // Do not modify this. It is optimized for execution
    791         // speed and flexibility!
    792         //TObject *o = n<fData.GetSize() && fData.UncheckedAt(n) ? fData.UncheckedAt(n) : fData.New(n);
    793 
    794         // Now we read a single cherenkov bunch
    795         //const Int_t rc = static_cast<MPhotonData*>(o)->ReadRflEvt(fin);
    796         const Int_t rc = Add(n).ReadRflEvt(fin, i);
    797 
    798         // Evaluate result from reading event
    799         switch (rc)
    800         {
    801         case kCONTINUE:  continue;        // No data in this bunch... skip it.
    802         case kFALSE:     return kFALSE;   // End of stream
    803         case kERROR:     return kERROR;   // Error occured
    804         }
    805 
    806         // Now increase the number of entries which are kept,
    807         // i.e. keep this photon(s)
    808         n++;
    809     }
    810 
    811     Shrink(n);
    812 
    813     SetReadyToSave();
    814 
    815     // *fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
    816     return kTRUE;
    817 }
    818 */
     509Int_t MPhotonEvent::ReadEventIoEvt(MCorsikaFormat *fInFormat)
     510{
     511   Int_t  bunchHeader[3];
     512   fInFormat->Read(bunchHeader, 3 * sizeof(Int_t));
     513
     514   Int_t n = 0;
     515
     516   for (int bunch = 0; bunch < bunchHeader[2]; bunch++)
     517      {
     518      Float_t buffer[8];
     519      fInFormat->Read(buffer, 8 * sizeof(Float_t));
     520
     521      if (Add(n).FillEventIO(buffer))
     522         n++;
     523      }
     524
     525   Resize(n);
     526   fData.UnSort();
     527
     528   SetReadyToSave();
     529
     530   //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
     531   return kTRUE;
     532
     533}
     534
     535Int_t MPhotonEvent::ReadCorsikaEvt(Float_t * data, Int_t numEvents, Int_t arrayIdx)
     536{
     537   Int_t n = 0;
     538
     539   for (Int_t event = 0; event < numEvents; event++)
     540      {
     541      const Int_t rc = Add(n).FillCorsika(data + 7 * event, arrayIdx);
     542
     543      switch (rc)
     544      {
     545      case kCONTINUE:  continue;        // No data in this bunch... skip it.
     546      case kERROR:     return kERROR;   // Error occured
     547      }
     548
     549      // This is a photon we would like to keep later.
     550      // Increase the counter by one
     551      n++;
     552      }
     553
     554   Resize(n);
     555   fData.UnSort();
     556
     557   SetReadyToSave();
     558
     559   //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
     560   return kTRUE;
     561
     562}
    819563
    820564// --------------------------------------------------------------------------

  • trunk/Mars/msim/MPhotonEvent.h

    r9941 r10060  
    6060
    6161    // I/O
    62     Int_t ReadCorsikaEvt(MCorsikaFormat *fInFormat, Int_t i);
    63     Int_t ReadCorsikaEvt(istream &fin, Int_t i);
    64     //Int_t ReadRflEvt(istream &fin, Int_t i);
     62    Int_t ReadEventIoEvt(MCorsikaFormat *fInFormat);
     63    Int_t ReadCorsikaEvt(Float_t * data, Int_t numEvents, Int_t arrayIdx);
    6564
    6665    // TObject

  • trunk/Mars/readcorsika.cc

    r9949 r10060  
    3838    gLog << all << endl;
    3939    gLog << "Sorry the usage is:" << endl;
    40     gLog << "   readcorsika [-h] [-?] [-vn] [-dec] [-a0] inputfile[.raw] [outputfile[.root]]" << endl << endl;
    41     gLog << "     inputfile:   corsika raw file" << endl;
    42     gLog << "     outputfile:  root file" << endl;
     40    gLog << "   readcorsika [-h] [-?] [-vn] [-dec] [-a0] [-A=arrayNo] [-T=telescopeNo] inputfile[.raw] [outputfile[.root]]" << endl << endl;
     41    gLog << "     inputfile:     corsika raw file or eventio file" << endl;
     42    gLog << "     outputfile:    root file" << endl;
     43    gLog << "   -A=arrayNo:      use data only of array number arrayNo" << endl;
     44    gLog << "   -T=telescopeNo:  use data only of telescope number telescopeNo (used only for eventio input file)" << endl;
    4345    gLog << "   -ff                       Force reading of file even if problems occur" << endl;
    4446    gLog.Usage();
    … …  
    7173
    7274    const Int_t  kCompLvl = arg.GetIntAndRemove("--comp=", 1);
     75    const Int_t  kArrayNo = arg.GetIntAndRemove("-A=", -1);
     76    const Int_t  kTelNo   = arg.GetIntAndRemove("-T=", -1);
    7377    const Bool_t kForce   = arg.HasOnlyAndRemove("-f");
    7478    const Bool_t kForceRd = arg.HasOnlyAndRemove("-ff");
     79
    7580
    7681    //
    … …  
    149154    MCorsikaRead read(kNamein);
    150155    read.SetForceMode(kForceRd);
     156    read.SetArrayIdx(kArrayNo);
     157    read.SetTelescopeIdx(kTelNo);
    151158    tasks.AddToList(&read);
    152159
Note: See TracChangeset for help on using the changeset viewer.