source: branches/Corsika7500Compatibility/mcorsika/MCorsikaRead.cc@ 19004

Last change on this file since 19004 was 18526, checked in by dbaack, 8 years ago
Bugfixes and gcc-5 compability changes
File size: 21.4 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): Thomas Bretz 11/2008 <mailto:tbretz@astro.uni-wuerzburg.de>
19 Reiner Rohlfs 11/2010 <mailto:Reiner.Rohlfs@unige.ch>
20!
21! Copyright: Software Development, 2000-2008
22!
23!
24\* ======================================================================== */
25
26//////////////////////////////////////////////////////////////////////////////
27//
28// MCorsikaRead
29//
30// Input Containers:
31// -/-
32//
33// Output Containers:
34// MCorsikaRunHeader
35// MCorsikaEvtHeader
36// MPhotonEvent
37//
38//////////////////////////////////////////////////////////////////////////////
39#include "MCorsikaRead.h"
40
41#include <errno.h>
42#include <fstream>
43
44#include <TSystem.h>
45
46#include "MLog.h"
47#include "MLogManip.h"
48
49#include "MParList.h"
50#include "MStatusDisplay.h"
51
52#include "MCorsikaFormat.h"
53#include "MCorsikaRunHeader.h"
54#include "MCorsikaEvtHeader.h"
55
56#include "MPhotonEvent.h"
57
58ClassImp(MCorsikaRead);
59
60using namespace std;
61
62
63/* ----------- please don't delete and don't care about (Thomas) ------------
64#define kBUFSZ 64 //1024*1024*64
65#include <iomanip.h>
66class bifstream : public istream, public streambuf
67{
68private:
69 char fBuffer[kBUFSZ]; //!
70 FILE *fd;
71
72 int sync()
73 {
74 memset(fBuffer, 0, kBUFSZ);
75 return 0;
76 }
77 int underflow()
78 {
79 int sz=fread(fBuffer, kBUFSZ, 1, fd);
80 //int sz=fread(fBuffer, 1, kBUFSZ, fd);
81 setg(fBuffer, fBuffer, fBuffer+kBUFSZ);
82
83 return sz==1 ? *(unsigned char*)fBuffer : EOF;//EOF;
84 //return sz==kBUFSZ ? *(unsigned char*)fBuffer : EOF;//EOF;
85 }
86public:
87 bifstream(const char *name) : istream(this)
88 {
89 fd = fopen(name, "rb");
90 setbuf(fBuffer, kBUFSZ);
91 }
92};
93*/
94
95// --------------------------------------------------------------------------
96//
97// Default constructor. It tries to open the given file.
98//
99MCorsikaRead::MCorsikaRead(const char *fname, const char *name, const char *title)
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)
103{
104 fName = name ? name : "MRead";
105 fTitle = title ? title : "Read task to read DAQ binary files";
106
107 fFileNames = new TList;
108 fFileNames->SetOwner();
109
110 if (fname!=NULL)
111 AddFile(fname);
112}
113
114// --------------------------------------------------------------------------
115//
116// Destructor. Delete input stream.
117//
118MCorsikaRead::~MCorsikaRead()
119{
120 delete fFileNames;
121 if (fInFormat)
122 delete fInFormat;
123}
124
125/*
126Byte_t MCorsikaRead::IsFileValid(const char *name)
127{
128 MZlib fin(name);
129 if (!fin)
130 return 0;
131
132 Byte_t c[4];
133 fin.read((char*)c, 4);
134 if (!fin)
135 return 0;
136
137 if (c[0]!=0xc0)
138 return 0;
139
140 if (c[1]==0xc0)
141 return 1;
142
143 if (c[1]==0xc1)
144 return 2;
145
146 return 0;
147}
148*/
149
150// --------------------------------------------------------------------------
151//
152// Add a new file to a list of files to be processed, Returns the number
153// of files added. (We can enhance this with a existance chack and
154// wildcard support)
155//
156Int_t MCorsikaRead::AddFile(const char *fname, Int_t entries)
157{
158 TNamed *name = new TNamed(fname, "");
159 fFileNames->AddLast(name);
160 return 1;
161
162}
163// --------------------------------------------------------------------------
164//
165// This opens the next file in the list and deletes its name from the list.
166//
167Int_t MCorsikaRead::OpenNextFile(Bool_t print)
168{
169
170 //
171 // open the input stream and check if it is really open (file exists?)
172 //
173 if (fInFormat)
174 delete fInFormat;
175 fInFormat = NULL;
176
177 //
178 // Check for the existance of a next file to read
179 //
180 TObject *file = fFileNames->At(fNumFile);
181 if (!file)
182 return kFALSE;
183
184 //
185 // open the file which is the first one in the chain
186 //
187 const char *name = file->GetName();
188
189 const char *expname = gSystem->ExpandPathName(name);
190 fInFormat = MCorsikaFormat::CorsikaFormatFactory(expname);
191 delete [] expname;
192
193 if (fInFormat == NULL)
194 return kERROR;
195
196 *fLog << inf << "Open file: '" << name << "'" << endl;
197
198 if (fDisplay)
199 {
200 // Show the number of the last event after
201 // which we now open a new file
202 TString txt = GetFileName();
203 txt += " @ ";
204 txt += GetNumExecutions()-1;
205 fDisplay->SetStatusLine2(txt);
206 }
207
208 fNumFile++;
209
210 if (!fParList)
211 return kTRUE;
212
213 //
214 // Search for MTaskList
215 //
216 MTask *tlist = (MTask*)fParList->FindObject("MTaskList");
217 if (!tlist)
218 {
219 *fLog << err << dbginf << "MTaskList not found... abort." << endl;
220 return kERROR;
221 }
222
223 //
224 // A new file has been opened and new headers have been read.
225 // --> ReInit tasklist
226 //
227 return tlist->ReInit(fParList) ? kTRUE : kERROR;
228}
229
230// --------------------------------------------------------------------------
231//
232// Return file name of current file.
233//
234TString MCorsikaRead::GetFullFileName() const
235{
236 const TObject *file = fFileNames->At(fNumFile-1);
237 return file ? file->GetName() : "";
238}
239
240// --------------------------------------------------------------------------
241//
242// Restart with the first file
243//
244Bool_t MCorsikaRead::Rewind()
245{
246 fNumFile=0;
247 fNumEvents=0;
248 return OpenNextFile()==kTRUE;
249}
250
251// --------------------------------------------------------------------------
252//
253// Calculates the total number of events from all input files and stores
254// the number in fNumTotalEvents.
255//
256Bool_t MCorsikaRead::CalcNumTotalEvents()
257{
258 fNumTotalEvents = 0;
259
260 Bool_t rc = kTRUE;
261
262 while (1)
263 {
264 switch (OpenNextFile(kFALSE))
265 {
266 case kFALSE:
267 break;
268 case kERROR:
269 rc = kFALSE;
270 break;
271 case kTRUE:
272
273 // read run header(1200), telescope position(1201) and
274 // first event header(1202)
275 Bool_t status = kTRUE;
276 Int_t blockType, blockVersion, blockIdentifier, blockLength;
277 while (status &&
278 fInFormat->NextBlock(fReadState/* == 4*/, blockType, blockVersion,
279 blockIdentifier, blockLength))
280 {
281 if (blockType == 1200)
282 status = fRunHeader->ReadEvt(fInFormat);
283
284 else if(blockType == 1201)
285 status = ReadTelescopePosition();
286
287 else if (blockType == 1202)
288 {
289 Float_t buffer[272];
290 status = fInFormat->Read(buffer, 272 * sizeof(Float_t));
291 status = fRunHeader->ReadEventHeader(buffer);
292 break;
293 }
294 else
295 fInFormat->Seek(blockLength);
296 }
297
298 if (status != kTRUE)
299 return status;
300
301 if (!fInFormat->SeekEvtEnd())
302 {
303 *fLog << (fForceMode?warn:err) << "Error: RUNE section not found in file." << endl;
304 if (!fForceMode)
305 return fForceMode ? kTRUE : kFALSE;
306 }
307
308 if (!fRunHeader->ReadEvtEnd(fInFormat, kTRUE))
309 {
310 *fLog << (fForceMode?warn:err) << "Error: Reading RUNE section failed." << endl;
311 if (!fForceMode)
312 return kFALSE;
313 }
314
315 fNumTotalEvents += fRunHeader->GetNumEvents()*fRunHeader->GetNumReuse()*
316 (fTelescopeIdx<0 ? fNumTelescopes : 1);
317 continue;
318 }
319 break;
320 }
321
322 return rc;
323}
324
325// --------------------------------------------------------------------------
326//
327// Reads the the position of all telescopes in one array
328//
329Int_t MCorsikaRead::ReadTelescopePosition()
330{
331 if (!fInFormat->Read(&fNumTelescopes, 4)) return kERROR;
332
333 if (fTelescopeIdx>=fNumTelescopes)
334 {
335 *fLog << err << "ERROR - Requested telescope index " << fTelescopeIdx <<
336 " exceeds number of telescopes " << fNumTelescopes << " in file." << endl;
337 return kERROR;
338 }
339
340 fTelescopeX.Set(fNumTelescopes);
341 fTelescopeY.Set(fNumTelescopes);
342 fTelescopeZ.Set(fNumTelescopes);
343 fTelescopeR.Set(fNumTelescopes);
344
345 if (!fInFormat->Read(fTelescopeX.GetArray(), 4*fNumTelescopes)) return kERROR;
346 if (!fInFormat->Read(fTelescopeY.GetArray(), 4*fNumTelescopes)) return kERROR;
347 if (!fInFormat->Read(fTelescopeZ.GetArray(), 4*fNumTelescopes)) return kERROR;
348 if (!fInFormat->Read(fTelescopeR.GetArray(), 4*fNumTelescopes)) return kERROR;
349
350 *fLog << all;
351 *fLog << "Number of telescopes: " << fNumTelescopes << " (using ";
352 if (fTelescopeIdx>=0)
353 *fLog << "telescope " << fTelescopeIdx;
354 else
355 *fLog << "all telescopes";
356 *fLog << ")" << endl;
357
358 const Int_t lo = fTelescopeIdx<0 ? 0 : fTelescopeIdx;
359 const Int_t up = fTelescopeIdx<0 ? fNumTelescopes : fTelescopeIdx+1;
360
361 for (int i=lo; i<up; i++)
362 {
363 *fLog << all << "Telescope #" << setw(4) << i << " [X/Y/Z (R)]: ";
364 *fLog << setw(7) << fTelescopeX[i] << "/";
365 *fLog << setw(7) << fTelescopeY[i] << "/";
366 *fLog << setw(7) << fTelescopeZ[i] << " (R=" << setw(4) << fTelescopeR[i] << ")" << endl;
367 }
368
369 fNumTelescope = 0;
370
371 return kTRUE;
372}
373
374// --------------------------------------------------------------------------
375//
376// Reads the header of the next block. Depending on the current fReadState
377// and the block-type of the new header the methods decides wether
378// to continue the reading (kCONTINEU), to exit the Process() method (kTRUE)
379// or to try to read a new file (kFALSE).
380// Return codes:
381// - kFALSE : end of file was detected and no new header was read
382// - kTRU: A EventEnd was already found (fReadStatus == 3) and
383// the current header is not the RunEnd
384// - kCONTINUE: all other cases.
385Int_t MCorsikaRead::ReadNextBlockHeader()
386{
387 if (fInFormat->NextBlock(fReadState/* == 4*/, fBlockType, fBlockVersion,
388 fBlockIdentifier, fBlockLength) == kFALSE)
389 // end of file
390 return kFALSE;
391/*
392 gLog << "Next fBlock: type=" << fBlockType << " version=" << fBlockVersion;
393 gLog << " identifier=" << fBlockIdentifier << " length=" << fBlockLength;
394 gLog << " readState= " << fReadState << endl;
395*/
396 if (fReadState == 3 && fBlockType != 1210)
397 // fReadState == 3 means we have read the event end
398 // most of the time we can now save our data. BUT:
399 // bBlockType != 1210 means the next block is not RUNE,
400 // we want to read the RUNE block, before we
401 // save everything (last events and RUNE)
402 return kTRUE;
403
404 if (fBlockType == 1204 && (fReadState != 2 && fReadState != 15))
405 // next is a new set of telescope arrays, we store the previous ones
406 // but not if this is the first one (fReadState != 2)
407 return kTRUE;
408
409 return kCONTINUE;
410
411}
412
413// --------------------------------------------------------------------------
414//
415// The PreProcess of this task checks for the following containers in the
416// list:
417// MCorsikaRunHeader <output> if not found it is created
418// MCorsikaEvtHeader <output> if not found it is created
419// MCorsikaEvtData <output> if not found it is created
420// MCorsikaCrateArray <output> if not found it is created
421// MCorsikaEvtTime <output> if not found it is created (MTime)
422//
423// If all containers are found or created the run header is read from the
424// binary file and printed. If the Magic-Number (file identification)
425// doesn't match we stop the eventloop.
426//
427// Now the EvtHeader and EvtData containers are initialized.
428//
429Int_t MCorsikaRead::PreProcess(MParList *pList)
430{
431 //
432 // open the input stream
433 // first of all check if opening the file in the constructor was
434 // successfull
435 //
436 fParList = 0;
437
438 //
439 // check if all necessary containers exist in the Parameter list.
440 // if not create one and add them to the list
441 //
442 fRunHeader = (MCorsikaRunHeader*)pList->FindCreateObj("MCorsikaRunHeader");
443 if (!fRunHeader)
444 return kFALSE;
445
446 fEvtHeader = (MCorsikaEvtHeader*)pList->FindCreateObj("MCorsikaEvtHeader");
447 if (!fEvtHeader)
448 return kFALSE;
449
450 fEvent = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent");
451 if (!fEvent)
452 return kFALSE;
453
454 *fLog << inf << "Calculating number of total events..." << flush;
455 if (!CalcNumTotalEvents())
456 return kFALSE;
457 *fLog << inf << " " << fNumTotalEvents << " found." << endl;
458
459 fNumFile=0;
460 fNumEvents=0;
461
462 fParList = pList;
463
464 return kTRUE;
465}
466
467// --------------------------------------------------------------------------
468//
469// The Process reads one event from the binary file:
470// - The event header is read
471// - the run header is read (only once per file)
472// - the raw data information of one event is read
473//
474Int_t MCorsikaRead::Process()
475{
476 //fEvent->Clear();
477 //fEvent->Reset();
478 fEvent->Resize(0);
479 while (1) // loop for multiple input files
480 {
481 if (fInFormat)
482 {
483
484 while (1) // loop per input block
485 {
486
487 if (fReadState == 4)
488 {
489 fTopBlockLength -= fBlockLength + 12;
490 if (fTopBlockLength <= 0)
491 // all data of a top block are read, go back to normal state
492 fReadState = 15; // not used
493 }
494
495
496 Int_t status = kTRUE;
497 switch (fBlockType)
498 {
499 case 1200: // the run header
500 status = fRunHeader->ReadEvt(fInFormat);
501 fReadState = 1; // RUNH is read
502 break;
503
504 case 1201: // telescope position
505 status = ReadTelescopePosition();
506 break;
507
508 case 1202: // the event header
509 Float_t buffer[272];
510 if (!fInFormat->Read(buffer, 272 * sizeof(Float_t)))
511 return kFALSE;
512
513 if (fReadState == 1) // first event after RUN header
514 {
515 fRunHeader->ReadEventHeader(buffer);
516 fRunHeader->Print();
517 }
518
519 status = fEvtHeader->ReadEvt(buffer);
520 if (fArrayIdx >= (Int_t)fEvtHeader->GetTotReuse())
521 {
522 *fLog << err << "ERROR - Requested array index " << fArrayIdx <<
523 " exceeds number of arrays " << fEvtHeader->GetTotReuse() <<
524 " in file." << endl;
525 return kERROR;
526 }
527
528 fReadState = 2;
529 break;
530
531 case 1204: // top level block for one array (only for eventio data)
532 if (fArrayIdx < 0 || fArrayIdx == fBlockIdentifier)
533 {
534 fReadState = 4;
535 fTopBlockLength = fBlockLength;
536 }
537 else
538 // skip this array of telescopes
539 fInFormat->Seek(fBlockLength);
540
541 break;
542
543
544 case 1205: // eventio data
545 {
546 Int_t telIdx = fBlockIdentifier % 1000;
547 if (fBlockVersion == 0 &&
548 (fTelescopeIdx < 0 || fTelescopeIdx == telIdx) )
549 {
550 status = fEvent->ReadEventIoEvt(fInFormat);
551
552 Int_t arrayIdx = fBlockIdentifier / 1000;
553 Float_t xArrOff, yArrOff;
554 fEvtHeader->GetArrayOffset(arrayIdx, xArrOff, yArrOff);
555 fEvtHeader->SetTelescopeOffset(arrayIdx,
556 xArrOff + fTelescopeY[telIdx],
557 yArrOff - fTelescopeX[telIdx] );
558 fEvent->AddXY(xArrOff + fTelescopeY[telIdx],
559 yArrOff - fTelescopeX[telIdx]);
560 fEvent->SimWavelength(fRunHeader->GetWavelengthMin(),
561 fRunHeader->GetWavelengthMax());
562 }
563 else if (fBlockVersion == 1000 && (fTelescopeIdx < 0 || fTelescopeIdx == telIdx) )
564 {
565 status = fEvent->ReadEventIoEvtCompact(fInFormat);
566
567 Int_t arrayIdx = fBlockIdentifier / 1000;
568 Float_t xArrOff, yArrOff;
569 fEvtHeader->GetArrayOffset(arrayIdx, xArrOff, yArrOff);
570 fEvtHeader->SetTelescopeOffset(arrayIdx,
571 xArrOff + fTelescopeY[telIdx],
572 yArrOff - fTelescopeX[telIdx] );
573 fEvent->AddXY(xArrOff + fTelescopeY[telIdx],
574 yArrOff - fTelescopeX[telIdx]);
575 fEvent->SimWavelength(fRunHeader->GetWavelengthMin(),
576 fRunHeader->GetWavelengthMax());
577 }
578 else
579 // skip this telescope
580 fInFormat->Seek(fBlockLength);
581 }
582 break;
583
584 case 1209: // the event end
585 status = fEvtHeader->ReadEvtEnd(fInFormat);
586
587 if (fReadState == 10 || fReadState == 2)
588 {
589 // corsika events
590 fEvtHeader->ResetNumReuse();
591 fEvtHeader->InitXY();
592 fBlockType = 1109; // save corsika events
593 continue;
594 }
595
596 fReadState = 3; // event closed, run still open
597 break;
598
599
600 case 1210: // the run end
601 status = fRunHeader->ReadEvtEnd(fInFormat, kTRUE);
602 fNumEvents += fRunHeader->GetNumEvents();
603 //fRunHeader->SetReadyToSave();
604 fReadState = 5; // back to starting point
605 // this may be the last block in the file.
606 // We exit to write the header into the file
607 fBlockLength = 0;
608 fBlockType = 0; // not available type
609 return kTRUE;
610 break;
611
612
613 case 1105: // event block of raw format
614 if (fReadState == 2 || fReadState == 10)
615 {
616 Int_t oldSize = fRawEvemtBuffer.size();
617 fRawEvemtBuffer.resize(oldSize + fBlockLength / sizeof(Float_t));
618 status = fInFormat->Read(&fRawEvemtBuffer[oldSize], fBlockLength);
619 fReadState = 10;
620 }
621 else
622 fInFormat->Seek(fBlockLength);
623 break;
624
625
626 case 1109: // save corsika events
627 fEvtHeader->InitXY();
628 status = fEvent->ReadCorsikaEvt(&fRawEvemtBuffer[0],
629 fRawEvemtBuffer.size() / 7,
630 fEvtHeader->GetNumReuse()+1);
631 fEvtHeader->IncNumReuse();
632
633 if (fEvtHeader->GetNumReuse() == fRunHeader->GetNumReuse())
634 {
635 // this was the last reuse. Set fBlockType to EVTE to save
636 // it the next time.
637 fRawEvemtBuffer.resize(0);
638
639 fReadState = 3;
640 fBlockType = 1209;
641 }
642 else
643 // save this reuse
644 return status;
645
646 break;
647
648 default:
649 // unknown block, we skip it
650 fReadState = 15;
651 fInFormat->Seek(fBlockLength);
652
653 }
654
655 if (status != kTRUE)
656 // there was an error while data were read
657 return status;
658
659 Int_t headerStatus = ReadNextBlockHeader();
660 if (headerStatus == kFALSE)
661 // end of file
662 break;
663 if (headerStatus == kTRUE)
664 // we shall quit this loop
665 return kTRUE;
666
667 // else continue
668 }
669
670 }
671
672 //
673 // If an event could not be read from file try to open new file
674 //
675 const Int_t rc = OpenNextFile();
676 if (rc!=kTRUE)
677 return rc;
678
679 if (ReadNextBlockHeader() == kFALSE)
680 return kFALSE;
681 }
682 return kTRUE;
683}
684
685
686// --------------------------------------------------------------------------
687//
688// Close the file. Check whether the number of read events differs from
689// the number the file should containe (MCorsikaRunHeader). Prints a warning
690// if it doesn't match.
691//
692Int_t MCorsikaRead::PostProcess()
693{
694
695 const UInt_t n = fNumEvents*fRunHeader->GetNumReuse()*(fTelescopeIdx<0?fNumTelescopes:1);
696
697 //
698 // Sanity check for the number of events
699 //
700 if (n==GetNumExecutions()-1 || GetNumExecutions()==0)
701 return kTRUE;
702
703 *fLog << warn << dec;
704 *fLog << "Warning - number of read events (" << GetNumExecutions()-1;
705 *fLog << ") doesn't match number in run header(s) (";
706 *fLog << fRunHeader->GetNumReuse() << "*" << fNumEvents << "=" << n << ")." << endl;
707
708
709 return kTRUE;
710}
711
Note: See TracBrowser for help on using the repository browser.