source: branches/Corsika7405Compatibility/mcorsika/MCorsikaRead.cc@ 18795

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