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

Last change on this file since 10009 was 9943, checked in by tbretz, 14 years ago
Implemented a working version of eventio which reads all telescope's data and all array's data.
File size: 15.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): Thomas Bretz 11/2008 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: Software Development, 2000-2008
21!
22!
23\* ======================================================================== */
24
25//////////////////////////////////////////////////////////////////////////////
26//
27// MCorsikaRead
28//
29// Input Containers:
30// -/-
31//
32// Output Containers:
33// MCorsikaRunHeader
34// MCorsikaEvtHeader
35// MPhotonEvent
36//
37//////////////////////////////////////////////////////////////////////////////
38#include "MCorsikaRead.h"
39
40#include <errno.h>
41#include <fstream>
42
43#include <TSystem.h>
44
45#include "MLog.h"
46#include "MLogManip.h"
47
48#include "MParList.h"
49#include "MStatusDisplay.h"
50
51#include "MCorsikaFormat.h"
52#include "MCorsikaRunHeader.h"
53#include "MCorsikaEvtHeader.h"
54
55#include "MPhotonEvent.h"
56
57ClassImp(MCorsikaRead);
58
59using namespace std;
60
61/* ----------- please don't delete and don't care about (Thomas) ------------
62#define kBUFSZ 64 //1024*1024*64
63#include <iomanip.h>
64class bifstream : public istream, public streambuf
65{
66private:
67 char fBuffer[kBUFSZ]; //!
68 FILE *fd;
69
70 int sync()
71 {
72 memset(fBuffer, 0, kBUFSZ);
73 return 0;
74 }
75 int underflow()
76 {
77 int sz=fread(fBuffer, kBUFSZ, 1, fd);
78 //int sz=fread(fBuffer, 1, kBUFSZ, fd);
79 setg(fBuffer, fBuffer, fBuffer+kBUFSZ);
80
81 return sz==1 ? *(unsigned char*)fBuffer : EOF;//EOF;
82 //return sz==kBUFSZ ? *(unsigned char*)fBuffer : EOF;//EOF;
83 }
84public:
85 bifstream(const char *name) : istream(this)
86 {
87 fd = fopen(name, "rb");
88 setbuf(fBuffer, kBUFSZ);
89 }
90};
91*/
92
93// --------------------------------------------------------------------------
94//
95// Default constructor. It tries to open the given file.
96//
97MCorsikaRead::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)
101{
102 fName = name ? name : "MRead";
103 fTitle = title ? title : "Read task to read DAQ binary files";
104
105 fFileNames = new TList;
106 fFileNames->SetOwner();
107
108 if (fname!=NULL)
109 AddFile(fname);
110}
111
112// --------------------------------------------------------------------------
113//
114// Destructor. Delete input stream.
115//
116MCorsikaRead::~MCorsikaRead()
117{
118 delete fFileNames;
119// if (fIn)
120// delete fIn;
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
164Bool_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
183// --------------------------------------------------------------------------
184//
185// This opens the next file in the list and deletes its name from the list.
186//
187Int_t MCorsikaRead::OpenNextFile(Bool_t print)
188{
189
190 //
191 // open the input stream and check if it is really open (file exists?)
192 //
193/* if (fIn)
194 delete fIn;
195 fIn = NULL;
196*/
197 if (fInFormat)
198 delete fInFormat;
199 fInFormat = NULL;
200
201 //
202 // Check for the existance of a next file to read
203 //
204 TObject *file = fFileNames->At(fNumFile);
205 if (!file)
206 return kFALSE;
207
208 //
209 // open the file which is the first one in the chain
210 //
211 const char *name = file->GetName();
212
213 const char *expname = gSystem->ExpandPathName(name);
214 fInFormat = MCorsikaFormat::CorsikaFormatFactory(expname);
215 delete [] expname;
216
217 if (fInFormat == NULL)
218 return kERROR;
219
220 *fLog << inf << "Open file: '" << name << "'" << endl;
221
222 if (fDisplay)
223 {
224 // Show the number of the last event after
225 // which we now open a new file
226 TString txt = GetFileName();
227 txt += " @ ";
228 txt += GetNumExecutions()-1;
229 fDisplay->SetStatusLine2(txt);
230 }
231
232 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();
312
313 if (!fParList)
314 return kTRUE;
315
316 //
317 // Search for MTaskList
318 //
319 MTask *tlist = (MTask*)fParList->FindObject("MTaskList");
320 if (!tlist)
321 {
322 *fLog << err << dbginf << "MTaskList not found... abort." << endl;
323 return kERROR;
324 }
325
326 //
327 // A new file has been opened and new headers have been read.
328 // --> ReInit tasklist
329 //
330 return tlist->ReInit(fParList) ? kTRUE : kERROR;
331}
332
333// --------------------------------------------------------------------------
334//
335// Return file name of current file.
336//
337TString MCorsikaRead::GetFullFileName() const
338{
339 const TObject *file = fFileNames->At(fNumFile-1);
340 return file ? file->GetName() : "";
341}
342
343// --------------------------------------------------------------------------
344//
345// Restart with the first file
346//
347Bool_t MCorsikaRead::Rewind()
348{
349 fNumFile=0;
350 fNumEvents=0;
351 return OpenNextFile()==kTRUE;
352}
353
354Bool_t MCorsikaRead::CalcNumTotalEvents()
355{
356 fNumTotalEvents = 0;
357
358 Bool_t rc = kTRUE;
359
360 while (1)
361 {
362 switch (OpenNextFile(kFALSE))
363 {
364 case kFALSE:
365 break;
366 case kERROR:
367 rc = kFALSE;
368 break;
369 case kTRUE:
370 if (!ReadEvtEnd())
371 {
372 rc = kFALSE;
373 break;
374 }
375
376 fNumTotalEvents += fRunHeader->GetNumEvents()*fRunHeader->GetNumReuse()*
377 (fTelescopeIdx<0 ? fNumTelescopes : 1);
378 continue;
379 }
380 break;
381 }
382/*
383 if (fIn)
384 delete fIn;
385 fIn = NULL;
386*/
387 return rc;
388}
389
390// --------------------------------------------------------------------------
391//
392// The PreProcess of this task checks for the following containers in the
393// list:
394// MCorsikaRunHeader <output> if not found it is created
395// MCorsikaEvtHeader <output> if not found it is created
396// MCorsikaEvtData <output> if not found it is created
397// MCorsikaCrateArray <output> if not found it is created
398// MCorsikaEvtTime <output> if not found it is created (MTime)
399//
400// If all containers are found or created the run header is read from the
401// binary file and printed. If the Magic-Number (file identification)
402// doesn't match we stop the eventloop.
403//
404// Now the EvtHeader and EvtData containers are initialized.
405//
406Int_t MCorsikaRead::PreProcess(MParList *pList)
407{
408 //
409 // open the input stream
410 // first of all check if opening the file in the constructor was
411 // successfull
412 //
413 fParList = 0;
414/*
415 if (!OpenStream())
416 return kFALSE;
417*/
418 //
419 // check if all necessary containers exist in the Parameter list.
420 // if not create one and add them to the list
421 //
422 fRunHeader = (MCorsikaRunHeader*)pList->FindCreateObj("MCorsikaRunHeader");
423 if (!fRunHeader)
424 return kFALSE;
425
426 fEvtHeader = (MCorsikaEvtHeader*)pList->FindCreateObj("MCorsikaEvtHeader");
427 if (!fEvtHeader)
428 return kFALSE;
429
430 fEvent = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent");
431 if (!fEvent)
432 return kFALSE;
433
434 *fLog << inf << "Calculating number of total events..." << flush;
435 if (!CalcNumTotalEvents())
436 return kFALSE;
437 *fLog << inf << " " << fNumTotalEvents << " found." << endl;
438
439 fNumFile=0;
440 fNumEvents=0;
441
442 fParList = pList;
443
444 return kTRUE;
445}
446
447// --------------------------------------------------------------------------
448//
449// Read a single event from the stream
450//
451Int_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}
521
522// --------------------------------------------------------------------------
523//
524// The Process reads one event from the binary file:
525// - The event header is read
526// - the run header is read
527// - all crate information is read
528// - the raw data information of one event is read
529//
530Int_t MCorsikaRead::Process()
531{
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;
558}
559
560// --------------------------------------------------------------------------
561//
562// Close the file. Check whether the number of read events differs from
563// the number the file should containe (MCorsikaRunHeader). Prints a warning
564// if it doesn't match.
565//
566Int_t MCorsikaRead::PostProcess()
567{
568 const UInt_t n = fNumEvents*fRunHeader->GetNumReuse()*(fTelescopeIdx<0?fNumTelescopes:1);
569
570 //
571 // Sanity check for the number of events
572 //
573 if (n==GetNumExecutions()-1 || GetNumExecutions()==0)
574 return kTRUE;
575
576 *fLog << warn << dec;
577 *fLog << "Warning - number of read events (" << GetNumExecutions()-1;
578 *fLog << ") doesn't match number in run header(s) (";
579 *fLog << fRunHeader->GetNumReuse() << "*" << fNumEvents << "=" << n << ")." << endl;
580
581 return kTRUE;
582}
583
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//
591Int_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}
Note: See TracBrowser for help on using the repository browser.