source: trunk/Mars/mcorsika/MCorsikaFormat.cc @ 10060

Last change on this file since 10060 was 10060, checked in by rohlfs, 10 years ago
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.
File size: 6.7 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18!   Author(s): Reiner Rohlfs 2010
19!   Author(s): Thomas Bretz  2010 <mailto:thomas.bretz@epfl.ch>
20!
21!   Copyright: Software Development, 2000-2010
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27//
28// MCorsikaFormat
29//
30//////////////////////////////////////////////////////////////////////////////
31#include "MCorsikaFormat.h"
32
33#include <errno.h>
34#include <fstream>
35
36#include "MLogManip.h"
37
38
39using namespace std;
40
41
42const unsigned int MCorsikaFormat::kSyncMarker = 0xd41f8a37;
43
44// --------------------------------------------------------------------------
45//
46MCorsikaFormat *MCorsikaFormat::CorsikaFormatFactory(const char * fileName)
47{
48    ifstream * fileIn = new ifstream(fileName);
49
50    const Bool_t noexist = !(*fileIn);
51    if (noexist)
52    {
53        gLog << err << "Cannot open file " << fileName << ": ";
54        gLog << (errno!=0?strerror(errno):"Insufficient memory for decompression") << endl;
55        delete fileIn;
56        return NULL;
57    }
58
59    char *buffer = new char[5];
60    memset(buffer, 0, 5);
61    fileIn->read(buffer, 4);
62    fileIn->seekg(-4, ios::cur);
63
64    if (strcmp(buffer, "RUNH") == 0)
65    {
66        delete [] buffer;
67        return new MCorsikaFormatRaw(fileIn);
68    }
69
70    if (*reinterpret_cast<unsigned int*>(buffer) == kSyncMarker)
71    {
72        delete [] buffer;
73        return new MCorsikaFormatEventIO(fileIn);
74    }
75
76    gLog << err << "File " << fileName <<
77            " is neither a CORSIKA raw nor EventIO file" << endl;
78    delete fileIn;
79    delete [] buffer;
80
81    return NULL;
82}
83
84Bool_t MCorsikaFormat::Read(void *ptr, Int_t i) const
85{
86   fIn->read((char*)ptr, i);
87   return !fIn->fail();
88
89}
90// --------------------------------------------------------------------------
91//
92Bool_t MCorsikaFormat::Eof() const
93{
94    return fIn->eof();
95}
96
97// --------------------------------------------------------------------------
98//
99MCorsikaFormat::~MCorsikaFormat() 
100{
101    delete fIn;
102}
103
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;
118    fIn->read(blockHeader, 4);
119    if (fIn->eof())
120        return kFALSE;
121
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}
143// --------------------------------------------------------------------------
144//
145Bool_t MCorsikaFormatRaw::SeekEvtEnd()
146{
147    // Search subblockwise backward (Block: 5733*4 = 21*273*4)
148    for (int i=1; i<22; i++)
149    {
150        fIn->seekg(-i*273*4, ios::end);
151
152        char runh[5]="\0\0\0\0";
153        fIn->read(runh, 4);
154
155        if (!strcmp(runh, "RUNE"))
156        {
157//            fIn->seekg(-4, ios::cur);
158            return kTRUE;
159        }
160    }
161   
162    return kTRUE;
163}
164///////////////////////////////////////////////////////////////////////////////
165///////////////////////////////////////////////////////////////////////////////
166
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      }
218    return kTRUE;
219}
220
221// --------------------------------------------------------------------------
222//
223Bool_t MCorsikaFormatEventIO::SeekEvtEnd()
224{
225
226    // the RUNE block it at the very end of the file.
227    fIn->seekg(-32, ios::end);
228
229    unsigned int blockHeader[4];
230    fIn->read((char*)blockHeader, 4 * sizeof(int));
231
232    if ( blockHeader[0]               == kSyncMarker &&
233        (blockHeader[1] & 0xffff)     == 1210        &&
234        (blockHeader[3] & 0x3fffffff) == 16)
235    {
236        // this seams to be a RUNE (1210)  block
237        fIn->seekg( 8, ios::cur);
238        return kTRUE;
239    }
240
241    return kFALSE;
242}
243
Note: See TracBrowser for help on using the repository browser.