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

Last change on this file since 10147 was 10147, checked in by rohlfs, 9 years ago
Redesign of the Process() - method.
File size: 27.3 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)
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   while (1)  // loop for multiple input files
477   {
478      if (fInFormat)
479      {
480
481      while (1)  // loop per input block
482         {
483
484         if (fReadState == 4)
485            {
486            fTopBlockLength -= fBlockLength + 12;
487            if (fTopBlockLength <= 0)
488               // all data of a top block are read, go back to normal state
489               fReadState = 15; // not used
490            }
491
492
493         Int_t status = kTRUE;
494         switch (fBlockType)
495            {
496            case 1200:       // the run header
497               status = fRunHeader->ReadEvt(fInFormat);
498               fReadState = 1;  // RUNH is read
499               break;
500
501            case 1201:       // telescope position
502               status = ReadTelescopePosition();
503               break;
504
505            case 1202:  // the event header
506               Float_t buffer[272];
507               if (!fInFormat->Read(buffer, 272 * sizeof(Float_t))) 
508                  return kFALSE;
509
510               if (fReadState == 1)  // first event after RUN header
511                  {
512                  fRunHeader->ReadEventHeader(buffer);
513                  fRunHeader->Print();
514                  }
515
516               status = fEvtHeader->ReadEvt(buffer);
517               if (fArrayIdx >= (Int_t)fEvtHeader->GetTotReuse())
518                  {
519                  *fLog << err << "ERROR - Requested array index " << fArrayIdx << 
520                        " exceeds number of arrays " << fEvtHeader->GetTotReuse() << 
521                        " in file." << endl;
522                  return kERROR;
523                  }
524
525               fReadState = 2;
526               break;
527
528            case 1204: // top level block for one array (only for eventio data)
529               if (fArrayIdx < 0 || fArrayIdx == fBlockIdentifier)
530                  { 
531                  fReadState = 4;
532                  fTopBlockLength = fBlockLength;
533                  }
534               else
535                  // skip this array of telescopes
536                  fInFormat->Seek(fBlockLength);
537
538               break;
539
540
541            case 1205: // eventio data
542               {
543               Int_t telIdx   = fBlockIdentifier % 1000;
544               if (fBlockVersion == 0                               &&
545                     (fTelescopeIdx < 0 || fTelescopeIdx ==  telIdx)     )
546                  {
547                  status = fEvent->ReadEventIoEvt(fInFormat);
548
549                  Int_t arrayIdx = fBlockIdentifier / 1000;
550                  Float_t xArrOff, yArrOff;
551                  fEvtHeader->GetArrayOffset(arrayIdx, xArrOff, yArrOff);
552                  fEvtHeader->SetTelescopeOffset(arrayIdx, 
553                                                   xArrOff + fTelescopeY[telIdx], 
554                                                   yArrOff - fTelescopeX[telIdx] );
555                  fEvent->AddXY(xArrOff + fTelescopeY[telIdx], 
556                                 yArrOff - fTelescopeX[telIdx]);
557                  fEvent->SimWavelength(fRunHeader->GetWavelengthMin(), 
558                                          fRunHeader->GetWavelengthMax());
559                  }
560               else
561                  // skip this telescope
562                  fInFormat->Seek(fBlockLength);
563               }
564               break;
565
566            case 1209:  // the event end
567               status = fEvtHeader->ReadEvtEnd(fInFormat);
568               
569               if (fReadState == 10 || fReadState == 2)
570                  {
571                  // corsika events
572                  fEvtHeader->ResetNumReuse();
573                  fEvtHeader->InitXY();
574                  fBlockType = 1109;  // save corsika events
575                  continue;
576                  }
577
578               fReadState = 3;  // event closed, run still open
579               break;
580
581
582            case 1210:  // the run end
583               status = fRunHeader->ReadEvtEnd(fInFormat, kTRUE);
584               fNumEvents += fRunHeader->GetNumEvents();
585               fRunHeader->SetReadyToSave();
586               fReadState = 5;           // back to  starting point
587               // this may be the last block in the file.
588               // We exit to write the header into the file
589               fBlockLength = 0;
590               fBlockType   = 0; // not available type
591               return kTRUE;
592               break;
593
594
595            case 1105:  // event block of raw format
596               if (fReadState == 2 || fReadState == 10)
597                  {
598                  Int_t oldSize = fRawEvemtBuffer.size();
599                  fRawEvemtBuffer.resize(oldSize + fBlockLength / sizeof(Float_t));
600                  status = fInFormat->Read(&fRawEvemtBuffer[oldSize], fBlockLength); 
601                  fReadState = 10;
602                  }
603               else
604                  fInFormat->Seek(fBlockLength);
605               break;
606
607
608            case 1109:  // save corsika events
609               fEvtHeader->InitXY();
610               status = fEvent->ReadCorsikaEvt(&fRawEvemtBuffer[0], 
611                                               fRawEvemtBuffer.size() / 7, 
612                                               fEvtHeader->GetNumReuse()+1);
613               fEvtHeader->IncNumReuse();
614
615               if (fEvtHeader->GetNumReuse() == fRunHeader->GetNumReuse())
616                  {
617                  // this was the last reuse. Set fBlockType to EVTE to save
618                  // it the next time.
619                  fRawEvemtBuffer.resize(0);
620
621                  fReadState = 3;
622                  fBlockType = 1209;
623                  }
624               else
625                  // save this reuse
626                  return status;
627
628               break;               
629
630            default:
631               // unknown block, we skip it
632               fReadState = 15;
633               fInFormat->Seek(fBlockLength);
634
635            }
636
637         if (status != kTRUE)
638            // there was an error while data were read
639            return status;
640
641         Int_t headerStatus =  ReadNextBlockHeader();
642         if (headerStatus == kFALSE) 
643            // end of file
644            break;
645         if (headerStatus == kTRUE)
646            // we shall quit this loop
647            return kTRUE;
648
649            // else continue
650         }
651
652      }
653
654      //
655      // If an event could not be read from file try to open new file
656      //
657      const Int_t rc = OpenNextFile();
658      if (rc!=kTRUE)
659          return rc;
660
661      if (ReadNextBlockHeader() == kFALSE)
662         return kFALSE;
663   }
664   return kTRUE;
665}
666
667// --------------------------------------------------------------------------
668//
669// The Process reads one event from the binary file:
670//  - The event header is read
671//  - the run header is read
672//  - all crate information is read
673//  - the raw data information of one event is read
674//
675Int_t MCorsikaRead::Process_Old()
676{
677
678   if (fReadState == 11)
679      {
680      // we are currently saving the events of the raw format in the root file
681      if (fEvtHeader->GetNumReuse() == fRunHeader->GetNumReuse())
682         {
683         // all data are saved
684         fRawEvemtBuffer.resize(0);
685         fReadState = 3;
686         }
687      else
688         {
689         fEvtHeader->InitXY();
690         Int_t rc = fEvent->ReadCorsikaEvt(&fRawEvemtBuffer[0], 
691                                           fRawEvemtBuffer.size() / 7, 
692                                           fEvtHeader->GetNumReuse()+1);
693         fEvtHeader->IncNumReuse();
694         return rc;
695         }
696      }
697
698   while (1)
699   {
700      if (fInFormat)
701      {
702
703         while (fInFormat->NextBlock(fReadState == 4, fBlockType, fBlockVersion, 
704                                     fBlockIdentifier, fBlockLength))
705         {
706            gLog << "Next fBlock:  type=" << fBlockType << " version=" << fBlockVersion;
707            gLog << " identifier=" << fBlockIdentifier << " length=" << fBlockLength << endl;
708
709
710            if (fReadState == 4)
711               {
712               fTopBlockLength -= fBlockLength + 12;
713               if (fTopBlockLength <= 0)
714                  // all data of a top block are read, go back to normal state
715                  fReadState = 3;
716               }
717
718            Int_t status = kTRUE;
719            switch (fBlockType)
720               {
721               case 1200:       // the run header
722                  status = fRunHeader->ReadEvt(fInFormat);
723                  fReadState = 1;  // RUNH is read
724                  break;
725
726               case 1201:       // telescope position
727                  status = ReadTelescopePosition();
728sleep(5);
729                  break;
730
731               case 1202:  // the event header
732                  Float_t buffer[272];
733                  if (!fInFormat->Read(buffer, 272 * sizeof(Float_t))) 
734                     return kFALSE;
735
736                  if (fReadState == 1)  // first event after RUN header
737                     {
738                     fRunHeader->ReadEventHeader(buffer);
739                     fRunHeader->Print();
740                     }
741
742                  status = fEvtHeader->ReadEvt(buffer);
743                  if (fArrayIdx >= (Int_t)fEvtHeader->GetTotReuse())
744                     {
745                     *fLog << err << "ERROR - Requested array index " << fArrayIdx << 
746                           " exceeds number of arrays " << fEvtHeader->GetTotReuse() << 
747                           " in file." << endl;
748                     return kERROR;
749                     }
750
751
752                  fReadState = 2;
753                  break;
754
755               case 1204:
756                  if (fArrayIdx < 0 || fArrayIdx == fBlockIdentifier)
757                     { 
758                     fReadState = 4;
759                     fTopBlockLength = fBlockLength;
760                     }
761                  else
762                     // skip this array of telescopes
763                     fInFormat->Seek(fBlockLength);
764
765                  break;
766               
767               case 1205:
768                  {
769                  Int_t telIdx   = fBlockIdentifier % 1000;
770                  if (fBlockVersion == 0                               &&
771                      (fTelescopeIdx < 0 || fTelescopeIdx ==  telIdx)     )
772                     {
773                     status = fEvent->ReadEventIoEvt(fInFormat);
774
775                     Int_t arrayIdx = fBlockIdentifier / 1000;
776                     Float_t xArrOff, yArrOff;
777                     fEvtHeader->GetArrayOffset(arrayIdx, xArrOff, yArrOff);
778                     fEvtHeader->SetTelescopeOffset(arrayIdx, 
779                                                    xArrOff + fTelescopeY[telIdx], 
780                                                    yArrOff - fTelescopeX[telIdx] );
781                     fEvent->AddXY(xArrOff + fTelescopeY[telIdx], 
782                                   yArrOff - fTelescopeX[telIdx]);
783                     fEvent->SimWavelength(fRunHeader->GetWavelengthMin(), 
784                                           fRunHeader->GetWavelengthMax());
785
786                     if (status == kTRUE)
787                        // end of reading for one telescope in one array ==
788                        // end of this Process - step
789                        return status;
790                     }
791                  else
792                     // skip this telescope
793                     fInFormat->Seek(fBlockLength);
794                  }
795                  break;
796
797               case 1209:  // the event end
798                  status = fEvtHeader->ReadEvtEnd(fInFormat);
799                  if (fReadState == 10 || fReadState == 2)
800                     {
801                     // all particles of this event are read, now save them
802                     fReadState = 11;
803                     fEvtHeader->ResetNumReuse();
804
805                     fEvtHeader->InitXY();
806                     Int_t rc = fEvent->ReadCorsikaEvt(&fRawEvemtBuffer[0], 
807                                                       fRawEvemtBuffer.size() / 7, 
808                                                       fEvtHeader->GetNumReuse()+1);
809                     fEvtHeader->IncNumReuse();
810                     return rc;
811                     }
812                  else
813                     fReadState = 3;  // event closed, run still open
814                  break;
815
816               case 1210:  // the run end
817                  status = fRunHeader->ReadEvtEnd(fInFormat, kTRUE);
818                  fNumEvents += fRunHeader->GetNumEvents();
819                  fRunHeader->SetReadyToSave();
820                  fReadState = 5;           // back to  starting point
821                  return status;
822                  break;
823
824               case 1105:  // event block of raw format
825                  if (fReadState == 2 || fReadState == 10)
826                     {
827                     Int_t oldSize = fRawEvemtBuffer.size();
828                     fRawEvemtBuffer.resize(oldSize + fBlockLength / sizeof(Float_t));
829                     status = fInFormat->Read(&fRawEvemtBuffer[oldSize], fBlockLength); 
830                     fReadState = 10;
831                     }
832                  else
833                     fInFormat->Seek(fBlockLength);
834                  break;
835
836               default:
837                  // unknown block, we skip it
838                  fInFormat->Seek(fBlockLength);
839
840               }
841
842            if (status != kTRUE)
843                return status;
844           
845         }
846
847      }
848
849      //
850      // If an event could not be read from file try to open new file
851      //
852      const Int_t rc = OpenNextFile();
853      if (rc!=kTRUE)
854          return rc;
855   }
856   return kTRUE;
857}
858
859// --------------------------------------------------------------------------
860//
861//  Close the file. Check whether the number of read events differs from
862//  the number the file should containe (MCorsikaRunHeader). Prints a warning
863//  if it doesn't match.
864//
865Int_t MCorsikaRead::PostProcess()
866{
867
868    const UInt_t n = fNumEvents*fRunHeader->GetNumReuse()*(fTelescopeIdx<0?fNumTelescopes:1);
869
870    //
871    // Sanity check for the number of events
872    //
873    if (n==GetNumExecutions()-1 || GetNumExecutions()==0)
874        return kTRUE;
875
876    *fLog << warn << dec;
877    *fLog << "Warning - number of read events (" << GetNumExecutions()-1;
878    *fLog << ") doesn't match number in run header(s) (";
879    *fLog << fRunHeader->GetNumReuse() << "*" << fNumEvents << "=" << n << ")." << endl;
880
881
882    return kTRUE;
883}
884
Note: See TracBrowser for help on using the repository browser.