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

Last change on this file since 10060 was 10060, checked in by rohlfs, 10 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.