source: trunk/Mars/mcorsika/MCorsikaRead.cc@ 10107

Last change on this file since 10107 was 10060, checked in by rohlfs, 14 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: 18.8 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
251Bool_t MCorsikaRead::CalcNumTotalEvents()
252{
253 fNumTotalEvents = 0;
254
255 Bool_t rc = kTRUE;
256
257 while (1)
258 {
259 switch (OpenNextFile(kFALSE))
260 {
261 case kFALSE:
262 break;
263 case kERROR:
264 rc = kFALSE;
265 break;
266 case kTRUE:
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())
297 {
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;
308 }
309
310 fNumTotalEvents += fRunHeader->GetNumEvents()*fRunHeader->GetNumReuse()*
311 (fTelescopeIdx<0 ? fNumTelescopes : 1);
312 continue;
313 }
314 break;
315 }
316
317 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;
363}
364
365// --------------------------------------------------------------------------
366//
367// The PreProcess of this task checks for the following containers in the
368// list:
369// MCorsikaRunHeader <output> if not found it is created
370// MCorsikaEvtHeader <output> if not found it is created
371// MCorsikaEvtData <output> if not found it is created
372// MCorsikaCrateArray <output> if not found it is created
373// MCorsikaEvtTime <output> if not found it is created (MTime)
374//
375// If all containers are found or created the run header is read from the
376// binary file and printed. If the Magic-Number (file identification)
377// doesn't match we stop the eventloop.
378//
379// Now the EvtHeader and EvtData containers are initialized.
380//
381Int_t MCorsikaRead::PreProcess(MParList *pList)
382{
383 //
384 // open the input stream
385 // first of all check if opening the file in the constructor was
386 // successfull
387 //
388 fParList = 0;
389
390 //
391 // check if all necessary containers exist in the Parameter list.
392 // if not create one and add them to the list
393 //
394 fRunHeader = (MCorsikaRunHeader*)pList->FindCreateObj("MCorsikaRunHeader");
395 if (!fRunHeader)
396 return kFALSE;
397
398 fEvtHeader = (MCorsikaEvtHeader*)pList->FindCreateObj("MCorsikaEvtHeader");
399 if (!fEvtHeader)
400 return kFALSE;
401
402 fEvent = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent");
403 if (!fEvent)
404 return kFALSE;
405
406 *fLog << inf << "Calculating number of total events..." << flush;
407 if (!CalcNumTotalEvents())
408 return kFALSE;
409 *fLog << inf << " " << fNumTotalEvents << " found." << endl;
410
411 fNumFile=0;
412 fNumEvents=0;
413
414 fParList = pList;
415
416 return kTRUE;
417}
418
419
420// --------------------------------------------------------------------------
421//
422// The Process reads one event from the binary file:
423// - The event header is read
424// - the run header is read
425// - all crate information is read
426// - the raw data information of one event is read
427//
428Int_t MCorsikaRead::Process()
429{
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;
611}
612
613// --------------------------------------------------------------------------
614//
615// Close the file. Check whether the number of read events differs from
616// the number the file should containe (MCorsikaRunHeader). Prints a warning
617// if it doesn't match.
618//
619Int_t MCorsikaRead::PostProcess()
620{
621
622 const UInt_t n = fNumEvents*fRunHeader->GetNumReuse()*(fTelescopeIdx<0?fNumTelescopes:1);
623
624 //
625 // Sanity check for the number of events
626 //
627 if (n==GetNumExecutions()-1 || GetNumExecutions()==0)
628 return kTRUE;
629
630 *fLog << warn << dec;
631 *fLog << "Warning - number of read events (" << GetNumExecutions()-2;
632 *fLog << ") doesn't match number in run header(s) (";
633 *fLog << fRunHeader->GetNumReuse() << "*" << fNumEvents << "=" << n << ")." << endl;
634
635 return kTRUE;
636}
637
Note: See TracBrowser for help on using the repository browser.