Index: /trunk/Mars/mcorsika/MCorsikaRead.cc =================================================================== --- /trunk/Mars/mcorsika/MCorsikaRead.cc (revision 10146) +++ /trunk/Mars/mcorsika/MCorsikaRead.cc (revision 10147) @@ -249,4 +249,9 @@ } +// -------------------------------------------------------------------------- +// +// Calculates the total number of events from all input files and stores +// the number in fNumTotalEvents. +// Bool_t MCorsikaRead::CalcNumTotalEvents() { @@ -318,4 +323,8 @@ } +// -------------------------------------------------------------------------- +// +// Reads the the position of all telescopes in one array +// Int_t MCorsikaRead::ReadTelescopePosition() { @@ -352,5 +361,5 @@ for (int i=lo; iNextBlock(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,4 +465,203 @@ } +// -------------------------------------------------------------------------- +// +// 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,5 +673,5 @@ // - the raw data information of one event is read // -Int_t MCorsikaRead::Process() +Int_t MCorsikaRead::Process_Old() { @@ -453,16 +700,15 @@ 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,5 +717,5 @@ Int_t status = kTRUE; - switch (blockType) + switch (fBlockType) { case 1200: // the run header @@ -480,4 +726,5 @@ case 1201: // telescope position status = ReadTelescopePosition(); +sleep(5); break; @@ -507,12 +754,12 @@ 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,12 +767,11 @@ 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,5 +791,5 @@ else // skip this telescope - fInFormat->Seek(blockLength); + fInFormat->Seek(fBlockLength); } break; @@ -551,5 +797,5 @@ 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,15 +826,15 @@ { 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,8 +875,9 @@ *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; } Index: /trunk/Mars/mcorsika/MCorsikaRead.h =================================================================== --- /trunk/Mars/mcorsika/MCorsikaRead.h (revision 10146) +++ /trunk/Mars/mcorsika/MCorsikaRead.h (revision 10147) @@ -48,4 +48,9 @@ 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,5 +60,5 @@ // 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,8 +71,10 @@ 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();