Changeset 10147 for trunk/Mars/mcorsika

Timestamp:

02/09/11 11:07:20 (14 years ago)

Author:

rohlfs

Message:

Redesign of the Process() - method.

Location:

trunk/Mars/mcorsika

Files:

• MCorsikaRead.cc (modified) (15 diffs)
• MCorsikaRead.h (modified) (3 diffs)

trunk/Mars/mcorsika/MCorsikaRead.cc

@@ -249 +249 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Calculates the total number of events from all input files and stores
+// the number in fNumTotalEvents.                                       
+//
 Bool_t MCorsikaRead::CalcNumTotalEvents()
 {
@@ -318 +323 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Reads the the position of all telescopes in one array
+//
 Int_t MCorsikaRead::ReadTelescopePosition()
 {
@@ -352 +361 @@
    for (int i=lo; i<up; i++)
    {
-      *fLog << inf << "Telescope #" << setw(4) << i << " [X/Y/Z (R)]: ";
+      *fLog << all << "Telescope #" << setw(4) << i << " [X/Y/Z (R)]: ";
       *fLog << setw(7) << fTelescopeX[i] << "/";
       *fLog << setw(7) << fTelescopeY[i] << "/";
@@ -361 +370 @@
 
    return kTRUE;
+}
+
+// --------------------------------------------------------------------------
+//
+// Reads the header of the next block. Depending on the current fReadState
+// and the block-type of the new header the methods decides wether 
+// to continue the reading (kCONTINEU), to exit the Process() method (kTRUE)
+// or to try to read a new file (kFALSE).
+// Return codes:
+//   - kFALSE :   end of file was detected and no new header was read
+//   - kTRU:      A EventEnd was already found (fReadStatus == 3) and
+//                the current header is not the RunEnd
+//   - kCONTINUE: all other cases.
+Int_t MCorsikaRead::ReadNextBlockHeader()
+{
+   if (fInFormat->NextBlock(fReadState == 4, fBlockType, fBlockVersion, 
+                            fBlockIdentifier, fBlockLength)               == kFALSE)
+      // end of file
+      return kFALSE;
+
+    gLog << "Next fBlock:  type=" << fBlockType << " version=" << fBlockVersion;
+    gLog << " identifier=" << fBlockIdentifier << " length=" << fBlockLength; 
+    gLog << " readState= " << fReadState << endl;
+   
+   if (fReadState == 3 && fBlockType != 1210)
+      // fReadState == 3    means we have read the event end 
+      //                    most of the time we can now save our data. BUT:
+      // bBlockType != 1210 means the next block is not RUNE, 
+      //                    we want to read the RUNE block, before we 
+      //                    save everything (last events and RUNE) 
+      return kTRUE;
+
+   if (fBlockType == 1204 && fReadState != 2)
+      // next is a new set of telescope arrays, we store the previous ones
+      // but not if this is the first one (fReadState != 2)
+      return kTRUE;
+
+   return kCONTINUE;
+  
 }
 
@@ -417 +465 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// The Process reads one event from the binary file:
+//  - The event header is read
+//  - the run header is read (only once per file)
+//  - the raw data information of one event is read
+//
+Int_t MCorsikaRead::Process()
+{
+   while (1)  // loop for multiple input files
+   {
+      if (fInFormat)
+      {
+
+      while (1)  // loop per input block
+         {
+
+         if (fReadState == 4)
+            {
+            fTopBlockLength -= fBlockLength + 12;
+            if (fTopBlockLength <= 0)
+               // all data of a top block are read, go back to normal state
+               fReadState = 15; // not used
+            }
+
+
+         Int_t status = kTRUE;
+         switch (fBlockType)
+            {
+            case 1200:       // the run header
+               status = fRunHeader->ReadEvt(fInFormat);
+               fReadState = 1;  // RUNH is read
+               break;
+
+            case 1201:       // telescope position
+               status = ReadTelescopePosition();
+               break;
+
+            case 1202:  // the event header
+               Float_t buffer[272];
+               if (!fInFormat->Read(buffer, 272 * sizeof(Float_t))) 
+                  return kFALSE;
+
+               if (fReadState == 1)  // first event after RUN header
+                  {
+                  fRunHeader->ReadEventHeader(buffer);
+                  fRunHeader->Print();
+                  }
+
+               status = fEvtHeader->ReadEvt(buffer);
+               if (fArrayIdx >= (Int_t)fEvtHeader->GetTotReuse())
+                  {
+                  *fLog << err << "ERROR - Requested array index " << fArrayIdx << 
+                        " exceeds number of arrays " << fEvtHeader->GetTotReuse() << 
+                        " in file." << endl;
+                  return kERROR;
+                  }
+
+               fReadState = 2;
+               break;
+
+            case 1204: // top level block for one array (only for eventio data)
+               if (fArrayIdx < 0 || fArrayIdx == fBlockIdentifier)
+                  { 
+                  fReadState = 4;
+                  fTopBlockLength = fBlockLength;
+                  }
+               else
+                  // skip this array of telescopes
+                  fInFormat->Seek(fBlockLength);
+
+               break;
+
+
+            case 1205: // eventio data
+               {
+               Int_t telIdx   = fBlockIdentifier % 1000;
+               if (fBlockVersion == 0                               &&
+                     (fTelescopeIdx < 0 || fTelescopeIdx ==  telIdx)     )
+                  {
+                  status = fEvent->ReadEventIoEvt(fInFormat);
+
+                  Int_t arrayIdx = fBlockIdentifier / 1000;
+                  Float_t xArrOff, yArrOff;
+                  fEvtHeader->GetArrayOffset(arrayIdx, xArrOff, yArrOff);
+                  fEvtHeader->SetTelescopeOffset(arrayIdx, 
+                                                   xArrOff + fTelescopeY[telIdx], 
+                                                   yArrOff - fTelescopeX[telIdx] );
+                  fEvent->AddXY(xArrOff + fTelescopeY[telIdx], 
+                                 yArrOff - fTelescopeX[telIdx]);
+                  fEvent->SimWavelength(fRunHeader->GetWavelengthMin(), 
+                                          fRunHeader->GetWavelengthMax());
+                  }
+               else
+                  // skip this telescope
+                  fInFormat->Seek(fBlockLength);
+               }
+               break;
+
+            case 1209:  // the event end
+               status = fEvtHeader->ReadEvtEnd(fInFormat);
+               
+               if (fReadState == 10 || fReadState == 2)
+                  {
+                  // corsika events
+                  fEvtHeader->ResetNumReuse();
+                  fEvtHeader->InitXY();
+                  fBlockType = 1109;  // save corsika events
+                  continue;
+                  }
+
+               fReadState = 3;  // event closed, run still open
+               break;
+
+
+            case 1210:  // the run end
+               status = fRunHeader->ReadEvtEnd(fInFormat, kTRUE);
+               fNumEvents += fRunHeader->GetNumEvents();
+               fRunHeader->SetReadyToSave();
+               fReadState = 5;           // back to  starting point
+               // this may be the last block in the file. 
+               // We exit to write the header into the file
+               fBlockLength = 0;
+               fBlockType   = 0; // not available type 
+               return kTRUE;
+               break;
+
+
+            case 1105:  // event block of raw format
+               if (fReadState == 2 || fReadState == 10)
+                  {
+                  Int_t oldSize = fRawEvemtBuffer.size();
+                  fRawEvemtBuffer.resize(oldSize + fBlockLength / sizeof(Float_t));
+                  status = fInFormat->Read(&fRawEvemtBuffer[oldSize], fBlockLength); 
+                  fReadState = 10;
+                  }
+               else
+                  fInFormat->Seek(fBlockLength);
+               break;
+
+
+            case 1109:  // save corsika events
+               fEvtHeader->InitXY();
+               status = fEvent->ReadCorsikaEvt(&fRawEvemtBuffer[0], 
+                                               fRawEvemtBuffer.size() / 7, 
+                                               fEvtHeader->GetNumReuse()+1);
+               fEvtHeader->IncNumReuse();
+
+               if (fEvtHeader->GetNumReuse() == fRunHeader->GetNumReuse())
+                  {
+                  // this was the last reuse. Set fBlockType to EVTE to save
+                  // it the next time.
+                  fRawEvemtBuffer.resize(0);
+
+                  fReadState = 3;
+                  fBlockType = 1209;
+                  }
+               else
+                  // save this reuse
+                  return status;
+
+               break;               
+
+            default:
+               // unknown block, we skip it
+               fReadState = 15;
+               fInFormat->Seek(fBlockLength);
+
+            }
+
+         if (status != kTRUE)
+            // there was an error while data were read
+            return status;
+
+         Int_t headerStatus =  ReadNextBlockHeader();
+         if (headerStatus == kFALSE) 
+            // end of file
+            break;
+         if (headerStatus == kTRUE)
+            // we shall quit this loop
+            return kTRUE;
+
+            // else continue
+         }
+
+      }
+
+      //
+      // If an event could not be read from file try to open new file
+      //
+      const Int_t rc = OpenNextFile();
+      if (rc!=kTRUE)
+          return rc;
+
+      if (ReadNextBlockHeader() == kFALSE)
+         return kFALSE;
+   }
+   return kTRUE;
+}
 
 // --------------------------------------------------------------------------
@@ -426 +673 @@
 //  - the raw data information of one event is read
 //
-Int_t MCorsikaRead::Process()
+Int_t MCorsikaRead::Process_Old()
 {
 
@@ -453 +700 @@
       if (fInFormat)
       {
-         Int_t blockType, blockVersion, blockIdentifier, blockLength;
-
-         while (fInFormat->NextBlock(fReadState == 4, blockType, blockVersion, 
-                                     blockIdentifier, blockLength))
+
+         while (fInFormat->NextBlock(fReadState == 4, fBlockType, fBlockVersion, 
+                                     fBlockIdentifier, fBlockLength))
          {
-//            gLog << "Next block:  type=" << blockType << " version=" << blockVersion;
-//            gLog << " identifier=" << blockIdentifier << " length=" << blockLength << endl;
+            gLog << "Next fBlock:  type=" << fBlockType << " version=" << fBlockVersion;
+            gLog << " identifier=" << fBlockIdentifier << " length=" << fBlockLength << endl;
 
 
             if (fReadState == 4)
                {
-               fTopBlockLength -= blockLength + 12;
+               fTopBlockLength -= fBlockLength + 12;
                if (fTopBlockLength <= 0)
                   // all data of a top block are read, go back to normal state
@@ -471 +717 @@
 
             Int_t status = kTRUE;
-            switch (blockType)
+            switch (fBlockType)
                {
                case 1200:       // the run header
@@ -480 +726 @@
                case 1201:       // telescope position
                   status = ReadTelescopePosition();
+sleep(5);
                   break;
 
@@ -507 +754 @@
 
                case 1204:
-                  if (fArrayIdx < 0 || fArrayIdx == blockIdentifier)
+                  if (fArrayIdx < 0 || fArrayIdx == fBlockIdentifier)
                      { 
                      fReadState = 4;
-                     fTopBlockLength = blockLength;
+                     fTopBlockLength = fBlockLength;
                      }
                   else
                      // skip this array of telescopes
-                     fInFormat->Seek(blockLength);
+                     fInFormat->Seek(fBlockLength);
 
                   break;
@@ -520 +767 @@
                case 1205:
                   {
-                  Int_t telIdx   = blockIdentifier % 1000;
-                  if (blockVersion == 0                               &&
-                      (fTelescopeIdx < 0 || fTelescopeIdx ==  telIdx) &&
-                       blockLength > 12)
+                  Int_t telIdx   = fBlockIdentifier % 1000;
+                  if (fBlockVersion == 0                               &&
+                      (fTelescopeIdx < 0 || fTelescopeIdx ==  telIdx)     )
                      {
                      status = fEvent->ReadEventIoEvt(fInFormat);
 
-                     Int_t arrayIdx = blockIdentifier / 1000;
+                     Int_t arrayIdx = fBlockIdentifier / 1000;
                      Float_t xArrOff, yArrOff;
                      fEvtHeader->GetArrayOffset(arrayIdx, xArrOff, yArrOff);
@@ -545 +791 @@
                   else
                      // skip this telescope
-                     fInFormat->Seek(blockLength);
+                     fInFormat->Seek(fBlockLength);
                   }
                   break;
@@ -551 +797 @@
                case 1209:  // the event end
                   status = fEvtHeader->ReadEvtEnd(fInFormat);
-                  if (fReadState == 10)
+                  if (fReadState == 10 || fReadState == 2)
                      {
                      // all particles of this event are read, now save them
@@ -580 +826 @@
                      {
                      Int_t oldSize = fRawEvemtBuffer.size();
-                     fRawEvemtBuffer.resize(oldSize + blockLength / sizeof(Float_t));
-                     status = fInFormat->Read(&fRawEvemtBuffer[oldSize], blockLength); 
+                     fRawEvemtBuffer.resize(oldSize + fBlockLength / sizeof(Float_t));
+                     status = fInFormat->Read(&fRawEvemtBuffer[oldSize], fBlockLength); 
                      fReadState = 10;
                      }
                   else
-                     fInFormat->Seek(blockLength);
+                     fInFormat->Seek(fBlockLength);
                   break;
 
                default:
                   // unknown block, we skip it
-                  fInFormat->Seek(blockLength);
+                  fInFormat->Seek(fBlockLength);
 
                }
@@ -629 +875 @@
 
     *fLog << warn << dec;
-    *fLog << "Warning - number of read events (" << GetNumExecutions()-2;
+    *fLog << "Warning - number of read events (" << GetNumExecutions()-1;
     *fLog << ") doesn't match number in run header(s) (";
     *fLog << fRunHeader->GetNumReuse() << "*" << fNumEvents << "=" << n << ")." << endl;
 
+
     return kTRUE;
 }

trunk/Mars/mcorsika/MCorsikaRead.h

@@ -48 +48 @@
     TArrayF fTelescopeR;       //! Radii of spheres around tel. (unit: cm)
 
+    Int_t   fBlockType;        //! block type of the just read block header
+    Int_t   fBlockVersion;     //! block version of the just read block header 
+    Int_t   fBlockIdentifier;  //! block identifier of the just read block header
+    Int_t   fBlockLength;      //! block length from the just read block header
+
     Int_t   fReadState;        //! 0: nothing read yet 
                                //  1: RUNH is already read
@@ -55 +60 @@
                                //  5: RUNE is already read
                                // 10: raw events are currently read
-                               // 11: raw events are currently saved
+                               // 15: undefined status
     Int_t   fTopBlockLength;   // remaining length of the current top-level block 1204
 
@@ -66 +71 @@
     Bool_t CalcNumTotalEvents();
     Int_t  ReadTelescopePosition();
+    Int_t  ReadNextBlockHeader();
 
     // MTask
     Int_t PreProcess(MParList *pList);
     Int_t Process();
+    Int_t Process_Old();
     Int_t PostProcess();