Changeset 10147


Ignore:
Timestamp:
Feb 9, 2011, 11:07:20 AM (9 years ago)
Author:
rohlfs
Message:
Redesign of the Process() - method.
Location:
trunk/Mars/mcorsika
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/Mars/mcorsika/MCorsikaRead.cc

    r10060 r10147  
    249249}
    250250
     251// --------------------------------------------------------------------------
     252//
     253// Calculates the total number of events from all input files and stores
     254// the number in fNumTotalEvents.                                       
     255//
    251256Bool_t MCorsikaRead::CalcNumTotalEvents()
    252257{
     
    318323}
    319324
     325// --------------------------------------------------------------------------
     326//
     327// Reads the the position of all telescopes in one array
     328//
    320329Int_t MCorsikaRead::ReadTelescopePosition()
    321330{
     
    352361   for (int i=lo; i<up; i++)
    353362   {
    354       *fLog << inf << "Telescope #" << setw(4) << i << " [X/Y/Z (R)]: ";
     363      *fLog << all << "Telescope #" << setw(4) << i << " [X/Y/Z (R)]: ";
    355364      *fLog << setw(7) << fTelescopeX[i] << "/";
    356365      *fLog << setw(7) << fTelescopeY[i] << "/";
     
    361370
    362371   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 
    363411}
    364412
     
    417465}
    418466
     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}
    419666
    420667// --------------------------------------------------------------------------
     
    426673//  - the raw data information of one event is read
    427674//
    428 Int_t MCorsikaRead::Process()
     675Int_t MCorsikaRead::Process_Old()
    429676{
    430677
     
    453700      if (fInFormat)
    454701      {
    455          Int_t blockType, blockVersion, blockIdentifier, blockLength;
    456 
    457          while (fInFormat->NextBlock(fReadState == 4, blockType, blockVersion,
    458                                      blockIdentifier, blockLength))
     702
     703         while (fInFormat->NextBlock(fReadState == 4, fBlockType, fBlockVersion,
     704                                     fBlockIdentifier, fBlockLength))
    459705         {
    460 //            gLog << "Next block:  type=" << blockType << " version=" << blockVersion;
    461 //            gLog << " identifier=" << blockIdentifier << " length=" << blockLength << endl;
     706            gLog << "Next fBlock:  type=" << fBlockType << " version=" << fBlockVersion;
     707            gLog << " identifier=" << fBlockIdentifier << " length=" << fBlockLength << endl;
    462708
    463709
    464710            if (fReadState == 4)
    465711               {
    466                fTopBlockLength -= blockLength + 12;
     712               fTopBlockLength -= fBlockLength + 12;
    467713               if (fTopBlockLength <= 0)
    468714                  // all data of a top block are read, go back to normal state
     
    471717
    472718            Int_t status = kTRUE;
    473             switch (blockType)
     719            switch (fBlockType)
    474720               {
    475721               case 1200:       // the run header
     
    480726               case 1201:       // telescope position
    481727                  status = ReadTelescopePosition();
     728sleep(5);
    482729                  break;
    483730
     
    507754
    508755               case 1204:
    509                   if (fArrayIdx < 0 || fArrayIdx == blockIdentifier)
     756                  if (fArrayIdx < 0 || fArrayIdx == fBlockIdentifier)
    510757                     {
    511758                     fReadState = 4;
    512                      fTopBlockLength = blockLength;
     759                     fTopBlockLength = fBlockLength;
    513760                     }
    514761                  else
    515762                     // skip this array of telescopes
    516                      fInFormat->Seek(blockLength);
     763                     fInFormat->Seek(fBlockLength);
    517764
    518765                  break;
     
    520767               case 1205:
    521768                  {
    522                   Int_t telIdx   = blockIdentifier % 1000;
    523                   if (blockVersion == 0                               &&
    524                       (fTelescopeIdx < 0 || fTelescopeIdx ==  telIdx) &&
    525                        blockLength > 12)
     769                  Int_t telIdx   = fBlockIdentifier % 1000;
     770                  if (fBlockVersion == 0                               &&
     771                      (fTelescopeIdx < 0 || fTelescopeIdx ==  telIdx)     )
    526772                     {
    527773                     status = fEvent->ReadEventIoEvt(fInFormat);
    528774
    529                      Int_t arrayIdx = blockIdentifier / 1000;
     775                     Int_t arrayIdx = fBlockIdentifier / 1000;
    530776                     Float_t xArrOff, yArrOff;
    531777                     fEvtHeader->GetArrayOffset(arrayIdx, xArrOff, yArrOff);
     
    545791                  else
    546792                     // skip this telescope
    547                      fInFormat->Seek(blockLength);
     793                     fInFormat->Seek(fBlockLength);
    548794                  }
    549795                  break;
     
    551797               case 1209:  // the event end
    552798                  status = fEvtHeader->ReadEvtEnd(fInFormat);
    553                   if (fReadState == 10)
     799                  if (fReadState == 10 || fReadState == 2)
    554800                     {
    555801                     // all particles of this event are read, now save them
     
    580826                     {
    581827                     Int_t oldSize = fRawEvemtBuffer.size();
    582                      fRawEvemtBuffer.resize(oldSize + blockLength / sizeof(Float_t));
    583                      status = fInFormat->Read(&fRawEvemtBuffer[oldSize], blockLength);
     828                     fRawEvemtBuffer.resize(oldSize + fBlockLength / sizeof(Float_t));
     829                     status = fInFormat->Read(&fRawEvemtBuffer[oldSize], fBlockLength);
    584830                     fReadState = 10;
    585831                     }
    586832                  else
    587                      fInFormat->Seek(blockLength);
     833                     fInFormat->Seek(fBlockLength);
    588834                  break;
    589835
    590836               default:
    591837                  // unknown block, we skip it
    592                   fInFormat->Seek(blockLength);
     838                  fInFormat->Seek(fBlockLength);
    593839
    594840               }
     
    629875
    630876    *fLog << warn << dec;
    631     *fLog << "Warning - number of read events (" << GetNumExecutions()-2;
     877    *fLog << "Warning - number of read events (" << GetNumExecutions()-1;
    632878    *fLog << ") doesn't match number in run header(s) (";
    633879    *fLog << fRunHeader->GetNumReuse() << "*" << fNumEvents << "=" << n << ")." << endl;
    634880
     881
    635882    return kTRUE;
    636883}
  • trunk/Mars/mcorsika/MCorsikaRead.h

    r10060 r10147  
    4848    TArrayF fTelescopeR;       //! Radii of spheres around tel. (unit: cm)
    4949
     50    Int_t   fBlockType;        //! block type of the just read block header
     51    Int_t   fBlockVersion;     //! block version of the just read block header
     52    Int_t   fBlockIdentifier;  //! block identifier of the just read block header
     53    Int_t   fBlockLength;      //! block length from the just read block header
     54
    5055    Int_t   fReadState;        //! 0: nothing read yet
    5156                               //  1: RUNH is already read
     
    5560                               //  5: RUNE is already read
    5661                               // 10: raw events are currently read
    57                                // 11: raw events are currently saved
     62                               // 15: undefined status
    5863    Int_t   fTopBlockLength;   // remaining length of the current top-level block 1204
    5964
     
    6671    Bool_t CalcNumTotalEvents();
    6772    Int_t  ReadTelescopePosition();
     73    Int_t  ReadNextBlockHeader();
    6874
    6975    // MTask
    7076    Int_t PreProcess(MParList *pList);
    7177    Int_t Process();
     78    Int_t Process_Old();
    7279    Int_t PostProcess();
    7380
Note: See TracChangeset for help on using the changeset viewer.