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 | using namespace std;
|
---|
39 |
|
---|
40 | const unsigned int MCorsikaFormat::kSyncMarker = 0xd41f8a37;
|
---|
41 |
|
---|
42 | // --------------------------------------------------------------------------
|
---|
43 | //
|
---|
44 | MCorsikaFormat *MCorsikaFormat::CorsikaFormatFactory(const char * fileName)
|
---|
45 | {
|
---|
46 | ifstream * fileIn = new ifstream(fileName);
|
---|
47 |
|
---|
48 | const Bool_t noexist = !(*fileIn);
|
---|
49 | if (noexist)
|
---|
50 | {
|
---|
51 | gLog << err << "Cannot open file " << fileName << ": ";
|
---|
52 | gLog << (errno!=0?strerror(errno):"Insufficient memory for decompression") << endl;
|
---|
53 | delete fileIn;
|
---|
54 | return NULL;
|
---|
55 | }
|
---|
56 |
|
---|
57 | char buffer[5]="\0\0\0\0";
|
---|
58 | fileIn->read(buffer, 4);
|
---|
59 | fileIn->seekg(-4, ios::cur);
|
---|
60 |
|
---|
61 | if (strcmp(buffer, "RUNH") == 0)
|
---|
62 | return new MCorsikaFormatRaw(fileIn);
|
---|
63 |
|
---|
64 | if (*reinterpret_cast<unsigned int*>(buffer) == kSyncMarker)
|
---|
65 | return new MCorsikaFormatEventIO(fileIn);
|
---|
66 |
|
---|
67 | gLog << err << "File " << fileName <<
|
---|
68 | " is neither a CORSIKA raw nor EventIO file" << endl;
|
---|
69 | delete fileIn;
|
---|
70 |
|
---|
71 | return NULL;
|
---|
72 | }
|
---|
73 |
|
---|
74 | void MCorsikaFormat::Read(void *ptr, Int_t i) const
|
---|
75 | {
|
---|
76 | fIn->read((char*)ptr, i);
|
---|
77 | }
|
---|
78 |
|
---|
79 | // --------------------------------------------------------------------------
|
---|
80 | //
|
---|
81 | Bool_t MCorsikaFormat::Eof() const
|
---|
82 | {
|
---|
83 | return fIn->eof();
|
---|
84 | }
|
---|
85 |
|
---|
86 | // --------------------------------------------------------------------------
|
---|
87 | //
|
---|
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 | //
|
---|
138 | MCorsikaFormat::~MCorsikaFormat()
|
---|
139 | {
|
---|
140 | delete fIn;
|
---|
141 | }
|
---|
142 |
|
---|
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 | Bool_t MCorsikaFormatRaw::SeekNextBlock(const char * id, unsigned short type) const
|
---|
154 | {
|
---|
155 | char blockHeader[5]="\0\0\0\0";
|
---|
156 | 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(id, "EVTH")==0 && strcmp(blockHeader, "RUNE")==0)
|
---|
165 | return kTRUE;
|
---|
166 |
|
---|
167 | gLog << err << "ERROR - Wrong identifier: " << id << " expected.";
|
---|
168 | gLog << " But read " << blockHeader << " from file." << endl;
|
---|
169 |
|
---|
170 | return kFALSE;
|
---|
171 | }
|
---|
172 |
|
---|
173 | // --------------------------------------------------------------------------
|
---|
174 | //
|
---|
175 | void MCorsikaFormatRaw::UnreadLastHeader() const
|
---|
176 | {
|
---|
177 | fIn->seekg(-4, ios::cur);
|
---|
178 | }
|
---|
179 |
|
---|
180 | // --------------------------------------------------------------------------
|
---|
181 | //
|
---|
182 | Bool_t MCorsikaFormatRaw::SeekEvtEnd()
|
---|
183 | {
|
---|
184 |
|
---|
185 | // Search subblockwise backward (Block: 5733*4 = 21*273*4)
|
---|
186 | for (int i=1; i<22; i++)
|
---|
187 | {
|
---|
188 | fIn->seekg(-i*273*4, ios::end);
|
---|
189 |
|
---|
190 | char runh[5]="\0\0\0\0";
|
---|
191 | fIn->read(runh, 4);
|
---|
192 |
|
---|
193 | if (!strcmp(runh, "RUNE"))
|
---|
194 | {
|
---|
195 | fIn->seekg(-4, ios::cur);
|
---|
196 | return kTRUE;
|
---|
197 | }
|
---|
198 | }
|
---|
199 |
|
---|
200 | return kTRUE;
|
---|
201 | }
|
---|
202 |
|
---|
203 | // --------------------------------------------------------------------------
|
---|
204 | //
|
---|
205 | // Returns one event (7 Float values) after the other until the EventEnd
|
---|
206 | // block is found.
|
---|
207 | // If a read error occurred, the readError is set to kTRUE and the function
|
---|
208 | // returns kFALSE;
|
---|
209 | //
|
---|
210 | Int_t MCorsikaFormatRaw::GetNextEvent(Float_t **buffer, Int_t)
|
---|
211 | {
|
---|
212 | static Float_t data[273];
|
---|
213 | static Float_t * next = data + 273;
|
---|
214 |
|
---|
215 | if (next == data + 273)
|
---|
216 | {
|
---|
217 | // read next block of events
|
---|
218 | if (!ReadData(273, data))
|
---|
219 | return kERROR;
|
---|
220 |
|
---|
221 | if (!memcmp(data, "EVTE", 4))
|
---|
222 | {
|
---|
223 | // we found the end of list of events
|
---|
224 | UnreadLastData();
|
---|
225 | return kFALSE;
|
---|
226 | }
|
---|
227 |
|
---|
228 | next = data;
|
---|
229 | }
|
---|
230 |
|
---|
231 | *buffer = next;
|
---|
232 | next += 7;
|
---|
233 |
|
---|
234 | return kTRUE;
|
---|
235 | }
|
---|
236 |
|
---|
237 | ///////////////////////////////////////////////////////////////////////////////
|
---|
238 | ///////////////////////////////////////////////////////////////////////////////
|
---|
239 |
|
---|
240 | // --------------------------------------------------------------------------
|
---|
241 | //
|
---|
242 | // Jumps from one top level object to the next until if finds the object with
|
---|
243 | // correct type and correct id. The id is identifier of the raw corsika block,
|
---|
244 | // for example "RUNH", RUNE", "EVTH"
|
---|
245 | //
|
---|
246 | // Returns kTRUE if the functions finds the type / id
|
---|
247 | // kFALSE if the functions does not find the type / id.
|
---|
248 | //
|
---|
249 | // After this call the file position pointer points just after the 4 bytes
|
---|
250 | // of the id.
|
---|
251 | //
|
---|
252 | Bool_t MCorsikaFormatEventIO::SeekNextBlock(const char * id, unsigned short type) const
|
---|
253 | {
|
---|
254 | cout << "Seek " << type << endl;
|
---|
255 | while (1)
|
---|
256 | {
|
---|
257 | // we read - synchronisation marker
|
---|
258 | // - type / version field
|
---|
259 | // - identification field
|
---|
260 | // - length
|
---|
261 | // - unknown field
|
---|
262 | // - id (first 4 bytes of data field)
|
---|
263 | int blockHeader[4];
|
---|
264 | fIn->read((char*)blockHeader, 4 * sizeof(int));
|
---|
265 |
|
---|
266 | if (fIn->eof())
|
---|
267 | {
|
---|
268 | gLog << err << "MCorsikaFormatEventIO::SeekNextBlock: ERROR - Unexpected end-of-file seeking " << id << " (" << type << ")." << endl;
|
---|
269 | return kFALSE;
|
---|
270 | }
|
---|
271 |
|
---|
272 | const unsigned short objType = blockHeader[1] & 0xFFFF;
|
---|
273 | cout << "objType=" << objType << endl;
|
---|
274 | if (type == objType)
|
---|
275 | {
|
---|
276 | if (!id)
|
---|
277 | break;
|
---|
278 |
|
---|
279 | char c[9];
|
---|
280 | fIn->read(c, 8);
|
---|
281 |
|
---|
282 | if (memcmp(c+4, id, 4)==0)
|
---|
283 | break;
|
---|
284 |
|
---|
285 | c[8] = 0;
|
---|
286 | gLog << err << "MCorsikaFormatEventIO::SeekNextBlock: ERROR - '" << c+4 << "' doesn't match expected one '" << id << "' in object " << type << endl;
|
---|
287 | return kFALSE;
|
---|
288 | }
|
---|
289 |
|
---|
290 | // we are looking for a event header, but found a runEnd.
|
---|
291 | // This will happen at the end of a run. We stop here looking for
|
---|
292 | // the next event header
|
---|
293 | if (type == 1202 && objType == 1210)
|
---|
294 | return kFALSE;
|
---|
295 |
|
---|
296 | // a unknown block, we jump to the next one
|
---|
297 | const int length = blockHeader[3] & 0x3fffffff;
|
---|
298 | fIn->seekg(length, ios::cur);
|
---|
299 | }
|
---|
300 |
|
---|
301 |
|
---|
302 | return kTRUE;
|
---|
303 | }
|
---|
304 |
|
---|
305 | // --------------------------------------------------------------------------
|
---|
306 | //
|
---|
307 | void MCorsikaFormatEventIO::UnreadLastHeader() const
|
---|
308 | {
|
---|
309 | fIn->seekg( -6 * (Int_t)(sizeof(int)), ios::cur);
|
---|
310 | }
|
---|
311 |
|
---|
312 | // --------------------------------------------------------------------------
|
---|
313 | //
|
---|
314 | Bool_t MCorsikaFormatEventIO::SeekEvtEnd()
|
---|
315 | {
|
---|
316 | if (fRunePos != streampos(0))
|
---|
317 | {
|
---|
318 | fIn->seekg(fRunePos, ios::beg);
|
---|
319 | return kTRUE;
|
---|
320 | }
|
---|
321 |
|
---|
322 | // it is the first time we are looking for the RUNE block
|
---|
323 |
|
---|
324 | // is the RUNE block at the very end of the file?
|
---|
325 | const streampos currentPos = fIn->tellg();
|
---|
326 |
|
---|
327 | fIn->seekg(-32, ios::end);
|
---|
328 |
|
---|
329 | unsigned int blockHeader[4];
|
---|
330 | fIn->read((char*)blockHeader, 4 * sizeof(int));
|
---|
331 |
|
---|
332 | if ( blockHeader[0] == kSyncMarker &&
|
---|
333 | (blockHeader[1] & 0xffff) == 1210 &&
|
---|
334 | (blockHeader[3] & 0x3fffffff) == 16)
|
---|
335 | {
|
---|
336 | // this seams to be a RUNE (1210) block
|
---|
337 | fIn->seekg( -4 * (Int_t)(sizeof(int)), ios::cur);
|
---|
338 | fRunePos = fIn->tellg();
|
---|
339 | return kTRUE;
|
---|
340 | }
|
---|
341 |
|
---|
342 | // we do not find a RUNE block at the end of the file
|
---|
343 | // we have to search in the file
|
---|
344 | fIn->seekg(currentPos, ios::beg);
|
---|
345 | if (!SeekNextBlock("RUNE", 1210))
|
---|
346 | return kFALSE;
|
---|
347 |
|
---|
348 | UnreadLastHeader();
|
---|
349 | fRunePos = fIn->tellg();
|
---|
350 |
|
---|
351 | return kTRUE;
|
---|
352 | }
|
---|
353 |
|
---|
354 | // --------------------------------------------------------------------------
|
---|
355 | //
|
---|
356 | // Returns one event (7 Float values) after the other until the EventEnd
|
---|
357 | // block is found.
|
---|
358 | // If a read error occurred, the readError is set to kTRUE and the function
|
---|
359 | // returns kFALSE;
|
---|
360 | //
|
---|
361 | Int_t MCorsikaFormatEventIO::GetNextEvent(Float_t **buffer, Int_t telescope)
|
---|
362 | {
|
---|
363 | static Float_t *data = NULL;
|
---|
364 |
|
---|
365 | static int lengthPhotonData = 0;
|
---|
366 | static int lengthEvent = 0;
|
---|
367 |
|
---|
368 | if (lengthPhotonData>0)
|
---|
369 | {
|
---|
370 | cout << "Return Bunch l2=" << lengthPhotonData << endl;
|
---|
371 |
|
---|
372 | lengthPhotonData -= 8 * sizeof(Float_t);
|
---|
373 | *buffer += 8; // Return the pointer to the start of the current photon data
|
---|
374 | return kTRUE;
|
---|
375 | }
|
---|
376 |
|
---|
377 | if (data)
|
---|
378 | {
|
---|
379 | delete [] data;
|
---|
380 | data = NULL;
|
---|
381 | cout << "Return end-of-tel LE=" << lengthEvent << endl;
|
---|
382 | return kFALSE;
|
---|
383 | }
|
---|
384 |
|
---|
385 | if (lengthEvent==0)
|
---|
386 | {
|
---|
387 | cout << "Searching 1204... " << flush;
|
---|
388 | // The length of the block is stored and we use it to determine
|
---|
389 | // when the next top level block is expected
|
---|
390 | const Int_t rc = NextEventObject(lengthEvent);
|
---|
391 | if (rc==kERROR)
|
---|
392 | return kERROR;
|
---|
393 | if (!rc)
|
---|
394 | {
|
---|
395 | cout << "skip EVTE." << endl;
|
---|
396 | return kFALSE;
|
---|
397 | }
|
---|
398 |
|
---|
399 | cout << "found new array." << endl;
|
---|
400 | }
|
---|
401 |
|
---|
402 | cout << "Searching 1205..." << flush;
|
---|
403 |
|
---|
404 | // Look for the next event photon bunch (object type 1205)
|
---|
405 | const Int_t tel = NextPhotonObject(lengthPhotonData);
|
---|
406 | if (tel<0)
|
---|
407 | return kERROR;
|
---|
408 |
|
---|
409 | lengthEvent -= lengthPhotonData+12; // Length of data + Length of header
|
---|
410 |
|
---|
411 | lengthPhotonData -= 12; // Length of data before the photon bunches
|
---|
412 | fIn->seekg(12, ios::cur); // Skip this data
|
---|
413 |
|
---|
414 | cout << " found (tel=" << tel << ",LEN=" << lengthEvent << ",len=" << lengthPhotonData << ")... " << flush;
|
---|
415 |
|
---|
416 | if (lengthPhotonData==0)
|
---|
417 | {
|
---|
418 | cout << "empty!" << endl;
|
---|
419 | return kFALSE;
|
---|
420 | }
|
---|
421 |
|
---|
422 | const UInt_t cnt = lengthPhotonData / sizeof(Float_t);
|
---|
423 | // Read next object (photon bunch)
|
---|
424 | data = new Float_t[cnt];
|
---|
425 | if (!ReadData(cnt, data, 0))
|
---|
426 | {
|
---|
427 | delete [] data;
|
---|
428 | data = NULL;
|
---|
429 | return kERROR;
|
---|
430 | }
|
---|
431 |
|
---|
432 | cout << "return!" << endl;
|
---|
433 |
|
---|
434 | lengthPhotonData -= 8 * sizeof(Float_t);
|
---|
435 | *buffer = data; // Return the pointer to the start of the current photon data
|
---|
436 |
|
---|
437 | return kTRUE;
|
---|
438 | }
|
---|
439 |
|
---|
440 | // --------------------------------------------------------------------------
|
---|
441 | //
|
---|
442 | // Looks for the next Block with type 1204 and return kTRUE.
|
---|
443 | // The function also stops moving forward in the file, if it finds a
|
---|
444 | // EventEnd block (1209). In this case kFALSE is returned
|
---|
445 | //
|
---|
446 | Int_t MCorsikaFormatEventIO::NextEventObject(Int_t &length) const
|
---|
447 | {
|
---|
448 | while (1)
|
---|
449 | {
|
---|
450 | // we read - synchronisation marker
|
---|
451 | // - type / version field
|
---|
452 | // - identification field
|
---|
453 | // - length
|
---|
454 | UInt_t blockHeader[4];
|
---|
455 | fIn->read((char*)blockHeader, 4 * sizeof(Int_t));
|
---|
456 |
|
---|
457 | if (fIn->eof())
|
---|
458 | {
|
---|
459 | gLog << err << "MCorsikaFormatEventIO::NextEventObject: ERROR - Unexpected end-of-file." << endl;
|
---|
460 | return kERROR;
|
---|
461 | }
|
---|
462 |
|
---|
463 | // For speed reasons we skip the check of the synchronization marker
|
---|
464 |
|
---|
465 | // Decode the object type
|
---|
466 | const unsigned short objType = blockHeader[1] & 0xFFFF;
|
---|
467 |
|
---|
468 | cout << "objType=" << objType << endl;
|
---|
469 |
|
---|
470 | // Decode length of block
|
---|
471 | length = blockHeader[3] & 0x3fffffff;
|
---|
472 |
|
---|
473 | // Check if we found the next array (reuse / impact parameter)
|
---|
474 | // blockHeader[2] == array number (reuse)
|
---|
475 | if (objType==1204)
|
---|
476 | return kTRUE;
|
---|
477 |
|
---|
478 | // we found an eventEnd block, reset file pointer
|
---|
479 | if (objType==1209)
|
---|
480 | break;
|
---|
481 |
|
---|
482 | // a unknown block, we jump to the next one
|
---|
483 | fIn->seekg(length, ios::cur);
|
---|
484 | }
|
---|
485 |
|
---|
486 | fIn->seekg( -4 * (Int_t)(sizeof(Int_t)), ios::cur);
|
---|
487 | length = 0;
|
---|
488 |
|
---|
489 | return kFALSE;
|
---|
490 | }
|
---|
491 |
|
---|
492 | // --------------------------------------------------------------------------
|
---|
493 | //
|
---|
494 | Int_t MCorsikaFormatEventIO::NextPhotonObject(Int_t &length) const
|
---|
495 | {
|
---|
496 | UInt_t blockHeader[3];
|
---|
497 |
|
---|
498 | // we read - synchronisation marker
|
---|
499 | // - type / version field
|
---|
500 | // - identification field
|
---|
501 | // - length
|
---|
502 | fIn->read((char*)blockHeader, 3 * sizeof(UInt_t));
|
---|
503 |
|
---|
504 | if (fIn->eof())
|
---|
505 | {
|
---|
506 | gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected end-of-file.." << endl;
|
---|
507 | return -1;
|
---|
508 | }
|
---|
509 |
|
---|
510 | const unsigned short objType = blockHeader[0] & 0xFFFF;
|
---|
511 | if (objType != 1205)
|
---|
512 | {
|
---|
513 | gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected EventIO object type " << objType << " (expected 1205)." << endl;
|
---|
514 | return -1;
|
---|
515 | }
|
---|
516 |
|
---|
517 | const unsigned short version = (blockHeader[0] >> 20) & 0x0FFF;
|
---|
518 | if (version != 0)
|
---|
519 | {
|
---|
520 | gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected EventIO object version " << version << " (expected 0 for object type 1205)." << endl;
|
---|
521 | return -1;
|
---|
522 | }
|
---|
523 |
|
---|
524 | length = blockHeader[2] & 0x3fffffff;
|
---|
525 |
|
---|
526 | // blockHeader[1] == 1000 x array number (reuse) + telescope number
|
---|
527 | return blockHeader[1] % 1000;
|
---|
528 | }
|
---|